lsqr

Solve a rectangular linear system using lsqr with default settings, and then adjust the tolerance and number of iterations used in the solution process.

Create a random sparse matrix A with 50% density. Also create a random vector b for the right-hand side of $\mathrm{Ax}=\mathit{b}$.

rng default
A = sprand(400,300,.5);
b = rand(400,1);

Solve $\mathrm{Ax}=\mathit{b}$ using lsqr. The output display includes the value of the relative residual error $\frac{\Vert \mathit{b}-\mathrm{Ax}\Vert }{\Vert \mathit{b}\Vert }$.

x = lsqr(A,b);

lsqr stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 0.26.

By default lsqr uses 20 iterations and a tolerance of 1e-6, but the algorithm is unable to converge in those 20 iterations for this matrix. Since the residual is still large, it is a good indicator that more iterations (or a preconditioner matrix) are needed. You also can use a larger tolerance to make it easier for the algorithm to converge.

Solve the system again using a tolerance of 1e-4 and 70 iterations. Specify six outputs to return the relative residual relres of the calculated solution, as well as the residual history resvec and the least-squares residual history lsvec.

[x,flag,relres,iter,resvec,lsvec] = lsqr(A,b,1e-4,70);
flag

flag = 
0

Since flag is 0, the algorithm was able to meet the desired error tolerance in the specified number of iterations. You can generally adjust the tolerance and number of iterations together to make trade-offs between speed and precision in this manner.

Examine the relative residual and least-squares residual of the calculated solution.

relres

relres = 
0.2625

lsres = lsvec(end)

lsres = 
2.5375e-04

These residual norms indicate that x is a least-squares solution, because relres is not smaller than the specified tolerance of 1e-4. Since no consistent solution to the linear system exists, the best the solver can do is to make the least-squares residual satisfy the tolerance.

Plot the residual histories. The relative residual resvec quickly reaches a minimum and cannot make further progress, while the least-squares residual lsvec continues to be minimized on subsequent iterations.

N = length(resvec);
semilogy(0:N-1,lsvec,'--o',0:N-1,resvec,'-o')
legend("Least-squares residual","Relative residual")

Examine the effect of using a preconditioner matrix with lsqr to solve a linear system.

Load west0479, a real 479-by-479 nonsymmetric sparse matrix.

load west0479
A = west0479;

Define b so that the true solution to $\mathrm{Ax}=\mathit{b}$ is a vector of all ones.

b = sum(A,2);

Set the tolerance and maximum number of iterations.

tol = 1e-12;
maxit = 20;

Use lsqr to find a solution at the requested tolerance and number of iterations. Specify six outputs to return information about the solution process:

[x,fl,rr,it,rv,lsrv] = lsqr(A,b,tol,maxit);
fl

fl = 
1

rr

rr = 
0.0017

it

it = 
20

Since fl = 1, the algorithm did not converge to the specified tolerance within the maximum number of iterations.

To aid with the slow convergence, you can specify a preconditioner matrix. Since A is nonsymmetric, use ilu to generate the preconditioner $\mathit{M}=\mathit{L}\mathit{U}$ in factorized form. Specify a drop tolerance to ignore nondiagonal entries with values smaller than 1e-6. Solve the preconditioned system ${\mathrm{AM}}^{-1}\left(\mathit{M}\mathit{x}\right)=\mathit{b}$ for $\mathit{y}=\mathrm{Mx}$ by specifying L and U as the M1 and M2 inputs to lsqr.

setup = struct('type','ilutp','droptol',1e-6);
[L,U] = ilu(A,setup);
[x1,fl1,rr1,it1,rv1,lsrv1] = lsqr(A,b,tol,maxit,L,U);
fl1

fl1 = 
0

rr1

rr1 = 
7.1047e-14

it1

it1 = 
13

The use of an ilu preconditioner produces a relative residual less than the prescribed tolerance of 1e-12 at the 13th iteration. The output rv1(1) is norm(b), and the output rv1(end) is norm(b-A*x1).

You can follow the progress of lsqr by plotting the relative residuals at each iteration. Plot the residual history of each solution with a line for the specified tolerance.

