Matrix inverse in numpy/python not giving correct matrix?

Question

I have a nxn matrix C and use inv from numpy.linalg to take the inverse to get Cinverse. My Cmatrix has elements of order 10**4 but my Cinverse matrix has elements of order 10**12 and higher (not sure if thats correct). When I do numpyp.dot(C,Cinverse), I do not get the identity matrix. Why is this?

I have a vector x which I multiply by itself to get a matrix.

x=array([ 121.41191662,   74.22830468,   73.23156336,   75.48354975,
     79.89580817])
c=np.outer(xvector,xvector)

this is a 5x5 matrix.

then I get its inverse by

from numpy.linalg import inv
cinverse=inv(c)

then I want to see if I can get identity matrix back.

identity=np.dot(C00,C00inv)

However, I do not get the identity matrix. cinverse has very large matrix elements around 10**13 and higher while c has matrix elements around 10,000.

Please be more specific. Can you show us your input and output? Also, what exactly is the code you are using? — juanpa.arrivillaga
– juanpa.arrivillaga, Commented Feb 14, 2017 at 20:06
Could you provide some example code with random matrices to show what you are trying to do? That doesn't sound right- when you say has elements of order 10**4, does that mean that number of elements of that is the size of the number, i.e. 10000ish? Ah, juanpa, beat me to it ;) — RexFuzzle
– RexFuzzle, Commented Feb 14, 2017 at 20:06
"I do not get the identity matrix" - how close is it to the identity? You should not expect a matrix of exactly 1s and 0s because of rounding error. — ali_m
– ali_m, Commented Feb 14, 2017 at 21:08
To be more specific, my matrix elements are of order 10**4. For example, 8000,10000,7000.... etc. I edited original post with some code. — Nabeel Eh
– Nabeel Eh, Commented Feb 14, 2017 at 23:48
@RexFuzzle I have edited the original post for clarity. Yes it is for example a 5x5 matrix with the magnitude of the elements being around 10000ish. — Nabeel Eh
– Nabeel Eh, Commented Feb 15, 2017 at 0:07

Paul Panzer · Accepted Answer · 2017-02-15 04:25:12Z

1

The outer product of two vectors (be they the same or not) is not invertible. Since it is just a stack of scaled copies of the same vector its rank is one. Rank defective matrices cannot be inverted.

I'm surprised that numpy is not raising an exception or at least giving a warning.

answered Feb 15, 2017 at 4:25

Paul Panzer

53.3k3 gold badges60 silver badges103 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

Nabeel Eh Over a year ago

Why is this? If you look at eqn 14 of academic.oup.com/ptep/article/doi/10.1093/ptep/ptu065/1560725/… you'll see that I need an outer product of two vectors which creates a covariance matrix (I am just considering l=m=0 case). Eqn 13 requires me to take the inverse

Paul Panzer Over a year ago

@NabeelEh But you do get the mathematical argument, I trust? Whatever that paper says. If you have more than one dimension, the outer product is not invertible. That is a fact. Now I really cannot be expected to read your papers for you, but this once: What you have there (14) is not an outer product. It is a sum of outer products. Obviously, if you add enough rank 1 matrices you can get full rank.

Nabeel Eh Over a year ago

I see. Yes, I get the argument. Thanks Paul, appreciate the help!

RexFuzzle · Accepted Answer · 2017-02-15 03:53:36Z

0

So here is some code that generates the inverse matrix, and I will comment about it afterwards.

import numpy as np
x = np.random.rand(5,5)*10000 # makes a 5x5 matrix with elements around 10000
xin = np.linalg.inv(x)
iden = np.dot(x,xinv)

Now the first line of your iden matrix probably looks something like this: [ 1.00000000e+00, -2.05382445e-16, -5.61067365e-16, 1.99719718e-15, -2.12322957e-16] . Notice that the first element is exactly 1, as it should be, but there others are not exactly 0, however they are essentially zero and should be regarded as zero according to machine precision.

answered Feb 15, 2017 at 3:53

RexFuzzle

1,4522 gold badges18 silver badges30 bronze badges

1 Comment

Nabeel Eh Over a year ago

Hi @RexFuzzle, I would be happy if those were my results. My vector is x=array([ 4140.56172552, 2561.7498602 , 2529.56826173, 2607.35090233, 2759.5614813 ]). I get my matrix by using c=np.outer(x,x) and then inverse using cinv=inv(c). Then using np.dot(c,cinv) my first line is [ 6.01898362e-01, 5.53081985e+00, 1.79258876e+00, -5.86243807e-01, 1.50883747e+01] which is not identity

Collectives™ on Stack Overflow

Matrix inverse in numpy/python not giving correct matrix?

2 Answers 2

3 Comments

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related