---
title: 'DiffPhD: Differentiable Elastodynamics Solver'
url: https://www.emergentmind.com/topics/diffphd
type: topic
---

# DiffPhD: Differentiable Elastodynamics Solver

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 [2605.14526]. 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 [2605.14526].

## 1. Problem setting and design objective

DiffPhD targets differentiable simulation of soft bodies for system identification, trajectory optimization, and Real2Sim transfer [2605.14526]. 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 [2605.14526].

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 [2605.14526]. 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 [2605.14526].

## 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 [2605.14526].

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
$$
\bm{A} = \frac{(1+\alpha h)\bm{M}}{h^{2}} + \sum_e \Bigl(w_e+\frac{\beta_e}{h}\Bigr)V_e \bm{G}_e^T\bm{G}_e.
$$

This separation is explicit in the method description. The local proximal operator uses mesh-wide scalar means \((\bar\mu,\bar\lambda,\bar k)\), while heterogeneity is encoded in the global system through \(w_e\) [2605.14526]. 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 [2605.14526]. 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,
$$
\bm{p}_{e}^{*} = \arg\min_{\bm{p}\in\mathbb{R}^{d\times d}} \frac{\bar{k}}{2}\|\bm{F}_e-\bm{p}\|_F^2 + V_e\psi_e(\bm{p}).
$$
Using the singular value decomposition \(\bm{F}_e=\bm{U}\operatorname{diag}(\bm{\sigma}_F)\bm{V}^T\), the minimization decouples into stretch space:
$$
\bigl(\bm{H}_{\psi}(\bm{\sigma}^{*})+\bar{k}\bm{I}\bigr)\bm{\sigma}^{*} = \bar{k}\,\bm{\sigma}_F.
$$
The stretch-space Hessian is
$$
\bigl[\bm{H}_{\psi}\bigr]_{ii} = \bar{\mu}\left(1+\sigma_i^{-4}\right) + \bar{\lambda}\,\frac{1-\ln\!\prod_j \sigma_j}{\sigma_i^2}, \qquad \bigl[\bm{H}_{\psi}\bigr]_{ij} = \frac{\bar{\lambda}}{\sigma_i\sigma_j}\quad (i\neq j).
$$
The converged projection is
$$
\bm{p}_e^{*}=\bm{U}\,\operatorname{diag}(\bm{\sigma}^{*})\,\bm{V}^{T}.
$$

The paper also describes an optional logarithmic volume-barrier variant,
$$
\phi^{\log}_e(\bm{\sigma}) = -\mu_e\sum_i \ln \sigma_i + \frac{\lambda_e}{2}\left(\sum_j \ln \sigma_j\right)^2,
$$
used when inversion avoidance is paramount [2605.14526].

## 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 [2605.14526]. To realize this, the method computes a single sparse factorization using nested dissection:
$$
\bm{P}\bm{A}\bm{P}^{T}=\bm{L}\bm{D}\bm{L}^{T},\qquad \bm{S}=\mathrm{diag}(d_i^{-1/2})\bm{L}^{-1}\bm{P},\qquad \bm{A}^{-1}=\bm{S}^{T}\bm{S}.
$$
Thus,
$$
\bm{A}^{-1}\bm{v}=\bm{S}^{T}(\bm{S}\bm{v}).
$$
The factor pair \((\bm{S},\bm{S}^T)\) is stored in persistent GPU buffers and reused everywhere unless topology, material parameters, or damping change [2605.14526].

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
$$
\Delta\bm{G} = \bigl[\bm{g}^{1}-\bm{g}^{0},\ldots,\bm{g}^{k-1}-\bm{g}^{k-2}\bigr],\qquad \Delta\bm{Q} = \bigl[\bm{q}^{1}-\bm{q}^{0},\ldots,\bm{q}^{k-1}-\bm{q}^{k-2}\bigr],
$$
with \(\bm{g}^k=\hat{\bm{q}}-\bm{q}^k\). The mixing coefficients are
$$
\bm{\gamma}^{*} = \bigl(\Delta\bm{G}^{T}\Delta\bm{G}+\rho_{\mathrm{aa}}\bm{I}\bigr)^{-1}\Delta\bm{G}^{T}\bm{g}^{k},
$$
and the extrapolated iterate is
$$
\bm{q}^{\mathrm{AA}} = (\bm{q}^{k-1}+\bm{g}^{k})-(\Delta\bm{Q}+\Delta\bm{G})\bm{\gamma}^{*}.
$$

Two safeguards are specified. If \(\|\bm{\gamma}^{*}\|_2>10\), the AA history is discarded and the solver falls back to the non-accelerated iterate; and the history window is reduced to \(m=1\) on heterogeneous meshes, versus \(m=5\) on homogeneous ones [2605.14526].

