Iterative algorithms for least squares problems with multiple right-hand sides

Rasmus Munk Larsen

HEPL Annex A206, Stanford University, Stanford, CA 94305-4085, USA


We discuss iterative algorithms for solving linear least squares problems with multiple right-hand sides when the coefficient matrix is large, ill-conditioned and sparse or structured.

The algorithms presented are all based on the Lanczos bidiagonalization process, but different implementations allow significant trade-off between the number of operations, required amount of storage, parallel efficiency and fraction of BLAS level 3 operations involved. We study the performance of three algorithms:

  1. The (block-) conjugate gradient algorithm applied to the normal equations (CGLS) augmented with a Galerkin projection step.
  2. A projection algorithm based on (block-) Lanczos bidiagonalization with partial reorthogonalization.
  3. Solving for each right-hand side independently by means of the LSQR algorithm.

As an example, ill-conditioned least squares problems with multiple right-hand sides arise in image deblurring when restoring multiple frames distorted by the same point spread function or when implementing the Backus-Gilbert (BG) inversion method. The BG method is used extensively in seismology, helioseismology, astronomy and other areas involving inverse problems. We present numerical results for examples from image deblurring and helioseismology, where BG implementations based on the three algorithms above were used to solve the 2D inverse problem of the solar rotation.