Linear Least Squares with Linear Equality Constraints - Iterative Solver

Question

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?

Answer 1

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$.

Answer 2

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||$).

Answer 3

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.