perf(preprocessing): whiten via the covariance matrix in Whitener::pca

[mysma-9403]

WhiteningMethod::Pca takes the SVD of the whole centered (nsamples x nfeatures) data matrix, while the Zca and Cholesky arms immediately next to it first form the (nfeatures x nfeatures) covariance and decompose that instead.

For PCA whitening the two are equivalent: the right singular vectors of the centered data are the eigenvectors of its covariance, and the singular values relate to the eigenvalues by lambda = s² / (nsamples - 1). Forming the covariance first turns an SVD whose cost grows with the sample count into a fixed-size one, leaving a single GEMM as the only work proportional to nsamples.

Results

n	d	before	after	speedup
2 000	16	4.11 ms	0.15 ms	27×
8 000	32	83.28 ms	1.19 ms	70×
32 000	64	2051.71 ms	13.31 ms	154×
64 000	64	5259.08 ms	33.16 ms	159×
8 000	128	768.69 ms	22.92 ms	34×

Accuracy, measured in the same run (W is the whitening matrix, Y the whitened data):

n	d	max‖WᵀW − Σ⁻¹‖ before → after	max‖cov(Y) − I‖ before → after
2 000	16	7.37e-18 → 6.07e-18	3.55e-15 → 2.44e-15
32 000	64	1.26e-17 → 6.07e-18	6.66e-15 → 5.55e-15
64 000	64	2.52e-17 → 6.51e-18	1.24e-14 → 2.89e-15
8 000	128	1.87e-17 → 8.24e-18	6.00e-15 → 7.11e-15

Worth being explicit about the trade-off rather than reading too much into those numbers: forming XᵀX squares the condition number, so on badly conditioned inputs this formulation is in principle less accurate than an SVD of X, even though it comes out slightly ahead on the well-conditioned uniform-random data above. Zca and Cholesky already accept exactly that trade-off, so this makes Pca consistent with them rather than introducing a new compromise. If you would rather keep the data-matrix SVD for Pca specifically on numerical grounds, that is a reasonable call and I am happy to close this.

Details

• The epsilon floor is applied to the reconstructed singular values rather than directly to the eigenvalues, so max(s, 1e-8) keeps the meaning it had before.
• Eigenvalues are clamped at zero before sqrt, because rounding can push a numerically-zero one slightly negative and produce a NaN.
• Singular vectors are sign-arbitrary, so individual rows of the whitening matrix may come out negated relative to before. This does not affect the whitening property, nor WᵀW.

Fewer samples than features

With nsamples <= nfeatures the covariance is rank deficient, and the SVD of the data matrix returns a differently shaped factor — so switching formulations there would change the shape of the whitening matrix. That case is routed to the original code path, and two new tests pin the shape:

• test_pca_matrix_more_features_than_samples (20 × 50)
• test_pca_matrix_square_input (16 × 16)

Which shape that is turns out to depend on the backend, which I had missed on the first push — the first of those tests asserted only the default backend's answer and failed the BLAS CI jobs. linfa_linalg::svd is a compact SVD and returns the (nsamples, nfeatures) factor; ndarray_linalg::svd is a full one and returns (nfeatures, nfeatures). The test now asserts both. This split is not introduced here: the nsamples <= nfeatures arm is byte-for-byte the code that was there before, so both shapes are exactly what master produces today. The new covariance path is unaffected either way, since it decomposes a square matrix, on which the two agree.

While confirming that, one pre-existing wrinkle became visible, mentioned only for the record — it is untouched by this PR and I have not tried to fix it here. Under BLAS in that same regime, s has min(nsamples, nfeatures) entries while v_t has nfeatures rows, so the scaling loop leaves the trailing rows unscaled. They span the null space of the centered data, so the corresponding output columns come out identically zero rather than unit variance. The compact SVD has no trailing rows and so no such columns. Neither backend can genuinely whiten a rank-deficient input, so this is a difference in how the degenerate case is presented rather than in whether it works. Happy to open a separate issue if that is worth tracking.

Benchmark context

1. Plugged in (AC power, battery 100 %, charged)
2. Power saving mode: off (lowpowermode 0)
3. Machine otherwise idle
4. Not thermally throttled
5. MacBookPro15,3 — Intel Core i9-9980HK @ 2.40 GHz, 8 cores / 16 threads, 32 GB RAM, macOS 15.7.7, rustc 1.93.1
6. Default (pure-Rust linfa-linalg) backend, best of 3 Whitener::pca().fit(..) calls, measured with a throwaway harness that is not part of this PR; A/B taken on this branch with and without the change so the binary shape is identical

Checks

• cargo fmt --all -- --check — clean
• cargo clippy --workspace --all-targets -- -D warnings — clean
• cargo test --release --workspace — no failures
• cargo test --release -p linfa-preprocessing --features blas,linfa/netlib-static --lib — 55 passed, 0 failed, and cargo clippy clean under the same features (MKL still will not link locally on macOS, and OpenBLAS trips over an unrelated openblas-build/ureq TLS mismatch, so netlib-static is what I could actually run; that is a real LAPACK, so it exercises the same ndarray-linalg code path CI does)

🤖 Generated with Claude Code

https://claude.ai/code/session_01L4uUK6s1CWFwx6u7Rx28go

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 77.92%. Comparing base (7fe5c86) to head (70287fe).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #453      +/-   ##
==========================================
+ Coverage   77.53%   77.92%   +0.39%     
==========================================
  Files         106      104       -2     
  Lines        7585     7535      -50     
==========================================
- Hits         5881     5872       -9     
+ Misses       1704     1663      -41     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.