Introduction

Projected Gauss-Seidel, often shortened to PGS, is a simple iterative method for solving linear complementarity problems. It is widely used in contact mechanics and rigid-body simulation because it is easy to implement and often good enough for large sparse systems where exact solvers are too expensive.

What Problem PGS Solves

A linear complementarity problem asks for vectors x and w such that:

• 'w = Mx + q'
• every entry of x is non-negative
• every entry of w is non-negative
• for each index, at least one of x[i] or w[i] is zero

That last rule is the complementarity condition. In contact simulation, it matches the physical idea that a contact can either apply force or have a separating gap, but not both in the active sense at the same time.

The Core PGS Update

Ordinary Gauss-Seidel updates one variable at a time using the newest available values. Projected Gauss-Seidel adds one extra step: after computing the unconstrained update, it projects the value back onto the feasible region by clamping it to be non-negative.

For each component, the update is conceptually:

1. compute the residual using the current iterate
2. divide by the diagonal entry of M
3. clamp the result with max(0, value)

That projection is what makes the method suitable for an LCP.

A Runnable Python Implementation

The following NumPy example solves a small LCP and checks the complementarity conditions:

1import numpy as np
2
3
4def projected_gauss_seidel(M, q, max_iter=1000, tol=1e-8):
5    n = len(q)
6    x = np.zeros(n)
7
8    for _ in range(max_iter):
9        previous = x.copy()
10
11        for i in range(n):
12            sigma = M[i, :] @ x - M[i, i] * x[i]
13            x[i] = max(0.0, (-q[i] - sigma) / M[i, i])
14
15        if np.linalg.norm(x - previous, ord=np.inf) < tol:
16            break
17
18    w = M @ x + q
19    return x, w
20
21
22M = np.array([[2.0, -1.0], [-1.0, 2.0]])
23q = np.array([-1.0, -2.0])
24
25x, w = projected_gauss_seidel(M, q)
26print("x =", x)
27print("w =", w)
28print("x >= 0:", np.all(x >= -1e-8))
29print("w >= 0:", np.all(w >= -1e-8))
30print("complementarity:", np.allclose(x * w, 0.0, atol=1e-6))

This code is intentionally small, but it captures the standard PGS pattern used in larger solvers.

Why PGS Is Popular

PGS has three practical advantages.

First, it is cheap per iteration. You do not need to factor the whole matrix, which matters when the system changes every frame in a physics engine.

Second, it works well with sparse matrices. Many LCPs from contacts or constraints are sparse and structured, so each sweep is fast.

Third, it is easy to warm start. If yesterday's or the previous frame's solution is close to today's, you can initialize with that guess and reduce the iteration count.

These advantages explain why PGS appears so often in real-time simulation code even when more sophisticated methods have better theoretical guarantees.

Where It Struggles

PGS is still only an iterative approximation. Convergence can be slow when the matrix is poorly conditioned or when the coupling between variables is strong. It can also behave badly if diagonal entries are zero or too small, because each update divides by M[i, i].

In practice, convergence depends heavily on problem structure. Well-scaled systems with reasonable diagonal dominance tend to behave far better than arbitrary dense matrices.

Some solvers improve PGS with relaxation factors, ordering heuristics, or block updates. Those changes can help, but they do not remove the basic limitations.

Common Pitfalls

The first pitfall is applying PGS to an LCP that does not meet the assumptions your implementation relies on, especially positive diagonal entries.

The second is checking only the iterate difference and ignoring feasibility. A small change in x does not automatically mean the complementarity conditions are satisfied well enough.

The third is expecting high-accuracy solutions from very few iterations. In many real-time applications PGS is chosen because approximate answers are acceptable, not because it is the most accurate solver.

Summary

• Projected Gauss-Seidel solves LCPs by combining Gauss-Seidel updates with a non-negativity projection.
• It is popular in contact and constraint problems because each iteration is cheap and easy to implement.
• Sparse, well-scaled systems are where PGS tends to work best.
• Convergence can be slow or unreliable on poorly conditioned problems.
• Always verify feasibility and complementarity, not just whether the iterate stopped moving.