CODEWITHSUNDEEP
pythonnumpylinear-algebra
Glostas
Glostas

Reputation: 1180

Python: Calculating the inverse of a pseudo inverse matrix

I am trying to calculate the pseudo inverse of a matrix which should be not very difficult. The problem is inverting the matrix.

I am using the following code:

A=numpy.random.random_sample((4,5,))
A=mat(A)
B=pseudoinverse(A)

def pseudoinverse(A):
    helper=A.T*A
    print helper*helper.I
    PI=helper.I*A.T
    return PI`

to test this I included the print line. helper*helper.I should give unity. The output I get from this is:

 [[ 2.    -1.     0.     0.     3.   ]

 [ 0.     2.     0.     0.     3.5  ]

 [ 0.    -0.5    1.125 -1.     2.25 ]

 [ 2.     0.     0.25  -2.     3.   ]

 [ 0.     0.     0.5   -2.     4.   ]]

which is clearly not unity. I don't know what I did wrong and really would like to know.

Upvotes: 1

Views: 3319

Answers (1)

MB-F
MB-F

Reputation: 23637

Your matrix A does not have full column rank. In consequence helper is singular and not invertible (If you print helper.I you will see some very large numbers).

The solution is to compute the right inverse instead of the left inverse:

helper = A * A.T
PI = A.T * helper.I

See Wikipedia for more details.

Unless you are doing this for exercise, you could also use numpy's built in implementation of the pseudeinverse.

edit

>>> numpy.random.seed(42)
>>> a = mat(numpy.random.random_sample((3, 4)))  # smaller matrix for nicer output
>>> h = a * a.T
>>> h * h.I
matrix([[  1.00000000e+00,   1.33226763e-15,   0.00000000e+00],
        [ -1.77635684e-15,   1.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   1.33226763e-15,   1.00000000e+00]])

Up to numeric precision this looks pretty much like an identity matrix to me.

The problem in your code is that A.T * A is not invertible. If you try to invert such a matrix you get wrong results. In contrast, A * A.T is invertible.

You have two options:

  1. change the direction of multiplication
  2. call pseudoinverse(A.T)

Upvotes: 1

Related Questions