Gradients of Functions of Large Matrices
Nicholas Krämer, Pablo Moreno-Muñoz, Hrittik Roy, Søren Hauberg
Technical University of Denmark
Kongens Lyngby, Denmark{pekra, pabmo, hroy, sohau}@dtu.dk

Abstract
Tuning scientific and probabilistic machine learning models – for example, partial differential equations, Gaussian processes, or Bayesian neural networks –often relies on evaluating functions of matrices whose size grows with the data
set or the number of parameters. While the state-of-the-art for evaluating thesequantities is almost always based on Lanczos and Arnoldi iterations, the present
work is the first to explain how to differentiate these workhorses of numerical linear algebra efficiently. To get there, we derive previously unknown adjoint systems for Lanczos and Arnoldi iterations, implement them in JAX, andshow that the resulting code can compete with Diffrax when it comes to differentiating PDEs, GPyTorch for selecting Gaussian process models and beats
standard factorisation methods for calibrating Bayesian neural networks. All
this is achieved without any problem-specific code optimisation. Find the code
at https://github.com/pnkraemer/experiments-lanczos-adjoints andinstall the library with pip install matfree.

1

Introduction

Automatic differentiation has dramatically altered the development of machine learning models by
allowing us to forego laborious, application-dependent gradient derivations. The essence of thisautomation is to evaluate Jacobian-vector and vector-Jacobian products without ever instantiating
the full Jacobian matrix, whose column count would match the number of parameters of the neuralnetwork. Nowadays, everyone can build algorithms around matrices of unprecedented sizes by
exploiting this matrix-free implementation. However, differentiable linear algebra for Jacobian-vector
products and similar operations has remained largely unexplored to this day. We introduce a new
matrix-free method for automatically differentiating functions of matrices. Our algorithm yields theexact gradients of the forward pass, all gradients are obtained with the same code, and said code runs
in linear time- and memory-complexity.
For a parametrised matrix A = A(θ) ∈ RN ×N and an analytic function f : R → R, we call f (A) afunction of the matrix (different properties of A imply different definitions of f (A); one of them isapplying f to each eigenvalue of A if A is diagonalisable; see [1]). However, we assume that A is theJacobian of a large neural network or a matrix of similar size and never materialise f (A). Instead, weonly care about the values and gradients of the matrix-function-vector product(θ, v) 7→ f [A(θ)]v(1)

assuming that A is only accessed via differentiable matrix-vector products. Table 1 lists examples.
Evaluating Equation 1 is crucial for building large machine learning models, e.g., Bayesian neuralnets: A common hyperparameter-calibration loss of a (Laplace-approximated) Bayesian neuralnetwork involves the log-determinant of the generalised Gauss–Newton matrix [13]XA(α) :=[Dθ g](xi )⊤ [Dg2 ρ](yi , g(xi ))[Dθ g](xi ) + α2 I,(2)(xi ,yi )∈data

Preprint. Under review.

Table 1: Some applications for functions of matrices. Log-determinants apply by combining
log det(A) = trace (log(A)) with stochastic trace estimation, which is why most vectors in thistable are Rademacher samples. “PDE” / “ODE” = “Partial/Ordinary differential equation”.
Function f

Application
PDEs & flows [2–5]
Gaussian process [6–8]
Invert. ResNets [9, 10]
Gaussian sampler [11]Neural ODE [12]λelog(λ)
log(1 + λ)√λλ2

Matrix A

Vector v

Parameter θ

PDE discret.
Kernel matrix
Jacobian matrix
Covariance matrix
Jacobian matrix

PDE initial valuev ∼ Rademacherv ∼ Rademacherv ∼ N (0, I)v ∼ Rademacher

PDE
KernelNetwork
Covariance
Network

where Dθ g is the parameter-Jacobian of the neural network g, Dg2 ρ is the Hessian of the loss function
ρ with respect to g(xi ), and α is a to-be-tuned parameter. The matrix A(α) in Equation 2 has as manyrows and columns as the network has parameters, which makes traditional, cubic-complexity linear
algebra routines for log-determinant estimation entirely unfeasible. To compute this log-determinant,one chooses between either (i) simplifying the problem by pretending that the Hessian matrix ismore structured than it actually is, e.g., diagonal [14]; or (ii) approximating log det(A) by combining
stochastic trace estimation [15]
L


1X ⊤
vℓ Avℓ ,trace (A) = E v ⊤ Av ≈
L

for



E vv ⊤ = I,

(3)

ℓ=1

with a Lanczos iteration A(θ) ≈ QHQ⊤ [16], to reduce the log-determinant to [17, 18]
log det(A) = trace (log A) ≈

L

L

ℓ=1

ℓ=1

1X ⊤
1X ⊤
vℓ log(A)vℓ ≈vℓ Q log(H)Q⊤ vℓ .
L
L

(4)

The matrix H in A ≈ QHQ⊤ has as many rows/columns as we are willing to evaluate matrix-vector
products with A; thus, it is small enough to evaluate the matrix-logarithm log(H) in cubic complexity.
Contributions This article explains how to differentiate
not just log-determinants but any Lanczos and Arnoldi
iteration so we can build loss functions for large models
with such matrix-free algorithms (thereby completingthe pipeline in Figure 1). This kind of functionality has
been sorely missing from the toolbox of differentiable
programming until now, even though the demand forfunctions of matrices is high in all of probabilistic andscientific machine learning [e.g. 2–12, 19–31].

Parameter θLarge matrix-vector productArnoldi, Lanczos
Small factorisation“Cheap” linalg

2

Related work

Optimise θ

Function of large matrix∇
Contribution∇(known)∇Here, we focus on applications in machine learning andillustrate how prior work avoids differentiating matrixLoss(known)free decomposition methods like the Lanczos and Arnoldi
iterations. Golub and Meurant [32] discuss applications Figure 1: Values (down) and gradientsoutside machine learning.
(up) of functions of large matrices.
Generative models, e.g., normalising flows [27, 33], rely on the change-of-variables formula, which
involves the log-determinant of the Jacobian matrix of a neural network. Behrmann et al. [9] and Chenet al. [10] combine stochastic trace estimation with a Taylor-series expansion for the matrix logarithm.
Ramesh and LeCun [34] use Chebyshev expansions instead of Taylor expansions. That said, Ubaru
et al. [17] demonstrate how both methods converge more slowly than the Lanczos iteration whencombined with stochastic trace estimation.
Gaussian process model selection requires values and gradients of log-probability density functionsof Gaussian distributions (which involve log-determinants), where the covariance matrix A(θ) has as
many rows and columns as there are data points [35]. Recent work [6, 7, 11, 19, 23, 36] all uses some2

combination of stochastic trace estimation with the Lanczos iteration, and unanimously identifies
gradients of log-determinants as (“d” shall be an infinitesimal perturbation; see Section 4)µ := log det(K(θ)), dµ = trace K(θ)−1 dK(θ) .
(5)Another round of stochastic trace estimation then estimates dµ [6, 23, 36]. In contrast, our contribution
is more fundamental: not only do we derive the exact gradients of the forward pass, but our formulation
also applies to, say, matrix exponentials, whereas Equation 5 only works for log-determinants.
Section 5 shows how our black-box gradients match state-of-the-art code for Equation 5 [6].
Laplace approximations and neural tangent kernels face the same problem of computing derivativesof log-determinants but with the generalised Gauss–Newton (GGN) matrix from Equation 2. In
contrast to the Gaussian process literature, prior work on Laplace approximations prefers structured
approximations of the GGN by considering subsets of network weights [37–39], or algebraic approximations of the GGN via diagonal, KFAC, or low-rank factors [31, 40–44]. All such approximations
imply simple expressions for log-determinants, which are straightforward to differentiate automatically. Unfortunately, these approximations discard valuable information about the correlation between
weights, so a linear-algebra-based approach leads to superior likelihood calibration (Section 7).
Linear differential equations, for instance ẏ(t) = Ay(t), y(0) = y0 are solved by matrix exponentials,
y(t) = exp(At)y0 . By this relation, matrix exponentials have frequent applications not just for thesimulation of differential equations [e.g. 2, 45], but also for the construction of exponential integrators[3, 26, 29], state-space models [5, 46], and in generative modelling [4, 26, 28, 47]. There are many
ways of computing matrix exponentials [48, 49], but only Al-Mohy and Higham [50] consider the
problem of differentiating it and only in forward mode. In contrast, differential equations have a richhistory of adjoint methods [e.g. 51, 52] with high-performance open-source libraries [53–56]. Still,the (now differentiable) Arnoldi iteration can compete with state-of-the-art solvers in JAX (Section 6).

3

Problem statement

Recall A = A(θ) ∈ RN ×N from Section 1. The focus of this paper is on matrices that are too large
to store in memory, like Jacobians of neural networks or discretised partial differential equations:
Assumption 3.1. A(θ) is only accessed via differentiable matrix-vector products (θ, v) 7→ A(θ)v.
The de-facto standard for linear algebra under Assumption 3.1 are matrix-free algorithms [e.g.
57, Chapters 10 & 11], like the conjugate gradient method for solving large sparse linear systems [58]. But there is more to matrix-free linear algebra than conjugate gradient solvers:
Matrix-free implementations of matrix decompositionseK
H
usually revolve around variations of the Arnoldi iteraNtion [59], which takes an initial vector v ∈ R and a
Q = Q
A(θ)
+ r
prescribed number of iterations K ∈ N and produces a
N ×K
K×Kcolumn-orthogonal Q ∈ R, structured H ∈ R,
residual vector r ∈ RN , and length c ∈ R such that
AQ = QH + r(eK )⊤ ,and Qe1 = cv

(6)

Figure 2: Lanczos/Arnoldi iteration.

Khold (Figure 2; e1 , eK ∈ R are the first and last unit vectors). If A is symmetric, H is tridiagonal,
and the Arnoldi iteration becomes the Lanczos iteration [16]. Both iterations are popular for implementing matrix-function-vector products in a matrix-free fashion [1, 57], because the decomposition
in Equation 6 implies A ≈ QHQ⊤ , thus(θ, v) 7→ f (A(θ))v ≈ Qf (H)Q⊤ v = c−1 Qf (H)e1 .

(7)

The last step, Q⊤ v = c−1 e1 , is due to the orthogonality of Q. Since the number of matrix-vector
products K rarely exceeds a few hundreds or thousands, the following Assumption 3.2 is mild:
Assumption 3.2. The map H 7→ f (H)e1 is differentiable, and Q fits into memory.
In summary, we evaluate functions of large matrices by firstly decomposing a large matrix
into a product of small matrices (with Lanczos or Arnoldi) and, secondly, using conventional
linear algebra to evaluate functions of small matrices. Functions of small matrices can already be differentiated efficiently [60–62]. This work contributes gradients of the Lanczos andArnoldi iteration under Assumptions 3.1 and 3.2, and thereby makes matrix-free implementations of matrix decompositions and functions of large matrices (reverse-mode) differentiable.
3

0.6
0.5

Wall time (sec)Automatic differentiation, i.e. “backpropagating through”
the matrix decomposition, is far too inefficient to be a viable option (Figure 3; setup in Appendix A). Our approach
via implicit differentiation or the adjoint method, respectively, leads to gradients that inherit the linear runtime andmemory-complexity of the forward pass.

SuiteSparse: 'bcsstk18'

0.4
0.3

N = 11948 rows

0.2

Backprop

Adjoint

Limitations and future work The landscape of Lanczos0.1
and Arnoldi-style matrix decompositions is vast, and some
Forward passadjacent problems cannot be solved by this single article:
0.0(i) Forward-mode derivatives would require a derivation
0
50 100 150 200 250
separate from what comes next. Yet, since functions of
Krylov-space depth
matrices map many to few parameters (matrices to vectors), reverse-mode is superior to forward-mode anyway Figure 3: Backpropagation vs our ad[66, p. 153]. (ii) We only consider real-valued matrices joint method on a sparse matrix [63–65].
(for their relevance to machine learning), even though the decompositions generalise to complex
arithmetic with applications in physics [67]. (iii) We assume Q fits into memory, which relates tocombining Arnoldi/Lanczos with full reorthogonalisation [57, 68–71]. Relaxing this assumption
requires gradients of partial reorthogonalisation (among other things), which we leave to future work.

4

The method: Adjoints of the Lanczos and Arnoldi iterations

Numerical algorithms are rarely differentiated automatically, and usually, some form of what is known
as “implicit differentiation” [66, 72] applies. The same is true for the Lanczos and Arnoldi iterations.
However, and perhaps surprisingly, we differentiate the iterations like a dynamical system using the“adjoint method” [51, 73, 74], a variation of implicit differentiation that uses Lagrange multipliers[66], and not like a linear algebra routine [60, 61, 75]. To clarify this distinction, we briefly review
implicit differentiation before the core contributions of this work in Sections 4.1 and 4.2.
Notation Let dx be an infinitesimal perturbation of some x. D is the Jacobian operator, and ⟨·, ·⟩ the
Euclidean inner product between two equally-sized inputs. For a loss ρ ∈ R that depends on some x,
the linearisation dρ = Dx ρ dx and the gradient identity dρ = ⟨∇x ρ, dx⟩ will be important [76, 77].
Implicit differentiation Let a : θ 7→ x be a numerical algorithm that computes some x from someθ. Assume that the input and output of a(·) satisfy the constraint c(θ, a(θ)) = 0. For instance, ifa(·) solves Ax = b, the constraint is c(A, b; x) = Ax − b with θ := {A, b}. We can use c(·) in
combination with the chain rule to find the derivatives of a(·),(c = 0 implies dc = 0)
0 = dc(θ, x) = Dx c(θ, x)dx + Dθ c(θ, x)dθ.
(8)
In other words, we “linearise” the constraint c(θ, x) = 0. The adjoint method [78] proceeds by“transposing” this linearisation as follows. Let ρ be a loss that depends on y with gradient ∇y ρ. Then,for all Lagrange multipliers λ with the same shape as the outputs of c(·), we knowdρ = ⟨∇x ρ, dx⟩ + ⟨λ, dc⟩ = ⟨∇x ρ + (Dx c)⊤ λ, dx⟩ + ⟨(Dθ c)⊤ λ, dθ⟩.
By matching Equation 9 to dρ = ⟨∇θ ρ, dθ⟩, we conclude that if λ solves the adjoint system∇x ρ + (Dx c)⊤ λ = 0,
⊤

