Implementation of Highams Algorithm 10.20 (scaling and squaring algorithm). This algorithm evaluates the matrix exponential X = e^A of A ∈ C^{n×n} using the scaling and squaring method.

scale_and_square.py

import numpy as np
import scipy.linalg as sla

theta_m = {3: 1.5e-2, 5: 5.4e-1, 7: 9.5e-1, 9: 2.1e0, 13: 5.4e0}

pade_coefficients = {
    3: [120, 60, 12, 1],
    5: [30240, 15120, 3360, 420, 30, 1],
    7: [17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1],
    9: [17643225600, 8821612800, 2075673600, 302702400, 30270240, 2162160, 110880, 3960, 90, 1],
    13: [64764752532480000, 32382376266240000, 7771770303897600, 1187353796428800, 129060195264000, 10559470521600, 670442572800, 33522128640, 1323241920, 40840800, 960960, 16380, 182, 1],
}

def pade_approximant(A, m):
    n, _ = A.shape
    I = np.eye(n)
    coeffs = pade_coefficients[m][::-1]
    # Precompute A^2 to speed up Horners rule
    A2 = A @ A
    U = I * coeffs[1]
    V = I * coeffs[0]
    for i in range(1, m//2+1):
        U = U @ A2 + I * coeffs[2*i+1]
        V = V @ A2 + I * coeffs[2*i]
    V = V @ A
    return U, V

def scale_and_square(A):
    # Determine the order m of the Pade approximant to use
    norm_A = np.linalg.norm(A, ord=1)
    m = max(m_val for m_val, theta in theta_m.items() if norm_A < theta)
    # Calculate the minimum integer s such that ||A/2^s|| < theta_m
    s = max(0, int(np.ceil(np.log2(norm_A / theta_m[m]))))
    # Scale down A by 2^s
    A_scaled = A / (2**s)
    # Calculate the Pade approximant of the scaled matrix
    U, V = pade_approximant(A_scaled, m)
    # X = (U + V)(U - V)^(-1)
    X = np.linalg.solve(U - V, U + V)
    # Square X repeatedly to undo the scaling
    for _ in range(s):
        X = X @ X
    return X

test_matrix = np.random.rand(3, 3)
matrix_exponential = scale_and_square(test_matrix)
print(matrix_exponential)
print(sla.expm(test_matrix))