function xB = LowRankLS(F,G,b)
% F is m-by-r, G is n-by-r, and b is m-by-1.
% Assume rank(F) = rank(G) = r < n.
% xB is the basic solution to min norm(F*G'*x - b)

[m,r] = size(F);
[n,r] = size(G);
%
% Suppose F = QR and Q = [Q1 Q2] where Q1 has r columns. It follows that
% norm(FG'x - b,2)^2 = norm(RG'x - Q'b)^2
% = norm(R(1:r,1:r)G'x - Q1'^b)^2 + norm(Q2'*b)^2
% so we need the basic solution to
% min norm(R(1:r,1:r)G'x - Q1'^b)

[Q1,R1] = qr(F,0);
c = Q1'*b;
[QG,RG,pVec] = qr(R1*G',0);
xB = zeros(n,1);
xB(pVec) = [RG(1:r,1:r)\(QG'*c);zeros(n-r,1)];

% This also gives a basic solution, but not necessarily the
% same one as above...
[Q1,R1] = qr(F,0);
z = R1\(Q1'*b); % this z minimizes norm(Fz-b)
% Now get a basic solution to G'x = z G'x = QRP'x = b
[QG,RG,pVec] = qr(G',0);
xB0 = zeros(n,1);
xB0(pVec) = [RG(1:r,1:r)\(QG'*z);zeros(n-r,1)];

ERROR = norm(xB0-xB)