function [x_ls, residual] = rdlstsq(A,b)
%rdlstsq Compute the minimum norm solution to the
%linear least-squares problem Ax = b.
%
%   x_ls = rdlstsq(A,b) uses the SVD to find the minimum norm solution
%   to the least-squares problem Ax = b, where A is a real m x n matrix,
%   and b is an m x 1 vector, m >= n. Use when A does not have full rank.

[~, n] = size(A);

[U, S, V] = svd(A);
sigma = diag(S);

% compute rank(A)
tol = max(size(A))*eps(max(sigma));
r = sum(sigma > tol);

x_ls = zeros(n,1);
for i = 1:r
   x_ls = x_ls + ((U(:,i)'*b)/sigma(i))*V(:,i);
end

residual = 0;
for i = r+1:size(A,1)
	residual = residual + (b'*U(:,i))^2;
end
residual = sqrt(residual);
