Linear Least Squares with Linear Equality Constraints - Iterative Solver I am looking for iterative procedures for solution of the linear least squares problems with linear equality constraints.
The Problem:
$$ \arg \min_{x} \frac{1}{2} \left\| A x - b \right\|_{2}^{2}, \quad \text{subject to} \quad B x = d $$
How can best the two systems can be combined so that iterative procedures can be applied on it?
 A: The easiest and most straight forward iterative method would be the Projected Gradient Descent.
Projection onto Linear Equality Equation
The Projected Gradient Descent requires a projection step onto the Linear Equation constraint.
The problem is given by:
$$ \arg \min_{x} \frac{1}{2} \left\| x - y \right\|_{2}^{2}, \quad \text{subject to} \quad B x = d $$
From the KKT Conditions one could see the result is given by:
$$ \hat{x} = y - {B}^{T} \left( B {B}^{T} \right)^{-1} \left( B y - d \right) $$
Full derivation is given at Projection of $ z $ onto the Affine Set $ \left\{ x \mid A x = b \right\} $.
Now, all needed is to integrate that into the Projected Gradient Descent framework.
%% Solution by Projected Gradient Descent

hObjFun     = @(vX) (0.5 * sum((mA * vX - vB) .^ 2));
hProjFun    = @(vY) vY - (mB.' * ((mB * mB.') \ (mB * vY - vD)));
vObjVal = zeros([numIterations, 1]);

mAA = mA.' * mA;
vAb = mA.' * vB;

vX          = mB \ vD; %<! Initialization by the Least Squares Solution of the Constraint
vX          = hProjFun(vX);
vObjVal(1)  = hObjFun(vX);

for ii = 2:numIterations

    vX = vX - (stepSize * ((mAA * vX) - vAb));
    vX = hProjFun(vX);

    vObjVal(ii) = hObjFun(vX);
end


The full code, including validation using CVX, can be found in my StackExchange Mathematics Q1596362 GitHub Repository.
A: Let $A\in M_{m,n}$ with $rank(A)=n\leq m$, that is, $A$ is one to one. Let $B\in M_{k,n}$ with $0<r=dim(\ker(B))<m$; here $x\in \mathbb{R}^n,d\in\mathbb{R}^k$. Thus $x\in \Pi$, an affine subspace of dimension $r$ of $\mathbb{R}^n$ and $Ax\in A(\Pi)$, an affine subspace of $\mathbb{R}^m$ of dimension $r$. If $x_0$ is the required solution, then $Ax_0$ is the orthogonal projection of $b$ on $A(\Pi)$.
Method 1. We seek a basis $(e_1,\cdots,e_r)$ of $\ker(B)$; then $(Ae_1,\cdots,Ae_r)$ is a basis of $A(\ker(B))$; we deduce an orthonormal basis (by Schmidt process) of $A(ker(B))$ ....
Method 2. Use the Lagrange method. The unknowns are $x\in \mathbb{R}^n$ and $\Lambda\in M_{1,k}$ and the $n+k$ linear equations are $Bx=d,2(Ax-b)^TA+\Lambda B=0$. 
A: A quick and dirty way that I have seen work well in practice is to minimize $$||Ax-b||+\lambda||Bx-d||$$ using iterative methods such as Levenberg-Marquardt. In your case, I would set a higher value for $\lambda$. In case finding a correct $\lambda$ is difficult (because of scale differences, etc), you can minimize $$\frac{1}{s-||Bx-d||}+||Ax-b||$$ where s is a threshold that you do not want ||Bx-d|| to exceed (in that case, you might want to initialize your minimization algorithm with the $x$ that minimizes $||Bx-d||$).   
A: I would propose reweighted linear least squares.
This way you can get away with only a simple least squares solver and very little overhead in updating the weight matrix.
$$+w\|Bx-d\|_2^2$$
Where the weight $w$ is updated to enforce the equality more and more in between iterations. You can even try a very big $w$ at the start. Maybe you don't even need much of reweighting.
