Pith. sign in

REVIEW 3 major objections 4 minor 17 references

Differentiate the Solver, Not the Equation: Reverse-Sweep Adjoints for Block Implicit Simulation

T0 review · 3 major / 4 minor · reviewed 2026-08-14 · deepseek-v4-flash

Pith's one-line read A reverse sweep of local 3x3 adjoint solves gives machine-precision gradients through the exact finite-depth Vertex Block Descent solver without forming any global system.

desk verdict Solid solver-level adjoint with machine-precision validation in the non-contact regime; the contact-scale experiment sits outside the theorem's hypotheses and needs an exactness check. read the letter →

arxiv 2608.08559 v1 pith:G52AZY2O submitted 2026-08-09 cs.GR cs.LGcs.NAmath.NA

classification cs.GRcs.LGcs.NAmath.NA
keywords solveradjointdifferentiationdifferentiableforwardimplicitsolvesautomatic
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Many physics simulators advance a deformable object by repeatedly solving small local problems: for each vertex, find the tiny displacement that lowers the energy, apply it, then move to the next vertex. To train a model or optimize a controller, a simulator must also provide gradients, telling how the final positions respond to small changes in material parameters, initial velocities, or forces. Two standard gradient routes exist. Unrolled automatic differentiation traces every arithmetic operation and backpropagates through the whole tape; it is exact but stores a huge graph. Equation-level implicit differentiation instead assumes the solver has converged and solves one global linear system; it uses little memory but returns the gradient of a converged solution, not of the short truncated solver that actually ran.

This paper takes a third route. It differentiates the executed solver itself, block by block. For any solver built from small implicit updates, the gradient is obtained by walking the updates in reverse order and, at each block, solving a tiny local adjoint problem. For Vertex Block Descent, that is a 3x3 solve per vertex, with the same graph coloring and parallel schedule as the forward pass. No global matrix is ever assembled.

The authors derive the exact conditions for this to work at finite solver depth and identify three terms that are easy to omit: the Jacobian of update damping, a Hessian-tangent term, and the derivatives of validity gates. Omitting any of them leaves gradients that look right at convergence but are wrong after few sweeps. They verify against autograd through the identical forward trajectory to machine precision, report large speed and memory gains, and demonstrate a one-GPU optimization over roughly 8 million vertices.

Extended reading notes

Core claim

The load-bearing claim is Theorem 1 and its instantiation: the exact discrete adjoint of an executed block implicit sweep is a reverse sweep of constant-size local adjoint solves, so the backward matches automatic differentiation through the identical executed forward to machine precision at every solver depth (Table 1: worst case 3.6e-15) while forming no global system. If correct, gradient-based learning, control, and inverse problems can use the solver structure itself as the backward pass, at a cost proportional to the forward and with memory linear in K times |V| at about 4 KB per vertex per sweep.

Load-bearing premise

The construction is exact only if the backward replays the same local linearization states as the forward: 'Provided the backward replays the same local linearization states as the forward, the adjoint of one sweep is a reverse block sweep of constant-size local solves' (Section 3.4). In production code this requires every non-smooth safeguard (saturation damping, regularization floors, validity gates, contact set changes) to be reconstructible and differentiable; the paper differentiates saturation and gates, but if a forward branch such as broad-phase pair changes or active-set contact cannot be replayed identically, the machine-precision claim fails.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The paper proposes differentiating the finite executed block-implicit solver rather than the converged equation, instantiating the idea on Vertex Block Descent. The central claim (Theorem 1) is that the exact discrete adjoint of a finite composition of block implicit updates is a reverse sweep of constant-size local adjoint solves, with no global linear-algebra object. The authors prove the theorem under C1-differentiability and disjoint same-color stencil assumptions, give closed-form local adjoints including safeguard terms (saturation Jacobian, Hessian tangent, validity gates), and validate the construction against unrolled automatic differentiation through the identical executed forward at solver depths K=1, 5, 20 on non-contact scenes, reporting machine-precision agreement. They also report ablations showing the error laws of omitted safeguard terms, a controlled cost comparison against unrolled AD and global equation-level adjoints, exactness transfers to XPBD and projective dynamics, inverse-problem results, and a large-scale demonstration with 10^6 contact-coupled bodies.

