%
%
% Solve rank deficient linear least squaeres problems.
%

example = 2;

if( example == 1 )
    % Generate problem data
    % Compute A
    t = 10.^(0:-1:-10)';
    A = [ ones(size(t))  t   t.^2  t.^3];

    % compute exact data
    xex   = ones(4,1);
    b     = A*xex;

    % try to fit   x1 + x2*t + x3*t^2 + x4*t^3 + x5*(1-t)  to data

    A = [ ones(size(t))  t   t.^2  t.^3 (1-t)];
    
else % example == 2

    % Load the Filip data set
    % (see http://www.itl.nist.gov/div898/strd/lls/data/Filip.shtml)
    load Filip.dat
    b = Filip(:,1);
    t = Filip(:,2);

    % try to fit   x1 + x2*t + x3*t^2 + ... + x11*t^10  to data
    A = ones(size(t,1),11);
    for i = 2:11
        A(:,i) = A(:,i-1).*t;
    end
end

[m,n] = size(A);




% Solve using QR decomposition:
[Q,R,P] = qr(A);
figure(2)
semilogy(abs(diag(R)),'*')
title('absolute values of diagonals of R')

% determine the effective rank
r = 1;
while( r < size(A,2) & abs(R(r+1,r+1)) >= max(size(A))*eps*abs(R(1,1)) )
       r = r+1;
end
disp([' The effective rank of A is ',int2str(r)])


% compute minimum norm solution
d  = Q'*b;

% z2 solves
%   min || [ R(1:r,1:r)\R(1:r,r+1:n) ]        [ R(1:r,1:r) \ d(1:r) ] ||
%       || [          I              ] * z2 - [      0              ] ||
      
     
B  = [ R(1:r,1:r)\R(1:r,r+1:n); eye( n-r ) ];
z2 = B \ [R(1:r,1:r) \ d(1:r);
           zeros(n-r,1) ];
          
z  = [- R(1:r,1:r) \ (R(1:r,r+1:n)*z2 - d(1:r));
              z2  ];
x  = P*z;
disp('Minimum norm solution computed using the QR decomposition')
disp(x)
disp([' || x ||_2 = ',num2str(norm(x))])




disp('Hit RETURN to continue....'); pause

% compute a particular solutiuon
d  = Q'*b;

% set z2 = 0
          
z  = [ R(1:r,1:r) \ d(1:r);
          zeros(n-r,1)  ];
x  = P*z;
disp('Least squares solution computed using the QR decomposition')
disp(x)
disp([' || x ||_2 = ',num2str(norm(x))])


disp('Least squares solution computed using A\b')
x = A\b;
disp(x)
disp([' || x ||_2 = ',num2str(norm(x))])




disp('Hit RETURN to continue....'); pause

% Solve using SVD:
[U,S,V] = svd(A);
s       = diag(S);
figure(2)
semilogy(s,'*')
title('singular values of A')
% determine the effective rank using singular values
r = 1;
while( r < size(A,2) & s(r+1) >= max(size(A))*eps*s(1) )
       r = r+1;
end
disp([' The effective rank of A is ',int2str(r)])

d  = U'*b;
x  = V* ( [d(1:r)./s(1:r);
            zeros(n-r,1) ] );
disp('Minimum norm solution computed using using the SVD')
disp(x)
disp([' || x ||_2 = ',num2str(norm(x))])