perf(preprocessing): whiten via the covariance matrix in Whitener::pca - #453
perf(preprocessing): whiten via the covariance matrix in Whitener::pca#453mysma-9403 wants to merge 3 commits into
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #453 +/- ##
==========================================
+ Coverage 77.93% 77.97% +0.03%
==========================================
Files 104 104
Lines 7547 7558 +11
==========================================
+ Hits 5882 5893 +11
Misses 1665 1665 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
70287fe to
7c6b880
Compare
|
First, I generally appreciate your interest in contributing to
but even the exact operation measured is unclear (was it a Performance is also generally tricky to measure, using tools like Also, for what it's worth, please remember that there's a person, not an agent, on the other side of this screen. It takes time and energy to review PRs like this. One of the key ways we evaluate the quality contributions is by getting a feel for the level of good-faith effort someone has put into their PR. When the content is clearly and primarily LLM-generated, the default assumption must be that it wasn't thoughtfully produced, which makes (at least me) much less likely to accept them. |
The existing `whitening` group varies the feature count at a fixed sample count and times `fit` together with `transform`. `transform` is linear in the sample count whichever way `fit` is implemented, so that group cannot show how a whitening method scales with the number of samples. Add a `whitening_fit` group that holds the feature count fixed and sweeps the sample count, timing `fit` on its own, for all three methods. A single `fit` there runs into the seconds, so the group lowers the shared 200-sample default rather than take tens of minutes.
`WhiteningMethod::Pca` took the SVD of the whole centered
(nsamples x nfeatures) data matrix, while the `Zca` and `Cholesky` arms
next to it first form the (nfeatures x nfeatures) covariance and decompose
that. The two are equivalent for PCA whitening: 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^2 / (nsamples - 1)`.
Forming the covariance first makes the decomposition's cost independent of
the sample count, leaving a single GEMM as the only work proportional to
`nsamples` - the same work `Zca` and `Cholesky` already did.
From the `whitening_fit` group added in the previous commit, timing `fit`
alone at nfeatures = 64, medians over 20 samples, p = 0.00 throughout:
nsamples before after change
2500 53.27 ms 2.06 ms -96.1%
5000 97.53 ms 2.97 ms -96.9%
10000 257.03 ms 4.60 ms -98.2%
20000 545.34 ms 8.06 ms -98.5%
40000 2.554 s 16.91 ms -99.3%
At 40000 samples `Cholesky` and `Zca` take 16.02 ms and 17.03 ms, so `Pca`
now sits alongside the two methods that already formed the covariance
instead of growing super-linearly away from them.
The epsilon floor is applied to the reconstructed singular values rather
than to the eigenvalues, so it keeps the meaning it had before, and the
eigenvalues are clamped at zero first because rounding can push a
numerically-zero one slightly negative and turn `sqrt` into a NaN.
Forming the covariance squares the condition number, so this is in
principle less accurate on badly conditioned inputs than an SVD of the data
matrix. `Zca` and `Cholesky` already make that trade; this brings `Pca` in
line with them rather than introducing a new compromise.
With at most as many samples as features the covariance is rank deficient,
and the SVD of the data matrix yields a differently shaped factor. Routing
that case to the original formulation keeps the shape of the whitening
matrix unchanged.
When `nsamples <= nfeatures` the whitening matrix keeps the shape it gets from the SVD of the data matrix, and that shape depends on the linalg backend: `linfa_linalg::svd` is compact and returns `nsamples x nfeatures`, `ndarray_linalg::svd` is full and returns `nfeatures x nfeatures`. That divergence is not introduced by the previous commit; it is existing behaviour that no test covered. Assert both shapes explicitly in `test_pca_matrix_more_features_than_samples` so it cannot drift unnoticed, and note in `test_pca_matrix_square_input` that the square case is the one where the two backends agree.
7c6b880 to
8c7211b
Compare
|
All fair, and the reproducibility point especially. What changed: Measurement. The crate already ships a criterion bench for this ( Timing Description. Rewritten - third person, only what bears on the diff. The dead session link and the co-author trailer are gone, from the description and from the commits. Commits. Three now, in the order that makes the progression followable: the bench, then the change, then the tests. Each one builds and passes tests on its own. One thing worth raising here rather than leaving buried in the diff: |
2 participants
WhiteningMethod::Pcatakes the SVD of the whole centered(nsamples, nfeatures)data matrix. TheZcaandCholeskyarms immediately next to it instead form the(nfeatures, nfeatures)covariance and decompose that.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^2 / (nsamples - 1). Forming the covariance first makes the decomposition's cost independent of the sample count, leaving a single GEMM as the only work proportional tonsamples.This only applies when
nsamples > nfeatures. Otherwise the covariance is rank deficient and the SVD of the data matrix yields a differently shaped factor, so that case keeps the original formulation.Benchmark
Both groups below are
criterion, live in the crate'swhitening_bench, and are part of this PR. Numbers come from criterion's own change detection.One of the two groups is added by this PR, so it does not exist on
master. To keep the harness identical on both sides of the comparison, the baseline is taken by checking out this branch and reverting only the implementation file:fitscaling with the sample countwhitening_fit, added here.nfeatures = 64,fittimed on its own, since that is the call this change touches and the sample count is the axis it is about.pcamasterpcathis branchcholeskyzcaMedians over 20 samples; every
pcarow is p = 0.00, with the 97% interval on the change spanning under half a percentage point.On
master,pcagrows super-linearly with the sample count whilecholeskyandzcagrow linearly. After the change all three are linear andpcasits alongside them, which is the point of the change: what remains proportional tonsamplesis the GEMM the other two methods already did.The
choleskyandzcacolumns are unchanged code, included as a reference line. Their apparent run-to-run deltas (up to a third at64x10000) are measurement noise at this sample count, and are a fair indication of the noise floor for the small measurements here — thepcaeffect is two orders of magnitude clear of it.Pre-existing group, varying the feature count
whitening, already in the crate.nsamples = 10000, timingfitfollowed bytransform.transformis untouched by this change and is linear in the sample count either way, so it dilutes the effect; included because it is the crate's existing measurement of this code.Medians over 200 samples, 97% confidence intervals within ±3.5% of each point estimate.
Both runs on the same machine, back to back: Intel Core i9-9980HK, 32 GB, macOS 15.7.7,
rustc1.93.1, criterion 0.5, default (linfa-linalg) backend.Details
max(s, 1e-8)keeps the meaning it had before.sqrt, because rounding can push a numerically-zero one slightly negative and produce aNaN.WᵀW.XᵀXsquares the condition number, so on badly conditioned inputs this is in principle less accurate than an SVD ofX.ZcaandCholeskyalready make that trade, so this makesPcaconsistent with them rather than introducing a new compromise.Shape when
nsamples <= nfeaturesThe fallback keeps
master's behaviour, which differs by backend:linfa_linalg::svdis a compact SVD and returns the(nsamples, nfeatures)factor,ndarray_linalg::svdis a full one and returns(nfeatures, nfeatures).test_pca_matrix_more_features_than_samples(20 × 50) asserts both;test_pca_matrix_square_input(16 × 16) covers the square case, where they agree.That backend split is not introduced here — it is
master's existing behaviour, pinned by the test so it cannot drift unnoticed.Checks
Run locally on the machine noted above, in addition to CI. The
blasruns usenetlib-staticrather than CI'sintel-mkl-static, since MKL is not available on this host:cargo test --release --workspace— no failurescargo test --release -p linfa-preprocessing --features blas,linfa/netlib-static --lib— 55 passed, 0 failedcargo clippy --workspace --all-targets --features benchmarks -- -D warnings— cleancargo clippy -p linfa-preprocessing --all-targets --features blas,linfa/netlib-static,linfa/benchmarks -- -D warnings— cleancargo fmt --all -- --check— cleanEach of the three commits was also built and tested on its own, so the branch is bisectable.