qdldl: replay a precomputed schedule on numeric refactorization

[batterseapower]

What

After the first factorization the sparsity pattern of L is fixed, but _factor_inner re-derives the per-column update lists on every refactor: an elimination-tree walk with marker arrays, a work buffer, and a reversal step. Interior-point use refactors the same pattern once per IP iteration, so that control flow is pure overhead.

This PR records the control flow once — a flat list of (cidx, pos) update steps per column, exactly as _factor_inner executes them — and replays it on subsequent refactorizations. The update range of a step is Lp[cidx]..pos because a column holds exactly its rows < k when step (k, cidx) runs, and pos (the write position of L[k,cidx]) is a fixed index once the pattern is fixed. This is the standard symbolic/numeric phase separation of sparse direct solvers (T. Davis, Direct Methods for Sparse Linear Systems, SIAM 2006, ch. 4); QDLDL's single-phase design favours one-shot factorization, which Clarabel's refactor-heavy usage inverts.

Why it is correct

The replay performs the identical floating-point operations in the identical order as _factor_inner(..., logical_factor=false) — same A-column scatter, same update loop bounds, same pivot regularization tests on the same values. So L, D, Dinv, the inertia count and the dynamic-regularization decisions are bit-for-bit identical.

Defensive design: the schedule is built lazily on the first refactor() (one-shot factorizations pay nothing), the build self-verifies every recorded position against the existing pattern (Li[pos] == k per step, total = nnz(L)), and any mismatch — or a pattern too large for the u32 step storage — falls back to the original path.

No tunable constants. The only magic number is the u32 index width, which is a storage choice guarded by an explicit overflow check and fallback.

Validation suites used throughout this series

All timing is on one Apple M-series core, single-threaded, --release with
debug symbols, and every comparison is interleaved: the two binaries are
run alternately problem-by-problem within each round, so thermal drift and
machine noise affect both sides equally. Reported times are the minimum over
rounds. (Session-to-session noise on this machine is large — up to 20% between
identical binaries at different times — so non-interleaved comparisons are not
trustworthy and none are quoted.)

1. Unit / integration tests. cargo test — 20 test binaries, all passing,
plus the specific new tests listed per PR below. cargo clippy clean on the
touched files.

2. Portfolio-rebalance set (9 problems, the motivating workload).
Real conic problems from a production portfolio-rebalance backtest:
n = 9,448–18,424, m = 22,329–43,273, nnz(A) = 271k–538k, diagonal P, two
large nonnegative cones and one second-order cone of dimension 3,279–6,271
(a factor-model risk constraint: ~156 dense-ish factor rows plus a diagonal
idiosyncratic block). Used for the headline timings.

3. In-sample public corpus (224 problems). Every tiny/small/medium
non-PSD problem in a locally converted corpus: 101 Maros–Mészáros QPs,
85 Netlib LPs, 24 structured conic problems (SOCP/EXP/POW built from UCI
data), 12 Netlib-Kennington LPs, 2 Mittelmann LPs. Each has an
independently verified reference objective. Used as the status/objective
regression gate for every change.

4. Out-of-sample suites (92 problems, fetched fresh from the internet).
Because suite 3 became a tuning gate, two further suites were added that
were never used to guide any decision:

• MIPLIB 2017 benchmark set as LP relaxations — 45 instances. Downloaded
  from miplib.zib.de; the MPS reader ignores INTORG/INTEND markers, so
  reading a MIPLIB file yields exactly the continuous relaxation. These are
  substantially harder for an interior-point method than the portfolio
  problems. 90s time limit per solve.

• CBLIB (Conic Benchmark Library) — 47 instances, sampled across families
  from cblib.zib.de, including the DIMACS classics nb, nb_L1, nb_L2,
  nql30/60/180, qssp30/60, sched_*. Converted with a new CBF reader;
  integer instances are taken as continuous relaxations, and rotated
  quadratic cones are mapped to second-order cones by the orthogonal
  rotation u=(x1+x2)/√2, v=(x1−x2)/√2.

  Conversion validated independently: nb solves to −5.0703094644e-2,
  matching its published DIMACS optimum (−0.05070309), and HiGHS (installed
  for the purpose) agrees on the LP-only conversions, e.g. gen_ip054
  6765.209042728 vs Clarabel 6765.2090428.

Tests specific to this PR

