perf(rescore): fit the folds one at a time, and the logistic regression in two parallel phases - #119
Merged
Merged
Conversation
`logreg_fit` was one serial loop over training rows per epoch, and serial all the way down: `z` is a left-fold into a single f64, which rustc can neither reassociate nor vectorise, so each of the d terms costs a full scalar-add latency. The only source of threads in the whole native rescorer was `percolator_lite`'s parallel map over `folds`, which defaults to 3, so the shipped default classifier ran ~2,000 full passes per fold (num_iter 10 x epochs 200) on three cores of whatever machine it was given. The epoch is now two phases. Phase 1 computes the per-row residual in parallel over rows; each row keeps its own strict ascending-j fold for `z`, so `z`, `p` and `err` are bit-for-bit unchanged. Phase 2 folds the residuals in: `grad[0]` stays one serial left-fold over the rows in ascending order, and `grad[j+1]` is an independent left-fold per column, so it parallelises over COLUMN BANDS of 16 (64 B = one cache line of f32, so the second pass moves its payload rather than 16x it). Every column sees the same terms in the same order. No sum is reordered anywhere. The band accumulators are a local array copied back once: accumulating straight into the `grad` slice shares cache lines between neighbouring bands, and that false sharing alone made the first version 3.4x SLOWER than the serial fit. Also parallelised, all per-row maps with no reduction in them: the standardised training matrix fill (`std_row_into` now writes a caller slice so the whole matrix can be filled by `par_chunks_mut`), the per-iteration rescoring of the training rows, and the held-out test-fold scoring. Numbers, `logreg_fit_epoch_cost` (new, #[ignore]d; both arms share the prebuilt row pointers and each allocates its own per-epoch working set inside the timer). i9-13900KS, 32 rayon threads, release, 387 features: 5.3x at 4,000 rows, 4.4x at 16,000, 2.7-4.3x at 60,000, 3.2x at 300,000 (serial 104 ms/epoch, two-phase 32). The gain falls as the slice grows because both passes become DRAM streams: 2 x 464 MB in 32 ms is ~29 GB/s, about all this desktop's two channels have. The survey projected ~12x from a 64-core server model; that is NOT resolvable on this machine and I am not claiming it. What the fixture does corroborate is the size of the problem: 104 ms/epoch at 300,000 x 387 scales to ~415 ms at a real ~1.2M-row fold, against the survey's 466 ms estimate, i.e. ~14 min per fold. Equality: byte-identical. No float sum changes order in either phase. `percolator_lite_reference` in the test module now calls a transcribed `logreg_fit_reference` (the old one-pass body) rather than the shared kernel, so the existing end-to-end test pins the OLD behaviour instead of comparing the new code with itself. New `two_phase_gradient_fits_bit_identical_weights` compares fitted weights bit for bit across widths that straddle the band (1, 5, 15, 16, 17, 32, 33, 387), row counts 1..301, an empty row set, a `y` shorter than `rows` (the old `zip` truncation) and zero epochs. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…at once `percolator_lite` mapped over folds with `into_par_iter()`, which was the only source of threads in the native rescorer, and the price was that every fold's standardised training copy was live simultaneously. On a six-run Astral pool (3,133,636 PSMs x 387 features) that is the 4.85 GB matrix plus three 3.23 GB copies = 14.55 GB; at 11.6M PSMs, 17.95 + 3 x 11.97 = 53.8 GB. Sequentially it is the matrix plus one copy: 8.08 GB and 29.9 GB, saving 6.47 and 23.9 GB. `train_idx`/`test_idx`, `mean`/`std`, `train_scores` and the `rows` pointer vector become single-copy at the same time (~280 MB at the Astral pool, ~1.0 GB at 11.6M), and `fold_of` and the two index vectors are now `u32` rather than `usize` (~60 MB / ~220 MB more). This is only affordable because the previous commit made one fold saturate the machine; before it, this would have been 3x the wall. What it still costs is the serial remainder that no longer overlaps across folds -- the tied-block walk in `target_decoy_q`, the positive-set build, the bias fold -- which is seconds per fold against a fit measured in minutes. Second change, the same call path: `fdr::target_decoy_q_split(&[f64], &[bool])` is the same kernel from two parallel columns instead of a slice of pairs. `percolator_lite` held a `Vec<(f64, bool)>` allocated OUTSIDE the `num_iter` loop and merely refilled inside it, so 16 bytes per training row stayed resident for the whole fit purely to pair two columns the fold already had: 3 x 2.089M x 16 = 100 MB concurrent at the Astral pool, 371 MB at 11.6M, now a 2.1 MB decoy mask gathered once. `entrapment_q` already took its columns separately, so the two kernels are now consistent. `target_decoy_q` keeps its signature and is a one-line wrapper over the shared body, so no caller outside this subsystem changes. Equality: byte-identical. Folds are genuinely independent -- each fits its own scaler on its own train rows and writes only its own disjoint test indices -- so the merge is a scatter, not a reduction, and the order the folds are consumed in cannot move a value; scattering inside the loop is the same scatter. The split q form reaches the same kernel with the same `(score, is_decoy)` per row. `total_d` is now counted off the sorted records instead of the input, which is the same integer over the same multiset. Validation: SMOKE_OK, 144 assertions. `cargo test --workspace` 314+ pass. `flat_training_matrix_scores_bit_identically_to_the_per_row_layout` now compares against a reference that is the complete OLD implementation -- parallel folds, `usize` bookkeeping, `Vec<Vec<f32>>` training rows, the paired-staging q call and the one-pass fitter -- via a transcribed `fit_standardizer_reference`, so it pins the behaviour being replaced rather than the shape replacing it. New `the_split_and_paired_q_forms_agree_bit_for_bit` and a should_panic test for the mismatched-column guard. Not done, blocked: the error text at stages/rescore.rs:250-256 still tells the user the native peak is "roughly (1 + folds) times" the matrix. After this it is 2x regardless of `folds`. That file is outside this task's file list. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…ed to `target_decoy_q` was rewritten to rank one 16-byte sortable record per row and to monotonize inside the same backward pass that computes each block's FDR. `entrapment_q` was left behind on exactly the form that rewrite documents as having been removed: a `key` column, a `usize` index permutation sorted with a comparator that dereferences it, a forward tied-block walk into a third n-vector, and a fourth for `q`. Same treatment, line for line. Four n-vectors (32 bytes per row) become two (`ranked` 16 + `q` 8 = 24): 75 against 100 MB at 3.13M rows, 278 against 371 MB at 11.6M. The serial `sort_by` over `usize` with two random loads per comparison becomes `par_sort_unstable_by` over a cache-resident record, and the tied-block walk is a sequential scan instead of two indirections per row. This kernel runs five or more times per entrapment rescore -- pooled PSM q, once per source for `run_psm_q`, three times through `grouped_q` -- over the whole scored table. Size, stated honestly: the memory term is ~25 MB / ~93 MB, which on its own is immaterial. The sort is the real term and I cannot put a wall-clock number on it here, because no entrapment-arm timing exists for this tree and the smoke fixture does not reach the kernel at all (`QMode::Entrapment` only). What earns the change is that FDR validity is a first-class objective, the rewrite is a transcription of one that already exists and is already proven in-repo, and two divergent copies of one kernel is how the next reader picks the wrong one. Equality: byte-identical. The stable sort on key-descending over an ascending index vector is the same permutation as an unstable sort on (key desc, row asc), and no two records compare equal. In the fused pass, `ne` and `nr` are integer counts at or above the block -- totals minus what is strictly below -- which is exactly what the forward walk had accumulated by the end of that block, so `(ratio * ne + 1) / max(1, nr)` is the same float from the same integers and no float reordering occurs at all. `qmin` accumulates worst-to-best in both, and one min per block rather than one per rank is the same value because every rank in a block held the same `f`. Entrapment still wins over real where a row claims both, which is the precedence the forward walk's `else if` had. Tests: `entrapment_q_reference` transcribes the old body and `entrapment_q_is_bit_identical_to_the_previous_kernel` compares bit patterns over heavy ties, signed zeros, NaN and both infinities, rows that are neither class, rows that are both, all-entrapment and all-real populations where `max(1, nr)` bites, the empty input, and three library-size ratios. Validation across all three commits: cargo fmt --check, clippy -D warnings, cargo test --workspace (315 pass), and SMOKE_OK with all 104 parquet and TSV artifacts hashing identically to the ones the base commit a816768 produced. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…re they lose The previous two commits only make sense as a pair, and the pair is not a free win: the old arrangement got `folds`-way parallelism for nothing, so running folds sequentially is faster only where the per-fold speedup beats `folds`. `percolator_lite_fold_cost` (#[ignore]d) measures exactly that, both arms handed the same `FeatureMatrix` and each allocating its own standardised copies, index vectors and per-epoch working set inside the timer; it also asserts the two arms agree bit for bit, so it is a second equality check. Measured, i9-13900KS, 32 rayon threads, release, 387 features, 200,000 PSMs (training slices 206-310 MB, nothing cache-resident): 1.51x at the default 3 folds, 0.91x at 5. So at the shipped default this is a wall-clock gain as well as a 1.8x smaller peak; under the opt-in `folds: 5` sensitivity recipe it is wall parity at best on a desktop and the memory is the whole reason. The same fixture at 30,000 PSMs reads 5.00x / 0.97x / 0.24x at 2 / 3 / 5 folds. That 0.24x is NOT a fold-count effect: the training slices there are 23 / 31 / 37 MB against this chip's 36 MB L3, so the largest one thrashes and the new fitter, which is bandwidth-bound, pays for it while the old serial one hides the misses behind its dependency chain. It is recorded because it is the sort of number that gets quoted out of context, and because it is the reason the fixture worth believing is the large one. Both comments in `percolator_lite` and `logreg_fit` now carry these numbers rather than a projection. No production code changes in this commit. Equality: unchanged, byte-identical; this commit adds a test and comments only. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…nds check An adversarial review extracted `fdr.rs` and `rescoring.rs` into a standalone crate, transcribed the PARENT commit's production `percolator_lite`, and measured against that. Its verdict: the equality claims all hold and are strong (bit identical at 50,000x387, 200,003x97 and 1,000,001x8, across rayon pools of 1/2/3/5/13 threads; `target_decoy_q_split` and `entrapment_q` likewise), the memory claim holds (300,000x387, folds=3 1,385 -> 767 MB, folds=5 2,308 -> 829 MB), and the measurement claims do not. It was right about the measurement and about a behavioural regression. WHAT WAS WRONG WITH THE MEASUREMENT `percolator_lite_fold_cost` timed `percolator_lite_reference`, which is the layout from TWO commits back: it rebuilds `Vec<Vec<f32>>` per training row, which the parent had already replaced with a flat `xtr`. That arm is now measurable (`MUMDIA_BENCH_BASELINE=per_row`) and costs 1.03x the parent's wall at 3 folds but 1.56x at 5 (211.6 s against 135.8 s at 200,000x387x10 iters), so it inflated exactly the configuration where this change is weakest. It also perturbs the arm it is compared against: at 5 folds it keeps ~800,000 live heap blocks, and the interleaved sequential arm then varied 124-212 s where it is 104-105 s against the parent. Second, `MUMDIA_BENCH_ITERS` defaulted to 2 against a shipped `rescore.num_iter` of 10. Third, every number was a single shot on a hybrid desktop CLAUDE.md warns about. CORRECTED NUMBERS The benchmark now times `percolator_lite_parent`, a line-for-line transcription of a816768's production `percolator_lite` (parallel fold map, `usize` bookkeeping, flat `xtr` filled serially, the `Vec<(f64, bool)>` staging buffer outside the `num_iter` loop, the serial one-pass fitter, `score_std_row` for the test fold). `MUMDIA_BENCH_ITERS` defaults to 10. `MUMDIA_BENCH_REPS` (default 3) repeats with the arms interleaved and the order alternating; the summary is best-of with the spread printed. i9-13900KS, 32 rayon threads, release, 387 features, 200,000 PSMs: folds num_iter 2 num_iter 10 (shipped) 3 1.91x 2.03x 5 1.31x 1.30x and at 30,000 PSMs, num_iter 10: 6.16x / 3.82x / 2.17x at 2 / 3 / 5 folds (at num_iter 2: 5.82x / 3.46x / 1.89x). Repeat spread 1.00-1.16x on the parent arm. The implied per-fold speedup is consistent (107.4 s of parent against 3 x 17.6 s, 135.8 against 5 x 20.9: about 6.1x and 6.5x) and beats both fold counts. So the honest answer is a win everywhere measured, not the "1.51x at 3 folds, 0.91x at 5" the commit claimed, and not parity. WITHDRAWN: that pair, and the "0.24x at 5 folds" on the 30,000-PSM fixture along with the L3-thrashing story invented to explain it -- repeated and interleaved, that exact configuration reads 1.89x, so the story was explaining a noise artifact. The review's mechanism for `num_iter` is real but small here: the per-fold serial remainder costs relatively more when the fit is short, worth 6% at 3 folds and nothing at 5, not a change of sign. `logreg_fit_epoch_cost` also gained repeats, and its old numbers were also single-shot low: 6.1x at 4,000 rows, 8.5x at 16,000, 5.1x at 60,000, 5.1x at 300,000 (serial 103 ms/epoch, two-phase 20) and 4.9x at 1,000,000, against 5.3x / 4.4x / 2.7-4.3x / 3.2x before. The serial arm reproduces to 1% (104 against 103 ms/epoch), so it was the two-phase arm that was mistimed. The bandwidth ceiling is 53-55 GB/s on this desktop, measured twice, not the 29 GB/s recorded before. THE BEHAVIOURAL FIX `percolator_lite` drove the fold loop off `fold_of.iter().enumerate()` where the parent used `(0..n).filter(|&i| fold_of[i] != test_fold)`. With a `fold_key` shorter than the matrix the parent panicked; the new code returned, leaving every trailing PSM at its unrescored `init_score`, from where it enters the q population indistinguishable from a real score. `assert_eq!` on `fold_key`, `is_decoy` and `init_score` against `n`, with a test for two of them. `logreg_fit`'s `rows.len().min(y.len())` truncation is likewise now an `assert_eq!`: the one-pass form fitted on the shorter prefix while dividing the gradient by `rows.len()`, and callers push the two in lockstep. OTHER CORRECTIONS IN THESE TWO FILES - the old peak was `folds` x the matrix, not "(1 + folds) x": a training copy is (folds-1)/folds of the matrix, so `folds` of them is folds-1 matrices. The new peak is 1 + (folds-1)/folds, so 1.67x at the default and 1.80x at folds 5, approaching 2x and never above it, and no longer growing with `folds`; - the GRAD_BAND rationale was wrong at d = 387. 16 f32 = 64 B BOUNDS the over-fetch at 2x, it does not remove it: the row stride is 1548 B, 1548 mod 64 = 12, so one row in 16 starts line-aligned and the rest straddle two lines, ~1.94x payload. The real trade is over-fetch (~1 + 16/B) against task count (ceil(d/B), already 25 against 32 threads at B = 16). Swept at 300,000x387x5 epochs, best of 3, with the serial control at 0.500-0.517 s in all four builds: 0.144 s at B = 8, 0.101 at 16, 0.123 at 32, 0.160 at 64. 16 is the measured optimum, so the constant stands on the measurement rather than on the reasoning; - the per-epoch cost is 2.4x the serial form's DRAM traffic, not 2x: phase 2b re-reads the `rows` pointer vector and `errs` once per band, 25 times per epoch at d = 387, which is 1.25 GB on top of two 3.23 GB matrix passes at a 2.09M-row fold; - `percolator_lite_fold_cost` now states what it does NOT measure: it is a kernel on a synthetic population, not `mumdia rescore` on a real competed table, so no share-of-stage-runtime claim can rest on it. REJECTED - the review's finding 9 calls the two q entry points a divergence hazard. They are not: `target_decoy_q` and `target_decoy_q_split` are both one-line wrappers over `target_decoy_q_core`, so there is no second estimator to drift. The five remaining pair-form call sites are a saving left on the table (186 MB at 11.6M rows) and are now named in `fdr.rs` with that distinction stated; they are outside this file set; - GRAD_BAND is not changed. The reasoning was wrong, the constant is right. NOT REACHABLE FROM THIS FILE SET (five stale "(1 + folds)" statements, all wrong by the arithmetic above; two of them also still describe the matrix as `Vec<Vec<f64>>`, which it has not been since it became a flat f32 `FeatureMatrix`): stages/rescore.rs:252, mumdia-core/src/config.rs:1572, docs/11_compete_rescore_fdr.md:301, docs/24_config_reference.md:319, CLAUDE.md:405. Equality: unchanged, byte-identical. The equality test now pins against BOTH older layouts (parent and per-row), and the benchmark asserts bit equality against the parent on every repeat, which it did at 200,000x387 and 30,000x387 over six configurations. SMOKE_OK, 144 assertions, `peptides.tsv` and `proteins.tsv` hashes unchanged across the bounds-check restoration. `cargo fmt --check`, `cargo clippy --workspace --all-targets -D warnings`, `cargo test --workspace` 318+ pass. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`rescoring.rs` now fits the folds one at a time, so only one standardised training copy is live. Five places still described the parallel-fold peak as `(1 + folds)x` the matrix, and four of them additionally described the matrix as `Vec<Vec<f64>>`, which it has not been since it became a flat f32 `FeatureMatrix`. Both figures were wrong in the same direction. A training copy is `(folds - 1) / folds` of the matrix, so `folds` of them is `folds - 1` matrices and the OLD peak was `folds x` the matrix, not `(1 + folds) x`. The new peak is `1 + (folds - 1) / folds`: 1.67x at the default 3 folds, 1.80x at 5, approaching 2x and never above it. Measured at 300,000 x 387: folds 3 1,385 -> 767 MB, folds 5 2,308 -> 829 MB. The one that misled an operator is the `rescore.max_feature_matrix_gib` bail text, which read 1.8x too high at the default and 2.8x at 5 folds. Corrected there, in the config doc it is generated from, in docs/11, in the regenerated docs/24 and in CLAUDE.md. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The chunked `write_scored_table` from #118 moved every environment read below it in `rescore.rs` by ten lines, and `docs/24_config_reference.md` records each read by source line. No setting, default or variable changed. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
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
1 participant
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.
PR #117 was merged at 09:28 UTC, before the rescoring commits and their fix round were pushed to the same branch. These are the six that were left behind.
What is in it
perf(rescore): fit the native logistic regression in two parallel phases, and fit the cross-validation folds one at a time rather than all at once.perf(fdr): putentrapment_qon the layout its sibling already uses.fix(rescore): the correction round after an adversarial review of the above.docs: the(1 + folds)xmemory statements the change invalidates.The benchmark was wrong, and the corrected answer is better than the original claim
The original commit measured against a
Vec<Vec<f32>>straw man rather than the production code it replaced, which inflated exactly the configuration where the change is weakest. Rebuilt aspercolator_lite_parent, a line-for-line transcription of the parent commit's productionpercolator_lite, withMUMDIA_BENCH_ITERSdefaulting to the shippedrescore.num_iterof 10 and three interleaved repeats reported best-of:At 30,000 PSMs,
num_iter10: 6.16x / 3.82x / 2.17x at 2 / 3 / 5 folds.Three numbers are withdrawn in the code: the committed "1.51x at 3 folds, 0.91x at 5", and the "0.24x at 5 folds" on the 30,000-PSM fixture together with the cache-thrashing story invented to explain it — repeated and interleaved, that configuration reads 1.89x.
MUMDIA_BENCH_BASELINE=per_rowmakes the bad baseline measurable: it costs 1.03x the parent's wall at 3 folds but 1.56x at 5.Memory
A training copy is
(folds - 1) / foldsof the matrix, sofoldsof them isfolds - 1matrices and the OLD peak wasfolds xthe matrix, not the(1 + folds) xfive places said. Sequentially it is1 + (folds - 1) / folds: 1.67x at the default 3 folds, 1.80x at 5, approaching 2x and never above it, and no longer growing withfolds. Verified at 300,000 x 387: folds 3 1,385 -> 767 MB, folds 5 2,308 -> 829 MB.The statement that actually misled an operator is the
rescore.max_feature_matrix_gibbail text, which read 1.8x too high at the default and 2.8x at 5 folds. Corrected there, in the config doc it is generated from, indocs/11, in the regenerateddocs/24and in CLAUDE.md.Also restored
Bounds checks that the perf commit had dropped:
assert_eq!onfold_key,is_decoyandinit_scoreagainstninpercolator_lite(the silent-partial-rescore path), andassert_eq!(rows.len(), y.len())inlogreg_fitreplacing amintruncation. Two newshould_panictests.GRAD_BAND = 16was re-derived rather than changed: at 300,000x387x5 epochs, best of 3, 0.144 / 0.101 / 0.123 / 0.160 s at B = 8 / 16 / 32 / 64.Validation
cargo fmt --check,cargo clippy --workspace --all-targets -- -D warnings, rustdoc with-D warnings,cargo test --workspace(361 + 5 + 2 + 6 + 49 + 32 pass), andSMOKE_OKwithpeptides.tsvf0b5dc38... andproteins.tsv5a65d304... unchanged.🤖 Generated with Claude Code