Three more ways to choose a PCA component count, and retire the fourth - #554
Merged
Conversation
A survey of what other packages implement turned up two schemes worth having and one worth removing. Added, all on PCA.select_n_components: cv_scheme="ek" is Eastment and Krzanowski (1982). An element is predicted by a score from a model without its column and a loading from a model without its row, so it enters neither decomposition. This is what Simca-P reports and what pcaMethods computes by default, so it is the scheme to reach for when a number has to line up with either. Its one leak is documented: preprocessing is fitted on the whole matrix, as in the classical scheme, so a cell reaches the fold models through its column's centre and spread, an effect of order 1/n. cv_scheme="sacv" and "gcv" are the leverage corrections of Josse and Husson (2012), which get leave-one-cell-out without refitting anything: inflate each residual by the leverage of the cell that produced it, the device that turns a regression residual into its leave-one-out counterpart. One decomposition instead of 560, and on a rank-3 block "sacv" picks the same count as the element-wise scheme and stays within five points of its curve through the optimum. Both degrade as the components approach the variables, because a column's leverage goes to one and the divisor goes to zero with it. The docstring says so and a test pins where it breaks. Deprecated: cv_scheme="row_wise" now raises a DeprecationWarning and goes in 2.0. It measures how well a held-out row reproduces itself, which falls monotonically to zero once the components equal the variables, so it cannot select anything. The package that shipped the equivalent removed it and said publicly it had been a mistake. 3229 tests pass, coverage 94.19%. ruff, ruff format and mypy clean. The strict docs build succeeds. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LXthGpHLQFfGubBiKFtGAE
Codecov Report✅ All modified and coverable lines are covered by tests. 📢 Thoughts on this report? Let us know! |
…othing The coverage bot found the one line of the new code no test reached, and it was the guard worth reaching: a rank-a bilinear model spends a(n + p - a - 1) free parameters, and once that meets the (n - 1)p degrees of freedom available, the inflation factor's denominator is no longer positive. On a 20 by 6 block the sixth component lands exactly there. The criterion now has to come back as NaN, rather than as a number carrying the wrong sign. The block builder took its component sizes from a hardcoded triple, so it could only make rank-3 data. It slices a fixed decreasing sequence instead, which leaves the default rank of three on the values the other tests are pinned to. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LXthGpHLQFfGubBiKFtGAE
Checking the new schemes against the curve Simca-P actually reported, on the LDPE data the book already uses, showed that two of them do not behave the way the synthetic tests suggested. On 54 rows and 19 variables the element-wise schemes turn over at two components, which is where Simca-P turns over. Neither leverage scheme turns over at all within the first five: their curves rise at every count, so the number they return is the largest one evaluated. That is the same failure that makes the row-wise scheme useless, and it was being returned silently as though it were an answer. The cause is visible in the data. The LDPE fit reaches 99.98% by eleven components, so the residual these schemes inflate has almost nothing left in it, and the degrees-of-freedom penalty cannot keep up. So any scheme whose Q2 rises at every component count evaluated, and whose recommendation therefore lands on the largest one tried, now raises a SpecificationWarning saying exactly that. It fires for the two leverage schemes on this data and for the row-wise scheme, and stays quiet for the two that hold values out. The user guide records the LDPE comparison, since a caveat with a number behind it is worth more than a general warning about approximations. One existing test counted the warnings row_wise emits and expected one. Its intent was that the flag row_wise ignores earns no scaling warning, which still holds, so it now asserts that rather than a count. 3231 tests pass. ruff, ruff format and mypy clean. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LXthGpHLQFfGubBiKFtGAE
Found on the FMC quality block, where the schemes were run over the full component range for the first time. Once the components reach the variables every column leverage is one, so the guard in "sacv" drops every cell as undefined. np.nansum of an all-NaN array is 0.0, not NaN, so PRESS came back as zero and Q2 as exactly 1.000, which then won the selection outright: 6 components on a 46 x 6 block and 8 on a 27 x 8 one, the worst answers available. The generalised criterion already guarded its own version of this; the smoothing one did not. Report nothing there instead, and where only some cells are dropped, rescale so the error and the reference it is divided by cover the same cells. The old code compared an error over part of the block against a null over all of it, which flatters the model in proportion to how many cells were dropped. No cell is dropped early on real data, so no published number moves; a high-leverage outlying row reaches it. Also reword the no-optimum warning, which named both element-wise schemes as the remedy whatever was running, so a caller already using one was pointed back at it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LXthGpHLQFfGubBiKFtGAE
The three schemes that factorise the matrix directly refuse a block with missing cells, which is right, but the refusal said only that some exist. What a caller does next depends on where they sit, and the message left them to find that out themselves. On the FMC dryer's quality block, 18 of the 19 missing cells are in one attribute: drop that column and all 46 batches survive, where dropping the incomplete rows would cost 19 of them. A block whose gaps are scattered has no such repair. The message now gives the count and says which case it is, so the choice is visible from the error alone. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LXthGpHLQFfGubBiKFtGAE
The no-optimum warning built its remedy behind "if alternatives:", but the list is the two element-wise schemes less whichever is running, so one always survives: a scheme cannot be both. Coverage flagged the branch that never runs, which is the right complaint about a guard that reads as if the empty case were possible. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LXthGpHLQFfGubBiKFtGAE
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.
A survey of how other packages cross-validate a PCA turned up two schemes worth having and one worth removing. All three changes are on
PCA.select_n_components.Added
cv_scheme="ek", the two-model scheme of Eastment and Krzanowski (1982). An element is predicted by a score from a decomposition that dropped its column and a loading from one that dropped its row, so it enters neither. This is what Simca-P reports as its PCAQ2and whatpcaMethods::Q2computes by default in R, so it is the scheme to reach for when a number has to line up with either. Rows and columns are each split intocvgroups rather than deleted one at a time, costing2 * n_foldsdecompositions instead ofn + p, and every fold model's singular vectors are sign-aligned to the full-data fit.One leak is documented rather than hidden. Centring and scaling are fitted once on the whole matrix, as in the classical scheme and in the reference implementations, so a cell reaches the fold models through its own column's centre and spread even though it is absent from both decompositions. That path is of order
1/n. It is why this is an option and not the default.cv_scheme="sacv"andcv_scheme="gcv", the leverage corrections of Josse and Husson (2012). These get leave-one-cell-out without refitting anything: inflate each residual by the leverage of the cell that produced it,(1 - 1/n - a_i)(1 - b_j), the same device that turns a regression residual into its leave-one-out counterpart viae_i / (1 - h_ii)."gcv"replaces the two leverages with one averaged constant and is blunter. Both are the defaults inFactoMineRandmissMDA.The point is cost. On a 100 by 12 block,
"sacv"costs one decomposition against 560 for the element-wise scheme, picks the same component count, and stays within five points of its curve through the optimum:ekfsacvThe last row is the trap, and it is pinned by a test. Both leverage schemes are first-order approximations that degrade as the components approach the variables, because a column's leverage is the share of its variance inside the retained subspace and reaches one when the components reach the variables. The divisor goes to zero with it. The docstring and the user guide say to read them well below that ceiling, which is why
FactoMineRdefaults to five.All three factorise the matrix directly, so they raise a clear
ValueErroron a block with missing cells and point at"ekf", which imputes an unmeasured cell alongside the held-out ones and never scores it. The refusal names how many cells are missing and whether one column holds them, because that decides the repair: on the FMC dryer's quality block 18 of the 19 sit in one attribute, so dropping that column keeps all 46 batches where dropping the incomplete rows would cost 19 of them. A block whose gaps are scattered has no such repair, and the message says which case it is. The three schemes also ignorescale_inside_folds,n_repeats,n_iterandtol, and they report no per-fold spread, since nothing is held out in folds that could disagree. A fabricated zero there would let the 1-SE rule read a false certainty, so the field isNaNinstead.Deprecated
cv_scheme="row_wise"now emits aDeprecationWarningand will be removed in 2.0, followingdocs/development/deprecation_policy.rst. It measures how well a held-out row reproduces itself: the row supplies the scores that rebuild it, so PRESS falls monotonically and reaches exactly zero once the components equal the variables. That is a compression error, not a prediction error, and it cannot select a component count.mdatoolsremoved the equivalent from its own PCA and said in its release notes that it "was a bad idea from the beginning". The existingSpecificationWarningstill fires alongside.The FMC quality block shows it plainly. On its complete 46 by 6 sub-block the row-wise scheme climbs -0.11, 0.93, 0.95, 0.97, 0.999, 1.000 and recommends six components, the number of attributes;
"ekf"on the same block recommends one.Fixed
The leverage criterion could report a perfect score on an empty total. Found by running every scheme over the whole component range on the FMC quality block, which is narrow enough (six to eight attributes) to reach the ceiling. Once the components reach the variables every column leverage is one, so the guard in
"sacv"drops every cell as undefined.np.nansumof an all-NaN array is0.0, notNaN, so PRESS came back as zero andQ2as exactly1.000, which then won the selection outright: six components on a 46 by 6 block and eight on a 27 by 8 one, the worst answers available."gcv"already guarded its own version of this;"sacv"did not. It now reportsNaNthere, and recommends 1 and 2 on those two blocks.Where only some cells are dropped, the criterion and the reference it is divided by now cover the same cells. The previous code compared an error over part of the block against a null over all of it, which flatters the model in proportion to how many cells were dropped. Nothing is dropped early on real data, so no published number moves; a high-leverage outlying row is what reaches this path.
The no-optimum warning named the scheme it was warning about. Its remedy listed both element-wise schemes whatever was running, so a caller already using
"ekf"was told to use"ekf". It now names only the alternatives, and says how many component counts it could actually evaluate rather than how many were asked for. The guard that built that remedy conditionally has gone with it: the list is the two schemes less whichever is running, so one always survives.Why the default did not change
The likelihood route that scikit-learn takes is elegant and immune to the circularity, but it assumes isotropic noise. On a rank-3 block with the noise ramped from 0.3 to 3.0 across the variables, which is ordinary for a plant with sensors of differing quality, the cross-validated log-likelihood and Minka's estimate both take ten or eleven components across five seeds while the element-wise scheme takes three every time.
"ekf"stays the default.Two datasets since bore that out from the other side. On LDPE,
"ekf"and"ek"both turn over at two components as Simca-P does (correlation with the recorded Simca-P curve +0.86 and +0.65 over the first eight counts), while"sacv"and"gcv"do not turn over at all there, because the fit reaches 99.98% and the residual they inflate has nothing left in it. On the FMC quality block, whose 19 missing cells are the ordinary case for plant data,"ekf"is the only one of the four that will run at all.Verification
uv run pytest: 3235 passed, 3 skipped, coverage 94.20% against the 92% gate. (An earlier revision of this description said 3234; that was the count before the last test landed.)ruff check .,ruff format --check .andmypy src/process_improveall clean on the current head.-W) succeeds.CHANGELOG.md,pyproject.tomlandCITATION.cffall at 1.84.0, dated 2026-09-10, a MINOR bump for new features. The fixes above are to code this PR itself adds, so they are folded into that entry rather than bumping again.🤖 Generated with Claude Code
https://claude.ai/code/session_01LXthGpHLQFfGubBiKFtGAE
Generated by Claude Code