The last stage of mj_forward turns the unconstrained acceleration into the constrained one. MuJoCo does this by solving a convex optimization problem every step, so the result is unique and the solver can be any method that finds the minimum. Three are built in: PGS, CG and Newton. They solve the same problem, from two directions and at very different costs. This lesson states the problem, reads the solver’s own statistics, and measures four things. It measures how fast each solver converges from one state, what the warmstart saves, what an under-converged solve does to a box on a slope, and what each solver costs as the system grows.
[!established] The problem the solvers solve Let $a_0$ be
qacc_smooth, $J$ the constraint Jacobian (efc_J), $a_\text{ref}$ the reference acceleration of the soft constraints (efc_aref, Lesson 9.2), and $R$ their regularizer (efc_R). The primal problem, over accelerations, is \(\min_{a}\; \tfrac12 (a - a_0)^\top M (a - a_0) + s(Ja - a_\text{ref}),\) with $s$ a convex penalty set by the constraint types and their states (Todorov, ICRA 2014; MuJoCo documentation, Computation). The dual problem, over constraint forces $f$ in the friction cones, is \(\min_{f}\; \tfrac12 f^\top (A + R) f + f^\top (J a_0 - a_\text{ref}),\qquad A = J M^{-1} J^\top.\) In the 3.14.0 source, PGS iterates onefc_force(the dual) and maps the result to joint space (engine_solver.c, the PGS section from line 454). CG and Newton are the “primal solvers” (from line 979) and iterate onqacc. They differ in one step: Newton multiplies the gradient by the inverse Hessian of the cost, CG by $M^{-1}$ (line 1528).
There is no lab for this lesson. The Lesson 20.1 lab stops just before mj_fwdConstraint writes efc_force and qacc. The picture to keep is a bowl: Newton sees its curvature and steps to the bottom, while CG and PGS walk down its slopes. Each step makes PGS or CG cheaper than Newton, but they need many more steps when the bowl is narrow, as it is for stiff, coupled contacts.
A tower of ten 5 cm boxes standing on a floor (60 dofs, 40 contacts, 160 constraint rows, one island), taken after 1 s. Each solver starts cold, from $a_0$ with zero force, and stops after $k$ iterations. The table gives the share of the constraint correction still missing, $\lVert a_k - a^\ast\rVert / \lVert a_0 - a^\ast\rVert$, where $a^\ast$ is Newton’s solution run to convergence:
| $k$ | 1 | 2 | 5 | 10 | 20 | 50 | 100 |
|---|---|---|---|---|---|---|---|
| PGS | 7.4 | 6.0 | 4.3 | 1.9 | 0.40 | 0.011 | 2.1e-05 |
| CG | 0.95 | 0.89 | 0.63 | 4.6e-07 | 7.3e-15 | 7.1e-15 | 6.8e-15 |
| Newton | 0 | 6.5e-15 | 6.6e-15 | 6.7e-15 | 6.6e-15 | 6.1e-15 | 6.6e-15 |
Newton reaches the minimum in one iteration here: with the active set right, the cost is quadratic and one Newton step solves it exactly. CG converges after about ten iterations. PGS overshoots at first, by more than seven times the correction; a plausible cause is that a sweep through stacked contacts pushes each box against forces not yet updated (inference). After 100 iterations it still misses by 2e-5. The solver records each iteration in mjData.solver, one block of mjNSOLVER = 200 entries per island, with solver_niter giving the count. At the default tolerance (1e-8) Newton stops after one iteration with an improvement of 110 and a gradient of 1.7e-14.
Each step starts from the better of two guesses: the previous step’s acceleration (qacc_warmstart) or $a_0$, compared by cost (warmstart in engine_forward.c, from line 1072). Over 2 s of the tower, the mean number of iterations per step falls from 75.4 to 41.0 for PGS and from 37.2 to 22.3 for CG. For Newton it falls from 1.0 to 0.2, because most steps start at the minimum already.
A box resting on a 20° slope with friction 0.5 should stay put, since $\tan 20^\circ$ = 0.36 is below 0.5. In 5 s every converged solve lets it creep 8.99 mm, the soft-contact creep of Lesson 9.3. Cut the iterations and the result changes:
| Solver, iterations | PGS 2 | PGS 5 | PGS 20 | PGS 100 | CG 2 | CG 5 | Newton 1 |
|---|---|---|---|---|---|---|---|
| distance slid in 5 s | 465.350 mm | 155.934 mm | 8.991 mm | 8.991 mm | 8.995 mm | 8.990 mm | 8.989 mm |
With two PGS iterations per step, static friction is never reached and the box slides almost half a metre. An iteration limit is part of the model: a policy trained with it learns a world where things slip.
At the default settings, time per step and iterations per step:
| Scene | nv | nefc | PGS | CG | Newton |
|---|---|---|---|---|---|
| tower of 10 boxes | 60 | 160 | 0.332 ms (40.2) | 0.079 ms (21.1) | 0.028 ms (0.1) |
| heap of 50 boxes | 300 | 1176 | 2.077 ms (98.8) | 0.501 ms (68.0) | 0.365 ms (3.9) |
| heap of 200 boxes | 1200 | 4300 | 17.377 ms (110.0) | 4.471 ms (61.6) | 11.180 ms (5.3) |
Newton needs few iterations, but each one solves with the Hessian, whose cost plausibly grows faster with size than CG’s multiplications by $M^{-1}$ (inference from the timings). On the 200-box heap CG becomes the cheapest: 4.5 against 11.2 ms. PGS is the slowest here on every scene.
[!recommendation] Choosing solver, iterations and tolerance Keep Newton, the default, for robots and scenes up to a few hundred degrees of freedom; it converges in a handful of iterations and holds stacks and friction exactly. Try CG when the system is large (here from about a thousand degrees of freedom) and measure. PGS was the slowest on every scene measured here; use it only after measuring an advantage on your problem. Never cut iterations to save time without checking a quantity that depends on convergence, such as a static-friction test like the slope. Read
solver_niterduring a run: a solver that keeps hitting the iteration limit is not converging.
INPUT: a 10-box tower, a box on a slope, heaps of 50 and 200 boxes; MuJoCo 3.14.0 solver defaults
PROCESS: cold-started solves stopped after k iterations against a converged reference; iterations per step with and without warmstart; slip on a slope with limited iterations; time and iterations per step
OUTPUT: printed tables
```python file=examples/l20_2_constraint_solver.py “"”Lesson 20.2: the constraint solver, PGS, CG and Newton, measured (Level 20).
INPUT a tower of 10 boxes (5 cm cubes, 0.1 kg, stacked on a floor) and a heap of N boxes dropped into a bin; MuJoCo 3.14.0 defaults: Newton, 100 iterations, tolerance 1e-8, pyramidal cones, warmstart on PROCESS (1) one state of the tower after 1 s: each solver started cold (warmstart disabled) and stopped after k iterations; the error of qacc against Newton run to convergence, and mjData.solver’s statistics; (2) iterations per step with and without warmstart over 2 s of the tower; (3) a box resting on a 20 degree slope with friction 0.5, each solver limited to a few iterations: how far it slides in 5 s; (4) cost: milliseconds per step and iterations per step for the tower and for heaps of 50 and 200 boxes OUTPUT printed tables
Run: python examples/l20_2_constraint_solver.py (about 1 minute) “””
import math import os import time
import mujoco import numpy as np
FAST = os.environ.get(“MJC_FAST”) == “1” SOLVERS = {“PGS”: mujoco.mjtSolver.mjSOL_PGS, “CG”: mujoco.mjtSolver.mjSOL_CG, “Newton”: mujoco.mjtSolver.mjSOL_NEWTON} SIDE = 0.025 # box half-size
def tower(n: int = 10) -> mujoco.MjModel:
boxes = ““.join(f”””
<body pos="0 0 {SIDE + 2 * SIDE * i:.4f}">
</mujoco>”””)
def slope(angle: float = 20.0, mu: float = 0.5) -> mujoco.MjModel:
“"”A box resting on a plane tilted by angle degrees about y.”””
a = math.radians(angle)
return mujoco.MjModel.from_xml_string(f””“
</worldbody></mujoco>”””)
def heap(n: int) -> mujoco.MjModel:
rng = np.random.default_rng(0)
side = int(np.ceil(np.sqrt(n / 4)))
boxes = ““.join(f”””
<body pos="{(i % side - side / 2) * 0.06:.3f} {((i // side) % side - side / 2) * 0.06:.3f} {0.1 + 0.06 * (i // side ** 2):.3f}"
euler="{rng.uniform(0, 90):.0f} {rng.uniform(0, 90):.0f} 0">
</mujoco>”””)
def settle(model: mujoco.MjModel, seconds: float) -> mujoco.MjData: data = mujoco.MjData(model) while data.time < seconds: mujoco.mj_step(model, data) return data
def solve_once(model, state_data, solver: str, iterations: int, tolerance: float = 0.0): “"”qacc, qacc_smooth and solver statistics from one cold-started forward pass of the given state.””” model.opt.solver, model.opt.iterations, model.opt.tolerance = SOLVERS[solver], iterations, tolerance model.opt.disableflags |= mujoco.mjtDisableBit.mjDSBL_WARMSTART data = mujoco.MjData(model) data.qpos[:], data.qvel[:] = state_data.qpos, state_data.qvel mujoco.mj_forward(model, data) stats = [(data.solver[i].improvement, data.solver[i].gradient) for i in range(data.solver_niter[0])] model.opt.disableflags &= ~int(mujoco.mjtDisableBit.mjDSBL_WARMSTART) return data.qacc.copy(), data.solver_niter[0], stats, data.qacc_smooth.copy()
if name == “main”: model = tower() state = settle(model, 1.0) print(f”MuJoCo {mujoco.version}; tower of 10 boxes after 1 s: nv {model.nv}, ncon {state.ncon}, nefc {state.nefc}, “ f”islands {state.nisland}”) print(“(1) one state, cold start: the share of the constraint correction still missing after k iterations,”) print(“ |qacc - qacc| / |qacc_smooth - qacc|, with qacc* from Newton run to convergence”) reference, n_ref, _, smooth = solve_once(model, state, “Newton”, 200, 1e-14) scale = np.linalg.norm(smooth - reference) ks = (1, 2, 5, 10, 20, 50, 100) print(f” reference: Newton, {n_ref} iteration(s); |qacc_smooth - qacc*| = {scale:.2f}”) print(“ “ + f”{‘k’:>8}” + ““.join(f”{k:>10}” for k in ks)) for name in SOLVERS: errs = [np.linalg.norm(solve_once(model, state, name, k)[0] - reference) / scale for k in ks] print(f” {name:>8}” + ““.join(f”{e:>10.1e}” for e in errs)) _, n_newton, stats, _ = solve_once(model, state, “Newton”, 100, 1e-8) print(f” Newton at the default tolerance 1e-8 stops after {n_newton} iterations; mjData.solver per iteration “ “(improvement, gradient): “ + “; “.join(f”{i:.1e}, {g:.1e}” for i, g in stats)) print(“(2) iterations per step over 2 s of the tower, default tolerance”) for name in SOLVERS: cells = [] for warm in (True, False): m = tower() m.opt.solver = SOLVERS[name] if not warm: m.opt.disableflags |= mujoco.mjtDisableBit.mjDSBL_WARMSTART d = mujoco.MjData(m) iters = [] while d.time < 2.0: mujoco.mj_step(m, d) iters.append(d.solver_niter[0]) cells.append(f”{‘warmstart’ if warm else ‘cold start’} mean {np.mean(iters[250:]):5.1f}, max {max(iters[250:]):3d}”) print(f” {name:>8}: “ + “; “.join(cells)) print(“(3) a box on a 20 degree slope with friction 0.5 (tan 20 deg = 0.36: it should stay): distance slid in 5 s”) for name, iterations in ((“PGS”, 2), (“PGS”, 5), (“PGS”, 20), (“PGS”, 100), (“CG”, 2), (“CG”, 5), (“CG”, 100), (“Newton”, 1), (“Newton”, 100)): m = slope() m.opt.solver, m.opt.iterations = SOLVERS[name], iterations d = mujoco.MjData(m) mujoco.mj_forward(m, d) start = d.xpos[1].copy() while d.time < (2.0 if FAST else 5.0): mujoco.mj_step(m, d) print(f” {name:>8}, {iterations:>3} iterations: {1000 * np.linalg.norm(d.xpos[1] - start):8.3f} mm”) print(“(4) cost per step, default settings (iterations used, time)”) scenes = [(“tower of 10”, tower(), 1.0)] + [(f”heap of {n}”, heap(n), 1.5) for n in ((50,) if FAST else (50, 200))] for label, m0, seconds in scenes: cells = [] for name in SOLVERS: m0.opt.solver = SOLVERS[name] d = settle(m0, seconds) t0, iters, steps = time.perf_counter(), [], 300 for _ in range(steps): mujoco.mj_step(m0, d) iters.append(sum(d.solver_niter[:max(d.nisland, 1)])) cells.append(f”{name} {1000 * (time.perf_counter() - t0) / steps:.3f} ms ({np.mean(iters):.1f} it)”) print(f” {label:<12} nv {m0.nv:>4}, nefc {d.nefc:>5}: “ + “, “.join(cells))
Output:
```text
MuJoCo 3.14.0; tower of 10 boxes after 1 s: nv 60, ncon 40, nefc 160, islands 1
(1) one state, cold start: the share of the constraint correction still missing after k iterations,
|qacc - qacc*| / |qacc_smooth - qacc*|, with qacc* from Newton run to convergence
reference: Newton, 1 iteration(s); |qacc_smooth - qacc*| = 31.02
k 1 2 5 10 20 50 100
PGS 7.4e+00 6.0e+00 4.3e+00 1.9e+00 4.0e-01 1.1e-02 2.1e-05
CG 9.5e-01 8.9e-01 6.3e-01 4.6e-07 7.3e-15 7.1e-15 6.8e-15
Newton 0.0e+00 6.5e-15 6.6e-15 6.7e-15 6.6e-15 6.1e-15 6.6e-15
Newton at the default tolerance 1e-8 stops after 1 iterations; mjData.solver per iteration (improvement, gradient): 1.1e+02, 1.7e-14
(2) iterations per step over 2 s of the tower, default tolerance
PGS: warmstart mean 41.0, max 80; cold start mean 75.4, max 81
CG: warmstart mean 22.3, max 39; cold start mean 37.2, max 49
Newton: warmstart mean 0.2, max 1; cold start mean 1.0, max 1
(3) a box on a 20 degree slope with friction 0.5 (tan 20 deg = 0.36: it should stay): distance slid in 5 s
PGS, 2 iterations: 465.350 mm
PGS, 5 iterations: 155.934 mm
PGS, 20 iterations: 8.991 mm
PGS, 100 iterations: 8.991 mm
CG, 2 iterations: 8.995 mm
CG, 5 iterations: 8.990 mm
CG, 100 iterations: 8.990 mm
Newton, 1 iterations: 8.989 mm
Newton, 100 iterations: 8.991 mm
(4) cost per step, default settings (iterations used, time)
tower of 10 nv 60, nefc 160: PGS 0.332 ms (40.2 it), CG 0.079 ms (21.1 it), Newton 0.028 ms (0.1 it)
heap of 50 nv 300, nefc 1176: PGS 2.077 ms (98.8 it), CG 0.501 ms (68.0 it), Newton 0.365 ms (3.9 it)
heap of 200 nv 1200, nefc 4300: PGS 17.377 ms (110.0 it), CG 4.471 ms (61.6 it), Newton 11.180 ms (5.3 it)
Objects slide or sink when they should rest. The solver is not converging. Check solver_niter against opt.iterations, then raise the iterations, switch to Newton, or reduce the stiffness that makes the problem hard.
A simulation becomes slow as objects pile up. The constraint count grew. Compare solvers on the real scene; for large systems CG can beat Newton.
Results change when the solver or tolerance changes. Either the solve is not converged, or the scene is sensitive, like a perfectly aligned tower whose symmetry any residual can break. Test convergence on a quantity that should not depend on the solver.
Disabling warmstart changed the trajectory. It changed where each iterative solve stopped. Converged solves should agree to the tolerance.
Repeat the slope test with cone="elliptic" and with impratio="10". How many PGS and CG iterations does each need to hold the box, and how far does the converged box creep?
Build heaps of 100 to 800 boxes and plot time per step against nv for each solver. Find the size at which CG overtakes Newton on your machine, and explain it from the cost of one iteration of each.
[!research] Solvers in learning pipelines GPU simulators favour solvers with fixed, parallel work per iteration (inference). Read how MJX (
mujoco.mjx, Lesson 20.3) and MuJoCo Warp choose and bound their solver iterations. Then repeat this lesson’s slope test there, because a cap chosen for speed changes the physics a policy learns from (inference from part 3). Todorov (ICRA 2014, in the reading list) gives the convex formulation these solvers share.
{"id": "20.2-check", "title": "Knowledge check", "questions": [
{"kind": "mcq", "q": "Which solver works on the dual problem, over constraint forces?",
"options": ["Newton", "CG", "PGS", "All three"],
"answer": 2,
"explain": "<p>PGS iterates on efc_force and maps it to joint space; CG and Newton are primal and iterate on qacc.</p>"},
{"kind": "mcq", "q": "Why does Newton converge in one iteration on the tower?",
"options": ["It is cheating with warmstart", "With the active set right, the cost is quadratic and one Newton step solves it exactly", "The tower has no contacts", "Tolerance is loose"],
"answer": 1,
"explain": "<p>The table's runs were cold-started; Newton's error after one iteration is 0 to machine precision.</p>"},
{"kind": "numeric", "q": "With PGS limited to 2 iterations, how far (mm, roughly) did the box slide on the slope in 5 s?",
"answer": 465, "tol": 15, "unit": "mm",
"explain": "<p>465.35 mm, against 8.99 mm of creep when converged.</p>"},
{"kind": "mcq", "q": "On the 200-box heap (nv 1200), which solver was fastest per step?",
"options": ["PGS", "CG", "Newton"],
"answer": 1,
"explain": "<p>CG 4.5 ms against Newton 11.2 ms and PGS 17.4 ms.</p>"}
]}
Lesson 20.3 measures where the rest of the time goes, and how to run many simulations at once.