%% Iterative Solution to Linear System
% 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 \textrm{Ax}=b.

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

%
% Solve \textrm{Ax}=b using lsqr. The output display includes the value of the relative residual error \frac{\left\|b-\textrm{Ax}\right\|}{\left\|b\right\|}.

x = lsqr(A,b);

% 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

% 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
lsres = lsvec(end)

% 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")

%
% Copyright 2019 The MathWorks, Inc.