function P1()
% Tests svdLS

clc
randn('seed',0)
% Example 1.
disp('Example 1 ')
disp('                     ')
m = 100; n = 10; r = 2;
F = randn(m,r); G = randn(n,r); b = randn(m,1);
% Explicit SVD solution...
[U,S,V] = svd(F*G');
x0 = V(:,1:r)*((U(:,1:r)'*b)./diag(S(1:r,1:r)));
% compare with this...
x = svdLS(F,G,b);
disp('   Explicit     You')
disp('--------------------')
disp([x x0])

% Example 2
m = 10000; n = 5000; r = 10;
F = randn(m,r); G = randn(n,r); b = randn(m,1);
x0 = svdLSCVL(F,G,b);
x = svdLS(F,G,b);
disp('Example 2')
disp(' ')
fprintf('Eg2Error = %10.3e\n\n',norm(x - x0))


function x = svdLSCVL(F,G,b)
% F is m-by-r with rank r.
% G is n-by-r with rank r.
% b is m-by-1;
% x is the minimum norm solution to min norm((F*G')*x - b,2)

% Suppose F = QR  and G = ZT are full QR factorizations. It follows that
% norm((FG')x - b) = = norm( Q(R*T')Z'x - b) = norm(R*T'y - c) 
% where c = Q'*b and x = Zy. Note that
% R*T' is zero except for the upper left r-by-r submatrix and that 
% submatrix is given by R1*T1' where R1 = R(1:r,1:r) and T1 = T(1:r,1:r).
% Thus, the minimum norm solution of the transformed problem is 
% given by y(1:r) = (R*T')\c(1:r)  and y(r+1:n) = 0. 
% Thus, x = Zy = Z(:,1:r)*y(1:r) and we can
% live with thin QR factorizations...

[Q1,R1] = qr(F,0);
[Z1,T1] = qr(G,0);
ctilde = Q1'*b;                     % This is c(1:r)
[U,S,V] = svd(R1*T1');
ytilde = V*((U'*ctilde)./diag(S));  % This is y(1:r)
x = Z1*ytilde;

% Note that it is not necessary to compute Q1*U and Z1*V. Can get by with
% cheaper mat-vec prods. Didn't deduct points for this oversight.

% -2 if your solution has an m-by-m array, an n-by-n array, and/or
% an m-by-n array. This is an order of mag more than you need and the
% underlying amount of work is also an order of mag more than necessary.
% More points off for not understanding what the minimum norm solution is.



