perf: tree-sparse LᵀDL factorization + O(n·depth) multibody solves

[haixuanTao]

The multibody mass matrix has branch-induced sparsity: M[i,j] ≠ 0 only when DOFs i and j lie on the same root-to-leaf path. Eliminating leaves-to-root factors M = Lᵀ·D·L with zero fill-in (Featherstone, RBDA ch. 6 — the same factorization MuJoCo uses in mj_factorM): L is nonzero only along ancestor chains.

What changes

• Factorization: O(n³) → O(n·depth²). Every solve (gravity RHS + all constraint-init unit-RHS back-solves in joint_constraints/contact_constraints/impulse finalize): O(n²) → O(n·depth).
• The per-DOF parent array is built from links_static inside gpu_mb_gravity_and_lu and stored in the former lu_pivots buffer — no binding or host-side changes anywhere.
• Factors stay in the dense mass_matrices storage (strict lower = L, diagonal = D); lu_solve_in_place keeps its signature, so all consumers compile unchanged.
• The factor + gravity solve run serially on lane 0 with zero workgroup barriers — data-dependent loops without barriers are legal under WGSL uniformity rules; the dense path needed ~5 barriers per column × MAX_MB_DOFS fixed iterations.
• Partial pivoting is dropped: the augmented M is SPD (armature/damping keep the diagonal positive), where unpivoted elimination is Cholesky-class stable — MuJoCo makes the same choice.

Validation

On a 12-DOF + free-base biped with ground contacts (WebGPU + native-CUDA backends, RTX 5090): max position drift vs the dense pivoted-LU path is 4e-7 after 10 steps, saturating ~7e-4 over 300 contact-rich steps — rounding-order drift, not divergence. Single-env frame times improve ~4% (WebGPU) / ~3% (CUDA-graph capture); the asymptotic win grows with DOF count and env count since the constraint back-solves dominate at batch scale.

Related perf PRs: #9, #13, #14.

🤖 Generated with Claude Code

[haixuanTao]

Batch-scale follow-up (same hardware, native-CUDA backend, 2048–8192 parallel envs, 12-DOF + free-base biped): this is a net regression of 10–15% steps/s vs the dense-LU path, and the DOF-scaling penalty it targeted is unchanged. The root cause is occupancy, not FLOPs: the sparse factor/solve runs serially on lane 0 while the dense path spread its (larger) work across all 32 lanes; at nv ≈ 18–31 that's more serial steps per warp, and the batched regime is throughput-bound. The single-env win reported above is real but comes from fewer barriers in a latency-bound regime that doesn't transfer.

Closing for now. The serial per-env elimination here is exactly the inner loop needed for an env-per-lane layout (one warp = 32 envs, MJX/Isaac-style mapping), which is the follow-up I'd pursue — that's an indexing/layout restructure on top of this math, so I'll reopen or supersede with that version.