Convergence is checked by a dual-gate criterion,
$$
\frac{\|\bm{q}^{k+1}-\bm{q}^{k}\|} {\varepsilon_{\mathrm{rel}}\|\bm{q}^{k}\|+\varepsilon_{\mathrm{abs}}} \le 1
\quad\text{and}\quad
\frac{\|\bm{b}(\bm{q}^{k+1})-\bm{b}(\bm{q}^{k})\|} {\varepsilon_{\mathrm{rel}}\|\bm{b}(\bm{q}^{k})\|+\varepsilon_{\mathrm{abs}}} \le 1, \quad k\ge 1.
$$
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 [2605.14526].

## 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 [2605.14526]. In differentiable PD, the local proximal operator must be differentiated through the implicit function theorem. For Neo-Hookean materials, the relevant matrix is
$$
\bm{H}_{e}^{\mathrm{prox}} = \bm{H}_{\psi}(\bm{\sigma}^{*})+\bar{k}\bm{I},
$$
which may be indefinite.

The paper states that naively inverting this matrix can destroy descent directions in the adjoint computation and render gradients unreliable [2605.14526]. DiffPhD therefore introduces the filtered family
$$
\tilde{\bm{H}}_{e}^{\mathrm{prox}}(\tau) = (1-\tau)\bm{H}_{e}^{\mathrm{prox}} +\tau\lvert \bm{H}_{e}^{\mathrm{prox}}\rvert, \qquad \tau\in[0,1].
$$
The discrete choices are interpreted as \(\tau=0\) for the unfiltered inverse, \(\tau=\tfrac12\) for clamping-like behavior, and \(\tau=1\) for absolute-value filtering [2605.14526].

Selection is state-adaptive. Using the last converged PD increment
$$
\Delta\bm{q}^{*}=\bm{q}^{*}-\bm{q}^{*-1},
$$
the actual and model decreases are
$$
\Delta\Phi_{\mathrm{act}} = \Phi(\bm{q}^{*-1})-\Phi(\bm{q}^{*}), \qquad \Delta\Phi_{\mathrm{mod}} = \frac12\bigl|(\Delta\bm{q}^{*})^T\bm{A}\Delta\bm{q}^{*}\bigr|,
$$
giving the trust-region ratio
$$
\rho = \frac{\Delta\Phi_{\mathrm{act}}}{\Delta\Phi_{\mathrm{mod}}}, \qquad \tau^{*} = \begin{cases} \tfrac12, & |\rho-1|\le \varepsilon_{\mathrm{TR}},\\[4pt] 1, & \text{otherwise}, \end{cases}
$$
with \(\varepsilon_{\mathrm{TR}}=0.1\).

Once filtered, the local differential satisfies
$$
\tilde{\bm{H}}_{e}^{\mathrm{prox}}(\tau^*)\, d\bm{\sigma}^{*} = \bar{k}\,d\bm{\sigma}_F.
$$
This yields \(\partial\bm{p}_e^*/\partial\bm{F}_e\), which is inserted into
$$
\frac{\partial\bm{b}}{\partial\bm{q}}[\delta\bm{q}] = \sum_e w_e V_e \bm{G}_e^T \frac{\partial \bm{p}_e^*}{\partial \bm{F}_e} \bm{G}_e\,\delta\bm{q}.
$$

The paper presents this as the exact location at which stable hyperelastic gradients are obtained [2605.14526]. 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 [2605.14526]. With contact active, the framework solves a block KKT adjoint system whose \(\bm{q}\)-block uses the same sparse factor pair and whose \(\bm{\lambda}\)-block uses dense LDLT on a reduced contact system [2605.14526]. Gradients are propagated to states, forces, and material parameters, including projective weights,
$$
\frac{\partial\mathcal{L}}{\partial w_e} \mathrel{+} = V_e\bigl( \bm{\mu}^T\bm{G}_e^T\bm{p}_e^{*} - \bm{\mu}^T\bm{G}_e^T\bm{G}_e\bm{q}^{*} \bigr),
$$
and then to \(E_e\) and \(\nu_e\) by the chain rule [2605.14526].

## 5. Contact, damping, and unified reduced systems

DiffPhD incorporates contact through the Signorini–Coulomb complementarity formulation adopted from FBA [2605.14526]. For each contact, the normal complementarity residual uses a Fischer–Burmeister map,
$$
\varphi_n(\delta_{n,c},r_{n,c},\lambda_{n,c}) = \delta_{n,c}+r_{n,c}\lambda_{n,c} - \sqrt{\delta_{n,c}^{2}+r_{n,c}^{2}\lambda_{n,c}^{2}} =0,
$$
with
$$
r_{n,c}=h^2W_{n,cc},
$$
where the Delassus operator is
$$
\bm{W}=\bm{J}\bm{A}^{-1}\bm{J}^{T}.
$$

