Constraint Solvers: PGS, TGS, Newton, Convex

Collision detection handed the solver a list of contacts. Turning those into forces is, at heart, solving a linear system: find the contact impulses that stop bodies from interpenetrating. The choice of how to solve it is what separates a game engine from a robotics simulator.

The contact solve is a linear system

Stack the active constraints into a Jacobian $J$ that maps generalized velocity to relative velocity at the contacts. The impulses $\lambda$ that produce a target relative velocity $b$ satisfy

$$A\,\lambda = b, \qquad A = J\,M^{-1}J^{\top}.$$

$A$ is the Delassus operator: the inverse inertia seen in contact space. It is symmetric positive semidefinite, and its entry $A_{ij}$ measures how much an impulse at contact $j$ disturbs contact $i$. For a contact that can push but not pull, the impulse is one-sided, so the equality becomes a linear complementarity problem (LCP):

$$w = A\,\lambda - b, \qquad 0 \le \lambda \;\perp\; w \ge 0.$$

Either a contact carries force and its velocity error is zero, or it is separating and carries no force. A dense factorization of $A$ costs $O(n^{3})$ in the number of contacts, and $A$ changes every step. Engines never pay that; they iterate.

Splitting the matrix: Jacobi, Gauss-Seidel, PGS

Write $A = D + L + U$: diagonal, strictly lower, strictly upper. Each split gives a stationary iteration.

Jacobi updates every component from old values at once:

$$\lambda^{k+1} = D^{-1}\!\left(b - (L+U)\,\lambda^{k}\right).$$

Gauss-Seidel uses each freshly updated value immediately, which is the $(D+L)^{-1}$ split and converges faster:

$$\lambda_i \leftarrow \frac{1}{A_{ii}}\left(b_i - \sum_{j<i} A_{ij}\,\lambda_j - \sum_{j>i} A_{ij}\,\lambda_j\right).$$

Projected Gauss-Seidel (PGS) adds the $\lambda \ge 0$ projection after each component update, which is what keeps contacts from pulling:

$$\lambda_i \leftarrow \max\!\left(0,\; \lambda_i + \frac{1}{A_{ii}}\big(b_i - (A\,\lambda)_i\big)\right).$$

Applied one contact at a time, reusing the latest impulses, this is exactly the sequential impulse method that Box2D and Bullet run.¹²

A stack, solved by hand

Take a column of three unit boxes, $m=1$, $g=10$. Let $c_0$ be the ground contact, $c_1$ and $c_2$ the contacts above it, with a free surface $c_3 = 0$ over the top box. Static equilibrium of each box is $c_j = c_{j+1} + mg$: a contact holds up everything above it plus the box it touches. In matrix form,

$$\begin{bmatrix} 1 & -1 & 0 \\ 0 & 1 & -1 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} 10 \\ 10 \\ 10 \end{bmatrix}, \qquad \text{exact}\; \begin{bmatrix} c_0 \\ c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} 30 \\ 20 \\ 10 \end{bmatrix}.$$

Solving row $j$ for its diagonal unknown gives the PGS update

$$c_j \leftarrow \max\!\left(0,\; c_{j+1} + mg\right),$$

swept $j = 0, 1, 2$ over one shared array. The coupling is to $c_{j+1}$, the next index, which this sweep has not touched yet, so $c_{j+1}$ is the stale value from the previous iteration. Start from $\lambda = 0$:

Trace sweep 2: $c_0 \leftarrow \max(0, c_1 + 10) = 20$ using the stale $c_1 = 10$; then $c_1 \leftarrow \max(0, c_2 + 10) = 20$; then $c_2 \leftarrow \max(0, 0 + 10) = 10$. Three sweeps for three boxes. Because every row couples upward to a not-yet-updated neighbor, the true load propagates down exactly one contact per sweep, and the ground contact is the last to learn the full weight it carries.

PGS on a box stackinteractive

boxes N = 3sweeps K = 0

sweeps

0 of 3 needed

residual max|c − c*|

30.00

When one contact per sweep is too slow

That $N$-sweeps-per-$N$-stack behavior is the whole problem with PGS. A tall pile, a long chain, or a heavy body resting on a light one all couple many constraints, and information crawls through them one contact at a time. With a fixed iteration budget the solve stays soft and the stack sags, exactly as the demo shows below $K = N$.

Three families push past this. TGS, PhysX's temporal Gauss-Seidel, splits the timestep into substeps and re-solves each, which tightens convergence and handles large mass ratios far better than plain PGS at the same total cost.³ Newton solvers linearize the full nonlinear complementarity system and take a few large, well-scaled steps instead of many cheap ones. Convex solvers, such as MuJoCo's, recast the contact problem as a convex optimization and solve it with conjugate-gradient or Newton iterations, trading higher per-iteration cost for far fewer iterations and much better conditioning.² The engine's pick is a direct bet on how tightly its contacts couple: cheap sweeps where they do not, second-order steps where they do.

Footnotes

1. Catto, Iterative Dynamics with Temporal Coherence (GDC 2005). ↩

2. Bender, Erleben, Trinkle & Coumans, Interactive Simulation of Rigid Body Dynamics in Computer Graphics (CGF 2014). ↩ ↩²

3. NVIDIA PhysX 5.4.1 docs, Rigid Body Dynamics. ↩