DiffPhD: Differentiable Elastodynamics Solver

• DiffPhD is a differentiable elastodynamics solver that uses stiffness-aware projective assembly to reliably handle extreme stiffness contrasts and non-convex hyperelastic behavior.
• The framework reuses a single GPU-accelerated sparse factorization for forward simulation, contact computations, and backward adjoint solves, driving significant speedups.
• It integrates backward-pass trust-region eigenvalue filtering to stabilize hyperelastic gradient calculations, ensuring accurate and convergent differentiation in contact-rich settings.

DiffPhD is a differentiable elastodynamics solver built around Projective Dynamics (PD), formulated as a unified GPU-accelerated differentiable Projective Dynamics framework for heterogeneous materials under large deformation and contact-rich interaction regimes (Lai et al., 14 May 2026). It is designed for settings in which existing differentiable PD methods, particularly differentiable Projective Dynamics approaches such as DiffPD, become unreliable: heterogeneous materials with extreme stiffness contrasts, compressible Neo-Hookean hyperelasticity, and frictional contact. The framework’s central organizing principle is to preserve the core PD advantage of reusing a single sparse factorization while extending that reuse to forward solves, contact computations, and backward adjoint computations within one forward–backward pass (Lai et al., 14 May 2026).

1. Problem setting and design objective

DiffPhD targets differentiable simulation of soft bodies for system identification, trajectory optimization, and Real2Sim transfer (Lai et al., 14 May 2026). The motivating claim is that three failure modes—material heterogeneity, hyperelastic non-convexity, and contact—are coupled rather than separable. A solver may remain stable on homogeneous materials, or on hyperelasticity alone, or on contact alone, yet fail when all three appear simultaneously (Lai et al., 14 May 2026).

The implicit Euler elastodynamics step is written as

$$\frac{\bm{M}}{h^{2}}(\bm{q}_{t+h}-\tilde{\bm{q}}) = \bm{f}_{\mathrm{ela}}(\bm{q}_{t+h}) - \bigl(\alpha\bm{M}+\bm{B}_{\bm{\beta}}\bigr)\bm{v}_{t+h} +\bm{J}_n^{T}\bm{\lambda}_n+\bm{J}_b^{T}\bm{\lambda}_b+\bm{J}_f^{T}\bm{\lambda}_f,$$

with free-fall target

$$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$

The elastic constitutive model is the heterogeneous compressible Neo-Hookean energy

$$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$

where

$$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$

This model captures large deformation and inversion resistance, but the paper emphasizes that it also introduces non-convexity and indefinite Hessians, which destabilize prior differentiable PD pipelines (Lai et al., 14 May 2026). DiffPhD’s stated objective is therefore not merely faster simulation, but strict gradient accuracy together with convergence under stiffness contrasts up to $100\times$, while retaining GPU efficiency through persistent sparse operators (Lai et al., 14 May 2026).

2. Heterogeneous Projective Dynamics formulation

The paper identifies the conditioning of the PD global operator as the central bottleneck under heterogeneity. In classic PD or DiffPD, the global matrix is

$$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$

When stiffness varies by $10\times$, $50\times$, or $100\times$, this operator becomes badly conditioned if stiffness variation is not embedded into the operator itself (Lai et al., 14 May 2026).

DiffPhD’s response is stiffness-aware projective assembly. Rather than placing per-element heterogeneity in the local Neo-Hookean proximal map, the framework routes heterogeneity into the global PD weights,

$$w_e=\bar{k}(\mu_e,\lambda_e,\varepsilon_\sigma)\propto \mu_e.$$

The assembled global matrix becomes

$$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$0

This separation is explicit in the method description. The local proximal operator uses mesh-wide scalar means $$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$1, while heterogeneity is encoded in the global system through $$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$2 (Lai et al., 14 May 2026). The stated rationale is that if the local proximal map varies by element, the fixed-point map becomes inconsistent across the mesh; soft elements can enter indefinite regimes and Anderson Acceleration can amplify the resulting oscillations (Lai et al., 14 May 2026). A plausible implication is that DiffPhD treats heterogeneity primarily as a global linear-algebraic conditioning problem rather than a purely local constitutive problem.

For each element, the Neo-Hookean local step is a proximal problem,

$$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$3

Using the singular value decomposition $$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$4, the minimization decouples into stretch space:

$$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$5

The stretch-space Hessian is

$$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$6

The converged projection is

$$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$7

The paper also describes an optional logarithmic volume-barrier variant,

$$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$8

used when inversion avoidance is paramount (Lai et al., 14 May 2026).

3. Forward iteration, persistent sparse factors, and GPU execution