(9)
(10)

then ∇θ ρ := (Dθ c) λ must be the gradient of ρ with respect to input θ. This is the adjoint method[66, Section 10.4]. In automatic differentiation frameworks like JAX [79], this gradient implementsa vector-Jacobian product with the Jacobian of a(·) – implicitly via the Lagrange multiplier λ,without differentiating “through” a(·) explicitly. In comparison to approaches that explicitly target
vector-Jacobian products with implicit differentiation [like 66, Proposition 10.1], the adjoint method
shines when applied to highly structured, non-vector-valued constraints, such as dynamical systems
or the Lanczos and Arnoldi iterations. The reason is that the adjoint method does not change ifc(·) becomes matrix- or function-space-valued, as long as we can define inner products and adjoint
operators, whereas other approaches (like what Blondel et al. [72] use for numerical optimisers)
would become increasingly laborious in these cases. In summary, to reverse-mode differentiate anumerical algorithm with the adjoint method, we need four steps: (i) find a constraint, (ii) linearise it,(iii) introduce Lagrange multipliers, and (iv) solve the resulting adjoint system. Carrying out thosefour steps for the Lanczos and Arnoldi iterations is the main contribution of the paper: Section 4.1
states both adjoint systems and Section 4.2 covers a matrix-free implementation.
44.1

Adjoint system of the Arnoldi and Lanczos iterations

Let ej be the jth unit vector. Denote by “◦” the element-wise matrix product, and define the matrices
I≤ := [δi≤j ]Ki,j=1 ,

I< := [δi<j ]Ki,j=1 ,

I≪ := [δi+1<j ]Ki,j=1 ,

(11)

so that for example, I≤ ◦ A extracts the lower triangular matrix of A (including the diagonal),and I≪ ◦ A = 0 enforces Hessenberg form [57]. The following two theorems do not require
Assumptions 3.1 and 3.2, which are only relevant for analysing the computational complexities.
Theorem 4.1 (Adjoint system of the Arnoldi iteration). Let K ∈ N, v ∈ R, and A ∈ RN ×N , and aloss ρ(·) ∈ R be given. If Q ∈ RN ×K , H ∈ RK×K , r ∈ RN , and c ∈ R solve the forward constraint
AQ = QH + r(eK )⊤ ,N

and if λ ∈ R , Λ ∈ R

N ×K

Qe1 = vc,
KI≤ ◦ [Q⊤ Q] = I,,γ ∈R ,Γ∈R

⊤

K×K

⊤

, and Σ ∈ R

I≪ ◦ H = 0,
K×K

Q⊤ r = 0,(12)satisfy the adjoint system⊤

0 = ∇Q ρ + A Λ − ΛH + λ(e1 ) + Q(I≤ ◦ Γ) + Q(I≤ ◦ Γ)⊤ + rγ ⊤

(13a)

⊤

0 = ∇H ρ − Q Λ + I≪ ◦ Σ
0 = ∇r ρ − ΛeK + Qγ

(13b)
(13c)

0 = ∇c ρ − v ⊤ λ,(13d)then the gradients of ρ with respect to A and v are∇A ρ := ΛQ⊤ ,∇v ρ := λc.

(14)Sketch of the proof. To derive the statement, start with Equation 12 as c(·). Apply the chain- andproduct rules liberally to get dc(·). Introduce Lagrange multipliers λ, Λ, γ, Γ, and Σ like in theprevious section, by adding Lagrange-multiplied constraints todρ = ⟨∇Q ρ, dQ⟩ + ⟨∇H ρ, dH⟩ + ⟨∇r ρ, dr⟩ + ⟨∇c ρ, dc⟩,

(15)

and rearrange the terms to see that Equation 13 implies Equation 14. Details are in Appendix B.
Theorem 4.2 (Adjoint system of the Lanczos iteration). Let a symmetric A ∈ RN ×N , as well
as v ∈ RN , K ∈ N, and a loss ρ be known. In the following equations, set b0 := 1 ∈ R,
x0 := 0 ∈ Rn , λK+1 := 0, µ0 := 0, ν0 := 0 to simplify the expressions. If x1 , ..., xK+1 ∈ RN , anda1 , ..., ak , b1 , ..., bk ∈ RK , satisfy the forward constraint
x1 − v/(v ⊤ v) = 0,−bk−1 xk−1 + (A − ak I)xk − bk xk+1 = 0,

k = 1, ..., K

(16a)
(16b)

x⊤
k+1 xk+1 − 1 = 0,

k = 1, ..., K,(16c)

k = 2, ..., K + 1

(16d)

x⊤
k−1 xk = 0,N

and if λ0 , ..., λK ∈ R , µ1 , ..., µK , ν1 , ..., νK ∈ R satisfy the adjoint system
0 = −λK bK + (∇xK+1 ρ + µK xK+1 + νK xK ),

(17a)

⊤

0 = −bk λk+1 + (A − ak I)λk − bk−1 λk−1 + (∇xk ρ + µk−1 xk + νk xk+1 + νk−1 xk−1 ),(17b)
0 = ∇ak ρ − λ⊤k xk ,(17c)⊤
0 = ∇bk ρ − λ⊤k+1 xk − λk xk+1 .

(17d)where all expressions involving k hold for all k = K, ..., 1, then∇v ρ :=

λ⊤0 x1x1 − λ0 ,v⊤ v∇A ρ :=

K
X

λk x⊤
k

(18)

k=1

are the gradients of ρ with respect to v and A.
Sketch of the proof. This theorem is proven similarly to that of Theorem 4.1, but instead of a few
equations involving matrices, we have many equations involving scalars because for symmetricmatrices, H must be tridiagonal [57], and we expand AQ = QH + r(eK+1 )⊤ column-wise. The
coefficients ak and bk are the tridiagonal elements in H. We rename qk from Arnoldi to xk for
Lanczos to make it easier to distinguish the two different sets of constraints. Details: Appendix C.
5

4.2

Matrix-free implementation

Solving the adjoint systems To compute ∇A ρ and ∇v ρ, we need to solve the adjoint systems.
When comparing the forward constraints to the adjoint systems, similarities emerge: for instance, theadjoint system of the Arnoldi iteration follows the same A(⊤) H − ΛH (⊤) + rest = 0 structure as theforward constraint. This structure suggests deriving a recursion for the backward pass that mirrorsthat of the forward pass. Appendix E contains this derivation and contrasts the resulting algorithmwith that of the forward pass. The main observation is that the complexity of the adjoint passes for
Lanczos and Arnoldi mirrors that of the forward passes. Gradients can be implemented purely withmatrix-vector products, which explains the significant performance gains observed in Figure 3.
Theorems 4.1 and 4.2’s expressions for ∇A ρ are not directly applicable when we only have matrixvector products with A. Fortunately, parameter-gradients emerge from matrix-gradients:
Corollary 4.3 (Parameter gradients). Under Assumption 3.1 and the assumptions of Theorem 4.1,
and if A is parametrised by some θ, the gradients of ρ with respect to θ are∇θ ρ =

KX∇ θ 7→ (ek )⊤ Q⊤ A(θ)⊤ Λek

(19)

k=1

which can be assembled online during the backward pass. For the Lanczos iteration, we assume theconditions of Theorem 4.2 instead of Theorem 4.1, replace Qek and Λek with xk and λk , let the sum
run from k = 0 to k = K, and the rest of this statement remains true.
Sketch of the proof. The proof of this identity combines the expression(s) for ∇A ρ from Theorems 4.1
and 4.2 with dA = Dθ Adθ. The derivations are lengthy and therefore relegated to Appendix D.
Reorthogonalisation It is well known that the Lanczos Table 2: Accuracy loss when differentiand Arnoldi iterations suffer from a loss of orthogonality ating the Arnoldi iteration on a Hilbert
and that reorthogonalisation of the columns in Q is often matrix in double precision (Appendix F).
necessary [68–71]. Reorthogonalisation does not affectthe forward constraints, so the adjoint systems remain the
Loss of accuracy
same with and without reorthogonalisation. But adjoint
Adjoint w/o proj.
5.83 · 10−3systems also suffer from a loss of orthogonality: The
equivalent of orthogonality for the adjoint system is the
Adjoint w/ proj.
1.17 · 10−10
projection constraint in Equation 13b, which constrains
1.17 · 10−10
Backprop.
the Lagrange multipliers Λ to a hyperplane defined by Q
and other known quantities. The constraint can – and should (Table 2) – be used whenever the forward
pass requires reorthogonalisation.1 In the case studies below, we always use full reorthogonalisationon the forward and adjoint pass, also for the Arnoldi iteration [71, Table 7.1], even though this isslightly less common than for the Lanczos iteration.
Summary (before the case studies) The main takeaway from Sections 4.1 and 4.2 is that now,
we do not only have closed-form expressions for the gradients of Arnoldi and Lanczos iterations(Theorems 4.1 and 4.2), but that we can compute them in the same complexity as the forward
pass, in a numerically stable way, and evaluate parameter-gradients in linear time- and spacecomplexity (Corollary 4.3). While some of the derivations are somewhat technical, the overall
approach follows the general template for the adjoint method relatively closely. The resulting
algorithm beats backpropagation “through” the iterations by a margin in terms of speed (Figure 3)
and enjoys the same stability gains from reorthogonalisation as the forward pass (Table 2). Our
open-source implementation of reverse-mode differentiable Lanczos and Arnoldi iterations can beinstalled via pip install [redacted]. Next, we put this code to the test on three challenging
machine-learning problems centred around functions of matrices to see how it fares against state-ofthe-art differentiable implementations of exact Gaussian processes (Section 5), differential equation
solvers (Section 6), and Bayesian neural networks (Section 7).
1
The adjoint system of the Lanczos iteration does not admit this projection constraint, but we can implement
re-orthogonalised Lanczos via calling the Arnoldi code. This induces only minimal overhead because fully
reorthogonalised Lanczos code has roughly the same complexity as Arnoldi code.

6

Table 3: Our method yields the same root-mean-square errors (RMSEs) as GPyTorch. It reaches
lower training losses but is ≈ 20× slower per epoch due to different matrix-vector-product backends(see Appendix G). Three runs, significant improvements in bold. We use an 80/20 train/test split.
Dataset

Size

Dim.

MethodRMSE ↓

Final training loss ↓Runtime (s/epoch) ↓elevators

16,59918

Adjoints
GPyTorch

0.09 ± 0.002
0.09 ± 0.003-0.91 ± 0.025-0.63 ± 0.062

1.69 ± 0.000
0.10 ± 0.004

protein

45,7309Adjoints
GPyTorch

0.39 ± 0.005
0.39 ± 0.005

0.73 ± 0.300
0.73 ± 0.075

12.62 ± 0.172
0.73 ± 0.039

kin40k

40,0008Adjoints
GPyTorch

0.12 ± 0.004
0.10 ± 0.010-0.30 ± 0.078-0.26 ± 0.094

8.27 ± 0.004
0.26 ± 0.024

kegg dir

48,82720

Adjoints
GPyTorch

0.12 ± 0.002
0.12 ± 0.005-0.59 ± 0.295-0.41 ± 0.054

13.25 ± 0.005
0.62 ± 0.262

kegg undir 63,608

26

Adjoints
GPyTorch

0.12 ± 0.002
0.12 ± 0.003-0.69 ± 0.263-0.40 ± 0.039

24.20 ± 0.004
1.67 ± 0.532

5Case study: Exact Gaussian processes

