# Algorithm for solving sparse equality-constrained least squares

I have a diagonal, positive-definite inner product matrix $M$ and want to find a minimizer of

$$\min_q \frac{1}{2} \|q-q_0\|_M^2\qquad \text{s.t.}\qquad C^Tq+c_0 = 0,$$

where $q_0, c_0$, and $C$ are given and $C$ is very large and sparse.

What is the best algorithm for efficiently solving this problem? The straightforward approach is to solve the normal equations $$\left[\begin{array}{cc}M & C\\C^T & 0\end{array}\right]\left[\begin{array}{c}q\\ \lambda\end{array}\right] = \left[\begin{array}{c} Mq_0 \\-c_0\end{array}\right]$$ for $q$; the matrix here is sparse and symmetric, but indefinite and possibly singular, so sparse iterative solvers like good old Conjugate Gradient cannot be used.

A little manipulation of the above gives an alternative two-step approach where you solve first for the Lagrange multipliers: $$C^T M^{-1} C\lambda = C^Tq_0 + c_0$$ and then substitute $\lambda$ into $$q = q_0 - M^{-1} C\lambda.$$ Since $M$ is diagonal it is trivial to invert, and $C^T M^{-1} C$ is symmetric positive-(semi)definite, but is dense and so also is inefficient to solve when the number of constraints is very large.

What is the standard approach to solving this problem?

• "the matrix here is sparse and symmetric, but indefinite and possibly singular, so sparse iterative solvers like good old Conjugate Gradient cannot be used" - right, that's where things like MINRES(-QLP) come in... – J. M. isn't a mathematician Oct 18 '11 at 14:56
• Thanks, that looks promising – user7530 Oct 18 '11 at 15:50

J.M. is correct that the difficulty is that C does not necessarily have full column rank, causing your augmented system to be singular. The standard approach in this case would be to regularize your problem, i.e., replace the constraint $C^T q + c_0 = 0$ with $C^T q + \delta r + c_0 = 0$ for some small value of $\delta > 0$ (e.g., between $10^{-6}$ and $10^{-8}$) and a new variable $r$ representing the residual of your constraints. Your constraint matrix is now $$\begin{bmatrix} C^T & \delta I \end{bmatrix}$$ which always has full row rank. At the same time, you'll want to add a term $\delta^2 \|r\|^2$ to your objective so as to minimize this residual.
There is an alternative if you have a means to identify efficiently a particular solution $\bar{q}_0$ of $C^T q + c_0 = 0$ (by factorization or using some special structure of the problem). You can "shift" your linear system by defining $\bar{q} := q - \bar{q}_0$: $$\begin{bmatrix} M & C \\ C^T & 0 \end{bmatrix} \begin{bmatrix} \bar{q} \\ \lambda \end{bmatrix} = \begin{bmatrix} M(q_0 - \bar{q}_0) \\ 0 \end{bmatrix}.$$ These are the optimality conditions of the (unconstrained) linear least-squares problem $$\min_{\lambda} \ \tfrac{1}{2} \| C \lambda - M(q_0 - \bar{q}_0) \|_{M^{-1}}^2.$$ NOW you can use the conjugate gradients, or better yet, LSQR (still from Michael Saunders' page). LSQR is a variant of CG which will be more accurate if $C$ is ill conditioned. Once you have found $\lambda$, you can recover $\bar{q} = q_0 - \bar{q}_0 - M^{-1} C \lambda$ and your solution is $q = \bar{q} + \bar{q}_0$.
For more general $M$ (non diagonal), there is a nice hybrid alternative. It consists in interpreting the shifted linear systems as the optimality conditions of the quadratic program $$\min_{\bar{q}} \ -(q_0 - \bar{q}_0)^T M \bar{q} + \tfrac{1}{2} \bar{q}^T M \bar{q} \quad \text{s.t.} \ C^T \bar{q} = 0.$$ The method is a hybrid between a factorization and CG. It needs to factorize the matrix $$\begin{bmatrix} I & C \\ C^T & 0 \end{bmatrix}$$ and applies CG to $M$. At each iteration, the search direction is projected into the nullspace of $C^T$, ensuring that all iterates satisfy $C^T \bar{q} = 0$. The main reference is http://epubs.siam.org/sisc/resource/1/sjoce3/v23/i4/p1376_s1
Clearly, for diagonal $M$, there's no advantage versus a direct factorization of your original system. But in many applications, replacing $M$ by $I$ leads to a very sparse matrix. Also $M$ need not be positive definite any more but only positive definite on the nullspace of $C^T$.