magnus.expmkernels

expmkernels.py

Compiled Cayley-Hamilton kernels for \(\exp(-iK)\), K Hermitian.

This module is the 'numba' backend of magnus.magnus._expm_stack. It computes the matrix exponential of a stack of small Hermitian matrices without an eigenvector solver, by applying to \(K\) the polynomial that interpolates \(\exp(-i\lambda)\) on the spectrum of \(K\):

\[\exp(-iK) = a_0 I + a_1 K + a_2 K^2 .\]

Cayley-Hamilton guarantees such a polynomial exists (degree \(d-1\) for a \(d \times d\) matrix); the eigenvalues are obtained in closed form.

Why this is worth a compiled kernel

np.linalg.eigh costs about 1.25 us per 3x3 regardless of stack size – 1, 108 or 4096 matrices, the per-matrix cost is flat – because it loops over LAPACK internally rather than vectorising over the stack. On a 108-slab Magnus pass that single call is roughly a quarter of the total.

The same algebra written in pure numpy does not help below stacks of about a hundred: it is some twenty numpy calls, each paying dispatch overhead on arithmetic that is otherwise trivial, so it loses at small stacks and wins only mildly at large ones. Only a compiled kernel removes the dispatch, which is why numba is used here and why a numpy version of these formulas is not offered as a third backend.

Dimensions 2 and 3 only

There is no practical closed form for the eigenvalues of a 4x4 or 5x5 Hermitian matrix, so 4nu and 5nu keep the eigh path. They stay correct; they do not get faster. supports_dim() is the single place that decides this.

Why the interpolation form is safe at a degeneracy

Coincident eigenvalues are the whole numerical risk in a Cayley-Hamilton scheme, because the interpolation coefficients divide by eigenvalue differences. Two facts remove it here, and both are load-bearing enough to state:

A Hermitian matrix is never defective. Its minimal polynomial has simple roots even when its characteristic polynomial does not, so a polynomial matching \(\exp(-i\lambda)\) on the distinct eigenvalues already reproduces the function exactly. The confluent (Hermite) form, which matches derivatives as well and is unavoidable for a general matrix, is not needed for this one. Nothing in this module differentiates anything.

The ill-conditioned coefficient multiplies a correspondingly small matrix. Write the interpolant in Newton form on eigenvalues sorted ascending, with the spectrum shifted so the median eigenvalue sits at zero:

\[\exp(-iZ) = f[z_0] I + f[z_0, z_1](Z - z_0 I) + f[z_0, z_1, z_2](Z - z_0 I)(Z - z_1 I) ,\]

with \(Z = K - \lambda_1 I\) and \(z_0 \le z_1 = 0 \le z_2\). The first divided difference is evaluated as \(f[a, b] = -i e^{-i(a+b)/2}\, \mathrm{sinc}((a-b)/2)\), which is cancellation-free for every pair including \(a = b\) (see _sinc). The second, \((f[z_0,z_1] - f[z_1,z_2])/(z_0 - z_2)\), does lose digits as the nodes coalesce – its absolute error grows like \(\epsilon/(z_2 - z_0)\) – but the matrix it multiplies has norm at most \((z_2 - z_0)^2\), so the product’s error is bounded by \(\epsilon\,(z_2 - z_0)\) and vanishes with the gap. Sorting is what makes this true: it is what guarantees that a small \(z_2 - z_0\) means all three eigenvalues are close, rather than one unlucky pair out of three.

So there is no tolerance, no crossover, and no near-degenerate branch to place correctly. The only guard is for \(z_0 = z_2\) exactly, where the term is multiplied by the zero matrix and is simply dropped.

Measured against scipy.linalg.expm, the error is 1e-16 at splittings of 1e-2, 1e-6, 1e-10, 1e-14 and exactly zero alike. The closed-form eigenvalues are much worse than that near a degeneracy – they degrade to ~1e-9, because \(\arccos\) has infinite derivative at the ends of its range, which is where a repeated root sits – and at \(\lVert K \rVert \sim 1\) it does not matter: the interpolation error is second order in the displacement of a coalescing node, so a node that is 1e-9 off contributes 1e-18. A test asserts both halves of that.

That argument has a range of validity, and an earlier version of this paragraph did not say so. The eigenvalue error scales with the norm, so the second-order suppression is fighting a term that grows: where a clustered spectrum meets a large norm the closed form reaches 2.7e-07 against eigh’s 3.0e-11, a factor of 7440. Neither a sweep over norms at generic separation (ratios 0.4-2.5) nor a sweep over separations at norm 1 (0.4-1.0) visits that corner, which is how the unqualified claim came to be written and believed. SEV_TOL is the gate that keeps it out of reach, and tests/test_expm_backend.py now crosses the two axes so the corner cannot go unmeasured again.

A note on the determinant

The cross term of \(\det X\) for Hermitian X is \(2\,\mathrm{Re}(X_{01} X_{12} \overline{X_{02}})\). Moving that conjugate to either of the other two factors produces matrices that are still nearly unitary and still plausible, and wrong by O(1). An earlier prototype of this kernel shipped that exact transposition and reported a 6x speed-up while returning garbage; _ch3_core carries the term explicitly and tests/test_expm_backend.py pins it against np.linalg.det.