Significance. If the central claim is correct, the contribution is significant: it replaces a global adjoint system with a block-local reverse sweep that inherits the forward solver's parallelism and memory behavior, and it makes finite-solver-depth exactness precise rather than asymptotic. The paper's strengths include a self-contained proof of the core theorem, machine-precision validation against an external autograd reference on the same executed forward, explicit identification and measurement of non-smooth safeguard terms that prior hand-derived adjoints omit, a same-codebase controlled cost comparison, and a reproducible harness with single-command experiments. The main weakness is that the contact-coupled experiments, including the headline 10^6-body scale demo, sit outside the theorem's hypotheses and are not validated against an autograd reference, so the exactness claim is not yet established for the regime the paper itself presents as its largest demonstration.

major comments (3)
  1. [§3.4, §5, Appendix A] The contact terms of §5 fall outside the hypotheses of Theorem 1. The proof in Appendix A requires each F_i to be C1 in a neighborhood of the executed linearization state and same-color updates to touch disjoint stencils. The IPC log barrier b(d) = -kappa (d-dhat)^2 log(d/dhat) has a third-derivative discontinuity at the activation boundary d = dhat, and the active set can change across sweeps; the Hessian-tangent term of Eq. (10) needs the third derivative of the barrier, so the C1 condition is not met in a neighborhood of the boundary. In addition, the center-of-mass penalty between separate bodies couples vertices through body centers and broad-phase pair lists, so two same-color vertices can have overlapping stencils through the shared body center; the disjoint-stencil step of Appendix A then does not apply. This is load-bearing because §6.6 uses exactly these mechanisms. Please either prove an extended exactness statement for these terms, or provide a small-scale autograd-vs-reverse-sweep validation with contact activation/deactivation and changing broad-phase pairs at multiple K, and restrict the exactness claims accordingly.
  2. [§6.6, Table 1, Appendix B] The scale experiment in §6.6 reports only task success (loss reduction to 87%) and never reports a gradient error against an autograd reference through the identical executed forward. Table 1 covers meshes up to |V|=631 with no contact, and Appendix B covers XPBD and projective dynamics on spring systems with no contact. Since the contact-coupled center-of-mass penalty and the IPC barrier are claimed to be differentiated by the same reverse sweep, the finite-depth exactness of the contact path is unverified. This matters precisely because the paper's thesis is that at small K the reverse sweep differs from equation-level differentiation; without an autograd check on a contact scene, there is no evidence that the contact backward pass is the exact adjoint of the executed forward. Please add such a check, including scenes where the contact active set changes between sweeps or steps.
  3. [§5, §H (Contact)] The statement that 'pair indices carry no gradient' does not by itself resolve the branch-replay problem. The broad-phase pair list is a discrete function of the current positions; if the reverse sweep does not store or exactly replay which pairs were active in each forward sweep, the backward pass may not correspond to the function the forward actually executed. The replay buffer stores only the applied update and the six symmetric Hessian entries per vertex, not the pair list or activation flags. The paper should specify how discrete branch decisions (active IPC contact pairs, broad-phase candidate pairs, saturation/regularization branch selections) are replayed exactly, or store them, and verify this replay in a contact experiment. Without this, the 'provided the backward replays the same local linearization states' qualifier in §3.4 is not satisfied for the contact machinery.
minor comments (4)
  1. [Abstract, §6.1] The abstract states that the backward matches automatic differentiation 'to machine precision at every solver depth'; the measured support in Table 1 is for non-contact scenes with |V| up to 631. Please qualify the abstract and §6.1 to indicate the validated regime, given that the contact-coupled experiments are not covered by the machine-precision measurements.
  2. [Appendix H] The text says the reverse pass 'undoes Δx_i to recover the pre-update state'; in floating-point arithmetic, x := x + Δx followed by x := x - Δx is not guaranteed to return the original x bit-exactly. The paper reports machine-precision agreement, so this apparently does not affect the measurements, but the mechanism should be explained or the reconstruction should be specified (e.g., storing the pre-update value or showing the error is below the reported tolerance).
  3. [§6.4] The cross-paper comparison with DiffIPC ('two orders of magnitude less memory') is not controlled and the caveat is stated; consider moving this comparison to a clearly labeled 'published numbers' paragraph and referring to Table 3 as the only controlled memory comparison.
  4. [§3.5 (iii)] In the description of validity gates, the quantity '8.2×10^-3' is reported without units; please state whether this is an absolute or relative gradient error and in what norm.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: the reverse-sweep adjoint is a direct chain-rule/implicit-function identity, validated against an external autograd reference on the same executed forward.

