qdldl: replay factor updates over consecutive-row runs - #9
Conversation
Public-corpus timing, measured per trackAdded because the earlier corpus timings in this description used a narrow 7-problem set. Re-measured on 14 timing-relevant public problems (Netlib-Kennington
All 14 problems individually faster with the full stack (−6.9% to −56.3%; 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. |
|
Submitted upstream as oxfordcontrol#230 (both commits of this track together, since a cross-repo PR must target upstream |
The dominant loop of the numeric factorization scatters into the sparse
accumulator, one indexed load and one indexed store per entry:
for j in range { y[Li[j]] -= Lx[j] * y_c }
Row indices within a column of L are often long blocks of consecutive
integers -- a set of dense-ish rows, which AMD orders last, produces a
trailing trapezoid shared by many columns. Replaying such a block as a
contiguous slice update drops the per-entry index load and lets the
compiler vectorize:
y[rs..rs+len] -= Lx[p..p+len] * y_c
That is worth up to 55% of solve time on factors built from long blocks.
But it costs bookkeeping per *block* -- the bounds, the clip against the
end of the update range, and the inner loop's own setup -- and on factors
whose blocks are two or three entries long that is not amortized: measured
across public problems, replaying every column block-wise 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 the same work, and the cost is
the block loop itself.)
Since the trade is a property of each column rather than of the matrix, 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, Sum(len^2) / Sum(len) -- and columns
failing the test are replayed entry-wise. A column with no recorded blocks
takes the entry-wise path, so the decision doubles as its own storage and
needs no extra array. Factors of mixed structure, which are the common
case, then take the contiguous path on exactly the columns that benefit.
RUN_MIN_LEN is 16, a vectorization threshold: eight iterations of a
two-wide double-precision loop, comfortably past its prologue and epilogue,
whereas at two or four entries the setup dominates.
Both paths perform the same operations on the same values in the same
order, so the factorization is bit-identical whichever is chosen. Two
tests assert that against a fresh factorization on matrices sitting on
opposite sides of the decision, and a third pins that the decision does go
both ways, so neither path can rot untested.
Measured on 80 public problems (Netlib, Netlib-Kennington, Maros-Meszaros,
Mittelmann, SDPLIB, structured conic), interleaved against the schedule
replay alone:
67 small and medium problems: -29.2% total, median -10.0%, 54 faster,
12 unchanged, one slower by 2.3%
13 large problems (>10s), held out from the calibration of RUN_MIN_LEN:
-35.7% total, all 13 faster
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CqGdm2vZZ6HA5HT8DsRePa
adfcb7c to
096ae97
Compare
|
Superseded by a broader benchmark; see the upstream PR for the current state. Summary of what changed: Timing this on the 14 slowest corpus problems was a biased sample. On all 80 small/medium/large public problems (with the Cause: the block-wise update trades one indexed load per entry for bookkeeping per block, which a two- or three-entry block does not amortize. Fixed by deciding per column — a column's blocks are recorded only if the block containing a typical entry reaches 16 entries — which removes the regressions and, on mixed-structure problems, beats replaying everything block-wise ( Final: −29.2% across the 67 small/medium problems (one +2.3% outlier) and −35.7% across the 13 large ones (all faster), bit-identical throughout. |
1 participant
The observation
The dominant loop of the numeric factorization scatters into the sparse accumulator:
one indexed load and one indexed store per element — the compiler cannot vectorize it. But KKT factors from block-structured problems place most of their arithmetic in columns whose row indices are long blocks of consecutive integers: a set of dense-ish rows, which AMD orders last, produces a trailing trapezoid shared by thousands of columns.
Measured on the portfolio KKT systems (flop-weighted, i.e. each entry weighted by its column's height, since a column's update range is walked once per dependent column):
The change
Record each L column's decomposition into maximal consecutive-row runs alongside the schedule, and replay the update range run by run:
Same independent elementwise operations, same values, same order ⇒ bit-identical
L,D,Dinv, inertia and regularization decisions. No per-element index load, contiguous slices, so the compiler vectorizes. Factorization only — triangular solves untouched (see "what was tried and rejected"). Runs stored asu32pairs with the same overflow fallback as #6.No tunable constants. The run decomposition is exact and data-derived; the "≥ 8" above appears only in the measurement that motivated the change, never in the code.
Why this instead of a bordered/Schur factorization
The natural plan for this problem class was an OOPS-style block factorization exploiting diagonal
Pplus the dense border (Gondzio & Grothey, EJOR 181(3):1019–1029, 2007). Flop accounting first: AMD already achieves zero fill-in on these systems (nnz(L) = 592,846 vs nnz(K) = 611,653) and orders the dense border last, so a bespoke block factorization would perform the same ~42M flops. The opportunity was never fewer flops — it was executing them at vector rate instead of scalar-gather rate, which this does in ~70 lines instead of a new KKT solver. Generic supernodal codes (faer) miss it because the ~2000 tall columns have differing sparse upper patterns and so do not form supernodes; only their dense-border tails coincide, which is exactly what the run decomposition captures.Validation suites used throughout this series
All timing is on one Apple M-series core, single-threaded,
--releasewithdebug 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 clippyclean on thetouched 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, twolarge 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, soreading 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 classicsnb,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:
nbsolves 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_ip0546765.209042728 vs Clarabel 6765.2090428.
Tests specific to this PR
No new tests: #6's
test_refactor_matches_fresh_factor_exactlyand..._with_regularizationalready assert bitwise equality of a refactorization against a fresh factorization, which is exactly the invariant this change must preserve, and they cover the new code path unchanged.Performance: before → after
Interleaved against #6's branch, min of 3 rounds × 2 reps. Objectives bitwise-equal and iteration counts identical on every problem below.
Portfolio-rebalance set (9 problems)
Public corpus (7 problems)
The public-corpus spread is expected and honest: the gain tracks how much of a factor's work sits in long runs, and two problems regress slightly where it does not. Since results are bit-identical, a run-density gate could suppress those, but I did not add one — the threshold would be fitted to this corpus rather than derived, and the aggregate is clearly positive.
Status/objective regression gate
224-problem in-sample corpus: 0 status changes, 0 iteration changes, 0 objective changes (guaranteed by bit-identity, verified anyway).
What was tried and rejected
Extending the same run decomposition to the triangular solves (
_lsolve/_dltsolve, ~21% of wall). It is also bit-identical there — the runs partition each column in position order and are traversed in increasing position, so even theLᵀdot product accumulates in the original order — and it gave −2.2% on the portfolio set. But it was +5.4% slower on the public corpus (mitt_rail582+17%), so it is not included. The reason is a clean asymmetry: factor updates are flop-weighted, dominated by tall long-run columns, whereas the solves are entry-weighted (every column is visited exactly once), so plain average run length governs — 5.64 on the portfolio problems versus 1.37 on rail582, and below ~3–4 entries per run the per-run loop overhead exceeds the one indexed load it saves.🤖 Generated with Claude Code
https://claude.ai/code/session_01CqGdm2vZZ6HA5HT8DsRePa