Add technical report + switch base IP solver to KLU

.gitignore

@@ -1,4 +1,5 @@
 *.mp4
+*.DS_Store
 
 # Files generated by invoking Julia with --code-coverage
 *.jl.cov

Project.toml

@@ -1,6 +1,6 @@
 name = "MixedComplementarityProblems"
 uuid = "6c9e26cb-9263-41b8-a6c6-f4ca104ccdcd"
-version = "0.2.3"
+version = "0.2.4"
 authors = ["David Fridovich-Keil <dfk@utexas.edu>"]
 
 [deps]

README.md

@@ -7,6 +7,16 @@ This package provides an easily-customizable interface for expressing mixed comp
 
 As of `v0.2.2`, `MixedComplementarityProblems.jl` implements CPU multithreading and GPU-parallelized solvers as well, enabled via `KernelAbstractions.jl`. Check out the benchmarking README [here](https://github.com/CLeARoboticsLab/MixedComplementarityProblems.jl/blob/main/benchmark/README.md) for more details.
 
+If you find this project useful in your work, please cite the accompanying [paper](https://arxiv.org/pdf/2608.00959):
+```
+@article{fridovich2026mcps,
+    title={MixedComplementarityProblems.jl: A Fast, Batched, Open-Source Interior Point Solver for Mixed Complementarity Problems},
+    author={David Fridovich-Keil},
+    year={2026},
+    journal={arXiv preprint arXiv:2608.00959}
+}
+```
+
 ## What are MCPs?
 
 Mixed complementarity problems (MCPs) are a class of mathematical program, and they arise in a wide variety of application problems. In particular, one way they can arise is via the KKT conditions of nonlinear programs and noncooperative games. This package provides a utility for constructing MCPs from (parameterized) games, cf. `src/game.jl` for further details. To see the connection between KKT conditions and MCPs, read the next section.
@@ -214,12 +224,17 @@ The returned `x`, `y`, `s` live on the GPU (as `CuArray`s); bring them back with
 > pattern the benchmarks use. `Adapt` is not a dependency of this package, so you would need
 > to add it to your own project (`] add Adapt`).
 
-> **Performance status (as of `v0.2.3`).** GPU now beats a many-threaded CPU run by a
-> consistent 2.5-2.8x on large-enough batched problems (e.g. the trajectory game once
-> `horizon ≳ 30`), though CPU remains faster or roughly at parity for smaller
-> per-instance problems — both handily beat PATH regardless. See the [benchmarking
+> **Performance status.** Both batched backends clear a batch far faster than sequential
+> `PATH` (the trajectory game runs ~60–90× faster on the multithreaded CPU, ~20–55× on the
+> GPU). Between the two backends the story is regime-dependent: the GPU wins end-to-end only
+> on large *dense* per-instance systems (e.g. randomly generated QPs with ≳128 primal
+> variables, 2–3×) or large batches of them, while the **CPU is faster on the trajectory
+> game** at all tested horizons. The GPU's batched sparse factorization does become cheaper
+> *per instance* than threaded KLU as the per-instance system grows, but the CPU's active-set
+> skip — it factorizes only the still-active instances each Newton step, whereas cuDSS always
+> processes the whole batch — keeps the CPU ahead end-to-end on the game. See the [benchmarking
 > README](https://github.com/CLeARoboticsLab/MixedComplementarityProblems.jl/blob/main/benchmark/README.md)
-> and PRs #54/#55 for the full breakdown and up-to-date numbers.
+> for the full breakdown and numbers.
 
 ## A fancier demo
 

benchmark/Project.toml

@@ -1,9 +1,14 @@
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MixedComplementarityProblems = "6c9e26cb-9263-41b8-a6c6-f4ca104ccdcd"

benchmark/README.md

@@ -115,11 +115,17 @@ julia> data = problem_size_scaling_benchmark();
 julia> problem_size_scaling_summary(data)
 ```
 
-**Current status (RTX 4090, as of 2026-07-04): GPU beats 32 CPU threads by a consistent
-2.5-2.8x on the trajectory game once the per-instance problem is large enough (`horizon
-≳ 30`, KKT dimension `d ≳ 2000`); below that, and for small dense QPs at small batch
-sizes, CPU is still faster or roughly at parity.** The headline number depends heavily on
-*which* regime you're in — see the breakdown below rather than quoting a single ratio.
+**Current status (RTX 4090, 32 threads, updated 2026-07-30): the CPU is faster than the GPU
+end-to-end on the trajectory game at every tested horizon; the GPU wins end-to-end only on
+large *dense* per-instance systems (≳128-primal QPs, 2–3×) or large batches of them.** Both
+batched backends clear a batch far faster than sequential `PATH` (game: CPU ~60–90×, GPU
+~20–55×). Raw data and analysis scripts are in [`benchmark/results/`](../results/).
+
+> **Correction.** An earlier version of this section (and PRs #54/#55) claimed the GPU *beats*
+> the CPU 2.5–2.8× on the game. That number was the GPU/CPU wall-clock *ratio* with the GPU in
+> the numerator — i.e. the GPU is ~2.5× **slower** — mislabeled as GPU-favorable. Re-measured
+> here (median-of-N), the raw ratios reproduce, but the direction is the opposite of the old
+> headline.
 
 - **cuDSS tuning:** `factorization_alg = "algo1"` (set manually in
   `ext/MixedComplementarityProblemsCUDSSExt.jl`, since `LinearAlgebra.lu()`'s convenience
@@ -128,31 +134,27 @@ sizes, CPU is still faster or roughly at parity.** The headline number depends h
   `algo2`-`algo5` are unsupported or slower; `reordering_alg`/`use_superpanels` don't help.
 - **Stall detection** (`max_stall_rounds` in `src/batched_solver.jl`) cut wall-clock
   2.7-7x on both devices by ending instances that neither converge nor diverge instead of
-  dragging every solve to `max_outer_iters` — necessary groundwork, but on its own it
-  *shrank* GPU's apparent advantage rather than growing it (GPU was previously winning
-  partly by handling wasted iterations better, not via a genuine linear-algebra edge).
-- **Why GPU wins at larger `horizon`:** raw per-instance `jacobian!+factorize!` cost
-  crosses over in GPU's favor right around `d ≈ 3500` (`horizon ≈ 50`) — confirmed by
-  isolated, apples-to-apples timing at fixed batch size, independent of solved fraction
-  or iteration count. `outer_iters` stays flat (9-16) across `horizon = 10..100`, so this
-  is **not** explained by `BatchedInteriorPoint`'s CPU-only active-set skip (which would,
-  if anything, favor CPU *more* as harder instances drop out early) — it's cuDSS's batched
-  factorization getting relatively cheaper per instance as the matrix grows, amortizing
-  its kernel-launch/occupancy overhead better than CPU's per-thread KLU factorization.
-- **Trajectory game, with the sampling/warm-start fix above** (`N = 1024`, `height = 50`):
-  GPU/CPU wall-clock ratio is 2.68x (`horizon=30`), 2.50x (`horizon=50`), 2.48x
-  (`horizon=70`), 2.81x (`horizon=100`) — consistently GPU-favorable. Solved fraction is
-  92% (`horizon=30`), 57% (`horizon=50`), 34-36% (`horizon=70`), 23-25% (`horizon=100`) —
-  a large improvement over the pre-fix collapse (as low as 0% at `horizon=100`), but still
-  short of a "realistic, high-confidence" benchmark at `horizon ≳ 50`; further tuning
-  (larger `lead_offset`/`lead_vy_boost`, or a genuinely better warm start) is still open.
-- At `horizon ≤ 20` (small `d`), CPU remains faster (GPU 2-2.5x slower) — GPU only pulls
-  ahead once there's enough per-instance work to amortize its overhead.
-- QP (small, dense, `num_primals = 32, num_inequalities = 16`): GPU/CPU ratio close to
-  parity (0.5-1.4x) across batch sizes; problem-size sweep (32/64/128 primals) is
-  non-monotonic (GPU wins at 128 primals, up to 2.9x, but loses at 64, 0.6-0.8x) —
-  confounded by the random QP generator's solved-fraction changing sharply with
-  `num_primals` (41% → 96% → 100%), not yet a clean isolated comparison.
-- `num_samples = 16384`+ currently OOMs for the trajectory game with an unexplained low
-  reported memory usage (~10%) at failure — not yet root-caused; dropped from the sweep
-  for now.
+  dragging every solve to `max_outer_iters`.
+- **The GPU factorization kernel *does* cross over — but end-to-end the CPU still wins.**
+  Isolated per-call timing (`percall_timing.csv`, all instances active) shows the GPU's
+  `jacobian!+factorize!` crossing over CPU's around `d ≈ 2500` (jac+fac GPU/CPU 1.92× at
+  `d=700` → 0.72× at `d=4900`) and `ldiv!` crossing even earlier (~`d≈1800`). But that
+  per-call number assumes the *whole* batch is factorized every step, which only holds on the
+  first Newton iteration. The real driver of end-to-end cost is `BatchedInteriorPoint`'s
+  **active-set skip**: on CPU each Newton step factorizes only the still-active instances
+  (cost ∝ active count), while cuDSS always factorizes the whole batch (flat cost — `active`
+  is a no-op on GPU). At `d = 3500` (`active_fraction_T50.csv`) the `jac+fac` GPU/CPU ratio
+  goes from 0.77× with all 1024 active to 8.1× with only 32 active. Since a real solve's active
+  set collapses fast as instances sub-converge, the CPU spends most of the solve in the regime
+  where it dominates. (This reverses the earlier claim here that the active-set skip was *not*
+  the explanation.)
+- **Trajectory game, end-to-end** (`N = 1024`, warm, median-of-5): `horizon=30` CPU 4.15s vs
+  GPU 8.01s (GPU 1.93× slower); `horizon=50` CPU 12.7s vs GPU 16.3s. Solved fraction (warm) is
+  93% (`horizon=30`), 57% (`horizon=50`) — matching earlier runs; cold-start collapses it
+  (≈17% at `horizon=30`), so cold high-horizon GPU/CPU ratios compare two mostly-failing
+  backends and aren't meaningful. GPU timings are also higher-variance at large horizon.
+- **QP** (`num_primals=32, num_inequalities=16`): the GPU pulls ahead as the batch grows
+  (GPU/CPU ~0.6× at `N=4096`, i.e. GPU ~1.7× faster), while the CPU is faster at small batch
+  sizes. Problem-size sweep (32/64/128 primals) is non-monotonic — GPU wins clearly at 128
+  primals (2–3×) but loses at 64 — confounded by the QP generator's solved fraction changing
+  sharply with `num_primals` (41% → 96% → 100%); not a clean isolated comparison.

benchmark/batched_benchmark.jl

@@ -33,6 +33,7 @@ function benchmark_throughput(
     run_batched = true,
     run_sequential_ip = true,
     run_path = true,
+    use_initial_guess = true,
 )
     @info "Generating random problems..."
     problem = generate_test_problem(benchmark_type; problem_kwargs...)
@@ -121,8 +122,10 @@ function benchmark_throughput(
     # Per-`benchmark_type` initial guess for the batched solve (e.g. the trajectory
     # game's zero-input rollout — see `generate_initial_guess` in
     # `trajectory_game_benchmark.jl`); `nothing` (cold zero start) for benchmark types
-    # that don't override it, matching prior behavior.
-    X₀ = run_batched ?
+    # that don't override it, matching prior behavior. Set `use_initial_guess = false` to
+    # force the batched solve to cold-start too, matching the (cold) sequential IP and
+    # PATH baselines for an apples-to-apples comparison.
+    X₀ = (run_batched && use_initial_guess) ?
         generate_initial_guess(benchmark_type, batched_mcp, Θ, device; problem_kwargs...) :
         nothing
 
@@ -144,7 +147,12 @@ function benchmark_throughput(
         tol,
         regularize_linear_solve,
     )
-    run_path && ParametricMCPs.solve(path_mcp, θs[1]; warn_on_convergence_failure = false)
+    run_path && ParametricMCPs.solve(
+        path_mcp,
+        θs[1];
+        convergence_tolerance = tol,
+        warn_on_convergence_failure = false,
+    )
 
     # --- Batched IP: one call over the whole batch, on `device`. ---
     batched = if run_batched
@@ -184,6 +192,7 @@ function benchmark_throughput(
             ParametricMCPs.solve(
                 path_mcp,
                 θ;
+                convergence_tolerance = tol,
                 warn_on_convergence_failure = false,
             ).status == PATHSolver.MCP_Solved
         end

benchmark/path.jl

@@ -81,8 +81,16 @@ function benchmark(
         ip_kwargs...,
     )
 
+    # Hold PATH to the same convergence tolerance as the IP solver for a fair comparison.
+    path_tol = get(ip_kwargs, :tol, 1e-6)
+
     @info "Warming up PATH solver..."
-    ParametricMCPs.solve(path_mcp, first(θs); warn_on_convergence_failure = false)
+    ParametricMCPs.solve(
+        path_mcp,
+        first(θs);
+        convergence_tolerance = path_tol,
+        warn_on_convergence_failure = false,
+    )
 
     # Solve and time.
     ip_data = @showprogress desc = "Solving IP MCPs..." map(θs) do θ
@@ -98,8 +106,12 @@ function benchmark(
 
     path_data = @showprogress desc = "Solving PATH MCPs..." map(θs) do θ
         # Solve and time.
-        elapsed_time = @elapsed sol =
-            ParametricMCPs.solve(path_mcp, θ; warn_on_convergence_failure = false)
+        elapsed_time = @elapsed sol = ParametricMCPs.solve(
+            path_mcp,
+            θ;
+            convergence_tolerance = path_tol,
+            warn_on_convergence_failure = false,
+        )
 
         (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
     end

benchmark/results/README.md

@@ -0,0 +1,75 @@
+# Experiments data (technical report §VII)
+
+Benchmark data behind the Experiments section of `docs/paper/main.tex`, collected 2026-07-30.
+Hardware/software: see `run_metadata.txt` (AMD Ryzen 9 7950X, 16 physical / 32 logical cores;
+NVIDIA RTX 4090 24 GB; Julia 1.12.6). All runs used `julia -t 32 --project=benchmark/gpu`.
+
+## Methodology
+- **Cold vs warm start.** We report both: *cold* (all solvers from their default start) is the
+  fair, apples-to-apples comparison against PATH and the unbatched solver, which are always cold;
+  *warm* (the batched solve seeded with `generate_initial_guess`, the trajectory game's zero-input
+  rollout) reflects realistic receding-horizon robotics use. Warm only affects the game (the QP has
+  no initial-guess override). Controlled by the `use_initial_guess` kwarg on
+  `benchmark_throughput` and the `WARMSTART` env var in `scripts/run_benchmarks.jl`.
+- **Distributions, not point estimates.** The QP's random draws are often infeasible by
+  construction, making per-instance solve times bimodal (fast feasible / slow-to-diverge
+  infeasible). Report mean±std, or violins when bimodal, from the raw per-sample CSVs — not a
+  single number. Batched sub-second timings are also noisy (GPU especially at large horizon), so
+  we sample many repetitions.
+- **Matched tolerance for PATH.** All solvers, including `PATH`, are held to the same
+  convergence tolerance `1e-4`. `PATH` runs at its built-in default `convergence_tolerance = 1e-6`
+  unless the option is set explicitly, which would hold it to a 100× tighter bar than our solver;
+  every PATH call now passes `convergence_tolerance = tol`. In practice this changes PATH's
+  timings and solved counts negligibly (quadratic local convergence drives the residual below
+  both thresholds in the same iteration) — verified by re-running: identical solved counts, times
+  within run-to-run noise. `scripts/rerun_path_tol.jl` re-timed PATH at `1e-4` and patched the
+  PATH rows of every CSV in place, leaving all batched/IP/GPU rows untouched (the pre-fix PATH
+  numbers, at `1e-6`, are in git history).
+
+## Key findings
+- **Batched solver vs PATH.** Both batched backends clear a batch far faster than sequential PATH
+  (game: CPU ~60–90×, GPU ~20–55×).
+- **CPU vs GPU is regime-dependent, and the CPU wins on the game.** The GPU only wins end-to-end on
+  large *dense* per-instance systems (128-primal QP: 2–3×) or large QP batches. This **corrects
+  PR #55's headline** ("GPU beats CPU 2.5–2.8× on the game"), which inverted the ratio direction /
+  over-generalized a per-call result. The GPU *factorization kernel* does cross over CPU at large d
+  (`percall_timing.csv`), but the CPU's active-set skip — it factorizes only the shrinking active
+  subset each Newton step, while cuDSS always processes the whole batch — keeps the CPU ahead
+  end-to-end (`active_fraction_T50.csv`: GPU/CPU 0.77×→8.1× as the active set shrinks 1024→32).
+- **KLU default.** Switching the unbatched interior-point solver's linear solve from UMFPACK to KLU
+  (in-place refactor reusing the symbolic analysis) is ~4–6× faster per solve at identical
+  reliability (`table1_summary.txt`), bringing the unbatched solver to ≈ parity with PATH on the QP
+  and ~50× faster on the game. Now the package default.
+
+## Files
+| file | what |
+|---|---|
+| `experiments.csv` | single-shot cold results: throughput, gpu_scaling, problem_size, horizon (cold+warm). `warm_start` column. |
+| `experiments_warm_reference.csv` | earlier full warm sweep (old schema, no `warm_start` col; QP rows cold-valid, game rows warm). |
+| `per_instance.csv` | raw per-instance solve time + status for PATH & unbatched IP (N=1024, QP+game) — for violins / bimodality. |
+| `per_instance_table1.csv` | Table I raw: per-instance PATH vs IP(UMFPACK) vs IP(KLU). `table1_summary.txt` has the summary. |
+| `per_rep.csv` | raw per-repetition full-batch wall-clock for batched CPU/GPU (+ PATH totals). Covers both the batch-size sweep (base configs `p32i16`/`T10` across `B`) and the problem-size sweep (`p64i32`, `p128i64`; game `T20`/`T30`/`T40`/`T50` at `B=1024`, cold+warm, median-of-5). |
+| `throughput_klu.csv` | clean single-run throughput (PATH vs sequential IP vs batched CPU) with the KLU default — the source for the report's Table II. |
+| `thread_scaling.csv` | batched throughput vs CPU thread count (1–32), QP and game — the source for the thread-scaling figure. |
+| `percall_timing.csv` | per-call `jacobian!`/`factorize!`/`ldiv!` CPU-vs-GPU vs horizon (d) — the kernel crossover. |
+| `active_fraction_T50.csv` | GPU-vs-CPU vs active-set size at fixed d=3500 — the mechanism behind CPU winning end-to-end. |
+| `confirm_game_T30_median.txt` | median-of-5 game T=30 CPU-vs-GPU (cold+warm). |
+| `run_metadata.txt` | hardware / git / Julia snapshot. |
+| `scripts/` | the driver + analysis scripts (repo root auto-resolved, or set `REPO_ROOT`). |
+
+## Reproduce
+From the repo root, one GPU stage at a time (single GPU — never overlap):
+```
+WARMSTART=0 STAGE=throughput   julia -t 32 --project=benchmark/gpu benchmark/results/scripts/run_benchmarks.jl
+#   STAGE ∈ {throughput, gpu_scaling, problem_size, horizon}
+julia -t 32 --project=benchmark/gpu benchmark/results/scripts/percall.jl        # kernel crossover
+julia -t 32 --project=benchmark/gpu benchmark/results/scripts/active_frac.jl     # active-set mechanism
+julia -t 32 --project=benchmark/gpu benchmark/results/scripts/table1_klu.jl      # Table I (KLU vs UMFPACK)
+STAGE=throughput_raw    julia -t 32 --project=benchmark/gpu benchmark/results/scripts/raw_data.jl   # STAGE ∈ {throughput_raw, scaling_raw, problem_size_raw}
+julia -t 32 --project=benchmark/gpu benchmark/results/scripts/rerun_path_tol.jl  # re-time PATH at 1e-4, patch PATH rows in every CSV
+```
+
+`problem_size_raw` (in `raw_data.jl`) fills the median-of-5 problem-size reps in `per_rep.csv`
+(QP `p64i32`/`p128i64`; game `T20`–`T50`); the base configs `p32i16`/`T10` already have ≥30
+reps from `throughput_raw`. `rerun_path_tol.jl` is CPU-only and must be run *after* any stage
+that regenerates a PATH-containing CSV, to restore matched-tolerance PATH numbers.

benchmark/results/active_fraction_T50.csv

@@ -0,0 +1,7 @@
+active_k,cpu_jacfac_ms,gpu_jacfac_ms,jacfac_GPU/CPU,cpu_ldiv_ms,gpu_ldiv_ms,ldiv_GPU/CPU
+1024,39.73,30.62,0.77,25.31,11.87,0.47
+512,20.31,29.25,1.44,16.45,11.77,0.72
+256,11.65,29.28,2.51,6.65,11.76,1.77
+128,7.44,27.91,3.75,4.56,11.77,2.58
+64,4.42,27.34,6.19,3.96,11.77,2.97
+32,3.36,27.14,8.08,3.84,11.77,3.07