SOL Logo

Systems Optimization Laboratory

Stanford University
Dept of Management Science and Engineering (MS&E)

Huang Engineering Center

Stanford, CA 94305-4121  USA

CGLS: CG method for Ax = b and Least Squares

  • AUTHOR: Michael Saunders
    CONTRIBUTORS: Per Christian Hansen, Folkert Bleichrodt, Christopher Fougner

    A MATLAB implementation of CGLS, the Conjugate Gradient method for unsymmetric linear equations and least squares problems: \begin{align*} \text{Solve } & Ax=b \\ \text{or minimize } & \|Ax-b\|^2 \\ \text{or solve } & (A^T A + sI)x = A^T b, \end{align*} where the matrix \(A\) may be square or rectangular (represented by an M-file for computing \(Ax\) and \(A^Tx\)) and \(s\) is a scalar (positive or negative).
    The method is stable if \(s = 0\) or \(s > 0\). More generally, it should be stable if \(A^T A + sI\) is positive definite. Otherwise it may be unstable.

    Documented in the MATLAB file below.
    Special feature: This is a simple CG-type code for unsymmetric equations and least squares, with the option of a negative shift. However, use with caution if \(s < 0\).

    If \(s \ge 0\), we recommend LSQR or LSMR. Do not use a symmetric CG code on the normal equations, even though it is equivalent to LSQR in exact arithmetic.

    With \(s = 0\), the method is due to Hestenes and Stiefel (1952). See the first LSQR paper, Paige and Saunders (1982a), ACM TOMS 8(1), 43-71. Trivial changes are needed to allow an arbitrary shift.

    01 Sep 1999: First version, implemented at DTU, Denmark.
    15 Apr 2013: More polished version contributed by Folkert Bleichrodt, Centrum Wiskunde & Informatica (CWI), Amsterdam, The Netherlands.
    02 Sep 2014: CUDA implementation contributed by Christopher Fougner, ICME, Stanford.
    08 Feb 2015: CUDA version can treat \(A\) as a matrix or an operator.
    09 Feb 2015: CUDA version allows \(A,b,x\) to be complex (but \(s\) is real).