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?

  • $\begingroup$ Can't you reduce to regular optimization with equality constraints ? $\endgroup$ – Claude Leibovici Jan 2 '16 at 8:27
  • $\begingroup$ Would a gradient based method be OK? Though there is also a direct solver. $\endgroup$ – Royi Aug 27 '17 at 16:35

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

  • $\begingroup$ a quick and non dirtay way is to notice that when $\lambda \to \infty$ minimizing $||Ax-b||^2+\lambda||Bx-d||^2$ is equivalent to minimize $||A ( B^+ d + B^{\perp} y) - b||^2$ and setting $x^* = ( B^+ d + B^{\perp} y^*)$ yields the desired solution $\endgroup$ – reuns Jan 21 '16 at 10:12
  • $\begingroup$ This is not obvious to me... Could you please explain (or post an answer, I'll be the first to upvote)? $\endgroup$ – Ash Jan 21 '16 at 18:02
  • $\begingroup$ when $\lambda \to \infty$ : $Bx^* \to d$ $\endgroup$ – reuns Jan 21 '16 at 19:08

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);

enter image description here

The full code, including validation using CVX, can be found in my StackExchange Mathematics Q1596362 GitHub Repository.


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.


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.