Sparser, Better, Faster, Stronger: Sparsity Detection for Efficient Automatic Differentiation
Authors: Adrian Hill
Guillaume Dalle
Venue: Accepted by TMLR ()

Abstract: From implicit differentiation to probabilistic modeling, Jacobian and Hessian matrices have many potential use cases in Machine Learning (ML), but they are viewed as computationally prohibitive. Fortunately, these matrices often exhibit sparsity, which can be leveraged to speed up the process of Automatic Differentiation (AD).
This paper presents advances in sparsity detection, previously the performance bottleneck of Automatic Sparse Differentiation (ASD). Our implementation of sparsity detection is based on operator overloading, able to detect both local and global sparsity patterns, and supports flexible index set representations. It is fully automatic and requires no modification of user code, making it compatible with existing ML codebases.
Most importantly, it is highly performant, unlocking Jacobians and Hessians at scales where they were considered too expensive to compute. On real-world problems from scientific ML, graph neural networks and optimization, we show significant speed-ups of up to three orders of magnitude. Notably, using our sparsity detection system, ASD outperforms standard AD for one-off computations, without amortization of either sparsity detection or matrix coloring.

Published in Transactions on Machine Learning Research (05/2025)
Sparser, Better, Faster, Stronger: Sparsity Detection for Efficient Automatic Differentiation
Adrian Hill"
hill@tt-berlin.de
BIFOLD -Brllin Institute for the Foundations of Learning and Data, Berlin, Germany Machine Learning Group, Technicol University of Berlin Berlin, Germany
Guillaume Dalle
guilauume.alleeenppy.
LVMT, ENPC, Institut Polgtechnique de Paris, Undv Custave Biffel, Marne-la- Morne-la-Vallée, France
Reviewed on OpenReview:
Abstract
From implicit differentiation to probabilistic modeling, Jacobian and Hessian matrices have many potential use cases in Machine Learning (MI), but they are viewed as computationally prohibitive. Fortunately, these matrices often exhibit sparsity, which can be leveraged to speed up the process of Automatic Differentiation (AD). This paper presents advances in sporsity detection, previously the performance bottleneck of Automatic Sparse Differentia- tion (ASD). Our implementation of sparsity detection is based operator overloading, able to detect both local and global spansity patterns, and supports flexible index set representa- lions. It is fully automatic and requires modification of user code, making it compatible with existing ML codebases. Most importantly, it is highly performant, unlocking Jacobians and Hessians at scales where they were considered too expensive to computc. On real-world problems from scientific ML, graph neural networks and optimization, we show significant speed-ups of up three orders of magnitude. Notably, using our sparsity detection system, ASD outperforms standard AD for one-off computations, without amortization of either sparsity detection matrix coloring.
1 Introduction
1.1 Motivation
Machine Learning (ML) has witnessed incredible progress in the last decade, a lot of which was driven by Automatic Differentiation (AD) (Griewank and Walther 2008 Baydin et al. 2018 Blondel and Roulet, 2024p. Thanks to AD, working out gradients by hand is no longer a requirement for training differentiable ble models. User-friendly software packages like TensorFlow (Abadi al. 2015), PyTorch (Paszke et al. 2019) and JAX (Bradbury et al. 2018) allow practitioners quickly experiment with different models and architectures, resting assured that gradients will be computed efficiently and correctly] without human intervention
However, while gradient-based optimization has become ubiquilous within ML, the practical use of Jacobians and Hessians remains scarce. Conventional wisdom tells us that for realistic applications, these matrices are too large to handle, since we cannot afford to store n² coefficients memory when the number of parameters n reaches millions. A common workaround is manipulate matrices in form of so-called lazy linear operators BBondel and Roulet 2024). which defined only by their action vectors.
'Both authors have contributed equally this work. While this generally holds true; AD has several pitfalls. Refer Hückelheim al. 2024 for taxonomy.

---

Published in Transactions on Machine Learning Research (05/2025)
(a) AD code transformation
(b) Standard AD Jacobian computation
VJP(x, ei)
f(x)
VJP(x,v)
VJP(x, e2)
VJP(x, e3)
VJP(x, e4) =
(c) ASD Jacobian computation
0
0
VJP(x, e1 + t.)=
f
0
0
VJP(x, +2
Sparsity detection
Coloring
Matrix-vector products
Decompression
Figure 1: Comparison of moverse-mode AD and ASD
(a) Given a function f, reverse-mode AD return a function (æ, v) - v Of())) computing vector-Jacobian products (VJPs). (b) AD computes Jacobians row-by-row by evaluating VJPs with all basis vectors. (c) ASD reduces the number of VJP evaluations by first detecting a sparsity pattern of non-zero values, then coloring orthogonal rows the pattern, and simultaneously evaluating VJPs orthogonal rows. The first step of sparsity detection the performance bottleneck of ASD and the focal point of this paper. The concepts shown in this figure directly translate to forward-mode, which computes Jacobians column-by-column instead of row-by-row.
Luckily, in numerous applications within MT., most notably in the sciences, Jacobians and Hessians exhibit sparsity characteristic that has remained largely ignored by current ML software. By automatically detecting and leveraging the sparsity of Jacobians and Hessians (see Figure 1), their automatic computation be sped up by orders of magnitude in high-dimensional settings. Furthermore, when materialize these matrices instead of representing them as lazy lincar operators, computational ad advantages such as factorizations and direct linear solves are suddenly unlocked.
1.2 Applications
We enumerate concrete scenarios where Jacobians or Hessians appear naturally.
Newton's method [Nocedal and Wright, 2006 Chapter 3) is a fast root-finding and optimization algorithm,
which is easier to implement with AD. To find zero of the vector-to-vector function f: R - Rm, Newton's method performs the following iteration:
a=1 1) = =x(1)= 0ʃ (a (1)) (()..
To minimize the vector-to-scalar function f: R" + R without constraints, which amounts to finding a zero of the gradient (=))), Newton's 's method turns into:
at + = xt)) - (/
In both cases, need to solve linear system of equations, involving either Jacobian matrix 0f(x) Hossian matrix 2 fxx) (which is the Jacobian of the gradient). Specifically for optimization, a lot of research effort went into quasi-Newlon methods and their limited-memory variants (Nocedal and Wright, 2006, Chapters 6 and 7), which leverage cheap approximations
the inverse Hessian. Still, evaluating the exact Hessian can prove beneficial, for instance study its spectral properties. In deep learning, the maximum eigenvalue of the training loss Hessian provides insights into the dynamics and stability of gradient descent. (Cohen et. al. 2020
2

---

Published in Transactions on Machine Learning Research (05/2025)
Implicit differentiation has become more prevalent in ML with the rise of implicit layers <Koltrr et al, 2020). When a vector-to-vector function f(x is defined implicitly by conditions of the form g( 0) O, the implicit function theorem lets us recover the Jacobian of f by solving yet another linear system, this time with partial Jacobians: 0f(0) = - 0)- (2g(f(0))
()0 0). For unconstrained optimization ((0) = arg miny cly I)) the optimality criterion is g( 0) = Vic(s((), 0) = 0, and so we obtain a partial Hessian to invert. The recent survey by [Blondel al. (2022) gives more insights and examples on implicit differentiation and its connection to AD.
Further applications to probabilistic modeling are discussed in Appendix A In all of the scenarios mentioned above, we observe that (1) Jacobians and Hessians are useful objects, (2) they are often large and computed with AD, and (3) exact computation is deemed intractable, which seemingly justifies approximations or lazy operators. When the matrices in question exhibit sparsity, we claim that the last item should be examined more closely, and possibly refuted.
1.3 Contributions
Despite it being a well-researched area (Griewank and Walther 2008), sparse differentiation is not widely used in the ML community. We identify two likely reasons for this situation. First, the lack of tooling for aulomatic sparsily delection, which forces potential users to work out the sparsity pattern by hand (a tedious and error-prone process). Second and most importantly, the separation between the AD and ML scientific commmmities (in terms of research groups, publication venues and programming languages), which means that AD advances do not percolate easily into ML circles. For instance, many of the existing sparse differentiation libraries are written in compiled languages like or Fortran. While powerful, these languages lack the flexibility and iteration speed required by modern ML workflows, which favor dynamic alternatives like Python, R and Julia.
Our contributions are meant to fill these gaps. On the theoretical side, we present operator overloading for sparsity detection in a new light, reformulating existing techniques from the AD literature as a binarization of the chain rule. The algorithms we describe have been known for at least two decades, but we hope our new presentation will be more natural for people unfamiliar with the AD literature. Our perspective abstracts away implementation details like computational graphs and the data structures used for bookkoeping. On the practical side, we introduce SparseConnect wtw.eeoooonntiiiitttaaeteee viityrracer jil (SCT), highly performant sparsily detection code written in the open-source Julia programming language (Bezanson et al., 2017p. SCT allows users to detect both local and global sparsity patterns, naturally handles dead ends which can occur in traditional graph-based approaches, allows for a flexible selection of data structures for sparsity pattern ions,, and enables novel tensor-level overloads that can drastically reduce the amount of operations required to compute sparsity patterns. Drawing on Julia's multiple dispatch paradigm allows it to generate highly performant machine code for each overloaded operator.
We also provide in-depth benchmarks of SCT part of a full Automatic Sparse Differentiation (ASD) system in Julia for the computation of Jacobians and Hessians. Thanks to sparsity detection, ASD aulomatically leverage complicated Jacobian and Hessian sparsity patterns without any human involvement. Our benchmarks demonstrate that this ASD implementation outperforms AD at scales previously considered impractical. Notably, it is performant enough to enable spoed-ups in onc-off computations of Jacobians and Hessians. Thus, while amortizing the computational cost of sparsity detection over several Jacobian and Hessian computations can help, is not always necessary.
The approach we present docs not require any rewriting of user functions, making it compatible with main- stream packages in the Julia ecosystem (such as those for deep learning and scientific ML). We hope that it will provide a blueprint for adaptation in other languages and frameworks, especially in Python which currently lacks sparsity detection and ASD) tooling
Reler Appendix B for an overview of ASD implementations other high-level programming languages.
3

---

Published in Transactions on Machine Learning Research (05/2025)
1.4 Notations
Scalar quantities are denoted by lowercase letters I, vectors by boldface lowercase letters 2, and matrices by boldface uppercase letters A. Given a vector æ; its coefficients are written Hú. Given a matrix A, its columns are written A.,, its rows are written A.; and its coefficients are written A.,- We use the word "tensor" when want to refer to either vectors matrices. The centered dot stands for product wetween two scalars, or the product between scalar and tensor. Integer ranges are denoted by (n....,m.. For a vector-to-scale function f: R" + R. we write ff()) E R" for its gradient vector and f²f(x) E R=x" for its Hessian matrix. For a pector-to-vector function f: IR" + R"" we write (f()) E Rmx= for its Jacobian matrix. The partial derivative (or partial Jacobian) of a function with respect to its k-th argument is denoted by 3u: while the total derivative with respect to a variable v is denoted by dg. Unless otherwise specified, all functions considered here are sufficiently differentiable the point of interest. We will also be interested in the sparsity patterns of gradients, Jacobians and Hessians. The "one" function is defined on numbers as ||x| = 1 if x # 0 and [ = 0 otherwise. Since the sparsity pattern of a generic tensor T is just the one function applied to every element, we write it as [T). We use V lo denote the binary OR operation V b. Note that multiplication has priority over the OR operation.
1.5 Outline
2 gives summary of AD and ASD techniques well the role of sparsity detection in ASD. contains updated formulation of sparsity detection. Section describes the software imple- mentation of our sparsity detection code SCT. Section 5 numerical experiments for Jacobian and Hessian sparsity detection, as well as resulting ASD benchmarks. Section 6 concludes with future rescarch perspectives.
Section Section 8
2
Background
2.1 Automatic differentiation
As highlighted in the survey by Baydin et al. (2018) AD is method for computing derivatives that is neither numeric (based on finite difference approximations) nor symbolic (based on algebraic manipulations of expressions). Instead, AD works with non-standard interpretation of the source code, allowing it to carry derivatives along with primal values. The reference textbook on the subject is the one by Griewank and Walther (2008), while a more recent treatment is given by Blondel and Roulet (2024). Here we briefly recap the complexities of the main AD modes: forward, reverse, and forward-over-revense.
Let us consider vector-to-vector function f: R. + IR and an input æ E R We also fix a perturbation along the input E I=" (tangent) and a perturbation along the output v E Rm (cotangent). We call 7 the unit time complexity of evaluating f(x) (which may scale with the input size n or output size m). Forward mode AD can compute (x, u) + Of()u, called the Jacobian- Vector Product (JVP), in time O(T). Symmetrically, reverse mode AD can compute (ax, v) + e (f)), called the Vector-Jacobian Product (VJP), in the same order of time ((r). In particular, the case m implies that gradients are cheap to compute in reverse mode, observed by Baur and Strassen 1983)..
Now let us consider a vector-to-scalat function f: R + R. still with input x and unit time complexity т. The gradient = fxx) can be computed by reverse mode AD. Forward mode AD can then be applied to the computation of g- Thus, forward-oven-reveesee mode AD can compute (ær, u) + f(x),, called the Viessian-Vector Hessian- Vector Product (HVP), in time (r)). This observation was first made by Pearlmutter (1994) and revisited by LeCun al. (2012) and Dagréou et al. (2024). Note that some AD systems support botched evaluation that is, the joint application of JVPs, VJPs or HVPs to a vector (or batch) of seeds (ui...., Up) all at once.
Thhi variant is also commonly called vector mode. Given that the seeds themselves can also be vectors, and that the word "mode" already refers to forward reverse, we hope that our choice of terminalogy will be less confusing.
4

---

Published in Transactions on Machine Learning Research (05/2025)
2.2 Lazy products are not always enough
The routines mentioned above compute matrix-vector products Au involving Jacobians Hessians, for a cost that is small multiple of the cost of f. But often enough, need more complex quantities like matrix-matrix products AU or solutions of linear systems A-1b. In such cases, solution based purely on matrix-vector products may be suboptimal, and having access to the full matrix A can yield accelerations.
A matrix-matrix product AU can be computed from n lazy matrix-voctor products AU possibly batchod. But given the materialized matrix A, more efficient procedures exist, for instance in implementations of BLAS Level III (Lawson et al. 1979 Blackford et al. 2002). Similarly, a linear system Au b can be solved using only matrix-voctor products if we resort iterative methods like the Conjugate Gradient (CG) (Hestenes) and Stiefel 1952) GMRES (Saad and Schultz 1986). The precision of these methods is tied to the number of ilerations, each of which costs around the same as one function call. On the other hand, given access to A, we can use direct factorization-bass solver such as those in LAPACK (Anderson et al. 1999).
When the materialized matrix A is encoded in a sparse format, like Compressed Sparse Column (CSC), Compressed Sparse Row (CSR) or COOrdinate (COO), these conclusions still hold. In particular, complex linear algebra operations can be executed even faster thanks to dedicated libraries. A prominent example is 4 Davis 2024), which includes sparse matrix factorizations and direct linear solvers. Since many high-dimensional Jacobians and Hessians are naturally sparse, can leverage this property to speed up computations if we are able to reconstruct them efficiently.
the SuiteSparse ecosystem
Numerical linear algebra is not only faster when working on materialized matrices (dense or sparse), it can also be more robust. Iterative solvers shine most when the matrix A is well-conditioned, or when a good proconditioner is known (Stewart, 2022). Outside of these conditions, the gain with respect to direct solvers is less obvious. Finally, it is worth noting that some nonlinear optimization libraries only accept Jacobian/ Ceoobinn/Hessian Hessian matrices (usually in sparse formats) for their linear system subroutines, and cannot work with lazy matrix-vector products. Examples include the popular Ipopt (Wächter and Biegler] 2006) and Knitro (Byrd et al. 2006). All of this suggests that it may be worth paying an initial fee to matcrialize the matrix, which we recoup as subsequent operations are sped up. This is supported by our experiments in Appendix G and Isubsection 5.2,
2.3 Reconstructing (sparse) matrices from products
One can reconstruct a dense matrix A by taking its products A., = Aeli) (resp. A.: = (T all the basis vectors of the input (resp. output) space. For a function f : RD - Imm the Jacobian is either built columm-by-column with n JVPs, or row-by-row with m VJPs, as shown in Figure lp. For the Hessian, both options are equivalent due to symmetry. With batched AD, several seeds can be provided at once for the products. In the extreme case where (1)) ...., (')) are all supplied together, batched AD amounts to a matrix-matrix product with the identity AI, The complexity of this operation scales with the number of basis vectors.
A) with
When the matrix A is known to be sparse, this painstaking reconstruction can be greatly accelerated. One option is to perform sparse batched AD (Griewank and Walther 2008 Chapter 7), essentially computing AI in one pass while dynamically exploiting sparsity inside the function. This approach only applies to a few AD systems because it requires sparsity-aware differentiation of each elementary operation, which is not always implemented. On the other hand, matrix-vector products are always available as the lowest-level primitive of any AD system. To leverage these products generically, thus focus on compressed evaluation of the matrix (Griewank and Walther 2008 Chapter 8), which is the standard method for ASD because it can be implemented on top of any eristing AD backend.
The core idea behind ASD is that, if columns A.j, and A. are structurally orthogonal (they do not share a non-zero coefficient), we can evaluate them together with a single matrix-vector product A(elh) + eli)) and then decompress the sum in a unique fashion. An illustration of this procedure is shown in Figure Finding large sets of uralyy independent columns or rows helps lower the number of products necessary
htsps://peope en....r.eee.e................. tamu. edu/ davi s/svi tas.pts..ppeeeeee.. -se.
5