The core systems concept in DiffPhD is that one sparse symmetric positive definite PD operator should serve three roles: the forward global PD solve, the contact Delassus operator, and the backward adjoint solve (Lai et al., 14 May 2026). To realize this, the method computes a single sparse factorization using nested dissection:

$$\tilde{\bm{q}} = \bm{q}_t+h\bm{v}_t +h^{2}\bm{M}^{-1}\bigl(\bm{f}_{\mathrm{ext}}+\bm{f}_{\mathrm{state}}\bigr).$$9

Thus,

$$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$0

The factor pair $$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$1 is stored in persistent GPU buffers and reused everywhere unless topology, material parameters, or damping change (Lai et al., 14 May 2026).

Forward fixed-point iteration is accelerated by type-II Anderson Acceleration (AA), but the method deliberately bounds the history window because long AA histories are reported to be unstable on heterogeneous meshes. The residual and history matrices are

$$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$2

with $$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$3. The mixing coefficients are

$$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$4

and the extrapolated iterate is

$$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$5

Two safeguards are specified. If $$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$6, the AA history is discarded and the solver falls back to the non-accelerated iterate; and the history window is reduced to $$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$7 on heterogeneous meshes, versus $$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$8 on homogeneous ones (Lai et al., 14 May 2026).

Convergence is checked by a dual-gate criterion,

$$\psi_e(\bm{F}_e) = \frac{\mu_e}{2}\bigl(\|\bm{F}_e\|_F^{2}-d\bigr) -\mu_e\ln J_e +\frac{\lambda_e}{2}(\ln J_e)^2,$$9

The paper motivates this by noting that a step-size-only criterion can falsely declare convergence when the first step is small or when contact state changes produce a nearly stationary iterate that is not actually a fixed point (Lai et al., 14 May 2026).

4. Backward-pass stabilization and exact differentiability

A principal technical contribution of DiffPhD is backward-pass trust-region eigenvalue filtering lifted to the prox-map Hessian (Lai et al., 14 May 2026). In differentiable PD, the local proximal operator must be differentiated through the implicit function theorem. For Neo-Hookean materials, the relevant matrix is

$$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$0

which may be indefinite.

The paper states that naively inverting this matrix can destroy descent directions in the adjoint computation and render gradients unreliable (Lai et al., 14 May 2026). DiffPhD therefore introduces the filtered family

$$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$1

The discrete choices are interpreted as $$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$2 for the unfiltered inverse, $$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$3 for clamping-like behavior, and $$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$4 for absolute-value filtering (Lai et al., 14 May 2026).

Selection is state-adaptive. Using the last converged PD increment

$$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$5

the actual and model decreases are

$$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$6

giving the trust-region ratio

$$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$7

with $$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$8.

Once filtered, the local differential satisfies

$$\mu_e=\frac{E_e}{2(1+\nu)},\qquad \lambda_e=\frac{E_e\nu}{(1+\nu)(1-2\nu)}.$$9

This yields $100\times$0, which is inserted into

$100\times$1

The paper presents this as the exact location at which stable hyperelastic gradients are obtained (Lai et al., 14 May 2026). This suggests that DiffPhD’s differentiability is not an automatic consequence of differentiable primitives, but a carefully regularized implicit differentiation scheme tailored to the structure of the Neo-Hookean prox map.

With contact inactive, the backward solve reduces to one GPU two-SpMV application through the persistent factor (Lai et al., 14 May 2026). With contact active, the framework solves a block KKT adjoint system whose $100\times$2-block uses the same sparse factor pair and whose $100\times$3-block uses dense LDLT on a reduced contact system (Lai et al., 14 May 2026). Gradients are propagated to states, forces, and material parameters, including projective weights,

$100\times$4

and then to $100\times$5 and $100\times$6 by the chain rule (Lai et al., 14 May 2026).

5. Contact, damping, and unified reduced systems

DiffPhD incorporates contact through the Signorini–Coulomb complementarity formulation adopted from FBA (Lai et al., 14 May 2026). For each contact, the normal complementarity residual uses a Fischer–Burmeister map,

$100\times$7

with

$100\times$8

where the Delassus operator is

$100\times$9

The reduced linear contact system is

$$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$0

with

$$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$1

A notable systems detail is that the Delassus matrix is assembled via batched sparse-dense multiplications,

$$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$2

again reusing the persistent factor (Lai et al., 14 May 2026). Position updates reuse cached columns $$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$3:

$$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$4

The same strategy governs damping. DiffPhD uses Rayleigh damping

$$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$5

which folds directly into the global operator

$$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$6

The elementwise damping coefficient is chosen proportional to stiffness,

$$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$7

