Matrix exponential in Python

Question

Given the equation,

$$\begin{equation*}\exp(Ax)=\sum_{n=0}^{\infty}\frac{x^n}{n!}A^n\end{equation*}$$

implement in Python a function which calculates $\exp(Ax)$ for a given matrix $A$ and scalar $x$ using the equation directly. The only part that is troubleieng me is the equation. How can I implement an infinite sum in Python?

I'm allowed to test my result with the expm from scipy.linalg but I have to directly use the equation.

I have written the following code with help from previous answers.

import math
import numpy as np
from scipy.linalg import expm

# Scalar x (will later on be for user input)
x = 1

matrix = np.array([[-5, 2, 3], [2, -6, 4], [4, 5, -9]])

# Using scipy to compute the matrix exponential (for comparison)
    
B = expm(matrix)
print(B)    

# Defining the equation
def equation(n):
    y = ((pow(x, n) * np.linalg.matrix_power(matrix, n)) / int(math.factorial(n)))
    return y


# Summing the equation with finite iterations
result = sum([equation(n) for n in range(0, 1000)])
print(result)

But it doesn't seem to give me the right answer. After hours of staring I can't seem to figure out the problem.

Answer 1

Stepping through some calls to other functions, the crucial part of the source code is here. @zwim provided a hint of how matrix exponentiation can be reduced to exponentiating scalars, but either way the basic answer, as @saulspatz noted, is that you can just add terms until new ones are so small they can be neglected.

For what it's worth, the not-so-basic answer (probably beyond the scope of what you need to implement in this exercise), which you'll see hinted at in the code I linked to, is to first check whether Padé approximants are an adequate, more efficient alternative than partial sums of the exponential map. The algorithmic details get quite thorny beyond that point - in fact, I suspect some research into whether we can do something even better is still ongoing - but their motive is quite simple: the obvious choice of a computable iteration converging to a given limit may not be the one that converges fastest.

Right, so how do we calculate $\exp(xA)$? One might exploit repeated squaring to get away with fewer terms, e.g.

def exp(A, x):
    Ax = A*x
    if max(np.absolute(Ax))>1e-4:
        sqrt = exp(A, x/2)
        return sqrt @ sqrt
    sq = Ax @ Ax
    return 1 + Ax + sq / 2 + sq @ Ax / 6

Note the above code doesn't explicitly check new terms are very small; it checks the matrix being exponentiated is very small. There's much about this code that has room for optimizations, but I leave that to you.

Answer 2

Here's a naive implementation:

def my_expm(A):
    m,n = np.shape(A)
    norm = np.sqrt(np.trace([email protected]))
    k = 0
    fact = 1
    matpow = np.eye(n)
    expmat = np.matrix(np.zeros([n,n]))
    lastnorm = 1
    while lastnorm > 1e-8:
        expmat += matpow/fact
        lastnorm = np.max(abs(matpow/fact))
        k += 1
        fact *= k
        matpow = A@matpow
    return expmat

A = np.matrix('0, -1; 1 0')*np.pi
print(my_expm(A))