Model selection for Gaussian processes has arguably been the strongest proponent of the Lanczos
iteration and similar matrix-free algorithms in recent years [6, 7, 11, 19, 20, 23, 80], and most ofthese efforts have been bundled up in the GPyTorch library [6]. For example, GPyTorch defaults
to choosing a Lanczos iteration over a Cholesky decomposition as soon as the dataset exceeds 800
data points.2 Calibrating hyperparameters of Gaussian process models involves optimising logmarginal-likelihoods of the regression targets, which requires computing x⊤ A−1 x and log det(A)for a covariance matrix A with as many rows and columns as there are data points. Recent works [6,7, 11, 19, 20, 23, 80] unanimously
suggest to differentiate log-determinants via µ := trace (log(A))
and dµ = trace A−1 dA (Equation 5). Since we seem to be the first to take a different path,benchmarking Gaussian processes in comparison to GPyTorch is a good first testbed for our gradients.
Setup: Like GPyTorch’s defaults We mimic recent suggestions for scalable Gaussian processmodels [7, 19]: we implement a pivoted Cholesky preconditioner [81] and combine it with conjugate gradient solvers for x⊤ K −1 x (which can be differentiated efficiently). We estimate thelog-determinant stochastically via log det(A) = trace (log(A)) = E[v ⊤ log(A)v], and computelog(A)v via the Lanczos iteration. While all of the above is common for “exact” Gaussian processes[6, 7, 19] (“exact” as opposed to variational approaches, which are not relevant for this comparison),there are three key differences between our code and GPyTorch’s: (i) GPyTorch is in Pytorch anduses KeOps [82] for efficient kernel-matrix-vector products. We use JAX and must build our own
low-memory matrix-vector products (Appendix G). (ii) GPyTorch runs all algorithms adaptively (wespecify tolerances and maximum iterations as much as possible). We use adaptive conjugate gradientsolvers and fixed ranks for everything else. (iii) GPytorch differentiates the log-determinant witha tailored approximation of Equation 5 [6]; we embed our gradients of the Lanczos iteration into
automatic differentiation. To keep the benchmark objective, we mimic the parametersuggestions
from GPyTorch’s default settings, and optimise hyperparameters of a Matérn 32 model on UCI
datasets with the Adam optimiser [83]. Appendix H lists parameters and discusses the datasets.
Analysis: Trains like GPyTorch; large scale only limited by matrix-vector-product backends In
this benchmark, we are looking for low reconstruction errors, fast runtimes and well-behaved loss
functions. Table 3 shows that this is the case for both implementations: the reconstruction errors are
essentially the same, and both methods converge well (we achieve lower training losses). This resultshows that by taking a numerically exact gradient of the Lanczos iteration, and leaving everything elseto automatic differentiation, matches the performance of state-of-the-art solvers. Larger datasets are
only limited by the efficiency of our matrix-vector products (in comparison to KeOps); Appendix G
discusses this in detail. Overall, this result strengthens the democratisation of exact Gaussian processes
because it reveils a simple yet effective alternative to GPyTorch’s domain-specific gradients.
2

Parameter max cholesky size: https://docs.gpytorch.ai/en/stable/settings.html.

7

Forward error

A2.

Arnoldi (ours)
Dopri5
Tsit5

100

10 2

10 4

10 4

No. Matvecs

102

Training loss

B.

Arnoldi (ours)
Dopri5
Tsit5

100

10 2

101

Gradient error

Arnoldi (ours), 3000 epochs
Dopri5, 3000 epochs
Tsit5, 3000 epochs

Loss

RMSE (relative)

A1.

10 1
10 2

101

No. Matvecs

0102

100 200 300 400 500 600 700 800

Time (sec)Figure 5: Arnoldi’s superior convergence on the forward pass (A1) is inherited by the gradients (A2)and ultimately leads to fast training (B). For training, Arnoldi uses ten matrix-vector products, andthe other two use 15 matrix-vector products (so they have equal error ≈ 10−4 in A1 and A2.)
Table 4: All three methods reconstruct the parameter well (std.-deviations exceed differences for testloss and RMSE), but Arnoldi and Dopri5 are faster than Tsit5. Dopri5 uses the BacksolveAdjoint,
and Tsit5 the RecursiveCheckpointAdjoint in Diffrax [54]. We contribute Arnoldi’s adjoints.

Loss on test set
Parameter RMSE
Runtime per epoch6

Arnoldi (adjoints; ours)

Dopri5 (diff. → disc.)

Tsit5 (disc. → diff.)

6.1e-03 ± 3.3e-04
2.9e-04± 4.4e-05
7.7e-02± 1.8e-05

6.3e-03±5.7e-04
2.6e-04 ±5.0e-05
7.2e-02 ±3.4e-05

5.9e-03 ±2.2e-04
2.7e-04± 5.2e-05
2.7e-01 ±1.1e-05

Case study: Physics-informed machine learning with PDEs

Much of the paper thus far discusses functions of matrices in the context of log-determinants. So,
in order to demonstrate performance for (i) a problem that is not a log-determinant and (ii) for anon-symmetric matrix which requires Arnoldi instead of Lanczos, we learn the coefficient field ω of
 2

∂2
∂
∂2
2
u(t;
x
,
x
)
=
ω(x
,
x
)
u(t;
x
,
x
)
+
u(t;
x
,
x
)
(20)
1
2
1
2
1
2
1
2
∂t2
∂x21
∂x22
subject to Neumann boundary conditions. We discretise this equation on a 128 × 128 grid in space
and transform the resulting 1282 -dimensional second-order ordinary differential equation into afirst-order differential equation, ẇ = Aw, w(0) = w0 , with solution operator w(t) = exp(At)w0 .
The system matrix A is sparse, asymmetric, and has 32, 768 rows and columns. We sample a true ωfrom a Gaussian process with a square exponential kernel and generate data by sampling 256 initialconditions and solving the equation numerically with high precision. Details are in Appendix I.
Setup: Arnoldi vs Diffrax’s Runge-Kutta methods for a 250k parameter MLP We learn
ω with a multi-layer perceptron (MLP) with approximately 250,000 parameters. We had similarreconstructions with fewer parameters but use 250,000 to display how gradients of the Arnoldi
iteration scale to many parameters. We compare an implementation of the solution operator (θ, w0 ) 7→
exp(A(θ))w0 with the Arnoldi iteration to Diffrax’s [54] implementation of “Dopri5” [84, 85]with a differentiate-then-discretise adjoint [86] as well
Solver: Arnoldi (ours)
Truth
0.014 1.0
0.014
as “Tsit5” [87] with a discretise-then-differentiate ad- 1.0
0.012
0.012
joint (recommended by [53, 54]). All methods receive
0.010
0.010
0.008
0.008
equal matrix-vector products per simulation.
0.5
0.5
0.006
0.006
Analysis: All methods train, but Arnoldi is more accurate for fixed matrix-vector-product budgets We
evaluate the approximation errors in computing thevalues and gradients of a mean-squared error loss for
all three solvers and then use the solvers to train the
MLP. We are looking for low approximation errors
for few matrix-vector products and for a good reconstruction of the truth. Figure 5 shows the results.
The Arnoldi iteration has the lowest forward-pass andgradient error, but Table 4 demonstrates how all ap8

0.0

0.0

1.0

0.5

Solver: Dopri5

1.0

0.5
0.0

0.0

0.5

1.0

0.004
0.002
0.000 0.0
0.014 1.0
0.012
0.010
0.008 0.5
0.006
0.004
0.002
0.000 0.0

0.0

0.0

0.5

Solver: Tsit5

0.5

1.0

1.0

Figure 4: All methods find the truth.

0.004
0.002
0.000
0.014
0.012
0.010
0.008
0.006
0.004
0.002
0.000

Table 5: Lanczos outperforms the diagonal approximation for calibrating a Bayesian version of an ImageNet pre-trained VAN. One training run; calibration estimated with 30 samples (sampling “Lanczos”
and “diagonal” with another Lanczos/diagonal approximation; see Appendix J). places365 [91] isthe out-of-distribution data; mean and std.-deviations of 3 runs. “MVs”: matrix-vector products.
Lanczos (50 MVs)

Lanczos (150 MVs)

Diagonal (150 MVs)

Runtime (sec)Marginal likelihood (log)

↓
↑

950
-192,757

4674-154,475

4314-876,444

Joint likelihood: train (log)
Joint likelihood: test (log)
Expected calibration error
AUROC (out-of-dist.)

↑
↑
↓
↑

-81.2 ± 15.4
-66.5 ± 11.2
0.5 ± 0.01
0.9 ± 0.03

-81.2 ± 15.4
-66.5 ± 11.2
0.5 ± 0.01
0.9 ± 0.03

-5,669.2 ± 124.1-5,260.9 ± 290.2
0.2 ± 0.003
0.5 ± 0.01

proaches lead to low errors on ω as well as on a test set (a held-back percentage of the training data);
see also Figure 4. The adjoints of the Arnoldi iteration match the efficiency of the differentiate-thendiscretise adjoint [86], and both outperform the discretise-then-differentiate adjoint by a margin.
This shows how linear-algebra solutions to matrix exponentials can compete with highly optimised
differential equation solvers. We anticipate ample opportunities of using the now-differentiable
Arnoldi iteration for physics-based machine learning.

7

Case study: Calibrating Bayesian neural networks

Next, we differentiate a function of a matrix on a problem that is native to machine learning: marginal
likelihood optimisation of a Bayesian neural network (a high-level introduction is in Appendix J).
Setup: Laplace-approximation of a VAN pre-trained on ImageNet We consider as gθ (x) a“Visual Attention Network” [88] with 4,105,800 parameters, pre-trained on ImageNet [89]. We
assume p(θ) = N (0, α−2 I), and Laplace-approximate the log-marginal likelihood of the data aslog p(y | x) ≈ log p(y, θ | x) −1
log det(A(α)) + const
2(21)

Neg. Log-marg.

where A(α) is the generalised Gauss–Newton matrix (GGN) from Section 1 (recall Equation 2).
We optimise α via Equation 21, implementing the log-determinant
×106
via stochastic trace estimation in combination with a Lanczos it1.0
eration (like in Section 5). Contemporary works [31, 37–43] rely
0.8
on sparse approximations of the GGN (such as diagonal or KFAC
approximations), so we compare our implementation to a diagonal
Diagonal
0.6
Lanczos
approximation of the GGN matrix, which yields closed-form logdeterminants. The exact diagonal of the 4-million-column GGN
0.4
matrix would require 4 million GGN-vector products with unit vectors, and like Deng et al. [90], we find this too expensive and resort
0.2
to stochastic diagonal approximation (similar to trace estimation;
0
100 200 300 400
all details are in Appendix J). We give both the stochastic diagoEpochs
nal approximation and our Lanczos-based estimator exactly 150
matrix-vector products to approximate Equation 21. We compare Figure 6: Lanczos vs diagonal
the evolution of the loss function over time and various uncertainty approx. for a Bayesian VAN.
calibration metrics. Figure 6 demonstrates training and Table 5 shows results.
Analysis: Lanczos uses matrix-vector products better (by a margin) The results suggest how, fora fixed matrix-vector-product budget, Lanczos achieves a drastically better likelihood at a similar
computational budget and already shows significant improvement with a much smaller budget.
Lanczos outperforms the diagonal approximation on all metrics except ECE. The subpar performanceof the diagonal approximation matches the observations of Ritter et al. [31]; see also [14]. The maintakeaway from this study is that differentiable matrix-free linear algebra unlocks new techniques for
Laplace approximations and allows further advances for Bayesian neural networks in general.
9Acknowledgments and Disclosure of Funding
This work was supported by a research grant (42062) from VILLUM FONDEN. The work was partly
funded by the Novo Nordisk Foundation through the Center for Basic Machine Learning Research inLife Science (NNF20OC0062606)

References[1] Nicholas J Higham. Functions of Matrices: Theory and Computation. SIAM, 2008.
[2] Efstratios Gallopoulos and Yousef Saad. Efficient solution of parabolic equations by Krylov
approximation methods. SIAM Journal on Scientific and Statistical Computing, 13(5):1236–
1264, 1992.
[3] Marlis Hochbruck and Alexander Ostermann. Exponential integrators. Acta Numerica, 19:
209–286, 2010.
[4] Changyi Xiao and Ligang Liu. Generative flows with matrix exponential. In International
Conference on Machine Learning, pages 10452–10461. PMLR, 2020.
[5] Patrik Axelsson and Fredrik Gustafsson. Discrete-time solutions to the continuous-timedifferential Lyapunov equation with applications to Kalman filtering. IEEE Transactions on
Automatic Control, 60(3):632–643, 2014.
[6] Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson.
GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. In
Advances in Neural Information Processing Systems, 2018.
[7] Jonathan Wenger, Geoff Pleiss, Philipp Hennig, John Cunningham, and Jacob Gardner. Preconditioning for scalable Gaussian process hyperparameter optimization. In International
Conference on Machine Learning. PMLR, 2022.
[8] Alexander Immer, Tycho FA Van Der Ouderaa, Mark Van Der Wilk, Gunnar Ratsch, andBernhard Schölkopf. Stochastic marginal likelihood gradients using neural tangent kernels. In
International Conference on Machine Learning, pages 14333–14352. PMLR, 2023.
[9] Jens Behrmann, Will Grathwohl, Ricky TQ Chen, David Duvenaud, and Jörn-Henrik Jacobsen.
Invertible residual networks. In International Conference on Machine Learning, pages 573–
582. PMLR, 2019.
[10] Ricky TQ Chen, Jens Behrmann, David K Duvenaud, and Jörn-Henrik Jacobsen. Residual
flows for invertible generative modeling. In Advances in Neural Information Processing
Systems, volume 32, 2019.
[11] Geoff Pleiss, Jacob Gardner, Kilian Weinberger, and Andrew Gordon Wilson. Constant-time
predictive distributions for Gaussian processes. In International Conference on Machine
Learning, pages 4114–4123. PMLR, 2018.
[12] Chris Finlay, Jörn-Henrik Jacobsen, Levon Nurbekyan, and Adam Oberman. How to train yourneural ODE: the world of Jacobian and kinetic regularization. In International Conference on
Machine Learning, pages 3154–3164. PMLR, 2020.
[13] Nicol N Schraudolph. Fast curvature matrix-vector products for second-order gradient descent.
Neural Computation, 14(7):1723–1738, 2002.
[14] Erik Daxberger, Agustinus Kristiadi, Alexander Immer, Runa Eschenhagen, Matthias Bauer,and Philipp Hennig. Laplace redux – effortless Bayesian deep learning. Advances in Neural
Information Processing Systems, 34:20089–20103, 2021.
[15] Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplaciansmoothing splines. Communications in Statistics-Simulation and Computation, 18(3):1059–
1076, 1989.
10

