magnus.expmkernels

expmkernels.py

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

This module is the 'numba' backend of magnus.magnus._expm_stack. For 2x2 and 3x3 matrices 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 vectorizing 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 4 and 5: Jacobi, not a closed form

There is no practical closed form for the eigenvalues of a 4x4 or 5x5 Hermitian matrix, and for a long time that sentence ended “so 4nu and 5nu keep the eigh path”. The conclusion did not follow: what made those dimensions slow was never the missing closed form but eigh’s fixed per-matrix LAPACK overhead, about 2.3 us on a 4x4 – two thirds of a whole d=4 Magnus pass. So 4x4 and 5x5 stacks go to _jacobi_expm_core instead, a batched cyclic Jacobi eigensolver that warm-starts each matrix from its predecessor’s eigenvectors and re-orthonormalizes that basis at every step; see its docstring for the scheme and for which of its details are load-bearing. Unlike the closed forms it is iterative, so this backend replacement is not bit-identical to what it replaces – it is held to the same accuracy class as eigh instead (within 5.5x at every norm, clustering and degeneracy measured, against the same yardstick that admits the 3x3 closed form at up to 10x). supports_dim() is the single place that decides the routing.

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 characteriztic 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 compiled 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 behave, and eigh answers beyond it. What the gate admits, worst over two spectrum families x 6 separations x 40 random bases per rung (docs/dev/calibrate_sev_tol.py):

m         worst |closed - expm|    eigh, same cells
4e1     1.0e-14                  1.4e-15
0e2     8.1e-14                  2.0e-14
1e3     2.1e-13                  3.9e-14
4e3     8.5e-13                  7.8e-14
0e4     2.0e-12                  6.1e-14

The worst measured is 2.0e-12 just under the gate and 2.3e-13 in the \(m \le 1.1\times10^3\) corner, so the guarantee is stated with headroom at 5e-12 absolute across everything the gate admits, and 5e-13 in that corner. The headroom is deliberate: these are worst-over-random-bases quantities and more sampling keeps finding slightly worse ones, which is exactly how the previous claim came to be false. Both bounds are far below any tolerance this package is asked for, and far below where the closed form runs away past the gate: cells at \(m \ge 1.1\times10^5\) reach 131x eigh, and at \(m \sim 4\times10^9\), 7440x.

An earlier version of this docstring claimed 2e-13 across the admitted range. That was measured on one corner of it and is not true of the rest; at \(m = 1.1\times10^3\) itself about 1% of random bases exceed it (11 of 1200, worst 2.3e-13). No value of this constant could have rescued that claim, because test_sev_tol_sits_inside_its_calibrated_window pins \(m = 1.1\times10^3\) as a cell that must stay on the kernel, so the gate cannot be lowered past the point where the claim already fails. The number was corrected instead.

Read m, never “spectral scale”. Two calibrations of this constant appeared to contradict each other – one finding the first unsafe cell at \(1.1\times10^5\), the other at \(4.4\times10^3\) – purely because they used different spectrum families and both called the result “scale \(10^2\)”. \([-s, -s(1-d), s]\) spans \(2s\) and gives \(m \simeq 0.44\,s^2\); \([0, d, S]\) spans \(S\) and gives \(m \simeq 0.11\,S^2\) – a factor of four in \(m\) at the same nominal scale. Compared at equal \(m\) the two families agree to within their sampling scatter, and the disagreement dissolves. \(m = \mathrm{tr}(X^2)/6\) is a spectral invariant; the word “scale” is not, and is what made this look like a contradiction.

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 compiled kernel.

True for 2 through 5. Dimensions 2 and 3 have Cayley-Hamilton closed forms; 4 and 5 go to the batched Jacobi eigensolver _jacobi_expm_core instead. An earlier version of this docstring reasoned that a 4x4 or 5x5 Hermitian eigenproblem has no practical closed form and concluded that 4nu and 5nu stay on eigh. The premise stands; the conclusion did not follow from it, because the missing closed form was never what made those dimensions slow – eigh’s fixed per-matrix LAPACK overhead (~2.3 us on a 4x4, two thirds of a d=4 Magnus pass) was, and an iterative solver with no such overhead removes it without any closed form. 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 through 5 (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 could not be answered at full accuracy – too ill-conditioned for the 3x3 closed form (see _ch3_core), or at the Jacobi sweep cap for a 4x4/5x5 (see _jacobi_expm_core) – and the caller should recompute with eigh.
Return type:: tuple
Raises:: ValueError – If d is not 2 through 5. Without this the else below handed unsupported 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 uninitialized 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.