Magnus Expansion Terms to Any Order

The numerical core evaluates the Magnus expansion with the coefficient of every commutator group written out explicitly in Python. That is fast, and it keeps the low-order expressions checkable by eye, but it also means the highest available order is whatever was typed in. magnus.expansionterms derives the same terms from the recursion itself, in exact rational arithmetic, at any order.

That serves two purposes. It lets the hard-coded coefficients be checked rather than trusted – the test suite regenerates them and compares – and it lets you inspect an order beyond the implemented ceiling without running anything.

The recursion

Writing \(A(t) = -i H(t)\), the Magnus expansion sums terms \(\Omega = \sum_n \Omega_n\) generated by the Bernoulli-number recursion [BCOR09], in the \(B_1 = -1/2\) convention:

\[\begin{split}\Omega_1(t) &= \int_0^t A(s)\, ds \\ \Omega_n(t) &= \sum_{j=1}^{n-1} \frac{B_j}{j!} \int_0^t S_n^{(j)}(s)\, ds ,\end{split}\]

with the \(S_n^{(j)}\) defined recursively in terms of the lower-order \(\Omega_m\):

\[\begin{split}S_n^{(1)} &= [\Omega_{n-1}, A] \\ S_n^{(j)} &= \sum_{m=1}^{n-j} [\Omega_m, S_{n-m}^{(j-1)}] , \qquad 2 \leq j \leq n-1 .\end{split}\]

Two facts about this recursion matter in practice.

Odd Bernoulli numbers vanish. \(B_j = 0\) for every odd \(j \geq 3\), so whole commutator groups drop out and only \(j = 1, 2, 4, 6, 8, \ldots\) contribute. The surviving coefficients are the ones the numerical core hard-codes:

\[\frac{B_1}{1!} = -\frac{1}{2} , \quad \frac{B_2}{2!} = \frac{1}{12} , \quad \frac{B_4}{4!} = -\frac{1}{720} , \quad \frac{B_6}{6!} = \frac{1}{30240} , \quad \frac{B_8}{8!} = -\frac{1}{1209600} .\]

Every term is a right-nested chain. Unrolling the recursion, each term of \(\Omega_n\) has the form

\[[\Omega_{m_1}, [\Omega_{m_2}, \ldots [\Omega_{m_j}, A] \ldots ]] , \qquad m_1 + m_2 + \cdots + m_j = n - 1 ,\]

carrying the coefficient \(B_j/j!\). The terms of the \(j\)-th group are therefore indexed by the compositions of \(n-1\) into \(j\) positive parts, of which there are \(\binom{n-2}{j-1}\). So

