How to calculate inverse using cramer's rule in python?
Question: I'm trying to generate the inverse matrix using numpy package in python.
Unfortunately , I'm not getting the answers I expected.

Original matrix:

`([17 17 5] [21 18 21] [2 2 19])`

Inverting the original matrix by Cramer's rule gives:

`([4 9 15] [15 17 6] [24 0 17])`

Apparently using `numpy.linalg.inv()` gives

`-3.19488818e-01,3.80191693e-01,-6.38977636e-03, 3.33333333e-01,
-3.33333333e-01, 2.26123699e-18, -2.84345048e-01, 2.68370607e-01,
5.43130990e-02n`

I expected that multiplying the original matrix and the inverse would have
given an identity matrix but as you can see I give a matrix filled with
floating points.

Where can be the issue?


Answer: I think you may have made a mistake when inverting the matrix by hand.

When I perform the following

    
    
    import numpy as np
    
    a = np.array([[17, 17, 5], [21, 18, 21], [2, 2, 19]], dtype=np.float)
    inv = np.linalg.inv(a)
    
    print np.dot(inv, a)
    

I get

    
    
    array([[  1.00000000e+00,   0.00000000e+00,   1.05471187e-15],
           [  1.11022302e-16,   1.00000000e+00,  -7.21644966e-16],
           [  1.38235777e-17,   5.65818009e-18,   1.00000000e+00]])
    

Which is fine, note that all the off diagonal elements are approximately zero
to machine precision so it looks like numpy is doing an ok job with this one!

Remember that floating point numbers don't work like real numbers and you can
expect small rounding errors to creep into your calculations unless you are
careful.

* * *

If you want to do this exactly have a look at
[sympy](http://www.sympy.org/en/index.html) which will be able to do the
calculation with exact maths (at the expense of it being a little slower).

    
    
    import sympy as sp
    
    a = sp.Matrix([[17, 17,5],[21,18,21],[2,2,19]])
    inv = a.inv()
    
    print inv
    print a * inv
    

yielding the exact inverse

    
    
    Matrix([
    [-100/313,  1/3, -89/313],
    [ 119/313, -1/3,  84/313],
    [  -2/313,    0,  17/313]])
    

which when multiplied by the original matrix gives the exact identity as you
would expect

    
    
    Matrix([
    [1, 0, 0],
    [0, 1, 0],
    [0, 0, 1]])
    