full rationale

The load-bearing result, Theorem 1, is a mathematical identity rather than a fitted or self-referential claim. Equation (5) states that the transpose of a composition is the reverse composition of transposes, and the local vector-Jacobian products of Eqs. (7)-(10) follow from differentiating the local optimality condition Fi = 0 and applying the implicit function theorem under the stated C1 and invertibility assumptions. No fitted constants, exclusions, or post hoc adjustments enter this derivation. The machine-precision validation is against PyTorch autograd through the identical executed forward, with forward trajectories verified identical to 1e-16 before comparing gradients, so the reference is external to the paper's own derived predictions. The ablations in Section 3.5 are explicitly omissions from the authors' own implementation, not redefinitions of the target result. The paper builds on Vertex Block Descent from Chen et al. (2024a), which includes a co-author of the present paper, but VBD is used as an external forward solver to instantiate the construction, not as evidence for the adjoint identity; Theorem 1 applies to any block implicit solver satisfying its hypotheses and is proved in Appendix A from the chain rule alone. The contact and scale experiments that may fall outside the theorem's C1 or disjoint-stencil conditions are acknowledged limitations and scoping statements, and any correctness concern there is a matter of scope or rigor, not circularity. Accordingly, no specific reduction of the derivation to its own inputs can be exhibited, and the appropriate finding is no significant circularity.

Assumptions & free parameters 3 free parameters · 5 assumptions · 0 invented entities

The method uses no fitted constants in its derivation and introduces no new physical or mathematical entities. It relies on standard chain rule and implicit function theorem results, plus domain assumptions about block-solver structure and replay fidelity. Implementation constants such as damping c and regularization tau are part of the executed forward and are differentiated, not tuned to produce the headline result.

free parameters (3)
  • saturation damping constant c = 0.1 m
    Forward safeguard used in the executed solver (App. H). The adjoint differentiates through it; the numerical value is a simulation choice, not fitted to the exactness result.
  • Hessian regularization tau = 1e-6 in h2-scaled system
    Added to the per-vertex block (App. H). Differentiated rather than treated as constant; the value does not influence the central exactness claim.
  • IPC barrier stiffness kappa and activation distance dhat = kappa=400, dhat=2 cm in Cover the Spot; other tasks vary
    Contact model constants for experiments. They set the contact force but are not fitted to the adjoint-exactness claim.
assumptions (5)
  • standard math Chain rule for composition: (dS)^T = (dS1)^T ... (dSn)^T for S = Sn composed with ... composed with S1
    Invoked in Section 3.2, Eq. (5), and Appendix A; forms the basis of the reverse-sweep construction.
  • standard math Implicit function theorem applies to each local block solve with invertible local system
    Theorem 1 and Eq. (7), Appendix A; requires C1 Fi and invertible dFi/dDelta_xi, or for the affine instance, invertible Hi.
  • domain assumption Distance-1 graph coloring makes same-color vertex updates touch disjoint stencils, so the parallel color update equals composition in any order
    Section 3.4 and Appendix A; needed for the reverse sweep over colors to be exact. If two same-color blocks shared an element, simultaneous and sequential updates would differ.
  • domain assumption The stored 72-byte replay buffer (Delta_x and symmetric H per vertex per sweep) is sufficient to reconstruct all forward linearization states and safeguard Jacobians
    Section 4 and App. H; the strongest operational premise. Verified empirically to 1e-16 but not proven for arbitrary data-dependent branches.
  • domain assumption All forward safeguards are differentiable and their Jacobians are included in the backward pass
    Section 3.5; identifies saturation, Hessian tangent, and validity gates. If a safeguard is non-differentiable or omitted, finite-depth exactness fails.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Differentiate the Solver, Not the Equation: Reverse-Sweep Adjoints for Block Implicit Simulation." pith.science (2026). https://pith.science/paper/G52AZY2O

