compute the two residual products in one pass over A - #11
Open
batterseapower wants to merge 1 commit into
Open
Conversation
The interior-point residual update needs both sparse products on every
iteration:
rx_inf = -A' z
rz_inf = s + A x
Each is a separate gemv, and each streams all of A's row indices and
values. Since both traverse the same nonzeros, one pass can serve both:
for each nonzero (i,j,v) the same loads of v and rowval[k] feed the
column accumulation of A'z and the scatter of A x.
Both existing kernels visit columns in order and entries within a column
in order, and the fused loop does the same, so each output accumulates
in exactly the sequence it did before and the results are bit-identical.
A test asserts that against the two separate calls.
Measured (interleaved A/B, min of 3 rounds x 2 reps, objectives
bitwise-equal throughout):
nine portfolio-rebalance SOCPs: -3.7% total wall clock, every problem
faster (-2.5% to -5.8%)
seven Netlib-Kennington / Mittelmann / Maros-Meszaros / conic
problems: +0.6% total, between -2.0% and +1.6%
The difference is cache behaviour, not arithmetic: fusing trades one
traversal of A for two indirect accesses to different vectors inside the
same loop, which pays when z and rz_inf stay resident and is a wash when
they do not. It is bit-identical either way, so the choice is purely
about speed.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CqGdm2vZZ6HA5HT8DsRePa
This was referenced Jul 21, 2026
Owner
Author
|
Held back from the upstream submission. Public profiling puts the two products at roughly 1% of runtime on public problems (0.73% on |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Independent of the rest of the series — this touches
residuals.rsand the CSC kernels, which no other open PR modifies, so it can land in any order.What
The interior-point residual update needs both sparse products on every iteration:
They are two separate
gemvcalls, and each one streams all of A's row indices and values. Since both traverse the same nonzeros, one pass can serve both: for each nonzero(i,j,v)the same loads ofvandrowval[k]feed the column accumulation ofA'zand the scatter ofA x.On the portfolio problems
nnz(A)is 271k–538k and this runs once per iteration, so it was a measurable share of runtime — profiling with#[inline(never)]put_csc_axpby_Tat 5.95% and_csc_axpby_Nat 1.31% of wall clock, together the largest remaining cost outside the factorization, the refinement residual and the triangular solves.Why it is correct
Both existing kernels visit columns in order and entries within a column in order. The fused loop does the same, so each output accumulates in exactly the sequence it did before — the results are bit-identical, not merely equivalent.
test_csc_neg_At_and_A_matches_separate_gemvsasserts that against the two separate calls on a matrix with empty and multi-entry columns.No tunable constants.
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
test_csc_neg_At_and_A_matches_separate_gemvs— the fused kernel vs_csc_axpby_T+_csc_axpby_N, compared bitwise on both outputs.Performance: before → after
Interleaved, min of 3 rounds × 2 reps. Objectives bitwise-equal, and statuses and iteration counts identical, on every problem — the solver follows exactly the same trajectory.
Portfolio-rebalance set (9 problems)
Public corpus
Measured twice, on a narrow and then a wider set, because the first result was small enough to be sampling noise.
14 timing-relevant problems, 3 rounds × 2 reps (the trustworthy measurement):
7-problem subset, 2 rounds × 2 reps (reported for completeness): 16.153s → 16.255s, +0.6% — three faster, four marginally slower. On the wider set with more rounds the aggregate is a small win rather than a small loss, which is the honest summary: the effect on the public corpus is real but close to noise, roughly −1%, while on the portfolio class it is a consistent −3.7%.
The mechanism is cache behaviour, not arithmetic — fusing trades one traversal of A for two indirect accesses to different vectors inside the same loop, which pays when
zandrz_infstay resident. Since the result is bit-identical either way, there is no numerical risk in either direction, and I did not add a size heuristic to choose between them because any threshold would be fitted to these corpora rather than derived.Status/objective regression gate
224-problem in-sample corpus: 0 status changes, 0 iteration changes, 0 objective changes — guaranteed by bit-identity, verified anyway.
🤖 Generated with Claude Code
https://claude.ai/code/session_01CqGdm2vZZ6HA5HT8DsRePa