qdldl: replay numeric refactorization from a precomputed schedule

[batterseapower]

Two commits, reviewable independently. Both are bit-identical to current main: they change how a numeric refactorization is executed, never what it computes.

Revised after a wider benchmark. An earlier version of this PR replayed every column block-wise, which cost 10–18% on small sparse LPs; the second commit now makes that choice per column. The numbers below are from the full corpus, not a fast subset.

Motivation

The numeric factorization dominates solve time. Profiling stock Clarabel with #[inline(never)] on the candidate kernels, so inlined frames cannot be misattributed:

	mm_exdata	kennington_cre_b
_factor_inner	86.4%	67.7%
refinement residual	7.5%	7.5%
triangular solves	2.9%	15.1%

Two properties of interior-point use make it improvable without touching the arithmetic.

Commit 1 — replay a precomputed schedule

The sparsity pattern of L is fixed after the first factorization, yet _factor_inner re-derives, for every column, the list of prior columns that update it: an elimination-tree walk with marker arrays, a work buffer and a reversal step. Clarabel refactors the same pattern once per iteration, so that control flow is pure overhead.

This records the discovered control flow once as a flat list of (source column, write position) steps and replays it. The update range of a step is Lp[cidx]..pos, because a column holds exactly its rows < k when that step runs, and pos is fixed once the pattern is. It 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 suits one-shot factorization, which interior-point use inverts.

Built lazily on the first refactor(), so single-factorization users pay nothing. The build self-verifies every recorded position against the existing pattern of L (Li[pos] == k, and the step count must equal nnz(L)); on mismatch, or if the pattern exceeds the u32 index storage, refactorization falls back to _factor_inner.

Measured alone on 67 public problems: −2.1% total, median −3.5%, 61 of 67 faster, one slower by 2.5%.

Commit 2 — replay updates block-wise on the columns where it pays

Row indices within a column are often long blocks of consecutive integers, so replaying a block as a contiguous slice update drops the per-entry index load and vectorizes:

y[rs..rs+len] -= Lx[p..p+len] * y_c        // instead of  y[Li[j]] -= Lx[j] * y_c

That is worth up to 55% on factors built from long blocks, but it costs bookkeeping per block — bounds, the clip against the end of the update range, the inner loop's setup — which a two- or three-entry block does not amortize. Applied unconditionally it costs 10–18% on small sparse LPs. (A per-block fast path does not help: for a short block the contiguous and indexed inner loops do identical work; the cost is the block loop itself.)

The trade is a property of each column, not of the matrix, so it is decided per column: a column's blocks are recorded only if the block containing a typical entry of that column reaches RUN_MIN_LEN entries — entry-weighted mean block length, Σlen²/Σlen — and the remaining columns are replayed entry-wise. A column with no recorded blocks takes the entry path, so the decision doubles as its own storage and needs no extra array. Factors of mixed structure, which are the common case, take the contiguous path on exactly the columns that benefit.

RUN_MIN_LEN is 16: eight iterations of a two-wide double-precision loop, comfortably past its prologue and epilogue, where two or four entries would be dominated by setup. It is the only constant introduced, and deciding per column is what removed the need for any whole-matrix threshold.

Correctness

Both commits perform the same floating-point operations in the same order as _factor_inner, and both replay paths in commit 2 do likewise, so L, D, Dinv, the inertia count and every dynamic-regularization decision are bit-for-bit unchanged. Asserted, not argued:

• test_refactor_matches_fresh_factor_exactly — a refactorization after updating every value equals a fresh factorization bitwise (L.nzval, D, Dinv, positive_inertia, regularize_count), checked on the first replay and a later one, and asserting the replay path was engaged rather than silently falling back.
• test_refactor_matches_fresh_factor_with_regularization — the same on a matrix built so dynamic regularization fires (asserted non-zero).
• test_both_replay_paths_agree_with_a_fresh_factor — the same for both paths of commit 2, on a banded matrix and a sparse quasidefinite one, chosen to sit on opposite sides of the per-column decision.
• test_run_selection_follows_block_structure — pins that the decision goes both ways, so neither path can rot untested.

Full suite green (cargo test, 20 binaries); cargo clippy clean on the touched file.

Benchmarks

Interleaved A/B — the two binaries run alternately problem-by-problem within each round, so machine drift affects both equally — minimum over rounds, single-threaded release build, one machine. Problems are every small-, medium- and large-band problem of a locally converted public corpus, i.e. selected by size rather than by where the change was expected to help.

67 small and medium problems (Netlib 19, Netlib-Kennington 14, Maros–Mészáros 13, SDPLIB 19, structured conic 11, Mittelmann 2), 2 rounds × 2 reps:

	total	median	faster	flat	slower
both commits	77.6s → 54.9s, −29.2%	−10.0%	54	12	1 (+2.3%, conic_pnorm15_winered_pow)

Largest gains mm_exdata −52.3%, sdplib_qap7 −46.7%, sdplib_theta1 −45.3%.

13 large problems (>10s), held out from the calibration of RUN_MIN_LEN, 1 rep:

