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?


2 Answers 2


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


$$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
  • $\begingroup$ Where in this solution does it handle the constraint $Cx=d$?? $\endgroup$
    – whuber
    Mar 7 at 2:05
  • $\begingroup$ I totally missed the part to first stack all equality constraints together and then find the (best possible) solution to the augmented system. $\endgroup$
    – dipetkov
    Mar 8 at 9:54
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber
    Mar 8 at 19:00

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.