% 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