---

Published in Transactions on Machine Learning Research (05/2025)
to recover A. As shown in the review by Gebremedhin et al. (2005), this matrix partitioning problem is equivalent to a graph coloring problem, where the graph is constructed based on the rows and columns of A. Overlapping columns or rows must get different colors, and the goal is to find an assignment which minimizes the total number of colors used. Denoting by Ca the number of colors in the columns and by cm the number of colors in the rows; we see that only en (resp. Cm) products are needed to build the matrix column-by-column (resp. row-by-row). For typical sparse matrices, Cn < n and C. & m, which makes ASD a huge improvement in complexity. For instance, a forward-mode sparse Jacobian can be computed with cost (c,,T) instead of O(Tr).
Crucially, ASD through compressod evaluation requires priori knowledge of the sparsity pattern (where the structural zeros are located). In some special cases, this pattern can be described manually: diagonal and banded matrices are common examples. However, more sophisticated sparsity patterns can emerge from complex code, which makes sparsity detection a key component of AST If the sparsity pattern does not depend on the input z, it can be reused across several AD calls at different points. The same goes for the result of coloring. Therefore, runtime measurements usually do not include the cost. of this "preparation" step, which is amortized in the long run.
Finally, note that ASD remains accurate even when the sparsity pattern is overestimated. If we predict that a coefficient can sometimes be non-zero, but it is in fact always zero, the result. will still be correct. Of course, if the coloring involves more colors than strictly necessary, this overestimation makes differentiation slower. Still, this tradeoff might be interesting to save time on sparsity detection.
2.4 Related work
The literature on sparse differentiation dates back 50 years. Curtis et al. 1974) first notice that, when computing sparse Jacobians, one can save time by evaluating fewer matrix-vector products. [Pwell and Toint (1979) exlend that insight to sparse Hessians. In the following years, the connection to graph coloring is discovered and several heuristic algorithms are proposed, see Gebremedhin et al. (2003) and references therein. While carly works expect the user to provide the sparsity pattern, automated sparsity detection quickly becomes a topic of research. Approaches to sparsity detection can be either dynamic (run-time) or static (compile-tine), like for AD itself. In way, sparsily delection is equivalent to boolean AD.
Dixon et al. 1990) propose an operator overloading method based on "doublets" and "triplets" that encode sparse gradients or Hessians. Griewank and Reese (1991) offer an alternative point of view by describing elimination of intermediate vertices in the linearized computational graph. Instead of derivative values, prop- agating only the sparsity patterns is often more efficient, given that binary information can be encoded into bit vectors (Geitner et al. 1995). Bischof et al. (1996) show that depending on the problem at hand, differ- ent sparse storage techniques may be preferred. Griewank and Mitev (2002) suggest a Bayesian criterion to select clever basis combinations and reduce the number of function calls even further. [iiering and Kaminski (2006) describe a static transformation of the source code, with rules that echo the aforementioned operator overloading, the bit vector encoding being present as well. Lohoff and Neftci (2024) propose an alternative approach to Jacobian accumulation that is based on reinforcement. learning and cross-country elimination GGriewank and Walther 2008 Chapter 9) and falls outside of forward and reverse mode. Thoir mothod leverages sparsity, but lacks support for dynamic control flow, needs to be retrained for each computat ional graph and is tailored to highly vectorizable functions. Specifically for Hessian sparsity patterns, Walther (2008) extends the operator overloading approach to recognize nonlinear interactions. Walther (2012) dis- cusses a faster variant of their initial algorithm, as well as the choice of the underlying sparse data structures. Meanwhile, Gower and Mello 2202)) introduce the adge-pushing algorithm which directly computes sparse Hessian with its values, bypassing the need for detection, coloring and compressed differentiation. While Walther 2008)) requires only forward propagation, Gower and Mello (0022) leverage a reverse pass to increase efficiency: a comparison can be found in Gower and Mello (2014). . The edge pushing algorithm is further improved by |Wann et al. (2016a); Petra et al. (2018), and shown to be equivalent the vertex elimination rule of Griewank and Reese (1911) by Wang et al. 2016b))
5 Even when sparsity patterns can be worked out by hand, helps to compare the result with automated sparsity detection.
6

---

Published in Transactions on Machine Learning Research (05/2025)
In terms of software, most existing sparse differentiation systems are implemented in low-level language like Fortran or C/C++. Prominent examples include ADIFOR (Bischof et al. 1992) and ADOL-C (Griewank
et al.,
1996 Walther 2009), along with the ColPack package for coloring (Gebremedhin et al., 2013
The closed-source MATLAB language also boasts a couple of implementations (Coleman and Verma 2000
Forth 2006; Weinstein and Rao 2017). Furthermore, several algebraic modeling languages for mathemalical programming and optimization include some form of sparse differentiation. is the case at least. for AMPL
(Fourer et a 1990), CasADi (Andersson et al., 2019) and JuMP, j1 (Dunning al. 2017), as well as the more recent and GPU-compatible ExaModels. ji (Shin et al.) 2024).
Nonetheless, lot of scientific and statistical code is developed directly in open-source, high-level languages like Python, R or Julia. Thus, there is a clear need for fully automatic sparse differentiation libraries which can differentiate user code without the translation layer of a modeling language. In Julia, the curent state of the art for sparsity detection relies on a package called Symbol cimbolliis.. ics. .jl for sparsity detection (Gowda ct al., 2019 2022). As ssection 5) demonstrates, our contributions inside SCT give rise to a much faster ASD pipeline. A survey of ASD implementation efforts in other high-level programming languages is given in appendix B
3 Detecting sparsity via operator overloading
We now present revised viewpoint on sparsity detection, exposing the principles behind first-order tracing. In contrast to standard literature (Griewank and Walther 2008 Walther 2008), we propose exposition that does not rely on the notion of computational graph and exploits local sparsity. It also provides a blueprint for easy implementation using a classification of operators, such as the one described in Section For second-order tracing, refer to Appendix D
3.1 Principle
Let f: R" + RT be a vector-to-vector function, and x E R". The Jacobian matrix 0f(z) E R= is obtained by stacking m gradient vectors, since its i-th row is the gradient (x(x)) E Rn of the scalar output component fz)). Thus, we can recover the Jacobian sparsity pattern 1[0f(x)] {0, mmmxa we know the gradient sparsity pattern [[V[((==]] (0, 1) of each component.
To achieve this, use a tracer, a number type which contains both primal value y(œ) E R (the actual number) and its gradient sparsity pattern 1[V)(=)] E {0, 1}" (a binary vector). Importantly, the gradient in question is taken with respect to the input vector æ. We initialize the procedure by turning every in the input into ez), where e, is the j-th basis vector. Using operator overloading, every intermediate scalar quantity involved in our function f replaced with a tracer. Thus, at the end of the computation, we recover a tracer (fi(z), ([Vf((=)) containing both the primal output and the desired gradient sparsity pattern. All have left. to do is write down the rules how two such numbers are combined, defining how gradient sparsity patterns propagate.
Remark 1. The tracer type is related to dual numbers, which are classic ingredient of forward-mode AD. More precisely, it encodes the sporsity pattern of a batched dual number, containing one directional derivative for each input z.j- Such batched dual numbers are a way to implement batched forward-mode AD (Revels
et al. 2016). One can also see it as the binary version of the sparse doublet in Dixon et al. 1990
3.2 Propagation rules
Let a(x) and ((x) be two intermediate scalar quantities in the computational graph of the function f(x). We compute new scalar 7(x) ((a(x), ((x) by applying two-argument operator io to a(x) and ((x). Our goal is to express the gradient sparsity pattern 1[V)(()] as a function of the intermediate sparsity patterns 1 [V)(x)) and [VV)(=)] By the chain rule, for any input index j E (1,....,n) the derivative with respect to input zj can be expressed as follows:
(x)(x)= day p(a(()) 8 ())))))) -)a(z)) @s)((a(z), 0,3(z)

---

Published in Transactions on Machine Learning Research (05/2025)
From now on, omit the dependence the input to lighten notations, but we keep in mind that everything is evaluated at point x:
0yr = de, 0 0 824 0,3
(1)
Bringing the indices together shows that Vy 01 Va + 024 VB. From there, the sparsity pattern emerges, if we extend V to represent the elementwise OR.
1 - 1 1 V 1[22]| 1
(2)
In other words, the propagation of gradient sparsity palterns through the operator only depends on two binary values: y] and [0]4]]. These binary values tell us whether the operator ip depends on each of its arguments at the first order.
Remark 2. Note the use of instead of in Equation 2 It is necessary because 1 [Va] and [VB forget about actual values, and thus remain blind to accidental cancellations. For instance, this method will always overestimale the sporsily pollern of se(ar, a- 3 whenever ax(x)) B(x). Fortumalely, as discussed above, these overestimates still give rise to correct derivatives inside ASD.
3.3 First-order operator classification
To implement Equation 2 we need to classify every elementary operator in our programming language depending on whether its partial derivative with respoct to each argument is zero. There are two ways to perform this classification: local (accurate) or global (conservative). Local classification takes inlo account the current value of a and B, while global classification considers every possible value. In the global case, we effectively replace 1[0 3)] with max.,9 a 100r 3)] in Equation 2
Tabbeeggive some examples. The max operator is especially interesting since it comes up in neural networks with activation function ReLU(x) max(x, 0). It is well-known that max (and therefore ReLU) induces local sparsity because it only depends on one of its two arguments: the first one whenever 2 B, and the second one whenever 8 2 a. Global sparsity will overlook this subtlety, because there exists a part of the space where a 2 3 and there exists a part of the space where 3 a, so that both arguments can influence the output at the first order. Local sparsity allows us figure out that only one of them will. As far we are aware, local sparsity has rarely been considered in previous works. An example of global and local sparsity patterns
convolutional layer can be found in Figures 5a and 5b in the appendix.
Global sparsity is still relevant, since it does not require propagating primal values through the computational graph, making it much cheaper to compute. Additionally, it yields a sparsity pattern that is valid over the entire input domain. The cost of the sparsity detection can therefore be amortized over the computation of multiple Jacobians and Hessians at different input points.
Local
Global
Operator p(a, B) exp, log sin, cos round, floor, ceil
1
122
1
1 1
1
1 a.c.
1
0 a.e.
0
1
1
1
max
a 2 8
8 2 a
1
1
min
a
$ 8
3 < a
Table 1: First-order classification of operators
Unary operators have no second argument. "a.e." means "almost everywhere" for the Lebesgue measure
4 Software implementation
We implement the tracer number types described in the previous section in an open-source software package called ObarrseConnectivii www.ttttetttyttttttttttoo (SCT). . For global sparsity detection, SCT implements two types
8

---

