matrix_exp and matrix_exp_frechet#

The matrix exponential extends e**x to square matrices: exp(A) = I + A + A**2/2! + A**3/3! + .... It solves systems of linear differential equations (x' = A x has the solution x(t) = exp(t A) x(0)), turns rotation generators into rotations, and gives the transition probabilities of continuous-time Markov chains. Its Frechet derivative says how exp(A) changes when A moves a little in a direction E, which is what a gradient through the matrix exponential needs. Implementations differ widely: truncated Taylor or Pade series of any degree, with any number of scaling and squaring steps.

In ExactRef, er.matrix_exp(A), also exactref.linalg.matrix_function.matrix_exp, is the exponential of each square matrix in A, as in torch.linalg.matrix_exp. er.matrix_exp_frechet(A, E), also exactref.linalg.matrix_function.matrix_exp_frechet, is its Frechet derivative at A in the direction E. Both take real or complex inputs.

Their tolerance is a backward error, because no rounding step is shared by all algorithms: a correct result is exp(A + dA) for a small dA, rounded to the compute format and then to the output format. They have no normal answer. Guarantees and limits lists implementations known to exceed the tolerance.

Arguments#

Argument

Meaning

Default

A

square matrices, shape (..., n, n), real or complex

required

E (matrix_exp_frechet only)

the direction, with A's shape

required

dtype

the output format, of the same kind as the inputs

the input format

rounding

the rounding of the reference and of the final step into the output format

er.Rounding.NEAREST_EVEN

overflow

what a result past the output format's range becomes: "ieee", "saturate" or "nan" for floats; "wrap", "error" or "saturate" for integers

each format's own rule; integers wrap

Examples#

Code
import math
import exactref as er

# A 2x2 float64 matrix
A = er.asarray([[0.0, 1.0], [-1.0, 0.0]], "float64")
# Its exponential, each entry correctly rounded to float64
m = er.matrix_exp(A)

# exp(A) is a rotation by 1 radian: [[cos 1, sin 1], [-sin 1, cos 1]]
print("rounded:")
for row in m.rounded.tolist():
    print(row)

# How far a correct float64 output can be from rounded, first row
print("tolerance, row 0:", m.tolerance().tolist()[0])
# ulps=4 allows an input error dA 4 times as large
print("tolerance(ulps=4), row 0:", m.tolerance(ulps=4).tolist()[0])
Output
rounded:
[0.5403023058681398, 0.8414709848078965]
[-0.8414709848078965, 0.5403023058681398]
tolerance, row 0: [4.218847493575595e-15, 4.107825191113079e-15]
tolerance(ulps=4), row 0: [1.6653345369377348e-14, 1.6653345369377348e-14]

ulps (default 1) scales the backward error. Both standards give the same tolerance.

Batches and complex matrices#

Code
# A batch of two 2x2 float64 matrices: zero and a diagonal matrix
zero, diagonal = [[0.0, 0.0], [0.0, 0.0]], [[1.0, 0.0], [0.0, 2.0]]
batch = er.asarray([zero, diagonal], "float64")

# The exponential of each matrix in the batch
exp_zero, exp_diagonal = er.matrix_exp(batch).rounded.tolist()
# exp of the zero matrix is the identity
print("exp(zero):", exp_zero)
# exp of a diagonal matrix: e**d for each diagonal entry d
print("exp(diagonal):", exp_diagonal)

# The 1x1 complex128 matrix [[i * pi]]
z = er.asarray([[1j * math.pi]], "complex128")
# Its exponential, close to exp(i * pi) = -1
print("exp(i * pi):", er.matrix_exp(z).rounded.tolist())
Output
exp(zero): [[1.0, 0.0], [0.0, 1.0]]
exp(diagonal): [[2.718281828459045, 0.0], [0.0, 7.38905609893065]]
exp(i * pi): [[(-1+1.2246467991473532e-16j)]]

The imaginary part is not 0 because the input is the float64 value nearest pi, not pi itself.

Special values#

Any non-finite entry makes the whole matrix NaN:

Code
# A 2x2 float64 matrix with one infinite entry
with_inf = er.asarray([[1.0, float("inf")], [0.0, 1.0]], "float64")
# Its exponential, rounded to float64
print("rounded:", er.matrix_exp(with_inf).rounded.tolist())
Output
rounded: [[nan, nan], [nan, nan]]

Checking PyTorch#

Code
import random
import torch

# Seeded, so every run uses the same values
rng = random.Random(0)
# A 4x4 float64 tensor of normal values, mean 0 and standard deviation 0.5
values = [[rng.gauss(0, 0.5) for _ in range(4)] for _ in range(4)]
T = torch.tensor(values, dtype=torch.float64)

# Could a correct implementation produce every element of PyTorch's exp(T)?
print("covered:", er.matrix_exp(T).covers(torch.linalg.matrix_exp(T)).all())
Output
covered: True

Backward#

er.matrix_exp.backward(A, grad) returns grad_A = matrix_exp_frechet(A^H, grad), which is what PyTorch computes.

Code
# An output gradient: 1 in the top-left entry, 0 elsewhere
E = er.asarray([[1.0, 0.0], [0.0, 0.0]], "float64")

# The gradient at A for the output gradient E
print("grad_A:")
for row in er.matrix_exp.backward(A, E).rounded.tolist():
    print(row)
Output
grad_A:
[0.6908866453380181, -0.42073549240394825]
[0.42073549240394825, -0.1505843394698784]

matrix_exp_frechet#

matrix_exp_frechet(A, E) is the top-right block of exp([[A, E], [0, A]]):

Code
# The Frechet derivative at A in the direction E
print("matrix_exp_frechet(A, E):")
for row in er.matrix_exp_frechet(A, E).rounded.tolist():
    print(row)

# The 4x4 block matrix [[A, E], [0, A]]
block = er.asarray([[0.0, 1.0, 1.0, 0.0], [-1.0, 0.0, 0.0, 0.0],
                    [0.0, 0.0, 0.0, 1.0], [0.0, 0.0, -1.0, 0.0]], "float64")

# Its exponential's top-right 2x2 block: rows 0 and 1, columns 2 and 3
print("top-right block of matrix_exp(block):")
for row in er.matrix_exp(block).rounded.tolist()[:2]:
    print(row[2:])
Output
matrix_exp_frechet(A, E):
[0.6908866453380181, 0.42073549240394825]
[-0.42073549240394825, -0.1505843394698784]
top-right block of matrix_exp(block):
[0.6908866453380181, 0.42073549240394825]
[-0.42073549240394825, -0.1505843394698784]

It has no backward of its own.

Full specification#

er.matrix_exp#

matrix_exp(A, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, overflow: str | None = None) -> Result#

Exponential of each square matrix in A (shape (..., n, n)), real or complex.

Every entry is the exact value correctly rounded to dtype. Entries are rounded from certified enclosures computed with python-flint, refined until both ends round alike. Each part of an entry is rational or transcendental (Lindemann-Weierstrass). A transcendental part is never on a rounding boundary, so refining it ends; a rational part can be, so parts still unrounded after the first refinements are tested for being rational with exact arithmetic, and rational parts are rounded from their exact values. An exact zero is +0, or -0 when rounding toward -inf. Any non-finite entry makes the whole matrix NaN.

tolerance() uses a backward-error model, because no internal boundary is shared by all algorithms (Taylor or Pade, any degree, any number of squarings): the computed result is exp(A + dA) with ||dA||_inf <= ulps * (n + 16) * u * max(||A||_inf, 1) (ulps defaults to 1, u is the compute format's unit roundoff; n and 16 count the roundings of an algorithm's inner products and its others, and the floor of 1 covers the rounding of exp(A)'s identity part, which every algorithm makes), rounded to the compute format and then to dtype. This is an assumption, measured on random matrices (n = 2 to 64, ||A||_inf 0.03 to 32): Higham's 2005 scaling and squaring with its degree-13 Pade approximant, and PyTorch's matrix_exp, stay within half of it. PyTorch's own results at small norms (float64 near 0.03, float32 up to 1) and scipy.linalg.expm 1.18.1 for some float64 matrices from ||A||_inf 4 lie outside it; there the reference matches mpmath at high precision. A may also be an earlier result, as Result describes: the model then holds for every A within its tolerance.

er.matrix_exp.backward#

backward(A, grad, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, overflow: str | None = None) -> Result#

Gradient of matrix_exp(A) for the output gradient grad.

grad_A = matrix_exp_frechet(A^H, grad) (^H: conjugate transpose), which is what PyTorch computes. grad must have A's shape and dtype. A and grad may also be earlier results.

er.matrix_exp_frechet#

matrix_exp_frechet(A, E, *, dtype=None, rounding: Rounding = Rounding.NEAREST_EVEN, overflow: str | None = None) -> Result#

Frechet derivative of the matrix exponential at A in the direction E.

L = d/dt exp(A + t E) at t = 0, for each square matrix in A and E (shape (..., n, n), the same for both), real or complex. L is the top-right block of exp([[A, E], [0, A]]), so this is matrix_exp of that 2n x 2n block matrix, restricted to the block (a proxy): every entry is correctly rounded once, and exact and tolerance() are matrix_exp's for those entries. The tolerance is matrix_exp's backward-error model on the block matrix, which is how PyTorch computes matrix_exp's backward. A and E may also be earlier results.

See also#