There has been recently a trend to study linear system identification with high order finite impulse response (FIR) models using the regularized least-squares approach. One key of this approach is to solve the hyper-parameter estimation problem that is usually nonconvex. Our goal here is to investigate implementation of algorithms for solving the hyper-parameter estimation problem that can deal with both large data sets and possibly ill-conditioned computations. In particular, a QR factorization based matrix-inversion-free algorithm is proposed to evaluate the cost function in an efficient and accurate way. It is also shown that the gradient and Hessian of the cost function can be computed based on the same QR factorization. Finally, the proposed algorithm and ideas are verified by Monte-Carlo simulations on a large data-bank of test systems and data sets. (C) 2013 Elsevier Ltd. All rights reserved.