5

I'm using sympy to find a matrix's inverse. I've the next problem. When I compute the inverse of matrix A and I want prove it, I got a matrix with fractions; I mean

>> import sympy
>> from sympy import pprint
>> from sympy.abc import *
>> import sys
>> sys.displayhook = pprint
>> from sympy.matrices import *
>> A = Matrix([[a, b],[c, d]])
>> B = A.inv()
>> B
>> [1       b*c           -b     ]
>> [- + ------------  -----------]
>> [a    2 /    b*c\    /    b*c\]
>> [    a *|d - ---|  a*|d - ---|]
>> [       \     a /    \     a /]
>> [                             ]
>> [      -c               1     ]
>> [  -----------       -------  ]
>> [    /    b*c\           b*c  ]
>> [  a*|d - ---|       d - ---  ]
>> [    \     a /            a   ]
>> B*A
>> [  /1       b*c     \       b*c        /1       b*c     \       b*d    ]
>> [a*|- + ------------| - -----------  b*|- + ------------| - -----------]
>> [  |a    2 /    b*c\|     /    b*c\    |a    2 /    b*c\|     /    b*c\]
>> [  |    a *|d - ---||   a*|d - ---|    |    a *|d - ---||   a*|d - ---|]
>> [  \       \     a //     \     a /    \       \     a //     \     a /]
>> [                                                                      ]
>> [                                             d          b*c           ]
>> [                0                         ------- - -----------       ]
>> [                                              b*c     /    b*c\       ]
>> [                                          d - ---   a*|d - ---|       ]
>> [                                               a      \     a /       ]

And I wanna get the next matrix

>> I = Matrix([
>> [1, 0],
>> [0, 1]])

My problem is the matrix A*B or B*A. Really I want to simplify the matrix A*B to get I. I tried simplify() but doesn't work.

Yakz
  • 81
  • 1
  • 4
  • what does it mean simplifying? is it some technical term or what? my second question - what are you going to prove? – marmeladze Oct 02 '15 at 22:57
  • 1
    Ok, B is the inverse of A and their product sould be the matrix identity I. In the product of the code I obtein a monster matrix A*B (or B*A) but I want the entrates of this matrix be zero (0). sympy doesn't simplify the entrates (1, 2), (1,1) and (2, 1). You can see the result of this product. Excuse my english, I don't speak it so good. Tranks for your time. – Yakz Oct 03 '15 at 03:17
  • I get the same result as you for A*B, but simplify(A*B) gives the identity matrix. I use notebook with the latest Anaconda 64 on windows 7. – brian Oct 03 '15 at 13:55

2 Answers2

6

You can apply the simplify function to each cell of the matrix with applyfunc, like this:

>>> (B*A).applyfunc(simplify)
[1  0]
[    ]
[0  1]
Adrien
  • 469
  • 3
  • 11
0

Forget python and sympy for a minute. Focus on finding an inverse of matrix with paper and pen.

For a A = [[a, b], [c,d]] matrix, we calculate inverse A^-1 as,

(1/D)*[[d, -b],[-c, a]]. Here D is determinant of A matrix (1/ad-bc)

This (A^-1) is equal to [[d/D, -b/D][-c/D, a/D]]

Let's take the first element from first row and follow the operations I've made. For me they actually make no sense, but this is the way how sympy does :) Then apply this procedure to other Matrix elements. 

=> d/D 
d/(a*d-b*c) 
a*d/(d*a^2 - a*b*c)
(a*d-b*c+b*c)/a^2*(d-b*c/a) 
(a*d - a*b*c/a + b*c)/a^2*(d-b*c/a)
(a*(d-b*c/a) + b*c)/a^2*(d-b*c/a)
a*(d-b*c/a)/a^2*(d-b*c/a) + b*c/a^2*(d-b*c/a)
1/a + b*c/a^2*(d-b*c/a) [this is how sympy outputs]


>>> A = Matrix([[a,b],[c,d]])
>>> B = A**-1 #same as B = A.inv()
>>> B[0]
1/a + b*c/(a**2*(d - b*c/a))

Now, let's have a look up what is sympy A*B output. 

>>> N = A*B
>>> N
Matrix([
[a*(1/a + b*c/(a**2*(d - b*c/a))) - b*c/(a*(d - b*c/a)),                                   0],
[c*(1/a + b*c/(a**2*(d - b*c/a))) - c*d/(a*(d - b*c/a)), d/(d - b*c/a) - b*c/(a*(d - b*c/a))]])
>>> pprint(N)
⎡  ⎛1       b⋅c     ⎞       b⋅c                           ⎤
⎢a⋅⎜─ + ────────────⎟ - ───────────            0          ⎥
⎢  ⎜a    2 ⎛    b⋅c⎞⎟     ⎛    b⋅c⎞                       ⎥
⎢  ⎜    a ⋅⎜d - ───⎟⎟   a⋅⎜d - ───⎟                       ⎥
⎢  ⎝       ⎝     a ⎠⎠     ⎝     a ⎠                       ⎥
⎢                                                         ⎥
⎢  ⎛1       b⋅c     ⎞       c⋅d         d          b⋅c    ⎥
⎢c⋅⎜─ + ────────────⎟ - ───────────  ─────── - ───────────⎥
⎢  ⎜a    2 ⎛    b⋅c⎞⎟     ⎛    b⋅c⎞      b⋅c     ⎛    b⋅c⎞⎥
⎢  ⎜    a ⋅⎜d - ───⎟⎟   a⋅⎜d - ───⎟  d - ───   a⋅⎜d - ───⎟⎥
⎣  ⎝       ⎝     a ⎠⎠     ⎝     a ⎠       a      ⎝     a ⎠⎦

It doesn't evaluate it to direct eye(2) but if you take elements, and simplify them, you'll see that they this messy matrix is actually and 2x2 identity matrix.

A pythonic way to check that (knowing given):

>>> N[0]
a*(1/a + b*c/(a**2*(d - b*c/a))) - b*c/(a*(d - b*c/a))
>>> N[1]
0
>>> N[3]
d/(d - b*c/a) - b*c/(a*(d - b*c/a))
>>> N[2]
c*(1/a + b*c/(a**2*(d - b*c/a))) - c*d/(a*(d - b*c/a))

>>> def will_evaluate_one(a,b,c,d):
...    return a*(1/a + b*c/(a**2*(d - b*c/a))) - b*c/(a*(d - b*c/a)) #N[0]     
...
>>> will_evaluate_one(1,2,3,9)
1
>>> will_evaluate_one(1,2,3,19)
1
>>> will_evaluate_one(1,2,23,19)
1
>>> will_evaluate_one(1,12,23,19)
1

>>> def will_also_evaluate_one(a,b,c,d):
...     return d/(d - b*c/a) - b*c/(a*(d - b*c/a)) #N[1] 
... 
>>> will_also_evaluate_one(2,4,5,6)
1
>>> will_also_evaluate_one(2,4,15,6)
1
>>> will_also_evaluate_one(2,14,15,6)
1
>>> will_also_evaluate_one(12,14,15,6)
1

Note: I've just realised that, sympy uses anlaytic inversion formula. See here: https://en.wikipedia.org/wiki/Helmert%E2%80%93Wolf_blocking

marmeladze
  • 6,468
  • 3
  • 24
  • 45