Computing pseudo inverse of a matrix using sympy

Question

How should I compute the pseudo-inverse of a matrix using sympy (not using numpy, because the matrix has symbolic constants and I want the inverse also in symbolic). The normal inv() does not work for a non-square matrix in sympy. For example if M = Matrix(2,3, [1,2,3,4,5,6]), pinv(M) should give

-0.9444   0.4444
-0.1111   0.1111
 0.7222  -0.2222

Answer 1

I think since this is all symbolic it should be OK to use the text-book formulas taught in a linear algebra class (e.g. see the list of special cases in the Wikipedia article on the Moore–Penrose pseudoinverse). For numerical evaluation pinv uses the singular value decomposition (svd) instead.

You have linearly independent rows (full row rank), so you can use the formula for a 'right' inverse:

>>> import sympy as sy
>>> M = sy.Matrix(2,3, [1,2,3,4,5,6])

>>> N = M.H * (M * M.H) ** -1

>>> N.evalf(4)
[-0.9444,  0.4444]
[-0.1111,  0.1111]
[ 0.7222, -0.2222]
>>> M * N
[1, 0]
[0, 1]

For full column rank, replace M with M.H, transpose the result, and simplify to get the following formula for the 'left' inverse:

>>> M = sy.Matrix(3, 2, [1,2,3,4,5,6])

>>> N = (M.H * M) ** -1 * M.H

>>> N.evalf(4)
[-1.333, -0.3333,  0.6667]
[ 1.083,  0.3333, -0.4167]
>>> N * M
[1, 0]
[0, 1]