`lsqr`

Solves a least-squares or linear system iteratively.

Syntax

x = lsqr(A, b)

Solves A*x = b, in the least-squares sense when A is rectangular.

x = lsqr(A, b, tol, maxit)

Sets the tolerance and iteration limit.

x = lsqr(A, b, tol, maxit, M1, M2, x0)

Adds a preconditioner and an initial guess.

[x, flag, relres, iter, resvec, lsvec] = lsqr(___)

Also returns the convergence flag, residual information and history vectors.

Description

LSQR applies the Lanczos bidiagonalization to the normal equations A'*A*x = A'*b without ever forming A'*A, which would square the condition number. Each iteration costs one product with A and one with A'.

It is intended for large sparse least-squares problems. Because the default maxit is only 20 for large matrices, raise it when the flag is 1.

Input arguments

A: Coefficient matrix, not necessarily square, or a function handle.
b: Right-hand side column vector.
tol: Tolerance; default 1e-6.
maxit: Maximum iterations; default min(size(A,1), size(A,2), 20).

Output arguments

x: The solution, a column vector with size(A,2) entries.
flag: 0 if converged, 1 if maxit was reached, 2 for an ill-conditioned preconditioner, 3 if stagnated.
relres: Relative residual norm(b - A*x)/norm(b).

Example

Square system

A = sparse(diag(4*ones(5,1)) + diag(ones(4,1),1) + diag(ones(4,1),-1));
b = (1:5)';
[x, flag, relres] = lsqr(A, b, 1e-10, 100);
fprintf('flag = %d, relres < 1e-10: %d\n', flag, relres < 1e-10)
disp(x')

Output

flag = 0, relres < 1e-10: 1
    0.1679    0.3282    0.5192    0.5949    1.1013

Notes

Rectangular systems may return a solution of the wrong size in this version; use square matrices, or A\b, for least squares.

lsqr