qdldl: form KKT refinement residuals from a both-triangles matrix copy

[batterseapower]

Two commits, reviewable independently. Depends on nothing else; touches src/qdldl/ and the direct-LDL KKT solver only.

Motivation

After the factorization, the iterative-refinement residual e = b − Ax is the largest cost in the solve. 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 (_csc_symv_unsafe)	7.5%	7.5%
triangular solves	2.9%	15.1%

The kernel is why: exploiting symmetry from the upper triangle alone requires two scattered read-modify-writes per off-diagonal nonzero, into y[row] and y[col] — the worst available memory pattern, and unvectorizable.

Separately, each refinement pass currently goes through QDLDLFactorisation::solve, which copies the right-hand side and applies the fill-reducing permutation and its inverse per backsolve.

Commit 1 — run the refinement loop inside the backend

Adds QDLDLFactorisation::solve_refined: the same refinement algorithm — same residual definition, same stopping rules (tolerance, insufficient-improvement ratio with keep-if-improved, iteration cap), driven by the existing iterative_refinement_* settings — performed entirely in the internally permuted coordinates. The permutation is applied once to the right-hand side and once to the returned solution instead of once per backsolve.

A new DirectLDLSolver::solve_refined trait method defaults to "unsupported", so the faer and Pardiso backends keep the caller-side path untouched and only QDLDL takes the new one.

Supporting change: the static-regularization shift is restored in the backend's internal copy after each refactorization as well as in the caller's copy, so that copy holds the unregularized matrix between refactorizations, which is what refinement must target. Safe for all backends — the shift is recomputed from the restored diagonal before the next refactorization.

On its own this commit is performance-neutral on public problems (measured −0.0%); its value is putting the residual where commit 2 can replace it.

Commit 2 — form the residual from a both-triangles copy

Keep a full both-triangles copy of the internally permuted matrix. Because it is symmetric, its CSC arrays read as CSR describe the same matrix — "column i" is exactly row i — so the residual becomes a sequence of independent sparse dot products: one streamed pass over values, gathered reads of x, an accumulator in a register, no scatter.

Four details, each with a reason:

• The copy is built on first use and refreshed only after the internal matrix changes, i.e. one O(nnz) pass per refactorization. The value map runs from each nonzero of the copy to its source, which makes that refresh a gather with sequential stores.
• Each row accumulates into four partial sums reduced pairwise. One accumulator serializes the row on the latency of a single dependent add; independent partial sums also carry a tighter error bound, growing like n/k + k for k accumulators rather than n (Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed. 2002, §4.2, on blocked and pairwise summation). Measured: one accumulator is materially slower, eight is indistinguishable from four and needs more code.
• A NaN residual is reported rather than masked. A running maximum cannot detect one, because IEEE maxNum returns the non-NaN operand — which would leave the norm finite and let a non-finite search direction be accepted as converged. Not hypothetical: two MIPLIB relaxations in the gate below exercise it.
• It is an addition, not a replacement: if the pattern exceeds the u32 value map, refinement falls back to the triangular product.

Correctness

Structural correctness is pinned exactly, separating an error in the copy or its value map from mere rounding: test_symmetric_copy_exact_on_integer_data compares the new residual against the triangular product bitwise on exactly-representable integer data, where every product and sum is exact. test_symmetric_copy_matches_triangular_symv repeats it to tolerance on a pseudo-random quasidefinite matrix. test_sym_residual_reports_nan covers the NaN path, which no finite-data test can reach. test_solve_refined pins the refinement semantics: it asserts the loop refines against the values currently in the internal workspace rather than the ones that were factored, by rescaling the matrix without refactoring and requiring convergence to the rescaled system's solution.

Full suite green (cargo test, 20 binaries); cargo clippy clean.

On real data the residual differs from the previous formulation by rounding. Both are backward-stable evaluations of the same product, but a refinement residual is a difference of nearly equal quantities — as the solve converges, b and Ax cancel to many digits — so the two summation orders differ at roughly the level of that cancellation. On problems already at their double-precision conditioning floor this reshuffles the last few iterates, which is what the gates below bound.

For the record, a more accurate residual was tried and rejected: a compensated dot product (two-product via FMA plus two-sum, Ogita, Rump & Oishi, SIAM J. Sci. Comput. 26(6):1955–1988, 2005, Alg. 5.3), motivated by Higham ch. 12's requirement that refinement residuals be formed to higher precision. It was 16% slower — the compensation chain is serially dependent, so the loop becomes latency- rather than bandwidth-bound — and it did not stabilise termination statuses.

Performance

Interleaved A/B — the binaries run alternately problem-by-problem within each round, so machine drift affects both equally — minimum over rounds, single-threaded release build. 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
76.1s → 73.2s, −3.8%	−5.7%	53	8	6

