<p>A greedy orthogonal least squares identification algorithm based on the Householder transformation is proposed for multi-variable Hammerstein models with unknown time delays. Initially, the model is overparametrized into a sparse identification model with a high dimensional sparse vector by introducing the upper bounds of the nonlinear order and a data regression length. Subsequently, the QR decomposition of the high-dimensional information matrix is carried out by the Householder transformation, which circumvents ill-conditioned solutions and significantly reduces the computational burden. Then a modified stopping rule based on the Akaike information criterion is employed to determine the sparsity level. Finally, the time delays, orders and parameters estimation are estimated from the sparse parameter vector. Numerical simulations and a dynamic continuous stirred-tank reactor process demonstrate that, compared to other identification algorithms, the proposed algorithm offers higher accuracy and better anti-interference performance.</p>