@misc{pith2026260808559,
  author       = {Pith},
  title        = {Pith review of: Differentiate the Solver, Not the Equation: Reverse-Sweep Adjoints for Block Implicit Simulation},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/G52AZY2O}},
  note         = {Machine review of arXiv:2608.08559}
}
abstract

Differentiable simulation is a key component in learning, control, and inverse problems, where gradients through nonlinear implicit solvers are required. Existing approaches either rely on unrolled automatic differentiation, whose memory grows with solver depth, or on equation-level implicit differentiation, which assembles global Jacobians and solves large sparse adjoint systems, discarding the locality of the forward solver -- and differentiating the converged equation rather than the finite computation that actually ran. We propose solver-level differentiation, which differentiates the executed solver itself. When a solver is composed of block implicit updates, its discrete adjoint is obtained by applying the corresponding adjoint updates in reverse order, yielding a reverse-sweep formulation whose backward pass mirrors the forward solver. From an operator perspective, the forward pass realizes an approximate inverse through ordered local solves, and the backward applies its transpose through reverse local adjoint solves, constructing no global system. We instantiate this idea on Vertex Block Descent, yielding a differentiable solver whose reverse colored Gauss-Seidel sweeps are composed entirely of local $3\times 3$ adjoint solves. The backward matches automatic differentiation through the identical executed forward to machine precision at every solver depth, where the equation-level adjoint is off by 37% after one sweep; in a controlled same-codebase, same-GPU comparison it is 33x faster and uses 71x less memory than unrolled automatic differentiation; and the same construction is exact on projective dynamics and extended position-based dynamics. We scale differentiable elastodynamics to $10^6$ contact-coupled soft bodies (8M vertices) on one GPU. Overall, this work highlights solver structure as a practical organizing principle for efficient differentiable simulation.

Figures

Figures reproduced from arXiv: 2608.08559 by the authors.

