Left inverse in numpy or scipy?

Question

I am trying to obtain the left inverse of a non-square matrix in python using either numpy or scipy. How can I translate the following Matlab code to Python?

>> A = [0,1; 0,1; 1,0]

A =

     0     1
     0     1
     1     0

>> y = [2;2;1]

y =

     2
     2
     1

>> A\y

ans =

    1.0000
    2.0000

Is there a numpy or scipy equivalent of the left inverse \ operator in Matlab?

Answer 1

For those who wish to solve large sparse least squares problems:

I have added the LSQR algorithm to SciPy. With the next release, you'll be able to do:

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import lsqr
import numpy as np

A = csr_matrix([[0., 1], [0, 1], [1, 0]])
b = np.array([[2.], [2.], [1.]])

lsqr(A, b)

which returns the answer [1, 2].

If you'd like to use this new functionality without upgrading SciPy, you may download lsqr.py from the code repository at

http://projects.scipy.org/scipy/browser/trunk/scipy/sparse/linalg/isolve/lsqr.py

Answer 2

You can calculate the left inverse using matrix calculations:

import numpy as np

linv_A = np.linalg.solve(A.T.dot(A), A.T)

(Why? Because:

)

Test:

np.set_printoptions(suppress=True, precision=3)
np.random.seed(123)

A = np.random.randn(3, 2)
print('A\n', A)

A_linv = np.linalg.solve(A.T.dot(A), A.T)
print('A_linv.dot(A)\n', A_linv.dot(A))

Result:

A
 [[-1.086  0.997]
 [ 0.283 -1.506]
 [-0.579  1.651]]
A_linv.dot(A)
 [[ 1. -0.]
 [ 0.  1.]]

Answer 3

You can use lsqr from scipy.sparse.linalg to solve sparse matrix systems with least squares

Answer 4

Use linalg.lstsq(A,y) since A is not square. See here for details. You can use linalg.solve(A,y) if A is square, but not in your case.

Answer 5

You can also look for the equivalent of the pseudo-inverse function pinv in numpy/scipy, as an alternative to the other answers that is.

Answer 6

Here is a method that will work with sparse matrices (which from your comments is what you want) which uses the leastsq function from the optimize package

from numpy import *
from scipy.sparse import csr_matrix
from scipy.optimize import leastsq
from numpy.random import rand

A=csr_matrix([[0.,1.],[0.,1.],[1.,0.]])
b=array([[2.],[2.],[1.]])

def myfunc(x):
    x.shape = (2,1)
    return (A*x - b)[:,0]

print leastsq(myfunc,rand(2))[0]

generates

[ 1.  2.]

It is kind of ugly because of how I had to get the shapes to match up according to what leastsq wanted. Maybe someone else knows how to make this a little more tidy.

I have also tried to get something to work with the functions in scipy.sparse.linalg by using the LinearOperators, but to no avail. The problem is that all of those functions are made to handle square functions only. If anyone finds a way to do it that way, I would like to know as well.