DiffPhD

A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

Shih-Yu Lai1,2,7*, Sung-Han Tien1*, Jui-I Huang1*, Yen-Chen Tseng1, Yi-Ting Chiu1, Siyuan Luo3, Ziqiu Zeng3, Fan Shi3†, Peter Yichen Chen4, Tiantian Liu5, Yu-Lun Liu6, Bing-Yu Chen1†
1National Taiwan University 2MoonShine Animation Studio 3National University of Singapore 4The University of British Columbia 5Independent Researcher 6National Yang Ming Chiao Tung University 7Aalto University, Finland
* Equal contribution • † Corresponding author
Accepted to ACM Transactions on Graphics (TOG), Vol. 45, No. 6
(Proceedings of SIGGRAPH Asia 2026)
arXiv Paper Video Code

DiffPhD is a GPU-accelerated differentiable Projective Dynamics solver that handles heterogeneous materials, hyperelasticity under large deformations, and contact-rich interactions—within a single stable forward–backward pass.

Abstract

Differentiable simulation of soft bodies is a foundation for system identification, trajectory optimization, and Real2Sim transfer. Yet, existing methods such as the differentiable Projective Dynamics (DiffPD) struggle when faced with heterogeneous materials with extreme stiffness contrasts, hyperelasticity under large deformations, and contact-rich interactions, which are common scenarios in the real world. We present DiffPhD, a unified GPU-accelerated differentiable Projective Dynamics framework for heterogeneous materials that tackles these intertwined challenges simultaneously. Our key insight is a careful integration of: (i) stiffness-aware projective weights to embed heterogeneity into the global system; (ii) trust-region eigenvalue filtering lifted to the backward pass for stable hyperelastic gradients and a type-II Anderson Acceleration scheme with dual-gate convergence to stabilize forward iteration under large stiffness contrasts; and (iii) a unified GPU pipeline that reuses a single sparse factor across forward, backward, and contact computations, with stiffness-amplified Rayleigh damping folded into the same factor for heterogeneity-aware dissipation at zero recurring cost. DiffPhD delivers analytic adjoints—machine-precision on the elastic path, bounded and quantified under contact—with up to an order-of-magnitude speedup over prior differentiable solvers. Crucially, this speedup does not come at the cost of stability: DiffPhD remains convergent on stiffness contrasts up to 100× where prior PD solvers degrade. This unlocks end-to-end gradient-based optimization on regimes previously bottlenecked by either solver fragility or per-iteration cost—shell–joint composite creatures, soft characters wielding stiff weapons, and soft-gripper robotic manipulation—all handled within a single forward–backward pass.

Video

Method Overview

DiffPhD's design rests on one observation: Projective Dynamics' global stiffness matrix A is the shared backbone of every expensive operation—forward global step, Delassus contact compliance, backward adjoint—so the solver is organized around making A cheap to invert and robust under contrast. Heterogeneity is routed structurally through stiffness-aware projective weights baked directly into the global matrix; a trust-region filter is lifted onto the proximal-map Hessian in the backward pass to yield accurate Neo-Hookean gradients at high Poisson's ratios and large deformations; and a single persistent sparse-inverse factor is reused across forward, backward, and frictional-contact computations for an order-of-magnitude wall-clock improvement.

Figure 2: Two consequences of material heterogeneity in Neo-Hookean PD
Figure 2. Two consequences of material heterogeneity in Neo-Hookean PD. (a) The heterogeneous Neo-Hookean energy: a bar with stiff (red) and soft (blue) regions deforms qualitatively differently under stretch and twist, and the energy's logarithmic barrier prevents element inversion. (b) Consequence C1 (global). The spectrum of the global tangent stiffness widens directly with the stiffness contrast μmax/μmin—absorbing this contrast into the assembled operator A keeps its inversion robust. (c) Consequence C2 (local). Non-convexity is element-localized: a soft element near collapse develops an indefinite local Hessian, while its stiff neighbour stays positive definite. A uniform per-element filter under-serves both regimes, motivating a state-adaptive trust-region filter applied only in the backward pass.
Figure 4: Stiffness-aware projective assembly
Figure 4. Stiffness-aware projective assembly. (a) Heterogeneity is routed through the global operator A, not the local prox: per-element material parameters enter A only through projective weights we ∝ μe; the local prox sees only mesh-wide scalar means. (b) Why per-element local prox fails: mixing soft and stiff elements in the local Newton makes its projective Hessian element-varying, so per-element iterates ping-pong across PD iterations and Anderson Acceleration amplifies the drift into mesh-level divergence. (c) A-only routing preserves contraction: stiff weights enter A and widen its spectrum with contrast, but every element still runs the same Newton-on-stretches, so the PD fixed-point stays contractive.
Figure 5: Trust-region filter on the prox-map Hessian, applied in the backward pass only
Figure 5. Trust-region filter on the prox-map Hessian, applied in the backward pass only. (a) The IFT operator can be indefinite. Under high deformation or near inversion, the per-element prox-map Hessian picks up negative eigenvalues, making vanilla inversion through the backward IFT unstable. (b) Trust-region blend. A scalar τ ∈ [0, 1] interpolates between three canonical filters: unprojected (τ = 0), eigenvalue clamping (τ = 1/2), and absolute-value filtering (τ = 1). A state-adaptive rule selects τ* from the trust-region ratio evaluated on the last PD increment. (c) Backward-only, A reused. The forward fixed-point is left untouched—any forward modification would perturb the AA mixing and destabilise iteration. The filtered Hessian feeds the adjoint chain, and the trust-region ratio's quadratic model reuses the persistent factor of A at zero extra cost.
Figure 6: One persistent factor, shared across forward, backward, and contact
Figure 6. One persistent factor, shared across forward, backward, and contact. (a) One persistent factor. The SPD operator A is factorised once via METIS nested dissection into A−1 = STS; the factor pair sits in persistent GPU buffers. (b) Two-SpMV inverse application. Any A−1v is evaluated as two sparse multiplications ST(Sv); the dense inverse is never materialised. (c) Shared across all stages. The same factor pair serves the forward global solve, the contact Delassus operator, and the backward adjoint—refactorisation is triggered only on topology, material, or damping change.

