Stability of conjugate gradient and Lanczos methods for linear least squares problems.

*(English)*Zbl 0924.65035A thorough study of two different implementations of the conjugate gradient and the Lanczos methods is presented, in order to find numerical solutions to linear least squares problems. Although these implementations are mathematically equivalent, their computational, and most important, numerical properties are quite different. For the conjugate gradient method the implementations considered are the original ones, as given by M. R. Hestenes and E. Stiefel in their classic paper [(CGLS1) J. Res. Nat. Bur. Standards 49, 409-436 (1952; Zbl 0048.09901)] and a variant of it (CGLS2) obtained by recurring on the residuals of the normal equations, instead of doing so on the original residuals. For the Lanczos method, the two implementations LSQR and LSCG differ in that LSQR employs Givens rotations in the bidiagonalization process.

By adapting a result obtained by A. Greenbaum [SIAM J. Matrix Anal. App. 18, No. 3, 535-551 (1997; Zbl 0873.65027)] for linear problems with square matrices, bounds are obtained for the errors in the computed residuals for both CGLS1 and LSQR. These bounds show that both implementations are normwise backward stable. A heuristic explanation is also given for the failure of both CGLS2 and LSCG to yield accurate results. As for storage and operational count, the unstable implementations are less costly than their stable counterparts. But for the tests reported in the paper, which deal with ill-conditioned problems, the difference in accuracy is overwhelmingly in favour of the stable implementations.

Another important point to bear in mind is that the unstable variants tend to converge faster (to inaccurate results) than CGLS1, whereas the situation is even when compared with LSQR, which almost always converges faster than CGLS1. On the other hand, regarding accuracy, CGLS1 beats them all, with a slightly higher a priori computational cost than CGLS2, and a significantly smaller one than both LSQR and LSCG. These are perhaps the reasons for the authors’ preference for CGLS1.

The paper does not report on effective computing times that would justify even more this preference, especially with respect to LSQR, although such justification is suggested by the included graphs. Another very nice graph shows for one example that the norm of the residuals of the normal equations decrease for CGLS1 to about machine precision, while for CGLS2, in spite of very poor convergence, they decrease to zero in the limit. The last conclusion to point out is that for the inconsistent case, CGLS1 and LSQR achieve even better accuracy than a normwise backward stable method.

By adapting a result obtained by A. Greenbaum [SIAM J. Matrix Anal. App. 18, No. 3, 535-551 (1997; Zbl 0873.65027)] for linear problems with square matrices, bounds are obtained for the errors in the computed residuals for both CGLS1 and LSQR. These bounds show that both implementations are normwise backward stable. A heuristic explanation is also given for the failure of both CGLS2 and LSCG to yield accurate results. As for storage and operational count, the unstable implementations are less costly than their stable counterparts. But for the tests reported in the paper, which deal with ill-conditioned problems, the difference in accuracy is overwhelmingly in favour of the stable implementations.

Another important point to bear in mind is that the unstable variants tend to converge faster (to inaccurate results) than CGLS1, whereas the situation is even when compared with LSQR, which almost always converges faster than CGLS1. On the other hand, regarding accuracy, CGLS1 beats them all, with a slightly higher a priori computational cost than CGLS2, and a significantly smaller one than both LSQR and LSCG. These are perhaps the reasons for the authors’ preference for CGLS1.

The paper does not report on effective computing times that would justify even more this preference, especially with respect to LSQR, although such justification is suggested by the included graphs. Another very nice graph shows for one example that the norm of the residuals of the normal equations decrease for CGLS1 to about machine precision, while for CGLS2, in spite of very poor convergence, they decrease to zero in the limit. The last conclusion to point out is that for the inconsistent case, CGLS1 and LSQR achieve even better accuracy than a normwise backward stable method.

Reviewer: Juan Pedro Milaszewicz (Ioannina)

##### MSC:

65F20 | Numerical solutions to overdetermined systems, pseudoinverses |

65F10 | Iterative numerical methods for linear systems |