Figure 1
Figure 1. Vertex Block Descent on a tetrahedral mesh. Left: distance-1 graph coloring; vertices that share an element never share a color, so each color class updates in parallel. Right: one local 3×3 Newton solve, touching only the elements incident to the updated vertex. 2 BACKGROUND AND NOTATION Notation. Bold italic denotes vectors (x, v, y, gi , pi , ∆xi); roman capitals denote matrices (M, Hi , I); plain italic denotes … view at source ↗
Figure 2
Figure 2. Data flow around one block update. The forward advances the block state; the reverse consumes the incoming adjoint x¯ + i (accumulated by reverse updates already executed) and produces the outgoing adjoint x¯ − i . Superscripts denote position in the sweep, not before/after in time. of its executed linearization state and that ∂Fi/∂∆xi is invertible there; for the affine instance Fi = Hi∆xi + gi it suffices that Hi … view at source ↗
Figure 3
Figure 3. Cube center of mass trajectory under the initial parameters versus the optimized ones, side and top view. The optimizer, not the initialization, drives the target hit. Shear modulus from a scalar sag target. A 484-vertex square cloth (side 0.6 m, extruded to a thin three-tet shell) is pinned at its four corners and sags into a hammock under gravity. Stiffer µ gives shallower center sag, so the inverse problem L(µ)=(… view at source ↗
Figures from the paper (10 more)
Figure 4
Figure 4. Figure 4: Beam µ/λ identification, target deformation (µ, λ)=(800, 400) under a 60◦ angular kick. Blue ghost: target equilibrium. Both parameters are recovered to < 0.3% [PITH_FULL_IMAGE:figures/full_fig_p015_4.png]
Figure 5
Figure 5. Figure 5: Dress µ/λ identification. Left: initial stiff configuration (µ=5000). Remaining panels: target draping progression (µ=1500). error is below 10−6 and stable across the ϵ sweep (so that it reflects the gradient, not FD truncation or round-off). All 11 checks pass; [PITH…
Figure 6
Figure 6. Figure 6: Bunny initial-state optimization: start, mid-bounce, and at the target [PITH_FULL_IMAGE:figures/full_fig_p016_6.png]
Figure 7
Figure 7. Figure 7: Cover the Spot. Three soft-body disks are jointly optimized by backpropagation through the full DiffVBD rollout to land such that their union covers a painted target spot. Top: top-down view of the optimized landing (left) and optimization curves (right). Bottom: 3D sn…
Figure 12
Figure 12. Figure 12: Beam µ identification at three twist stages. Blue: initial (µinit=105 ). Coral: target (µtrue=106 ). Green: recov￾ered. Material identification. Given a target rollout produced by unknown Lame parameters, recover them. On a ´ 1029- vertex cantilever twisted 180◦ and r…
Figure 8
Figure 8. Figure 8: Shear modulus recovery on a four-corner-pinned cape. [PITH_FULL_IMAGE:figures/full_fig_p018_8.png]
Figure 9
Figure 9. Figure 9: Fitting a procedural dress to a body by optimizing a rest-shape scale. [PITH_FULL_IMAGE:figures/full_fig_p018_9.png]
Figure 14
Figure 14. Figure 14: Flag under the recovered time￾varying wind, at three times in the rollout. Wind control. Given a target flag trajectory, recover a spatially-uniform, time-varying wind on a 442-vertex flag pinned along its top edge, parameterized as piecewise￾constant over 8 equal tim…
Figure 10
Figure 10. Figure 10: Wind-force recovery on a top-pinned cape. [PITH_FULL_IMAGE:figures/full_fig_p019_10.png]
Figure 11
Figure 11. Figure 11: Optimized cube center of mass trajectories at three solver depths under identical 40-iteration Adam budgets, with the uncontrolled trajectory under the initial parameters overlaid (gray dashed). Every optimized trajectory passes close to the target; the ordering betwe…

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

17 extracted references · 7 canonical work pages

  1. [8]

    Yifei Li, Tao Du, Kui Wu, Jie Xu, and Wojciech Matusik

    doi: 10.1145/3386569.3392425. Yifei Li, Tao Du, Kui Wu, Jie Xu, and Wojciech Matusik. DiffCloth: Differentiable cloth simulation with dry frictional contact.ACM Transactions on Graphics, 42(1):2:1–2:20,

  2. [16]

    Both are float64 end to end

    kernels; the second is a fused CUDA C++ extension compiled on first use through PyTorch’s load inline and cached thereafter. Both are float64 end to end. They agree to 10−16 on forward state and on gradients, which is what lets us use either as the reference for the other; the fused path is 7–22× faster because one kernel launch per color performs assembl...

  3. [17]

    I INVERSE-PROBLEM TASKS IN DETAIL These are the per-task descriptions for §6.5; the headline numbers are Table 4 in the main text, and the reproduction protocols are Appendix J

    every5 steps, bodies are packed into a single mega-mesh processed in chunks of 65,536, gradients are clipped at norm 5, and the material is optimized through a sigmoid parameterization that confines µ to a bounded range instead of relying on projection. I INVERSE-PROBLEM TASKS IN DETAIL These are the per-task descriptions for §6.5; the headline numbers ar...

  4. [1964]

    doi: 10.1214/aoms/1177703732. Krishna Murthy Jatavallabhula, Miles Macklin, Florian Golemo, Vikram V oleti, Linda Petrini, Martin Weiss, Breandan Considine, J ´erˆome Parent-L´evesque, Kevin Xie, Kenny Erleben, Liam Paull, Florian Shkurti, Derek Nowrouzezahrai, and Sanja Fidler. gradSim: Differentiable simulation for system identification and visuomotor c...

  5. [1990]

    Keenon Werling, Dalton Omens, Jeongseok Lee, Ioannis Exarchos, and C

    doi: 10.1109/5.58337. Keenon Werling, Dalton Omens, Jeongseok Lee, Ioannis Exarchos, and C. Karen Liu. Fast and feature-complete differentiable physics for articulated rigid bodies with contact. InRobotics: Science and Systems (RSS),

  6. [2004]

    Daniel Freeman, Samuel S

    Luke Metz, C. Daniel Freeman, Samuel S. Schoenholz, and Tal Kachman. Gradients are not all you need.arXiv preprint arXiv:2111.05803,

  7. [2007]

    jvcir.2007.01.005

    doi: 10.1016/j. jvcir.2007.01.005. Georgios Pavlakos, Vasileios Choutas, Nima Ghorbani, Timo Bolkart, Ahmed A. A. Osman, Dimitrios Tzionas, and Michael J. Black. Expressive body capture: 3D hands, face, and body from a single image. InIEEE Conference on Computer Vision and Pattern Recognition (CVPR),

  8. [2016]

    Yixin Chen, David I. W. Levin, and Timothy R. Langlois. Fluid control with Laplacian eigenfunctions. InACM SIGGRAPH 2024 Conference Papers, 2024b. Jonas Degrave, Michiel Hermans, Joni Dambre, and Francis Wyffels. A differentiable physics engine for deep learning in robotics.Frontiers in Neurorobotics, 13:6,

Show all 17 references
  1. [2018]

    Tuur Stuyck and Hsiao-yu Chen

    doi: 10.1145/3180491. Tuur Stuyck and Hsiao-yu Chen. DiffXPBD: Differentiable position-based simulation of compliant constraint dynamics.Proceedings of the ACM on Computer Graphics and Interactive Techniques, 6(3):51:1–51:14,

  2. [2019]

    neckline (y=1.45 m) and drapes through 35 DiffVBD steps. A single scalarscale s multiplies all rest-pose radii, so the optimized quantity enters through the rest metric rather than through the material or the initial state; the target is the hem heighty ⋆ hem=0.654m under L(s)...

  3. [2020]

    10 Michael B

    doi: 10.1145/3414685.3417766. 10 Michael B. Giles. On the use of Runge-Kutta time-marching and multigrid for the solution of steady adjoint equations. Technical Report NA-00/10, Oxford University Computing Laboratory,

  4. [2021]

    For a single block, the C 1 hypothesis on Fi together with invertibility of ∂Fi/∂∆xi at the executed linearization state gives Eq

    A PROOF OFTHEOREM1 Equation (5) is the chain rule applied to S=S n ◦ · · · ◦S1, so the adjoint visits blocks in reverse order. For a single block, the C 1 hypothesis on Fi together with invertibility of ∂Fi/∂∆xi at the executed linearization state gives Eq. (7) by the implicit...

  5. [2022]

    Moritz Geilinger, David Hahn, Jonas Zehnder, Moritz B¨acher, Bernhard Thomaszewski, and Stelian Coros

    doi: 10.1145/3490168. Moritz Geilinger, David Hahn, Jonas Zehnder, Moritz B¨acher, Bernhard Thomaszewski, and Stelian Coros. ADD: Analytically differentiable dynamics for multi-body systems with frictional contact. ACM Transactions on Graphics, 39(6):190:1–190:15,

  6. [2023]

    Hyung Ju Terry Suh, Max Simchowitz, Kaiqing Zhang, and Russ Tedrake

    doi: 10.1145/3606923. Hyung Ju Terry Suh, Max Simchowitz, Kaiqing Zhang, and Russ Tedrake. Do differentiable simulators give better policy gradients? InInternational Conference on Machine Learning (ICML), pp. 20668–20696,

  7. [2024]

    doi: 10.1145/3657648. Peter J. Huber. Robust estimation of a location parameter.The Annals of Mathematical Statistics, 35 (1):73–101,

  8. [2025]

    doi: 10.2312/pg.20251265

    ISBN 978-3-03868-295-0. doi: 10.2312/pg.20251265. Anka He Chen, Ziheng Liu, Yin Yang, and Cem Yuksel. Vertex block descent.ACM Transactions on Graphics, 43(4):116:1–116:16, 2024a. doi: 10.1145/3658179. Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Ne...

  9. [2026]

    Training deep nets with sublinear memory cost.arXiv preprint arXiv:1604.06174,

    Tianqi Chen, Bing Xu, Chiyuan Zhang, and Carlos Guestrin. Training deep nets with sublinear memory cost.arXiv preprint arXiv:1604.06174,

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.