LCP vs PGS Solvers in Game Development

In real-time game physics, rigid body contact and constraint resolution are commonly formulated as a Linear Complementarity Problem (LCP). While Projected Gauss-Seidel (PGS) is technically an iterative algorithm used to solve an LCP, the game development industry typically distinguishes "LCP solvers" as direct, active-set pivoting methods (such as Lemke’s or Dantzig’s algorithms) and "PGS solvers" as iterative, matrix-splitting methods. This article examines the core structural differences between these two approaches, focusing on matrix assembly, memory management, execution predictability, and handling of frictional constraints.

Mathematical Formulation and Solution Strategy

Both approaches resolve the velocity-level constraint equations subject to complementarity conditions:

\[w = A\lambda + b \ge 0, \quad \lambda \ge 0, \quad \lambda^T w = 0\]

Where \(A = J M^{-1} J^T\) is the Delassus operator, \(J\) is the constraint Jacobian, \(M\) is the system mass matrix, \(b\) is the constraint violation/bias term, and \(\lambda\) represents the constraint impulses.

• Direct LCP Solvers (Pivoting Methods): Use combinatorial search methods such as Lemke’s algorithm or Dantzig’s principal pivoting method. Structurally, they partition constraints into active and inactive sets, iteratively swapping indices between sets and solving linear sub-systems until the complementarity conditions are satisfied exactly.
• PGS Solvers (Iterative Splitting): Decompose the system matrix \(A\) into lower, diagonal, and upper components (\(A = L + D + U\)). The solver updates each constraint impulse \(\lambda_i\) sequentially while immediately clamping (projecting) the result onto the admissible set \([0, \infty)\) or friction limits:

\[\lambda_i^{k+1} = \max\left(0, \lambda_i^k - \frac{1}{A_{ii}} \left( b_i + \sum_{j < i} A_{ij}\lambda_j^{k+1} + \sum_{j > i} A_{ij}\lambda_j^k \right)\right)\]

Matrix Assembly vs. Matrix-Free Architectures

A fundamental structural difference lies in how data is stored and accessed during execution.

Direct LCP solvers require explicit construction of the Delassus matrix \(A\). Assembling \(A = J M^{-1} J^T\) requires sparse matrix-matrix multiplications and dynamic memory allocation to manage pivoting steps and matrix factorizations (such as Cholesky or LU updates). This introduces significant cache thrashing and memory overhead when handling large constraint islands.

PGS solvers are almost universally implemented as matrix-free algorithms, often referred to in game engines as Sequential Impulses (SI). Because PGS updates constraints one row at a time, the engine does not need to build \(A\) globally. Instead, it computes the scalar diagonal \(J_i M^{-1} J_i^T\) once per constraint, reads body velocities directly from memory, calculates the impulse, and applies the delta directly back to body velocities. This keeps data localized in CPU L1/L2 caches and makes vectorization (SIMD) straightforward.

Time Complexity and Frame-Budget Predictability

Game engines operate under strict per-frame time budgets (typically 16.6ms or 8.3ms).

• Direct LCP Solvers: Have a worst-case computational complexity of \(O(2^n)\) or \(O(n^3)\) operations, where \(n\) is the number of constraint rows. If the pivot search hits degenerate conditions, computation time spikes drastically, leading to visible frame drops. Furthermore, direct methods cannot be terminated halfway; stopping early yields an incomplete pivot state with invalid velocities.
• PGS Solvers: Exhibit \(O(k \cdot n)\) complexity, where \(k\) is a fixed number of iterations (typically 4 to 10 in modern titles). PGS is an "anytime" algorithm: it can be terminated unconditionally when the frame budget expires. While early termination results in softer constraints or slight penetration, it avoids simulation crashes and maintains frame-rate stability.

Warm Starting and Temporal Coherence

Physics engines exploit temporal coherence across consecutive frames:

• In PGS Solvers: Warm starting is trivial. The impulse \(\lambda\) calculated in frame \(t\) is cached and applied as the initial guess \(\lambda^0\) in frame \(t+1\). This dramatically reduces the number of iterations required to reach visual stability for resting stacks of bodies.
• In Direct LCP Solvers: Caching the active set from the previous frame can accelerate initial pivoting, but changes in contact topology require rebuilding factorizations and updating pivot bases, reducing the performance advantage gained from temporal coherence.

Friction and Non-Linear Limits

Under Coulomb friction, the friction limits depend dynamically on the normal contact force: \(|\lambda_t| \le \mu \lambda_n\).

Direct LCP solvers struggle with this coupling because it turns the system into an NP-hard non-monotone or non-linear complementarity problem. Engines using direct LCPs typically linearize the friction cone into a polyhedral pyramid, adding multiple constraint rows per contact and further bloating matrix dimensions.

PGS naturally handles coupled friction limits through sequential projection. The normal impulse \(\lambda_n\) is computed and clamped first; immediately following, the tangential impulses \(\lambda_t\) are evaluated and clamped against the newly calculated \([-\mu \lambda_n, \mu \lambda_n]\) bound within the same iteration loop. While this Gauss-Seidel-style decoupling is an approximation that can cause minor shear inaccuracies, it is computationally robust and visually indistinguishable in real-time applications.