%% Specify Tolerance to Reduce Impact of Noisy Data
% Show how specifying a tolerance for the rank computation in lsqminnorm can help define the scale of the problem so that random noise does not corrupt the solution.
%
% Create a low-rank matrix of rank 5 and a right-hand side vector b.

rng default % for reproducibility
U = randn(200,5);
V = randn(100,5);
A = U*V';
b = U*randn(5,1) + 1e-4*randn(200,1);

%
% Solve the linear system \textrm{Ax}=b using lsqminnorm. Compute the norms of A*x-b and x to check the quality of the solution.

x = lsqminnorm(A,b);
norm(A*x-b)
norm(x)

%
% Now add a small amount of noise to the matrix A and solve the linear system again. The noise affects the solution vector x of the linear system disproportionately.

Anoise = A + 1e-12*randn(200,100);
xnoise = lsqminnorm(Anoise,b);
norm(Anoise*xnoise - b)
norm(xnoise)

%
% The reason for the big difference in the solutions is that the noise affects the low-rank approximation of A. In other words, lsqminnorm is treating small values on the diagonal of the R matrix in the QR decomposition of A as being more important than they are. Ideally, these small values on the diagonal of R should be treated as zeros.
% Plot the diagonal elements of the R matrix in the QR decomposition of Anoise. A large number of the diagonal elements are on the order of 1e-10.

[Q,R,p] = qr(Anoise,0);
semilogy(abs(diag(R)),'o')

%
% The solution to this issue is to increase the tolerance used by lsqminnorm so that a low-rank approximation of Anoise with error less than 1e-8 is used in the calculation. This makes the result much less susceptible to the noise. The solution using a tolerance is very close to the original solution x.

xnoise = lsqminnorm(Anoise, b, 1e-8);
norm(Anoise*xnoise - b)
norm(xnoise)
norm(x - xnoise)

% Copyright 2012 The MathWorks, Inc.