  For the large-scale linear discrete ill-posed problem $\min\|Ax-b\|$ or
$Ax=b$ with $b$ contaminated by a white noise, the Lanczos bidiagonalization
based LSQR method and its mathematically equivalent Conjugate Gradient (CG)
method for $A^TAx=A^Tb$ are most commonly used. They have intrinsic
regularizing effects, where the number $k$ of iterations plays the role of
regularization parameter. However, there has been no answer to the
long-standing fundamental concern by Bj\"{o}rck and Eld\'{e}n in 1979: 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 standard-form Tikhonov regularization. In
this paper, assuming that the singular values of $A$ are simple, we analyze the
regularization of LSQR for severely, moderately and mildly ill-posed problems.
We establish accurate estimates for the 2-norm distance between the underlying
$k$-dimensional Krylov subspace and the $k$-dimensional dominant right singular
subspace of $A$. For the first two kinds of problems, we then prove that LSQR
finds a best possible regularized solution at semi-convergence occurring at
iteration $k_0$ and that, for $k=1,2,\ldots,k_0$, (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
in natural order; (iii) the $k$-step LSQR always captures the $k$ dominant SVD
components of $A$. For the third kind of problem, we prove that LSQR generally
cannot find a best possible regularized solution. Numerical experiments confirm
our theory.
