expm/benchmarks/pure_python_vanilla.py at main · john-science/expm · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
import numpy as np
from numpy.linalg import norm, solve
from numpy import eye, dot, ceil, log2


def expm(A):
    """ Compute the matrix exponential using Pade approximation
    https://github.com/rngantner/Pade_PyCpp/blob/master/src/expm.py

    Args:
        A: Matrix (shape(M,M)) to be exponentiated
    Returns:
        np.array: Matrix (shape(M,M)) exponential of A
    """
    A_L1 = norm(A, 1)
    n_squarings = 0

    if A_L1 < 1.495585217958292e-2:
        U, V = _pade3(A)
    elif A_L1 < 2.539398330063230e-1:
        U, V = _pade5(A)
    elif A_L1 < 9.504178996162932e-1:
        U, V = _pade7(A)
    elif A_L1 < 2.097847961257068:
        U, V = _pade9(A)
    else:
        n_squarings = max(0, int(ceil(log2(A_L1 / 5.371920351148152))))
        A = A / 2 ** n_squarings
        U, V = _pade13(A)

    P = U + V  # p_m(A) : numerator
    Q = -U + V  # q_m(A) : denominator
    R = solve(Q, P)

    # squaring step to undo scaling
    for _ in range(n_squarings):
        R = dot(R, R)

    return R


def _pade3(A):
    """ Helper method for expm

    Args:
        A (np.array): input Matrix
    Returns:
        tuple: Mystery Components
    """
    b = tuple([120., 60., 12., 1.])
    shape = (len(A[:, 0]), len(A[0, :]))
    ident = eye(*shape, dtype='float64')
    A2 = dot(A, A)
    U = dot(A, (b[3] * A2 + b[1] * ident))
    V = b[2] * A2 + b[0] * ident
    return U, V


def _pade5(A):
    """ Helper method for expm

    Args:
        A (np.array): input Matrix
    Returns:
        tuple: Mystery Components
    """
    b = tuple([30240., 15120., 3360., 420., 30., 1.])
    shape = (len(A[:, 0]), len(A[0, :]))
    ident = eye(*shape, dtype='float64')
    A2 = dot(A, A)
    A4 = dot(A2, A2)
    U = dot(A, b[5] * A4 + b[3] * A2 + b[1] * ident)
    V = b[4] * A4 + b[2] * A2 + b[0] * ident
    return U, V


def _pade7(A):
    """ Helper method for expm

    Args:
        A (np.array): input Matrix
    Returns:
        tuple: Mystery Components
    """
    b = tuple([17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.])
    shape = (len(A[:, 0]), len(A[0, :]))
    ident = eye(*shape, dtype='float64')
    A2 = dot(A, A)
    A4 = dot(A2, A2)
    A6 = dot(A4, A2)
    U = dot(A, b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * ident)
    V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * ident
    return U, V


def _pade9(A):
    """ Helper method for expm

    Args:
        A (np.array): input Matrix
    Returns:
        tuple: Mystery Components
    """
    b = tuple([17643225600., 8821612800., 2075673600., 302702400., 30270240., 2162160., 110880., 3960., 90., 1.])
    shape = (len(A[:, 0]), len(A[0, :]))
    ident = eye(*shape, dtype='float64')
    A2 = dot(A, A)
    A4 = dot(A2, A2)
    A6 = dot(A4, A2)
    A8 = dot(A6, A2)
    U = dot(A, b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * ident)
    V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * ident
    return U, V


def _pade13(A):
    """ Helper method for expm

    Args:
        A (np.array): input Matrix
    Returns:
        tuple: Mystery Components
    """
    b = tuple([64764752532480000., 32382376266240000., 7771770303897600., 1187353796428800., 129060195264000.,
               10559470521600., 670442572800., 33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.])
    shape = (len(A[:, 0]), len(A[0, :]))
    ident = eye(*shape, dtype='float64')
    A2 = dot(A, A)
    A4 = dot(A2, A2)
    A6 = dot(A4, A2)
    U = dot(A, dot(A6, b[13] * A6 + b[11] * A4 + b[9] * A2) + b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * ident)
    V = dot(A6, b[12] * A6 + b[10] * A4 + b[8] * A2) + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * ident
    return U, V