Exponential Algorithms

Geodesic is the exponential map, and evaluating a matrix exponential is a numerical problem rather than a formula one writes down. Since 0.2.0 the retraction therefore carries an algorithm saying how that exponential is evaluated, and there are five of them. They compute the same map. They differ in accuracy at a large step, in cost, and in which backends they run on.

This page addresses that numerical problem: why the exponential needs an algorithm at all, the two independent choices any such algorithm makes, what each of the five does, what they measure, and which one to reach for. The retractions themselves — what they are, and how Cayley and Geodesic differ as maps — are on the Retractions page, and the optimizer that uses them on Optimization on Homogeneous Spaces.

The exponential needs an algorithm

The notation this page uses

Everything below is stated in the notation of Both retractions factor the lift on the Retractions page, which is where it is derived. Briefly:

• $\bar{B} \in \mathbb{R}^{N\times{}N}$ is the horizontal lift, a skew-symmetric matrix of the same rank the dimension of the manifold (see StiefelLieAlgHorMatrix and GrassmannLieAlgHorMatrix). It is the direction an OptimizerMethod produces and the argument the retraction is applied to.
• $B'$ and $B''$ are the two $N\times{}2n$ thin factors of that lift, so that $\bar{B} = B'(B'')^T$. lift_factors returns them; for the Stiefel manifold they are the two matrices in

  \[\begin{bmatrix} A & -B^T \\ B & \mathbb{O} \end{bmatrix} = \underbrace{\begin{bmatrix} \tfrac{1}{2}A & \mathbb{I} \\ B & \mathbb{O} \end{bmatrix}}_{B'} \underbrace{\begin{bmatrix} \mathbb{I} & \mathbb{O} \\ \tfrac{1}{2}A & -B^T \end{bmatrix}}_{(B'')^T},\]

  and for a GrassmannLieAlgHorMatrix the same expression with $A \equiv \mathbb{O}$.
• $X := (B'')^TB' \in \mathbb{R}^{2n\times{}2n}$ is the reduced matrix. It is small even when $N$ is large, and it is the argument every matrix function on this page is evaluated on.

The reduction

The geodesic requires $\exp(\bar{B})$, but forming and exponentiating the full $N\times{}N$ lift would throw away its low-rank structure. Writing $\bar{B}=B'(B'')^T$ and $X=(B'')^TB'$ gives, for every $k\geq 1$,

\[\left(B'(B'')^T\right)^k = B'X^{k-1}(B'')^T.\]

Substitution into the exponential series therefore gives the exact identity

\[\exp\left(B'(B'')^T\right) = \mathbb{I} + B'\,\mathfrak{A}(X)\,(B'')^T, \qquad \mathfrak{A}(X) = \sum_{k=1}^\infty \frac{X^{k-1}}{k!},\]

so computing $\mathfrak{A}(X)$ is the central numerical task: once it is available, the $N\times{}N$ exponential is assembled using only the two thin factors, and the only matrix function evaluated anywhere is $\mathfrak{A}$ on a $2n\times{}2n$ argument.^[1] That is what preserves the cost advantage when $n\ll N$, and it is why this page is about $\mathfrak{A}$ throughout rather than about $\exp$.

There are two separate algorithmic choices to make in evaluating it:

1. how to approximate $\mathfrak{A}$ at a small argument — the kernel; and
2. whether to scale a large argument down and recover the original value afterwards — the recovery.

Taylor and Padé are kernels. Scaling and modified squaring is a framework around a kernel, not a competing approximation, so the two choices are independent. Classical dense matrix-exponential routines normally combine a Padé kernel with scaling and squaring [18, 19]; Taylor kernels can also be effective when scaling keeps the argument small [20].

All five algorithms in this package are:

Note that the two algorithms a reader is most likely to reach for, ScaledSquaring and NativePade, differ only in the kernel column. Both scale the argument down and both recover it with the same modified-squaring recurrence; the names are historical and unfortunate in this respect, since NativePade scales and squares exactly as much as ScaledSquaring does.

The Backend column says what the algorithm needs, not what it is known to run on: the first three require only matrix products and reductions, so whether they run on a particular accelerator depends on that backend's support for those operations, while the last two call into dense LAPACK and are therefore CPU-only.

Geodesic makes the choice visible:

Geodesic(ScaledSquaring())   # the default, and `Geodesic()`
Geodesic(NativePade())
Geodesic(AugmentedPade())    # the dense-CPU reference
Geodesic(ProjectedSkew())
Geodesic(TaylorSeries())     # the pre-0.2.0 behaviour; not a usable retraction

All five return the same exponential map, so the one-parameter subgroup property holds for every one of them. Two of the five are not ordinary choices, and both are documented here for a reason beyond being constructible:

• TaylorSeries is the small-argument kernel that ScaledSquaring uses. The section on it is not a deprecation notice but the case for scaling: it is what the kernel does when the argument is not small, and therefore what the framework of the next section exists to prevent.
• AugmentedPade is the route that hands $\mathfrak{A}$ to a native LAPACK exponential via one embedding. That route is available if wanted, and is strongly discouraged: it exponentiates a matrix four times the size needed, discards three quarters of the result, and gives up the backend portability of the direct algorithms for no gain in accuracy.

1. Direct Taylor series

TaylorSeries evaluates $\mathfrak{A}$ directly from its defining series, without scaling, and stops when a term falls below eps. The series converges for every matrix, but convergence alone does not guarantee an accurate floating-point sum.

Why the reduced argument is the hard case

The reduction of the previous section bought a small matrix, and it did so at a price: $X$ is a considerably worse argument for a power series than $\bar{B}$ is. Two facts about the pair say why, and both are short enough to derive here, because everything the rest of the page does about scaling follows from them.

$X$ and $\bar{B}$ have the same nonzero eigenvalues

This is the standard fact that $MK$ and $KM$ share their nonzero spectrum [22], applied to $\bar{B} = B'(B'')^T$ and $X = (B'')^TB'$. In one direction, if $\bar{B}v = \lambda{}v$ with $\lambda\neq0$, put $w := (B'')^Tv$; then

\[Xw = (B'')^TB'(B'')^Tv = (B'')^T\bar{B}v = \lambda\,(B'')^Tv = \lambda{}w,\]

and $w\neq0$ because $B'w = \bar{B}v = \lambda{}v \neq 0$. In the other, if $Xw = \lambda{}w$ with $\lambda\neq0$, put $v := B'w$; then $\bar{B}v = B'(B'')^TB'w = B'Xw = \lambda{}v$, and $v\neq0$ because $(B'')^Tv = Xw = \lambda{}w \neq 0$. The two constructions are mutually inverse, so the nonzero eigenvalues agree. Since $X$ is $2n\times{}2n$ and $\bar{B}$ has rank at most $2n$, that is the entire spectrum of both, apart from the zeros $\bar{B}$ carries on its null space.

That pins the spectral radius exactly. $\bar{B}$ is real skew-symmetric and therefore normal, $\bar{B}^T\bar{B} = -\bar{B}^2 = \bar{B}\bar{B}^T$, so its singular values are the moduli of its eigenvalues and its spectral radius is its two-norm. Writing $\rho$ for the spectral radius, $\rho(X) := \max_i|\lambda_i(X)|$, the paragraph above gives an equality rather than an estimate:

\[\rho(X) = \rho(\bar{B}) = \|\bar{B}\|_2.\]

The norm of $X$ is quadratic in the norm of $\bar{B}$

Multiplying out the two thin factors of The notation this page uses gives $X$ in closed form:

\[X = (B'')^TB' = \begin{bmatrix} \mathbb{I} & \mathbb{O} \\ \tfrac{1}{2}A & -B^T \end{bmatrix} \begin{bmatrix} \tfrac{1}{2}A & \mathbb{I} \\ B & \mathbb{O} \end{bmatrix} = \begin{bmatrix} \tfrac{1}{2}A & \mathbb{I} \\ \tfrac{1}{4}A^2 - B^TB & \tfrac{1}{2}A \end{bmatrix}.\]

Every block is linear in the entries of the lift except the lower-left one, which is quadratic. Replacing $(A, B)$ by $(tA, tB)$ therefore multiplies $\bar{B}$ by $t$ and that block by $t^2$, so $\|X\|$ grows like $\|\bar{B}\|_2^2$. The constant depends on how the lift's weight is split between its two blocks, and the two pure cases bracket it:

• with $B = \mathbb{O}$ the block is $\tfrac{1}{4}A^2$ and $\|\bar{B}\|_2 = \|A\|_2$, so $\|X\|_2 \to \tfrac{1}{4}\|\bar{B}\|_2^2$;
• with $A = \mathbb{O}$ — which is every GrassmannLieAlgHorMatrix — it is $-B^TB$ and $\|\bar{B}\|_2 = \|B\|_2$, so $\|X\|_2 \to \|\bar{B}\|_2^2$.

A lift with a full $B$ block sits at the top of that range, $B$ being $(N-n)\times{}n$ against $A$'s $n\times{}n$. Both claims are cheap to check, so the sweep checks them — $\rho(X)$ against $\|\bar{B}\|_2$ column by column, the eigenvalues themselves against each other, and the quadratic ratio:

# The 2n eigenvalues of B̄ of largest modulus, in the same order as `eigvals(X)`: the rest are the
# zeros on its null space.
byᵢₘ(λs) = sort(λs, by = λ -> (imag(λ), real(λ)))
top(λs, k) = byᵢₘ(λs[partialsortperm(abs.(λs), 1:k, rev = true)])

rows = map(lifts) do B
    B′, B′′ = lift_factors(B)
    X = B′′' * B′
    λX = byᵢₘ(eigvals(X))
    λB = top(eigvals(Matrix(B)), length(λX))
    nrm = opnorm(Matrix(B), 2)
    [fixed(nrm), fixed(maximum(abs, λX)), sci(maximum(abs, λX - λB)),
     sci(opnorm₁(X)), fixed(opnorm(X, 2) / nrm^2)]
end

table(["‖B̄‖₂", "ρ(X)", "‖λ(X) − λ(B̄)‖∞", "‖X‖₁", "‖X‖₂/‖B̄‖₂²"], rows)

The first two columns are the same number to every digit printed — the equality above — and the third says the whole spectrum agrees and not just its largest modulus. The last column is the quadratic constant, and from the second row on it sits at the $B$-dominant end of $[\tfrac{1}{4}, 1]$, as it should for a lift whose $B$ block is the larger one. The first row is the one place the asymptotic statement does not apply: $X$ carries an identity block, so $\|X\| \geq 1$ however small the lift is, and a quadratic in $\|\bar{B}\|_2$ cannot describe a lift of norm $0.44$. That regime is harmless — it is the regime in which the series needs no help at all.

Why that combination is bad for a power series

The two facts together are the difficulty. $X$ carries the spectrum of a normal matrix while its norm is the square of that spectrum's radius, so

\[\frac{\rho(X)}{\|X\|_2} \approx \frac{1}{\|\bar{B}\|_2} \longrightarrow 0 \qquad\text{as}\qquad \|\bar{B}\|_2 \to \infty:\]

the reduction takes a normal matrix and returns one that is as far from normal as its size allows. That is the sense in which this page and the docstrings call $X$ "strongly non-normal". The label matters less than its consequence, which is the subject of [23]: for a non-normal matrix the transient behaviour of powers and of the exponential is governed by the norm while the asymptotics are governed by the spectrum, and the gap between the two can be arbitrarily large.

A truncated power series is precisely what that gap breaks. The kernel forms the partial sums

\[S_m := \sum_{k=1}^{m}\frac{X^{k-1}}{k!}, \qquad m = 1, 2, \ldots, M,\]

and returns $S_M$, where $M$ is the first index at which a term falls below eps. The individual terms are sized by $\|X\|$, while the limit $\mathfrak{A}(X)$ is of moderate size because it is sized by the spectrum. The sequence $S_1, S_2, \ldots, S_M$ therefore rises to a peak and cancels back down to its limit. Rounding error is committed at the scale of the largest partial sum and survives into an answer of size $\|\mathfrak{A}(X)\|$, so the decimal digits lost to that cancellation are about

\[\text{digits lost} \approx \log_{10} \frac{\max_{1\leq m\leq M}\|S_m\|}{\|\mathfrak{A}(X)\|}.\]

The maximum is over $m$: it is the largest of the $M$ partial sums the summation actually forms — a peak along a finite sequence, not a limit. The sequence $S_m$ converges, and the ratio measures how far above its own limit it travels on the way. Where the norm and the spectral radius are comparable — a normal argument, or any small one — the terms decrease from the first, the largest partial sum is $S_M$ itself, and the ratio is one: nothing cancels and nothing is lost. The denominator is the true $\mathfrak{A}(X)$, taken from AugmentedPade, rather than the computed $S_M$, which is the quantity under suspicion.

Summing the series exactly as 𝔄 does — the same recurrence, the same termination test — and recording the peak on the way:

"""
The peak of the partial sums against the value they converge to. The recurrence divides before it
multiplies, as the kernel's does, so no *term* overflows where the sum itself does not.
"""
function cancellation(X)
    term = unit_matrix(X)               # k = 1, i.e. X⁰/1!
    S = copy(term)
    peak = norm(S)
    k = 1
    while norm(term) > eps(eltype(X))
        k += 1
        term = term * X / k             # Xᵏ⁻¹/k!
        S = S + term
        peak = max(peak, norm(S))
    end
    (; terms = k, peak = peak, computed = norm(S))
end

rows = map(lifts) do B
    B′, B′′ = lift_factors(B)
    X = B′′' * B′
    c = cancellation(X)
    reference = 𝔄(X, AugmentedPade())
    [fixed(norm(Matrix(B))), string(c.terms), sci(c.peak), fixed(norm(reference)),
     fixed(log10(c.peak / norm(reference))),
     sci(norm(𝔄(X, TaylorSeries()) - reference) / norm(reference))]
end

table(["‖B̄‖", "terms", "peak ‖Sₘ‖", "‖𝔄(X)‖", "digits lost", "relative error"], rows)

The estimate and the measurement track each other: at the fourth lift about eight and a half digits are predicted lost and the computed series is wrong in the eighth, while by the fifth the predicted loss passes the sixteen digits Float64 has and the relative error is above one — no correct digits at all, rather than a lost digit or two. The number of terms is the number of chances to accumulate that loss, and it grows with the lift. The termination test cannot see any of it: a term below eps bounds the truncation error, and what has gone wrong is the rounding error, which is of size $\varepsilon\cdot\max_m\|S_m\|$ where the test is calibrated for $\varepsilon\cdot\|\mathfrak{A}(X)\|$. That the two coincide at a small argument is the entire reason the kernel is usable there and nowhere else.

That the direct series is not a method for the matrix exponential has been known since [24] — restated with twenty-five more years of evidence in [25]. What is specific here is the direction of the effect: the factorisation that made the problem cheap also made its argument worse than the matrix one started with, so this package meets the failure earlier than a dense implementation would.

What that forward error does to the retraction is the number an optimizer actually cares about. check(geodesic(B, TaylorSeries())), i.e. $\|Y^TY - \mathbb{I}\|$ of the retracted point, over the same sweep:

table(["‖B̄‖", "`check`"],
      [[fixed(norm(Matrix(B))), sci(check(geodesic(B, TaylorSeries())))] for B in lifts])

By $\|\bar{B}\| \approx 79$ the "retracted" point is not on the Stiefel manifold in any sense, and at the top of the sweep the series has overflowed outright. The Cayley and ScaledSquaring columns of the same sweep, for comparison, are in Staying on the manifold below and stay at round-off throughout.

Two things about the failure are worth recording rather than merely deprecating. First, it is silent: an optimizer using it takes a step, gets a matrix back, and nothing anywhere reports that the matrix is no longer on the manifold. No optimizer in this package currently verifies that a retracted point is still on its manifold, which is a gap rather than a decision — see #76. Second, the obvious first fix does not work: making the termination test relative to the partial sum rather than absolute was measured to change none of the numbers above, at any lift norm, for the reason the displayed ratio gives. The loss is in the cancellation inside the sum, not in the point at which the summation stops. Making the argument small is what helps, and that is the next section.

2. Scaling and modified squaring

The previous section blames the size of the argument, so the remedy is to evaluate the kernel somewhere else. Scaling and squaring computes a matrix function at $X/2^s$ and then reconstructs its value at $X$ — for $\exp$ by repeated squaring, $\exp(X) = (\exp(X/2^s))^{2^s}$. For $\mathfrak{A}$ the recovery formula is not squaring but a modified squaring [20], and it is worth deriving because it is the one step of these algorithms that has no counterpart in the classical dense routine.

The doubling identity for $\mathfrak{A}$

Everything follows from one identity. From $\exp(2Y) - \mathbb{I} = (\exp(Y) - \mathbb{I})(\exp(Y) + \mathbb{I})$ and $\exp(Y) = \mathbb{I} + Y\mathfrak{A}(Y)$,

\[\mathfrak{A}(2Y) = \frac{\exp(2Y) - \mathbb{I}}{2Y} = \mathfrak{A}(Y)\,\frac{\exp(Y) + \mathbb{I}}{2} = \mathfrak{A}(Y)\,\frac{2\mathbb{I} + Y\mathfrak{A}(Y)}{2},\]

that is

\[\mathfrak{A}(2Y) = \mathfrak{A}(Y) + \tfrac{1}{2}Y\mathfrak{A}(Y)^2.\]

The middle step divides by $Y$, but the two ends do not: both sides are power series in $Y$, so the identity holds as an identity of power series and therefore at every $Y$, singular or not. It is an exact algebraic identity, not an approximation — the $\mathfrak{A}$ case of the modified squaring of [20].

Now put

\[W_k := \frac{\mathfrak{A}(X/2^k)}{2^k},\]

which is the object both algorithms below actually carry: not $\mathfrak{A}$ of anything, but $\mathfrak{A}$ at the $k$-times-halved argument, divided by $2^k$. That scaling is what makes the recovery step clean. Applying the doubling identity at $Y = X/2^k$ gives

\[W_{k-1} = \frac{\mathfrak{A}(X/2^{k-1})}{2^{k-1}} = \frac{1}{2^{k-1}}\left(\mathfrak{A}(Y) + \tfrac{1}{2}Y\mathfrak{A}(Y)^2\right) = 2W_k + W_kXW_k,\]

where the last equality uses $Y = X/2^k$ to turn $Y\mathfrak{A}(Y)^2/2^k$ into $W_kXW_k$. So one recovery step is one application of

\[W \longleftarrow 2W + WXW,\]

at $2n\times{}2n$, with the original $X$ and not the scaled one. Starting from $W_s$ and applying it $s$ times lands on $W_0 = \mathfrak{A}(X)$, which is what was wanted.^[2]

The framework, and the two algorithms in it

Scaling and modified squaring is a framework: it says nothing about how $W_s$ is obtained, only how to get from $W_s$ to $W_0$. Filling in a kernel gives an algorithm, and this package has two:

They differ in one column — the kernel. In particular NativePade scales and squares too, with the same rule for $s$ and the same recurrence; the algorithm names suggest a contrast that does not exist.

$\theta$ is also the same quantity in both, and means the same thing: the ceiling the scaled argument is brought below, entering each algorithm through the one formula in the scaling column and nowhere else. Only the admissible range differs, and it does so because of the kernel rather than because of $\theta$: the Taylor series sums until its terms vanish, so it copes with whatever it is handed, whereas the Padé kernel does a fixed number of steps that were sized for $\|Y\|_1\leq1/2$ and is inaccurate above it. Both constructors default to 0.5, so the same $\theta$ and the same $s$ are what the comparisons on this page hold fixed; see What a large $\theta$ costs NativePade for the measurement behind the restriction.

The names are a historical record rather than a taxonomy, and they are worth reading carefully because of what they do not offer. ScaledSquaring means scaling with a Taylor kernel and NativePade means scaling with a Padé kernel: the two axes are welded together, so there is no way to ask for the framework with a kernel of your choosing, and in particular ScaledSquaring has no Padé support of any kind. Separating the two — a kernel argument to one scaling-and-squaring algorithm, with Taylor the default — is the natural shape for this code and is #63, along with the question of which kernel is the better choice, which this page does not settle.

What scaling buys the Taylor kernel is exactly what 1. Direct Taylor series measured, and the same cancellation function measures it — this time on $X/2^s$, which is the argument the kernel is actually handed:

θ = 0.5   # the default of both algorithms

rows = map(lifts) do B
    B′, B′′ = lift_factors(B)
    X = B′′' * B′
    s = max(0, ceil(Int, log2(opnorm₁(X) / θ)))
    c = cancellation(X / 2^s)
    reference = 𝔄(X, AugmentedPade())
    [fixed(norm(Matrix(B))), string(s), fixed(opnorm₁(X / 2^s)), string(c.terms),
     fixed(c.peak / c.computed),
     sci(norm(𝔄(X, ScaledSquaring()) - reference) / norm(reference))]
end

table(["‖B̄‖", "s", "‖X/2ˢ‖₁", "terms", "cancellation ratio", "relative error"], rows)

Every row is summed below $\theta$, where the terms decrease from the first one: the cancellation ratio is one to four decimals instead of the $10^{20}$ of the unscaled table, six to eleven terms suffice instead of hundreds, and the relative error of the whole algorithm stays at round-off across a sweep on which the bare kernel loses every digit it has. In this way the framework makes an otherwise unusable kernel usable.

The complete ScaledSquaring algorithm

Given the reduced matrix $X = (B'')^TB'$ and the threshold $\theta$:

1. Choose the scaling. Set $s = \max(0, \lceil\log_2(\|X\|_1/\theta)\rceil)$, so that $\|X/2^s\|_1 \leq \theta$.
2. Evaluate the kernel. Sum the Taylor series of 1. Direct Taylor series at $X/2^s$ and divide by $2^s$, giving $W_s = \mathfrak{A}(X/2^s)/2^s$.
3. Recover. Apply $W \leftarrow 2W + WXW$ exactly $s$ times, with the original $X$.
4. Return $W_0 = \mathfrak{A}(X)$.

The division by $2^s$ in step 2 is part of the definition of $W_s$ and not an approximation, and neither is step 3. Every operation is $2n\times{}2n$; nothing $N\times{}N$ is squared, inverted or exponentiated anywhere.^[2]

The threshold $\theta$

θ defaults to 0.5. A smaller value performs more scaling steps; a larger value asks the kernel to handle a larger argument. That default is not taken from the literature: the $\theta_m$ tables of [18, 19] are derived for $\exp$ rather than for $\mathfrak{A}$, and they bound a backward error in $\|X\|$ — the least informative norm available here, for the reason Why the reduced argument is the hard case gives. How 0.5 was in fact arrived at differs between the two algorithms:

• For NativePade it is forced. Its five Newton–Schulz steps are only accurate while $\|I - q_6(Y)\|_1 < 1$, and $\theta \leq 1/2$ is what guarantees that with room to spare. See Applying the denominator without a matrix solve for the bound and What a large $\theta$ costs NativePade for the measured failures above it.
• For ScaledSquaring it is an empirical choice, and a nearly free one, because the algorithm turns out to be insensitive to it. The measurement is in Sensitivity to the threshold θ below, recomputed at build time rather than quoted here: over the 32-fold range $\theta \in [0.125, 4]$ both check and the forward error move by less than a factor of six, and neither is monotone in $\theta$. Nothing in that table singles out 0.5; it is a reasonable value rather than a tuned one, and it is one lift and one seed rather than a general error bound.

Advantages. It is the cheapest of the algorithms on the $\mathfrak{A}$ call itself, and it comes as close to exp(Matrix(B)) as AugmentedPade does, which is as close as anything here gets; the two measurements are What they cost and Agreeing with the exponential below. It needs only matrix products, norms and a kernel-written identity, so it has no dense-LAPACK dependency. It is the default because it is the cheapest of the algorithms that have none.^[3]

Disadvantages. Its orthogonality is the outcome of an arithmetic cancellation rather than a structural property, so check drifts upwards with the size of the lift — by two to three orders of magnitude over Staying on the manifold, and considerably further in Float32. Only ProjectedSkew avoids that drift. And it takes about twice the squarings it needs, for the reason in the note below.

3. Padé approximation

Scaling is settled. The previous section made the argument small and gave an exact recovery for the original one, and it did so without committing to any particular kernel. What it did not settle is which approximation to evaluate at the small argument. ScaledSquaring evaluates a Taylor series there because that is what it already had; NativePade evaluates a rational function instead, and this section is why that is a different and defensible choice.

A degree-$m$ Taylor polynomial simply keeps the first $m+1$ terms,

\[T_m(z)=\sum_{k=0}^m c_kz^k, \qquad f(z)-T_m(z)=O(z^{m+1}).\]

Its coefficients are fixed one by one by the series, so matching more terms requires increasing the polynomial degree. Padé takes a different route: it chooses a numerator and a denominator together, so that their quotient matches more series coefficients without raising the numerator to the same degree. It is therefore an alternative small-argument kernel, not an alternative to scaling.

To see how this works, first consider a scalar analytic function $f$. Its $[m/n]$ Padé approximant is the rational function

\[R_{m,n}(z) = \frac{P_m(z)}{Q_n(z)}, \qquad \deg P_m\leq m,\quad \deg Q_n\leq n,\quad Q_n(0)=1,\]

and is defined by the matching condition

\[Q_n(z)f(z)-P_m(z)=O\left(z^{m+n+1}\right).\]

Equating powers of $z$ determines the coefficients of $P_m$ and $Q_n$. The denominator is what makes this different from truncating a Taylor series: when the quotient is expanded again as a power series, even low-degree $P_m$ and $Q_n$ generate infinitely many powers. For example,

\[\frac{1+z/2}{1-z/2}=1+z+\frac{z^2}{2}+\frac{z^3}{4}+\cdots\]

matches $e^z$ through degree two, whereas a polynomial with numerator degree one can match only $1+z$. More generally, a $[m/n]$ Padé approximant matches through degree $m+n$ whenever its matching condition is solvable — for $e^z$ it always is, as the next section shows. This is why a rational kernel can capture much more local information than a Taylor polynomial with a similar polynomial degree [18, 26]. The scalar variable $z$ serves only to fix the coefficients; a matrix takes its place once they are known.

Deriving the coefficients for $\mathfrak{A}$

Nothing is fitted to $\mathfrak{A}$, and nothing is fitted numerically. The coefficients come from the matching condition above, applied to $f=\exp$ at $m=7$ and $n=6$; one subtraction then converts that approximant into one for $\mathfrak{A}$. Write its numerator and denominator as $P^{\exp}_7$ and $Q^{\exp}_6$:

\[e^z = \frac{P^{\exp}_7(z)}{Q^{\exp}_6(z)} + O(z^{14}), \qquad P^{\exp}_7(0)=Q^{\exp}_6(0)=1.\]

Keep $m$ and $n$ general for the derivation, since one formula covers every degree:

\[P_m(z)=\sum_{k=0}^m a_kz^k, \qquad Q_n(z)=\sum_{j=0}^n b_jz^j, \qquad b_0=1.\]

Because $e^z=\sum_{r\geq0}z^r/r!$, the $z^k$ coefficient of $Q_n(z)e^z$ is the convolution at degree $k$,

\[c_k := \sum_{j=0}^{n}\frac{b_j}{(k-j)!},\]

reading $1/(k-j)!$ as zero when $j>k$. That single quantity is what the rest of the derivation is about. The matching condition $Q_ne^z-P_m=O(z^{m+n+1})$ says that the $z^k$ coefficient of $Q_ne^z-P_m$ vanishes for every $k\leq m+n$, and since $P_m$ has no term above degree $m$, it says two different things on the two ranges of $k$:

\[a_k=c_k \quad (k=0,\ldots,m), \qquad c_k=0 \quad (k=m+1,\ldots,m+n).\]

Only the right-hand block is a system to be solved: $n$ equations in the $n$ unknowns $b_1,\ldots,b_n$. The left-hand block constrains nothing, since each $a_k$ occurs in one equation and nowhere else — with the denominator known it reads the numerator off. The construction is a statement about the denominator alone. For $m=7$ and $n=6$ the ranges are the degrees $0,\ldots,7$ and $8,\ldots,13$.

The solution is

\[a_k=\frac{(m+n-k)!}{(m+n)!}\binom{m}{k}, \qquad b_k=(-1)^k\frac{(m+n-k)!}{(m+n)!}\binom{n}{k},\]

and one binomial identity settles both blocks at once. Substituting this $b_j$ into $c_k$, clearing $(m+n)!$, and writing $d := m+n-k$ turns the convolution into $(m+n)!\,c_k = d!\sum_{j}(-1)^j\binom{n}{j}\binom{m+n-j}{d}$, using $(m+n-j)!/(k-j)! = d!\binom{m+n-j}{d}$ — a product of $d$ consecutive integers, since the exponent difference does not depend on $j$. What remains is the claim

\[\sum_{j=0}^n(-1)^j\binom{n}{j}\binom{m+n-j}{d}=\binom{m}{d-n},\]

which is the coefficient of $x^d$ in $\sum_j(-1)^j\binom{n}{j}(1+x)^{m+n-j} = (1+x)^m\bigl((1+x)-1\bigr)^n = x^n(1+x)^m$, where $x$ is a bookkeeping variable unrelated to the $z$ of the approximant. The claim splits exactly where the matching condition splits:

• if $k\geq m+1$, then $d-n=m-k<0$, and $x^n(1+x)^m$ has no $x^d$ term at all, so $c_k=0$ — all $n$ denominator equations hold;
• if $k\leq m$, the coefficient is $\binom{m}{m-k}=\binom{m}{k}$, and restoring the $d!/(m+n)!$ gives $c_k = \frac{(m+n-k)!}{(m+n)!}\binom{m}{k}$, which is the claimed $a_k$.

The two formulas are one expression with $m$ and $n$ interchanged and a sign, $a_k(m,n)=(-1)^kb_k(n,m)$: the reflection that $e^{-z}=1/e^z$ induces on the table. They are classical [18, 26]; the derivation appears here so the numbers below are traceable rather than quoted. At $m=7$ and $n=6$:

\[\begin{array}{r|ccc} k & a_k & b_k & a_k-b_k\\\hline 0 & 1 & 1 & 0\\ 1 & \tfrac{7}{13} & -\tfrac{6}{13} & 1\\ 2 & \tfrac{7}{52} & \tfrac{5}{52} & \tfrac{1}{26}\\ 3 & \tfrac{35}{1716} & -\tfrac{5}{429} & \tfrac{5}{156}\\ 4 & \tfrac{7}{3432} & \tfrac{1}{1144} & \tfrac{1}{858}\\ 5 & \tfrac{7}{51480} & -\tfrac{1}{25740} & \tfrac{1}{5720}\\ 6 & \tfrac{7}{1235520} & \tfrac{1}{1235520} & \tfrac{1}{205920}\\ 7 & \tfrac{1}{8648640} & 0 & \tfrac{1}{8648640} \end{array}\]

The last column is the promised subtraction. Since $\mathfrak{A}(z)=(e^z-1)/z$,

\[\begin{aligned} \mathfrak{A}(z) &= \frac{e^z-1}{z}\\ &= \frac{P^{\exp}_7(z)-Q^{\exp}_6(z)}{z\,Q^{\exp}_6(z)} + O(z^{13})\\ &= \frac{p_6(z)}{q_6(z)} + O(z^{13}), \end{aligned}\]

where

\[p_6(z)=\frac{P^{\exp}_7(z)-Q^{\exp}_6(z)}{z}, \qquad q_6(z)=Q^{\exp}_6(z).\]

The $k=0$ entry of that column vanishes, both constant terms being one, so the difference is divisible by $z$ and $p_6$ is a polynomial. Its coefficients are the remaining entries shifted down one degree, $(p_6)_k=a_{k+1}-b_{k+1}$ for $k=0,\ldots,6$ with $b_7=0$. Numerator and denominator both have degree six, so this is a $[6/6]$ approximant of $\mathfrak{A}$, and it agrees with

\[\mathfrak{A}(z)=1+\frac{z}{2!}+\frac{z^2}{3!}+\cdots\]

through the $z^{12}$ term, where truncating that series after $z^6$ agrees only through $z^6$. The first term missed is small: in exact arithmetic

\[q_6(z)\mathfrak{A}(z)-p_6(z)=\frac{z^{13}}{149597947699200}+O(z^{14}),\]

so at the scaled argument NativePade evaluates, $|z|\leq1/2$, that term is about $8\cdot10^{-19}$ — the same order as the inverse-iteration residual derived below, and far under the Float64 rounding of a result of size one. The polynomials the implementation uses are therefore

\[\begin{aligned} p_6(z)={}&1+\frac{z}{26}+\frac{5z^2}{156}+\frac{z^3}{858} +\frac{z^4}{5720}+\frac{z^5}{205920}+\frac{z^6}{8648640},\\ q_6(z)={}&1-\frac{6z}{13}+\frac{5z^2}{52}-\frac{5z^3}{429} +\frac{z^4}{1144}-\frac{z^5}{25740}+\frac{z^6}{1235520}. \end{aligned}\]

From scalar Padé to matrix Padé

The scalar calculation above has done its only job: it has determined the coefficients. Matrix functions defined by power series use those same scalar coefficients, with each scalar power $z^k$ replaced by the matrix power $Y^k$. Thus $p_6(z)$ and $q_6(z)$ become the matrix polynomials $p_6(Y)$ and $q_6(Y)$ displayed above with $z$ replaced by $Y$ and $1$ replaced by $I$. The matching statement also transfers directly: the expansion of the rational matrix function agrees with

\[\mathfrak{A}(Y)=I+\frac{Y}{2!}+\frac{Y^2}{3!}+\cdots\]

through the $Y^{12}$ term.

The quotient is not taken entry by entry. For scalars, $r=p/q$ means $qr=p$. For matrices, the corresponding operation is therefore to solve

\[q_6(Y)R=p_6(Y),\]

or, when $q_6(Y)$ is invertible, to write

\[R=q_6(Y)^{-1}p_6(Y)\approx\mathfrak{A}(Y).\]

Because $p_6(Y)$ and $q_6(Y)$ are polynomials in the same matrix, they commute; the order shown is the one used by the implementation. Notice that no inverse of $Y$ occurs. Unlike the expression $(e^Y-I)Y^{-1}$, the Padé formula remains directly usable when $Y$ is singular. In _native_pade_polynomials, the code first forms $Y^2$ and $Y^4$ and groups the displayed polynomials around those shared powers. This evaluates both degree-6 polynomials with matrix products instead of constructing the powers one by one.

Applying the denominator without a matrix solve

A conventional evaluation of this rational matrix function would solve

\[q_6(Y)W=p_6(Y).\]

To derive the iteration [27], analysed in [26, §7.2], temporarily write $A=q_6(Y)$ and seek a matrix $Z$ satisfying $Z^{-1}-A=0$. The Fréchet derivative of matrix inversion is

\[D(Z^{-1})[H]=-Z^{-1}HZ^{-1}.\]

One Newton correction therefore solves

\[-Z_j^{-1}H_jZ_j^{-1}=A-Z_j^{-1},\]

which gives $H_j=Z_j-Z_jAZ_j$ and hence

\[Z_0=I, \qquad Z_{j+1}=Z_j+H_j=Z_j\left(2I-AZ_j\right).\]

This rearrangement is Newton's method using no factorization or linear solve. It is not safe from an arbitrary starting point. Writing $E_j := I - AZ_j$ for the inverse residual at step $j$, the condition is that the initial one be a contraction in some submultiplicative norm:

\[\|E_0\| = \|I - AZ_0\| < 1.\]

For the present choice $A=q_6(Y)$ and $Z_0=I$ that reads $\|I-q_6(Y)\|<1$, which is exactly what scaling will deliver below. Since $AZ_j=I-E_j$ and $2I-AZ_j=I+E_j$, the update gives

\[AZ_{j+1}=(I-E_j)(I+E_j)=I-E_j^2, \qquad E_{j+1}=E_j^2.\]

Consequently $E_j=E_0^{2^j}$ in exact arithmetic and $\|E_j\|\leq\|E_0\|^{2^j}$: once $\|E_0\|<1$, the iteration converges quadratically, squaring the residual at every step — the ordinary local behaviour of Newton's method [26, §7.2], with $\|E_0\|<1$ identifying the basin rather than supplying a contraction constant. The code writes the first step explicitly as $Z_1=2I-q_6(Y)$ and performs four more, giving $E_5=(I-q_6(Y))^{32}$. Scaling ensures $\|Y\|_1\leq 1/2$; from the displayed coefficients,

\[\|I-q_6(Y)\|_1 \leq \sum_{k=1}^6 |(q_6)_k|\,\|Y\|_1^k \leq \sum_{k=1}^6 \frac{|(q_6)_k|}{2^k} = 0.2563204283419127\ldots <0.257,\]

so $\|E_5\|_1 \leq 0.257^{32}<1.3\cdot10^{-19}$. This explains both the fixed five Newton-Schulz steps and the constructor restriction $0<\theta\leq 1/2$: together they make the solve-free inverse accurate to approximately Float64 precision using matrix multiplication alone. It is also the same order as the $8\cdot10^{-19}$ kernel truncation error derived above, so neither half of the accuracy argument is the weaker one.

The complete NativePade algorithm

Given the reduced matrix $X=(B'')^TB'$:

1. Choose $s=\max(0,\lceil\log_2(\|X\|_1/\theta)\rceil)$ and set $\alpha=2^s$.

2. Set $Y=X/\alpha$, so $\|Y\|_1\leq\theta\leq1/2$, and evaluate $p_6(Y)$ and $q_6(Y)$.

3. Compute $Z_5\approx q_6(Y)^{-1}$ with the five Newton-Schulz steps above.

4. Form $W_s=Z_5p_6(Y)/\alpha$. The division by $\alpha$ is required because

   \[\exp\left(B'(B'')^T/\alpha\right) =I+B'\left[\mathfrak{A}(X/\alpha)/\alpha\right](B'')^T.\]

5. Apply $W\leftarrow2W+WXW$ exactly $s$ times — the low-rank modified-squaring steps of the previous section, unchanged and for the same reason. After the last one, $W\approx\mathfrak{A}(X)$.

Padé is therefore the small-argument kernel in NativePade, not a replacement for scaling and squaring. Steps 1 and 5 are shared with ScaledSquaring verbatim; steps 2 to 4 are the only difference between the two algorithms. Scaling makes the rational approximation and its solve-free denominator application reliable; modified squaring transports that small-argument result back to the original $X$.

What scaling buys the Padé kernel

1. Direct Taylor series is what the Taylor kernel does at an argument that is not small. The same measurement for the Padé kernel, over the same eight lifts, makes the point that the two kernels are in the same position: neither is usable on its own.

Two unscaled variants are worth separating, because an unscaled Padé kernel fails in two independent ways. The first evaluates $q_6(X)^{-1}p_6(X)$ with an exact dense solve, so its error is the approximation error of the $[6/6]$ rational function alone — the best the kernel could possibly do without scaling, and not something NativePade could run on a GPU backend anyway. The second is what NativePade would actually compute with $s=0$: the same rational function with the five Newton–Schulz steps, which lose their residual bound the moment $\|Y\|_1$ exceeds $1/2$. Relative error against AugmentedPade:

# The kernel of the previous subsections on its own, with an exact denominator solve — the best an
# unscaled `[6/6]` can do. `\` is a dense LAPACK solve and is used only for this comparison.
function pade_exact_solve(X)
    p, q = _native_pade_polynomials(X, unit_matrix(X))
    q \ p
end

# The same kernel with the denominator applied the way `NativePade` applies it.
function pade_newton_schulz(X)
    𝕀 = unit_matrix(X)
    p, q = _native_pade_polynomials(X, 𝕀)
    q⁻¹ = 2 * 𝕀 - q
    for _ in 1:4
        q⁻¹ = q⁻¹ * (2 * 𝕀 - q * q⁻¹)
    end
    q⁻¹ * p
end

rows = map(lifts) do B
    B′, B′′ = lift_factors(B)
    X = B′′' * B′
    reference = 𝔄(X, AugmentedPade())
    relative(W) = norm(W - reference) / norm(reference)
    [fixed(norm(Matrix(B))), fixed(opnorm₁(X)),
     sci(relative(pade_exact_solve(X))), sci(relative(pade_newton_schulz(X))),
     sci(relative(𝔄(X, NativePade())))]
end

table(["‖B̄‖", "‖X‖₁", "kernel, exact solve", "kernel, Newton–Schulz", "`NativePade`"], rows)

The kernel column with the exact solve is accurate at the first lift, has lost half its digits by $\|X\|_1 \approx 16$, and is wrong by $O(1)$ — a relative error above one, meaning no correct digits at all — from $\|X\|_1 \approx 600$ upwards. It does not overflow the way the Taylor kernel does, because a rational function stays bounded where a truncated series does not; it simply approximates a different function. The Newton–Schulz column is worse and worse for the separate reason of the previous subsection: past $\|Y\|_1 = 1/2$ the residual no longer contracts, so five squarings of a quantity larger than one amplify rather than converge.

The last column is the same kernel with steps 1 and 5 of the algorithm around it, and it holds $10^{-14}$ across the sweep. What the framework supplies is the guarantee that the kernel is only ever asked for $\|Y\|_1 \leq \theta$, where both failure modes are absent.

What a large $\theta$ costs NativePade

$\theta$ is the same quantity in the two algorithms — the same ceiling on $\|X/2^s\|_1$, reached through the same formula for $s$, as The framework, and the two algorithms in it sets out. What differs is how much of it the kernel tolerates, and that is the reason NativePade(θ) rejects $\theta > 1/2$ where ScaledSquaring(θ) accepts any positive value. ScaledSquaring sums its series until the terms vanish, so a larger $\theta$ only asks it to sum more terms; Sensitivity to the threshold θ shows what that costs, which is almost nothing over a 32-fold range. NativePade does a fixed five Newton–Schulz steps against a fixed rational approximant, and both were sized for $\|Y\|_1 \leq 1/2$. Above that, the iteration count is simply too small, and no quantity the algorithm computes says so.

Worst relative error of the kernel against AugmentedPade over 400 random $6\times6$ arguments of one-norm exactly $\theta$ — that is, the kernel asked for exactly what a threshold of $\theta$ would hand it:

Random.seed!(2024)  # nothing else on this page draws after `lifts`, so this is self-contained

worst(θ) = maximum(1:400) do _
    Y = randn(6, 6)
    Y = θ * Y / opnorm₁(Y)
    reference = 𝔄(Y, AugmentedPade())
    norm(pade_newton_schulz(Y) - reference) / norm(reference)
end

θs = [0.5, 1.0, 1.5, 2.0, 3.0]
table(vcat("``\\theta``", string.(θs)), [vcat("worst relative error", sci.(worst.(θs)))])

Only the pattern is reproducible, not the digits: these are maxima over a random draw, and the two rightmost entries move by orders of magnitude between seeds. What does not move is where the transition happens — full accuracy up to $\theta = 1$, several digits gone by $3/2$, none left by $2$. The default sits a factor of two below the first of those, which is the room the residual bound of Applying the denominator without a matrix solve quantifies.

Lowering $\theta$ below $1/2$ is safe and only adds modified-squaring steps. Raising it is not available, and that asymmetry is deliberate. It is a restriction on the argument, not a redefinition of it: a NativePade(θ) and a ScaledSquaring(θ) built with the same $\theta$ scale identically and hand their kernels the same matrix, which is why every comparison on this page can hold $\theta$ fixed at the shared default.

Advantages. It never forms a matrix larger than $2n\times{}2n$ and uses only reductions, matrix products and a kernel-written identity, so it is the algorithm ScaledSquaring can be checked against on a backend that forbids scalar indexing. The two share their scaling rule and their recovery step, but nothing in the kernel: a truncated series against a rational function, and a division by $2^s$ against a Newton–Schulz inverse. Agreement between them is therefore evidence about the two kernels, though not about the framework they have in common. Its accuracy is indistinguishable from ScaledSquaring's and AugmentedPade's in Float64 throughout Agreeing with the exponential.

Disadvantages. Its fixed rational evaluation does more small matrix products than ScaledSquaring — around 1.8× the $\mathfrak{A}$ call in What they cost, still below AugmentedPade — and it allocates the most of the three algorithms that evaluate $\mathfrak{A}$ directly. An allocation count is the figure least likely to stay a constant factor on a backend where allocating can cost a synchronisation rather than a malloc, and the ordering of the three is itself unexplained; both are #77.

The Float64 indifference does not carry to Float32, where it has the worst check of the three while its forward error there is no outlier at all. What degrades is the orthogonality of the retracted point rather than the agreement with the exponential; both are in Float32 below.

4. AugmentedPade

AugmentedPade is the one algorithm here that chooses neither a kernel nor a recovery: it delegates both. For the $4n\times{}4n$ augmented matrix,

\[\exp\begin{pmatrix} X & \mathbb{I} \\ \mathbb{O} & \mathbb{O} \end{pmatrix} = \begin{pmatrix} \exp(X) & \mathfrak{A}(X) \\ \mathbb{O} & \mathbb{I} \end{pmatrix},\]

which is the standard device for recovering $\mathfrak{A}$ from a routine that computes only $\exp$ [26, 28]. One call to Base.exp therefore returns $\mathfrak{A}(X)$ in its upper-right block, and the numerics are Julia's own dense matrix exponential — a Padé approximant with its own scaling and squaring [18, 19]. This is the package algorithm that corresponds directly to the conventional Padé description of scaling and squaring, and the one place a reader can compare the two sections above against a textbook implementation.

Advantages. It introduces no package-specific approximation at all. Everything delicate is done by the most heavily exercised matrix-exponential implementation available, which is why it is the reference every accuracy table on this page is measured against. Its accuracy matches ScaledSquaring's throughout Agreeing with the exponential, and it is, unexpectedly, the lightest allocator of the three algorithms that evaluate $\mathfrak{A}$ — see What they cost — because Base.exp works in a few reused buffers where both native algorithms produce a fresh $2n\times{}2n$ temporary per operation.

Disadvantages. It exponentiates a matrix four times the size needed and discards three quarters of the result, which makes the $\mathfrak{A}$ call the most expensive of the three, though much less so once the $N\times{}N$ assembly around it is counted. And Base.exp on a dense matrix needs LAPACK, so it does not run on a GPU backend.

5. ProjectedSkew

ProjectedSkew does not go through $\mathfrak{A}$ at all. It nonetheless exploits the same rank-$2n$ structure as the other four, and nothing $N\times{}N$ is exponentiated here either; it simply reaches that structure by a change of basis rather than by The reduction. It exponentiates the lift in a basis of the lift's own range, where the lift is a small skew-symmetric matrix.

$\bar{B}$ is skew-symmetric of rank at most $2n$, so its range and its row space coincide and both sit inside the range of $B'$. A thin QR of $B'$ gives an $N\times{}2n$ orthonormal $Q$ with $\bar{B} = QMQ^T$, where $M = Q^T\bar{B}Q$ is skew-symmetric and $2n\times{}2n$, and

\[\exp(\bar{B}) = \mathbb{I} + Q\left(\exp(M) - \mathbb{I}\right)Q^T.\]

$\exp(M)$ is then formed from an eigendecomposition rather than from a series: $iM$ is Hermitian for real skew $M$, so $M = -iV\Lambda{}V^*$ with $V$ unitary and $\Lambda$ real, and

\[\exp(M) = \Re\left(V e^{-i\Lambda} V^*\right),\]

which is orthogonal by construction — a product of a unitary matrix, a diagonal of unit-modulus numbers, and a unitary matrix — rather than by the cancellation of a series. That is the property the whole algorithm exists for, and it is what the measurements below show: the orthogonality of the result does not depend on how accurately anything was summed.

Unlike the other four, ProjectedSkew bypasses $\mathfrak{A}$ and specialises geodesic directly. It therefore has no 𝔄(X, ProjectedSkew()) method — worth knowing if you call $\mathfrak{A}$ yourself.

Advantages. It is the only algorithm whose check does not degrade with the size of the lift: its column in Staying on the manifold is level where the other four climb by two orders of magnitude. The gap is widest in Float32, where the others are at the mercy of the format and this one is not. That is the case for choosing it: a long Float32 run, where the departure from the manifold accumulates over thousands of steps and staying on the manifold matters more than agreeing with the exponential to the last bit.

Disadvantages. It usually has the largest forward error against exp(Matrix(B)), in both formats and at all but the top of the sweep — see Agreeing with the exponential. It needs a qr and an eigen instead of matrix products, which costs more than ScaledSquaring per retraction in What they cost and rules out a GPU backend.

The trade is a real exchange rather than a ranking: orthogonality is structural here and accuracy against the exponential is not, and the two tables below disagree with each other for exactly that reason.

Using them

None of these types is exported, so import the ones you use:

using GeometricOptimizers
using GeometricOptimizers: Geodesic, Cayley, ScaledSquaring, NativePade, AugmentedPade, ProjectedSkew, check
import Random
Random.seed!(123)

Y = rand(StiefelManifold, 5, 3)
B = GeometricOptimizers.global_rep(GlobalSection(Y), rand(5, 3))

check(Geodesic()(B)) < 1e-14, check(Cayley()(B)) < 1e-14

# output

(true, true)

A retraction is passed to Optimizer as a keyword argument, and the algorithm travels inside it:

optimizer = Optimizer(ps, L; algorithm  = Adam(),
                             retraction = Geodesic(ProjectedSkew()))

Geodesic() is Geodesic(ScaledSquaring()), and ScaledSquaring(θ) accepts a scaling threshold; the default is 0.5:

Geodesic().algorithm == ScaledSquaring(0.5)

# output

true

The algorithm reaches the lower-level entry points too. geodesic takes it as an optional last argument, both for a tangent vector at a point and for a horizontal lift:

using GeometricOptimizers: geodesic, rgrad

Δ = rgrad(Y, rand(5, 3))

check(geodesic(Y, 300 * Δ, ProjectedSkew())) < 1e-13

# output

true

and GeometricOptimizers.𝔄 can be called on a bare matrix, which is the level at which four of the five algorithms are implemented — ProjectedSkew being the exception, as above:

using GeometricOptimizers: 𝔄
import Random
Random.seed!(123)

X = randn(6, 6)

isapprox(𝔄(X, ScaledSquaring()), 𝔄(X, NativePade()); rtol = 1e-12) &&
    isapprox(𝔄(X, NativePade()), 𝔄(X, AugmentedPade()); rtol = 1e-12)

# output

true

If what you want is the exponential itself rather than $\mathfrak{A}$, it is one line — $\exp(B'(B'')^T) = \mathbb{I} + B'\mathfrak{A}(B', B'')(B'')^T$, at a cost still set by $n$ rather than by $N$, since the only matrix function evaluated is $\mathfrak{A}$ on the $2n\times{}2n$ product. It is what geodesic computes before wrapping the result in a Manifold, with ScaledSquaring as geodesic's default:

using GeometricOptimizers: lift_factors
using LinearAlgebra: I
import Random
Random.seed!(1234)

B = 60 * rand(StiefelLieAlgHorMatrix, 20, 3)      # ‖B̄‖ ≈ 393
B̂, B̄ = lift_factors(B)

isapprox(I + B̂ * 𝔄(B̂, B̄, ScaledSquaring()) * B̄', exp(Matrix(B)); rtol = 1e-10)

# output

true

The difference between scaled and unscaled Taylor evaluation becomes visible for a large lift:

using GeometricOptimizers
using GeometricOptimizers: Geodesic, TaylorSeries, check
import Random
Random.seed!(1234)

B = 60 * rand(StiefelLieAlgHorMatrix, 20, 3)      # ‖B̄‖ ≈ 393

check(Geodesic()(B)) < 1e-12, check(Geodesic(TaylorSeries())(B)) < 1e-12

# output

(true, false)

6. Numerical comparison

Everything in this section other than the timings is recomputed when this page is built, so the figures are those of the version of the package the documentation was built from rather than a quote that can go stale. scripts/retraction_accuracy.jl produces the same tables — including the timings — from the command line.

The figures quoted in the algorithm sections above are read off one build of these tables. They are there to make each section's argument concrete, and the tables here are what to trust: round-off-level numbers depend on the BLAS and the machine, so the last digit of a quotation will not always match the table below it. What does not move is the shape of each column, and that is what the sections claim.

Staying on the manifold

check of the retracted point, $\|Y^TY - \mathbb{I}\|$, on a random StiefelLieAlgHorMatrix(20, 3) scaled up. This is the quantity a retraction is supposed to keep at round-off, and it is the one the test suite asserts on. Cayley is in the last column as the reference, since it evaluates no matrix function at all.

table(["‖B̄‖", "`ScaledSquaring`", "`NativePade`", "`AugmentedPade`", "`ProjectedSkew`", "`TaylorSeries`", "`Cayley`"],
      [[fixed(norm(Matrix(B))),
        sci(check(geodesic(B, ScaledSquaring()))),
        sci(check(geodesic(B, NativePade()))),
        sci(check(geodesic(B, AugmentedPade()))),
        sci(check(geodesic(B, ProjectedSkew()))),
        sci(check(geodesic(B, TaylorSeries()))),
        sci(check(cayley(B)))] for B in lifts])

TaylorSeries fails by many orders of magnitude and is the column the !!! danger above quotes. Every other column stays at round-off across the whole Float64 sweep, and only ProjectedSkew's is level. The other four — Cayley included, and it is the one that drifts furthest — grow by two to three orders of magnitude from one end to the other, because their orthogonality is an arithmetic outcome where ProjectedSkew's is structural. Round-off at $\|\bar{B}\| \approx 770$ is still round-off, so this separates the algorithms without condemning any of the four usable ones.

Agreeing with the exponential

Relative distance to exp(Matrix(B)), i.e. to the exponential of the full $N\times{}N$ lift. This is a different question from the one above — a retraction that re-orthonormalised its result would have a perfect check and be wrong here — and the test suite asserts both.

table(["‖B̄‖", "`ScaledSquaring`", "`NativePade`", "`AugmentedPade`", "`ProjectedSkew`"],
      [begin
           reference = exp(Matrix(B))
           err(algorithm) = norm(Matrix(geodesic(B, algorithm)) - reference) / norm(reference)
           [fixed(norm(Matrix(B))), sci(err(ScaledSquaring())), sci(err(NativePade())),
            sci(err(AugmentedPade())), sci(err(ProjectedSkew()))]
       end for B in lifts])

All four grow slowly with the norm of the lift, and the ordering is roughly the reverse of the previous table: ProjectedSkew is the furthest from the exponential at all but the largest of these norms. ScaledSquaring, NativePade and AugmentedPade sit within a factor of 1.6 of each other on every row, so the trade is between that group and ProjectedSkew: one is orthogonal by construction, the others agree with exp more closely. The four converge again at the top of the sweep, where the reference exp(Matrix(B)) is itself no more accurate than what is being measured against it. The table reports an experiment, not an error bound; the ordering can depend on the matrix and the floating-point format, which is what the next section is for.

Float32

The same orthogonality residual in Float32, the format the MNIST example in Optimization on Homogeneous Spaces actually runs in. Nothing can do better than about $10^{-6}$ here, but the four do not degrade alike.

table(["‖B̄‖", "`ScaledSquaring`", "`NativePade`", "`AugmentedPade`", "`ProjectedSkew`"],
      [[fixed(norm(Matrix(B))),
        sci(check(geodesic(B, ScaledSquaring()))),
        sci(check(geodesic(B, NativePade()))),
        sci(check(geodesic(B, AugmentedPade()))),
        sci(check(geodesic(B, ProjectedSkew())))] for B in sweep(Float32)])

ProjectedSkew is flat here too, and by the top of the sweep it is one to two orders of magnitude below the others. Which of the other three is worst depends on the lift over most of the range, but not at the top: there NativePade is the worst of them by a factor of about $2.5$, which is the one place the Float64 indifference of the previous table does not carry over.

And the forward error in the same format, which gives a different ranking and is worth having next to it. The reference is exp of the lift promoted to Float64, with the difference taken there as well: exp of a Float32 matrix is itself only Float32-accurate, so comparing against it would measure the reference as much as the algorithm.

table(["‖B̄‖", "`ScaledSquaring`", "`NativePade`", "`AugmentedPade`", "`ProjectedSkew`"],
      [begin
           reference = exp(Matrix{Float64}(Matrix(B)))
           err(algorithm) = norm(Matrix{Float64}(Matrix(geodesic(B, algorithm))) - reference) /
                            norm(reference)
           [fixed(norm(Matrix(B))), sci(err(ScaledSquaring())), sci(err(NativePade())),
            sci(err(AugmentedPade())), sci(err(ProjectedSkew()))]
       end for B in sweep(Float32)])

Here ProjectedSkew is the worst of the four at every norm — by 1.4× to 3.9× over most of the sweep and by 23× at the smallest lift, where it is the only one not at Float32 round-off — and the three $\mathfrak{A}$ algorithms are within 1.7× of each other throughout, NativePade included, its check outlier above notwithstanding. Taken together the two tables say what the trade actually is in Float32: ProjectedSkew buys orthogonality at the price of agreement, and it is the only one of the four for which that is a structural exchange rather than an accident of the arithmetic. Whether it matters in an optimization run depends on how the two errors accumulate in that application.

Sensitivity to the threshold θ

ScaledSquaring's only parameter, swept over a 32-fold range on one lift:

Random.seed!(99)
B = 30 * rand(StiefelLieAlgHorMatrix{Float64}, 20, 3)
reference = exp(Matrix(B))

table(["θ", "`check`", "error vs `exp`"],
      [begin
           Y = geodesic(B, ScaledSquaring(θ))
           [string(θ), sci(check(Y)), sci(norm(Matrix(Y) - reference) / norm(reference))]
       end for θ in (0.125, 0.25, 0.5, 1.0, 2.0, 4.0)])

Both columns move by less than a factor of six across the whole 32-fold range, and not monotonically. The default 0.5 sits inside that band and nothing in the measurement singles it out, which is the point: it needs no tuning because no value in the range does appreciably better. This is one lift and one seed, so it supports 0.5 as a reasonable default rather than establishing an optimal threshold, and it says nothing about NativePade, whose θ is a ceiling and not a preference.

What they cost

Unlike everything above, these are timings and therefore machine-dependent, so they are quoted rather than measured at build time. minimum of 50 repetitions, a single BLAS thread, on an Apple M-series laptop; julia --project=. scripts/retraction_accuracy.jl reproduces them on yours, and every figure in this section comes from one run of it. They are illustrative, and the only sound way to use them is to remeasure on the target hardware.

The $\mathfrak{A}$ call on its own, at $N = 200$, $n = 10$, is what separates the three algorithms that evaluate it:

The allocation row is the one figure here that is not machine-dependent — @allocated is exact — and it does not rank the three the way runtime does. AugmentedPade, which builds a $4n\times{}4n$ matrix and discards three quarters of the result, allocates the least, because Base.exp works in a few reused buffers where both native algorithms produce a fresh $2n\times{}2n$ temporary per operation. On a CPU that is a detail. On a backend where an allocation costs a synchronisation it may not be, which is worth knowing about the algorithm whose whole purpose is to be portable. Neither native algorithm reuses buffers, and both could; that is #77.

One whole retraction, which adds the $N\times{}N$ assembly they all share, in milliseconds:

and the same in KiB allocated:

ScaledSquaring and NativePade are level on runtime to within the run-to-run noise at every size but one; the isolated call above is where the extra rational work shows. The exception is $n = 50$, the one column where the $2n\times{}2n$ argument is large enough for the difference to survive the assembly: 3.33 against 3.03, with AugmentedPade's 5.98 as it pays for the $4n\times{}4n$ embedding. Allocations separate them everywhere and by more, up to 1.5× at $n = 50$, since that is the metric the shared assembly dilutes least at small $N$. ProjectedSkew stays close on both — a QR and an eigendecomposition of a $2n\times{}2n$ matrix are not expensive things, and it is the lightest of the five at $n = 50$. Cayley is level with the exponential up to $N \approx 50$ and loses by a factor of 15 by $N = 1000$, which is the $O(N^3)$ against $O(N^2n)$ of the Retractions page.

Choosing one

Cayley remains the package default. It requires only a small linear solve and is robust for large steps. Geodesic computes the exponential map and has the one-parameter subgroup property. In the measurements above it also becomes cheaper once the ambient dimension is large relative to the manifold dimension. Choose between them according to which map the algorithm needs, then benchmark representative problem sizes if cost matters.

For the algorithm: ScaledSquaring, i.e. the default, unless one of the alternatives has the property you specifically need.

ScaledSquaring, NativePade and TaylorSeries avoid dense LAPACK and scalar indexing in package code; AugmentedPade does not, and ProjectedSkew needs the backend's own qr and eigen. Whether any of the first three runs on a particular accelerator depends on that backend's matrix-multiplication and reduction support. ScaledSquaring is the default among them because it does less work, and NativePade is the independent cross-check that was previously missing there.

Adding one

A new algorithm is a subtype of AbstractExponentialAlgorithm that supplies

GeometricOptimizers.𝔄(X::AbstractMatrix, ::NewAlgorithm)

Everything above it — geodesic for lifts and for tangent vectors, retraction, update_section!, the optimizer — then follows, for both manifolds, with no further methods. An algorithm that does not evaluate $\mathfrak{A}$ at all supplies

GeometricOptimizers.geodesic(B::AbstractLieAlgHorMatrix, ::NewAlgorithm)

instead, which is what ProjectedSkew does. A new retraction is a subtype of AbstractRetraction supplying retraction(::NewRetraction, x); the callable form R(x) comes for free.

Whichever it is, test/retractions/exponential_accuracy.jl is where it earns its place: the sweep there asserts that a retracted point stays on the manifold at every lift norm and that it still agrees with exp(Matrix(B)), which together rule out both of the ways an exponential can be wrong here.

Reference