  For the large-scale linear discrete ill-posed problem $\min\|Ax-b\|$ or
$Ax=b$ with $b$ contaminated by white noise, the Golub-Kahan bidiagonalization
based LSQR method and its mathematically equivalent CGLS, the Conjugate
Gradient (CG) method applied to $A^TAx=A^Tb$, are most commonly used. They have
intrinsic regularizing effects, where the iteration number $k$ plays the role
of regularization parameter. The long-standing fundamental question is: {\em
Can LSQR and CGLS find 2-norm filtering best possible regularized solutions}?
The author has given definitive answers to this question for severely and
moderately ill-posed problems when the singular values of $A$ are simple. This
paper extends the results to the multiple singular value case, and studies the
approximation accuracy of Krylov subspaces, the quality of low rank
approximations generated by Golub-Kahan bidiagonalization and the convergence
properties of Ritz values. For the two kinds of problems, we prove that LSQR
finds 2-norm filtering best possible regularized solutions at semi-convergence.
Particularly, we consider some important and untouched issues on best, near
best and general rank $k$ approximations to $A$ for the ill-posed problems with
the singular values $\sigma_k=\mathcal{O}(k^{-\alpha})$ with $\alpha>0$, and
the relationships between them and their nonzero singular values. Numerical
experiments confirm our theory. The results on general rank $k$ approximations
and the properties of their nonzero singular values apply to several Krylov
solvers, including LSQR, CGME, MINRES, MR-II, GMRES and RRGMRES.
