Source-linked AI summary
DiffPD: Differentiable Projective Dynamics
Tao Du, Kui Wu, Pingchuan Ma, Sebastien Wah, Andrew Spielberg, Daniela Rus, Wojciech Matusik
TL;DR
Differentiable soft-body simulation must balance stable implicit integration against expensive gradient computation. DiffPD combines Projective Dynamics with differentiable contact handling and reuses a prefactorized Cholesky system for backpropagation. It is 4-19 times faster than Newton’s method across several applications and reconstructs a real-world collision event from video.
Problem
Existing explicit differentiable simulators require tiny time steps, while implicit methods typically incur expensive linearized dynamics during backpropagation.
Method
DiffPD is an implicit FEM simulator based on Projective Dynamics that reuses a prefactorized Cholesky decomposition for backpropagation and supports penalty-based and complementarity-based contact models.
Results
4-19 times faster than Newton’s method across system identification, inverse design, motion planning, and closed-loop control applications, with a real-to-sim collision reconstruction demonstration.
Takeaways & Limitations
DiffPD provides differentiable soft-body simulation for learning, design, planning, control, and real-to-sim reconstruction across applications with nearly 30,000 DoFs.
Takeaways & Limitations
DiffPD assumes a restricted material-energy form and does not fully resolve differentiable complementarity-based contact and friction.
Abstract
from arXiv · showhide
We present a novel, fast differentiable simulator for soft-body learning and control applications. Existing differentiable soft-body simulators can be classified into two categories based on their time integration methods: Simulators using explicit time-stepping schemes require tiny time steps to avoid numerical instabilities in gradient computation, and simulators using implicit time integration typically compute gradients by employing the adjoint method and solving the expensive linearized dynamics. Inspired by Projective Dynamics (PD), we present Differentiable Projective Dynamics (DiffPD), an efficient differentiable soft-body simulator based on PD with implicit time integration. The key idea in DiffPD is to speed up backpropagation by exploiting the prefactorized Cholesky decomposition in forward PD simulation. In terms of contact handling, DiffPD supports two types of contacts: a penalty-based model describing contact and friction forces and a complementarity-based model enforcing non-penetration conditions and static friction. We evaluate the performance of DiffPD and observe it is 4-19 times faster compared with the standard Newton's method in various applications including system identification, inverse design problems, trajectory optimization, and closed-loop control. We also apply DiffPD in a reality-to-simulation (real-to-sim) example with contact and collisions and show its capability of reconstructing a digital twin of real-world scenes.
1 INTRODUCTION
DiffPD addresses the bottleneck in differentiable soft-body simulation by extending Projective Dynamics to efficient backpropagation, while supporting differentiable contact handling and broad applications.
- Differentiable soft-body simulators support optimization tasks such as system identification, motion planning, controller design, and inverse design.
- Tiny time steps limit explicit differentiable simulators, whereas implicit methods require expensive dynamics linearization during backpropagation.
- DiffPD extends Projective Dynamics to differentiable simulation by decoupling nonlinear dynamics and reusing a prefactorized global linear system during backpropagation.
- DiffPD is 4-19 times faster than Newton’s method across applications and achieves 9-11 times lower wall-clock time than deep reinforcement learning in closed-loop controller training.
- DiffPD supports penalty-based contact and friction models alongside complementarity-based non-penetration contact and static friction.
2 RELATED WORK
Related work spans differentiable soft-body simulation, Projective Dynamics, contact handling, and soft-robot design and control, with DiffPD combining PD-style efficiency and gradient computation.
- Differentiable simulation: Differentiable simulators have been developed for rigid bodies, soft bodies, cloth, and fluids, with explicit soft-body methods requiring small stable time steps.
- Projective Dynamics: Projective Dynamics offers an implicit alternative to Newton’s method and has been extended with material, rigid-body, conservation, and solver improvements.
- Projective Dynamics: DiffPD inherits the PD framework and augments it with gradient computation, while leaving GPU acceleration as future work.
- Contact handling: PD contact methods commonly use soft constraints or fictitious collision energies, which introduce task-dependent artificial stiffness requiring careful tuning.
- Soft robot design and control: Gradient-based soft-robot design and control methods use differentiability to improve efficiency, while some partially differentiable FEM methods trade accuracy for speed.
3 BACKGROUND
The background formulates implicit FEM simulation as nonlinear optimization, explains PD’s local-global solver, and identifies prefactorized linear systems as its efficiency source.
- Implicit time integration: Implicit time integration seeks positions x satisfying a nonlinear system derived from mass, external forces, and potential energy.
- Implicit time integration: Newton’s method repeatedly linearizes and factorizes the dynamics, creating a computational bottleneck for large-scale problems.
- Backpropagation with implicit time integration: Implicit backpropagation differentiates the equilibrium condition and solves an adjoint linear system, but reconstructing and refactorizing the Hessian remains expensive.
- Projective Dynamics: PD decouples nonlinear energies through local projections onto constraint manifolds and a global quadratic minimization over positions.
- Projective Dynamics: Each PD step is non-increasing in the surrogate objective, and convergence yields a critical point of the implicit integration objective.
- Projective Dynamics: The global PD matrix is constant and symmetric positive definite, so its Cholesky factorization can be precomputed and reused through back-substitution.
4 DIFFERENTIABLE PROJECTIVE DYNAMICS
DiffPD accelerates implicit backpropagation by splitting the Hessian into a constant global matrix and a parallelizable nonlinear component, while characterizing convergence through spectral radius.
- DiffPD reformulates the backpropagation Hessian into a constant global matrix and a local nonlinear component, paralleling PD’s forward-simulation decomposition.
- The local-global backpropagation solver computes nonlinear contributions locally and reuses the forward simulation’s Cholesky factorization for global back-substitution.
- Convergence is guaranteed from any initial guess when ρ(A^-1ΔA) < 1, although stronger theoretical results depend on the specific energy functions.
- Experiments report no convergence issues, suggesting the spectral-radius condition is likely satisfied in the tested applications despite the lack of a theoretical guarantee.
- The solver can be interpreted as a simplified quasi-Newton method with a constant Hessian approximation, while BFGS and line search offer further acceleration or safeguards.
5 CONTACT HANDLING
DiffPD supports penalty-based and complementarity-based contact handling within its differentiable PD framework. The complementarity approach enforces non-penetration and static friction while preserving efficient forward and backward computation through Cholesky-based linear algebra.
- Contact models: DiffPD offers penalty-based contact with static and dynamic friction, or complementarity-based contact enforcing non-penetration and static friction.The penalty model is simpler to integrate, while the complementarity model avoids scene-dependent parameters but supports static friction only.
- Penalty-based contact: Penalty-based contact adds a fictitious contact energy for penetration forces and incorporates friction forces directly into the external-force term.The friction treatment is needed because friction depends on nodal velocities rather than positions.
- Complementarity-based contact: Complementarity-based contact requires either zero penetration with no normal force or contact with a nonnegative normal force.The signed-distance function identifies whether a node is outside the obstacle or on the contact surface.
- Complementarity-based contact: The complementarity solver uses an active set of candidate contacts, solves a reduced system, then adds or removes contacts when force or penetration conditions are violated.Candidate nodes are fixed to target locations during the reduced solve, and the resulting forces and signed distances determine active-set updates.
- Efficient differentiation: Contact constraints modify the local-global PD steps while retaining convergence of the surrogate solver and compatibility with the prefactorized system.Because contact nodes are usually a small subset, row and column deletions are formulated as a low-rank update.
- Efficient differentiation: O(n^2c) to O(nc^2) reduces the overhead of applying the complementarity update beyond one linear solve.The reduced complexity applies to computing A^-1PU and the remaining operations associated with the update.
- Trade-offs: The complementarity contact algorithm trades exact Coulomb friction for differentiability, physical plausibility, and compatibility with DiffPD.The method enforces non-penetration and infinitely large static friction, while more physically accurate friction models exist.
6 EVALUATION
DiffPD is evaluated against explicit and implicit differentiable simulators across time-stepping, speed, contact handling, and optimization tests. The experiments show faster computation under supported conditions, while complementarity-based contact acceleration depends on a relatively small active contact set.
- Comparisons with Explicit Method: Implicit integration permits larger time steps and fewer stored frames than explicit integration, reducing backpropagation memory requirements.The evaluation describes a 10-millisecond implicit step and expects a 20× memory increase for the explicit method in the example.
- Comparisons with Explicit Method: Fig. 2 compares loss and gradient-norm changes for explicit and DiffPD simulations along 16 random directions near the initial nodal positions.Five directions are shown individually, while the aggregate plots report means and standard deviations across all 16 directions.
- Comparisons with Other Implicit Methods: At a 1e-4 convergence threshold, DiffPD agrees numerically with Newton’s method and is 12-16 times faster forward and 6.5-9 times faster in backpropagation.These measurements use 8 threads on the Cantilever example; the threshold is selected by comparison with Newton-Cholesky.
- Comparisons with Other Implicit Methods: For the Rolling sphere contact test, DiffPD is 10 times faster than Newton’s method in both forward and backward modes when the threshold reaches 1e-6.The low-rank update avoids refactoring, while caching A^-1_II×C adds 25% forward and 44% backward speedups with 8 threads.
- Ablation Study: Penalty-based contact produces qualitatively different duck sliding distances as friction decreases, and gradient optimization substantially improves landing at a target position.The tests use a 16,776-DoF duck on slopes and a curved surface with trainable friction coefficients.
- Ablation Study: Complementarity-based acceleration assumes a relatively small active contact set; this assumption may fail for co-dimensional objects such as ropes or cloth.The reported 3D contact applications had no more than 6% active contact nodes, while the Napkin study probes larger contact-set ratios.
7 APPLICATIONS
DiffPD is applied across system identification, optimization, control, and real-to-sim tasks, achieving comparable solutions with substantial computational savings. Contact-rich tasks reduce the speedup, while the method still improves optimization efficiency and reconstructs motions more closely to observed scenes.
- DiffPD achieves comparable application results while reducing time cost by almost an order of magnitude across five task categories.The categories are system identification, inverse design, trajectory optimization, closed-loop control, and real-to-sim applications.
- 7.1 System Identification: 9× speedup accelerates system identification while producing material parameters close to the reference plant.DiffPD estimates Young’s modulus at 1.00 MPa and Poisson’s ratio at 0.4, compared with the reference values of 1 MPa and 0.4.
- 7.2 Initial State Optimization: 9 times overall speedup improves initial-state optimization for the Bunny, although its optimized trajectory still misses the arbitrarily chosen target.The target is not guaranteed to be reachable, and DiffPD’s final loss is larger than Newton’s method’s.
- 7.3 Trajectory Optimization: 9× speedup solves the routing-tendon task almost perfectly, while contact-rich locomotion tasks achieve 4–6 times speedups.Locomotion optimization remains computationally difficult because it involves significant contact and can take hours with Newton’s method.
- 7.4 Closed-Loop Control: 8–19 times speedup accelerates closed-loop neural-controller optimization, which achieves comparable results to PPO with one or two orders of magnitude greater sampling efficiency.Wall-clock speedups against PPO are 9–11 times for both reported examples.
- 7.5 A Real-to-Sim Experiment: Optimizing image-space ball-position loss refines camera, initial-state, and contact-model estimates and produces motions that more closely mimic the video.The authors note that accounting for lens distortion or using a more accurate collision model could further improve the result.
8 LIMITATIONS AND FUTURE WORK
The paper identifies limitations in material-model support, contact realism, scalability, optimization, and contact-rich locomotion performance, and proposes corresponding future directions.
- Material models: DiffPD’s derivation does not straightforwardly support general hyperelastic materials because backpropagation requires reassembling the Hessian at each time step.The authors report no speedup over a direct solver for hyperelastic materials and leave acceleration as future work.
- Contact handling: DiffPD’s contact methods do not fully resolve differentiable complementarity-based contact and friction.The authors propose more realistic static and sliding friction models while retaining efficiency and differentiability.
- Scalability: The current largest examples contain thousands of elements and tens of thousands of DoFs, motivating scaling by one or two orders of magnitude.GPU implementations are suggested as a way to unlock larger applications.
- Locomotion and contact-rich applications: DiffPD is 4-6 times faster than standard Newton’s method in the reported locomotion examples, a less significant speedup than when PD assumptions hold.The authors suspect contact slows both forward simulation and backpropagation, especially under the sparse-contact assumption.
- Optimization: Gradient-based continuous optimization can terminate prematurely at local minima because the loss landscape is non-convex.The paper suggests global optimization or combinations with gradient-free methods such as genetic algorithms or reinforcement learning.
DiffPD: Differentiable Projective Dynamics (Supplemental Material)
The supplemental material identifies DiffPD as a work on Projective Dynamics and differentiable simulation, with the paper and supplement attributed to the listed authors.
- Bibliographic information: The paper’s listed authors are Tao Du, Kui Wu, Pingchuan Ma, Sebastien Wah, Andrew Spielberg, Daniela Rus, and Wojciech Matusik.The authors are affiliated with MIT CSAIL in the supplied material.
- Bibliographic information: The work is categorized under physical simulation in the ACM Computing Classification System.The supplied classification lists “Computing methodologies → Physical simulation.”
- Bibliographic information: The paper’s additional keywords are Projective Dynamics and differentiable simulation.These keywords identify the paper’s central technical area.
- Bibliographic information: The supplemental material is titled “DiffPD: Differentiable Projective Dynamics (Supplemental Material)” and is listed as a four-page October 2021 publication.The supplied citation includes the authors and DOI placeholder.
- Bibliographic information: The supplied preprint metadata identifies arXiv:2101.05917v3, dated 11 October 2021.This is the version and date shown in the supplemental material.
1 ENERGY DEFINITIONS
The supplemental material defines the quadratic elastic and muscle energies used in the examples and explains how their projections and gradients are handled during differentiation.
- Energy definitions: The supplemental material reviews quadratic energies used in the examples and connects each definition to its backpropagation gradient computation.This section provides energy-specific differentiation details rather than a single universal energy formula.
- Background elasticity: The corotated energy uses a deformation gradient F and the closest rotation matrix R, with gradients obtained through the polar decomposition.The discretized energy density uses the norm of F−R and is summed over Gaussian quadratures.
- Background elasticity: The supplemental material states that the background elastic model combines corotated and volume-preserving terms.These terms provide the elastic energy used in the examples described in the supplement.
- Volume-preserving energy: The volume-preserving energy projects F onto the manifold M = {D ∈ R3×3 : |D| = 1}.The projection is computed using an SVD of F and an optimization problem, after which D is reconstructed from U, ΣD, and V⊤.
- Volume-preserving energy: Gradients of the volume-preserving projection are derived by differentiating its KKT condition with respect to the singular values of F.Perturbing the singular values one at a time yields the corresponding changes in the optimal d∗.
- Muscle fibers: The muscle-fiber energy projects Fm onto a sphere whose radius is the actuation signal r.The muscle direction m is predefined per element, and gradients with respect to x and r follow by the chain rule.
2 EXPERIMENT DETAILS
The supplemental experiments apply optimization to material, initial-state, actuation, and controller variables across multiple soft-body examples, with figures tracking variables and loss over function evaluations.
- Plant: The Plant example optimizes Young’s modulus and Poisson’s ratio after deforming a fixed-root plant and releasing the external force.The ground-truth values are Young’s modulus 1e6 and Poisson’s ratio 0.4.
- Bouncing ball: The Bouncing ball example optimizes Young’s modulus and Poisson’s ratio from a specified initial height and velocity under gravity.The ground-truth values are Young’s modulus 2e6 and Poisson’s ratio 0.4; Fig. 2 tracks both variables and loss.
- Bunny: The Bunny example optimizes nine variables describing initial orientation, position, and velocity under bounds that prevent trivial target placement.Newton’s methods achieve lower final loss than DiffPD, but none of the methods solves the task perfectly.
- Routing tendon: All three methods produce satisfying Routing tendon solutions with very similar muscle-actuation magnitudes.Fig. 4 plots actuation magnitude and loss against function evaluations.
- Torus: The Torus task is sensitive to the initial guess, indicating a challenging loss landscape with more local minima than earlier examples.The torus uses muscle-driven mass redistribution and momentum to roll forward under the stated contact model.
- Quadruped: Quadruped locomotion is easier to optimize than the Torus task, possibly because the body is stable and sinusoidal controls reduce dimensionality.The example optimizes an open-loop sinusoidal control sequence.
- Cow: The Cow example uses a setup similar to the Quadruped example but with a different mesh model.Its optimization progress is shown in Fig. 7 through actuation magnitude and loss.
- Closed-loop control: The closed-loop controller maps state inputs to muscle actions through a neural network during forward simulation and backpropagation with DiffPD.For marine creatures, inputs include centroid and limb velocities plus relative limb positions; Fig. 8 depicts the forward and backward pipeline.