semilogy(0:length(rv)-1,rv/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Examine the effect of supplying lsqr with an initial guess of the solution.

Create a random rectangular sparse matrix. Use the sum of each row as the vector for the right-hand side of $\mathrm{Ax}=\mathit{b}$ so that the expected solution for $\mathit{x}$ is a vector of ones.

A = sprand(700,900,0.1);
b = sum(A,2);

Use lsqr to solve $\mathrm{Ax}=\mathit{b}$ twice: one time with the default initial guess, and one time with a good initial guess of the solution. Use 75 iterations and the default tolerance for both solutions. Specify the initial guess in the second solution as a vector with all elements equal to 0.99.

maxit = 75;
x1 = lsqr(A,b,[],maxit);

lsqr converged at iteration 64 to a solution with relative residual 8.7e-07.

x0 = 0.99*ones(size(A,2),1);
x2 = lsqr(A,b,[],maxit,[],[],x0);

lsqr converged at iteration 26 to a solution with relative residual 9.6e-07.

With an initial guess close to the expected solution, lsqr is able to converge in fewer iterations.

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = lsqr(A,b,tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

Solve a linear system by providing lsqr with a function handle that computes A*x and A'*x in place of the coefficient matrix A.

Create a nonsymmetric tridiagonal matrix. Preview the matrix.

A = gallery('wilk',21) + diag(ones(20,1),1)

A = 21×21

    10     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     2     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     2     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     2     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     2     0     0     0     0     0     0     0     0     0     0
      ⋮

Since this tridiagonal matrix has a special structure, you can represent the operation A*x with a function handle. When A multiplies a vector, most of the elements in the resulting vector are zeros. The nonzero elements in the result correspond with the nonzero tridiagonal elements of A.

The expression $\mathit{A}\mathit{x}$ becomes:

$\mathit{A}\mathit{x}=\left[\begin{array}{ccccccc}10& 2& 0& \cdots & & \cdots & 0\\ 1& 9& 2& 0& & & ⋮\\ 0& 1& \ddots & 2& 0& & \\ ⋮& 0& 1& 0& \ddots & \ddots & ⋮\\ & & 0& \ddots & 1& \ddots & 0\\ ⋮& & & \ddots & \ddots & \ddots & 2\\ 0& \cdots & & \cdots & 0& 1& 10\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_3\\ ⋮\\ ⋮\\ {\mathit{x}}_{21}\end{array}\right]=\left[\begin{array}{c}{10\mathit{x}}_1+2{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+2{\mathit{x}}_3\\ ⋮\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+2{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$.

The resulting vector can be written as the sum of three vectors:

$\mathit{A}\mathit{x}=\left[\begin{array}{c}{10\mathit{x}}_1+2{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+2{\mathit{x}}_3\\ ⋮\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+2{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$=$\left[\begin{array}{c}0\\ {\mathit{x}}_1\\ {\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}{10\mathit{x}}_1\\ {9\mathit{x}}_2\\ ⋮\\ 9{\mathit{x}}_{20}\\ 10{\mathit{x}}_{21}\end{array}\right]+2\cdot \left[\begin{array}{c}{\mathit{x}}_2\\ {\mathit{x}}_3\\ ⋮\\ {\mathit{x}}_{21}\\ 0\end{array}\right]$.

Likewise, the expression for ${\mathit{A}}^{\mathit{T}}\mathit{x}$ becomes:

${\mathit{A}}^{\mathit{T}}\mathit{x}=\left[\begin{array}{ccccccc}10& 1& 0& \cdots & & \cdots & 0\\ 2& 9& 1& 0& & & ⋮\\ 0& 2& \ddots & 1& 0& & \\ ⋮& 0& 2& 0& \ddots & \ddots & ⋮\\ & & 0& \ddots & 1& \ddots & 0\\ ⋮& & & \ddots & \ddots & \ddots & 1\\ 0& \cdots & & \cdots & 0& 2& 10\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_3\\ ⋮\\ ⋮\\ {\mathit{x}}_{21}\end{array}\right]=\left[\begin{array}{c}{10\mathit{x}}_1+{\mathit{x}}_2\\ {2\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ ⋮\\ ⋮\\ {2\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {2\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$.

${\mathit{A}}^{\mathit{T}}\mathit{x}=\left[\begin{array}{c}{10\mathit{x}}_1+{\mathit{x}}_2\\ {2\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ ⋮\\ ⋮\\ {2\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {2\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]=2\cdot \left[\begin{array}{c}0\\ {\mathit{x}}_1\\ {\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}{10\mathit{x}}_1\\ {9\mathit{x}}_2\\ ⋮\\ 9{\mathit{x}}_{20}\\ 10{\mathit{x}}_{21}\end{array}\right]+\left[\begin{array}{c}{\mathit{x}}_2\\ {\mathit{x}}_3\\ ⋮\\ {\mathit{x}}_{21}\\ 0\end{array}\right]$.

In MATLAB®, write a function that creates these vectors and adds them together, thus giving the value of A*x or A'*x, depending on the flag input:

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end

(This function is saved as a local function at the end of the example.)

Now, solve the linear system $\mathrm{Ax}=\mathit{b}$ by providing lsqr with the function handle that calculates A*x and A'*x. Use a tolerance of 1e-6 and 25 iterations. Specify $\mathit{b}$ as the row sums of $\mathit{A}$ so that the true solution for $\mathit{x}$ is a vector of ones.

b = full(sum(A,2));
tol = 1e-6;  
maxit = 25;
x1 = lsqr(@afun,b,tol,maxit)

lsqr converged at iteration 21 to a solution with relative residual 4.5e-12.

x1 = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end