[16] Cornelius Lanczos. An iteration method for the solution of the eigenvalue problem of linear
differential and integral operators. 1950.
[17] Shashanka Ubaru, Jie Chen, and Yousef Saad. Fast estimation of tr(f (A)) via stochastic
Lanczos quadrature. SIAM Journal on Matrix Analysis and Applications, 38(4):1075–1099,
2017.
[18] Tyler Chen and Eric Hallman. Krylov-aware stochastic trace estimation. SIAM Journal on
Matrix Analysis and Applications, 44(3):1218–1244, 2023.
[19] Ke Wang, Geoff Pleiss, Jacob Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gordon Wilson. Exact Gaussian processes on a million data points. Advances in Neural Information
Processing Systems, 32, 2019.
[20] Kaiwen Wu, Jonathan Wenger, Haydn Jones, Geoff Pleiss, and Jacob R Gardner. Largescale Gaussian processes via alternating projection. International Conference on Artificial
Intelligence and Statistics, 2024.
[21] Josie König, Max Pfeffer, and Martin Stoll. Efficient training of Gaussian processes withtensor product structure. Preprint on ArXiv:2312.15305, 2023.
[22] Alexander James Davies. Effective implementation of Gaussian process regression for machinelearning. PhD thesis, University of Cambridge, 2015.
[23] Kun Dong, David Eriksson, Hannes Nickisch, David Bindel, and Andrew G Wilson. Scalable
log determinants for Gaussian process kernel learning. In Advances in Neural Information
Processing Systems, 2017.
[24] Geoff Pleiss, Martin Jankowiak, David Eriksson, Anil Damle, and Jacob Gardner. Fast matrix
square roots with applications to Gaussian processes and Bayesian optimization. Advances in
Neural Information Processing Systems, 33:22268–22281, 2020.
[25] Hao Wang. The Krylov subspace methods for the computation of matrix exponentials. 2015.
[26] Qinsheng Zhang and Yongxin Chen. Fast sampling of diffusion models with exponential
integrator. International Conference on Learning Representations, 2023.
[27] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In
International Conference on Machine Learning, pages 1530–1538. PMLR, 2015.
[28] Emiel Hoogeboom, Victor Garcia Satorras, Jakub Tomczak, and Max Welling. The convolution
exponential and generalized Sylvester flows. In Advances in Neural Information Processing
Systems, volume 33, pages 18249–18260, 2020.
[29] Marlis Hochbruck, Christian Lubich, and Hubert Selhofer. Exponential integrators for large
systems of differential equations. SIAM Journal on Scientific Computing, 19(5):1552–1574,
1998.
[30] Alexander Immer, Matthias Bauer, Vincent Fortuin, Gunnar Rätsch, and Khan Mohammad
Emtiyaz. Scalable marginal likelihood estimation for model selection in deep learning. In
International Conference on Machine Learning, 2021.
[31] Hippolyt Ritter, Aleksandar Botev, and David Barber. A scalable Laplace approximation forneural networks. In International Conference on Learning Representations, volume 6, 2018.
[32] Gene H Golub and Gérard Meurant. Matrices, Moments and Quadrature with Applications,volume 30. Princeton University Press, 2009.
[33] Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow
matching for generative modeling. In International Conference on Learning Representations,2023.
[34] Aditya Ramesh and Yann LeCun. Backpropagation for implicit spectral densities. In Advancesin Neural Information Processing Systems, 2018.
11

[35] Christopher KI Williams and Carl Edward Rasmussen. Gaussian Processes for Machine
Learning, volume 2. MIT Press Cambridge, MA, 2006.
[36] Kurt Cutajar, Michael Osborne, John Cunningham, and Maurizio Filippone. Preconditioning
kernel matrices. In International Conference on Machine Learning, pages 2529–2538. PMLR,
2016.
[37] Erik Daxberger, Eric Nalisnick, James U Allingham, Javier Antorán, and José Miguel
Hernández-Lobato. Bayesian deep learning via subnetwork inference. In International
Conference on Machine Learning, pages 2510–2521. PMLR, 2021.
[38] Agustinus Kristiadi, Matthias Hein, and Philipp Hennig. Being Bayesian, even just a bit, fixes
overconfidence in ReLu networks. In International Conference on Machine Learning, pages
5436–5446. PMLR, 2020.
[39] Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Mr Prabhat, and Ryan Adams. Scalable Bayesian optimization using
deep neural networks. In International Conference on Machine Learning, pages 2171–2180.
PMLR, 2015.
[40] John Denker and Yann LeCun. Transforming neural-net output levels to probability distributions. Advances in neural information processing systems, 3, 1990.
[41] Hippolyt Ritter, Aleksandar Botev, and David Barber. Online structured Laplace approximations for overcoming catastrophic forgetting. Advances in Neural Information Processing
Systems, 31, 2018.
[42] Wesley J Maddox, Gregory Benton, and Andrew Gordon Wilson. Rethinking parametercounting in deep models: Effective dimensionality revisited. arXiv preprint arXiv:2003.02139,
2020.
[43] Apoorva Sharma, Navid Azizan, and Marco Pavone. Sketching curvature for efficient out-ofdistribution detection for deep neural networks. In Uncertainty in Artificial Intelligence, pages
1958–1967. PMLR, 2021.
[44] Marco Miani, Frederik Warburg, Pablo Moreno-Muñoz, Nicki Skafte Detlefsen, and Søren
Hauberg. Laplacian autoencoders for learning stochastic representations. In Advances in
Neural Information Processing Systems, 2022.
[45] Marlis Hochbruck and Christian Lubich. On Krylov subspace approximations to the matrix
exponential operator. SIAM Journal on Numerical Analysis, 34(5):1911–1925, 1997.
[46] Simo Särkkä and Arno Solin. Applied Stochastic Differential Equations, volume 10. Cambridge
University Press, 2019.
[47] Severi Rissanen, Markus Heinonen, and Arno Solin. Generative modelling with inverse heatdissipation. In International Conference on Learning Representations, 2023.
[48] Cleve Moler and Charles Van Loan. Nineteen dubious ways to compute the exponential of amatrix. SIAM Review, 20(4):801–836, 1978.
[49] Cleve Moler and Charles Van Loan. Nineteen dubious ways to compute the exponential of amatrix, twenty-five years later. SIAM Review, 45(1):3–49, 2003.
[50] Awad H Al-Mohy and Nicholas J Higham. Computing the Fréchet derivative of the matrix
exponential, with an application to condition number estimation. SIAM Journal on Matrix
Analysis and Applications, 30(4):1639–1657, 2009.
[51] Yang Cao, Shengtai Li, Linda Petzold, and Radu Serban. Adjoint sensitivity analysis fordifferential-algebraic equations: The adjoint DAE system and its numerical solution. SIAM
Journal on Scientific Computing, 24(3):1076–1089, 2003.
[52] Moritz Diehl and Sébastien Gros. Numerical optimal control. Optimization in Engineering
Center (OPTEC), 2011.
12

[53] Christopher Rackauckas and Qing Nie. DifferentialEquations.jl–a performant and feature-richecosystem for solving differential equations in Julia. Journal of Open Research Software, 5(1):15–15, 2017.
[54] Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
[55] Chris Rackauckas, Mike Innes, Yingbo Ma, Jesse Bettencourt, Lyndon White, and Vaibhav Dixit. DiffEqFlux.jl–a julia library for neural differential equations. arXiv preprint
arXiv:1902.02376, 2019.
[56] Yingbo Ma, Vaibhav Dixit, Michael J Innes, Xingjian Guo, and Chris Rackauckas. A comparison of automatic differentiation and continuous sensitivity analysis for derivatives ofdifferential equation solutions. In 2021 IEEE High Performance Extreme Computing Conference (HPEC), pages 1–9. IEEE, 2021.
[57] Gene H Golub and Charles F Van Loan. Matrix Computations. JHU press, 2013.
[58] Magnus Rudolph Hestenes, Eduard Stiefel, et al. Methods of conjugate gradients for solving
linear systems, volume 49. NBS Washington, DC, 1952.
[59] Walter Edwin Arnoldi. The principle of minimized iterations in the solution of the matrix
eigenvalue problem. Quarterly of Applied Mathematics, 9(1):17–29, 1951.
[60] Matthias Seeger, Asmus Hetzel, Zhenwen Dai, Eric Meissner, and Neil D Lawrence. Autodifferentiating linear algebra. Preprint on ArXiv:1710.08717, 2017.
[61] Sebastian F Walter and Lutz Lehmann. Algorithmic differentiation of linear algebra functions with application in optimum experimental design (extended version). Preprint on
ArXiv:1001.1654, 2010.
[62] Andreas Griewank and Andrea Walther. Evaluating Derivatives: Principles and Techniques of
Algorithmic Differentiation. SIAM, 2008.
[63] Scott P Kolodziej, Mohsen Aznaveh, Matthew Bullock, Jarrett David, Timothy A Davis,
Matthew Henderson, Yifan Hu, and Read Sandstrom. The SuiteSparse matrix collection
website interface. Journal of Open Source Software, 4(35):1244, 2019.
[64] Timothy A Davis and Yifan Hu. The University of Florida sparse matrix collection. ACM
Transactions on Mathematical Software (TOMS), 38(1):1–25, 2011.
[65] Iain S Duff, Roger G Grimes, and John G Lewis. Sparse matrix test problems. ACM
Transactions on Mathematical Software (TOMS), 15(1):1–14, 1989.
[66] Mathieu Blondel and Vincent Roulet. The elements of differentiable programming. arXiv
preprint arXiv:2403.14606, 2024.
[67] Gerrit C Groenenboom and Henk M Buck. Solving the discretized time-independent
Schrödinger equation with the Lanczos procedure. The Journal of Chemical Physics, 92(7):4374–4379, 1990.
[68] Christopher Conway Paige. The computation of eigenvalues and eigenvectors of very large
sparse matrices. PhD thesis, University of London, 1971.
[69] Christopher C Paige. Computational variants of the Lanczos method for the eigenproblem.
IMA Journal of Applied Mathematics, 10(3):373–381, 1972.
[70] Christopher C Paige. Error analysis of the Lanczos algorithm for tridiagonalizing a symmetricmatrix. IMA Journal of Applied Mathematics, 18(3):341–349, 1976.
[71] Steffen Börm and Christian Mehl. Numerical methods for Eigenvalue Problems. Walter deGruyter, 2012.
[72] Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe LlinaresLópez, Fabian Pedregosa, and Jean-Philippe Vert. Efficient and modular implicit differentiation.
Advances in Neural Information Processing Systems, 35:5230–5242, 2022.
13

[73] Michael Betancourt, Charles C Margossian, and Vianey Leos-Barajas. The discrete adjoint method: Efficient derivatives for functions of discrete sequences. Preprint on
ArXiv:2002.00326, 2020.
[74] Brian Kha Tran and Melvin Leok. Geometric methods for adjoint systems. Journal of
Nonlinear Science, 34(1):25, 2024.
[75] Denisa AO Roberts and Lucas R Roberts. QR and LQ decomposition matrix backpropagation algorithms for square, wide, and deep–real or complex–matrices and their software
implementation. Preprint on ArXiv:2009.10071, 2020.
[76] Thomas P Minka.
Old and new matrix algebra useful for statistics.
www.stat.cmu.edu/minka/papers/matrix.html, 4, 2000.

See[77] Mike Giles. An extended collection of matrix derivative results for forward and reverse mode
automatic differentiation. 2008.
[78] Jean Céa. Conception optimale ou identification de formes, calcul rapide de la dérivéedirectionnelle de la fonction coût. ESAIM: Modélisation mathématique et analyse numérique,20(3):371–402, 1986.
[79] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal
Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, andQiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL
http://github.com/google/jax.
[80] Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian
processes (KISS-GP). In International Conference on Machine Learning, pages 1775–1784.
PMLR, 2015.
[81] Helmut Harbrecht, Michael Peters, and Reinhold Schneider. On the low-rank approximationby the pivoted Cholesky decomposition. Applied Numerical Mathematics, 62(4):428–440,
2012.
[82] Benjamin Charlier, Jean Feydy, Joan Alexis Glaunes, François-David Collin, and Ghislain
Durif. Kernel operations on the GPU, with autodiff, without memory overflows. Journal of
Machine Learning Research, 22(74):1–6, 2021.
[83] Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. San Diega,
CA, USA, 2015.
[84] J. R. Dormand and P. J. Prince. A family of embedded Runge–Kutta formulae. J. Comp. Appl.
Math, 6:19–26, 1980.
[85] Lawrence F. Shampine. Some practical Runge-Kutta formulas. Mathematics of Computation,
46(173):135–150, 1986. doi: https://doi.org/10.2307/2008219.
[86] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary
differential equations. Advances in Neural Information Processing Systems, 31, 2018.
[87] Charalampos Tsitouras. Runge–Kutta pairs of order 5 (4) satisfying only the first column
simplifying assumption. Computers & Mathematics with Applications, 62(2):770–775, 2011.
[88] Meng-Hao Guo, Cheng-Ze Lu, Zheng-Ning Liu, Ming-Ming Cheng, and Shi-Min Hu. Visual
attention network. Computational Visual Media, 9(4):733–752, 2023.
[89] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. ImageNet: A large-scale
hierarchical image database. In IEEE Conference on Computer Vision and Pattern Recognition,
pages 248–255. IEEE, 2009.
[90] Zhijie Deng, Feng Zhou, and Jun Zhu. Accelerated linearized Laplace approximation forBayesian deep learning. Advances in Neural Information Processing Systems, 35:2695–2708,
2022.
14

