Matrix multiplication gives unsual result in Python (SciPy/PyLab)

Question

I'm new to Python, and a bit rusty with my linear algebra, so perhaps this is a simple question. I'm trying to implement a Taylor Series expansion on a Matrix to compute exp(A), where A is just a simple 3x3 matrix. The formula, BTW for this expansion is sum( A^n / n! ).

My routine works alright up to n=9, but at n=10, the numbers in the Matrix suddenly become negative. This is the problem.

A**9 matrix([[ 250130371, 506767656, 688136342], [ 159014912, 322268681, 437167840], [ 382552652, 775012944, 1052574077]])

A**10 matrix([[-1655028929, 1053671123, -1327424345], [ 1677887954, -895075635, 319718665], [ -257240602, -409489685, -1776533068]])

Intuitively A^9 * A should produce larger numbers for each member of the matrix, but as you can see, A^10 isn't giving that result.

Any ideas?

from scipy import *
from numpy import *
from numpy.linalg import *
#the matrix I will use to implement exp(A)
A = mat('[1 3 5; 2 5 1; 2 3 8]')
#identity matrix
I = mat('[1 0 0; 0 1 0; 0 0 1]')
#first step in Taylor Expansion (n=0)
B = I
#second step in Taylor Expansion (n=1)
B += A
#start the while loop in the 2nd step
n = 2
x=0
while x<10:
    C = (A**n)/factorial(n)
    print C
    print " "
    n+=1
    B+= C
    print B
    x+=1

print B

Thanks for any help you can give!

Answer 1

Your matrix is created with elements of type int32 (32-bit integer). You can see this by printing the value of A.dtype. 32-bit integers can only hold values up to about 2 billion, so after that they will wrap around to negative values.

If 64-bit integers are large enough, you can use them instead:

A = mat('[1 3 5; 2 5 1; 2 3 8]', dtype=numpy.int64)

Otherwise, you can use floating point numbers. They have a much larger maximum value, but limited precision, so there may be some inaccuracies.

A = mat('[1 3 5; 2 5 1; 2 3 8]', dtype=float)

In this case floating-point is probably the best choice, since you don't want your results to be integers after dividing by n!.

Answer 2

I don't know much about scientific python but I do know what is going wrong. It seems that the matrix elements are represented as 32 bit signed integers. That means that they are restricted to the range -2^31 <= x < 2^31. At A^10, the numbers become too big and they "wrap around". I checked and the top-left coefficient actually is 2639938267, which, when wrapped around, gives 2639938267 - 2^32 = -1655028929.

I don't know how to set the data type in python, so I don't know how you can solve this. I'm sure it's possible, though.

(I could also suggest you try sage: www.sagemath.org, which is mathematics software based on python. It automatically uses infinite precision. It's how I checked these numbers just now).

Good luck! Timo