[SciPy-User] Faster approach to calculate the matrix exponential

[Pauli Virtanen]

26.06.2014 14:07, Troels Emtekær Linnet kirjoitti:
> Using the scripy method does not sound viable.
> 
> It was discussed here,  post 2
> 
> http://thread.gmane.org/gmane.science.nmr.relax.devel/6446/match=eliminating+86+bottleneck+numeric+r1rho

That link doesn't load currently.

Nevertheless, using the Scipy implementation inside for loops is faster
than using the eig() based approach, with a crossover (on my machine)
when the inner matrix size is bigger than 16x16. This comparison applies
to Scipy 0.14.0 --- the implementation in earlier Scipy versions is less
efficient.

Benchmark:
"""
from scipy.linalg import expm

def matrix_exponential_scipy(A, dtype=None):
    B = np.empty(A.shape, dtype=A.dtype)
    for i in range(A.shape[0]):
        for j in range(A.shape[1]):
            for k in range(A.shape[2]):
                for l in range(A.shape[3]):
                    for m in range(A.shape[4]):
                        B[i,j,k,l,m] = expm(A[i,j,k,l,m])
    return B

import numpy as np

N = 18
A = np.random.rand(3, 4, 5, 3, 2, N, N)

start_1 = time.time()
B_1 = matrix_exponential_scipy(A)
end_1 = time.time()

start_2 = time.time()
B_2 = matrix_exponential_rank_NE_NS_NM_NO_ND_x_x(A)
end_2 = time.time()

print("Scipy:", end_1 - start_1)
print("Numpy:", end_2 - start_2)
"""
#-> ('Scipy:', 0.3015120029449463)
#-> ('Numpy:', 0.31020593643188477)


-- 
Pauli Virtanen