Equality-constrained least-squares when the matrix is singular

Question

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?

Answer 1

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

Answer 2

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