REVIEW 3 major objections 5 minor 1 cited by
Scalable Analysis and Design Using Automatic Differentiation
T0 review · 3 major / 5 minor · reviewed 2026-08-07 · deepseek-v4-flash
Pith's one-line read Localizing automatic differentiation to quadrature points gives exact finite-element Jacobians at nearly hand-coded speed.
desk verdict The core AD-localization idea and benchmark are solid, but the blanket claim that B is coordinate-independent overreaches for shape optimization, and the experimental evidence is narrower than the abstract promises. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing mechanism is the Finite Element Operator Decomposition. It splits the assembled operator into $P$ (global to subdomain), $G$ (subdomain to element), $B$ (element to quadrature points), and $D$ (pointwise evaluation of constitutive relations). The decomposition isolates the only nonlinear, problem-dependent part into $D$, so AD can be applied there alone. The resulting chain-rule identity $J_p = P^T G^T B^T J_D BGP$ is the formula that lets a solver compute exact tangent operators without differentiating through mesh transfer or parallel data exchange.
What would settle it
Run the integration-point-only Jacobian on a curved isoparametric element or a moving mesh, where $B$ depends on coordinates or design parameters, and compare it with a finite-difference Jacobian of the full residual. A mismatch would falsify the claim that localizing AD to $D$ captures the exact derivative.
Extended reading notes
Core claim
The authors claim that exact Jacobians for nonlinear finite element systems can be obtained by localizing AD at the quadrature point. Writing the residual action as $A_p(u)=P^T G^T B^T D(BGPu)$ and differentiating gives $J_p(u)=P^T G^T B^T J_D(BGPu)BGP$, so AD only ever differentiates the pointwise operator $D$. Since $P$, $G$, and $B$ are linear and fixed, the identity is exact and requires no changes to parallel communication. The paper backs the claim with timings: for third-order tetrahedra, forward-mode AD at the integration point is 10-15% slower than a hand-coded tangent assembly, whereas element-level AD, especially reverse mode, is markedly slower and more memory-hungry. This is what makes derivative computation both automatic and scalable.
Load-bearing premise
The argument assumes $P$, $G$, and $B$ are fixed linear operators that do not depend on the solution, the mesh coordinates, or the design parameters; if any of them does, derivatives through that operator are missed and the integration-point-only Jacobian is no longer exact.
Editorial extensions
If this is right
- Tangent matrices for new nonlinear constitutive models can be generated automatically, eliminating hand-derived derivative formulas.
- Parallel scalability is preserved because the AD loop never touches global transfer or communication code.
- Adjoint loads for topology optimization and other PDE-constrained design can be computed from pointwise code alone, letting new objectives and constraints be added quickly.
- The approach is independent of the AD engine chosen, and a native dual-number implementation matches the performance of external AD libraries.
Reading between the lines
- A likely extension is to time-dependent and multiphysics problems, where only the pointwise operator $D$ changes; the paper argues this but demonstrates only a steady scalar problem.
- For shape optimization, moving meshes, or isoparametric elements, $B$ becomes design-dependent, so the “differentiate only $D$” recipe would need to be extended through $B$; this limit is not tested in the paper.
- The 10-15% overhead figure is from one CPU tetrahedral benchmark; on GPUs or with reverse mode, the balance between time and memory may shift, and the method's advantage may become memory efficiency rather than speed.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes to localize automatic differentiation (AD) to the integration-point level within the finite element operator decomposition (FEOD), expressing the nonlinear finite element operator as Ap(u)=P^T G^T B^T D(BGPu) and its Jacobian as Jp(u)=P^T G^T B^T J_D(uq) BGP. The authors argue that this is computationally effective, scalable, automatic, and non-intrusive for existing parallel solvers, and they demonstrate it on a p-Laplacian problem in MFEM, comparing integration-point-level AD against element-level AD and a hand-coded implementation, with a topology optimization example as a design application.
Significance. The central chain-rule factorization is simple and correct under the stated assumption that the finite element operators P, G, and B are constant with respect to the differentiated quantities. The integration-point-level localization is a practical and potentially useful contribution for fixed-mesh problems, especially for high-order elements; the benchmark against a hand-coded implementation and the use of the open-source MFEM library are strengths. However, the experimental support is considerably narrower than the claims in the title and abstract: only one scalar problem, one mesh, a single process count, and no scaling study are presented, and the design component is only a single figure without quantitative details. The paper would be strengthened by explicitly scoping its claims and providing more evidence for scalability and design applicability.
major comments (3)
- [AUTOMATIC DIFFERENTIATION IN FINITE ELEMENT ANALYSIS, Eq. (2)] The factorization Jp(u)=P^T G^T B^T J_D(uq) BGP is exact only when P, G, and B are independent of the variable with respect to which one differentiates. The manuscript states that B does not depend on physical coordinates or design parameters, but this is false for isoparametric elements, moving meshes, ALE descriptions, or shape optimization, where B includes the inverse element-map Jacobian and therefore depends on nodal coordinates/design variables. In those cases, differentiating BGPu introduces an extra term (∂B/∂ρ)GPu, so the resulting Jacobian or sensitivity is not exact. The authors should either restrict the claims to fixed-reference-geometry problems with design parameters entering only through the pointwise operator D, or extend the formulation to differentiate through B.
- [Table 1 and the 'scalable' claim] The performance demonstration rests on a single benchmark: one p-Laplacian problem, one cube mesh with 200K elements, and 12 MPI processes, with no variance reported despite averaging over 100 runs. No scaling study is provided, so the abstract's claim of 'scalable' and 'large-scale' is not supported by the experimental data. Additionally, the statement that 'forward AD is only 10-15% slower than the highly optimized hand-coded implementation' is true only for the third-order element (Tet3: 2.86 s vs 2.53 s); for Tet1 and Tet2 the forward RES times are about 17% and 31% slower, respectively. The performance claims should be qualified accordingly.
- [Figure 2 and the 'design' contribution] The topology optimization example is presented as a single image with no accompanying problem formulation, objective and constraint definitions, optimization algorithm, convergence history, or performance measurements. Since 'Design' appears in the title and abstract, this is insufficient support for the design-related claims. The authors should either add quantitative details about the optimization setup and results or narrow the title and abstract to focus on analysis/Jacobian evaluation.
minor comments (5)
- [Introduction] There are typographical issues such as 'di fferentiation' and 'e fficient' with stray spaces; additionally, references [4] and [8] use 'at al.' and should be 'et al.'.
- [AUTOMATIC DIFFERENTIATION IN FINITE ELEMENT ANALYSIS] The sentence beginning 'Reverse AD mode introduces additional overhead for memory management and for relations with vector length uq equal to the output of the vector function D(uq) forward AD mode will be preferable' is grammatically incomplete and unclear; it should be rewritten.
- [Table 1] The units and meaning of 'Ke [s]' should be clarified: it is not clear whether the times are total assembly times for the whole mesh, average times per element, or something else. Reporting standard deviations or error bars would also help, given that the results are averaged over 100 runs.
- [References] Reference [11] (Hascoet and Pascual) appears in the reference list but is not cited anywhere in the text.
- [Equations (1)-(2)] The notation uq = BGPu is used in Eq. (2) but is defined only implicitly; it should be defined immediately after Eq. (1).
Circularity Check
No significant circularity; the AD-at-quadrature-point Jacobian is a direct chain-rule application, benchmarks are external, and self-citations are tool references rather than load-bearing premises.
full rationale
The paper's central derivation, Equations (1) and (2), is an ordinary chain-rule differentiation of the FEOD factorization Ap = P^T G^T B^T D(BGPu). No fitted parameter is renamed as a prediction, and no quantity is defined in terms of the result it is supposed to establish. The only nontrivial assumption is that P, G, and B are linear and independent of the solution, physical coordinates, and design parameters; this is stated explicitly as a hypothesis of the decomposition, not derived from the target result. The skeptical concern about isoparametric elements, moving meshes, and shape optimization is a validity limitation outside the fixed-reference-geometry regime, not a circularity, because the Jacobian formula is still exact whenever the stated assumption holds. The performance claims are benchmarked against an external hand-coded implementation and against independent AD tools (CoDiPack and Enzyme) in addition to MFEM's native dual numbers, so they are not internally forced. Self-citations to the MFEM FEOD work provide context and implementation machinery, but the conclusion that the proposed localization yields exact Jacobians and low overhead does not reduce to those citations; it rests on the displayed equations and the reported measurements. No circular step is present.
Assumptions & free parameters
assumptions (4)
- domain assumption The finite element operators P, G, and B are linear and independent of the solution, physical coordinates, and design parameters.
- domain assumption The pointwise operator D is differentiable with respect to quadrature-point values.
- standard math Automatic differentiation tools (CoDiPack, Enzyme, MFEM dual numbers) compute exact derivatives of the local operator.
- domain assumption Performance measurements via Caliper and PAPI on 12 MPI processes with 200K elements are representative of parallel scalability.
Cite this review
Pith. "Pith review of Scalable Analysis and Design Using Automatic Differentiation." pith.science (2026). https://pith.science/paper/OTYO5TQI
@misc{pith2026250600746,
author = {Pith},
title = {Pith review of: Scalable Analysis and Design Using Automatic Differentiation},
year = {2026},
howpublished = {\url{https://pith.science/paper/OTYO5TQI}},
note = {Machine review of arXiv:2506.00746}
}
read the original abstract
This article aims to demonstrate and discuss the applications of automatic differentiation (AD) for finding derivatives in PDE-constrained optimization problems and Jacobians in non-linear finite element analysis. The main idea is to localize the application of AD at the integration point level by combining it with the so-called Finite Element Operator Decomposition. The proposed methods are computationally effective, scalable, automatic, and non-intrusive, making them ideal for existing serial and parallel solvers and complex multiphysics applications. The performance is demonstrated on large-scale steady-state non-linear scalar problems. The chosen testbed, the MFEM library, is free and open-source finite element discretization library with proven scalability to thousands of parallel processes and state-of-the-art high-order discretization techniques.
Figures
Forward citations
Cited by 1 Pith paper
-
PDE-Constrained High-Order Mesh Optimization
A PDE-constrained optimization framework moves high-order mesh nodes to minimize a weighted sum of solution error and mesh distortion, cutting discretization error by up to 10x in Poisson and linear elasticity tests.
Reference graph
Works this paper leans on
-
[1]
A. Griewank and A. Walther, Evaluating derivatives : principles and techniques of algorithmic differentiation (SIAM, 2008)
work page 2008
-
[2]
Nørgaard et al., Structural and Multidisciplinary Optimization 56, 1135–1146 (2017)
S. Nørgaard et al., Structural and Multidisciplinary Optimization 56, 1135–1146 (2017)
work page 2017
-
[3]
R. Anderson et al., Computers & Mathematics with Applications 81, 42–74 (2021), development and Appli- cation of Open-source Software for Problems with Numerical PDEs
work page 2021
-
[4]
Andrej at al., The International Journal of High Performance Computing Applications p
J. Andrej at al., The International Journal of High Performance Computing Applications p. 10943420241261981 (2024), https://doi.org/10.1177/10943420241261981
-
[5]
Reverse-mode automatic di fferentiation and optimization of GPU kernels via enzyme,
W. Moses et al., “Reverse-mode automatic di fferentiation and optimization of GPU kernels via enzyme,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’21 (Association for Computing Machinery, New York, NY , USA, 2021)
work page 2021
-
[6]
Sagebaum et al., ACM Transactions on Mathematical Software (TOMS) 45, p
M. Sagebaum et al., ACM Transactions on Mathematical Software (TOMS) 45, p. 38 (2019)
work page 2019
-
[7]
I. Toulopoulos and T. Wick, SIAM Journal on Scientific Computing 39, A681–A710 (2017), https://doi.org/10.1137/16M1067792
-
[8]
Caliper: Performance introspection for HPC software stacks,
D. Boehme at al., “Caliper: Performance introspection for HPC software stacks,” in SC ’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (2016), pp. 550–560
work page 2016
Show all 11 references
-
[9]
Jagode, A
H. Jagode, A. Danalis, H. Anzt, and J. Dongarra, The International Journal of High Performance Computing Applications 33, 1113–1127 (2019)
2019
-
[10]
Duswald at al., Computer Methods in Applied Mechanics and Engineering 429, p
T. Duswald at al., Computer Methods in Applied Mechanics and Engineering 429, p. 117146 (2024)
2024
-
[11]
Hascoet and V
L. Hascoet and V . Pascual, ACM Trans. Math. Softw. 39may (2013), 10.1145/2450153.2450158
2013
Reviewed August 7, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.