problem	before	after	delta
sdplib_qap8	35.0s	17.1s	−51.2%
sdplib_qap9	135.2s	66.8s	−50.6%
sdplib_theta2	161.0s	79.7s	−50.5%
sdplib_control7	93.6s	48.4s	−48.3%
sdplib_control5	12.0s	6.3s	−47.8%
sdplib_control6	36.6s	19.2s	−47.6%
sdplib_truss8	199.4s	104.5s	−47.6%
sdplib_arch0	29.9s	17.7s	−40.9%
kennington_pds_20	253.6s	182.4s	−28.1%
mitt_rail4284	407.7s	299.1s	−26.6%
kennington_pds_10	19.4s	14.7s	−24.1%
mitt_rail2586	108.2s	97.9s	−9.5%
kennington_osa_60	19.6s	18.0s	−8.1%
total	1511.1s	971.8s	−35.7%

All 13 faster, none slower.

Regression gate: across all 80 problems, zero status changes, zero iteration-count changes and zero objective changes — required by bit-identity, and verified rather than assumed. The single timing regression is conic_pnorm15_winered_pow at +2.3%, close to this machine's noise floor.

Cost

nnz(L)×8B for the step schedule plus (n+1)×4B, and a block list for the qualifying columns only, allocated on first refactor. About 5MB on the largest problem measured here.

[batterseapower]

Correction to the benchmark evidence above, before you spend time on this.

The 14 problems I timed were selected as the slowest band of a larger corpus, and that selection flattered the result. Re-running on 67 public problems (every small- and medium-band problem across Netlib, Netlib-Kennington, Maros–Mészáros, Mittelmann, SDPLIB and structured conic — SDPLIB was missing before because my harness lacked the sdp feature), interleaved, 2 rounds × 2 reps:

	total	median	faster	flat	slower
this PR (both commits)	−27.7%	−4.6%	39	11	17
first commit only (schedule replay)	−2.1%	−3.5%	61	5	1
second commit only (run kernels)	—	—	17	11	26

So the aggregate is a large win, but it is not uniform: the second commit is bimodal. It is worth −51% on mm_exdata and −20% to −24% on mm_qshell, mm_stcqp1, mm_stcqp2, mitt_cont1 and kennington_cre_b, while costing +11% to +18% on 26 problems, mostly small Netlib LPs — netlib_nesm +18.5%, kennington_cre_c +16.1%, netlib_maros +16.0%, kennington_cre_a +14.3%.

In wall-clock terms the 17 problems that regress in this PR as a whole total 2.1s against a 21s saving, so nothing here is large in absolute terms. But they are reproducible at 4 rounds × 3 reps, they are bit-identical in output (so this is pure kernel overhead, not a different iteration path), and 12 of the 17 are Netlib LPs that anyone evaluating this would run.

I also want to retract the explanation I would have offered: the effect does not track run length. Problems regress at flop-weighted mean run lengths up to 12.2 and improve from 4.0, so there is no clean threshold on that statistic and I cannot currently gate the kernel on it honestly.

The first commit, by contrast, looks like what a bit-identical overhead removal should look like: 61 of 67 faster, one problem 2.5% slower, and it changes no floating-point result.

I would suggest treating the first commit as the reviewable change here and letting me withdraw the second until I understand its regression properly — happy to restructure the PR that way, or to close and resubmit, whichever you prefer. Apologies for putting an under-powered benchmark in front of you.

[batterseapower]

Revised. The description above is now the complete picture; this note records what changed since the correction I posted earlier, so the thread reads in order.

The regressions are gone, and they came from one specific decision. Isolating each commit on the 67-problem set showed the schedule replay was a near-universal win (61 of 67 faster, one 2.5% slower) while replaying every column block-wise was bimodal: worth up to 55% on long-block factors, costing 10–18% on small sparse LPs.

A whole-matrix threshold was the wrong fix. I first gated the block path on a matrix-level statistic — ≥80% of flop-weighted update work in blocks of ≥16 rows — which looked convincing on the 67 (−27.1%, zero regressions, and the boundary sat in a gap in the data). The held-out large problems then showed why it was wrong: they have mixed structure, scored below the threshold, and lost their speedups entirely — kennington_pds_20 −22.7% → −1.2%, kennington_pds_10 −7.5% → +1.4%. A single number per matrix cannot express "this matrix has both kinds of column", which is the common case.

Deciding per column fixes both ends and, on the mixed problems, beats replaying everything block-wise:

	all columns block-wise	matrix-level gate	per column
kennington_pds_20	−22.7%	−1.2%	−28.1%
kennington_pds_10	−7.5%	+1.4%	−24.1%
mitt_rail2586	−10.9%	−1.7%	−9.5%

because within one matrix the short-block columns now take the cheaper path too. It also removed the fitted 80% constant: what remains is RUN_MIN_LEN = 16, a vectorization threshold rather than a corpus-fitted one.

Final: −29.2% across the 67 small/medium problems with one +2.3% outlier, and −35.7% across the 13 large problems with all 13 faster. Bit-identical output throughout, so no status, iteration count or objective changes anywhere in the 80.