[91] Bolei Zhou, Agata Lapedriza, Jianxiong Xiao, Antonio Torralba, and Aude Oliva. Learning
deep features for scene recognition using places database. Advances in Neural Information
Processing Systems, 27, 2014.
[92] Tianlin Liu, Jose Antonio Lara Benitez, Florian Faucher, AmirEhsan Khorashadizadeh,
Maarten V de Hoop, and Ivan Dokmanić. WaveBench: Benchmarks datasets for modeling linear wave propagation PDEs. Transactions on Machine Learning Research, 2023.
[93] Peter Deuflhard, Andreas Hohmann, Peter Deuflhard, and Andreas Hohmann. Three-termrecurrence relations. Numerical Analysis in Modern Scientific Computing: An Introduction,
pages 151–178, 2003.
[94] Peter Deuflhard. On algorithms for the summation of certain special functions. Computing,
17:37–48, 1976.
[95] Simon Bartels, Kristoffer Stensbo-Smidt, Pablo Moreno-Muñoz, Wouter Boomsma, Jes
Frellsen, and Soren Hauberg. Adaptive Cholesky Gaussian processes. In International
Conference on Artificial Intelligence and Statistics, pages 408–452. PMLR, 2023.
[96] Rui Camachol. Inducing models of human control skills. In Machine Learning: ECML-98:
10th European Conference on Machine Learning Chemnitz, Germany, April 21–23, 1998
Proceedings 10, pages 107–118. Springer, 1998.
[97] Anton Schwaighofer and Volker Tresp. Transductive and inductive methods for approximate
Gaussian process regression. Advances in Neural Information Processing Systems, 15, 2002.
[98] Paul Shannon, Andrew Markiel, Owen Ozier, Nitin S Baliga, Jonathan T Wang, Daniel Ramage,
Nada Amin, Benno Schwikowski, and Trey Ideker. Cytoscape: a software environment for
integrated models of biomolecular interaction networks. Genome research, 13(11):2498–2504,
2003.
[99] P. J Prince and J. R. Dormand. High order embedded Runge–Kutta formulae. Journal ofComputational and Applied Mathematics, 7(1):67–75, 1981.
[100] Antonio Stanziola, Simon R Arridge, Ben T Cox, and Bradley E Treeby. j-Wave: An opensource differentiable wave simulator. SoftwareX, 22:101338, 2023.
[101] David JC MacKay. Bayesian interpolation. Neural computation, 4(3):415–447, 1992.
[102] Alexander Immer, Maciej Korzepa, and Matthias Bauer. Improving predictions of Bayesianneural nets via local linearization. In International Conference on Artificial Intelligence andStatistics, pages 703–711. PMLR, 2021.
[103] Vardan Papyan. Traces of class/cross-class structure pervade deep learning spectra. Journal of
Machine Learning Research, 21(252):1–64, 2020.
[104] Costas Bekas, Effrosyni Kokiopoulou, and Yousef Saad. An estimator for the diagonal of amatrix. Applied Numerical Mathematics, 57(11-12):1214–1229, 2007.
[105] Mahdi Pakdaman Naeini, Gregory Cooper, and Milos Hauskrecht. Obtaining well calibrated
probabilities using Bayesian binning. In Proceedings of the AAAI Conference on Artificial
Intelligence, volume 29, 2015.

15

Appendix: Overview
Some of the results in the main paper promised detailed information about setups, data, compute, oradditional proofs. For example, the case study about partial differential equations involves a data
generation process which will receive further explanation in this supplement.
The appendix provides the following details: Appendices A and F elaborate on Figure 3 and Table 2
respectively; Appendices B to E contain proofs for the main results; and Appendices G to J describe
the setup used for the case studies. Notably, Appendix G describes how we implement low-memory
matrix-vector products in JAX (to replicate what makes libraries like KeOps [82] so efficient), andAppendix I outlines a PDE data set similar to that by Liu et al. [92].
Code Most of the contributions of this paper pertain to differentiable implementations of numericalalgorithms – at the heart of it are our reverse-mode differentiable Lanczos and Arnoldi iterations. We
provide JAX code to reproduce all experiments at the URL
https://github.com/pnkraemer/experiments-lanczos-adjoints
and have packaged all numerical methods in a JAX library that can be installed via
pip install matfree
Next to Lanczos and Arnoldi, this includes variants of conjugate gradient methods [58] and pivoted
Cholesky preconditioners [81] (which we used for Gaussian processes), low-memory kernel-matrixvector products (also Gaussian processes), efficient GGN-vector products (used for Bayesian neuralnetworks), and matrix-free sampling algorithms to sample from Gaussian processes and Laplaceapproximated Bayesian neural networks (used for Gaussian processes and Bayesian neural networks).
These methods are known in some form or another, but until now, they have all lacked a software
implementation in the current JAX ecosystem (with some exceptions relating to conjugate gradients).
Compute All experiments before the case studies were run on CPU. The Gaussian process anddifferential equation case studies run on a V100 GPU, the Bayesian neural network one on a P100
GPU. The Gaussian process and Bayesian neural network studies run in a few hours, all other code
finishes in a few minutes.

A

Additional context for Figure 3

To create Figure 3, we load the “bcsstk18” matrix from the SuiteSparse matrix collection [63–65]. This matrix is symmetric and has 11,948 rows/columns and 149,090 nonzero entries. We
implement matrix-vector products with jax.experimental.sparse, and use a Lanczos iteration
without reorthogonalisation. We time the execution of the forward pass, as well as the backward
pass with and without implementing a custom vector-Jacobian product (the custom vector-Jacobian
product involves the adjoint). The results were shown in Figure 3, and displayed how rapidly
the computational complexity of automatic differentiation increases, whereas the computational
complexity of our custom gradient mirrors that of the forward pass. This figure showed how without
our proposed gradients, differentiating through the Lanczos iteration is unfeasible.
To add to this benchmark, we repeat the same for the Arnoldi iteration and show both compilation andruntime for Lanczos and Arnoldi in Figure 7 (this includes the curves from Figure 3 again). We see
Figure 7: Run and compilation times for Lanczos and Arnoldi. The “backprop” curves were stopped
at 100, because for higher values we encountered memory issues. All experiments run on CPU.

0.0

ForwardAdjoint
Backprop050

100 150 200 250

Krylov-space depth

Arnoldi: Run time

ForwardAdjoint
Backprop

10
5
0

50

100 150 200 250

Krylov-space depth

ForwardAdjoint
Backprop

1.5
1.0
0.5
0.0
0

50

100 150 200 250

Krylov-space depth

Arnoldi: Compile time

15

Wall time (sec)

0.1

15

Wall time (sec)Wall time (sec)0.2

Lanczos: Compile time
Wall time (sec)Lanczos: Run time

0.3

ForwardAdjoint
Backprop

10
5
0

50

100 150 200 250

Krylov-space depth

different behaviour for Lanczos and Arnoldi. Whereas for back-propagation through Lanczos without
16

the custom gradient, the compilation times remain constant for increasing Krylov-space depth K andruntimes increase rapidly, the reverse is true for the Arnoldi iteration. The adjoint method mirrorsthat of the forward pass in both benchmarks. In either case, the increasing memory requirements forbackpropagation through Lanczos and Arnoldi without our proposed adjoints becomes apparent.

B

Proof of Theorem 4.1

Here is how we derive the adjoints of the Arnoldi system. The structure is the usual: take the forward
constraint, differentiate, add Lagrange multipliers (“transpose”), and identify the adjoint system.
B.1

Linearisation

The linearisation of Equation 12 is(dA)Q + AdQ − (dQ)H − QdH − (dr)(eK )⊤ = 0(dQ)e1 − (dv)c − vdc = 0∈ RN ×K ,
∈R

N ×1

(22a)

,

(22b)

I≤ ◦ [⟨Qei , (dQ)ej ⟩ + ⟨Qej , (dQ)ei ⟩]K
i,j=1 = 0

∈ RK×K ,(22c)

I≪ ◦ dH = 0

K×K

(22d)

∈R

[⟨r, (dQ)ej ⟩ + ⟨Qej , dr⟩]K
j=1 = 0

,

K×1

∈R
.
To see this, apply the chain- and product rules to the original constraint in Equation 12.
B.2

(22e)

Transposition

Let ρ = ρ(Q, r, H, c) ∈ R be a scalar function of the outputs. In the following, interpret vectors as
N × 1 matrices and scalars as 1 × 1 matrices. The values of inner products and the realisations of ρ
are the only scalars.
For all Λ ∈ RN ×K , λ ∈ RN ×1 , Γ, Σ ∈ RK×K , and γ ∈ RK×1 , we havedρ = ⟨∇Q ρ, dQ⟩ + ⟨∇H ρ, dH⟩ + ⟨∇r ρ, dr⟩ + ⟨∇c ρ, dc⟩= ⟨∇Q ρ, dQ⟩ + ⟨∇H ρ, dH⟩ + ⟨∇r ρ, dr⟩ + ⟨∇c ρ, dc⟩(23a)+ ⟨Λ, (dA)Q + AdQ − (dQ)H − QdH − (dr)(eK )⊤ ⟩+ ⟨λ, (dQ)e1 − (dv)c − vdc⟩
+ ⟨Γ, I≤ ◦ [⟨Qei , (dQ)ej ⟩ + ⟨Qej , (dQ)ei ⟩]K
i,j=1 ⟩+ ⟨Σ, I≪ ◦ dH⟩+ ⟨γ, [⟨r, (dQ)ej ⟩ + ⟨Qej , dr⟩]K
j=1 ⟩

(23b)

= ⟨∇Q ρ, dQ⟩ + ⟨∇H ρ, dH⟩ + ⟨∇r ρ, dr⟩ + ⟨∇c ρ, dc⟩+ ⟨ΛQ⊤ , dA⟩ + ⟨A⊤ Λ, dQ⟩ − ⟨ΛH ⊤ , dQ⟩ − ⟨Q⊤ Λ, dH⟩ − ⟨ΛeK , dr⟩+ ⟨λ(e1 )⊤ , dQ⟩ − ⟨λc⊤ , dv⟩ − ⟨v ⊤ λ, dc⟩+ ⟨Q(I≤ ◦ Γ), dQ⟩ + ⟨Q(I≤ ◦ Γ)⊤ , dQ⟩+ ⟨I≪ ◦ Σ, dH⟩+ ⟨rγ ⊤ , dQ⟩ + ⟨Qγ, dr⟩

(23c)
⊤

⊤

=: ⟨ZQ , dQ⟩ + ⟨ZH , dH⟩ + ⟨Zr , dr⟩ + ⟨Zc , dc⟩ + ⟨ΛQ , dA⟩ + ⟨λc , dv⟩
with the constraints
ZQ := ∇Q ρ + A⊤ Λ − ΛH ⊤ + λ(e1 )⊤ + Q(I≤ ◦ Γ) + Q(I≤ ◦ Γ)⊤ + rγ ⊤ ∈ RN ×K⊤

ZH := ∇H ρ − Q Λ + I≪ ◦ Σ
Zr := ∇r ρ − ΛeK + Qγ
⊤

(23d)
(24a)

∈R

K×K

(24b)

∈R

N ×1

(24c)

1×1

Zc := ∇c ρ − v λ∈ R . (24d)
Solving the adjoint system, ZQ = 0, ZH = 0, Zr = 0, and Zc = 0 as a function of Λ, λ, Γ, Σ, and γ,
yields the desired ∇A ρ = ΛQ⊤ and ∇v f = λc⊤ . Theorem 4.1 is complete.
17

C

Proof of Theorem 4.2

Recall that ρ = ρ(x1 , ..., xK+1 ; a1 , ..., aK ; b1 , ..., bK ) shall be a scalar/loss that depends on the
output of the algorithm. Denote by ∇xk ρ the gradient of ρ with respect to each Lanczos vector xk ,and by ∇ak ρ and ∇bk ρ the gradients with respect to ak and bk respectively.
The differential of normalisation, i.e., the operation s 7→ h = s/(s⊤ s) is1
dh = ⊤ I − hh⊤ ds.
(25)
s s
The next steps are the usual ones: we start with the forward constraint, linearise, add Lagrange
multipliers, and identify the adjoint system. The order of the middle two steps (linearise, multipliers)
is interchangeable; while for Arnoldi, we linearise first and then add Lagrange multipliers, for
Lanczos, we go the other way.
n
K
Define Lagrange multipliers {λk }K
k=0 ⊆ R and {µk /2, νk }j,k=1 ⊆ R,

ρ=ρ+

K
X

λ⊤
k (−bk−1 xk−1 + (A − ak I)xk − bk xk+1 )(b0 = 1, x0 = 0)

k=1

− λ⊤
0


x1 − √



v
v⊤ v

+

K

K

k=1

k=1

X
1X
νk x⊤
µk (x⊤
k xk+1 .
k+1 xk+1 − 1) +2

(26)

Differentiate and use that the forward constraint must be satisfied,
dρ =

K+1
X

(∇xk ρ)⊤ dxk +

k=1

+

K
X

(∇ak ρ)⊤ dak +

k=1

K
X

K
X

(∇bk ρ)⊤ dbk

k=1

λ⊤k [(dA)xk − (dak )xk − (dbk−1 )xk−1 − (dbk )xk+1 ](db0 = 1, dx0 = 0)

k=1

+

K
X

λ⊤
k [Adxk − ak dxk − bk−1 dxk−1 − bk dxk+1 ](b0 = 1, x0 = 0)

k=1
⊤
⊤
− λ⊤
0 (dx1 − (I − x1 x1 )/(v v)dv)

+

K
X

µk x⊤
k+1 dxk+1

k=1

+

K
X

⊤
νk (x⊤
k+1 dxk + xk dxk+1 ).

(27)

k=1

Sort all terms by differential,
dρ =

K
X

Zak dak +

k=1

K
X

Zbk dbk +

k=1

K+1
X

Zxk dxk + Zv dv + trace (ZA dA) .

(28)

k=1