The reduced linear contact system is
$$
\bm{M}_{\mathrm{sys}}\Delta\bm{\lambda} = \bm{h}_{\mathrm{vec}} - \bar{\bm{J}\bm{A}^{-1}\bm{r}_{\mathrm{con}}}, \qquad \bm{M}_{\mathrm{sys}} = \bm{\Omega}\bm{W}\bm{\Omega}^{T}+\bm{E},
$$
with
$$
\bm{r}_{\mathrm{con}}=\bm{b}+\bm{J}^{T}(\bm{\Omega}\bm{\lambda}^k).
$$

A notable systems detail is that the Delassus matrix is assembled via batched sparse-dense multiplications,
$$
\bm{Y}=\bm{S}\bm{J}^T,\qquad \bm{W}=\bm{Y}^T\bm{Y},
$$
again reusing the persistent factor [2605.14526]. Position updates reuse cached columns \(\bm{a}_c=\bm{A}^{-1}\bm{j}_c\):
$$
\bm{q}^{k+1} = \bm{S}^{T}(\bm{S}\bm{b}) + \sum_{c=1}^{K}\omega_c\lambda_c^{k+1}\bm{a}_c.
$$

The same strategy governs damping. DiffPhD uses Rayleigh damping
$$
\bm{f}_d = -\bigl(\alpha\bm{M}+\bm{B}_{\bm{\beta}}\bigr)\bm{v}, \qquad \bm{B}_{\bm{\beta}} = \sum_e \beta_e V_e \bm{G}_e^T\bm{G}_e,
$$
which folds directly into the global operator
$$
\bm{A} = \frac{(1+\alpha h)\bm{M}}{h^2} + \sum_e \Bigl(w_e+\frac{\beta_e}{h}\Bigr)V_e\bm{G}_e^T\bm{G}_e.
$$
The elementwise damping coefficient is chosen proportional to stiffness,
$$
\beta_e=\beta_0\,\frac{\mu_e}{\mu_{\mathrm{ref}}}, \qquad \mu_{\mathrm{ref}}=\max_e \mu_e.
$$

The paper explicitly notes that this is opposite to the more standard strategy of damping soft regions more heavily [2605.14526]. 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 [2605.14526]. 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 [2605.14526]. The reported benchmarks include Cantilever, Armadillo twist, Crab, Gatorman, Napkin, Bouncing Ball, Plant, Bunny, Routing Tendon, Torus, Oreo, and Dice [2605.14526].

Before the summary table, two empirical claims organize the paper’s interpretation. First, DiffPhD is reported to remain convergent on stiffness contrasts up to \(100\times\), whereas prior PD solvers degrade in that regime [2605.14526]. 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 [2605.14526].

| Benchmark group | Examples | Reported outcomes |
|---|---|---|
| Heterogeneous forward simulation | Cantilever, Armadillo twist, Crab | Convergent up to \(100\times\) stiffness contrast; on Crab, up to about \(8.69\times\) forward speedup over DiffPD and about \(12.96\times\) over MAS |
| Contact-rich simulation | Gatorman, Napkin | On heterogeneous Napkin at \(50\times50\), up to about \(23.53\times\) backward speedup over DiffPD |
| Inverse problems and Real2Sim | Plant, Bunny, Routing Tendon, Torus, Dice | Plant loss \(0.029\) vs. \(0.692\) for DiffPD/MAS; Bunny heterogeneous loss \(0.132\) vs. \(0.403\) and \(0.374\); Routing Tendon \(0.441\) while DiffPD and MAS stall around \(201.674\); 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 \(0.029\) for DiffPhD versus \(0.692\) for DiffPD and MAS, **Bunny (heterogeneous)** at \(0.132\) versus \(0.403\) and \(0.374\), and **Routing Tendon** at \(0.441\) while DiffPD and MAS stall around \(201.674\) [2605.14526]. On **Torus**, the objective is matched or slightly improved, with a reported \(14.41\times\) speedup over MAS, though the paper notes that the benefit is smaller because contact itself dominates [2605.14526].

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 [2605.14526]. The **Oreo** example is presented as a forward-only capability check with extreme stiffness contrast across a multi-body gripper or object system [2605.14526].

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 [2605.14526]. 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 [2605.14526]. 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 [2605.14526].

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 [2605.14526].

Source: https://www.emergentmind.com/topics/diffphd