\[\#\,\Omega_n = \sum_{j\,:\,B_j \neq 0} \binom{n-2}{j-1} ,\]

which is what the implementation for orders 7 and above iterates over directly, rather than storing 129 expressions.

How many terms there are

import magnus.expansionterms as et

{n: et.count_terms(n) for n in range(1, 13)}

{1: 1,
 2: 1,
 3: 2,
 4: 3,
 5: 5,
 6: 9,
 7: 17,
 8: 33,
 9: 65,
 10: 129,
 11: 257,
 12: 513}

The count roughly doubles per order. This is why the implemented ceiling (magnus.magnus.MAGNUS_EXP_ORDER_MAX, currently 10) is a deliberate choice: the terms remain easy to generate far beyond it, but the work per slab grows with their number.

The expansion, written out

Orders 1 through 6, which the numerical core spells out inline:

et.print_magnus_terms(6)

Omega_1  (1 term)
    int A

Omega_2  (1 term)
    int -1/2 [Om_1, A]

Omega_3  (2 terms)
    int -1/2 [Om_2, A]
    int +1/12 [Om_1, [Om_1, A]]

Omega_4  (3 terms)
    int -1/2 [Om_3, A]
    int +1/12 [Om_1, [Om_2, A]]
    int +1/12 [Om_2, [Om_1, A]]

Omega_5  (5 terms)
    int -1/2 [Om_4, A]
    int +1/12 [Om_1, [Om_3, A]]
    int +1/12 [Om_2, [Om_2, A]]
    int +1/12 [Om_3, [Om_1, A]]
    int -1/720 [Om_1, [Om_1, [Om_1, [Om_1, A]]]]

Omega_6  (9 terms)
    int -1/2 [Om_5, A]
    int +1/12 [Om_1, [Om_4, A]]
    int +1/12 [Om_2, [Om_3, A]]
    int +1/12 [Om_3, [Om_2, A]]
    int +1/12 [Om_4, [Om_1, A]]
    int -1/720 [Om_1, [Om_1, [Om_1, [Om_2, A]]]]
    int -1/720 [Om_1, [Om_1, [Om_2, [Om_1, A]]]]
    int -1/720 [Om_1, [Om_2, [Om_1, [Om_1, A]]]]
    int -1/720 [Om_2, [Om_1, [Om_1, [Om_1, A]]]]

Orders 7 through 10 are generated from the recursion rather than typed in. Here is order 7 in full, and the sizes of the rest:

for term in et.omega_terms(7):
    print("   int " + et.format_term(term))

   int -1/2 [Om_6, A]
   int +1/12 [Om_1, [Om_5, A]]
   int +1/12 [Om_2, [Om_4, A]]
   int +1/12 [Om_3, [Om_3, A]]
   int +1/12 [Om_4, [Om_2, A]]
   int +1/12 [Om_5, [Om_1, A]]
   int -1/720 [Om_1, [Om_1, [Om_1, [Om_3, A]]]]
   int -1/720 [Om_1, [Om_1, [Om_2, [Om_2, A]]]]
   int -1/720 [Om_1, [Om_1, [Om_3, [Om_1, A]]]]
   int -1/720 [Om_1, [Om_2, [Om_1, [Om_2, A]]]]
   int -1/720 [Om_1, [Om_2, [Om_2, [Om_1, A]]]]
   int -1/720 [Om_1, [Om_3, [Om_1, [Om_1, A]]]]
   int -1/720 [Om_2, [Om_1, [Om_1, [Om_2, A]]]]
   int -1/720 [Om_2, [Om_1, [Om_2, [Om_1, A]]]]
   int -1/720 [Om_2, [Om_2, [Om_1, [Om_1, A]]]]
   int -1/720 [Om_3, [Om_1, [Om_1, [Om_1, A]]]]
   int +1/30240 [Om_1, [Om_1, [Om_1, [Om_1, [Om_1, [Om_1, A]]]]]]

{n: et.count_terms(n) for n in (8, 9, 10)}

{8: 33, 9: 65, 10: 129}

Going beyond the implemented ceiling

The generator has no ceiling of its own. Order 12, which the numerical core does not implement, is still inspectable:

terms_12 = et.omega_terms(12)
print(f"Omega_12 has {len(terms_12)} terms; the first three are")
for term in terms_12[:3]:
    print("   int " + et.format_term(term))

Omega_12 has 513 terms; the first three are
   int -1/2 [Om_11, A]
   int +1/12 [Om_1, [Om_10, A]]
   int +1/12 [Om_2, [Om_9, A]]

Coefficients come back as fractions.Fraction, so they are exact at any order and can be compared without a tolerance:

from fractions import Fraction

et.bernoulli_factor(8) == Fraction(-1, 1209600)

True

Checking the implementation against the derivation

The reason this module is worth having is that it shares no code with the numerical core, so the two agreeing is evidence rather than a tautology. The hard-coded group factors:

import magnus.magnus as mg

{'F1': (mg.F1, float(et.bernoulli_factor(2))),
 'F2': (mg.F2, float(et.bernoulli_factor(4))),
 'F3': (mg.F3, float(et.bernoulli_factor(6))),
 'F4': (mg.F4, float(et.bernoulli_factor(8)))}

{'F1': (0.08333333333333333, 0.08333333333333333),
 'F2': (-0.001388888888888889, -0.001388888888888889),
 'F3': (3.306878306878307e-05, 3.306878306878307e-05),
 'F4': (-8.267195767195768e-07, -8.267195767195768e-07)}

tests/test_expansionterms.py goes further: it evaluates the generated terms numerically on a sampled \(A(t)\) and compares them, order by order, against what magnus.magnus.magnus_expansion() produces internally. Agreement is at machine precision for every order from 1 to 10, which covers both the hand-written low orders and the generated high ones in a single check.

Choosing an order

Higher order buys a genuinely faster convergence rate in the slab width. Measured against a tight ODE ground truth on a smooth Hamiltonian, a single slab:

Magnus order

Observed convergence rate

Cost per slab, relative to order 6

4

\(\sim h^{6.6}\)

0.25

6

\(\sim h^{8.3}\)

1

7

\(\sim h^{8.0}\)

2.7

8

\(\sim h^{10.2}\)

5.0

10

measurement floor

16.9

Two caveats keep this from being a simple “higher is better”.

First, the cost grows faster than the order does, and narrowing the slabs at order 4 or 6 often reaches a given accuracy for less total work than raising the order. Magνs warns (magnus.magnus.MagnusHighOrderCostWarning) when an order above 6 is requested, for exactly this reason.

Second, and more fundamentally, the Magnus series converges only while \(\int \lVert A \rVert\, dt < \pi\). Beyond that radius no order helps, and the answer is narrower slabs – see Methodology and magnus.magnus.MagnusConvergenceWarning.

Note

Orders above 6 require integration_method='trapezoid' or 'simpson'. The Gauss-Legendre commutator-free schemes ('gl', the default) are separately derived integrators [BCR00], not products of this recursion, and exist only up to order 6; requesting more from them raises a ValueError rather than quietly returning an order-6 result.

API

See magnus.expansionterms in the API reference for the full signatures.