Published in Transactions on Machine Learning Research (05/2025)
of tracers: gradient tracers, which hold a gradient sparsity pattern 1 [[y(2)), and hessian tracers, which additionally hold a hessian sparsity pattern 1 ([))] As outlined in subsectiio 3.3 local sparsity detection requires an additional wrapper type that holds and propagates the value of the primal computation. The choice of data structures used to represent sparsity patterns [[)( and [[V) is flexible and can be altered to fit. the problem at hand, as we will discuss in subsection 4.1
To add an operator overload, a given operator must be classified according to Table SCT then automat- ically generates performant code that implements [Equation 2] as well as correctness tests that check the classification against the forward mode AD system ForwaraDiff... j1 (Revels et al. 2016 6 ated operator overloads work on generic Julia code that supports Real numbers. By making use of Julia's 's multiple dispatch paradigm, external software packages automatically call SCT's overloads when they are evaluated with our tracer types, requiring no code modification. As a result, it is already used by the SciML ocosystem (Julia's 's equivalent of SciPy): c-g. for nonlinear root-finding (Pal et al. 2024) and constrained optimization (Dxxit and Rackauckas] 2023 Furthermore, we make SCT's code generation utilities available for third party packages, enabling custom operator overloads.
The gener-
4.1 Sparsity pattern representations
For the purpose of mathematical exposition, it was convenient to define gradient and Hessian sparsity patterns as sparse binary tensors 1 [[V)(=)] E [0, and [[V]v(=)] E {0, 1|mxn. In SCT, multiple data structures can be flexibly used to represent these sparsity patterns. One approach is to represent sparsity patterns as the sets of (pairs of) indices of non-zero entries: set [Vy(ar] == (i (1..... .n} such that ay((x) # 0} and [I....,n]? such that (x) + 0}. Every operation on sparse binary tensors has
sset [V>y(z)] := {(i.j)
equivalent on index sets. The clementwise OR [[Va] V | is union set [Va] U Ulset [V Mcanwhile,
the outer product OR 1 [a] V [ is a Cartesian product set [Va] x set [V We can thus translate
Equation 2 as 1 [T] = [0] 1ett [a] L 1[ (ase..VV)
To implement it, various data structures for sets, such as hash tables or bit vectors, can be used. In the end, the right choice of set data structure will depend on the performance of unions and ileration, as remarked by Walther (2012). For small problems up to a few hundreds of inputs, bit vectors tend to be the most efficient, as each index requires only one bit of memory. Additionally, unions are fast, amounting to an OR operation on bits. The downside of bit vectors is that their memory requirement is constant, regardless of sparsity, and grows with the number of inputs. For sparse computations with very large inputs, hash tables that store indices as integers can therefore be more efficient. Going beyond sets, vectors of indices can be used that allow for duplicated storage. The upside is that unions are just a vector concatenation, which makes them very fast, but a final de-duplication step is always needed. The optimal sparsity pattern representation depends on the dimensions of the problem, the sparsity level, and more generally the structure of the computational graph. Since there is no universal right answer, our generie implementation of tracer types allows users to select the sparsity pattern representation that best suits their task. All sparsity pattern representations, including data structures for set types, can additionally be extended by users through multiple dispatch.
4.2 Tensor-level overloads
For the detection of global sparsity patterns, overloads are not exclusively implemented on a scalar level, but also on tensors of tracers. This allows us to bypass the original scalar computational graph for increased computational performance without any loss of precision. For example, when multiplying two matrices of tracers, instead of falling back to elementwise multiplication and addition using scalar overloads, we can make use of the fact that both multiplication and addition have an identical first-order classification (sce Table 1). This allows us to reorder operations, first computing elementwise OR operations along rows and columns of both matrices respectively. For two matrices of size (n x p) and (p x m), this reduces the amount +)+m)(pp m) 1) + n m, leading to a significant increase in performance. of CR operations from n-m (2p- - 1) to (n +
An in-depth derivation of this is given in appendix
TTe same holds for second-order tracing using https::/scift /
Table for classification, generating code according to Equation
9

---

Published in Transactions on Machine Learning Research (05/2025)
Similar overloads are implemented for common matrix operators like matrix norms and determinants. For functions that depend non-trivially on all inputs, like the determinant, a conservative sparsity pattern can be computed by returning the union of all input sparsity patterns, thus bypassing expensive operations.
4.3 Limitations
While our local tracer types fully support control flow statements, our global tracer types only support very limited subse [ due to their lack of a primal value. Common boolean functions (e.g. iszero, isfinitta also aren't supported. This is by design: if a global tracer were to enter branch in code, the returned pattern wouldn't 't be guaranteed to be correct, as a conservative pattern requires the evaluation of all branches. While overestimates of sparsity patterns are acceptable (as described in subsection 2.3) underestimates can lead to an erroneous coloring of orthogonal rows or columns and subsequently wrong ASD computations. Therefore, it is preferable to have global tracers throw an error rather than have them silently return incorrect sparsity patterns. Fortunatcly, this issuc is casily cireumvented by loveraging operator overloads. Building on the example from ssubsection 3.3, the function ReLU(x) max(x; 0) normally requires a comparison of primal values, branching to return either the tracer x 0. By classifying the first-order derivative of ReLU as not globally sparse, SCT is able to generate code for global tracers that circumvents this comparison, directly returning r which yields the most conservative pattern estimate.
The second limitation of SCT is its current lack of GPU support. As described in subsection 4.1 sparsity pattern representations of large inputs currently make use of hash tables or bit vectors, which dynamically allocate memory. While statically sized index sets should in theory be able to run on GPUs, this work is still pending. On top of this issue, GPU parallelized computations in deep learning typically make use of vectorized code, which often results in block-wise sparsity patterns. Such patterns can benefit from specialized implementation of sparsity detection that does not work on the scalar level. Still, scalar-level implementations such as ours are very relevant for scientific ML, which can involve a lot of fine control flow and individual indexing.
5 Numerical experiments
To benchmark SCT's performance in sparsity detection, we compare the previous ttaoe-of-tho-art in Julie") which relies on Symbolics j1. When measuring the benefits of SCT in the broader com- putation of Jacobians and Hessians, make use of two additional open-source Julia packages inside a single ASD pipeline. The first is aparsatatrixcolorinee
11 (Montoison et al.] 2025), which imple-
ments and improves coloring algorithms from Gebremedhin et al. (2055 2007 2009 2013). The second is Di ifffernttia:: on wownnterfacce.... j1 (Dalle and Hill 2025a a common interface for AD and ASD, which allows users to switch the sparsity detection method between SCT and Symbolica j1l in one line code.
A brief introduction to sparse matrix coloring can be found in Appendix C A demonstration of the API and its usage is given in Appendix F.
In some benchmarks, we distinguish between prepared and unprepared AD or ASD. Preparation refers to the parts of the pipeline that can be amortized across several computations. For AD it is often negligible, but for ASD preparation includes the costly sparsity detection and coloring, whose results can only be reused if they are input-agnostic. Prepared benchmarks are representative for applications which require the computation of multiple Jacobians of the same function, unprepared benchmarks correspond to one-off computations or cases where the sparsity pattern could conceivably change between runs.
All experiments and benchmarks were using Julia 1.11 on an Apple M3 Pro CPU with 36 GB of RAM. The plots are rendered using Makie.jl (Danisch and Krumbiegel 2021 While the bench- marks in the main body of this paper focus on sparse Jacobians, further experiments on sparse Ja- cobians with linear solves, sparse Hessians, and implicit differentiation on Graph Neural Networks can
be found in appendices G and 5.2 respectively. To ensure reproducibility, we provide the complete
Forr example, Julia's felse function supported, since SCT overloads it, regular if-else blocks aren't. SSince sparsity detection is language-speciffe, benchmarking across languages difficult.
10

---

Published in Transactions on Machine Learning Research (05/2025)
source code and matching Julia environments in public repository https:///gttuu..cc///aaaate sparser-better- saterrr---fattttttrsssssse
Hill 2025
5.1 Jacobian computation: Brusselator
5.1.1 Jacobian sparsity detection
To compare the performance of SCT with Symbolics, j1, we benchmark on the same example as Gowda et al. (211)), the Brusselator semilinear partial differential equation (PDE). which describes the spatial evolution of an autocatalytic chemical reaction (Prigogine and Lefever 1968 The PDE discretized to Nx x coupled ordinary differential equations using finite differences. Selected sparsity patterns are shown in appendix [..] The Brusselator benchmark is representative for the application of ASD to the field of scientific machine learning, where it can be used to accelerate the computation of Jacobians in hybrid Neural ODEs (Chen et al. 2018). Instead of training neural networks to learn the full dynamics of dynamical system from data, Rackauckas et al. (2021) study the incorporation of mechanistic priors, which in turn exhibit sparsity in the resulting Jacobian.
Table measures the wall time Jacobian sparsity detection depending on the dimensionality N of the discretized Brusselator PDE. SCT outperforms the state of the by one order of magnitude on large problems and two orders of magnitude small problems.
Problem
Sparsity
Sparsity detection*
Inputs
Outputs
Zeros
Colors?
Symbolics
SCT3
6
72
72
91.67%
9
5.07 10-3
2.10
10
(241.5)
12
288
288
97.92%
10
2.12 10-2
8.85
10-6
(240.0)
24
1152
1152
99.48%
10
7.48-10-2
(190.8)
48
4608
4608
99.87%
10
-
1.96 10-3
(157.2)
96
18432
18432
99.97%
10
1.45
10
1.71-107
(84.5)
192
73728
73728
99.99%
10
7.19 10
2.44 10-1
(29.5)
Wall time in seconds.
Number colors resulting from greedy column coloring.
In parentheses: Wall time ratio compared to Symbolics.jls jll's' sparsity detection (bigher is better).
Table 2: Performance comparison of Jacobian sparsity detection the Brusselator PDE.
5.1.2 Jacobian computation
We now benchmark the full computation of dense and sparse Jacobians on the same Brusselator PDE. For both AD and ASD, the forward-mode backend Forwardiiff.j] (Revels et al. 2016) is used to evaluato JVPs. The wall time breakdown of [Fgure 2 shows that sparsity detection previously was significant performance bottleneck. For small to medium problems, sparsity detection using Symbolics. jl takes up more wall time than the full AD Jacobian computation, which is no longer true for SCT.
The same measurements give rise lo the comparison in Figure 3 For large problems, prepared ASD using SCT accelerates the computation Jacobians by more than three orders magnitude compared to classical AD. More surprisingly, on all but the smallest Brusselator problem (N = -6), one-off unprepared ASD using SCT also outperforms prepared AD. The performance of SCT therefore opens performance gains via one-off ASD to more settings. Detailed benchmark timings given in appendix G.2
5.2 Implicit differentiation: GNNs
Bocause it gives access to materialized Jacobian matrix, ASD unlocks the use of direct linear solvers instead of iterative ones. Among other applications, this speed up differentiation of implicitlyddefine neural
11

---

Published in Transactions on Machine Learning Research (05/2025)
AD ASD (8ym..)
N=6
(SCT)
0.000
0.002
0.004
ASD (5ym.)
N=12
ASD (SCT)
0.000
0.005
0.010
0.015
0.020
AD
(8yma.)
ASD (SCT)
0.02
0.04
0.06
AD
ASD (5ym.)
N=48
0.1
0.2
ASD (8ym.))
N=96
ASD (8C/))
0.0
0.5
1.0
1.5
Walltime [s]
Sparsity detection
Memory allocation & coloring
Jacobian-vector products
Figure 2: Performance decomposition of unprepared AD and ASD Jacobian computations on small
discretizations of the Brussclator PDE. Using Symbolica bisb j11 (Sym.), sparsity detection is the performance bottleneck of ASD
10000
AD, prepared
Wall time [s]
ASD, prepared
10-2.5
ASD, unprepared using SCT
ASD, umprepared using Symbolics
10 5.0
12
24
48
96
192
N
Figure 3: Performance comparison of AD and ASD Jacobian computations on the Brussclator PDE.
network layers. We demonstrate it here with implicit Graph Neural Networks (GNNs), as introduced by Gu et al. (020)).
Consider graph A represented by its adjacency matrix, and matrix of vertex features U. To compute vertex embeddings and perform predictions, standard GNNs apply the following update fixed number of times:
x (++1) = ( ( (0) A (
Here, @ nonlinear activation, while ( and ( denote layer-specific weights. Conversely, an implicit GNN layer solves the fixed- toint equation
X = +WXA+
using an unbounded number of iterations. Instead of backpropagating gradients through these itcrations, which is very costly, Gu et al. (2020) suggest lo use implicit differentiation for computing θw X and On X. They derive an efficient formula to solve the linear system inside the implicit function theorem, but in general one could also use generic itcrative or direct solvers. For direct solvers, leveraging sparsity paramount even in moderate dimensions, because most graphs have sparse adjacency structures. We implement and compare both variants using the package ImplicitDi werteiieeifffeeetteeeoe ffeernffeeeteeeee iation. j]1 (Dalle) 2022). This library allows the
12

---

Published in Transactions on Machine Learning Research (05/2025)
VJP (iterative solve)
Training loss
0.50
Sparse Jacobian (direct solve)
0.25
0.00
0.0
0.5
0.84
Test accuracy
0.82
0.80
0.78
0.76
0.0
0.5 Wall time
1.0
Figure Training an implicit GNN with implicit differentiation, using either iterative solves (with AD) or direct linear solves (with ASD)
use of iterative solvers (based on JVPs or VJPs) as well as direct solvers (based on dense or sparse matrices) in the linear system of equations given by the implicit function theorem. The specific iterative solver used there is gmres from Krylov. (Montoison and Orban 2023 with its default parameters.
Our test case is simplified versions of the chains dataset (Gu et al. 2020 Appendix E.1). We perform node classification on a chain of = 10 nodes with 100-dimensional features, where only the first feature of the first node of the chain carries meaningful information (0 or 1 depending on the class to predict). The architecture is an implicit GNN with 16 hidden neurons and ReLU activation, followed by dropout and a linear prediction layer with cross-entropy loss. For the sake of simplicity, some implementation details differ from Gu et al. (2020). The most important ones are listed here, we refer the reader our experimental code for the rest. We added some small noise to the 99 non-informative features, whereas they were uniformly zero in the original paper. We replaced the orthogonal projection onto the & ball with simple scaling of the weights matrix W based on the spectral radius of A. This also guarantees that the fixed-point iteration won't t diverge. We build the training set using one positive and one negative chain, instead of drawing and masking nodes at random from several chains.
We allow ourselves these approximations because our goal is not to compare accuracy values with the initial experiment. Instead, we seek measure training times for two variants of implicit differentiation on a simple but representative toy example. One variant uses reverse-mode VJPss inside an iterative solve, with Zygote.jl [Innes] 2019) the backend, while the other computes forward-mode sparse Jacobians for direct solve, with ForwardDiff. j1 [Revels et al.] 2016) as the backend. The results are presented on Figure 4 where each dot corresponds lo one training epoch. As we can see, the orange and blue curves follow nearly identical trajectories, with the exact same sequence of jumps and drops. This proves that for our setting, the direct method yiolds the same training results as the iterative one, which is not surprising when the linear system non-degenerate. However, the orange curve is nearly two times faster in terms of time per epoch. Even though it starts with slight delay (due to the initial sparsity detection), the direct
13

---

Published in Transactions on Machine Learning Research (05/2025)
method based on sparse Jacobian matrices quickly catches up with and overtakes the iterative method based VJPs.
Thus, sparse Jacobians can yield a sizeable training speedup. Here, part of it may be due to the number of matrix-vector product evaluations, and part of to the increased efficiency of forward mode compared to reverse mode. The achievable speedup also depends on the parameters of the implicit linear solver, for instance the required precision or the maximum number of iterations.
In this specific instance, the sparsity pattern would have been straightforward to guess manually: it is block version of the chain graph's 's adjacency matrix. However, realistic GNN architectures are much more convoluted: they may involve various feature channels, skip connections, graph rewiring, etc. With ASD, these complications are no longer a hurdle to the user, who can freely compute sparse Jacobians whenever implicit differentiation benefits from them.
6
Conclusion
ASD is an essential part of the scientific computing toolkit. While its core ideas have been known for decades, its adoption in high-level ML frameworks is still lagging. We presented a refreshed formulation of sparsity detection using operator overloading, and described an ellicient software implementation which is already used at scale. Our hope is that such advances can spark renewed interest in sparse differentiation in the ML community and enable the practical use of Jacobian and Hessian matrices in domains where they were previously considered too expensive to compute.
Still, there are numerous research avenues to pursue. On the theoretical side, sparsity detection could be generalized to encompass various kinds of structures and symmetries, for instance block structure. This in turn could lead to efficient decomposition techniques for large-scale problems in variely of domains. On the practical side, our operator overloading implementation only supports a limited amount of control flow when computing global sparsity patterns, for which some amount of program transformation is needed. The packages we developed were designed for CPUs, but deep learning applications will require GPU support and iln-cation-free routines, leveraging hareware-specific primitives. Finally, we plan to explore interoperability or adaptation for the Python language, which the default choice in modern ML workflows.
Acknowledgments
We would like to thank Alexis Montoison for his valuable feedback and his work on AparseMatriscoloinnss
We also thank Robert Müller, Andreas Ziehe, Niklas Schmitz, Stefan Gugler and Marcel Langer for insightful discussions and feedback. Furthermore, we gratefully acknowledge funding from the German Federal Ministry of Education and Research under the grant BIFOLDD2B..
References
M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, Shlens,
B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals,
P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng. TensorFlow: Large-scale machine learning heterogeneous systems, 2015. URL https: //wp::///www..// /ww. www.srfwss....oo.. www..enssrflow.orgg///
E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J.
Du Croz, A. Greenbaum, S. Ham- marling, A. McKenney, and D. Sorensen. LAPACK Users' Guide. Society for Industrial and Applied Mathematics, Philadclphia, PA, 3 edition, 1999. ISBN 0-89871-447-8 (paperback).
J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl. CasADi: software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1-36, doi: (010077/112332011-000
Mar. 2019. ISSN 1867-2957. s 512532-0188 112532-018-0139- 0139-41
URL httpt:::/ /ttts:://oi..oo 0ttp://0010000010000000077.7
14

---

Published in Transactions on Machine Learning Research (05/2025)
S. Babaeinejadsarookoleee, A. Birchfield, R. D. Christie, C. Coffrin, C. DeMarco, R. Diao, M. Ferris, S. Flis- counakis, S. Greene, R. Huang, C. Josz, R. Korab, B. Lesieutre, J. Maeght, T. W. K. Mak, D. K. Molzahn, T. J. Overbye, P. Panciatici, B. Park, J. Snodgrass, A. Tbaileh, P. V. Hentenryck, and R. Zim- merman. The Power Grid Library for Benchmarking AC Optimal Power Flow Algorithms, Jan. 2021.
URL
abs/1 908. 02788
W. Baur and V. Strassen. The complexity of partial derivatives. Theoretical Computer Science, 22 13):317-300, Feb. 1983. ISSN 0304-3975. doi: 10.1016 10.1016/00344339555/33 (00 10-X. URL httpp:://wwt sciencedi
A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind. Automatic Differentiation in Machine Learning: A Survey, Journal of Machine Learning Research, 18(153):1 43, 2018. ISSN 1533-7928, URL html
http:// /mmrr.mlr...///... wttp//pppss//////1--------
M. Betancourt. A Conceptual Introduction to Hamiltonian Monte Carlo, July 2018. URL http://arxiv. org/ htt//bbb1110.2233
J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. Julia: A Fresh Approach to Numerical Computing. doi: 10.1137/1410006711 URL https:
SIAM Review, 59(1):65 98, Jan. 2017. ISSN 0036- 1445, 1095-7200.
77 //epubs. siam. wtt/ 1ttp:///00011.113//111100000 1137/141000677
C. Bischof, A. Carle, G. Corliss, A. Griewank, and P. Hovland. ADIFOR Generating Derivative Codes from Fortran Programs. Scientific Programming, 111):717832, 1992. ISSN 1875-919X. doi: 10. 10.1155/1992/ 717832. URL pettps://onlinelliee
wi wttpc//cc/////iceee /ttmm///../ab//1..11////9///////
C. H. Bischof, P. M. Khademi, A. Buaricha, and C. Alan. Efficient computation of gradients and Jacobians by dynamic exploitation of sparsity in automatic differentiation. Oplimiization Methods and Sofhware, 7 (1):1-30 39, Jan. 1996. ISSN 1055-6788. doi: 10. 108/105557890805555
URL
L. S. Blackford, A. Petitet, R. Pozo, K. Remington, R. C. Whaley, J. Demmel, J. Dongarra, 1. Duff, S. Hammarling, and G. Henry. An updated set of basic linear algebra subprograms (BLAS). ACM Transactions on Mathematical Software, 28(2):135 151, June 2002. ISSN 0098-3500. doi: 10.1145/567800 567807.
M. Blondel and V. Roulet. The Flements of Differentiable Programming, July 2024. URL http://arxiv.. org/ abs/ /400. 14606
M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-Lopez, F. Pedregosa, and J.-P. Vert. Efffcientt and Modular Implicit Differentiation. In Advances in Neural Information Processing Systems, volume 35, pages 3230 5242, Dec. 2022. URL
J. Bradbury, R. Frostig, P. Hawkins, M. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. Van- derPlas, S. Wanderman-Milne, and Q Zhang. JAX: Composable transformations of Python+NumPy programs, 2018. URL http:://ttpttee
M. Braun. sparseHessian An Package for Estimating Sparse Hessian Matrices. Journal of Statistical Software, 82(10), 2017. ISSN 1548-7660. doi: 00.18637/jss.v082.. URL httt:://www http///www.sttttsftt.. org/ arg/v82/110
R. H. Byrd, J. Nocedal, and R. A. Waltz. Knitro: An Integrated Package for Nonlinear Optimization. In G. Di Pillo and M. Roma, editors, Large-Scale Nontinear Optimization, pages 35-59. Springer US Boston, MA, 2006. ISBN 978-0-387-30065-8 doi: 10.1007/0-387-30065 URL htttss://aoiorgg/
M. B. Cain, R. P. O'Neill, and A. Castillo. History of Optimal Power Flow and Formula- tions. Technical report, 2012. URL httpp:://ww.
acopf-1-hi stoor--form lat inntoon-ee instoonntttti.i..pd
15

---

Published in Transactions on Machine Learning Research (05/2025)
J. Chen, H. Stoppels, H. Ranocha, and J. López. IterativeSolvers.jl: Iterative algorithms for solving linear systems, eigensystems, and singular value problems. JuliaLinean LinearAlgebra, 2013. URL httpp:://github. com/ com/ http://llauierrreeerrseeeetetees
R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud. Neural Ordinary Dif- ferential Equations. In Advances Newral Information Processing Systems, volume 31. Cur- ran Associates, Inc., 2018. URL lattps:/pproeedinns neurips. ceuuipps.cc/ppppe///00//// 35388ffbblfaaae42222443333333333 Abb.:s33t999-bbbbtatt......... htmll
C. Coffrin and O. Dowson. Rosetta-opf: AC-OPF Implementations in Varions NLP Modeling Frameworks. advanced network science initiative, 2022.
J. Cohen, S. Kaur, Y. Li, J. Z. Kolter, and A. Talwalkar. Gradient Descent on Neural Networks Typically Occurs at the Edge of Stability. In International Conference on Learning Represenlations, Oct. 2020. URL httt::://openeveee
pet/forum?id- -
rTtvkGeM
T. F. Coleman and A. Verma. ADMAT: the SIAM Workshop on Object Oriented Methods for Inler-Operuie Scienlific and Engineering Compuling, SIAM, Philadelphio, PA, volume 2, 1998.
An automatic diffcrentiation toolbox for MATLAB. In Proceedings of
T. F. Coleman and A. Verma. ADMIT-1: Automatic differentiation and MATL.AB interface toolbox. ACM Trans. Math. Softw., 26(1):150 175, Mar. 2000. ISSN 0098-3500. doi: 10.1145/3447337.447779. URL
attps://al. acm. ort 11e.11...11181881111 /ttp///.011.111151/7888337777771
A. R. Curtis, M. J. D. Powell, and J. K. Reid. On the Estimation of Sparse Jacobian Matrices. IMA Journal of Applied Mathematics, 13(1):117-113, Feb. 1974. ISSN 0272-4960. doi: 00009//mmma///// URL
sttps://dans org/ 10. 1091/ 11..maaat11.........
M. Dagréou, P. Ablin, S. Vaiter, and T. Morean. How to compute Hessian-vector products? In The Third Blogpost Track at ICLR 2024, Feb. 2024.
G. Dalle. Machine Learning and Combinatorial Optimization Algorithms, with Applications to Railway Planning. PhD thesis, Ecole des Ponts ParisTech, Dec. 2022. URL
tel 040-3403332
G. Dallc and A. Hill. A Common Interface for Automatic Diffcrentiation, May 2025a. URL htttt://p/arii org aatlaab//22550 2505. 05542]
G. Dalle and A. Hill.
Differeniiationnnterfaee... Zenodo, Mar. 2025b. URL https://zenddo.org records/150613
G. Dalle and A. Montoison. SparseMaatrixCobringggj... Zenodo, Mar. 2025. URL https://zenodb.org
Tecords/14984422
S. Danisch and J. Krumbiegel. Makie.jl: Flexible high-performance data visualization for Julia. Journal of
Open Source Software, 6(65) Sept. 2021. ISSN 2475-9066. doi: 10.21105//jos.033449 URL ht https:
[hooss.hheh..rrrr/appe//////11.. 03349
T. Davis. hrt.moocaAA..eeeaaaa/aaeeeeepee.
May 2024.
URL https://github.co com/
V. K. Dixit and C. Opliniiaiin.. A Unified Optimization Package. Zenodo, Mar. 2023. URL attps://zenode: htt.:/dd..ogg/eecoo/////33//
L. C. W. Dixon, Z. Maany, and M. Mohseninia. Automatic differentiation of large sparse systems. Journal of Economic Dynamics and Control, 14(2):299-311 May 1990. ISSN 0165-1889. doi: 10. 1016 1016/0165- 1889 URL https: http://wwwwweeeneeeeeetccctecc wtt.e//wweeeeeeeeoooooe/eeeeraaa 01651889909002330
16

---

Published in Transactions on Machine Learning Research (05/2025)
I. Dunning, J. Huchette, and M. Lubin. JuMP: A Modeling Language for Mathematical Optimization. SIAM Revieu, 59(2):295 320, Jan. 2017. ISSN 0036-1445. doi: 10/gftshn. URL pttthttss://epuua s amm.rog doi/
J. A. Fike and J. J. Alonso. Automatic Diffcrentiation Through the Use of Hyper-Dual Numbers for Second Derivatives, In S. Forth, P. Hovland, E. Phipps, J. Utke, and A. Walther, editors, Recent Advences in Algorithmic Differentiation, pages 163-173, Berlin, Heidelberg, 2012. Springer. ISBN 978-3-642-30023-3. doi: 10.1007/978-3-6422 3-2-23333-33333 15.
R. A. Fisher. Theory of Stalistical Estimation.
Mathematical Proceedings of the Cam- ISSN 1469-8064, 0305-0041. doi:
bridge Philosophical Society, 22(5):700-77 725, July 1925. mat
10110779003050041000008888
URL httpss//www..cnbbbbrre...o./////
abs http://wweetttaaaaaaaaaaa///aaaa
S. A. Forth. An efficient overloaded implementation of forward mode automatic differentiation in MATLAB. ACM Transactions Mathematical Software, 32(2):195-222, June 2006. ISSN 0098-3500. doi: 10.1145/
14188551141888 URL TRLhhtpps: httph:tt///// /cmm:///ll..c..mm.....mm org. doi /ttp://011111111111111161113 11411
R. Fourer, D. M. Gay, and B. W. Kernighan. A Modeling Language for Mathematical Programming. Managementt Science, 36(5):519 554, May 1990. ISSN 0025-1909. doi: 10.1287 (..1287/mnsc.365.555 mmnsc.36.5.51.
URL
attps://pubson
informs.... org/ /ttp//ol//bb//1...227//moooo 36.5.
A. H. Gebremedhin, Manne, and A. Pothen. What Color Is Your Jacobian? Graph Coloring for Computing Derivatives. SIAM Review, 47(4):629 705, Jan. 2005. ISSN 0036-1445. doi: 10/cmwds4. URL https:
A. H. Gebremedhin, A. Tarafdar, F. Manne, and A. Pothen. New Acyclic and Star Coloring Algorithms with Application to Computing Hessians. SLAM Journal on Scientific Computing, 29(3):1042 1072, Jan.
2007. ISSN 1064-8275. doi: 10.1137/050639977 URL betps://epuue siam. httpa//...aaa..oggoaab//////////
I137/050639879
A. H. Gebremedhin, A. Tarafdar, A. Pothen, and A. Walther. Efficient Computation of Sparse Hossians Using Coloring and Automatic Differentiation. INFORMS Journal on Computing, 21(2):209 223, May 2009. ISSN 1091-9856. doi: 10.1287 10.1287/jjoc.00000 1080..000. doi /t/aab//1/ 10. 287/1 1080. 0286
URL http:://pubsnnine. informs org/
A. H. Gebremedhin, D. Nguyen, M. M. Patwary, and A. Pothen. ColPack: Software for graph coloring and related problems in scientific computing. ACM Transactions on Mathematicul Software, 40(1):1:1 1:31, Oct. 2013. ISSN 0098-3500. doi:
URL https://dl org/ doi/
10.145/2533005255117
C. Geitner, J. Utke, and A. Griewank. Automatic Computation of Sparse Jacobians by Applying the Method of Newsam and Ramsdell. 1995. URL. ftstpa:/wwir cem. org/
Comput of Sparse- by the Geitnef---itttteeeeeeee le22:388ff133332222222223ddddddd
R. Giering and T. Kaminski. Automatic Sparsity Detection Implemented as Source-to-Son Source Transfor- mation. V. N. Alexandrov, G. D. van Albada, M. A. Sloot, and J. Dongarra, editors, Computational Science - ICCS 2006, pages 591-598, Berlin, Hoidelberg, 2006. Springer. ISBN 978-3-540-34336-8. doi: 10.1007/11758549 81.
M. Girolami and B. Calderhead. Riemann Manifold Langevin and Hamiltonian Monte Carlo Methods, Jour- nal of the Royal Statistical Society Series B: Statistical Methodology, 73(2): Mar. 2011. ISSN 1369- 7412. doi: 10.1111/j.1467-9866 000.098888..00.00666.... URL https::/doi.org// http:///oo.oorr/////// htt.://////....1.../.........../ 1411111/..117988888
2010 00765.x
17

---

Published in Transactions on Machine Learning Research (05/2025)
S. Gowda, Y. Ma, V. Churavy, A. Edelman, and C. Rackauckas. Sparsity Programming: Automated Sparsity- Aware Optimizations in Differentiable Programming. In Program Transformations for ML. Workshop al NewrIPS 2019, Sept. 2019. URL http::////pppeeeeeeeeee/ee/rr/rr
S. Gowda, Y. Ma, A. Cheli, M. Gwóázdż, V. B. Shah, A. Edelman, and C. Rackauckas. High-performance symbolic-mumeria via multiple dispatch. ACM Commun, Comput. Algebra, 55(3):92 96, Jan. 2022, ISSN doi: 10.114/3551528.355
1932-2232.
URL https://dl. acm org/
B511535
R. M. Gower and M. P. Mello. A new framework for the computation of Hessians. Oplimization Methods and Software, (7(2):251122 Apr. 2012. ISSN 1055-6788. doi: 10.1080/ 1111000/106667882011155000 URL
nttps://doi... org. 10. 10001105567680088 201.1821111100099911
R. M. Gower and M. P. Mello. Computing the sparsity pattern of Hossians using automatic differentiation.
ACM Transactions on Mathematical Software, 40(2):10:1 10:15, Mar. 2014. ISSN 0098-3500. doi: 10.
11 1145/2490254. URL https://d/ http::s/ddlaa....o.m/mmmmm
acm. www..aammmmmo........ htogg/doi/100 11111/..111//299225555
Griewank and C. Mitev. Detecting Jacobian sparsity patterns by Bayesian probing. Mathematical Programming, 13-)) June 2002. ISSN 1436-4646. doi: 10.1007/111107110051.. URL https:
//doi.org/100 http:////...10//1111100100000000
A. Griewank and S. Reese. On the calculation of Jacobian matrices by the Markowitz rule. Technical Report ANI,/CP-75176; CONF-910180-4, Argomne National Lab., IL (United States), Dec. 1991. URL
nttps://www. osti http///b/oo////11100/0 https//www..ot...oo/bbbloool1113
A. Griewank and A. Walther. Evaluating Derivatives: Principles and Techniques of Algorithmic Differen- tiation. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2nd ed edition, 2008 ISBN 978-0-89871-655-7 URL httth::///ppbbba https//s.....orr/oo/o/////////10 /tt.://...1117117787887777777779
9780898717761
A. Griewank, D. Juedes, and J. Utke. Algorithm 755: ADOL-C: A package for the automatic differentiation of algorithms written C/C++. ACM Transactions on Mathematical Software, 22(2):131 167, June 1996. ISSN 0098-3500. doi: 00.114/229443.2252 229474
URL https: //t:://ll.. acm. wttp:://///......//.//..
F. Gu, H. Chang, W. Zhu, S. Sojoudi, and L. El Ghaoui. Implicit graph neural networks. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Newral Information Processing Systems, volume 33, pages 11984-11999 Curran Associates, Inc., 2020. URL https: httttpp:://prro htttps:://proceedingg. Patee-22a3---PPP-----------
M. Hestenes and E. Sticfel. Mothods of conjugate gradients for solving lincar systems. Journal of Research of the National Bureau of Standards, 49(6):409, Dec. 1952. ISSN 0091-0635. doi: 10.6028/jree 00tp08/jrrr.49.00000 049.04.. URL htttps:///vlppbbsb http:///wss//.ppoo.....a..//ssss htt.s//ppppppppss//////////////p
pdff
A. Hill. htchll//parreebbetee--ffsst-tsre zenodo. 15473177
v0.1.0, May 2025. URL https://dbi.ozg/10
J. Hückelheim, H. Menon, W. Moses, B. Christianson, P. Hovland, and L. Hascoët. A taxonomy of automatic
differentiation pitfalls. WIREs Data Mining and Knowledge Discovery, 14(6)=e1555, 2024. ISSN 1942-4795. doi: 10.1002/widm. 1555. URL https://onlinet httpa//aww..llll....ooo///////oo http:///00/aa/1..100//...o.
555
M. Innes. Flux: Elegant machine learning with Julia. Journal of Open Sowrce Softturr,, 6625 May 2018. ISSN 2475-9066.
doi: 10.21105//jss.000022 URL attp::/jjass http:///aa...eeo...cccc//pppp//p
joss. 00602
M. Innes Don't Unroll Adjoint: Differentiating SSA-Form Programs, Mar. 2019. URL http://arxiy org/abs/1e 08..9991111111771 07951
18

---

Published in Transactions on Machine Learning Research (05/2025)
M. Innes, E. Saba, K. Fischer, D. Gandhi, M. C. Rudilosso, N. M. Joy, T. Karmali, A. Pal, and V. Shah. Fashionable Modelling with Flux, Nov. 2018. URL http:///rxiv..vrt wtt.o///////////oo////////////./ .01477
M. Innes, A. Edelman, K. Fischer, C. Rackauckas, E. Saba, V. B. Shah, and W. Tebbutt. A Differentiable Programming System to Bridge Machine Learning and Scientific: Computing, July 2019. URL http:
//arxii..rrg/ab//19900 07587
JuliaDiff contributors.
Oct. 2024. URL htl.at/gs//ggttub...o////uaaaa/a
Sparsebiffrools.1
R. V. Kharche and S. A. Forth. Source Transformation for MATLAB Automatic Differentiation. In V. N. Alexandrov, G. D. van Albada, P. M. A. Sloot, and J. Dongarra, editors, Compulational Science ICCS 2006, pages 558 565, Berlin, Heidelberg, 2006. Springer. ISBN 978-3-540-34356-8, doi: 10.1007/11758545 77.
J. Z. Kolter, D. Duvenaud, and M. Johnson. Deep Implicit Layers - - ODEs, Deep Equilibirum Models, and Beyond, 2020. URL http:://impp//eeee layers tuutoriallooro
K. Kristensen, A. Nielsen, C. W. Berg, H. Skaug, and B. M. Bell. TMB: Automatic Differentiation and Laplace Approximation. Journal of Statistical Software, 70:1-21, Apr. 2016. ISSN 1548-7660. doi: 10. 18637/jss.v070..m URL https:://dot org/10...0.0 /637/$ss. (0063///...0000.00.
v070.
C. L. Lawson, R. J. Hanson, D. R. Kincaid, and F. T. Krogh. Basic Linear Algebra Subprograms for Fortran Usage. ACM Transactions Mathematical Software, 5(3):308 323, Sept. 1979. ISSN 0098-3500. doi: 00.114//355841.35
URL https:/ /dl http://lllccccc...cc
1ttpe//111111113555551 355844
Y. A. LeCun, L. Bottou, G. B. Orr, and K.-R. Müller. Efficient BackProp. In G. Montavon, G. B. Orr, and K.-R. Müller, editors, Neural Networks: Tricks of the Trade: Second Edition, pages 9-48. Springer, Berlin, Heidelberg, 2012. ISBN 978-3-642-3528898 doi: 10.1007/978-3-1 -52299--99-8 URL https:
J. Lohoff and E. Neftei. Optimizing Automatic Differentiation with Deep Reinforcement Learning. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, Nov. 2024. URL https:
//cpenreview net/forum?iddnnne 90aa:
A. Montoison and D. Orban. Krylov.jl: A Julia basket of hand-picked Krylov methods. Journal of Open Source Software, 8189):5187, Sept. 2023. ISSN 2475-9066. doi: 10.21105/foss.05187. URL https://
Joss. hoss.thaac...rr/ppppe:///...1110
Joss.05187
A. Montoison, G. Dalle, and A. Gebremedhin. Revisiting Sparse Matrix Coloring and Bicoloring, May 2025. URL
07308
W. Moses and V. Churavy, Instead of Rewriting Foreign Code for Machine Learning, Automatically Synthe- size Fast Gradients. In Advances in Neural Information Processing Systems, volume 33, pages 12472- 12485. Curran Associates, Inc., 2020. URL htttps://proceedd neur ps. cc/paper//020/
Lec///93222222444422551118222333 Absttatttttaabttatat.t... Abstract html
W. S. Moses. Automated Derivative Sparsity via Dead Code Elimination, 2023. URL https://c.
wsmoses. cow.:ppresentatiiosss/eeeooolll
W. S. Moses, V. Churavy, L. Paehler, J. Hückelheim, S. H. K. Narayanan, M. Schanen, and J. Doerfert.
Reverso-mode automatic differentiation and optimization of GPU kernels via Enzyme. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC '21, pages 1-16, New York, NY, USA, Nov. 2021. Association for Computing Machinery. ISBN 978-1- 8442-1. doi: 10.1455/4488817.3776666
URL
P. Nobel. Auto diff: An automatic differentiation package for Python. In Proceedings of the 2020 Spring Simulation Conference, SpringSim '20, pages 1-12, San Diego, CA, USA, May 2020. Society for Computer Simulation International. ISBN 978-1-7138-1288--3
19

---

Published in Transactions on Machine Learning Research (05/2025)
J. Nocedal and S. Wright. Numerical Optimiaation. Springer, New York, NY, 2nd edition edition, July 2006. ISBN 978-0-387-30303-11
A. Pal, F. Holtorf, A. Larsson, T. Loman, Utkarsh, F. Schäefer, Q. Qu, A. Edelman, and C. Rackauckas. NonlinearSolve..t High-
Mar. 2024. URL Perfommeroom and Robust Solvers for Systems of Nonlinear Equations in Julia,
org/ abs/abs20000 241121433..1111111 16341]
A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, 7. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala. PyTorch: An Imperative Style, High-Performance Deep Learning Library. In Advances in Neural Information Processing Systems, volume 32. Cur- ngs wwwriee..rccocoo.o hwwprrppp..ccppppe///1///////// ran Associates, Inc., 2019. URL Abstract html
B. A. Pearlmutter. Fast Exact Multiplication by the Hessian. Neural Computation, 6(1):147 160, Jan. 1994. ISSN 0899-7667. doi: 10.1162/nncco.1944...11 URL https:://eeexplooot Leee org/ wtt..//.ee...bt.ttttce..// abst document/67955
C. G. Petra, F. Qiang, M. Lubin, and J. Huchette. On efficient Hessian computation using the edge pushing al- gorithm in Julia. Optimization Methods and Softwwree, 33(4-))11010-0102 Nov. 2018. ISSN 1055-6788. doi: 10.1080/10556785 781.7880001...1466.. URL /ttp:///pp///....o//00000000///0
2018. 200.1668888.1111186662332 1480625
L. Piloto, S. Liguori, S. Madjiheurem, M. Zgubic, S. Lovett, H. Tomlinson, S. Elster, C. Apps, and S. With- crspoon. CANOS: A Fast and Scalable Neural AC-OPF Solver Robust To N-1 Perturbations, Mar. 2024. URL http://
org/ abs/ 2403. 146661
M. J. D. Powell and Ph. L. Toint. On the Estimation of Sparse Hessian Matrices. SIAM Journal on Numerical Analysis, 16(6):1000-1077 Dec. 1979. ISSN 0036-1429. doi: 10.1137/071600 URL https:
//epubs sit.ssiiamm. http:///oo/aaa//0.111177170
I. Prigogine and R. Lefever Symmetry Breaking Instabilities in Dissipative Systems. II. The Journal of Chemical Physics, 4894):1695-17700, Feb. 1968. ISSN 0021-9606. doi: 10.1063/1.1608890, URL https:
C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and
A. Edelman. Universal Differential Equations for Scientific Machine Learning, Nov. 2021. URL http:
//arxiv.org./ wtt.://rrgbabb//0000 04385]
J. Revels. Diff:seDiff.jl: Reverse Mode Automatic Diflerentiation for Julia. JuliaDiff, 2016. URL https: //github.com///n http:///eaaa.lfffeeeeeeereeeee
J. Revels, M. Lubin, and T. Papamarkou. Forward-Mode Automatic Differentiation in Julia, July 2016. URL http: http://aaxii.v// org/ abs///bb.///111/11
C. Robert and G. Casella, Monte Carlo Statistical Methods. Springer New York, Aug. 2005. ISBN 978-0- 387-21239-5.
Y. Saad and M. H. Schultz. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7(3):866 869, July 1986. ISSN 0196-5204. doi: 10.1137 001137/0900058. URL http:://ppp..ps..a..o.oooooaaaa/
M. Schubert. Mffcchueert/pparssgeee Sept. 2024. URL. http://tssggggttacccc/////dtttee
S. Shin, F. Pacaud, and M. Anitescu. Accelerating Optimal Power Flow with GPUs; STMD Abstraction of Nonlinear Programs and Condensed-Spacee Interior-Point Methods, Feb. 2024. URL http://arxiv.
org/abs/2301 16830
D. Simpson. Un garçon pas comme lcs autres (Bayes) An unexpected dctour into partially symbolic, sparsity-expoiting autodiff: or Lord won't you buy me a Laplace approximation, May 2024. URL https: 7/dansbl og. net 1 httptelrr....pppppps///000000000 a ace/ ace
20

---

Published in Transactions on Machine Learning Research (05/2025)
D. E. Stewart. Numerical Analysis: A Graduate Course, volume 4 of CMS/CAIMS Books in Mathematics. Springer International Publishing, Cham, 2022. ISBN 978-3-031-008100-4 978-3-031-00111--7 doi: 10.1007/ 978-3-001-08121-7. URL https://link
J.-W. de Meent, B. Paige, H. Yang, and F. Wood. An Introduction to Probabilistic Programming. 1rXi::1899100566 fcs, stat), Oct. 2021. URL http:://aa......rrraaa//////////
A. Wächter and L. T. Biegler. On the implementation of an interior-point filter line-search algorithm for largo-scale nonlincar programming. Mathematical Programmiing P Mar. 2006. ISSN 1436-4646. doi: 10. 10.1007/111107--0--065558 /ttpa//taaa///ooo..oooooo/oo///0
A. Walther. Computing sparse Hessians with automatic differentiation. ACM Transactions on Mathemalical Software, 34(1):3:113:1 Jan. 2008. ISSN 0098-3500.
doi:
URL https:
00009///|//////////////////////0 //dl. aal..cc..om wttp:////////.ooooo// http:////....11111222555555
A. Walther. Getting Started with ADOL-C. DROPS-IDN//
2009. doi: Documentt//...22///a3eeee......p 09061.1.1
URL
A. Walther. On the Efficient Computation of Sparsity Patterns for Hessians. In S. Forth, P. Hovland,
E. Phipps, J. Utke, and A. Walther, editors, Recent Advances in Algorithmic Differentiation, pages 139 149, Berlin, Heidelberg, 2012. Springer. ISBN 978-3-642-3002333, doi: 10. 10.1007/966-3-642230022
13.
M. Wang, A. Gebremedhin, and A. Pothen. Capitalizing live variables: New algorithms for efficient Hessian computation via automatic differentiation. Mathematical Programming Computation, 8(4):393 433, Dec. 2016a. ISSN 1867-2957. doi: (0.1007/s12533 016-0100-3. 1..10077/12532-006-0000333 URL /ttLtttt/////i.o///100/0
s s12532- 016 01000-
M. Wang, A. Pothen, and P. Hovland. Edge Pushing is Equivalent to Vertex Elimination for Computing Hessians. In 2016 Proceedings of the SIAM Workshop Combinatorial Scientific Computing (CSC), Proccedings, pages 102-111. Society for Industrial and Applied Mathematics, Jan. 2016b. doi: 10.1137/ L.9781611974690.chl. URL httts::/epubb: siam. wwtp:////o///....ooo.ooo 97111/7111.7111119445550 ch11
M. J. Weinstein and A. V. Rao. Algorithm 984: ADiGator, Toolbox for the Algorithmic Differentiation of Mathematical Functions in MATLAB Using Source Transformation via Operator Overloading. ACM Transactions on Mathematical Software, 21(2):211122:255, Aug. 2017. ISSN 0098-3500. doi: 10.1145/
3104990. URL nittps://d) acm. www./oog//oooooooooo 110.0.1155/111999000 04990
J. Willkomm, C. H. Bischof, and H. M. Bücker. A new user interface for ADiMat: Toward accurate and efficient derivatives of MATLAB programmes with case of use. International Journal of Computational Science and Engineering, 9(5-6):408-415, Jan. 2014. ISSN 1742-7185. doi: 10.1504/IJJSE. URL https:://// / http///widerrcceeceoonnoo// http:///wn....ccoa/oo/////////// http:/oi/aab//1..150////aaaa/ 2 06452
21

---

Published in Transactions on Machine Learning Research (05/2025)
A More applications
Probabilistic modeling makes frequent use of Hessian matrices. In frequentist statistics, the Fisher in- formation matrix (Fisher 1925) is defined as the expected Hessian of the negative log-likelihodd: I(0) -Ex-p(-10) 00 V [x O)]- its inverse gives estimate of the variance for asymplolically Gaussian estima- tors. The Bayesian counterpart of this notion is Laplace approximation, whereby the posterior distribution of an estimator is approximated with a Gaussian. The precision of that Gaussian distribution is taken as the observed Fisher information. When the exact Hessian of the log-density is intractable to compute or invert, diagonal approximations are a common workaround. As we witness a shift. from simple models to full-fledged probabilistic programs (van de Meent et al.) 2021), AD becomes key requirement to handle inscrutable log-density functions. Moreover, Markov-Chain Monte-Carlo is family of techniques that allow sampling from high-dimensional, unnormalized densities (Robert and Casella, 2005). To better exploit the geometry of distribution, Hamiltonian Monte-Carlo (Betancourt 2018) incorporates derivatives in the simulation. and those derivatives can be computed with AD. Its Riemannian extension (Girolami and Calderhead 2011 gives a central role to the Fisher information matrix when defining the metric tensor.
B ASD implementations in high-level programming languages
Among scientific computing languages, ASD has percolated most prominently into the MATLAB ecosystem, but the closed-source nature of the language hinders adoption in ML communities. As for Python and R, their ASD libraries are currently less developed than their Julia counterparts.
MATLAB.
The ADMIT package (Coleman and Verma, 2000 relics on an external (C or MATLAB) AD tool to compule sparse Jacobians and Hessians, augmenting it with coloring and compression. However, the chosen AD tool must also be able to detect Jacobian and Hessian sparsity patterns, as is the case for ADMAT (Coleman and Verma 1998) and ADOL-C (Griewank et al.) 1996). ADiMAT (Willkomm et al. 2014) has similar sparse functionality but requires user input for the sparsity pattern. The MAD (Forth 2006 library offers two options for sparse Jacobians and Hessians: either coloring and compression, or direct use of sparse derivative storage inside elementary operations. MSAD (Kharche and Forth 2006) enhances MAD by replacing
operator overloading with source transformation. Finally, ADiGator (Weinstein and Rao 2017 moves much complexity as possible to compile time, in order to lessen the runtime impact of ASD.
Python. The main AD frameworks in Python are TensorFlow (Abadi et al. 2015), PyTorch Paszkk et al.] 2019) and JAX (Bradbury el al. 2018). We are not aware of any ASD libraries for Tensorflow or PyTorch. In JAX, we mostly found sparsejac (Schubert) (024), which is limited to Jacobians and does not support sparsity detection. A blog post by [Simpson (2024) describes potential avenues for sparsity detection in JAX, but they are not fully implemented. The Graphax package (Loholf and Neftei 2024 implements the cross-country elimination technique described in insbsecction 2.4] Since it can exploit sparsity in Jacobian computations, it could potentially be used as an alternative to ASD when vectorized code is involved. Finally, the library auto_diff (Nobol 2020) computes sparse Jacobians of plain NumPy codc.
R. Our search for ASD in R turned up two packages: sparseHessianFD Braun (2017). which asks the user to provide a gradient, and TMB (Kristensen al. 2016), which is focused on one statistical application (Laplace approximation) and requires some of the code to be written in C++.
Julia. The state of the art for ASD used to involve a combination of Symboli i cs. ji Gooda et al. 2019 2022) with SparseDiffToo iftooosfffooo..11 ji (JuliaDiff' contributors, 2024). It was limited to a couple of AD backends, not optimized to support Hessian matrices, and (as our experiments show) rather slow for sparsity detection. As a result, the ecosystem is currently switching to a new, more modular pipeline based on Differentiatio ieaaiinInnteffacc j1 (Dalle and Hill 2025b) and wewareaatrixoooorriio and Montoison 2025).
j1 (Dalle
Remark 3. During preperation of this paper, we became aware of ongoing work around Spadina [Moses] 2023), which relies on dead code removal during compilation to compute only nonzero matriz entries. It is planned to interface with Enzyme [Moses and Churavy 2020 Moses et. al. 2021) and JAX.
22

---

Published in Transactions on Machine Learning Research (05/2025)
C
Coloring basics
brief primer coloring of sparse matrices for ASD, as implemented in The goal is to make this paper self-contained, but refer the reader to Montoison et al. (2025) for more details. In our experiments, we only use greedy column coloring with
Here give
SparseMatrixcoooringg
direct decompression (sce bclow).
C.1 Compression and decompression
As explained in the review by Gebremedhin al. (2005), coloring is a necessary step for matrix compression. Consider matrix A E R=x" and mapping from columns to colors [I....,ni) n) + { Co where is the largest color used. We construct the compressed matrix B E Rmxe. such that each column of B corresponds to a color k. The column B. is obtained by summing the columns of A which share that same color:
B = Σ A., Av
with
a (R)
Σ e
15
IS +0)=
o0)=k
If the coloring & induces injective mapping A + B, then decompression is possible: B contains all the necessary information to recover A exactly. Thus, it is possible to evaluate A using only c, < n matrix-vector products: one for each multi-basis vector (R).
In the general case, decompression from to A might require solving a linear system. avoid paying this cost, common to consider strouger condition called direct decompression: every coefficient A., should be directly readable from some coefficient B₁,-- The question then becomes: what constraints should we impose the coloring & to ensure that direct decompression is possible?
C.2 Coloring asymmetric matrices
The following explanation summarized from Gebremedhin et al. (2005 Section 3). In the asymmetric case, direct dccompression possible if ø induces structurally orthogonal partition of the columns of A: for every non-zero A,,j, column j should be the only one with a non-zero in row i among its color group. With such partition, when we the columns inside the group Φ(j), no other entry is added to A which implies that Biot)) Aij.
Let us represent the matrix A with bipartite graph Go = J.)). Its vertices include the set of rows I= [1,...m] and the set of columns I = {1, ",,,..,,],, while its edges link rows and columns which intersect at a non-zero element: E, = ((i.)) E x J A,, 03. The structural orthogonality condition satisfied whenever á is a partial distance-2 coloring of the columm vertices color to each column such that if j1 and je are distance-2 neighbors in the bipartite graph (there exists such that (i, & and (2) E &,), then they have different colors Φ(11) # 6(j2).
Gr. In other words, á should assign
To achieve this result, a simple greedy algorithm sulfices: iterate over columns in their native order, and give each one the smallest color that is still not used among its distance-2 neighbors. We use this very algorithm in our Jacobian computation experiments.
C.3 Coloring symmetric matrices
The following explanation is summarized from Gebremedhin al. (2005 Section 4). In the symmetric case, direct decompression is possible if induces symmetrically orthogonal partition of the columns of A : for every non-zero Ai,j, either column j should be the only one with non-zero in row i among its color group, or column i should be the only one with a non-zero in row j among its color group. This is a more relaxed notion of partition, which makes sense because we have more opportunities for recovery: the coefficient A,, = Ar can read either from Bi,()) or from Bj
Let us represent the matrix A with an adjacency graph Go Ea), whose vertices are the set of columns (or rows) and whose edges are defined by E. [(i.]) E J2 Aij # 0}. The symmetric orthogonality
23

---

Published in Transactions on Machine Learning Research (05/2025)
condition is satisfied whenever is star coloring of Ga- A star coloring is standard (distance-1) coloring with the additional constraint that every path on four vertices uses at least three colors. In star coloring, every subgraph induced by pair of colors is a collection of stars.
Again, a simple greedy algorithm suffices: iterate over columns in their native order, and give each one the smallest color that is still not used among neighbors, while avoiding the creation of four-vertex three- colored paths. As observed by Gebremedhin al. (2007), representing the graph as set of two-colored stars provides efficient implementation of this idca. We use optimized version of their algorithm in our Hessian computation experiments.
C.4
Other details
Both asymmetric and symmetric decompression routines optimized for Julia's default sparse matrix format, Compressed Sparse Column (CSC).
D Hessian tracing
Here extend the reasoning of section 3 to second-order cross-derivatives.
D.1 Principle
R be voctor-to-scalar function, and R". This time, want the sparsity pattern of the Hessian matrix (mx E ROxn Extending what we did for gradients, we define a second-order tracer lype which contains a primal value y(œ) E R, its gradient sparsity pattern 1[Vy(=)] E {0.1)" and its Hessian sparsity pattern 1[V2 y(x)] E {0, 1]mxx. Again, the Hessian in question is taken with respect to the imput æ, and we will replace every intermediate scalar quantity in with a tracer. We start by turning every Tj into (zy, , 0) (where the third term the initial empty Hessian sparsity pattern), and the end we recover This time, we need to describe how Hessian sparsity patterns propagate.
Let R" +
R"
(F(#), 1[))
Remark Our second-order tracer is related hyperdual numbers, which can be found in second-order fordard-mode AD Fike and Alonso 2012)) More precisely, it describes the sparsity pattern of batched hyperduul number, which could be used to implement second-order batched forword-mode AD. One can also see it as the binary version of the sparse triplet in Dizon et al. 1990).
D.2 Second-ordes propagation rules
We reuse the same framework as for Jacobian tracing in subsectioon 3.2 but this time we go step further. Differentiating Equation I once more with respect to F; gives us:
8/- - = dz, [01. 0,0) + da, 12...00 0,
Now we apply the product rule:
o (d 01p 0ga + 01p 0.0,a] + 02z 0,8 + 024 0,0,8]
We recognize second-order derivatives in the second term of each bracket:
0 = [dz, (81g) a a + 0. 0f + [dz, (02pp) 0,8 + O2p 03
(3)
For the first term of each bracket, we can once again apply Equation 1 but with the differentiated operators 0rip sp and Ozp instead of sp, and with the total derivative dz, instead of dzy;
du, (014) = 0 (01
0, 0 0,8
= 0f.g 0, 0 0,3
dz, (8xp) = 0 (82p) 0,a 02 0,8
0fgz - 0,a + 03
24

---

Published in Transactions on Machine Learning Research (05/2025)
Plugging these into Equation 3 get:
8 = [10fg. 0,a 82g 0,5) 0, + 01 0
+ [0Pg 0,a - 0ip 0.3) 0, 02p 0
And sorting by the operator derivatives involved, we conclude:
of 01p -. 024- 0
(first derivatives)
+ 8..-0000 0,a 8jp 0.3 0,3 + +22 0,0x 0,3 + 0220 - 0,0
(second derivatives)
(cross derivatives)
Bringing the indices logether with vector and matrix notation shows us:
v 31p- a + 02p V23
+ 0pp (Va) Vo) 83 (V))
+ 8 (Va)(V 822 (V
And once again, the sparsity pattern emerges, using V for elementwise OR between matrices. Like before, use < to emphasize the fact that the resulting sparsity patterns overestimates, which may be too conservative.
[0 1
V
1 1 3 1 [ [[()
M
v
( (V)
V
v
[0 [[(Va)(V)))
V
1 (0%44 (V()
Let us generalize V to also represent the outer product OR between two vectors, so that [ab [ This gives our final expression:
[0 1 a
V 1 [824] 1[V23
[[
S
V
1 10 Va V [
V g (1 V (
(4)
10 (1[Va] V (
V 1 10 (IV3
As we can see, the propagation of Hossian sparsity patterns through the operator to only depends on five values:
[] [022))
18f%),
[0 and 100244
These binary values tell whether the operator yo locally depends each of its arguments the first and second order.
Thanks to operator overloading at the scalar level, if the control flow reaches dead end (a value which is not reused for the function output), the corresponding dependencies will not appear in the computed sparsity pattern. This contrasts with the method of /Walther] 2008), where all intermediate values contribute to the final result, leading potential overestimation of the Hessian sparsity pattern.
D.3 Second-order operator classification
To implement [Eqation 4 we need to refine the classification from Table [ by considering second derivatives. Once again, the distinction between local and global sparsity plays key rolc. Some examples are given Table
E Tensor-level overload example
Using the standard matrix multiplication algorithm, propagating tracers through the matrix multiplication C = AB with A paxp and B E 12xx requires p multiplications and n (p 1) additions, since for all n m entries Cijz
Oi = Σ Ar..ra.
(5)
k=
25

---

Published in Transactions on Machine Learning Research (05/2025)
Local
Global
Operator ((a, ())
1
10f
12
1 (0
1 10
10
exp, log
1
+, max, min
0
0
0
0
0
a.c.
1
0
a.e.
a.e.
1
Table Sccond-order classification of operators
Unary operators have no second argument. "a.e." means "almost everywhere" for the Lebesgue measure
Applying the first-order propagation rule from Equation to the scalar multiplication operator
( (2) = 21
we obtain the global propagation rule
1 v 1 [8 1 V V 1[02g] V 1 V V
An identical propagation rule can also be derived for addition. Inserting into Equation 5 we obtain
1[V.,,] = V
which, if naively implemented, requires a total of n m - (2p- 1) elementwise OR operations to propagate tracers through the entire matrix multiplication. Rewriting this as the equivalent
[VC,s = V IVA.VV V Bx
[VA
reveals that we can instcad first compute intermodiate quantitics 1 and 1VB, by taking n (p 1) elementwise operations across rows of A and (p-1) operations across columns of B respectively. The total amount of operations is therefore reduced to (n mp + leading to significant increase in performance.
F
Code demonstration
We now showcase the API of the ASD pipeline used in our experiments section
F.1
Sparsity detection
To demonstrate the generality of SCT's tracer- approach to sparsity detection, we compute the global Jacobian sparsity patlern of a convolutional layer provided by the deep learning framework F Flux. jl
Innes
et al.
2018
Innes 2018
using sparseliemeertiiityrrreeee Flux
import requived packages
FFloatE, 10, 3, 1) layer = Convi (5, 5), => 11
create inpus
preate convolutional layer
detector Trace:9pars
specify global spareity pattern compute patzern
jaccbien_sparsizyy (layer, x, detector)
Listing 1: Detecting the Jacobian sparsity pattern of a convolutional layer using sparscoonectivvityyraeerr,,
26

---

Published in Transactions on Machine Learning Research (05/2025)
The full code is shown in Listing 1 We import the two required packages and sample random input tensor in the size of 10 10 image with 3 color channels and a batch size of 1. We then create convolutional layer Conv with a kernel of size 5 x 5, mapping the 3 input channels to a single output channel. To compute global sparsity patterns,
is chosen 10 The Jacobian sparsity pattern function. Note that SCT doesn't implement custom overloads for convolutional layers. Instead, Flux.j13s generic implementation of a convolution falls back to elementary operators like addition and multiplication, which are overloaded on SCT's tracer types. This is enabled by Julia's multiple dispatch paradigm and doesn't require writing any additional code.
is then computed by calling the FracersparsityDeteccr
The resulting sparsity pattern is shown in Figure 5aj with colors resulting from subsequent greedy column coloring. The banded structure of the matrix results from the size of the convolutional kernel as well as the number of input channels. Figure 5c shows the sparsity pattern that results from increasing the batch size from Since the convolutions of the two inputs in the batch are parallel and separable computations, the resulting sparsity pattern block diagonal matrix. Since this is the case for all separable parallel computations, sparsity detection be used tool to debug the automatic parallelization of computer programs, Figure 5d] additionally increases the number of output channels from 1 to 2. Note that while the size of the Jacobian sparsity pattern incrcases across all three figures, the number of of colors stays constant.
F.2 Computing Jacobians with AD and ASD
Listings 2 and 3 both show the computation of the Jacobian of a convolutional layer from Flux. jil the former using AD, the latter using AST). Since ASD is fully automatic, the only difference in code wetween the two computations lies in the package imports and the specification of the backend used to call to the jacobian function. Possible choices include ForwardDi wewfrddDff.... j1 (Revels et al. 2016). ReverseDi ff. i1
[Reels 2016). Zygote.ji (Innes 2019
Innes el al. 2019) and Enzyme. j1
(Moses and Churavy 2020
Moses et al. 2021
Impost required packages using iww.rrettaaiiinnnnteertet using ForwardDirt using Flux
common interface to AD backends forward-mode backend
deep learning framework
Create input Lensor and convolutional layer x = randa (Ploatt2,, 10, 10, 3, 1)
layer Conv (S, 5), =>
Specify backend
ad_backend = AuzoForwardDs
Compute Jacobian
jaccbian (layer, sackknd, x)
Listing 2: AD computation of the Jacobian of a convolutional layer using Differeentiatiinnnnteeeeaa
* Impost required packages
using maaselcnmeetivvityrraeee
sparsity
using ricooooinnnooooinn
sparsity pottern coloring
using Ditferentiat inttortteeraa
common interface AD backends
using ForwardDiff
forward mode AD backend
using Flux
deep learring framework
+ Create input sensor and corvoluticnal layer x randa randn (21oat32, 10, 10, 3, 1) layer Cont (5, 5),
Specify ASD backend aad bed_backend AutoSparae AutoFozvactDiEEE()
# from SparseConnect ivity
colaring_algoritnn calerieg.lpstttteereeeooesssssee
Foror local sparsity detection, TracerTioca 1 Sparsi used instead of Trasersparsitty2eeecte.
27

---

Published in Transactions on Machine Learning Research (05/2025)
(a) Global sparsity pattern of Conv ( ((5, 51, 3 => 1) using batch size 1
(b) Local sparsity pattern of Conv ((5, 5), 1) using batch size
(c) Global sparsity pattern of Conv (5, 5), using batch size 2
(d) Global sparsity pattern of Conv (5, 5), using batch size 2
Figure 5: Jacobian sparsity patterns of small convolutional layers from Flux. j1 applied to 10 10 image with three color channels. Squares correspond to non-zero entries in the Jacobian, with colors resulting from greedy column coloring.
28

---

Published in Transactions on Machine Learning Research (05/2025)
Compute dacobian jacebian (layer, asc_backend,
Listing 3: ASD computation of the Jacobian of a convolutional layer using j1 and
Di Differentiattionnnterra sparsetooeeiiiiiyyyrrrccc htt.:///weerraa....ppporrooeeeee
SpersetatrixColorings. j1
G More Jacobian experiments:
Brusselator
G.1 Sparsity detection
display the colored sparsity patterns of the Brusselator PDE (see section 5) for various discrctization levels. Writing explicit formulas for this kind of intricate structure is fairly tiresome and subject to human error, which justifies automated spansity detection approaches.
In Figure 6
G.2
Jacobian computation and linear solves
In Table 4
we present two sets of results:
comparison between AD and ASD for Jacobian computation, already discussed in section 5)
comparison between various methods for solving a Jacobian linear system, which we now focus on. Solving linear systems involving Jacobian matrices typically useful as part of Newton root- finding step. This happen within larger iteration, like backward Euler method for differential equations.
We distinguish between three different techniques for the Newton step. First, one needs to choose between direct solvers (which factorize the matrix) and iterative solvers. Second, one picks an iterative solver, there is still a choice between precomputing the matrix or using the lazy JVP operator. Our specific implementation uses Julia's standard library SparseArrays.ji for the direct solve, which in turn relies SuiteSparse. Mcanwhile, the iterative solver is the default recommendation BiCGstab (1) from ieerativeSolveee jl (Chen et al. 2013). In the rightmost part of Taale 4 we observe that the lazy version wins for small instances, but then the direct solver is actually faster for large enough matrices.
Problem
Jacobian computation
Newton step
AD (prepared)
ASD (prepared)
ASD
JVP (iterative)
Jacobian (iterative)
Jacobian (direct)
N
2
(unprepared)
6
1.64 10
1.97
10
(8.3)
3.36
10-5 (0.5)
2.07 10 5
2.19
10-
4.52
10
12
2.44-10-
8.67 10-6
(28.1)
10-0
(1.4)
1.61
10-4
2.42 10-4
24
10-3
3.43
10-5
(117.2)
1.31
(3.1)
1.04
1.34 10-3
1.24
10-3
48
10-2
1.68
10
(451.4)
1.70
(4.5)
8.34
1.17-
8.98
10-3
96
1.35
10
6.68
10
(2017.2)
2.16
10
(6.2)
7.56
10
1.11- 10-1
4.07
10
192
2.25 1
3.09
10-3
(7293.5)
4.62 10
(4.9)
1.05
10
1.07 10
2.30
10-1
Wall time seconds.
In parentheses: Wall time ratio compared to prepared (higher better).
Table 4: Performance comparison of Jacobian computation & linear solves the Brusselator PDE.
While this behavior might seem surprising, it is casily explained by the principles underlying ASD and iterative solvers. Assuming that the function at hand is very expensive to compute, pure lincar-algebraic manipulations will be much cheaper in comparison, sophisticated ones like Sui SuiteSparse linear solve.
29

---

Published in Transactions on Machine Learning Research (05/2025)
(a) =
(b) N
(c) N =
(d) N 6
Figure 6: Jacobian sparsity patterns of the discretized Brusselator PDE of size N x N x 2. Squares correspond to non-zero entries in the Jacobian, with colors resulting from greedy column coloring.
30

---

Published in Transactions on Machine Learning Research (05/2025)
Thus, it all comes down to how many JVPs need to be performed. Iterative solvers require one JVP per iteration, so the total complexity is dictated by the precision requirement and the conditioning number of the Jacobian. Meanwhile, materializing Jacobian with ASD requires one JVP per distinct color, so the total complexity is dictated by the sparsity pattern only. Thus, in cases where high precision is required, the matrix ill-conditioned and the coloring number is low, ASD with a direct solver can prevail over lavy iterative counterpart.
H Hessian experiments: ACOPF
H.1 Sparsity detection
At the core of power systems planning and power markets lics a nonlinear constrained optimization problem called the allernating current oplimal power flow (ACOPF) problem. As of today, the full ACOPF remains unsolved. While approximate solution techniques exist, they result in unnecessary emissions and spending, with potential anmual savings from optimal solution estimated in the tens of billions of dollars (Cain et al. 2012). The problem has also attracted the interest of the ML community with recent developments in neural ACOPF solvers (Piloto et al., 2024). To validate algorithms for the OPF problem, a suite of benchmarks called Power Grid Lib (PGLib) has been developed by Babxeiinejadsarookkaee et al. (2021). which can rum via the rosetta-opf (Coffrin and Dowson 2022) implementation of the ACOPF.
We benchmark the sparsity detection of Symbol i ics. j1 and SCT on the Hessian of the Lagrangian of several PGLib optimization problems. The results are summarized in Table 5 SCT outperforms Symbolics.
every problem, regardless of size and sparsity.
H.2 Hessian computation
We benchmark the computation of dense and sparse Hessians for the Lagrangian of several PGLib optimization problems. Our methodology with respect to prepared and unprepared computations mirrors that used in subsubscction 5.1.2 In both AD and ASD benchmarks, ForwardDiff. jl is used over the reverse-mode backend ReverseDif ff.jl to evaluate HVPs.
The results are summarized in Table 6 For large problems, prepared ASD provides an increase in perfor- mance of three orders of magnitude over AD. Even one-off unprepared ASD provides performance benefits over AD in all PGLib cases. Once again, performance gains of one-off ASD on small problems are largely due to the performance of SCT: timings of the 3 3_Imdd problem in Table 6 reveal that just the sparsity detection of Symbolics. jl alone was previously less performant than the full computation of the Hessian with AD.
31

---

Published in Transactions on Machine Learning Research (05/2025)
Problem
Sparsity
Sparsity detection'
Name
Inputs
Zeros
Colors
Symbolies
SCT
3
Imbd
24
91.15%
6
1.29
10-3
5.59
10-6
(23.1)
5 _jm
94.99%
8
2.49 1
1.19
10-4
(20.9)
14 jece
118
97. 97.84%
10
9.02-10-3
5.19
10-4
(17.4)
24 deee ieee_rts rts
266
99. 29.2
12
1.50 1.50-10-
(12.8)
30
236
98.89%
12
2.04 10-2
1.60 1.60-10-3
(12.7)
30 feee
236
98.89%
12
1.60 1.60-10-3
(12.7)
39 epri
282
99.10%
10
14-2
2.10-10
(11.7)
57 ieee
448
99.41%
14
4.68.10-2
(9.5)
60 C
518
99.56%
12
5.15-10-2
5.76-10-
(8.9)
73 jeee ris
824
99.74%
12
9.84-10-2
1.09-10-2
99.0)
89 pegase
1042
99.74%
26
1.80
2.20
(8.2)
118 ieee
1088
99.77%
12
1.57.10-1
2.11 10
(7.5)
162 icee dtc
1484
99.82%
16
2.99.10-1
3.33 -
99.0)
179 goc
1468
99.83%
14
2.59.10-1
3.09 10
88.4)
197 snem
1608
99.85%
14
3.02-10-1
-.57-10
88.5)
200 activ
1456
99.82%
12
2.59
2.92 10
88.9)
240 psere
2558
99.91%
16
6.72
7.32 10-2
99.2)
300 ieee
2382
99.80%
14
6.20
10-2 6.95
88.9)
500 goc
4254
99.94%
14
-0
1.40
10-0
(12.9)
588 sdet
4110
99.94%
14
1.71
10 1.40
(12.2)
793 goc
5432
99.95%
14
10"
2.65 10-1
(11.2)
1354 pegase
11192
99 .9.98%
18
-
4.14
(38.1)
1803 snem
15246
99.98%
16
3.02-
7.17
10-1
(42.1)
1888 rte
14480
99.98%
18
2.72
10 6.53
(41.7)
1951 rte
15018
99.98%
20
6.50
10-1
(47.7)
2000 goc
19008
99.99%
18
-
-
(59.5)
2312 9oc
17128
99.98%
16
4.43
8.69 -
(51.0)
2383up k
17004
99.98%
16
1
-
(51.7)
2736sp_4 k
19088
99.99%
14
-
1.02 10
(62.1)
2737sop
18988
99.99%
16
5.62
100
(55.1)
2742_900
24540
99.99%
14
1.37
10
-00
(122.8)
2746wop
19582
99.99%
16
6.61
10d
1.06-100
(62.5)
2746wp_
19520
99.99%
14
6.35 10
1.04-10"
(60.9)
2848 rie
21822
99.99%
20
8.57 10a
1.23 1.2
10
(69.5)
2853 sdet
23028
99.99%
26
1.04
8
(121.2)
2868 rte
22090
99.99%
20
1.10
-02
(86.7)
2869 pegase
25086
99.99 99.99%
28
-
1.06-10
(136.2)
301 21ppp__
21082
99.99%
14
8.
1.22
(68.4)
3022 goc
23238
99.99%
18
1.45
9.83
(147.8)
3120sp_ k
21608
99.99%
18
9.46
103
(72.1)
3375up k
24350
99.99%
18
1.27
10
9.82 10-1
(128.9)
Wall time sccouds.
Number colors resulting from greedy symmetric coloring.
parentheses: Wall time ratio compared to Symbolics.jlss sparsity detection (higher better).
Table 5: Performance comparison of Hessian sparsity detection the Lagrangian of PGLib optimization problems.
32

---

Published in Transactions on Machine Learning Research (05/2025)
Problem
Sparsity
Hessian computation'
Name
Inputs
Zeros
Colors
AD (prepared)
ASD (prepared)
ASD (unprepared) 3
3 Imbd
24
91.15%
6
1.82 10-4
8.29 3.29-10-6
(2.2)
1.45 10-4
(1.3)
5 pjm
44
94.99%
8
6.33 10
17
(3.7)
3.03 10
(2.1)
14_ feee
118
97.84%
10
5.38
4.84-10-4
(11.1)
1.12- 10-3
(4.8)
24_ deee_ rts
266
99.22%
12
2.56-10-* 10
1.04-10
(24.7)
10-3 2.74-10
(9.3)
30 as
236
98.89%
12
2.39.
1.10-10-3
(21.8)
2.84.10-
(8.4)
30 jeee
236
98.89%
12
-
1.09-10-*
(21.6)
2.87-10-3
(8.3)
39 19_err
282
99.10%
10
10
1.21-10
(27.1)
3.43
(9.6)
57 jeee
448
99.41%
14
8.
3.96-10-1
(22.2)
9.23-10-3
(9.5)
60
518
99.56%
12
1.15.10-1
2.36-10-3
(48.6)
8.61
(13.3
73 feee ris
824
99.74%
12
20
-.47-10-3
(79.1)
10-
(17.8)
89 pegase
1042
99.74%
26
5.61-10-
1.61-10
(34.8)
4.28-10-2
(13.1)
118 jeee
1088
99.77%
12
5.55 10
-
(105.8)
3.13.10-2
(17.7)
162 _ice_eee_ dtc
1484
99.82%
16
1.16 10
-
(75.7)
5.5
(00.9)
179 gue
1468
99.83%
14
1.08 10
1.3
(81.3)
50-
(21.4)
197 snem
1608
99.85%
14
1.34 10
1.46-10-
(92.2)
-0-
(23.0)
200 activ
1456
99.82%
12
1.02 10 10
6.94-10-3
(146.6)
3.88 3.88-10-2
(26.3)
240 psere
2558
99.91%
16
3.51 10"
2.50-10-
(140.2)
10-4
(33.6)
300 icee
2382
99.89%
14
3.00 10 10
2.14-10-2
(140.3)
9.67.
(31.1)
500 goc
4254
99.94%
14
-
3.85-10-2
(307.3)
2.20 10
(53.7)
588 sdet
4110
99.94%
14
1014
3.60-10
(316.1)
2.14.10-1
(53.3)
793 goe
5432
99.95%
14
10
4.91-10-2
(443.1)
3.33 10-1
(65.3)
1954 pegase
11192
99.98%
18
1.36 10
1.21-10-1
(1128.4)
6.21 10
(219.6)
1803 snem
15246
99.98%
16
2.09
1.66-10-1
(1259.5)
1.07 10"
(195.0)
1888 rte
14480
99.98%
18
8.15
1.43
(5706.7)
8.76
(930.4)
1951 rie
15018
99.98%
20
-
1054-10-1
(1293.4)
1.00
(199.1)
2000 goe
19008
99.99%
18
3.58 10a
2.15-10-1
(1669.5)
1.61 10t
(222.7)
2312 goe
17128
99.98%
16
2.75 102
10-
(1470.7)
1.35 10
(204.5)
2383up_
17004
99.98%
16
2.65 10a
1.80-10-1
(1468.2)
1.14 100
(231.4)
2736sp_
19088
99.99%
14
3.30 10
1.78-10-1
(1857.2)
1.40 10
(235.5)
2737sop
18988
99.99%
16
3.29 10a
2.02-10
(1629.8)
1.47 10
(223.0)
2742_ goc
24540
99.99%
14
6.50 10a
2.41
(2694.1)
1.78 10
(366.3)
2746wop k
19582
99.99%
16
3.64 10a
2.07-10-1
(1755.7)
1.54 10
(235.6)
2746wp_
19520
99.99%
14
3.53 2
1.77-10
(1991.4)
1.51 10
(234.5)
2848 rie
21822
99.99%
20
4.67 10
2024-10-1
(2083.5)
1.80 10
(259.7)
2853 sdet
23028
99.99%
26
5.38 102
3.62-10-
(1486.9)
1.68 10th
(320.6)
2868 rte
22090
99.99%
20
5.02 10a
2.35-10-1
(2137.9)
1.73 10
(290.0)
2869 pegase
25086
99.99%
28
5.08 10a
4.07-10-1
(1249.0)
1.99 100
(255.5)
3012wp k
21082
99.99%
14
4.33 102
1.96
(2208.3)
1.77 10
(245.1)
3022 goc
23238
99.99%
18
5.76 10a
2.51
(2296.9)
1.48 10
(390.7)
3120sp_
21608
99.99%
18
4.56
10
2.26-10
(2019.2)
1.90 10
(240.1)
3375wp_
24350
99.99%
18
6.25
10
2.54 2.54-10-1
(2463.9)
1.71 10
(365.1)
Wall time in seconds.
Number of colors resulting from greedy symmetric coloring.
In parentheses: Wall time ratio compared to prepared AD (higher is better).
Table Performance comparison of AD and ASD Hessian computation on the Lagrangian of PGLib optimization problems.
33


Reviews:
[{"review_id": "NgfePjWm8D", "reviewer": "Reviewer_PDA9", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The paper is very well presented with a great description of earlier work.\n- The presentation of sparsity detection patterns with operator overloading is very clear.\n- The experiments demonstrate the benefits of the approach in Julia\n\nWeaknesses:\n- The whole pipeline is not presented. In my opinion, the matrix coloring part should be presented to give a full overview. Since the contribution of the paper is the pipeline and not a new algorithm, an incomplete pipeline makes the paper simply incomplete.\n- Experiments can be made more extensive. As the authors mention at the beginning, \"matrix-free\" methods have been largely adopted to avoid full Jacobian/Hessian computations. It would be great to see the performance of the proposed package in downstream tasks like solving a linear system. By having access to the whole Jacobian, the authors may benefit from methods that are more numerically stable for example. It's unclear though it the cost of detecting the sparsity pattern could be amortized in this regime. Full experimental details on such questions would largely benefit the manuscript.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents a new pipeline to compute efficiently sparse Jacobians or Hessians in a differentiable programming framework. To be exact, the underlying algorithms are not new but the package itself, coded in Julia, is new, versatile and demonstrates improvements compared to existing libraries in Julia. \nThe authors explain in detail the benefits and applications of automatic differentiation libraries able to take advantage of sparse Jacobians/Hessians. They also provide a detailed related work of this old subject. \nThe authors focus, in this manuscript, on detecting sparsity patterns in the Jacobian and give a new presentation on the subject based on operator overloading. A second essential step (matrix coloring) is relegated to another manuscript. \nExperiments demonstrate that the new package outperforms previous packages by an order of magnitude in speed in terms of sparsity pattern detection. Experiments also demonstrate that the new package is competitive at computing the whole Jacobian, even without considering amortizing the cost of sparsity detection over several Jacobian computations.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Main changes:\n- I believe the manuscript would benefit from being fused with the ones that would describe the matrix coloring and the overall automatic differentiation interface. Otherwise the authors will end up with three incomplete papers with clear overlap of contributions and even maybe content. As much as I appreciate the current manuscript, I think such practice (dividing papers in smaller papers), may not be well regarded by the community.\n- Add experiments such as:\n  - Speed to make one optimization step in the ACOPF problem. More generally, comparison between \"lazy products\" and the proposed package to solve linear system like $Ha =b$ or $Ja = b$. I suppose the ACOPF problem solves such a linear system every step so it can be a starter. In optimization, the sparse Jacobians could particularly shine to solve nonlinear least-square problems where iterative algorithms use linear systems defined typically by the Jacobian. \n  - Comparison with closed source package if possible.\n\nClarification questions:\n- About the prepared ASD case, since the goal is to compute just one Jacobian/Hessian matrix how can the costs be amortized? For global sparsity patterns, it is clear that the global patterns can be reused to compute sparse Jacobians/Hessians at other points. But the local sparsity patterns, as presented by the authors, will depend each time on the point at which the Jacobian/Hessian is computed and cannot be amortized no?\n- \"batched evaluation\" refers to vmap in jax?\n\n\nDetails\n- There has also been a recent efforts about optimizing automatic differentiation engines through reinforcement learning (using vertex elimination and alpha zero) [1]. Comparing the present library is probably out of scope but worth mentioning nevertheless\n- The authors may consider avoiding phrasing that do not age well like \"most advanced automatic sparse differentiation system\" (just remove \"most\"). Similarly do not refer to \"future papers\"\n- page 3 top, \"Python currently lack such tooling\" what tooling exactly? Automatic sparse differentiation? or multiple dispatch paradygm?\n- Spectral analyses of Hessians or Jacobians could also largely benefit from having access to the actual matrix. In deep learning, there has been a growing interest in understanding the spectral properties of the Hessian of the objective along training for example [2]. It could be worth mentioning it.\n- \"It is worth mentioning that some optimization solvers only accept explicit Jacobians/Hessian matrices\" -> give references. Also is it a bottleneck in terms of code or algorithm?\n- \"optimization solvers often welcome matrices in CSC or CSR format\" -> give reference\n- While finding local sparsity patterns is clearly explained through operator overloading, finding global sparsity patterns is not really explained.\n- \"This allows us to bypass the original scalar computational graph... \" -> this sentence is not clear. More importantly the scope of the tensor-overloads part is not clear. It would be important for the reader to understand whether future work is needed to optimize the approach for convolutions.\n- Just before section 5: \"An demonstration\" -> \"A demonstration\"\n- Appendix A.5 \"In the end, the right choice of set implementation will depend on ...\" What set implementations the authors are talking about? Is it referring to various sparse encoding formats?\n\nSuggestions (feel free to ignore):\n- May consider a table summarizing the different approaches for automatic sparse differentiation and their main differences if there are.\n- Maybe mention in the notations that the operator $\\cdot$ has priority over the operator $\\vee$ in the diverse mathematical formulations.\n- The sentence \"We can apply the previous remarks on ...\" is somewhat weirdly formulated. Could be replaced by two sentences \"The gradient $g = \\nabla f(x)$ can be computed by reverse mode AD. Forward mode AD can then be applied on the computation of $g$\"\n- May put the coding block (page 23) in main text, since this manuscript is presenting the library.\n\n[1] Optimizing Automatic Differentiation with Deep Reinforcement Learning, Lohoff et al, NeurIPS 2024\n[2] Gradient Descent on Neural Networks Typically Occurs at the Edge of Stability, Cohen et al, ICLR 2021", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "x1gBExlMog", "reviewer": "Reviewer_ksEs", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "As the authors mention in the paper, the Automatic Differentiation and Machine Learning communities are quite separated. Therefore, it must be noted that my review is provided from the perspective of a member of the ML community. This implies that I might not be familiar with some methods developed by the AD community, and I approach these algorithms from the perspective of a user.\n\nStrengths:\n- The paper provides a clear exposition of the previous works and the proposed method. It is immediately clear where one can expect the acceleration of the computation and why.\n- The method demonstrates very strong empirical results on the selected tasks. That is, there is no question that the proposed method improves the efficiency of the computation and will do so for different hardware, problems, and hyperparameters. This simplifies comparison and serves as a good motivation for the usage of the proposed method.\n- The authors supplement an open-source implementation of the proposed method as a single Julia library. This is essential for achieving the declared goals of disseminating advanced AD techniques among the ML community.\n\nWeaknesses:\n- The Julia implementation is not the best or the fastest way to spread your method among the ML community. Although I don't have a clear understanding of this design choice, this is a clear barrier that will prevent many practitioners from using the proposed method.\n- The CPU-based implementation does not help achieve the declared goals. Indeed, as the authors mention in their work, CPU-based solutions are currently not the main focus of the ML community. Furthermore, the parallel nature of computations on GPU raises concerns about the possibility of implementing the proposed algorithm on GPU and the potential benefits from ASD.\n- The paper lacks code examples explaining to the reader how to use the proposed framework. Since the goal of the paper is, in part, educational, it would be more convincing if the authors demonstrated how one can easily use the proposed method instead of conventional AD methods.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a new exposition and implementation of the Automatic Sparse Differentiation (ASD) based on the operator overloading concept. In detail, the authors argue that there exists a certain barrier between two scientific communities: the community developing Automatic Differentiation (AD) algorithms and the Machine Learning (ML) community. In particular, Automatic Sparse Differentiation is still not widely adopted by the ML community despite being a very promising direction for optimizing the efficiency of ML algorithms. \n\nIn order to bridge the gap between the communities and make ASD more practical, the authors provide a minimalistic description of the operator overloading approach and unify several previous ASD methods in a single Julia library. They provide an open-source implementation of the library and benchmark it on several tasks relevant to the ML community. The benchmarks demonstrate the great utility of ASD with operator overloading in terms of the wall time required for the Jacobian and Hessian computations.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The major change I request is the addition of verbatim code snippets demonstrating the code difference between AD and ASD. This would greatly improve readability and help achieve the objective of promoting advanced AD techniques.\n\nThe list of minor comments:\n- abstract, page 3, the following phrase requires an explanation \"one-off computations of Jacobians\"\n- page 4, second to the last paragraph, typo in \"to a select few AD systems\";\n- page 6, the usage of color in eq. (1) and the one before is not clear, i.e. what exactly the colour show here and why is the colour not used in other equations?;\n- page 8, first line, typo in \"$x_j$ into $x_j$ into\";\n- page 9, last paragraph of Section 4, typo in \"Not only is it more generic with respect to the AD backend, it also much faster\" and \"An demonstration of the code\";\n- Figure 2, the resolution of the presented sparsity patterns is very low, in stark contrast with Figure 1;\n- page 11, the following fact is not evident from Table 3 as suggested \"sparsity detection alone took more wall time than the full AD Jacobian computation\";\n- Tables 3 and 4, it is not clear what \"pattern detection\" stands for. I would suggest being consistent with the text where the authors use \"sparsity detection\";\n- why does Table 3 have only pattern detection wall time and does not have a full computation wall time comparison?", "broader_impact_concerns": "The paper does not raise any broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "U0wJiWbBCS", "reviewer": "Reviewer_BWxp", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n* This paper brings Automatic Sparse Differentiation to the machine learning community. It can be convenient for problems that require the derivation of sparse Jacobian/Hessian matrices and for which computing some Jacobian/Hessian-vector products is not sufficient. \n\n* The authors propose an efficient implementation of the ASD pipeline in Julia.\n\n* An extensive benchmark of ASD against AD is provided, showing the efficiency of ASD in practice.\n\n* The paper is very well written and the presentation is clear.\n\n\n### Weaknesses\n* See the requested changes section for minor issues.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper deals with Automatic Sparse Differentiation (ASD), a set of techniques that leverage the sparsity of Jacobian/Hessian to speed up their computation. The paper presents propagation rules for sparsity pattern detection to approximate the sparsity pattern of Jacobian and Hessian matrices, thanks to a binarization of the chain rule. These propagation rules are the cornerstone of efficient sparsity detection and are implemented in a Julia package, as well as the matrix coloring and differentiation steps. The sparsity detection implementation for Jacobian and Hessian is benchmarked against another sparsity detection technique. Moreover, the full ASD pipeline is numerically compared with AD.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "### Minor\n* **p.9 and p.10**: Since the reference (Chen et al., 2018) is straddling on pages 9 and 10, the hyperlink to the reference covers the footnote of page 9 and the header of page 10. A solution is proposed [here](https://tex.stackexchange.com/questions/330669/cite-hyperlink-across-two-pages).\n\n\n### Clarity\n* **p.6 and p.7**: The unfamiliar reader can be confused by the appearance of $y$, without being defined. The author should clarify that $y(x)$ is an intermediate value in the computation of $f_i(x)$. Moreover, I don't get why when it is written \"$y\\in\\mathbb{R}$\" there is no argument $x$ and when it is written \"$1[\\nabla y(x)]\\in\\{0,1\\}^n$ the argument is written.\n\n* **Notation confusion in sparcity pattern**: I wonder if there is confusion in the notation between the actual sparsity pattern $1[\\nabla \\gamma]$ (as defined in the notation section) and the \"approximate\" sparsity pattern contained in the tracer. Indeed, if I understand correctly, the equation (2) tells us that $1[\\nabla \\gamma] = 0$ if and only if $1[\\partial_1 \\varphi]\\cdot 1[\\nabla \\alpha] = 1[\\partial_2 \\varphi]\\cdot 1[\\nabla \\beta] = 0$. This is not true if $\\varphi(\\alpha, \\beta) = \\alpha + \\beta$ and $\\alpha(x) = x = -\\beta(x)$. In this case, for $x\\neq 0$ we have $1[\\nabla \\gamma(x)] = 0$ and  $1[\\partial_1 \\varphi]\\cdot 1[\\nabla \\alpha] = 1[\\partial_2 \\varphi]\\cdot 1[\\nabla \\beta] = 1$. What is true is that\n$\n1[\\nabla \\gamma] \\leq 1[\\partial_1\\varphi]\\cdot 1[\\nabla \\alpha] \\vee 1[\\partial_2\\varphi]\\cdot 1[\\nabla \\beta]\n$. This tells us that the propagation rule in the right-hand-side is not exact but it does not approximate by 0 a non-zero coefficient of $1[\\nabla \\gamma]$, which is crucial.\nThus, I think that in equations (2) and (3) (for the Hessian), there should either be a $\\leq$ instead of $=$ or the notation should be changed to reflect that the sparsity pattern is not exact but an approximation.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "NdkmloXGTJ", "reviewer": "Reviewer_vWkC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**:\n\n- The SCT sparsity detector proposed by the authors is two to three orders of\n  magnitude faster on benchmark problems than previous software in Julia.\n\n- Experiments show that SCT reduces to the overhead of ASD, making it faster\n  than idealized dense automatic differentiation. \n\n- The description of forward/reverse mode differentiation and ASD is clear.\n\n\n**Weaknesses**:\n\n- The paper describes only the first part of the ASD pipeline, leaving matrix\n  coloring and the differentiation interface to \"future papers\".  \n\n- No experiments are given for machine learning (ML) applications. \n\n- The conceptual novelty and implementation difficulty for the authors ASD\n  approach is not clear.\n\n\n### Additional Details\n\nThis is an interesting paper which outlines a new software package for\nleveraging sparsity in automatic differentiation. \nThe paper is well-written and provides clear general description the\nASD package and its theoretical principles. \nThe package itself is written in Julia, a widely used language in the optimization and\nscientific computing communities.\nWhile the experimental evaluation is rather small in scope, it shows that SCT \noutperforms existing software by several orders of magnitude.\nTaking these items in to account, I think the package is of sufficient interest\nto the TMLR community. \nUnfortunately, the positive aspects of the paper are outweighed by my concerns\nregrading the paper's weaknesses.\nThese issues, which are discussed in detail below, must be reasonably addressed\nbefore I can support accepting this paper.\n\n**Small Scope**: \nThis paper only describes the SCT sparsity detector and explicitly leaves the\nrest of the ASD pipeline to future work.\nThis has three effects: (i) it makes it difficult to weigh the contributions\nof the sparsity detector compared to the other elements of the pipeline; (ii)\nit makes it seem like the authors are trying to \"slice their contributions\nthin\" and maximize the number of publications obtained from this work; and\n(iii) it harms the pedagogical value of the paper as the reader learns nothing\nabout matrix coloring or the differentiation interface.\nI strongly urge the authors to consider revising their paper to include\ndescriptions of all three software components. \nThese could easily be fit into the paper by paring down on the five pages of\nintroduction/related work and compressing verbose parts of the text. \n\n**ML Experiments**: \nIf the authors only wish to describe SCT, then I believe it is critical that\nclear applications of the package are developed for ML problems.\nAs the authors note, automatic differentiation is a key tool in ML and many ML\nproblems exhibit Jacobian/Hessian sparsity. \nExperiments showing that faster ASD has a positive down-stream affect on ML\ntasks would make an excellent argument for the utility of the software package.\nIt would also increase the papers contributions, which I consider to be quite\nthin.\nOne obvious problem the authors could consider is training ReLU networks; as\nthe text notes, ReLU activations induce Jacobian sparsity, so perhaps ASD would\nbe useful here.\n\n**Contribution Size**:\nAutomatic differentiation is not my area of expertise and I found it hard to\nquantify the novelty/interest of the software using the text.\nThe authors comment that they present a \"revised viewpoint\" on sparsity\ndetection, but this viewpoint seems quite standard and intuitive to me given \nfamiliarity with forward-mode differentiation.\nSimilarly, the authors make no comment on the actual implementation of the package.\nWas significant code/performance optimization necessary for SCT to be fast?\nWhat were the main challenges in the implementation and why are they interesting?\nAnswering these questions would clarify the core contributions of the paper.\n\n### Minor Issues:\n\n- Page 1: \"...resting assured that gradients will be computed efficiently and\n  correctly without human intervention.\" --- Actually, there are several\n  well-known bugs in automatic differentiation systems as result of iterative\n  application of sub-gradient computation rules. For example, functions like\n  $\\min\\{x, \\max\\{x, 0\\}\\} = x$ are often differentiated incorrectly in PyTorch.\n\n- Page 2: \"...most advanced automatic sparse differentiation system\" --- what\n  do you mean by \"most advanced\"? Does this mean fastest, best at detecting\n  sparsity, etc? I think this statement should be made based on clear metrics.\n\n- Page 5: \"Specifically for Hessian patterns, Walther (2008) extends the\n  operator overloading approach...\" --- it is more typical to use \"they\" as a\n  genderless singular pronoun when referring to papers with only one author.\n\n- Eq. (1): Equations are part of a sentence like any other textual element.\n  This equation stands alone as its own sentence, which is slightly odd.\n\n- Page 7: \"These binary values us whether the\" --- missing word?\n\n- Perhaps the authors should consider submitting this work to a more\n  software-focused venue, such as JMLR's [machine learning open source\n  software](https://jmlr.org/mloss/mloss-info.html) track.  While I don't think\n  this paper is necessarily a bad fit for TMLR, it may receive more attention\n  at software-focused venues.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a new automatic sparse differentiation (ASD) framework for\nthe Julia programming language. \nThe ASD framework has three components: a package for propagating\nJacobian/Hessian sparsity patterns during forward passes, a matrix coloring\npackage for determining the number of necessary vector-Jacobian products,\nand an interface for differentiation.\nThe authors primarily describe how to propagate sparsity patterns, leaving\ndescriptions of the other two packages to future work.\nDuring the forward pass, sparsity patterns are computed from parent nodes\nin the compute graph using local propagation rules, much like the \"seed\" in\nforward mode differentiation.\nThe output value of the forward procedure gives the sparsity pattern of the\nJacobian/Hessian and can be used for matrix coloring and compressed\ndifferentiation. \nNumerical simulations show that this approach is much faster than previous\nsoftware for automatic sparsity detection in Julia.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Questions**:\n\n- What is the motivation for presenting `SparseMatrixColorings.jl` and\n  `DifferentiationInterface.jl` in separate (future) papers?\n- Can SCT significantly improve training time for ML problems with sparse\n  Jacobians, such as training ReLU networks? \n- Was significant code/performance optimization necessary for SCT to be fast?\n- What were the main challenges in the implementation and why are they\n  interesting to the TMLR community?", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Overall, this paper is of great pedagogical value to the ML community, and while quite narrow, the experiments demonstrate the effectiveness of the proposed package.

The GNN experiment is the one that will speak the most to the ML community, so it should be part of the main text.

I share the reviewer's concern about splitting the contributions: the camera-ready version of the manuscript should include a description of the colouring algorithm; this paper must be self-contained.
Decision: Accept (poster)