Syntax

Description

x = lsqr(A,b)attempts to solve the system of linear equationsA*x=bforxifAis consistent, otherwise it attempts to solve the least squares solutionxthat minimizesnorm(b-A*x). Them-by-ncoefficient matrixAneed not be square but it should be large and sparse. The column vectorbmust have lengthm. You can specifyAas a function handle,afun, such thatafun(x,'notransp')returnsA*xandafun(x,'transp')returnsA'*x.

Parameterizing Functionsexplains how to provide additional parameters to the functionafun, as well as the preconditioner functionmfundescribed below, if necessary.

Iflsqrconverges, a message to that effect is displayed. Iflsqrfails to converge after the maximum number of iterations or halts for any reason, a warning message is printed displaying the relative residualnorm(b-A*x)/norm(b)and the iteration number at which the method stopped or failed.

lsqr(A,b,tol)specifies the tolerance of the method. Iftolis[], thenlsqruses the default,1e-6.

lsqr(A,b,tol,maxit)specifies the maximum number of iterations.

lsqr(A,b,tol,maxit,M)andlsqr(A,b,tol,maxit,M1,M2)usen-by-npreconditionerMorM = M1*M2and effectively solve the systemA*inv(M)*y = bfory, wherey = M*x. IfMis[]thenlsqrapplies no preconditioner.Mcan be a functionmfunsuch thatmfun(x,'notransp')returnsM\xandmfun(x,'transp')returnsM'\x.

lsqr(A,b,tol,maxit,M1,M2,x0)specifies then-by-1initial guess. Ifx0is[], thenlsqruses the default, an all zero vector.

[x,flag] = lsqr(A,b,tol,maxit,M1,M2,x0)also returns a convergence flag.

Wheneverflagis not0, the solutionxreturned is that with minimal norm residual computed over all the iterations. No messages are displayed if you specify theflagoutput.

[x,flag,relres] = lsqr(A,b,tol,maxit,M1,M2,x0)also returns an estimate of the relative residualnorm(b-A*x)/norm(b). Ifflagis0,relres <= tol.

[x,flag,relres,iter] = lsqr(A,b,tol,maxit,M1,M2,x0)also returns the iteration number at whichxwas computed, where0 <= iter <= maxit.

[x,flag,relres,iter,resvec] = lsqr(A,b,tol,maxit,M1,M2,x0)also returns a vector of the residual norm estimates at each iteration, includingnorm(b-A*x0).

[x,flag,relres,iter,resvec,lsvec] = lsqr(A,b,tol,maxit,M1,M2,x0)also returns a vector of estimates of the scaled normal equations residual at each iteration:norm((A*inv(M))'*(B-A*X))/norm(A*inv(M),'fro'). Note that the estimate ofnorm(A*inv(M),'fro')changes, and hopefully improves, at each iteration.

n = 100; on = ones(n,1); A = spdiags([-2*on 4*on -on],-1:1,n,n); b = sum(A,2); tol = 1e-8; maxit = 15; M1 = spdiags([on/(-2) on],-1:0,n,n); M2 = spdiags([4*on -on],0:1,n,n); x = lsqr(A,b,tol,maxit,M1,M2);

lsqr converged at iteration 11 to a solution with relative residual 3.5e-009

function x1 = run_lsqr n = 100; on = ones(n,1); A = spdiags([-2*on 4*on -on],-1:1,n,n); b = sum(A,2); tol = 1e-8; maxit = 15; M1 = spdiags([on/(-2) on],-1:0,n,n); M2 = spdiags([4*on -on],0:1,n,n); x1 = lsqr(@afun,b,tol,maxit,M1,M2); function y = afun(x,transp_flag) if strcmp(transp_flag,'transp') % y = A'*x y = 4 * x; y(1:n-1) = y(1:n-1) - 2 * x(2:n); y(2:n) = y(2:n) - x(1:n-1); elseif strcmp(transp_flag,'notransp') % y = A*x y = 4 * x; y(2:n) = y(2:n) - 2 * x(1:n-1); y(1:n-1) = y(1:n-1) - x(2:n); end end end

x1=run_lsqr;

lsqr converged at iteration 11 to a solution with relative residual 3.5e-009

References

See Also

bicg|bicgstab|cgs|gmres|minres|norm|pcg|qmr|symmlq