so that enforcing that all Zak , Zbk , Zxk terms are zero yields constraints for the multipliers fromwhich we can compute gradients with respect to v and A.
What are those terms? Let λK+1 = 0, µ0 = 0, and ν0 = 0 (to simplify notation below); then,
ZxK+1 = −λK bK + (∇xK+1 ρ + µK xK+1 + νK xK ),and for all k = K, ...1,(29)(recall x0 = 0, b0 = 1, µ0 = 0, ν0 = 0, λK+1 = 0)

Zxk = −bk λk+1 + (A⊤ − ak I)λk − bk−1 λk−1 + (∇xk ρ + µk−1 xk + νk xk+1 + νk−1 xk−1 ),(30a)
Zak = ∇ak ρ − λ⊤k xk ,(30b)⊤
Zbk = ∇bk ρ − λ⊤k+1 xk − λk xk+1 .

(30c)
18

The expressions are like the forward-Lanczos constraints, and the main differences are (i) that therecursions are run backwards in “time” and (ii) the existence of a nonzero bias term in the adjoints(marked by parentheses). Enforcing Zak , Zbk , Zxk to be zero identifies
"K
#
 ⊤

X
λ0 x1 ⊤∇v ρ = Zv :=x − λ⊤, ∇A ρ = ZA :=
xk λ⊤
(31)
0
k .
v⊤ v 1
k=1

Theorem 4.2 is complete.

D

Proof of Corollary 4.3

The following derivation covers only the case for Lanczos, i.e., we use the variables x1 , ..., xK+1
instead of Arnoldi’s q1 , ..., qK . But the derivation is the same for both methods.
The expression we manipulate is"KXdρ = trace

#
xk λ⊤
k

!
dA

+ const(32)

k=1

where all non-dA-related quantities are treated as some unimportant constants.
For a single parameter θj , dA = (∇θj A)⊤ dθj and we have
"K
#
!
"K
#
!
X
X
⊤
⊤
⊤
xk λk Dj Adθj = tracexk λk (∇θj A)
dθj
trace
k=1

(chain rule)k=1= trace

K
X

!
⊤
λ⊤
k (∇θj A) xkdθj (cyclic property of traces)

k=1

=

K
X

⊤
λ⊤
k (∇θj A) xk dθj(a scalar is its own trace)k=1= ∇θ j

"K
X

#
⊤
λ⊤
k A xkdθj(linearity of diff. & summation)dθj .

(symmetry of A)

k=1= ∇θ j

"K
X

#
λ⊤
k Axk

k=1

In conclusion, the derivative of ρ wrt θj is∇θ j ρ = ∇j

"K
X

#
⊤
x⊤
k A λk(33)

k=1

and stacking all of those partial derivatives on top of each other, we obtain
#
"
K
K
X
X


⊤
⊤
∇θ ρ = [∇θj ρ]j = ∇ θ 7→
xk A(θ) λk =∇ θ 7→ A(θ)⊤ λk xk .
k=1

(34)

k=1

We already compute A(θ)⊤ λk during the backward pass, so we are a single vector-Jacobian productwith xk away from a matrix-parameter-gradient instead of a matrix-gradient. This requires O(p + n)
storage, and is computed online, which makes the memory-complexity independent of K.

E

Solving the adjoint system

The upcoming section details how to solve the adjoint system for both, the Lanczos and the Arnoldi
iterations. It reuses notation from Appendices B and C.
We begin with Lanczos, because the solution is less technical, and because starting with Lanczos canprovide a template for solving Arnoldi’s adjoint system. All results in the present section (except for
those that explicitly point to Deuflhard et al. [93], Deuflhard [94], which are marked as such) are newand a contribution of this work.
19

E.1

LanczosK
The inputs to the adjoint system are the Lanczos vectors {xk }K+1
k=1 and the coefficients {ak }k=1 as
K+1K
well as {ak }k=1 from the forward pass, the corresponding input derivatives {∇xk ρ}k=1 , {∇ak ρ}K
k=1 ,and {∇ak ρ}K,andmatrixA
andinitialvectorv.
k=1

The overall strategy for solving the adjoint system of the Lanczos iteration (Theorem 4.2) is thefollowing: for every k = K, ..., 1, alternate the two steps:
1. Combine orthogonality with the Zak constraints to get νk , and combine it with the Zbk
constraints to get µk .
2. Once each µk and νk are available, solve for λk and repeat with the next lower k.
This results in the following procedure: To start, set ζK+1 = −(∇xK+1 ρ) and λK+1 = 0. Then, forall k = K, ..., 1, computeξk = ζk+1 /bk

(35a)

⊤
µ̃k = ∇bk ρ − λ⊤k+1 xk + xk+1 ξkν̃k = ∇ak ρ + x⊤
k ξk

(35b)
(35c)
(35d)

λk = −ξk + µ̃k · xk+1 + ν̃k · xk
ζk = −∇xk ρ − A⊤ λk + ak · λk + bk · λk+1 − bk · ν̃k · xk+1
Repeat with k = k − 1.

(35e)
(35f)

Finally, set λ0 = ζ1 . Only λ and ζ affect subsequent steps; µ, ν, and ξ are only needed for computing
λ and ζ. The k-th step depends on ak , bk , xk , xk+1 , ∇ak ρ, ∇bk ρ, ∇xk ρ.
The strategy above yields all {λk }K
k=0 . Finalise the gradients∇v ρ =

λ⊤0 x1x1 − λ0 ,v⊤ v∇A ρ =KX

λ k x⊤
k

(36)

k=1

which can be embedded into any reverse-mode algorithmic-differentiation-engine.
E.2

ArnoldiThe process for Arnoldi is similar to that for Lanczos, but the derivation is more technical. It shares
many similarities with deriving the Arnoldi iteration (i.e., the forward pass), so we begin by providing
a perspective on recurrence relations (which include the Arnoldi iteration) through the lens of linear
system solvers [93, 94] before we use this perspective to solve the adjoint system.
E.2.1

Solving the original system

At its core, finding the Arnoldi vectors amounts to solving−AQ + QH + r(eK )⊤ = 0,
Qe1 − cv = 0,

(37a)
(37b)

which is possible in closed form as follows: We rewrite the first two constraints as(e1 ⊗ I)vec (Q) = vec (cv)⊤−(I ⊗ A)vec (Q) + (H ⊗ I)vec (Q) + (eK ⊗ I)vec (r) = 0.
This expression is equivalent to
e1 ⊗ IH⊤ ⊗ I − I ⊗ A

0eK ⊗ I 
vec (Q)vec (cv)∈ RN (K+1)×1=
vec (r)
0

(38a)
(38b)

(39)

This is similar to the work by Deuflhard [94], who explain adjoints of three-term recurrence relations.
Since c2 = ⟨v, v⟩ holds (ie, c is known), the first row of Q is known. Then, the first row of Q together
with the orthogonality constraints yields the first row of H ⊤ , which then defines the next row of
the linear system thus the second row of Q. Alternating between deriving the next row of H ⊤ andsolving the lower triangular systems is then Arnoldi’s algorithm:
20

Algorithm
E.1 (Arnoldi’s forward pass; paraphrased). Assume that v and K are known. Computepc = ⟨v, v⟩. Then, for k = 1, ..., K, alternate the following tsteps:
1. Derive the next column of H using the orthogonality constraints.
2. Forward-substitute (“solve”) the block-lower-triangular system for the next column of Q
(respectively r at the last iteration).
Return all Q and H, as well as c and r.
The same principle applies to the adjoint system, and the only difference is that the notation is slightly
more complicated:
E.2.2

Solving the adjoint system

The constraints ZQ = 0 and Zr = 0 mirror those of AQ − QH − r(eK )⊤ = 0 and Qe1 − cv = 0;
the constraints ZH = 0 and Zc = 0 mirror the orthogonality constraints Q⊤ Q = I and Q⊤ r = 0.
Therefore, we start with ZQ = 0 and Zr = 0,
∇Q f + A⊤ Λ − ΛH ⊤ + λ(e1 )⊤ + Q(I≤ ◦ Γ) + Q(I≤ ◦ Γ)⊤ + rγ ⊤ = 0∇r f − ΛeK + Qγ = 0.

(40a)
(40b)

Introduce the auxiliary quantities
Ψ(Γ, γ) := ∇Q f + Q(I≤ ◦ Γ) + Q(I≤ ◦ Γ)⊤ + rγ ⊤ψ(γ) := ∇r f + Qγ

∈ RN ×K

(41a)

N ×1

(41b)

∈R

.

Ψ and ψ only serve the purpose of simplifying the notation in the coming part; there is no “meaning”
associated with them. Vectorise both expressions,(I ⊗ A⊤ )vec (Λ) − [(H ⊤ )⊤ ⊗ I]vec (Λ) + (e1 ⊗ I)vec (λ) = −vec (Ψ(Γ, γ))⊤[(eK ) ⊗ I]vec (Λ) = vec (ψ(γ))
and observe that this can be written as a linear system


 

e1 ⊗ I I ⊗ A⊤ − (H ⊤ )⊤ ⊗ I
vec (λ)−vec (Ψ(Γ, γ))=vec (Λ)vec (ψ(γ))
0(eK )⊤ ⊗ I

(42a)
(42b)

(43)

with a system matrix that is the transpose of the system matrix of the forward pass. The matrix
is upper triangular, and the equation can be solved with backward substitution provided ψ(γ) andΨ(Γ, γ) are known.
The defining quantities Ψ and ψ emerge by combining the adjoint recursion with the projection
constraints ZH = 0 (for Γ, which yields Ψ) and Zr = 0 (for γ, which yields ψ). We use Zc = 0 to
get a single element in Γ; more on this below. Summarise the adjoint pass:
Algorithm E.2 (Arnoldi’s adjoint pass; paraphrased). Assume Q, H, c, and r as well as the gradientsof f with respect to those quantities. Then, compute ψ via computing γ using Zr = 0. Then, fork = K, ..., 1, alternate the following two steps:
1. Derive the next row of Ψ by combining ZQ = 0 with the projection constraint ZH = 0
2. Backward-substitute (“solve”) for the next row of Λ (recall: we loop backwards)
Finally, use Zc = 0 to get the first row of Ψ and solve for λ. Then, return ∇A f = ΛQ⊤ and∇v f = λc⊤ .
The structure of the adjoint pass is similar to the forward pass (Table 6). In the following, we will
elaborate on each of those steps. We assume that the reader knows how to solve a lower triangularlinear system. We focus on constructing ψ and Ψ via Γ and γ.
21

Table 6: Forward versus adjoint (backward) pass

ForwardAdjoint

E.2.3

System matrix

Solve via

Recursively define

Using

Lower triangular
Upper triangular

Forward substitution
Backward substitution

System matrix
Right-hand side

Orthogonality
Projection: ZH = 0

Initialisation

Initialisation of the adjoint pass implies computing ψ. To get ψ, we need γ: Consider multiplying Zr
with Q⊤ ,0 = Q⊤ Zr = Q⊤ ∇r f − Q⊤ ΛeK + γ(use the definition of Zr )⊤(since ZH = 0)⊤(reorder)⊤((I≪ ◦ Σ)eK = 0)

= Q ∇r f − (∇H f + I≪ ◦ Σ) eK + γ= (Q ∇r f − ∇H f eK ) + (I≪ ◦ Σ)eK + γ= (Q ∇r f − ∇H f eK ) + γ

where we use that the last column in I≪ ◦ Σ consists entirely of zeros. All other quantities are known.
Therefore, γ is isolated and we identify
γ = ∇H f eK − Q⊤ ∇r f.

(44)

ψ(γ) = ∇r f + Qγ

(45)Next, use this γ to build ψ,and the initialisation step is complete.
E.2.4

Recursion

With ψ in place, we get the last column of Λ. To get the next column of Λ, we need to derive theright-hand side Ψ. To get Ψ, we need Γ.
Multiply Q⊤ ZQ to obtain
0 = Q⊤ ZQ
⊤

(46)
⊤

⊤

⊤

⊤

⊤

⊤

⊤

= Q ∇Q f + Q A Λ − Q ΛH + Q λ(e1 ) + I≤ ◦ Γ + (I≤ ◦ Γ)
(since Q⊤ Q = I and Q⊤ r = 0)
= Q⊤ ∇Q f + Q⊤ A⊤ Λ − (∇H f + I≪ ◦ Σ)H ⊤ + Q⊤ λ(e1 )⊤ + I≤ ◦ Γ + (I≤ ◦ Γ)⊤(since ZH = 0)
= (Q⊤ ∇Q f − ∇H f H ⊤ ) + Q⊤ A⊤ Λ + Q⊤ λ(e1 )⊤ + I≤ ◦ Γ + (I≤ ◦ Γ)⊤ + (I≪ ◦ Σ)H ⊤ .
(reorder the terms)
Since H is Hessenberg, (I≪ ◦ Σ)H ⊤ is strictly lower triangular. Therefore, multiplication with I≥
removes Σ from the expression,
0 = I≥ ◦ Q⊤ ∇Q f − ∇H f H ⊤ + Q⊤ A⊤ Λ + Q⊤ λ(e1 )⊤ + I≤ ◦ Γ + (I≤ ◦ Γ)⊤
(47a)
 ⊤

⊤
⊤ ⊤
⊤
⊤
⊤
= I≥ ◦ Q ∇Q f − ∇H f H + Q A Λ + I≥ ◦ [Q λ(e1 ) ] + I= ◦ Γ + I≥ ◦ Γ . (47b)
The term involving λ can be simplified as follows: Due to the presence of (e1 )⊤ , we know thatQ⊤ λ(e1 )⊤ is lower triangular (in fact, it has a single nonzero column). Thus, I≥ ◦ (Q⊤ λ(e1 )⊤ ) is
proportional to e1 (e1 )⊤ ,I≥ ◦ (Q⊤ λ(e1 )⊤ ) = [(Qe1 )⊤ λ] e1 (e1 )⊤
⊤