• test_refactor_matches_fresh_factor_exactly — refactor after value updates must equal a fresh factorization bitwise (L.nzval, D, Dinv, inertia, regularize count) on a 70×70 quasidefinite matrix with genuine fill-in, checked twice (first replay and a subsequent one), and asserts the replay path was actually engaged rather than silently falling back.
• test_refactor_matches_fresh_factor_with_regularization — the same, on a matrix constructed so dynamic regularization actually fires (asserted).

Performance: before → after

Interleaved, min of 3 rounds × 2 reps. Objectives are bitwise-equal and statuses and iteration counts identical on every problem below, which is the strongest possible quality result: the solver follows exactly the same trajectory, only faster.

Portfolio-rebalance set (9 problems)

problem	before	after	delta
20260101_full	1.625s	1.560s	−4.0%
20260102_full	1.716s	1.685s	−1.8%
20260102_screened	0.988s	0.967s	−2.2%
20260105_full	2.851s	2.817s	−1.2%
20260105_screened	0.785s	0.767s	−2.2%
20260115_full	1.629s	1.600s	−1.8%
20260115_screened	0.984s	0.950s	−3.5%
20260116_full	1.486s	1.465s	−1.5%
20260116_screened	0.848s	0.826s	−2.6%
total	12.911s	12.635s	−2.1%

Public corpus (7 timing-relevant problems)

problem	before	after	delta
kennington_osa_30	4.105s	3.883s	−5.4%
kennington_ken_13	0.968s	0.922s	−4.8%
mm_q25fv47	0.275s	0.264s	−3.9%
conic_logreg_winewhite_exp	0.330s	0.318s	−3.8%
kennington_pds_06	3.560s	3.428s	−3.7%
kennington_cre_b	5.217s	5.032s	−3.5%
mitt_rail582	2.200s	2.160s	−1.8%
total	16.654s	16.006s	−3.9%

All 16 problems individually faster. The win is modest because these factors are flop-bound (measured Σᵢlᵢ² ≈ 42M on the largest problem, dominated by a dense trapezoid from ~156 dense-ish rows); the removed overhead grows in relative terms the sparser the factor, which is why the sparse LPs gain more.

Status/objective regression gate

Bit-identity makes this vacuous but it was run anyway: 224-problem in-sample corpus — 0 status changes, 0 iteration-count changes, 0 objective changes.

Cost

nnz(L)×8B + (n+1)×4B for the schedule (≈4.7MB on the largest problem here), allocated on first refactor only.

Note for reviewers

#9 is stacked on this branch and adds the run-decomposition kernels; the combination is −12.7% on the portfolio set. This PR's test helper test_matrix_quasidef and #10's test_matrix_kkt_like are near-duplicate random-KKT generators on independent branches — worth deduping into one helper once both land.

🤖 Generated with Claude Code

https://claude.ai/code/session_01CqGdm2vZZ6HA5HT8DsRePa

[batterseapower]

Public-corpus timing, measured per track

Added because the earlier corpus timings in this description used a narrow 7-problem set. Re-measured on 14 timing-relevant public problems (Netlib-Kennington cre_b/cre_d/pds_06/osa_07/osa_14/osa_30/ken_11/ken_13, Mittelmann cont1/rail582, Maros–Mészáros exdata, Netlib pilot/d2q06c, structured conic pnorm15_winewhite_pow), interleaved, 3 rounds × 2 reps:

track	before	after	delta
#6 + #9 (factor: schedule + run kernels)	39.991s	31.277s	−21.8% — bit-identical throughout
#8 + #10 (refinement: backend loop + symmetric residual)	39.563s	37.663s	−4.8%
#11 (fused residual products)	39.517s	39.038s	−1.2%
all five together	39.530s	29.134s	−26.3%

All 14 problems individually faster with the full stack (−6.9% to −56.3%; mm_exdata 8.178s → 3.575s), 12 of 14 with bitwise-identical objectives.

Worth stating explicitly since the original motivation was a portfolio workload: this is not a workload-specific optimisation. The public corpus gains (−26.3%) are of the same order as the portfolio set (−42.6%), and the largest single component (#6 + #9, −21.8%) changes no floating-point result at all.

[batterseapower]

Submitted upstream as oxfordcontrol#230 (both commits of this track together, since a cross-repo PR must target upstream main and so cannot be stacked). Keeping this open as the fork-side record.