The paper explicitly notes that this is opposite to the more standard strategy of damping soft regions more heavily (Lai et al., 14 May 2026). Its stated argument is that high-frequency ringing in heterogeneous twist or impact scenarios is often seeded by stiff regions at the interface, so damping the stiff source suppresses ringing without overdamping the soft, visually important motion (Lai et al., 14 May 2026). This suggests a heterogeneity-aware dissipation model in which damping is aligned with stiffness-generated spectral pollution rather than with deformation amplitude alone.

6. Benchmarks, empirical claims, and scope

DiffPhD is evaluated on heterogeneous forward simulation, contact-rich forward simulation, inverse problems, and Real2Sim or robotics scenarios on a single RTX 4090 GPU plus an i7-12700 CPU (Lai et al., 14 May 2026). The reported benchmarks include Cantilever, Armadillo twist, Crab, Gatorman, Napkin, Bouncing Ball, Plant, Bunny, Routing Tendon, Torus, Oreo, and Dice (Lai et al., 14 May 2026).

Before the summary table, two empirical claims organize the paper’s interpretation. First, DiffPhD is reported to remain convergent on stiffness contrasts up to $$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$8, whereas prior PD solvers degrade in that regime (Lai et al., 14 May 2026). Second, the framework reports strict gradient accuracy together with up to an order-of-magnitude speedup over prior differentiable solvers, with some benchmark-specific forward or backward speedups exceeding that scale (Lai et al., 14 May 2026).

Benchmark group	Examples	Reported outcomes
Heterogeneous forward simulation	Cantilever, Armadillo twist, Crab	Convergent up to $$\bm{A} = \frac{\bm{M}}{h^{2}}+\sum_e w_e V_e \bm{G}_e^T\bm{G}_e.$$9 stiffness contrast; on Crab, up to about $10\times$0 forward speedup over DiffPD and about $10\times$1 over MAS
Contact-rich simulation	Gatorman, Napkin	On heterogeneous Napkin at $10\times$2, up to about $10\times$3 backward speedup over DiffPD
Inverse problems and Real2Sim	Plant, Bunny, Routing Tendon, Torus, Dice	Plant loss $10\times$4 vs. $10\times$5 for DiffPD/MAS; Bunny heterogeneous loss $10\times$6 vs. $10\times$7 and $10\times$8; Routing Tendon $10\times$9 while DiffPD and MAS stall around $50\times$0; Dice converges in 13 L-BFGS evaluations in roughly 13.5 minutes

The inverse-problem section is particularly important because it ties solver architecture to optimization quality rather than only wall-clock performance. Reported final losses include Plant at $50\times$1 for DiffPhD versus $50\times$2 for DiffPD and MAS, Bunny (heterogeneous) at $50\times$3 versus $50\times$4 and $50\times$5, and Routing Tendon at $50\times$6 while DiffPD and MAS stall around $50\times$7 (Lai et al., 14 May 2026). On Torus, the objective is matched or slightly improved, with a reported $50\times$8 speedup over MAS, though the paper notes that the benefit is smaller because contact itself dominates (Lai et al., 14 May 2026).

For Real2Sim and manipulation, the Dice experiment recovers probe contact location from real PokeFlex data and is reported to converge in 13 L-BFGS evaluations in roughly 13.5 minutes wall-clock (Lai et al., 14 May 2026). The Oreo example is presented as a forward-only capability check with extreme stiffness contrast across a multi-body gripper or object system (Lai et al., 14 May 2026).

The ablation study separates projection choice from GPU residency. Its conclusion is that projection choice governs gradient quality, while the unified GPU pipeline provides the speedup (Lai et al., 14 May 2026). This distinction is significant because it frames DiffPhD not as a single numerical trick, but as a layered design: heterogeneous PD assembly for conditioning, trust-region prox differentiation for gradient stability, and persistent sparse-factor reuse for throughput.

A common misconception would be to regard DiffPhD as merely a GPU port of Projective Dynamics. The paper does not support that reading. Its main contributions are stated as stiffness-aware projective assembly, backward-pass trust-region filtering for hyperelastic gradients, unified sparse-factor reuse across forward, backward, and contact, and stabilized forward iteration with bounded Anderson Acceleration and dual-gate convergence (Lai et al., 14 May 2026). Another misconception would be to treat contact, damping, and differentiation as add-on modules. The framework instead folds damping into the same operator, constructs contact through the same inverse applications, and solves adjoints through the same factorization (Lai et al., 14 May 2026).

Taken together, DiffPhD can be understood as a reformulation of differentiable Projective Dynamics in which the global operator becomes the central computational object: heterogeneity-aware, reused across all phases, and stabilized sufficiently to support end-to-end gradient-based optimization in shell–joint composite creatures, soft characters wielding stiff weapons, soft-gripper robotic manipulation, and related contact-rich regimes (Lai et al., 14 May 2026).