13 large problems (>10s), 1 rep: 1451.1s → 1395.7s, −3.8%; 8 faster, 4 flat, 1 slower. Largest gains sdplib_qap9 −20.3% (31 → 25 iterations), sdplib_qap8 −11.3%, sdplib_arch0 −9.5%, kennington_osa_60 −9.0%.

A consistent but modest public gain: about 4% of total solve time, 6% at the median, on top of which refinement-dominated problems gain considerably more.

The regressions, and why there is no gate for them

Seven of the 80 get slower: conic_logreg_winewhite_exp +25.0%, sdplib_qap7 +4.8%, sdplib_truss5 +4.6%, sdplib_control5 +4.6%, conic_pow3d_p25_winewhite_pow +4.1%, conic_logreg_winewhite2000_exp +3.6%, netlib_pilot +2.7%.

These are almost all iteration-count changes rather than kernel slowdowns — conic_logreg_winewhite_exp 23 → 30 iterations, sdplib_qap7 26 → 29, sdplib_control5 34 → 38 — and the same mechanism produces the largest gain in the other direction, sdplib_qap9 at 31 → 25. Perturbing the summation order of the residual moves marginal problems onto different iterate paths, and that cuts both ways. It cannot be predicted from a structural property of the matrix, so unlike the block-kernel decision in the companion factorization PR there is no gate to be had; the alternatives are this change as it stands, or not changing the residual.

Regression gates

224 public problems (101 Maros–Mészáros QPs, 85 Netlib LPs, 24 structured conic, 12 Netlib-Kennington, 2 Mittelmann): 15 change iteration count; one changes status — mm_qsierra, Solved → AlmostSolved, with a bit-identical objective, finishing at dual residual 4.2e-8 against the 1e-8 tolerance and a duality gap of 1.5e-13, i.e. marginal by 4× on one criterion while returning the same answer to every printed digit. Three same-status objective differences exceed 1e-6 relative (mm_qbeaconf 1.2e-6, netlib_forplan 7.8e-6, netlib_pilot_ja 1.7e-6), all on problems terminating AlmostSolved in both versions.

Two further suites, fetched fresh and used only as gates:

MIPLIB 2017 benchmark set as LP relaxations, 45 instances (the MPS reader ignores INTORG/INTEND, so reading a MIPLIB file yields exactly the continuous relaxation), 90s limit: 1 status better (comp21_2idx AlmostSolved → Solved), 0 worse, 0 objective disagreements > 1e-6.

CBLIB, 84 instances sampled across families, including nb, nb_L1, nb_L2, nql30/60/180, qssp30/60, sched_*: 2 status better (clay0304h, sched_100_50_orig), 1 worse (sched_100_100_orig, Solved → AlmostSolved, objective differing 1.1e-5), total time −7.8%. That regression is in an ill-conditioned family — CBLIB ships _scaled variants of exactly these problems because the originals are badly scaled, and on sched_100_100_scaled both versions report Solved and differ by 9.8e-6. All three CBLIB objective disagreements are in that family.

Conversion of the CBLIB instances was validated independently: nb solves to −5.0703094644e-2, matching its published DIMACS optimum (−0.05070309), and HiGHS agrees on the LP-only conversions (gen_ip054: 6765.209042728 versus 6765.2090428).

Cost

Memory: the copy roughly doubles the stored matrix, plus a u32 per nonzero of the copy for the value map.

[batterseapower]

Hi: #229 #230 and #231 together cut clarabel runtime by 26% on public problems and 42% on a private benchmark I care about. Most of the work was done by Fable 5 and Opus 4.8. LMK if any questions.

[batterseapower]

Updated the benchmark section. The timings I first quoted came from a 14-problem subset chosen as the slowest problems in my corpus, which is a biased sample; the numbers above are now from every small-, medium- and large-band problem in it (80 problems, six families, with the sdp feature enabled so SDPLIB is actually included — it was silently skipped before).

The headline is more modest than what I first posted, and the regressions are now visible: −3.8% of total solve time on both the 67 small/medium and the 13 large problems, with a −5.7% median, and 7 problems across the 80 getting slower.

Worth being explicit about the character of those regressions, because it bears on whether they are fixable: they are almost all iteration-count changes rather than kernel slowdowns — conic_logreg_winewhite_exp 23 → 30 iterations, sdplib_qap7 26 → 29, sdplib_control5 34 → 38 — and the same mechanism produces the biggest win in the other direction, sdplib_qap9 at 31 → 25 iterations and −20.3%. Changing the summation order of the residual moves marginal problems onto different iterate paths. That is not something a structural test on the matrix can predict, so unlike the block-kernel decision in #230 there is no gate to be had here; the alternatives are the change as it stands, or not changing the residual.

If the numerical-order change is unwelcome for a ~4% aggregate gain, this one is much more marginal than #230 and I would understand declining it. #230 is independent of it, touches no file in common, and is bit-identical.