Results & Demos

We validate DiffPhD on heterogeneous, hyperelastic, contact-rich benchmarks spanning forward simulation stability, differentiable inverse problems (system identification, initial-state and trajectory optimization), and Real2Sim robotic manipulation.

Heterogeneous Forward Simulation
Homogeneous · PD Energy
Heterogeneous · PD Energy
Homogeneous · Neo-Hookean
Heterogeneous · Neo-Hookean
Homogeneous vs. Heterogeneous
Armadillo (Twist)

A twisted Armadillo mesh height-partitioned into three sub-regions with stiffness ratios from 10× to 100×, compared under both the corotated PD energy and full Neo-Hookean hyperelasticity. DiffPhD remains stable across the entire range where DiffPD fails to converge beyond ~10×.

Heterogeneous Forward Simulation
Homogeneous
Heterogeneous (front)
Heterogeneous (back)
Homogeneous vs. Heterogeneous
Crab Pull-Up

A shell–joint composite crab with a stiff carapace (10×) coupled to soft joints is pulled upward by its carapace, cleanly lifting the soft legs off the ground—unlike the single-stiffness baseline.

Contact-Rich Forward Simulation
Mesh-to-mesh contact
Gatorman

A high-resolution Gatorman character wields a stiff sword (60× the body's reference stiffness) to strike a target sphere, exercising DiffPhD's unified Signorini–Coulomb contact pipeline under transient impact.

Contact-Rich Forward Simulation
Homo · 6%
Homo · 24%
Homo · 38%
Homo · 50%
Hetero · 6%
Hetero · 24%
Hetero · 38%
Hetero · 50%
Homogeneous vs. Heterogeneous
Napkin: Codimensional Contact

A one-layer napkin drapes over spherical obstacles of growing size, with contact sets from 6% to 50% of the mesh (columns, left to right). The heterogeneous variant (bottom row: 0.1× soft half vs. 1× stiff half) produces asymmetric, more pronounced draping than the homogeneous variant (top row).

Differentiable Inverse Problems
Initial Guess
Optimized
Ground Truth
System identification
Plant

System identification on an articulated potted-plant mesh whose slender branches induce strong geometric non-uniformity, disambiguating Young's modulus and Poisson's ratio from a single 100-frame oscillatory trajectory.

Differentiable Inverse Problems
Homogeneous · Initial Guess
Homogeneous · Optimized
Heterogeneous · Initial Guess
Heterogeneous · Optimized
Initial-state optimization
Bunny: Inverse Design

On the heterogeneous Bunny benchmark, DiffPhD optimizes initial position, orientation, and velocity so the bunny's center of mass reaches a target after 100 frames—converging where baseline solvers plateau on a suboptimal residual.

Differentiable Inverse Problems
Homogeneous · Initial Guess
Homogeneous · Optimized
Heterogeneous · Initial Guess
Heterogeneous · Optimized
Trajectory optimization
Routing Tendon

A vertical beam actuated through 16 muscle groups is optimized so its tip reaches a target at frame 100, contrasting homogeneous and heterogeneous muscle-fibre routing layouts through the same volume.

Real2Sim Robot Manipulation
Camera: Room In
Camera: Room Out
Extreme-contrast capability check
Oreo: Robot Manipulator

A Google Robot arm closes its 1013 Pa rigid claws on a five-body cookie stack (108 Pa ball, 4.55×105 Pa soft body), the only solver among those tested to remain numerically stable at this contrast.

Real2Sim Robot Manipulation
Initial Guess
Ground Truth
Optimized
Real2Sim recovery
Dice

Driven by real captured data of a UR5 arm pressing a soft die, DiffPhD recovers the in-plane contact location so the simulated deformation matches the captured video, closing the real-to-sim gap end-to-end.

BibTeX

@article{lai2026diffphd,
  author  = {Lai, Shih-Yu and Tien, Sung-Han and Huang, Jui-I and Tseng, Yen-Chen and Chiu, Yi-Ting and Luo, Siyuan and Zeng, Ziqiu and Shi, Fan and Chen, Peter Yichen and Liu, Tiantian and Liu, Yu-Lun and Chen, Bing-Yu},
  title   = {DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration},
  journal = {arXiv preprint arXiv:2605.14526},
  year    = {2026}
}