% Generate random 6x5 rank-deficient 
% matrix A for which columns 3 and 4 
% are almost dependent on columns 1 and 2.
%
% Then solve least squares problems with
% "true solutions" that have all components
% equal to 1 but include random errors in
% the data vector.
%
% Illustrate the fact that the solution vector
% is not well-determined, but that almost
% as small a residual can be obtained using
% 3 columns of A rather than all 5, and this
% reduction gives a solution vector that is
% not affected much by the random errors.
%
% Class demonstration 09/2008
% Dianne P. O'Leary

a = rand(6,3);
t= a(:,1:2)*rand(2,2)+rand(6,2)*1.e-3;
A = [a(:,1),a(:,2),t(:,1),t(:,2),a(:,3)];
[Qp,Rp,Pp] = qr(A)
[Q,R] = qr(A)

% Repeatedly generate new rhs data.

for k=1:100,

% Generate a right-hand side vector b
% which is almost in the range of A.

   b = A*ones(5,1) + rand(6,1)*1.e-3;

% Compute the least squares solution x
% to minimize r = b - A x over all
% vectors in the subspace spanned by
% the linearly independent columns of A.

   c = Qp'*b;
   y = Rp(1:3,1:3) \ c(1:3);
   
   xp = Pp*[y;0;0];

% Compute the least squares solution x
% to minimize r = b - A x.

   x=A\b;

% Compare the residuals.

   rp = b - A*xp;
   r = b - A*x;
   disp('The two solution vectors:')
   [xp,x]
   disp('The two residual norms:')
   [ norm(rp), norm(r)]
   disp('pausing')
   pause

end