# Equality-constrained least-squares when the matrix is singular

I have to solve the linear system $$Ax = b$$ in the least-square sense, where matrix $$A$$ is singular. To resolve this, I introduce $$Cx = d$$. However, even after this, I am not able to solve it using scipy.optimize.minimize as it keeps complaining that the matrix is singular. I am trying to solve it using the following:

optimum_func = lambda x: np.sum((Ax -b)**2)
constraint = ({'type': 'eq', 'func': lambda x: C@x-d})


What is the right way to solve this problem? And what theory did I miss that it does not work?

If it is possible to obtain a unique solution to $$Ax=b$$ subject to the constraints $$Cx =d,$$ then you can obtain it with Lagrange multipliers.

To be explicit, suppose $$A$$ is an $$n\times p$$ matrix, $$x$$ is a $$p$$-vector, $$b$$ is an $$n$$-vector, $$C$$ is an $$r\times p$$ matrix (the coefficients of $$r$$ linear constraints), and $$d$$ is an $$r$$ vector (the values of those constraints). The Lagrange multipliers $$\lambda$$ will be an $$r$$-vector.

The system that results from minimizing $$||Ax-b||^2 = x^\prime A^\prime A x - 2x^\prime A^\prime b + b^\prime b$$ subject to $$Cx-d = 0$$ is $$A^\prime A x + \lambda^\prime C = 0$$ and $$Cx = d,$$ to be solved for $$x$$ and $$\lambda.$$ In block matrix notation this can be written as a set of $$p+r$$ equations in the $$p+r$$ variables,

$$\pmatrix{A^\prime A & C^\prime \\ C & 0}\pmatrix{x \\ \lambda} = \pmatrix{A^\prime b\\d}$$ with solution

$$\pmatrix{x \\ \lambda} = \pmatrix{A^\prime A & C^\prime \\ C & 0}^{-1}\pmatrix{A^\prime b\\d}$$

This system can still be singular. That would mean imposing the constraints did not create an identifiable model.

### Important comment

This solution does not assume $$A^\prime A$$ is singular. Thus, it solves any least-squares problem subject to a set of linear constraints. You might have to understand the solution in the sense of a generalized (Moore-Penrose) inverse, as is always the case with ordinary least squares.

### Simple example

Let

$$A = \pmatrix{1&1\\1&1\\1&1};\quad b = \pmatrix{1\\2\\3};\quad C = \pmatrix{1&2};\quad d = \pmatrix {-1},$$

for which the dimensions are $$n=3,$$ $$p=2,$$ and $$r=1.$$ $$A$$ is obviously rank-deficient because its columns are linearly dependent.

This problem seeks to minimize $$(x_1+x_2-1)^2 + (x_1+x_2-2)^2 + (x_1+x_2-3)^2$$ subject to the constraint $$x_1 + 2x_2 = -1.$$ Writing

$$\pmatrix{A^\prime A & C^\prime \\ C & 0} = \pmatrix{3&3&1\\3&3&2\\1&2&0}\quad \text{and}\quad \pmatrix{A^\prime b\\d} = \pmatrix{6\\6\\-1},$$

the solution is

$$\pmatrix{x_1\\x_2\\\lambda} = \pmatrix{5\\-3\\0}.$$

To check the solution, you can easily check that $$Cx=d$$ holds. Notice that the constraint $$Cx=d$$ implies $$x_1+x_2=-1-x_2.$$ Plugging this into the least squares objective gives

\begin{aligned} (x_1+x_2-1)^2 + (x_1+x_2-2)^2 + (x_1+x_2-3)^2 \\ = (-x_2-1)^2 + (-x_2-2)^2 + (-x_2-3)^2\\ = 3x_2^2 + 18 x_2 + 29 \end{aligned}

whose derivative $$6x_2 + 18$$ has a unique critical point at $$x_2=-3,$$ whence $$x_1 = -1 - 2x_2 = 5,$$ verifying the solution really does minimize the squared distance between $$Ax$$ and $$b.$$

If Cx = d has a unique solution x0, then you can check whether x0 is a solution to Ax = b as well.

More generally, since you have two equality constraints for x, you can stack the matrices A & C and the vectors b & d (vertically) to construct a new system of equations that contains all information you have about the solution.

Let's call this system A_star x = b_star. It's possible that the bigger system is under-determined, over-determined or has a unique solution. It depends on the inputs A, C, b and d.

In any case A_star is not a square matrix: it has more rows than columns because we added the Cx = d constraints to the original system. So you can find the best solution to A_star x = b_star and check whether it is an exact solution to Ax = b.

Here are two approaches: one using the pseudo inverse and the other using least squares.

import numpy as np

A = np.array([[3, 4], [6, 8]])
b = np.array([[15], [30]])

C = np.array([[1, 0]])
d = np.array([[1]])

A_star = np.vstack([A, C])
b_star = np.vstack([b, d])

# We can't use np.linalg.solve because A_star is not square.
# If A.shape = (n, n) and C.shape = (m, n), then A_star.shape = (m + n, n).
# np.linalg.solve(A_star, b_star)

# Compute the pseudo inverse and solve
A_star_pinv = np.linalg.pinv(A_star)
x = A_star_pinv @ b_star
# Check if A*x = b, or how close the solution is
A @ x

# Or find the least squares solution
soln = np.linalg.lstsq(A_star, b_star)
x = soln[0]
# Check if A*x = b
A @ x

• Where in this solution does it handle the constraint $Cx=d$??
– whuber
Mar 7 at 2:05
• I totally missed the part to first stack all equality constraints together and then find the (best possible) solution to the augmented system. Mar 8 at 9:54
• Yours is an interesting approach. But the question appears to view the equations $Ax=b$ differently than the "constraint" $Cx=d.$ By calling the latter a constraint, the OP is saying it has to be satisfied exactly, not just in the least squares sense.
– whuber
Mar 8 at 19:00