Added in version 1.0.0.

Routine listings

HAVE_NUMBA - Whether the compiled kernels are available

SEV_TOL - Conditioning above which eigh answers instead

supports_dim - Whether a given matrix dimension has a kernel

expm_herm_stack - exp(-iK) and the eigenvalues of K, for a stack

Attributes

`HAVE_NUMBA`	Whether numba imported, and so whether the compiled kernels exist.
`SEV_TOL`	Above this \(m = \mathrm{tr}(X^2)/6\), a 3x3 is handed back to `eigh`.

Functions

`supports_dim`(→ bool)	Returns whether dimension `d` has a Cayley-Hamilton kernel.
`expm_herm_stack`(→ tuple)	Returns \(\exp(-iK)\) and the eigenvalues of K, for Hermitian K.

Module Contents

magnus.expmkernels.HAVE_NUMBA[source]

Whether numba imported, and so whether the compiled kernels exist.

False leaves every backend decision to eigh; nothing else in the package changes. numba is an optional dependency (pip install magnuspy[fast]).

Type:: bool

magnus.expmkernels.SEV_TOL = 10000.0[source]

Above this \(m = \mathrm{tr}(X^2)/6\), a 3x3 is handed back to eigh.

The closed-form solve is used only in the range of \(\lVert K \rVert\) where it is measured to match eigh, and eigh answers beyond it. Calibrated against the separation-by-scale grid in tests/test_expm_backend.py: no cell at spectral scale \(\le 10^2\) (\(m \le 1.1\times10^3\)) is worse than eigh by more than 2e-13 absolute, while cells at scale \(10^3\) (\(m \ge 1.1\times10^5\)) reach 131x and at \(10^5\) reach 7440x. This sits between them with a factor-100 margin either way.

Why the scale and not the clustering, which is the actual mechanism. The damage needs a clustered spectrum and a large norm together: \(\arccos\) has infinite derivative at \(u = \pm 1\), so clustering turns rounding in \(u\) into an eigenvalue error \(\sim\sqrt{\epsilon}\,\lVert K \rVert\), which only matters once the norm is large. But the clustering half cannot be gated on, because the danger is a band rather than a tail: at exact degeneracy the pair comes out bit-identical and the answer is fine (measured 0.3-1.8x of eigh at \(u = \pm 1\) exactly), the damage sits at intermediate separations, and \(1/(1-u^2)\) is largest exactly where there is no problem. A one-sided threshold on it therefore cannot work – verified by calibration, which found no separating value.

So this gate is deliberately conservative rather than tight: it also declines large-norm spectra that are not clustered and would have been fine (measured 0.6x of eigh). That costs speed on those, never accuracy, and it costs nothing where the speed comes from – a Magnus slab has \(\lVert\Omega\rVert \lesssim \pi\) by construction, so slab chains are never declined, and an ordinary 3nu constant-density or vacuum call measures \(\lVert K \rVert \approx 4\). What it does decline is the large accumulated phase of an eV-scale sterile splitting, where accuracy is worth more than the microsecond.

Added in version 1.0.0.

Type:: float

magnus.expmkernels.supports_dim(d: int) → bool[source]

Returns whether dimension d has a Cayley-Hamilton kernel.

True for 2 and 3. A 4x4 or 5x5 Hermitian eigenproblem has no practical closed form, so 4nu and 5nu stay on eigh: correct, and no faster than before. This is the one place that decision is made.

Parameters:: d (int) – Matrix dimension.
Returns:: Whether a kernel exists for that dimension.
Return type:: bool

magnus.expmkernels.expm_herm_stack(K: numpy.ndarray) → tuple[source]

Returns \(\exp(-iK)\) and the eigenvalues of K, for Hermitian K.

The drop-in replacement for the eigh half of magnus.magnus._expm_stack: it returns the eigenvalues alongside the exponential because the caller needs them anyway, for the slab-width convergence warning, and this way there is no second spectral computation.

Parameters:: K (np.ndarray) – Hermitian matrix or stack of them, shape (…, d, d), with d 2 or 3 (see supports_dim()).
Returns:: \(\exp(-iK)\), shape (…, d, d); the eigenvalues in ascending order, shape (…, d); and a float conditioning severity for the whole stack. A severity above SEV_TOL means at least one matrix was too ill-conditioned for the closed-form solve and the caller should recompute with eigh; see _ch3_core.
Return type:: tuple
Raises:: ValueError – If d is not 2 or 3. Without this the else below handed 4x4 and 5x5 input to the 3x3 kernel, which returned no exception, an error of 2.4 against scipy.linalg.expm, a unitarity violation of 11.3, and uninitialised memory in the fourth eigenvalue – and segfaulted at d=1 by indexing K[i,2,1] with numba’s bounds checking off. supports_dim() is documented as the single place that decides which dimensions are handled, and now actually is.