(48)

⊤

(since Qe1 = cv)⊤(since Zc = 0)

= cv λ e1 (e1 )= c ∇c f e1 (e1 )and all quantities are known; hence,

I= ◦ Γ + I≥ ◦ Γ⊤ = −I≥ ◦ [Q⊤ ∇Q f − ∇H f H ⊤ + Q⊤ A⊤ Λ] − c ∇c f e1 (e1 )⊤
22

(49)

Algorithm E.3 (Forward pass). Initialise k =1, Qe1 . Then, for k = 1, ..., K:
1. Use orthogonality for a new row inthe system matrix in Equation 39.
2. Solve for the next column of Q
3. Optional: re-enforce Q⊤ Q = I.
Solve for r and return Q, H, r, c.

Algorithm E.4 (Backward pass). Initialise
k = K, Λe1 . Then, for k = K, ..., 1:
1. Use projection for a new column ofthe right-hand side Ψ
2. Solve for the next column of Λ
3. Optional: re-enforce ZH = 0.
Solve for λ and return ∇θ ρ and ∇v ρ.

Figure 8: Forward and backward pass of the Arnoldi iteration (paraphrased)must hold.
Now, the most important observation is the following: the last column of I= ◦ Γ + I≥ ◦ Γ⊤ depends on
the last column of Q⊤ A⊤ Λ and known quantities; the penultimate column depends on the penultimate
column of Γ, and so on. But at the time of assembling the last column of Γ, the last column of Λ is
known! More generally, we always know one more column of Λ than of Γ, so we can recursively
assemble I= ◦ Γ + I≥ ◦ Γ⊤ :Letsym (M ) := I≥ ◦ M + (I> ◦ M )⊤(50)

be a symmetrisation operator. We define it for the sole purpose of reconstructing

sym I= ◦ Γ + I≥ ◦ Γ⊤ = I≤ ◦ Γ + (I≤ ◦ Γ)⊤ .

(51)

Let us use it:I≤ ◦ Γ + (I≤ ◦ Γ)⊤ = sym −I≥ ◦ [Q⊤ ∇Q f − ∇H f H ⊤ + Q⊤ A⊤ Λ] − c ∇c f e1 (e1 )⊤ (52a)= sym −I≥ ◦ [Q⊤ ∇Q f − ∇H f H ⊤ + Q⊤ A⊤ Λ] − c ∇c f e1 (e1 )⊤ . (52b)This yields the next row/column of Γ + Γ⊤ , and therefore the next row of Ψ. From there, we canassemble the next column of Λ and iterate. Figure 8 (respectively Algorithms E.3 and E.4) compare
pseudocode for forward and adjoint passes. Altogether, the implementation of the adjoint pass is very
similar to that of the forward pass.
At the final step, we obtain not the last column of Λ but λ, though this is a byproduct of solving thetriangular linear system. It does not need further explanation.
Remark E.5 (Σ). Like for gradients of QR decompositions [61, 75], we never solve for Σ.

F

Setup for Table 2
I : A 7→ (H, Q, r, c) 7→ QHQ⊤

(53)

where (H, Q, r, c) are the result of a full-rank Arnoldi
iteration (i.e. K = N ). For K = N , QHQ⊤ = A and I
must have an identity Jacobian; thus,ε := ∥IN 2 − I∥RMSE(54)

measures the loss of accuracy when differentiating theArnoldi iteration. A small ε is desirable.

Loss of accuracy

To create Table 2, we implement an operator
Adj. w/ re-proj.
Adj. w/o re-proj.
Backprop

10 2
10 6
10 10
10 14
10 18

2

4

6

Hilbert matrix size

8

Figure 9: Accuracy loss ε when differenThen, using double-precision, we construct a Hilbert
tiating I for a Hilbert matrix of increasN
N ×N
matrix A = [1/(i + j + 1)]i,j=1 ∈ R
which is a fa- ing size N . Uses double precision.
mously ill-conditioned matrix and a common test-bed forthe loss of orthogonality in methods like the Lanczos and Arnoldi iteration [e.g. 71, Table 7.1]. We
evaluate three algorithms, all of which rely on the Arnoldi iteration with full reorthogonalisation on
the forward-pass: One algorithm does not re-project on the adjoint constraints, another one does, and23

Figure 10: For matrices with at least 10,000 rows/columns, KeOps remains the state of the art. This
experiment uses a square-exponential kernel, on an artificial dataset with d = 3 dimensions.
for reference we compute ε when “backpropagating through” the re-orthogonalised Arnoldi iteration
as a third option. Figure 3 has demonstrated that the first two options beat the third one in terms ofspeed, but we consider numerical accuracy here.
We evaluate ε for N = 1, ..., 8 (see Figure 9), and show the values for N = 8 in Table 2. The
numerical accuracy of the re-projected adjoint method matches that of differentiating “through”
re-orthogonalisation, and outperforms not re-projecting by a margin.

G

Memory-efficient kernel-matrix-vector products in JAX

Matrix-free linear algebra requires efficient matrix-vector products. For kernel function k = k(x, x′ ),and input data x1 , ..., xN , Gaussian process covariance matrices are of the form A = [k(xi , xj )]Ni,j=1 .
Matrix-vector products with A thus look likeNNX
v 7→ Av = k(xi , xj )vj 
(55)
j=1

i=1

and can be assembled row-wise, either sequentially or parallely.
The more rows we assemble in parallel, the faster the runtime but also the higher the memory
requirements, so we follow Gardner et al. [6] and choose the largest number of rows of A that still
fit into memory, say r such rows, and assemble Av in blocks of r. In practice, we implement thisin JAX by combining jax.lax.map and jax.vmap, but care has to be taken with reverse-mode
automatic differentiation through (v, θ) 7→ A(θ)v because by default, reverse-mode differentiation
stores all intermediate results. To solve this problem, we place checkpoints around each such batch ofrows, which reduces the memory requirements but roughly doubles the runtime. (We place anothercheckpoint around each stochastic trace-estimation sample, which roughly doubles the runtime again.)An alternative to doing this manually is the KeOps library [82], which GPyTorch [6] builds on.
However, there currently exists no JAX-compatible interface to KeOps which is why we have toimplement the above solution.
Figure 10 compares the runtime of our approach to that of KeOps custom CUDA code. We see that
we are competitive, but roughly 5× slower for medium to large datasets. Multiplying this with the
4× increase due to the checkpoints discussed above explains the 20× increase in runtime comparedto GPyTorch. Being 20× slower than GPyTorch per epoch is only due to the matrix-vector products,
and has nothing to do with the algorithm contribution. Future work should explore closing this gap
with a KeOps-to-JAX interface.

H

Experiment configurations for the Gaussian process study

Data For the experiments we use the “Protein”, “KEGG (undirected”, “KEGG (directed)”, “Elevators”, and “Kin40k” datasets (Table 7, adapted from Bartels et al. [95]). All are part of the UCI data
24

Table 7: Datasets used in this study.
DatasetSourceProteinElevators
Kin40K
KEGG (undir)
KEGG (dir)Available here.3
Camachol [96]Schwaighofer and Tresp [97]Shannon et al. [98]Shannon et al. [98]repository, and accessible through there.
The data is subsampled to admit the train/test split of 80/20%, and to admit an even division into
the number of row partitions. More specifically, we use 10 partitions for the kernel-matrix vector
products. This way, we have to discard less than 1% of the data; e.g., on KEGG (undir), we use
63,600 instead of the original 63,608 points.
We calibrate a Matérn prior with smoothness ν = 1.5, using 10 matrix-vector products per Lanczos
iteration, conjugate gradients tolerance of ϵ = 1, a rank-15 pivoted Cholesky preconditioner, and
10 Rademacher samples. We evaluate all samples sequentially (rematerialising on the backward
pass to save memory, as discussed in Appendix G). The conjugate-gradients tolerances are taken to
be absolute (instead of relative), and the parametrisations of the Gaussian process models and loss
functions match that of GPyTorch.
For every model, we calibrate an independent lengthscale for each input dimension, as well as anscalar observation noise, scalar output-scale, and the value of a constant prior mean. All parameters
are initialised randomly. We use the Adam optimiser with learning rate 0.05 for 75 epochs. All
experiments are repeated for three different seeds.

I

Partial differential equation data

We generate data for the differential equations as follows: Recall the problem setup of a partial
differential equation
 2

∂2
∂
∂2
2
u(t; x1 , x2 ) = ω(x1 , x2 )u(t; x1 , x2 ) +u(t; x1 , x2 )
(56)
∂2t
∂x21
∂x22
with Neumann boundary conditions. The coefficient field ω is space- but not time-dependent.
First, we discretise the Laplacian operator with central differences on an equidistant, tensor-productmesh that consists of 128 points per dimension, which yields 1282 grid points. The resulting
second-order ordinary differential equationd2
w = ω 2 M w,
dt2

(57)

where M is the discretised Laplacian, is then transformed into a first-order differential equation
  

d w
0
I
=
w =: Aw.
(58)
ω2 M 0dt ẇ
This equation is solved by the matrix exponential, and the system matrix A is asymmtric (byconstruction), and highly sparse because M is. Matrix-vector products with A are cheap, because we
can implement them with jax.scipy.signal.convolve2d.
Then, we sample a true ω from a Gaussian process with a square-exponential covariance kernel, using
lengthscale softplus(−0.75) and output-scale softplus(−10). We sample from this process with
the Lanczos algorithm [11] using Krylov-depth K = 32.
Then, we use another Gaussian process with the same kernel, but lengthscale softplus(0) andoutput scale softplus(0), to sample 256 initial distributions – again with the Lanczos algorithm
25Inputs

1.0
0.5
0.0

Targets

1.0

0.0

0.5

1.0

0.5
0.0

0.0

0.5

1.0

1.4 1.0
1.2
1.0
0.8
0.6 0.5
0.4
0.2
0.0
0.2 0.0
3.0 1.0
2.5
2.0
1.5
1.0 0.5
0.5
0.0
0.5
1.0 0.0

0.0

0.0

0.5

0.5

1.0

1.0

1.00 1.0
0.75
0.50
0.25
0.00 0.5
0.25
0.50
0.75
1.00 0.0
0.5
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5

1.0

0.0

0.5

1.0

0.5
0.0

0.0

0.5

1.0

1.2
0.9
0.6
0.3
0.0
0.3
0.6
0.9
1.2
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
7.5

Figure 11: Three exemplary input/output pairs from the PDE dataset.
[11]. These 256 initial conditions are solved with Diffrax’s implementation of Dopri8 [99] using 128
timesteps. Some example input/output pairs are in Figure 11.
This setup is similar to that of the WaveBench dataset [92], with the main difference being that the
WaveBench dataset uses a slightly different formulation of the wave equation.4 We use the one above
because it lends itself more naturally to matrix exponentials, which are at the heart of this experiment.

J

Implementation details for the Bayesian neural network study

J.1

Bayesian neural networks with Laplace approximations

Another possible application of the gradients of matrix functions is marginal-likelihood-based
optimisation of Bayesian Neural Networks. Suppose gθ (x) is the output of a neural network withparameters θ ∈ RP . The choice of the model shall be denoted by M and consist of both continuousand discrete hyperparameters (such as network architecture, likelihood precision, prior precision,etc.). For some choice of prior given byp(θ | M)(59)p(y|x, θ, , M) = p(y | x, gθ (x), M)(60)

and likelihood
we can specify a Bayesian model. The posterior distribution is then given by:p(θ, y | x, M) ∝ p(y | x, gθ (x), M)p(θ | M)dθ.

(61)

The marginal likelihood is given by normalizing constant of this posterior, i.e.
Zp(y | x, M) = p(y | x, gθ (x), M)p(θ | M)dθ.

(62)

As suggested by MacKay [101], this marginal likelihood can be used for model selection in Bayesianneural networks. Immer et al. [30] use the Laplace approximation of the posterior to obtain access tothe marginal likelihood of the Bayesian neural network and its stochastic gradients.
The Laplace approximation of the marginal likelihood is given by:1
log p(y | x, M) ≈ log p(y, θMAP | x, M) − log det
2



1
HθMAP
2π


(63)

where HθMAP = −∇2θ log p(y, θMAP | x, M). Usual choices of the prior are N (0, α−1 I). Usually
this Hessian is approximated with the generalized Gauss-Newton (GGN) matrix [102]HθMAP ≈ A(α) :=JX[Dθ gθMAP ])(xj )⊤ [Dg2 ρ](yj , gθMAP (xj ))[Dθ gθMAP ])(xj )⊤ + α2 I

(64)

j=1
3

Link: http://archive.ics.uci.edu/dataset/265/physicochemical+properties+of+protein+tertiary+structure

4

To ensure radiating boundary conditions, Liu et al. [92] follow Stanziola et al. [100]’s model of simulating
the wave equations as a sytem of first-order equations.

26

where D2 ρ is the Hessian of the loss, and Dθ g the Jacobian of g (recall Equation 2). This objective
is used to optimize the prior precision of the model or any continuous model hyperparameters.
Matrix-vector products with the GGN matrix can be accessed through automatic differentiation using
Jacobian-vector and vector-Jacobian products. With these efficient matrix-vector products, one canestimate the log-determinant of GGN using matrix-free techniques like the Lanczos iteration.
To make predictions using the Laplace approximation of the posterior, we also need to sample fromthe normal distribution N (θMAP , A−1 ). Samples from this distribution can be written as:θ = θMAP + A−1/2 ϵ(65)

