Why do Mathematica and Python's answers differ when dealing with singular matrix equations?

Question

I have been dealing with linear algebra problems of the form A = Bx in Python and comparing this to a colleague's code in MATLAB and Mathematica. We have noticed differences between Python and the others when B is a singular matrix. When using numpy.linalg.solve() I throw a singular matrix error, so I've instead implemented .pinv() (the Moore Penrose pseudo inverse).

I understand that storing the inverse is computationally inefficient and am first of all curious if there's a better way of dealing with singular matrices in Python. However the thrust of my question lies in how Python chooses an answer from an infinite solution space, and why it chooses a different one than MATLAB and Mathematica do.

Here is my toy problem:

B = np.array([[2,4,6],[1,0,3],[0,7,0]])
A = np.array([[12],[4],[7]])
BI = linalg.pinv(B)
x = BI.dot(A)

The answer that Python outputs to me is:

[[ 0.4]
[ 1. ]
[ 1.2]]

While this is certainly a correct answer, it isn't the one I had intended: (1,1,1). Why does Python generate this particular solution? Is there a way to return the space of solutions rather than one possible solution? My colleague's code returned (1, 1, 1) - is there a reason that Python is different from Mathematica and MATLAB?

Why don't you print out BI instead for a more useful comparison? — Eric
"My colleague's code returned (1,1,1)" - was this the matlab, mathematica, or both? Can you show that code? — Eric
pinv(B) * A gives me the same result in matlab. So what is your matlab code? — Eric
np.linalg.lstsq is the singular version of solve, which gives the same asnwer as what you have. — Eric

jme jme · Accepted Answer · 2016-10-25T00:55:35

In short, your code (and apparently np.linalg.lstsq) uses the Moore-Penrose pseudoinverse, which is implemented in np.linalg.pinv. MATLAB and Mathematica likely use Gaussian elimination to solve the system. We can replicate this latter approach in Python using the LU decomposition:

B = np.array([[2,4,6],[1,0,3],[0,7,0]])
y = np.array([[12],[4],[7]])
P, L, U = scipy.linalg.lu(B)

This decomposes B as B = P L U, where U is now an upper-diagonal matrix, and P L is invertible. In particular, we find:

>>> U
array([[ 2.,  4.,  6.],
       [ 0.,  7.,  0.],
       [ 0.,  0.,  0.]])

and

>>> np.linalg.inv(P @ L) @ y
array([[ 12.],
       [  7.],
       [  0.]])

The goal is to solve this under-determined, transformed problem, U x = (P L)^{-1} y. The solution set is the same as our original problem. Let a solution be written as x = (x_1, x_2, x_3). Then we immediately see that any solution must have x_2 = 1. Then we must have 2 x_1 + 4 + 6 x_2 = 12. Solving for x_1, we get x_1 = 4 - 3 x_2. And so any solution is of the form (4 - 3 x_2, 1, x_2).

The easiest way to generate a solution for the above is to simply choose x_2 = 1. Then x_1 = 1, and you recover the solution that MATLAB gives you: (1, 1, 1).

On the other hand, np.linalg.pinv computes the Moore-Penrose pseudoinverse, which is the unique matrix satisfying the pseudionverse properties for B. The emphasis here is on unique. Therefore, when you say:

my question lies in how Python chooses an answer from an infinite solution space

the answer is that all of the choosing is actually done by you when you use the pseudoinverse, because np.linalg.pinv(B) is a unique matrix, and hence np.linalg.pinv(B) @ y is unique.

To generate the full set of solutions, see the comment above by @ali_m.

Why do Mathematica and Python's answers differ when dealing with singular matrix equations?

1 Answers