  function [xLS,rhoxLS,r] = LSviaSVD(A,b,tol)
% A is an mxn matrix (m>=n) and b is a column
% Uses the SVD to solve the 
% (nearly) rank-deficient least squares problem min || Ax - b ||.
% tol is a small positive number that "defines" zero. In particular,
% if a floating point number f satisfies |f| <= tol, then it is treated
% as if it were zero. 
% xLS is the minimizing x and rhoxLS = || A*xLS - b ||.
% r is an estimate of A's rank given the tolerance tol.

[m,n] = size(A);
[U,S,V] = svd(A,0);
% This means that A = USV' where U is mxn with orthonormal columns,
% V is nxn with orthonormal columns, and S is nxn and diagonal with
% S(1,1)>=S(2,2)>=...>=S(n,n)>=0
% Look for the first negligble diagonal entry in S
k = 0;
while k<n && S(k+1,k+1) > tol
    k = k+1;
end
% r is an estimate of A's rank where
r = k;
% Note that if [U,S,V] = svd(A), then 
%  || Ax - b ||^2 =  || (U'AV)(V'x) - U'b ||^2 
%                 =  || Sz - d ||^2
% where z = V'x and d = U'b. Since S has the property that 
% S(r+1,r+1),...,S(n,n) are negligible we have
%
%  || Ax - b ||^2 = || S(1:r,:)z - d(1:r)||^2 + || d(r+1:m) ||^2
%
% To obtain the "minimum norm" solution we set z(r+1:n) = 0 and solve 
% the nonsingular diagonal system S(1:r,1:r)z(1:r) = d(1:r).
% We can do all this even if we have the thin svd(A,0) factorization...
xLS = V(:,1:r)*((U(:,1:r)'*b)./diag(S(1:r,1:r)));
rhoxLS = norm(A*xLS - b);