where ϵ ∼ N (0, I). The main bottleneck in this computation is the inversion and matrix square rootof the GGN matrix, and we implement it with a Lanczos iteration using f (x) = x−1/2 . Since the
GGN is empirically known to have low-rank [103], doing a few Lanczos iterations can get us close toan accurate estimation.
J.2

Experiment setup

We estimate the diagonal of the GGN stochastically via (“◦” is the element-wise product) [104]L

diagonal (A) = E[v ◦ Av] ≈1X
vℓ ◦ Av,
L

E[vv ⊤ ] = I.

(66)

ℓ=1

We use 150 matrix-vector products for both diagonal calibration and our Lanczos-based estimation.
We use 30 Monte-Carlo samples to estimate the log-likelihoods for evaluating the test metrics, and weuse places365 [91] as an out-of-distribution dataset to compute OOD-AUROC. We also computethe expected calibration error (ECE) [105] of the model.
Data: We show scalability by doing Laplace approximation on Imagenet1k image classification [89].
The training set consists of approximately 1.2 million images, each belonging to one of 1000 classes.
We find that we can take small subsets of this dataset and still converge to the same prior precision.
Our computational budget allows us to use 10 percent of the samples for each class. However, evenfor very small subsamples of the data, we converge to a very similar prior precision.
Method: To optimize the prior precision we use the marginal likelihood as the objective. We use
the RMSprop optimizer with a learning rate of 0.01 for 100 epochs for optimizing both the diagonal
and Lanczos approximations of the GGN.

27

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the
paper’s contributions and scope?
Answer: [Yes]
Justification: The main contributions, adjoint systems for the Lanczos and Arnoldi iterations,are explained in Section 4. The three case studies are in Sections 5 to 7.
Guidelines:• The answer NA means that the abstract and introduction do not include the claims
made in the paper.
• The abstract and/or introduction should clearly state the claims made, including thecontributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
• The claims made should match theoretical and experimental results, and reflect how
much the results can be expected to generalize to other settings.
• It is fine to include aspirational goals as motivation as long as it is clear that these goals
are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The limitations are discussed in the paragraph titled “Limitations and futurework” on page 4. All assumptions (i.e., Assumptions 3.1 and 3.2) are contextualised in thesentences before and after they are have been introduced.
Guidelines:• The answer NA means that the paper has no limitation while the answer No means thatthe paper has limitations, but those are not discussed in the paper.
• The authors are encouraged to create a separate ”Limitations” section in their paper.
• The paper should point out any strong assumptions and how robust the results are toviolations of these assumptions (e.g., independence assumptions, noiseless settings,
model well-specification, asymptotic approximations only holding locally). The authors
should reflect on how these assumptions might be violated in practice and what the
implications would be.
• The authors should reflect on the scope of the claims made, e.g., if the approach was
only tested on a few datasets or with a few runs. In general, empirical results often
depend on implicit assumptions, which should be articulated.
• The authors should reflect on the factors that influence the performance of the approach.
For example, a facial recognition algorithm may perform poorly when image resolution
is low or images are taken in low lighting. Or a speech-to-text system might not be
used reliably to provide closed captions for online lectures because it fails to handle
technical jargon.
• The authors should discuss the computational efficiency of the proposed algorithms
and how they scale with dataset size.
• If applicable, the authors should discuss possible limitations of their approach toaddress problems of privacy and fairness.
• While the authors might fear that complete honesty about limitations might be used byreviewers as grounds for rejection, a worse outcome might be that reviewers discover
limitations that aren’t acknowledged in the paper. The authors should use their best
judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers
will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions anda complete (and correct) proof?
28

Answer: [Yes]
Justification: The main contributions, Theorems 4.1 and 4.2 and Corollary 4.3, are proven in
Appendices B to D. Appendix E discusses solving the adjoint system in full detail.
Guidelines:• The answer NA means that the paper does not include theoretical results.
• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced.
• All assumptions should be clearly stated or referenced in the statement of any theorems.
• The proofs can either appear in the main paper or the supplemental material, but ifthey appear in the supplemental material, the authors are encouraged to provide a shortproof sketch to provide intuition.
• Inversely, any informal proof provided in the core of the paper should be complemented
by formal proofs provided in appendix or supplemental material.
• Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusionsof the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: The most important information about the experiment setup is a part of themain paper in Sections 5 to 7; further information can be found in Appendices G to J.
Appendices A and F discuss the setup for Figure 3 and Table 2. Code will be publishedupon acceptance.
Guidelines:• The answer NA means that the paper does not include experiments.
• If the paper includes experiments, a No answer to this question will not be perceived
well by the reviewers: Making the paper reproducible is important, regardless ofwhether the code and data are provided or not.
• If the contribution is a dataset and/or model, the authors should describe the steps taken
to make their results reproducible or verifiable.
• Depending on the contribution, reproducibility can be accomplished in various ways.
For example, if the contribution is a novel architecture, describing the architecture fully
might suffice, or if the contribution is a specific model and empirical evaluation, it may
be necessary to either make it possible for others to replicate the model with the samedataset, or provide access to the model. In general. releasing code and data is often
one good way to accomplish this, but reproducibility can also be provided via detailed
instructions for how to replicate the results, access to a hosted model (e.g., in the caseof a large language model), releasing of a model checkpoint, or other means that areappropriate to the research performed.
• While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on thenature of the contribution. For example(a) If the contribution is primarily a new algorithm, the paper should make it clear how
to reproduce that algorithm.
(b) If the contribution is primarily a new model architecture, the paper should describe
the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should
either be a way to access this model for reproducing the results or a way to reproduce
the model (e.g., with an open-source dataset or instructions for how to construct
the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which caseauthors are welcome to describe the particular way they provide for reproducibility.
In the case of closed-source models, it may be that access to the model is limited insome way (e.g., to registered users), but it should be possible for other researchers
to have some path to reproducing or verifying the results.
29

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental
material?
Answer: [Yes]
Justification: Code has been submitted as a part of the supplementary material, and will be
published upon acceptance.
Guidelines:• The answer NA means that paper does not include experiments requiring code.
• Please see the NeurIPS code and data submission guidelines (https://nips.cc/
public/guides/CodeSubmissionPolicy) for more details.
• While we encourage the release of code and data, we understand that this might not be
possible, so “No” is an acceptable answer. Papers cannot be rejected simply for notincluding code, unless this is central to the contribution (e.g., for a new open-source
benchmark).
• The instructions should contain the exact command and environment needed to run toreproduce the results. See the NeurIPS code and data submission guidelines (https://nips.cc/public/guides/CodeSubmissionPolicy) for more details.
• The authors should provide instructions on data access and preparation, including how
to access the raw data, preprocessed data, intermediate data, and generated data, etc.
• The authors should provide scripts to reproduce all experimental results for the new
proposed method and baselines. If only a subset of experiments are reproducible, theyshould state which ones are omitted from the script and why.
• At submission time, to preserve anonymity, the authors should release anonymized
versions (if applicable).
• Providing as much information as possible in supplemental material (appended to thepaper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand theresults?
Answer: [Yes]
Justification: See the answer to “4. Experimental Result Reproducibility” above.
Guidelines:• The answer NA means that the paper does not include experiments.
• The experimental setting should be presented in the core of the paper to a level of detailthat is necessary to appreciate the results and make sense of them.
• The full details can be provided either with the code, in appendix, or as supplemental
material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate
information about the statistical significance of the experiments?
Answer: [Yes]
Justification: All case studies report mean and standard deviations of multiple runs. The
only exception is the Bayesian neural network example, which uses a single training run(but evaluates test metrics on multiple seeds).
Guidelines:• The answer NA means that the paper does not include experiments.
• The authors should answer ”Yes” if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support
the main claims of the paper.
30• The factors of variability that the error bars are capturing should be clearly stated (forexample, train/test split, initialization, random drawing of some parameter, or overall
run with given experimental conditions).
• The method for calculating the error bars should be explained (closed form formula,call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).
• It should be clear whether the error bar is the standard deviation or the standard errorof the mean.
• It is OK to report 1-sigma error bars, but one should state it. The authors should
preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesisof Normality of errors is not verified.
• For asymmetric distributions, the authors should be careful not to show in tables orfigures symmetric error bars that would yield results that are out of range (e.g. negative
error rates).
• If error bars are reported in tables or plots, The authors should explain in the text how
they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce
the experiments?
Answer: [Yes]
Justification: The introductory part of the appendix contains a paragraph titled “Compute”.
Guidelines:• The answer NA means that the paper does not include experiments.
• The paper should indicate the type of compute workers CPU or GPU, internal cluster,or cloud provider, including relevant memory and storage.
• The paper should provide the amount of compute required for each of the individual
experimental runs as well as estimate the total compute.
• The paper should disclose whether the full research project required more computethan the experiments reported in the paper (e.g., preliminary or failed experiments that
didn’t make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the
NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?
Answer: [Yes]
Justification: We use datasets that are either self-created (Section 6 and Appendix I), or
common test-cases for machine learning methods (UCI datasets, ImageNet) or numericalalgorithms (SuiteSparse matrix collection).
Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
• If the authors answer No, they should explain the special circumstances that require a
deviation from the Code of Ethics.
• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative
societal impacts of the work performed?
Answer: [No]Justification: This research provides a foundational algorithm for computational sciences,and societal impact is difficult if not impossible to predict.
Guidelines:31• The answer NA means that there is no societal impact of the work performed.
• If the authors answer NA or No, they should explain why their work has no societal
impact or why the paper does not address societal impact.
• Examples of negative societal impacts include potential malicious or unintended uses
(e.g., disinformation, generating fake profiles, surveillance), fairness considerations(e.g., deployment of technologies that could make decisions that unfairly impact specificgroups), privacy considerations, and security considerations.
• The conference expects that many papers will be foundational research and not tied
to particular applications, let alone deployments. However, if there is a direct path to
any negative applications, the authors should point it out. For example, it is legitimate
to point out that an improvement in the quality of generative models could be used togenerate deepfakes for disinformation. On the other hand, it is not needed to point out
that a generic algorithm for optimizing neural networks could enable people to trainmodels that generate Deepfakes faster.
• The authors should consider possible harms that could arise when the technology is
being used as intended and functioning correctly, harms that could arise when thetechnology is being used as intended but gives incorrect results, and harms following
from (intentional or unintentional) misuse of the technology.
• If there are negative societal impacts, the authors could also discuss possible mitigation
strategies (e.g., gated release of models, providing defenses in addition to attacks,
mechanisms for monitoring misuse, mechanisms to monitor how a system learns fromfeedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible
release of data or models that have a high risk for misuse (e.g., pretrained language models,
image generators, or scraped datasets)?
Answer: [NA]Justification: This paper does not contribute data or models that would require safeguards.
Guidelines:• The answer NA means that the paper poses no such risks.
• Released models that have a high risk for misuse or dual-use should be released withnecessary safeguards to allow for controlled use of the model, for example by requiringthat users adhere to usage guidelines or restrictions to access the model or implementing
safety filters.
• Datasets that have been scraped from the Internet could pose safety risks. The authors
should describe how they avoided releasing unsafe images.
• We recognize that providing effective safeguards is challenging, and many papers do
not require this, but we encourage authors to take this into account and make a best
faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in
the paper, properly credited and are the license and terms of use explicitly mentioned andproperly respected?
Answer: [Yes]
Justification: All assets used in the experiments have been cited. See also the answer to “9.
Code of Ethics”.
Guidelines:• The answer NA means that the paper does not use existing assets.
• The authors should cite the original paper that produced the code package or dataset.
• The authors should state which version of the asset is used and, if possible, include a
URL.
• The name of the license (e.g., CC-BY 4.0) should be included for each asset.
32• For scraped data from a particular source (e.g., website), the copyright and terms ofservice of that source should be provided.
• If assets are released, the license, copyright information, and terms of use in the
package should be provided. For popular datasets, paperswithcode.com/datasets
has curated licenses for some datasets. Their licensing guide can help determine thelicense of a dataset.
• For existing datasets that are re-packaged, both the original license and the license of
the derived asset (if it has changed) should be provided.
• If this information is not available online, the authors are encouraged to reach out tothe asset’s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation
provided alongside the assets?
Answer: [Yes]
Justification: The code (attached to this submission) is documented. The appendices containall other information.
Guidelines:• The answer NA means that the paper does not release new assets.
• Researchers should communicate the details of the dataset/code/model as part of theirsubmissions via structured templates. This includes details about training, license,
limitations, etc.
• The paper should discuss whether and how consent was obtained from people whose
asset is used.
• At submission time, remember to anonymize your assets (if applicable). You can either
create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper
include the full text of instructions given to participants and screenshots, if applicable, as
well as details about compensation (if any)?
Answer: [NA]Justification: This study does not involve crowdsourcing nor human subjects.
Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research withhuman subjects.
• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be
included in the main paper.
• According to the NeurIPS Code of Ethics, workers involved in data collection, curation,or other labor should be paid at least the minimum wage in the country of the datacollector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human
Subjects
Question: Does the paper describe potential risks incurred by study participants, whether
such risks were disclosed to the subjects, and whether Institutional Review Board (IRB)
approvals (or an equivalent approval/review based on the requirements of your country or
institution) were obtained?
Answer: [NA]Justification: This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research withhuman subjects.
33• Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you
should clearly state this in the paper.
• We recognize that the procedures for this may vary significantly between institutions
and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and theguidelines for their institution.
• For initial submissions, do not include any information that would break anonymity (ifapplicable), such as the institution conducting the review.

34