Differentiable Parameter Optimization for DAEs with State-Dependent Events

Neural surrogates for stiff differential-algebraic equations (DAEs) face two barriers: soft-constraint methods leave algebraic residuals that stiffness amplifies into errors, and hard-constraint methods require trajectory data from stiff integrators. We introduce an extended Newton implicit layer that enforces algebraic constraints exactly and reduces fast dynamics to their quasi-steady-state values in a single differentiable solve. Embedded in a physics-informed DeepONet, the layer recovers all fast and algebraic states exactly from slow-state predictions, removes the per-window stiffness-amplification pathway, and yields a stiffness-scaled Implicit Function Theorem gradient absent from penalty methods. Cascaded implicit layers extend this to multi-component systems with provable convergence. On a grid-forming inverter (stiffness ratio of about 4712), extended Newton attains 1.42% error versus 39.3% (penalty) and 57.0% (standard Newton); augmented Lagrangian and feedback linearization diverged. Two independently trained models compose without retraining (0.72% to 1.16% error, exact constraint satisfaction). Cross-domain validation on the Robertson stiff DAE (stiffness ratio up to $10^5$) confirms generalization. Conformal prediction provides 90% coverage with automatic out-of-distribution detection.

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.

Solving massive-scale optimization problems requires scalable first-order methods with low per-iteration cost. This tutorial highlights a shift in optimization: using differentiable programming not only to execute algorithms but to learn how to design them. Modern frameworks such as PyTorch, TensorFlow, and JAX enable this paradigm through efficient automatic differentiation. Embedding first-order methods within these systems allows end-to-end training that improves convergence and solution quality. Guided by Fenchel-Rockafellar duality, the tutorial demonstrates how duality-informed iterative schemes such as the alternating direction method of multipliers, and the primal-dual hybrid gradient can be learned and adapted through representative case studies.