  For the large-scale linear discrete ill-posed problem $\min\|Ax-b\|$ or
$Ax=b$ with $b$ contaminated by a white noise, Lanczos bidiagonalization based
LSQR and its mathematically equivalent CGLS are most commonly used. They have
intrinsic regularizing effects, where the number $k$ of iterations plays the
role of regularization parameter. However, hitherto there has been no answer to
the long-standing fundamental concern of Bj\"{o}rck and Eld\'{e}n in 1979: {\em
for which kinds of problems LSQR and CGLS can find best possible regularized
solutions}? Here a best possible regularized solution means that it is at least
as accurate as the best regularized solution obtained by the truncated singular
value decomposition (TSVD) method or by standard-form Tikhonov regularization
and cannot be improved under certain conditions. In this paper we make a
detailed analysis on the regularization of LSQR for severely, moderately and
mildly ill-posed problems. For the first two kinds of problems, we prove that
LSQR finds best possible solutions at semi-convergence and the following
results hold until semi-convergence: (i) the $k$-step Lanczos bidiagonalization
always generates a near best rank $k$ approximation to $A$; (ii) the $k$ Ritz
values always approximate the first $k$ large singular values of $A$ in natural
order; (iii) the $k$-step LSQR always captures the $k$ dominant SVD components
of $A$; (iv) the diagonals and subdiagonals of the bidiagonal matrices
generated by Lanczos bidiagonalization decay as fast as the singular values of
$A$. However, for the third kind of problem, the above results do not hold
generally. We also analyze the regularization of the other two Krylov solvers
LSMR and CGME, proving that LSMR has similar regularizing effects to LSQR for
each kind of problem and both are superior to CGME. Numerical experiments
confirm our theory on LSQR.
