Gradients of Functions of Large Matrices
Authors: Nicholas Krämer
Pablo Moreno-Muñoz
Hrittik Roy
Søren Hauberg
Venue: NeurIPS 2024 spotlight (2024)
Keywords: Automatic differentiation
numerical methods
linear algebra
implicit differentiation
adjoint methods
differential equations
Bayesian neural networks
Gaussian processes
Primary area: machine_learning_for_physical_sciences
TL;DR: We derive previously unknown gradients of Lanczos and Arnoldi iterations and use them for PDEs, Gaussian processes, and Bayesian neural networks.

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_ these quantities 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, and show 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 [link redacted] and install the library with *pip install [redacted]*.

Gradients of Functions of Large Matrices
Nicholas Krämer, Pablo Moreno-Muñoz, Hrittik Roy, Søren Hauberg
Technical University of Denmark Kongens Lyngby, Denmark
(pekra, pabmo, hroy, sohau sohau)@dtu..o @dtu. dk
Abstract
Tuning scientific and probabilistic machine learning models for example, par- tial 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 these quantities almost always based on Lanczos and Arnoldi iterations, the present work is the first to explain how to differentiate these workhorses of numeri- cal linear algebra efficiently. To get there, we derive previously unknown ad- joint systems for Lanczos and Amnoli iterations, implement them in JAX, and show that the resulting code can compete with Diffrax when it comes to dif- ferentiating PDEs, GPy'Torch for selecting Gaussian process models and beats standard factorisation methods for calibrating Bayesian neural networks. All this is achieved without any problem-specife code optimisation. Find the code at https: //gtthb..... install the library with pip install matfree.
httpp/aeeee/eppereeeeeeeee waannnosssaaansss-eeee lan-zos-aaaaaaaaaa- -adjointaaatt and
Introduction
Automatic diffcrentiation has dramatically altered the development of machine learning models by allowing to forego laborious, application-dependenl gradient derivations. The essence of this automation 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 neural network. 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 new matrix-free method for automatically differentiating functions of matrices. Our algorithm yields the exact gradients of the forward pass, all gradients are obtained with the same code, and said code runs in linear time- and memory comoorycomplexity.
For a parametrised matrix A = A(0) E NXN and an analytic function f TR. - R. we call ((A function of the matrix (different properties of A imply different definitions of f(A); diagonalisable; see [1]). However, we assume that A is the
of them is
applying ff to each eigenvalue of Jacobian of large neural network or matrix of similar size and never materialise f(A). Instead, product
only care about the values and gradients of the matrix-function-ved
(8, v) + f[A(O)]u
(1)
assuming that A is only accessed via differentiable matrix-vector products. Table lists examples.
Evaluating Equation is crucial for building large machine learning models, e.g., Bayesian neural
networks: A common hyperparameter-caliireee
loss of (Laplace-approxiiatel) Bayesian neural network involves the log-determinant of the generalised Gauss-Newton matrix [13]
A(a) =
[Dog(<x)) [t)) +22,
(2)
Edata
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Table 1: Some applications for functions of matrices. Log-determinants apply by combining log det(A)= = trace (Ao) with stochastic trace estimation, which is why most vectors in this table are Rademacher samples. "PDE" "ODE"" = "Partial/Ordinary differential equation".
Application
Function f
Matrix A
Vector v
Parameter 0
PDEs & flows [2-5]
e
PDE discret.
PDE initial value
PDE
Gaussian process [6-8]
log(A) log(1 +
Kernel matrix
V Rademacher
Kernel
Invert. ResNets [9, 10]
Jacobian matrix
Rademacher
Network
Gaussian sampler [11]
Covariance matrix
V ~~ N N(0, 7)
Covariance
Neural ODE [12]
X²
Jacobian matrix
V Rademacher
Network
where Dog is the parameter-Jacobian of the neural network g. D2p is the Hessian of the loss function with respect to ( and a is to-be-tuned parameter. The matrix A(a) in Equation 2 has as many rows and columns as the network has parameters, which makes traditional, cubie-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 is more structured than it actually is, e.g., diagonal [14]; (ii) approximating log det(A) by combining stochastic trace estimation [15]
trace (A) E Aul ~
oa Ave, for
E vv
I,
(3)
L
8=1
with Lanczos iteration A(0) QHQ T
[16), reduce the log-determinant to [17, 18]
1
1 Σ de oi log((()) L
log det(AA) = trace (log A) ~
Σ Uy log(A)ve ~
ve.
(4)
L
The matrix H in A =0 QHQT T has as many rows/columns as willing to evaluate matrix vector products with A; thus, it is small enough to evaluate the matrix lagrix-logarithm log(H) in cubic complexity.
Contributions
This article explains how to differentiate not just log-determinants but any Lanezos and Arnoldi iteration so we can build loss functions for large models with such matrix-free algorithms (thereby completing the pipeline in Figure 1). This kind of functionality has been sorely missing from the toolbox of differentiable programming until now, even though the demand for functions of matrices is high in all of probabilistic and scientific machine learning [e.g. 2-12, 19-31].
Parameter
Optimise
Large matrix-vector product
V
Arnoldi, Lanczox
Contribution
Small factorisation
"Chcap" linalg
(known)
2 Related work
Function large matrix
V
Here, we focus on applications in machine learning and illustrate how prior work avoids differentiating matrix- free decomposition methods like the Lanczos and Arnoldi iterations. Golub and Meurant [32] discuss applications outside machine learning.
Loss
(known)
Figure 1: Values (down) and gradients (up) of functions of large matrices.
Generative models, e.g., normalising flows [27, 33], rely the change-of-variable. formula, which involves the log-determinant of the Jacobian matrix of a neural network. Behrmann et al. [9] and Chen et al. [10] combine stochastic trace estimation with a Taylor-series expansion for the matrix logarithm. Ramesh and LeCun [34] use Chebyshev expansions instead Taylor expansions. That said, Ubaru et al. [17] demonstrate how both methods converge more slowly than the Lanczos iteration when combined with stochastic trace estimation.
Gaussian process model selection requires values and gradients of log- pog-robbability density functions of 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 some

---

combination of stochastic trace estimation with the Lanczos iteration, and unanimously identifies ("d" shall be an infinitesimal perturbation; see Section 4)
gradients of log-determinant
I == log det det(K(0)),
du trace ((())))a((0)))
(5)
Another round of stochastic trace estimation then estimates du [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 only works for log-determinants. Section 5 shows how our black-box gradients match state-of-the-ar code for Equation 5 [6].
Laplace approximations and neural tangent kernels face the same problem of computing derivatives of log-determinants but with the generalised Gauss-Newton (GGN) matrix from Equation 2. In contrast 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 approx- imations 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 automati- cally. Unfortunately, these approximations discard valuable information about the correlation between weights, so tmear-algebrrabasee approach leads to superior likelihood calibration (Section 7).
Linear differential equations, for instance (t) Ayt)), y(0) 3/o are solved by matrix exponentials, v(t = exp(At At) By this relation, matrix exponentials have frequent applications not just for the simulation 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, 47J. 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 rich history of adjoint methods [e.g. 51, 52] with high-performance open-soure source libraries [53-56]. Still, the (now differentiable) Arnoldi iteration can compete with state-of-the-an solvers in JAX (Section 6).
3 Problem statement
Recall A = A(0) RNXN from Section 1. The focus of this paper is matrices that too large to store in memory, like Jacobians of neural networks or discretised partial differential equations: Assumption 3.1. A(0) only accessed via differentiable matrix-vector products (0, v) A(0)..
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 sys- tems [58]. But there is more to matrix-free linear algebra than conjugate gradient solvers: Matrix-free implementations of matrix decompositions usually revolve around variations of the Arnoldi itera- tion [59], which takes an initial vector E R. and prescrihed number of iterations K E N and produces a column-orthogonal Q E RNNKK structured H R residual vector r E TR,. and length e E IR such that
x
A(0)
=
Q
AQ QH r(ek) and Qe
Figure 2: Lanczos/Armoldi iteration.
hold (Figure 2; C E 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 popular for imple- menting
products in a matrix-free fashion [1, 57), because the decomposition in Equation 6 implies A QHQ thus
(0, ~ (((())) =
The last step. Q tee. 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 H1 + f(H)e is differentiable, and Q fits into memory.
In summary, we evaluate functions of large matrices by firstly decomposing large matrix into product of small matrices (with Lanczos Amoldi) and, secondly, using conventional linear algebra to evaluate functions of small matrices. Functions of small matrices can al- ready be differentiated efficiently [60-62]. This work contributes gradients of the Lanczos and Arnoldi iteration under Assumptions 3.1 and 3.2, and thereby makes matrix-free implementa- tions of matrix decompositions and functions of large matrices (reverse-mode differentiable.
3

---

Automatic differentiation, i.e. "backpropagating through" the matrix decomposition, is far too inefficient to be vi- able option (Figure 3; setup in Appendix A). Our approach via implicit differentiation or the adjoint method, respec- tively, leads to gradients that inherit the linear runtime and mempry-complexity of the forward pass.
0.6
SuiteSarrse 'besstkII'
0.5
Wall time (sec)
0.3
N- 11948 rows
Adjaint
0.2
Backprop
Limitations and future work The landscape of Lanczos- and Arooddi style matrix decompositions is vast, and some adjacent problems cannot be solved by this single article: (i) Forward-mood derivatives would require derivation separate from what comes next. Yet, since functions of matrices map many to few parameters (matrices to vec- tors), reverse-mode is superior to forward-mode anyway [66, p. 153]. (ii) We only consider real-valued matrices (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 to combining Arnoldi/Lanczo with full reorthogonalisation [57, 68-71J. Relaxing this assumption requires gradients of partial reorthogonalisationn (among other things), which we leave to future work.
Forward pass
0.0
200
250
Krylov-space depth
Figure 3: Backpropagation our ad- joint method on sparse matrix [63-65].
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,
differentiate the iterations like dynamical system using the "adjoint method" [51. 73, 74], a variation of implicit differentiation that uses Lagrange multipliers [66], and not like 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 dæ be an infinitesimal perturbation of some æ. D is the Jacobian operator, and (:--) the Euclidean inner product between two equally-sized inputs. For a loss p E R that depends on some x, the linearisation dp- Dzp dr and the gradient identity dp <VzP, dx) will be important [76, 77].
Implicit differentiation Let a + x numerical algorithm that computes some æ from some 0. Assume that the input and output of a(.) satisfy the constraint c(0, a(0)) 0. For instance, if a(-) solves Ax b, the constraint cf b; Az 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 0= 0=c((0,x) x) = c,0 x) + Doc(O, x)00
(8)
In other words. we "lincarise" the constraint c(0,z) = 0. The adjoint method [78] proceeds by "transposing" this lincarisation follows. Let p be a loss that depends on y with gradient v and recall the gradient identity from the "Notation" paragraph above. Then, for all Lagrange multipliers A with the same shape as the outputs of c(-), know that since c 0 implies dc 0.
dp = (V.,, = (V.P, dæ) (c) = (V.p+ (Dac) A, + ((Doc) T X, d0)
must hold. By matching Equation 9 to dp (Vop, dl)) (this time, regarding p as function of 0, of x;; recall the "Notation" paragraph), we conclude that if X solves the adjoint system
Vzp+ (Dac)
0,
(10)
then Vep:= (Dac) must be the gradient of p with respect input 0. This is the adjoint method [66, Section 10.4]. In automatic differentiation frameworks like JAX [79], this gradient implements 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-vas valued constraints, such as dynamical systems or the Lanezos and Arnoldi iterations. The reason is that the adjoint method does not change c(-) becomes matrix- or function-spact spaeevi--spaee-llleddd as long as we can define inner products and adjoint operators, whereas other approaches (like what Blondel et [72] use for numerical optimisers) would become increasingly laborious in these cases. In summary, to reverse- mode differentiate

---

numerical algorithm with the adjoint method, we need four steps: (i) find constraint, (ii) linearise it, (iii) introduce Lagrange multipliers, and (iv) solve the resulting adjoint system. Carrying out those four 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.1
Adjoint system of the Arnoldi and Lanczos iterations
Let ej be the jth unit vector. Denote by "o" the wise matrix product, and define the matrices
K = [8igslig=1 L< == [6t2]lkz=1> L [64+1<5
(11)
so that for example, 1g o A extracts the lower triangular matrix of (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 E N, R. and A RNXN, and loss p()) E IR be given. IfQ E RNNK H E R T. RN, and c E R solve the forward constraint
AQ QH + t((e)), Q61 vC, I [ II
I o H =0 0, Q r 0,
(12)
and if X E RN, A ARRN YE RK I' RRaK and E E R K k satisfy the adjoint system
0 VQp+ A AH + +(e)) + Q(1 .T) Q rT) + +m7
(13a)
0 = VHP- VHP-Q"A- A l oF
(13b)
0 = V ÀCK + Q7
(13c)
0= VeP vT>,
(13d)
then the gradients of p with respect to A and are
AP = AQT, V = Ac
(14)
Sketch of the proof. To derive the statement, start with Equation 12 as c(-). Apply the chain- and product rules liberally to get dc(.). Introduce Lagrange multipliers X. A, I, and E like in the previous section, by adding Lagrange-multiplec constraints to
dp- = (Vop, dQ) + (VHP, d + (Vrp, dr) + (ep,ddc)
(15)
and rearrange the terms to see that Equation 13 implies Equation 14. Details are Appendix B. 0
Theorem 4.2 (Adjoint system of the Lanezos iteration). Let symmetric A N well
as E RN, K E N, and a loss be known. In the following equations, set bo == 1 E IR To 0 E R" AK+1 = pto 0, Mo 0 to simplify the expressions. y1, E R. and (1. ..., Ok. b1 ..... by E RRK satisfy the forward constraint
I1 v/(uT
0,
(16a)
-h-111 18 + (A - 0% I) 13 1x = 0,
k = 1,... K
(16b)
F 1 Tk+1 1 = 0,
k = 1, ,..., K
(16c)
k = 2, .... K
(16d)
and if Au, dK E R.N, Mi.... MK, E IR satisfy the adjoint system
0- = AxbK + (VZK-1P) P PK + vk
(17a)
0- buduui (AT 0n a1())= by 1--1 + VorP P- + 148- 13 PE + UR 10k-1)
(17b)
0 = VagP AI IT:
(17c)
0 Vo.P- XT1 XT 2k+1-
(17d)
where all expressions involving hold for all ls K, ..., 1, then
k
AJ Φ1 rı
Ao
V = Σ
A
(18)
R v
k=
the gradients of p with respect and A.
5

---

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 symmetric matrices, H must be tridiagonal [57], and we expand AQ = QH + r((x+1() column-wise. The coefficients ag and by are the tridiagonal elements in H. We rename ds from Arnoldi to 2% for
Lanczos to make it easier to distinguish the two different sets of constraints. Details: Appendix C. L
4.2 Matrix-free implementation
Solving the adjoint systems To compute Ap and Vep. we need solve the adjoint systems. When comparing the forward constraints to the adjoint systems. similarities emerge: for instance, the adjoint system of the Amoldi iteration follows the same A H - AH(T) + rest = 0 structure as the forward constraint. This structure suggests deriving recursion for the backward pass that mirrors that of the forward pass, Appendix E contains this derivation and contrasts the resulting algorithm with 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 with matrix-vector products, which is helpful because it makes our custom backward pass as matrix- free as backpropagation "through" the forward pass would be. This matrix muatrix-free free implementation in combination with the efficient recursions Theorems 4.1 and 4.2 explains the significant performance gains of our method compared to naive backpropagation, observed in Figure 3.
Theorems 4.1 and 4.2's expressions for V Ap are not directly applicable when we only have matrix- vector products with A. Fortunately, aurameter-gradient 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 0, the gradients of p with respect to 0 are
K (ex) ATO) AcE].
Vep
(19)
which can be assembled online during the backward pass. For the Lanczos iteration, assume the conditions of Theorem 4.2 instead of Theorem 4.1, replace Qek and Aeş with Ty and NR, let the sum run from k- 1 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 V Ap from Theorems 4. and 4.2 with dA Dg AdO. The derivations are lengthy and therefore relegated to Appendix D. I
Reorthogonallsation It is well known that the Lanczos and Arnoldi iterations suffer from a loss of orthogonality and that reorthogonalisation of the columns in Q is often necessary [68-71]. Reorthogonalisation does not affect the forward constraints, so the adjoint systems remain the same with and without reorthogonaliiation.. But adjoint systems also suffer from a loss of orthogonality: The equivalent of orthogonality for the adjoint system is the projection constraint in Equation 13b, which constrains the Lagrange multipliers A a hyperplane defined by Q and other known quantities. The constraint can and should (Table 2)- - be used whenever the forward pass requires reorthogonatisation. ! In the case studies below, we always use full reorthogonalisation on the forward and adjoint pass, also for the Arnoldi iteration [71, Table 7.1]. even though this is slightly less common than for the Lanezos iteration.
Table 2: Accuracy loss when differenti- ating the Arnoldi iteration on Hilbert matrix in double precision (ó decom- pose with a full-rank Arnoldi iteration, then reconstruct the original matrix; mea- sure |
IlI: details in Appendix F).
Loss of accuracy
Adjoint w/o proj. Adjoint w/ proj.
5.83 10-3
1.17 10 -10
Backprop.
1.17 10
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-gee -rradeeeg-grdettss... in linear time- and space- complexity (Corollary 4.3). While some of the derivations are somewhat technical, the overall
'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 Arnolii code.

---

Table 3: Our method yields the same root-mean- root-mean-square errors (RMSEs) as GPyTorch. reaches
lower training losses but is = x0x slower per epoch due to different matrix-vector-prodde backends
(see Appendix G). Three runs, significant improvements in bold. We use an 80/20 train/test split.
Dataset
Size
Dim.
Method
RMSE +
Final training loss +
Runtime (s/epoch)
Adjoints
0.09 + 0.002
-0.91 0.025
1.69 + 0.000
elevators
16,599
18
GPyTorch
0.09 = 0.003
-0.63 = 0.062
0.10 = 0.004
Adjoints
0.39 + 0.005
0.73 0.300
12.62 + 0.172
protein
45,730
GPyTorch
0.39 = 0.005
0.73 + 0.075
0.73 + 0.039
Adjoints
0.12 + 0.004
-0.30 + 0.078
8.27 + 0.004
kin40k
40,000
GPyTorch
0.10 + 0.010
-0.26 + 0.094
0.26 + 0.024
Adjoints
0.12 = 0.002
-0.59 + 0.295
13.25 + 0.005
kegg-dir
48,827
20
GPyTorch
0.12 + 0.005
-0.41 + 0.054
0.62 + 0.262
Adjoints
0.12 = 0.002
-0.69 + 0.263
24.20 + 0.004
kegg undir 63,608
26
GPyTorch
0.12 = 0.003
-0.40 + 0.039
1.67 + 0.532
approach follows the general template for the adjoint method relatively closely. The resulting algorithm beats backpropagation "through" the iterations by 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 be installed via "pip install matfree". Next, we put this code to the test on three challenging machine-learning problems centred around functions of matrices see how it fares against state-of- the-art differentiable implementations of exact Gaussian processes (Section 5), differential equation solvers (Section 6), and Bayesian neural networks (Section 7).
Case study: Exact Gaussian processes
Model selection for Gaussian processes has arguably been the strongest proponent the Lanczos
iteration and similar matrix-fres algorithms in recent years [6, 7, 11, 19, 20, 23, 80], and most of these efforts have been bundled up the GPyTorch library [6]. For example, GPyTorch defaults to choosing Lanczos iteration over Cholesky decomposition data points. ? Calibrating hyperparameters of Gaussian process models involves optimising log- marginal-likelhhoode likelihoods of the regression targets, which requires computing AA-'a 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 du trace dA) (Equation 5). Since we seem to be the first to take different path, benchmarking Gaussian processes comparison to GPyTorch is a good first testbed for our gradients.
the dataset exceeds 800
Setup: Like GPyTorch's defaults We mimic recent suggestions for scalable Gaussian process models [7. 19]: we implement a pivoted Cholesky preconditioner [81] and combine it with con- jugate gradient solvers for x: K-1a (which can be differentiated efficiently). We estimate the log-determinant stochastically via log (At( = trace (log(A)) E log(A)o), and compute log(A) (Aju via the Lanezos iteration. While all of the above is common for "exact" Gaussian processes [6, 7, 19] ("exact" as opposed variational approaches, which are not relevant for this comparison), there are three key differences between our code and GPyTorch's s: (i) GPyTorch is in Pytorch and uses KeOps [82] for efficient bernel-matrix-ve low-memory matrix-vecto products (Appendix G). (ii) GPyTorch runs all algorithms adaptively (we specify tolerances and maximum iterations as much as possible). We use adaptive conjugate gradient solvers and fixed ranks for everything else. (iti) GPytorch differentiates the log determiitnnt with a tailored approximation of Equation 5 [6]: we embed our gradients of the Lanczos iteration into automatic differentiation. To keep the benchmark objective, mimic the parameter suggestions from GPy Torch's 's default settings, and optimise hyperparameters of a Matérn ()) model on UCI datasets with the Adam optimiser [83]. Appendix H lists parameters and discusses the datasets.
products. We use JAX and must build our own
Parameter mww.cholesky.ctt.mo http::::iiiitttt///ts
attps:///oce wrw.occs..gttrrrc...cooooo wpw.yorc..../e/////ttaee//tttttt
7

---

A1.
Forward error
A2.
Gradient error
8.
Training loss
(rers) SEEE
RMSE (relative)
10"
Dears
10a
Depet
Deard, opochh
20-2
10-2
20.1
10
10-2
102
10a
10
10a
0
100
200
300
400
500
600
700
800
No. Matvecs
No. Matvecs
Time (sec)
Figure 5: Aroldi's superior convergence on the forward pass (A1) is inherited by the gradients (A2; mind the shared y-axis) and ultimately leads to fast training (B). For training, Arnoldi uses ten matrix-vector products, and the other two use 15 (so they have equal error R 10-1 in Al and A2.)
Analysis: Trains like GPy 3PyTorch; large scale only limited by matrix-vector- pronuturr-rruuttt 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 result shows that by taking a numerically exact gradient the Lanczos iteration, and leaving everything else to 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 GPy JPyyor'''s is domain- -ppaifc gradients.
6
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) problem that is not a log-determinant and (ii) for a non-symmetric matrix which requires Arnoldi instead of Lanczos, we learn the coefficient field of
02 D
0²
82
ultt 22) w(z), x22²²
u(t: x2, dzp
+
u(t: 21, 22)
(20)
subject to Neumann boundary conditions. We discretise this equation on a 128 x 128 grid in space and transform the resulting 1282-dimensions. second-order ordinary differential equation into a first-order differential equation, ú Aw, w(0) = win with solution operator w(t) = exp(Att (At)wi- The system matrix sparse, asymmetric, and has 32, 768 rows and columns. We sample a true w from Gaussian process with a square exponential kernel and generate data by sampling 256 initial conditions and solving the equation numerically with high precision. Details are in Appendix
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 similar reconstructions with fewer parameters but use 250,000 to display how gradients of the Amoldi iteration scale to many parameters. We compare an implementation of the solution operator (6, wp) H (xp(A(()) Jwo with the Arnoldi iteration to Diffrra''s 's [54] implementation of "Dopri5" [84, 85] with differentiate-the as "Tsit5" [87] with joint (recommended by [53, 54)). All methods receive
adjoint [86] as well ad-
Truth
Salver: Amoidi (ouns)
1.0
0.034 1.0
0.014
0.022
equal matrix-vector products per simulation.
5.006
0.008
0.5
0.038
0.036
Analysis: All methods train, but Arnoldi is more ac-
5.004
0.04
3,002
0.022
curate for fixed matrix-vector vector--oo-ot budgets We
2.5
1.0
evaluate the approximation errors in computing the values and gradients of a mean- -quuree error loss for all three vors and then use the solvers to train the
Solver: Dopris
Solver Tat5
0.012
0.012
5.030
0.010
MLP. We are looking for low approximation errors for few matrix-vector products and for a good recon- struction of the truth. Figure 5 shows the results. The Arnoldi iteration has the lowest forward-pass and gradient error, but Table 4 demonstrates how all ap-
0.5
0.000
3.008
0.5
0.030
0.036
0.092
0.032
0.0
0.000
9.0
0.000
Figure 4: All methods find the truth.
8

---

Table 4: All three methods reconstruct the parameter well deviaaiisss exceed differences for test- loss and RMSE), but Armoldi and DopriS are faster than Tsit5. Dopris uses the Backsolveddjoing and Tsit5 the Reeursivelhhciioonnktjjjooo in Diffrax [54]. We contribute Arnoldi's adjoints.
Arnoldi (adjoints; ours)
DopriS (diff. + disc.)
Tsit5 (disc. - diff.)
Loss on test set Parameter RMSE
6.1e-03 = 3.3e-04
6.3e-03 + 5.7e-04
5.9c-03 + 2.2c-04
2.9c-04 4.4e-05
2.6e-04 = 5.0c-05
2.7c-04 = 5.2c-05
Runtime per epoch
7.7e-02 + 1.8e-05
7.2c-02 + 3.4e-05
2.7e-01 + 1.1e-05
proaches lead to low errors see also Figure 4. The adjoints of the Arnolii iteration match the efficiency of the differentiate- discretise adjoint [86]. and both outperform the
well test set (a held-back percentage of the training data);
the---ifffefeeeetttteee adjoint by margin. This shows how lincar-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 function of a matrix on a problem that is native machine learning: marginal likelihood optimisation of a Bayesian neural network (a high-level introduction is in Appendix J).
Setup:
Laplace-aporoximation of VAN pre-trained on ImageNet We consider as (xx)
"Visual Attention Network" [88] with 4,105, 100 parameters, pre-trained on ImageNet [89]. We assume p(0) NN0,0 2T), and aapoxac-app the log-marginal likelihood of the data
log p(y x)) x log (ly, 0 x) log (et())()) + const
(21)
2
where A(ar) is the generalised
matrix (GGN) from Section 1 (recall Equation 2).
We optimise a via Equation 21, implementing the log-determinant
xI0
via stochastic trace estimation in combination with Lanczos it-
1.0
cration (like in Section 5). Contemporary works [31, 37-433 rely on sparse approximations of the GGN (such diagonal or KFAC approximations), so we compare our implementation to a diagonal approximation of the GGN matrix, which yields closed-form log-
Neg. Log Log-marg.
Diagonal Lanzzos
determinants. The exact diagonal of the 4-million-colle matrix would require 4 million GGN-vector products with unit vec- tors, and like Deng et al. [90], we find this too expensive and resort to stochastic diagonal approximation (similar to trace estimation; all details are in Appendix J). We give both the stochastic diago- nal approximation and our Lanczos-basse estimator exactly 150 matrix-vector products to approximate Equation 21. We compare the evolution of the loss function over time and various uncertainty calibration metrics. Figure 6 demonstrates training and Table 5 shows results.
GGN
0.4
0.2
100
Epochs
Figure 6: Lanczos vs diagonal approx. for a Bayesian VAN.
Analysis: Lanczos matrix-vecto products better (by a margin) The results suggest how, for a fixed sertix-veeoo---rooeee budget, Lanczos achieves a drastically better likelihood at a similar computational budget and already shows significant improvement with a much smaller budget. Lanezos outperforms the diagonal approximation on all metrics except ECE. The subpar performance of the diagonal approximation matches the observations of Ritter al. [31]; see also [14]. The main takeaway from this study is that differentiable matrix-free free linear algebra unlocks new techniques for Laplace approximations and allows further advances for Bayesian neural networks in general.

---

Table 5: Lanczos outperforms the diagonal approximation for calibrating Bayesian version of an Im-
ageNet 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] the out-of-distribution data; mean and -dd.ddvvatsoons of 3 runs. "MVs": matrix-vector products.
Lanezos (50 MVs)
Lanezos (150 MVs)
Diagonal (150 MVs)
Runtime (sec) Marginal likelihood (log)
950
4674
4314
-192,757
-15,,475
-876,444
Joint likelihood: train (log)
-81.2 +: 15,4
-81.2 +: 15,4
-5,669.2 = 124.1
Joint likelihood: (log)
-66.5 + 11.2
-66.5 + 11.2
-5,260.9 + 290.2
Expected calibration error
0.5 +
0.5
0.2 0.003
AUROC cout-of-dist)
T
0.9 + 0.03
0.9 = 0.03
0.5 = 0.010
Acknowledgments and Disclosure of Funding
This work was supported by research grant (42062) from VILLUM FONDEN, The work was partly funded by the Novo Nordisk Foundation through the Center for Basic Machine Learning Research in Life Science (NNF200C0062606), This project received funding from the European Research Council (ERC) under the European Union's 's Horizon programme (grant agreement 101 125993).
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-time differential Lyapunov equation with applications to Kalman filtering. IEEE Transactions on Automatic Control, 60(3):632-64 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. Pre- conditioning 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, and Bernhard Schölkopf. Stochastic marginal likelihood gradients using neural tangent kernels. In International Conference on Machine Learning, pages 14333-14352 PMLR, 2023.
Jens Behrmann, Will Grathwohl, Ricky TQ Chen, David Duvenaud, and Jörn-Henrik Jacobsen. Invertible residual networks. In International Conference on Machine Learning, pages 573- 582. PMLR, 2019.
[10] Ricky TQ Chen, Jens Behrmano, David K Duvenaud, and Jörn-Henrik Jacobsen. Residual flows for invertible generative modeling. In Advances Neural Information Processing Systems, volume 32, 2019.
10

---

[11] Geoff Pleiss, Jacob Gardner, Kilian Weinberger, and Andrew Gordon Wilson, Constant-tinn predictive distributions for Gaussian processes. International Conference on Machine Learning, pages 4124-4123, PMLR, 2018.
[12] Chris Finlay, Jörn-Henrik Jacobsen, Levon Nurbekyan, and Adam Oberman. How to train your neural 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- products for second-order gradient descent. Neural Computation, 14(7):1723-1735, 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-20111 2021.
[15] Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in StatisticsSSmmulatii and Computation, 18(3):1059 1076, 1989.
[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-122) 2023.
[19] Ke Wang, Geoff Pleiss, Jacob Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gor- don 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 Gardner. Large- scale Gaussian processes via alternating projection. International Conference on Artificial Intelligence and Statistics, 2024.
[21] Josie König, Max Pfeffer, and Martin Stoll. Efficient training of Gaussian processes with tensor product structure. Preprint on ArXiv:2312.15305 2023.
[22] Alexander James Davies. Effective implementation of Gaussian process regression for machine learning. 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 Gaussian processes and Bayesian optimization. Advances in Neural Information Processing Systems, 01-22268-22288 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 Tomezak, 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. 1965):152-1174, 1998.
[30] Alexander Immer, Matthias Bauer, Vincent Fortuin, Gunnar Rä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 for neural networks. In International Conference on Learning Representations, volume 6, 2018.
[32] Gene H Golub and Gérard Meurant. Matrices, Moments and Quadrature with Applications, volume 30. Princeton University Press, 2009.
[33] Yaron Lipman, Ricky TQ Chen, Heli Ben- Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In International Conference Learning Representations, 2023.
[34] Aditya Ramesh and Yann LeCun. Backpropagation for implicit spectral densities. In Advances in Neural Information Processing Systems, 2018.
[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 Allingham, Javier Antorán, and José Miguel Herná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 Sun- daram, 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 distribu- tions. Advances in Neural Information Processing Systems, 1990.
[41] Hippolyt Ritter, Aleksandar Botev, and David Barber. Online structured Laplace approxima- tions for overcoming catastrophic forgetting. Advances in Neural Information Processing Systems, 31, 2018.
[42] Wesley J Maddox, Gregory Benton, and Andrew Gordon Wilson. Rethinking parameter counting 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-of- distribution detection for deep neural networks. In Uncertainty in Artificial Intelligence, pages 1968-1967. PMLR, 2021.
[44] Marco Miani, Frederik Warburg, Pablo Moreno-Muñ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-192) 1997.
[46] Simo Särkkä and Arno Solin, Applied Stochastic Differential Equations, volume 10. Cambridge University Press, 2019.
12

---

[47] Severi Rissanen, Markus Heinonen, and Amo Solin. Generative modelling with inverse heat dissipation. In International Conference on Learning Representations. 2023.
[48] Cleve Moler and Charles Van Loan. Nineteen dubious ways to compute the exponential of matrix. SIAM Review, 20(4):801-836, 1978.
[49] Cleve Moler and Charles Van Loan. Nineteen dubious ways to compute the exponential of matrix, twenty-five years later. SIAM Review, 45(())3-49, 2003.
[50] Awad H Al-Mohy and Nicholas Higham. Computing the Fréchet derivative of the matrix exponential, with application 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 for unfferential-algebran equations: The adjoint DAE system and its numerical solution. SIAM Journal on Scientific Computing, 24(3):1076-108) 2003.
[52] Moritz Diehl and Sébastien Gros. Numerical optimal control. Optimization in Engineering Center (OPTEC), 2011.
[53] Christopher Rackauckas and Qing Nie. Differential&guationss.--- performant and feature-rich ecosystem 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 Vaib- hav Dixit. FiuxEl julia library for neural differential equations. arXiv preprint arXiv:1902.02376 2019.
[56] Yingbo Ma. Vaibhav Dixit, Michael Innes, Xingjian Guo, and Chris Rackauckas. A com- parison of automatic differentiation and continuous sensitivity analysis for derivatives of differential equation solutions. 2021 IEEE High Performance Extreme Computing Confer- ence (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 the solution of the matrix eigenvalue problem. Quarterly of Applied Mathematics, 91):17-29, 1951.
[60] Matthias Seeger, Asmus Hetzel, Zhenwen Dai, Eric Meissner, and Neil D Lawrence. Auto- differentiating linear algebra. Preprint on ArXivv1710.087717, 2017.
[61] Sebastian Walter and Lutz Lehmann. Algorithmic differentiation of linear algebra func- tions with application in optimum experimental design (extended version). Preprint 1rX1v:1001.1654 2010.
[62] Andreas Griewank and Andrea Walther. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiaaion. 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 Davis and Yifan Hu. The University of Florida sparse matrix collection. ACM Transactions on Mathematical Software (TOMS), (8(1):1-12))) 2011.
[65] lain 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 ar31 14606, 2024.
13

---

[67] Gerrit C Groenenboom and Henk M Buck. Solving the discretized time-independent Schrödinger equation with the Lanczos procedure. The Journal of Chemical Physics, 92 (7):4374 43799, 1990.
[68] Christopher Conway Paige. The computation of eigenvalues and eigenvectors of very large sparse malrices. 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-385 1972.
[70] Christopher C Paige. Error analysis of the Lanczos algorithm for tridiagonalizing symmetric matrix. IMA Journal of Applied Mathematicx, 18(3):341-34,, 1976.
[71] Steffen Börm and Christian Mehl. Numerical methods for Eigenvalue Problems. Walter de Gruyter, 2012.
[72] Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe Llinares- López, Fabian Pedregosa, and Jean-Philippe Vert. Efficient and modular implicit differentiation. Advances in Neural Information Processing Systems, 35:5230-5220 5242, 2022.
[73] Michael Betancourt, Charles C Margossian, and Vianey Leos-Barajas. The discrete ad- joint method: Efficient derivatives for functions 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 backpropa- gation algorithms for square, wide, and deep-real or complex-matrices and their software implementation. Preprint on ArXiv:2009.100 10071, 2020.
[76] Thomas P Minka. Old and new matrix algebra useful for statistics. www.:/aw.......ddd///iinaappppee
See
4, 2000.
[77] Mike Giles. An extended collection of matrix derivative results for forward and reverse mode automatic differentiation. 2008.
[78] Jean Céa. Conception optimale ou identification de formes, calcul rapide de dérivée directionnelle de la fonction coût. ESAIM: Modélisation mathématique et analyse numé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-Men Milne,, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL http:// gi thub. com/google/jer
[80] Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS-GP). In International Conference Machine Learning, pages 1775-1784. PMLR, 2015.
[81] Helmut Harbrecht, Michael Peters, and Reinhold Schneider. On the low-rank approximation by the pivoted Cholesky decomposition. Applied Numerical Mathematics, 62(4):428 440, 2012.
[82] Benjamin Charlier, Jean Feydy, Joan Alexis Glaunes, Franç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. family of embedded Runge-Kutta formulae. J. Comp. Appl. Math, 6:19-26, 1980.
14

---

[85] Lawrence F. Shampine. Some practical Runge-Kurtt Kutta formulas. Mathematics of Computation, 46(173):135-150, 1986. doi: http://di..rr/110223300000000000
[86] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in Neural Information Processing Systems, 31, 2018.
[87] Charalampos Tsitouras. Lunge-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-Nimg Liu, Ming-Ming Cheng, and Shi-Min Hu. Visual attention network. Computational Visual Media, 94):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 255. IEEE, 2009.
[90] Zhijie Deng, Feng Zhou, and Jun Zhu. Accelerated linearized Laplace approximation for Bayesian deep learning. Advances in Neural Information Processing Systems, 35:2695-2708 2022.
[91] Bolei Zhou, Agata Lapedriza, Jianxiong Xiao, Antonio Torralba, and Aude Oliva. Learning deep features for scene recognition using places database. Advances Neural Information Processing Systems, 27, 2014.
[92] Tianlin Liu, Jose Antonio Lara Benitez, Florian Faucher, AmirEhsan Khorashadizadeh, Maarten V de Hoop, and Ivan Dokmanić. WaveBench: Benchmarks datasets for model- ing linear wave propagation PDEs. Transactions on Machine Learning Research, 2023.
[93] Peter Deufthard, Andreas Hohmann, Peter Deufthard, and Andreas Hohmann. Three-term recurrence relations. Numerical Analysis Modern Scientific Computing: An Introduction, pages 151-178, 2003.
[94] Peter Deuftaard On algorithms for the summation of certain special functions. Computing, 17:37-48, 1976.
[95] Simon Bartels, Kristoffer Stensbo-Smidt, Pablo Moreno-Muñ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-998: 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: software environment for integrated models of biomolecular interaction networks. Genome research, 1344,)444882254, 2003.
[99] P. J Prince and R. Dormand. High order embedded Runge-Kutta Kutta formulae. Journal of Computational and Applied Mathematics, 7(1):67-75, 1981.
[100] Antonio Stanziola, Simon R Arridge, Ben T Cox, and Bradley E Treeby. j-Wave: An open- source 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 Bayesian neural nets via local linearization. In International Conference Artificial Intelligence and Statistics, pages 703-711. PMLR, 2021.
[103] Vardan Papyan. Traces class/cross-class structure pervade deep learning spectra. Journal of Machine Learning Research, 21(252);1 2020.
15

---

[104] Costas Bekas, Effrosyni Kokiopoulou, and Yousef Saad. An estimator for the diagonal of matrix. Applied Numerical Mathematics, 57(11-22)12144-222, 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.
16

---

Appendix: Overview
Some of the results in the main paper promised detailed information about setups, data, compute, or additional 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 respectively; Appendices B 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), and Appendix 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 numerical algorithms - at the heart it are our reverse-mode differentiable Lanczos and Amnldi iterations. We provide JAX code to reproduce all experiments at the URL
attps://githune
and have packaged all numerical methods in a JAX library that can be installed via
pip install matfree
Next to Lanczos and Amoldi, this includes variants of conjugate gradient methods [58] and pivoted Cholesky preconditioners [81] (which we used for Gaussian processes), low-memory kernel-matrix- vector products (also Gaussian processes), efficient GGN-vector products (used for Bayesian neural networks), and matrix-free sampling algorithms to sample from Gaussian processes and Laplace- approximated 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 software implementation in the current JAX ecosystem (with some exceptions relating conjugate gradients).
Compute All experiments before the case studies were run on CPU. The Gaussian process and differential 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 "besstk18" 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. www.eeeepemmetttee sparse, and use 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 and runtime for Lanczos and Arnoldi in Figure 7 (this includes the curves from Figure again). We see
Figure 7: Run and compilation times for Lanezos and Amnoldi The "backprop" curves were stopped at 100, because for higher values we encountered memory issues. All experiments run on CPU.
Lanczos: Run time
Lanczos: Compile time
Amoldi: Run time
Amoldi: Compile time
25
Wal time cred)
MOJ Umo Goc))
time (1NC))
WSY time crect
0.2
30
0.1
0.3
130
150
200
100
150
200
50
100
159
200
250
50
100
159
200
250
Kryfor-span Kyslor-spaca
Kryfov-spaca
different behaviour for Lanczos and Amolli. Whereas for back-propagation through Lanczos without
17

---

the custom gradient, the compilation times remain constant for increasing Krylov-space depth K and runtimes increase rapidly, the reverse is true for the Arnoldi iteration. The adjoint method mirrors that of the forward pass in both benchmarks. In either case, the increasing memory requirements for backpropagation 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 Lincarisation
The linearisation of Equation is
(dA)Q + AdQ (dQ) - QdH - (ex) T =0 VR
(22a)
(dQ)e1 (dv)c odc = 0 Nxl
(22b)
Ig (dQ)es + (Qej, (dQ)
ERKxK
(22c)
I< dH
E BR
(22d)
[rr, (dQ)ez) + (Qe,. dr> = To see this, apply the chain- and product rules the original constraint in Equation 12.
R
(22c)
B.2 Transposition
Let p = p(Q,r, r, H, e) E IR be a scalar function of the outputs. the following, interpret vectors as N x 1 matrices and scalars are the only scalars.
x matrices. The values of inner products and the realisations of p
For all A E R.N
A Rx1 T,E R and 7 E Rx11 we have
dp = (Voa, dQ) + (V dHf + , + (Vo,, de)
(23a)
(VQp, dQ) + (V |H] dr),dr) <Ve,, de)
+ (dA)Q AdQ
QdH (er)
+ -dd)) wdc)
(dr)c
(T, o [(Qe, (dQ)e,) (Qej, (2)));)
+ oHH)
+ (v. (dQ)e,) dr)
(23b)
-Voo.dQ) + (VHP, dll) dV),, + (Ve,,d)) de)
+ (AQQ, + (AT A,dQ) - (AH ,dQ) - (Q A, dH) - Ae),dr)
+ (d(e)),,dQ)) (AcT, ,dv) T
+ (0(1 ()), + (Q(1s o T, ,d))
+ o,
+ ( dQ) + (Q7,dr)
(23c)
= =ZQ (Zn, dlf) + (Zr, + (Ze de) + (AQT, ,dA) + (AcT, ,d)
(23d)
with the constraints
Zo = Vop + T AH + + + Q11 Q(1 o +r7
NxK
(24a)
Zμ = VEР VEP- - Q
1 1
E kx K
(24b)
2, = -ep-Aex + Q7
E B
(24c)
Z. T
R 1xl
(24d)
Solving the adjoint system, Zo = 0, Zh 0, Z, 0. and 2, 0 as function of A, A, T, Σ, and
yields the desired VAP AQ and V.f Ac Theorem 4.1 is complete.
18

---

C Proof of Theorem 4.2
B.r...,, agi bi....,kk' shall be scalar/loss that depends on the output of the algorithm. Denote by V the gradient of with respect to each Lanczos vector zk. and by Var and b. the gradients with respect to aj and by respectively.
Recall that p ,(2 IK+1i ai,...
The differential of normalisation, i.e., the operation 8 + h = 8/(88 T 8) is
dh S (1 h ds.
(25)
The next steps are the usual ones: 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.
Define Lagrange multipliers {hx C R and {ux/2, 11 S R,
p = p+
(- 1 (A axT) bL.r+11)
(by 1, 2o = 0)
K
K
Ao
I1
N& (x
#441
-
+
UR JR IR+1-
(26)
VuTu
2
k=1
k=1
Differentiate and use that the forward constraint must be satisfied,
K+1
K
K
dp
sx p)
d +
(Vaz))
da Σ (VDp)) Tdbx
das-
dbb
1=1
k=1
Ao [(dA)== - (dax)ar
(dbs 1, drg = 0)
+
Ao [Adrı axdæx - bb-iddk-1 - 1
(b) = 1, x = 0)
1=1
Ao (dæı I1- )/(0 w)dd))
Pk 2 dox+1
+
Σ
Uk (aIziidd dar + oht dz
(27)
1=1
Sort all terms by differential,
K+1
dp- Σ Zog dag
Σ Zo, db
Σ k=1
Zo,
dxx + Zodv + trace (ZadA).
(28)
so that enforcing that all Zau, Zbu, Zzu terms are zero yields constraints for the multipliers from which we can compute gradients with respect to v and A.
What are those terms? Let Ax+1 = 0, pD = 0, and in = 0 (to simplify notation below); then,
Zux+1 -xxb + (V V P + PK + #кナK),
(29)
and for all = (=KKK... K,..,,
(recall xo = 0, by = 1, 0, vo = 0, Ax+1 =
0)
Ze. = -bsde+1 + (AT - ax1)dx b + (VagP + pk-1 + 14 + 1
(30a)
Za. = Vaaf - -XT TK.
(30b)
Zn. = VoaP ATIIIE - Al Z2+1-
(30c)
19

---

The expressions are like the forward-Lanczo constraints, and the main differences are (i) that the recursions are rum backwards "time" and (ii) the existence nonzero bias term in the adjoints (marked by parentheses). Enforcing 7. Zou Zen, to be zero identifies
Ag 11 21
1o )
op Z
VAP ZA
J Σ A=1
(31)
Theorem 4.2 is complete.
D Proof of Corollary 4.3
The following derivation covers only the case for Lanczos, i.e., the variables T1,,,.;; 2K+1 instead of Arnoldi's 's 91, ЧK- But the derivation is the same for both methods.
The expression we manipulate is
dp trace
A
aA
+ const
(32)
k=1
where all non-dA-related quantities are treated as some unimportant constants.
For single parameter fg. dA (Vo, A) d0, and have
K
K
trace
Σ x
DjAd0,
trace
Σ
ARA
A)
d0,
(chain rule)
K
= trace
X (Vo, (Vo,A) A)TxA xx
d0j (cyclic property of traces)
k=1
Σ
No [o; A) 2gd0
(a scalar is its own trace)
K
Vo,
X AT
d8,
(linearity diff. & summation)
K
Va,
AT
Acr
d0j.
(symmetry of A)
In conclusion, the derivative of P wrt Oj is
Vo V,
x A
A
(33)
and stacking all of those partial derivatives on top of each other, we obtain
K
Vap V e, plj
V
0
Σ 2 A(O)
T A
EV A(0) An
(34)
A=1
1=1
We already compute A(0)T > during the backward pass, so we are single vector-Jacobian product with away from a matrix parameter-gres gradient instead of matrix saarix-grrdient.. This requires O(p + n) storage, and is computed online, which makes the memory- complexxyy 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 solution is technical, and because starting with Lanczos can provide template for solving Arnoldi's adjoint system. All results in the present section (except for those that explicitly point to Deuflhard et [93], Deufthard [94], which are marked as such) are new and contribution of this work.
20

---

E.1 Lanczos
The inputs the adjoint system are the Lanezos vectors {= and the coefficients {ax fk as well as {a k= from the forward pass, the corresponding input derivatives (VagP) K P) and {V VapPhem and matrix A and initial vector v.
The overall strategy for solving the adjoint system of the Lanczos iteration (Theorem 4.2) is the following: for every 1 K, ..,, alternate the two steps:
1. Combine orthogonality with the Za constraints to get Vk, and combine it with the Zo constraints to get ja-
2. Once each μa and are available, solve for d and repeat with the next lower k.
This results in the following procedure: To start, set (K+1 -(V xx and Ax+1 0. Then, for all K ...,T 1, compute
Ex Sk+1/bx
(35a)
Ab Vo,P 14 x
(35b)
Va, xx _/EC. &x
(35c)
A Ex + 1 Ik+1 + V/
(35d)
CA vnp A A + au Au + by AR+1 bb ix Tk-1
(35e)
Repeat with
(35f)
Finally, set Àg = Ç1. Only À and ( affect subsequent steps; μ. V, and & are only needed for computing À and Ç. The k-tth step depends on ax: bb, Tk, 2k+1, Vo, P, Vo P, V P.
The strategy above yields all {AR A Finalise the gradients
AT I1 U
Ao.
VAP=
A
(36)
which can be embedded into any reverse-mode mode algeritmii-iifferetitttionnn
E.2
Arnoldi
The 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 perspective on recurrence relations (which include the Arolli 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 Amoldi vectors amounts to solving
AQ + QH r((x) T =0,
(37a)
Qei cv = O,
(37b)
which is possible in closed form as follows: We rewrite the first two constraints as
(el @ T)vec (Q) vec (cv)
(38a)
- & A)vec (Q) (H" & L)vec (Q) (ex @ I)vec (r) 0.
(38b)
This expression equivalent
@HeI-1OA
(H 6 1
eK ex@1 @
0 (vee (Q) vec (r)
vvec (co) (00ec((vr))))
0
NIK+HI)
(39)
This is similar to the work by Deuflhard [94], who explain adjoints of three-term recurrence relations. Since ² = (v.e)) holds (ie, i known), the first row of is known. Then, the first row of Q together with the orthogonality constraints yields the first row of HT, which then defines the next row of the linear system thus the second row of Q. Alternating between deriving the next row of H and solving the lower triangular systems is then Amoldi's algorithm:
21

---

Algorithm E.1 (Amoldi's forward pass: paraphrased), Assume that e and K are known. Compute V(v, v). Then, for k = 1,... K, alternate the following tsteps:
1. Derive the next column of H using the orthogonality constraints.
2. Forward-substitte ("solve") the block-lower-triangglar system for the next column of Q (respectively r at the last iteration).
Return all Q and H, as well as 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 Zo = 0 and Z, 0 mirror those of AQ QuI r(ex)) = 0 and Qel - cv = 0;
the constraints Zи = 0 and Ze = 0 mirror the orthogonality constraints Q Q = =I and Q r 0. Therefore, we start with ZQ and Z, O,
Vos + A A + X(e)) + Q[1 oT) Q(1 oT)
T +ra =
(40a)
v. - Aex + Q7 = 0.
(40b)
Introduce the auxiliary quantities
(T, 7) Vof + Q41 oT) Ql1 r) + my (()) = v,f -
E R
(41a)
E Nx1
(41b)
I and &i only serve the purpose of simplifying the notation in the coming part; there is no "meaning" associated with them. Vectorise both expressions,
(1 x A )vec (A) - [(H T & llvec (A) + (el @ L)vec (A) = - vec (V((I,)))
(42a)
[(ex) @ L]vee (A) = vec (N((7)
(42b)
and observe that this can be written as a linear system
(1&T &I I - -TTTHHTT ( @ 2) 0
( (voc(())
vec (A)) (A)
(-voc (T 7))
- vec (v((1))
(43)
(ex) @
with a system matrix that is the transpose of the system matrix of the forward pass. The matrix upper triangular, and the equation can be solved with backward substitution provided v1()) and (T) are known.
The defining quantities V and vi emerge by combining the adjoint recursion with the projection constraints Z = 0 (for T which yields V) and Z, 0 (for 7, which yields v/). We use Z, to get single element in I'; more on this below. Summarise the adjoint pass:
Algorithm E.2 (Arnoldi's adjoint pass; paraphrased). Assume Q, H, c, and well as the gradients of f with respect to those quantities. Then, compute vú via computing using Z, = 0. Then, for k K..., 1, alternate the following two steps:
1. Derive the next row of V by combining Zo 0 with the projection constraint Zu =0
2. Backward-substiutee ("solve") for the next row A (recall: we loop backwards)
Finally, use Ze to get the first row of W and solve for Then, return Va = AQT and v. Ac
The structure of the adjoint pass is similar to the forward pass (Table 6). In the following, we will elaborate on each those steps. We assume that the reader knows how solve a lower triangular linear system. We focus constructing & and v via I' and 7.
22

---

Table 6: Forward versus adjoint (backward) pass
System matrix
Solve via
Recursively define
Using
Forward Adjoint
Lower triangular Upper triangular
Forward substitution Backward substitution
System matrix Right-hand side
Orthogonality Projection: Zи =
E.2.3
Initialisation
Initialisation of the adjoint pass implies computing vb. To get v we need y: Consider multiplying Z, with Q
0 = Q Z, = -'v +
(use the definition of Z,)
+7 = (Q'v.J-Vn.fee + (Iæ o E)EK +7
= =Q'v.. (Vv/ Le Ek +7
(since Zμ 0) (reorder)
= =QQTV, VuSex)++ +7
(U& Σ)ex = 0)
where we use that the last column in T& Σ consists entirely of zeros. All other quantities are known. Therefore, y isolated and we identify
7 Vnfex Q V
(44)
Next, use this y build 2b.
(67) = V,f+
(45)
and the initialisation step is complete.
E.2.4 Recursion
With v in place, we get the last column of A. To get the next column of A, we need derive the right-hand side V. To get V, we need l.
Multiply QT Zq to obtain
= QT Zz
(46)
= =Q + - +1o
(since QTQ I and Q
E)HT + (ee)(e))TT + o + o) QQ + A - (V
(since ZH = 0)
V +Q A M(r) + +go (Ig (I& H)HT (reorder the terms)
Since H is Hessenberg, (L& o H T is strictly lower triangular. Therefore, multiplication with 1> removes E from the expression,
Vol = 12 T - V + T A Q Men) +
(47a)
=1 Q Vof - T +Q A + +12 . |
(47b)
The term involving À can be simplified as follows: Due to the presence of (e we know that Q >(e)) is lower triangular (in fact, it has single nonzero column). Thus, I> o (Q X)er)T is proportional to Te
1 (Q ((Qe)))Tees (()((())T
(48)
cu -ee(ee
(since Qei = cv)
=cV.fer(er)
(since Z. 0)
and all quantities are known; hence,
I_ T + Ia TT -Ig 1Q WHSH -
(49)
23

---

Algorithm E.3 (Forward pass). Initialise 1, Qe1- Then, for 1..... K:
Algorithm E,4 (Backward pass). Initialise
k K, Aeg. Then, for k K,.., 1:
Use orthogonality for a new row in the system matrix in Equation 39,
Use projection for new column of the right-hand side y
2. Solve for the next column of Q
2. Solve for the next column of A
3. Optional: re-enforce QTQ
3. Optional: re-enforce ZH 0.
Solve for and return Q, r.
Solve for À and return Vop and V.p.
Figure 8: Forward and backward pass Arnoldi iteration (paraphrased)
must hold.
Now, the most important observation the following: the last column of -= o1' 12 o depends on the penultimate
the last column of Q' A A and known quantities; the penultimate column depends column of r, and so on. But at the time of assembling the last column of I, the last column of A is known! More generally, we always know one more column of than of r, so we can recursively assemble I= oo + +2 :
Let
sym (M) 12 a M + (I> o M)
(50)
be a symmetrisation operator. We define it for the sole purpose reconstructing
sym of = - (lgo))"
(51)
Let us use it:
1 of oT) = sym (-1- QQ Vol - VHfH + <e (52a) =ym (-Iz Q Vo V
(52b)
This yields the next row/column of PT and therefore the next row of V. From there, we can assemble the next column of A 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, obtain not the last column of A but X, though this byproduct of solving the triangular linear system. It does not need further explanation.
Remark E.5 (E). Like for gradients of QR decompositions [61, 75J, we never solve for E.
Setup for Table 2
To create Table 2, we implement operator
A (H,Q, QHQT
(53)
of accuracy
10-2
regroj.
Adi we re Baskorsp
where (H, Q,,,)),) are the result of full-rank Aroodi
10-6
iteration (i.e. K = N). For K = N, QHQ = and I must have identity Jacobian; thus,
10-10
10-14
8 = N T||RMSS
(54)
10-18
measures the loss of accuracy when differentiating the Arnoldi iteration. A small is desirable.
2
Hilbert matrix size
Figure Accuracy loss e when differen- tiating I for Hilbert matrix of increas- ing size N. Uses double precision.
Then, using double-precision, we construct Hilbert matrix A [1/(i + + D) - NxN which a fa- mously ill-conditioned matrix and a common test-bed for the loss of orthogonality in methods like the Lanczos and evaluate three algorithms, all of which rely on the Arnoldi iteration with full reorthogonalisation on the forward-pass: One algorithm does not re- re-project the adjoint constraints, another does, and
Arnoldi iteration [e.g. 71, Table 7.1]. We

---

Matvec run times (sec) vs. size of dataset
Glybrihh (oia Kosps)
Ours jeia mec!
10-1
Ount jeia vman) K-10 Oun jvia vman) K-100
Run time (sec)
10"h
10"3
10a
10a
Dataset size
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 dimensions.
for reference compute when "backpropagating through" the re-orthogonalised Amoldi iteration as a third option. Figure has demonstrated that the first two options beat the third one terms of speed, but consider numerical accuracy here.
We evaluate & for N 1, ..., 8 (see Figure 9), and show the values for N Table 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-veetor products in JAX
Matrix-fres free linear algebra requires efficient matrix-vecte products. For kernel function k: k(x, x²), and input data ø1, N, Gaussian process covariance matrices are of the form A = [k(I, 22 [2,]1, [))]"===> Matrix-vector products with A thus look like
N
N Σ
U - Av =
k(a,, uz
(55)
and can be assembled wise,ise, 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 al. [6] and choose the largest number of rows of A that still fit into memory, say such rows, and assemble Av in blocks of r. In practice, we implement this in JAX by combining jax. lax. map and jax. vmap, but care has to be taken with reverse-mode automatic differentiation through (v, 0) + A(0)c because by default, reverse-mode differentiation stores all intermediate results. To solve this problem, we place checkpoints around each such batch of rows, which reduces the memory requirements but roughly doubles the runtime. (We place another checkpoint around each stochastic trace eatemmsstittionn 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 com-cmmpaeelle interface to KeOps which why we have to implement the above solution.
Figure 10 compares the runtime of our approach to that of KeOps custom CUDA code. We that we are competitive, but roughly 5x slower for medium to large datasets. Multiplying this with the 4x increase due the checkpoints discussed above explains the 20x increase in runtime compared to GPyTorch. Being 20x slower than GPy 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 KeOps-to-JAX interface.
H Experiment configurations for the Gaussian process study
Data For the experiments we use the "Protein". "KBGG (undirected",, "KEGG (directed)". "Eleva- tors", and "Kin40k" datasets (Table 7, adapted from Bartels et al. [95]). All are part of the UCI data
25

---

Table Datasets used in this study.
Dataset
Source
Protein
Available here."
Elevators
Camachol [96]
Kin40K
Schwaighofer and Tresp [97]
KEGG (undir)
Shannon et al. [98]
KEGG (dir)
Shannon et al. [98]
repository, and accessible through there.
The data is subsampled admit the train/test split of 80/20%, and to admit 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., KEGG (undir), we use 63,600 instead the original 63,608 points.
We calibrate a Matérn prior with smoothness v/ = 1.5, using 10 matrix-vector products per Lanczos iteration, conjugate gradients tolerance of e 10 Rademacher samples. We evaluate all samples sequentially (rematerialising on the backward pass to save memory, as discussed in Appendix G). The conjugate-gradient tolerances are taken to be absolute (instead of relative), and the parametrisations of the Gaussian process models and loss functions match that of GPyTorch
rank-15 pivoted Cholesky preconditioner, and
For every model, we calibrate an independent lengthscale for each input dimension, as well as an scalar observation noise, scalar output-scale, and the value of 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
8
82 Fut
2²
821 u ((1zz1100 = x2),
(1, + u(t 11,
(56)
Ou
with Neumann boundary conditions. The coefficient field w is space- but not time-dependent.
First, we discretise the Laplacian operator with central differences on equidistant, tensor-product mesh that consists of 128 points per dimension, which yields 128" grid points. The resulting second-order ordinary differential equation
d²
w = w Mu,
(57)
dtz
where M is the discretised Laplacian, is then transformed into first-order differential equation
w =
0 w M
0) w Aw.
(58)
di
This equation is solved by the matrix exponential, and the system matrix A is asymmtric (by construction), and highly sparse because M is. Matrix-vector products with A are cheap, because can implement them with wessssscpy... wwips..iinall................... convolve2d.
Then, sample true from Gaussian process with square-exoonential covariance kernel, using lengthscale softplus( (rtplus--0075) (-0.75) and sutput-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(O) and output scale softplus(O), to sample 256 initial distributions again with the Lanczos algorithm
26

---

1.0
1.0
0.75
Inputs
1,0
0.50
0.8
0.25
0.3
0.5
0.5
0.00
0.5
0.0
-0.25
-0.3
-0.75
0.0
0.0
-1.00
0.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
1.0
3.0
1.0
0.5
1.0
6.0
Targets
3.0
0.5
0.5
-1.
-2.0
0.5
-1.5
-0.5
-6.0
0.0
0.0
0.0
0.5
1.0
0.5
1.0
0.5
1.0
Figure 11: Three exemplary input/output pairs from the PDE dataset.
[11]. These 256 initial conditions solved with Diffrax's '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 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 marginal-likelioodd- optimisation of Bayesian Neural Networks. Suppose go(x) is the output of neural network with parameters 0 E I The choice of the model shall be denoted by M and consist of both continuous and discrete hyperparameters (such network architecture, likelihood precision, prior precision, etc.). For some choice of prior given by
p(0 MM)
(59)
and likelihood
plylæ, 6, M) ply (.))
(60)
we can specify a Bayesian model. The posterior distribution then given by:
p(0,y M) ply M)p) MM)d.
(61)
The marginal likelihood is given by normalizing constant of this posterior, i.c.
plu x M) = p(u >. go (x)) M)p(0 MM)0.
(62)
As suggested by MacKay [101], this marginal likelihood can used for model selection in Bayesian neural networks. Immer et al. [30] use the Laplace approximation of the posterior to obtain access to the marginal likelihood the Bayesian neural network and its stochastic gradients.
The Laplace approximation of the marginal likelihood is given by:
1 logp(y I, M) = poy, dMAP x. M) 2 log det
1 2TT
Hosae
(63)
where Houee V log plu, OMAP M, Usual choices of the prior are (0,, a -11). Usually this Hessian is approximated with the generalized Gauss-Newton (GGN) matrix [102]
How ss A(a) = E (eDoooneD (W 9o (
+ a²1
(64)
1=1
Link http://archive... Acs.
"To ensure radiating boundary conditions, Liu et al. [92] follow Stanziola et al. [100]'s model of simulating the wave equations sytem of first-order equations.
27

---

where D² is the Hessian of the loss, and Dag the Jacobian of (recall Equation 2). This objective 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 can estimate the log-determinant of GGN using matrix-free techniques like the Lanezos iteration.
To make predictions using the Laplace approximation of the posterior, also need sample from the normal distribution N (OMAP, A )) Samples from this distribution can be written as:
0 MMAP A 1
(65)
where (0, I). The main bottleneck in this computation is the inversion and matrix square root of the GGN matrix, and we implement it with Lanczos iteration using fxx = x 1/2, Since the GGN is empirically known have low-rank [103], doing few Lanezos iterations can get us close to an accurate estimation.
J.2 Experiment setup
We estimate the diagonal of the GGN stochastically via (*o" is the element-wise product) [104]
diagonal (A) = Elu o Au]
Σ
Up o Au, Blue = I.
(66)
We use 150 matrix-vector products for both diagonal calibration and our Lanezos-ban -based estimation. We use 30 Monte-Carlo samples to estimate the log-likelihoods for evaluating the test metrics, and we use places365 [91] as an out-of-distribution dataset to compute OOD-AUROC. We also compute the expected calibration error (ECE) [105] of the model.
Data: We show scalability by doing Laplace approximation on Imagenet 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, even for 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 learning rate of 0.01 for 100 epochs for optimizing both the diagonal and Lanczos approximations of the GGN.
28

---

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 Arnolii 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 the contributions made 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 is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: The limitations are discussed in the paragraph titled "Limitations and future work" on page 4. All assumptions (i.e., Assumptions 3.1 and 3.2) are contextualised the sentences before and after they are have been introduced.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "L_imitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations 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 the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or 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 the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant 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 of assumptions and complete (and correct) proof?
29

---

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 cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof' can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proof' provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: The most important information about the experiment setup part of the main paper in Sections 5 to 7; further information can be found in Appendices G to J. Appendices and F discuss the setup for Figure and Table 2. Code will be published upon acceptance.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is 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 novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, 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 case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature 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) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results or 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 case authors are welcome 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 in some way (c.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
30

---

reproducible, they 5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Question:
Answer: [Yes]
Justification: Code has been submitted as 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 NeurlPS code and data submission guidelines (https://nip) cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: / /nips. wtc/public/guuueeeo
FrodeSubmiseioonoole
for more details.
The authors should provide instructions on data access and preparation, including how 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 subset of experiments should 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 supplemental material (appended to the paper) 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, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
See the answer to "4. Experimental Result Reproducibility" above.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, 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 the Bayesian neural network example, which uses single training run (but evaluates test metrics 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 emror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
31

---

The factors of variability that the error bars are capturing should be clearly stated (for example, 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 the standard error of the mean.
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 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars 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
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: The introductory part of the appendix contains 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 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 estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https:
Answer: [Yes]
Justification: We use datasets that are either self-created (Section and Appendix I), or common test-cases for machine learning methods (UCI datasets, ImageNet) or numerical algorithms (SuiteSparse matrix collection).
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations 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:
32

---

The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deploymentt of technologies that could make decisions that unfairly impact specific groups), 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 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 to generate 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 train models 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 the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (c.g., gated release models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
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 high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere 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 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 and properly 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 URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
33

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperswi thcode cow datasess has curated licenses for some datasets. Their licensing guide can help determine the license of 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 encouraged to reach out to the asset's '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) documented. The appendices contain all 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 part of their submissions 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 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 details about compensation (if any)?
Answer: [NA]
Justification: This study does involve crowdsourcing nor human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS 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 data collector.
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 crowdsoureing nor research with human subjects.
34

---

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 the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
35


Reviews:
[{"review_id": "RnNZm4U9GU", "reviewer": "Reviewer_RkAM", "summary": "This paper proposes a new approach to perform automatic differentiation for function of large matrices. Specifically, the paper outlines the backward computation of the matrix-vector product f(A(\\theta)) * v where A(\\theta) is the jacobian of a large NN that will not fit into memory. The proposed approach uses the Lanczos and Arnoldi iteration method to factorize the matrix-vector product, which permits inexpensive backpropagation, and then derive the backward computation for the Lanczos and Arnoldi iteration via the adjoint method. This technique is tested in three different scenarios that require the evaluation of such matrix-vector product: Gaussian processes (GP), Physics Informed ML with PDE, and Bayesian NN.", "questions": "I have no further question since the method is sound & clearly presented. My only problem is that the empirical results are really limited and do not show any significant improvement over existing methods.", "limitations": "No potential negative societal impacts", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- Theoretically interesting approach to auto diff of large Jacobian vector product\n- Potentially useful in scenarios where functions of large matrices are involved in the objective function\n- Interesting case studies detailing several such scenarios and how to apply the proposed method in those scenarios\n- Improved performance in majority of experiments\n- Generally clear and self-contained presentation, easy to read and follow", "weaknesses": "1. Motivation is somewhat unclear. The paper claims that the proposed method provides exact gradient of the forward pass, but it doesn't seem true when the Lanczos & Arnoldi iterations themselves do not yield an exact factorization (Eq. 7). This method, although theoretically quite interesting, should belong to the same class with other approximation methods. As such, I would like the authors to explicitly discuss the benefit of this approximation compared to previous work (which are currently positioned as being inferior to an exact gradient method)\n \n2. Empirical results are not convincing:\n- Table 3 seems to highlight the wrong improvement. It is moot to compare training loss when the loss landscapes are different. The real performance measure should be RMSE, in terms of which both methods converge to the same values (although the proposed Arnoldi method is 20x slower -- so I'm not sure what's the benefit here)\n- Case study 6 again shows the same performance between the Arnoldi method & Dopri5.\n- Bayesian NN is usually trained with Variational Inference. Diagonal method is a rather crude approximation, and it is not surprising that the proposed method outperforms it.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the review and the positive assessment! \n\nWe would like to briefly reply to the points you list as weaknesses:\n\n**Motivation:** \nYou are correct that the Lanczos and Arnoldi iterations yield approximate matrix-function-vector products. When we write \"exact gradient of the forward pass\", we mean that our adjoints yield the exact gradient of the Lanczos/Arnoldi approximation (the forward pass). We do not claim that it is the exact gradient of the matrix function.\n\n**Table 3:** \nWe agree that the RMSE is crucial. However, since neither method beats the other in terms of RMSE, Table 3 does not highlight any RMSE-winners. For context, kegg_dir also has no winner in the training loss. As for the interpretation of this experiment, we consider it a success that our new black-box algorithm matches the performance of a technique that (i) specialises in Gaussian-process log-determinants and (ii) has been honed for many years. We discuss this in the \"Analysis\" block on page 7. The 20x runtime difference is not due to our algorithm. Instead, it is because the reference code (GPyTorch) relies on hand-crafted CUDA kernels through the KeOps library. The 20x performance difference is due to the inability to use this library in JAX; see the discussion in Appendix G.\n\n**Table 6:**\nHere, too, is the main message that our black-box Arnoldi code can compete with state-of-the-art differential equation adjoints without specialising in ODEs/PDEs. We use the same algorithm for all three case studies, whereas all reference codes specialise in each respective task. For example, the differentiate-then-discretise adjoint for Dopri5 may be fast (as fast as the Arnoldi adjoint for the linear PDE) but only works for ODEs (respectively space-discretised, time-dependent PDEs). Please note how, in Figure A2, our adjoints are more efficient than Dopri5 with a differentiate-then-discretise adjoint in terms of the number of matrix-vector products.\n\n**BNN:**\nWe respectfully disagree with the \"usually\" in your statement that \"Bayesian NN is usually trained with Variational Inference\".\nThere are many strategies for setting up Bayesian neural networks, including the Laplace approximation and variational inference; see Papamarkou et al. (2024).\nHowever, unlike other methods, the Laplace approximation is particularly relevant for our work because marginal-likelihood calibration requires (the gradient of) a function of a large matrix: the Gauss--Newton matrix.\nOur paper does not demonstrate how Laplace approximations compare to other techniques.\nInstead, it outlines how differentiable linear algebra makes Laplace approximations more efficient.\nOur results show how this has been successful: Our unspecialised implementation of matrix functions beats specialised Laplace-approximation codes for the VAN model. We believe that this reinforces the importance of differentiable linear algebra for advances in probabilistic machine learning.\n\n\nThank you again for the review; we hope that we were able to clarify a few potential misconceptions. In either case, we are looking forward to your reply!\n\n**Reference:**\nPapamarkou, Theodore, et al. \"Position: Bayesian Deep Learning is Needed in the Age of Large-Scale AI.\" Forty-first International Conference on Machine Learning. 2024.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "D9mSAgNDdk", "reviewer": "Reviewer_3YyQ", "summary": "The paper proposes an adjoint method for functions of matrices that utilize Arnoldi/Lanczos iterations to compute gradients with respect to large dimensional variables and demonstrates their approach's utility on a variety of common compute/memory intensive tasks encountered in machine learning. They derive the adjoints by implicitly differentiating the solution of Arnoldi/Lanczos iterations with respect to the input matrix and vector and then show how this can be used to compute gradients with respect to parameters that the input matrix depends on. They then examine their method on three case studies: Gaussian process hyperparameter optimization, solving for parameters of a discretized two-dimensional (in space) PDE and post-training calibration of a Bayesian neural network. Their results demonstrate an improvement over the standard approaches on a range of metrics.", "questions": "- The last part of the explanation of implicit differentiation in section 4 was confusing. On the one hand $\\mathrm{d}\\rho$ is defined in eqn (9) but the form $\\mathrm{d}\\rho = \\langle \\nabla_{\\theta} \\rho, \\mathrm{d}\\theta \\rangle$ is given on the following line. It was not clear as to why these two should be equal. Being more concrete on the definition of $\\rho$ would help here\n- As a solution to the second weakness listed above an algorithm comparing adjoint Arnoldi/Lanczos with reverse-mode AD may help to make clear the differences between the methods.\n- On the sentence starting with \"Gradients ...\" on line 198, doesn't reverse-mode AD also only use matrix-vector products as well?\n- A definition of \"loss of accuracy\" in Table 2 should probably be given in the main text.\n- A bit of discussion as to why the JAX low-memory matrix-vector product implementations cannot compete with KeOps would be a worthwhile addition to the paper.", "limitations": "Limitations on the conceptual aspect of the approach are discussed at the end of section 3 but some discussion in light of the experimental results would also be of value", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- An adjoint method for the Arnoldi iteration is certainly a useful addition to the automatic differentiation toolbox given the techniques wide use and would thus be of significant interest to the Neurips community.\n- The paper was mostly clear to read and understand. The conciseness of the experiment section is also a strength as the key aspects of the setups and the results are succinctly presented.\n- The paper's results indicate the methods utility in obtaining better runtimes, more accurate gradients and/or lower train and test losses on a range of important tasks in machine learning. In particular, the BNN calibration experiment clearly demonstrates how their method can handle matrices that may be too large to fit on a GPU and so it doesn't have to resort to subsequent approximation which they show leads to better performance.", "weaknesses": "- Although the paper empirically shows how backprop through a matrix function is slower than their adjoint method for Arnoldi iterations in Fig 3, the fact that a sparse matrix (as opposed to a dense one) was utilized motivates some skepticism over whether we would expect to see this same runtime relationship for a dense matrix. In fact, the exact GP results suggest that implementing matrix-matrix (/-vector) operations more efficiently can be enough to make reverse-mode AD work quicker for dense matrices. More discussion on why the adjoint Arnoldi iteration is more appropriate would strengthen their argument for the improved efficiency of their method.\n- The description of the adjoint systems in section 4.1 is lacking some explanation to make clear the key aspects of the adjoint methods and why this is preferable over reverse-mode AD. Reverse-mode AD also only utilizes vector-Jacobian products and so is a kind of \"matrix-free\" method. Explanation of how their method differs from this would make more apparent their contribution.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the positive evaluation! \n\nBefore we answer your questions, we would like to reply to your points listed under \"Weaknesses\" briefly:\n\n**Matvecs, sparse/dense complexity:**\nThe complexity of our adjoint mirrors that of the forward pass of Lanczos/Arnoldi and depends almost entirely on the efficiency of the matrix-vector product. Therefore, Lanczos and Arnoldi are mainly used on sparse or highly structured matrices like PDE discretisations or Jacobians of neural networks, which is why we use a sparse matrix in Figure 3. \nWith dense matrices, an efficient matrix-vector product can quickly become the limiting factor already for the forward pass (which is what happens in the Gaussian process case study; see the discussion in Appendix G).\nFinally, unless we misunderstand your sentence \n\n    In fact, the exact GP results suggest that implementing matrix-matrix (/-vector) operations more efficiently can be enough to make reverse-mode AD work quicker for dense matrices\n\nwe would like to reinforce that the Gaussian process case study does not compare to 'backpropagation through the solver' but an alternative custom gradient operation (alternative to our adjoint), as implemented in GPyTorch; see Equation 5 and/or the introduction of Section 5.\n\n\n\n**Reverse-mode AD:**\nWhen you refer to reverse-mode AD in 'why this (the adjoint method) is preferable over reverse-mode AD', do you mean what we call 'backpropagation through the solver'? The adjoint method is one of many ways of implementing reverse-mode AD. If you are looking for a comparison of how backpropagation through the solver relates to our gradient implementation: we show in Figure 3 that our code is far more efficient. See also Appendix A for additional information.\n\n\n**Answers to questions:**\n1. Thank you for bringing that up. The definition is that $\\mathrm{d}x$ is an infinitesimal perturbation. The gradient identity is additional information (and critical for the following derivations). We will update the sentence accordingly.\n2. We agree that this is crucial; see Figure 3.\n3. Thank you for bringing that up. If we back-propagate ''through'' Lanczos/Arnoldi, the resulting (automatic) gradients would use reverse-mode derivatives of matrix-vector products. Since our algorithm replaces this step, the sentence in line 198 emphasises that our code is matrix-free, just like the forward pass; see also Corollary 4.3.\nWe will revise the sentence to make this distinction more clear.\n4. Equation 54 in Appendix F formally defines the loss of accuracy. We will link this information more clearly in the main text.\n5. We agree! We dedicate Appendix G to discussing JAX vs KeOps.\n\n\nThank you again for your review!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "c0pcZv6Spu", "reviewer": "Reviewer_X6aJ", "summary": "There are some useful iterative methods for calculating important matrix products, to wit, Lanczos and Arnoldi iterations, which apparently did not have known derivatives, until now. The paper provides a generic framework for calculating the derivatives of such iterative linear operator approximations. The core of the paper is interpreting the Lanczos resp Arnoldi iteration as an iterative system and using the implicit derivative trick", "questions": "I have no questions. This paper was well explained, and clear about its goals, and how it approached them.", "limitations": "There are many limitations to the methods divulged here, but they seem to be adequately articulated in the paper.", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "While these results are not *world* shaking, they seem to provide immediate quality-of-life improvements to linear algebra users in high-value problems. The results unify several known \"matrix free\" tricks, which is satisfying, and point the way to more. \n\nThe paper is well-written, compact, and easy to understand, which is a great relief this deep into the review process.", "weaknesses": "As far as I can see, few. The paper seems to present what I need to know.\n\nIt would not be impossible for these results to be known elsewhere in the literature, but if they are known, I have not seen them.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate your positive assessment!\nWe are happy that you acknowledge how our work \"provide(s) immediate quality-of-life improvements to linear algebra users in high-value problems\" because this was precisely our goal!\nWe hope that you continue to fight for this paper's acceptance.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jZ7etiMnRW", "reviewer": "Reviewer_FpA9", "summary": "This paper introduces a new matrix free method for automatically differentiating functions of matrices. The computational problem discussed in this paper is of interest because the matrix dimension scales with respect to the size of the dataset (e.g. Gaussian process regression, etc.). The authors' algorithm yields the exact gradients of the forward pass, all gradients are obtained with the same code, and said code runs in linear time- and memory-complexity. The proposed method is also matrix-free and does not form the matrix explicitly which is the key to scalability. The authors' method evaluate functions of large matrices by: 1) first decomposing the matrix into a product of small matrices (e.g. Lanczos, Arnoldi iterations); 2) then evaluating functions on the small matrix produced during the Lanczos/Arnoldi iteration. The authors note that functions of small matrices can be differentiated efficiently. The authors provide the method to compute gradients of the Lanczos/Arnoldi iteration so that backpropagation can be implemented efficiently. The authors use differentiation via adjoints of the Lanczos/Anroldi iterations (this is described in Theorem 4.1 and Theorem 4.2). The parameter gradient expression that can be used for any problem is provided in Corollary 4.3. The authors provide experimental results which provide impressive speedup over GPyTorch over problems of interest including: Gaussian process, physics-informed PDEs, Bayesian neural network tuning.", "questions": "1. There are parallelization techniques that may be useful for even scaling up the proposed method.\n2. How many samples do you take for stochastic estimation of log-determinant for the exact Gaussian process result?\n3. What norm is the gradient error measured in?", "limitations": "The authors describe several limitations - but as the authors mention, these extensions describe a more expanded write-up.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. This paper tackles an important computational bottleneck in tuning hyperparameters of potentially complex models that are used in scientific and machine learning applications.\n2. The proposed method is general and can accomodate any matrix functionals and do not need explicit derivations of the gradients.\n3. The authors provide experimental results on several problems of interest to machine learning practitioners.", "weaknesses": "There are no particular weaknesses I found on reading the paper. The paper is well written.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your positive evaluation! \n\nTo answer your questions:\n\n1. We agree that these directions are interesting for future research!\n2. We use 10 Rademacher samples to match GPyTorch's default settings. Appendix H explains all Gaussian-process-related parameters. \n3. The errors are relative mean-squared errors (using the square root of the machine epsilon in the current floating-point accuracy as a nugget to avoid dividing by zero). We accidentally omitted this information in the submission and will add it to the next version.\n\nThanks again for your review!\nWe hope that you continue to fight for this paper's acceptance.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all reviewers for their reviews and for assessing the paper so positively!
We are grateful that all reviewers praised the clarity of the contribution.
Reviewers FpA9 and X6Aj did not find any particular weaknesses, and we believe that the weaknesses listed by RkAM and 3YyQ might be easy to resolve.

We replied to all reviews in separate threads and look forward to the discussion.

Thank you again, and best wishes,

Authors

Meta-review: This paper presents a novel method for efficiently computing gradients for functions of large matrices using Lanczos and Arnoldi iterations, tackling challenges in automatic differentiation for tasks such as Gaussian processes, PDEs, and Bayesian neural networks. Implemented in JAX, the approach offers notable improvements in speed and scalability without the need for explicit matrix formation or problem-specific optimizations. Reviewers commend the paper for its strong theoretical contribution, clear presentation, and empirical validation, with only minor criticism regarding the empirical evaluation. I concur with the reviewers' consensus on acceptance.
Decision: Accept (spotlight)