function x = UnderD(A,b,xStar)
% A is an m-by-n matrix with rank(A) = m < n.
% b is m-by-1 and xStar is n-by-1.
% x minimizes norm(x - xStar,2) subject to the constraint that Ax = b

[m,n] = size(A);
[U,S,V] = svd(A);
btilde = U'*b;
y = btilde./diag(S);
x0 = V(:,1:m)*y;
% Note that A*x0 = b.
% Determine z so || (x0 + V(:,m+1:n)*z) - xStar || = min.
% Thus z solves the least square problem min || Cz - d || where
% C = V(:,m+1:n) and d = xStar - x0. Observe that C has orthonormal
% columns, so by normal equations, z = V(:,m+1:n)'*d
z = V(:,m+1:n)'*(xStar - x0);
x = x0 + V(:,m+1:n)*z;