function [x,mu] = constrained_ls(A,b,B)
% Minimize norm(Ax-b) subject to Bx = 0.

[Q,R] = qr(A,0); % The "economy" QR
C = B/R; % C = B*R^{-1}
d = Q'*b; % Project b onto range of A
mu = C'\d; % Minimize norm(C*mu-d)
r = d-C'*mu; % Residual of minimization
x = R\r; % Compute constrained minimizer