# (c) Kevin Dunn, 2010-2026. MIT License. Based on own private work over the years.
"""Principal Component Analysis (PCA) estimator (ENG-01).
The sklearn-compatible :class:`PCA` transformer, with NIPALS, SVD and TSR
fitting paths and full missing-data support. Diagnostics, confidence limits and
plotting are pulled in from the sibling submodules and bound as convenience
methods after ``fit()``.
"""
from __future__ import annotations
import itertools
import logging
import numbers
import time
import typing
import warnings
import numpy as np
import pandas as pd
from scipy.linalg import blas, eigh
from sklearn.base import BaseEstimator, TransformerMixin, _fit_context
from sklearn.decomposition import PCA as _SkPCA # noqa: N811
from sklearn.model_selection import BaseCrossValidator
from sklearn.utils import Bunch
from sklearn.utils.validation import check_is_fitted, validate_data
from ..univariate.metrics import detect_outliers_esd
from ._base import _LatentVariableModel, _LazyFrame
from ._common import (
Q2_MIN_INCREMENT,
DataMatrix,
DiagnoseResult,
NotEnoughVarianceError,
SelectionRule,
SpecificationWarning,
_align_to_fit_features,
_fit_asked_for_by_the_caller,
_float_frame,
_nz,
_reject_sparse,
_rows_with_data,
_select_n_components,
check_at_least_one,
check_x_varies,
components_to_fit,
constant_columns,
epsqrt,
refuse_empty_block,
resolve_loop_settings,
unobserved_columns,
user_stacklevel,
)
from ._diagnostics import _t2_score_sd
from ._limits import JudgedRows, judged_rows_note
from ._nipals import quick_regress, scores_converged, ssq, terminate_check
from ._preprocessing import MCUVScaler, _warn_scaling_traps
from ._projection import coerce_observed_mask, operator_for_pattern, project_rows
logger = logging.getLogger(__name__)
class _EkfPress(typing.NamedTuple):
"""What one element-wise k-fold pass measured. See :func:`_pca_ekf_press`."""
press: np.ndarray
per_fold_press: np.ndarray
per_column_press: np.ndarray
null_model_ss: float
per_column_null_ss: np.ndarray
press_input_units: np.ndarray
def _low_rank_reconstruction(X: np.ndarray, rank: int) -> np.ndarray:
"""Return the best rank-``rank`` approximation of ``X``, the truncated SVD's ``U_r S_r V_r^T``.
The projector onto the leading singular subspace comes from the top
eigenvectors of the smaller Gram matrix, ``X X^T`` or ``X^T X``, rather than
a full SVD. Only ``rank`` vectors are computed: on a 394 x 650 fold of the
tablet spectra this is about five times faster than ``np.linalg.svd``
natively, and about eight times faster under Pyodide, where numpy has no
BLAS and ``scipy.linalg`` brings OpenBLAS.
Forming the Gram matrix squares the singular values. A component whose
singular value is below about ``1e-8`` of the largest therefore loses its
digits, but it also contributes less than that to the reconstruction, which
is below what the EM loop that calls this can resolve.
Parameters
----------
X : np.ndarray of shape (n_samples, n_features)
Centred (and usually scaled) data, all cells finite.
rank : int
Number of components to keep; capped at ``min(n_samples, n_features)``.
Returns
-------
np.ndarray of shape (n_samples, n_features)
The reconstruction.
"""
n, p = X.shape
rank = min(rank, n, p)
wide = n <= p # pragma: no mutate - on a square X both Gram matrices give the same projector
# syrk reads X.T without a copy (it is Fortran-ordered) and fills the upper
# triangle only: X X^T when trans=1, X^T X when trans=0.
gram = blas.dsyrk(
1.0, # pragma: no mutate - any multiple of the Gram matrix has the same eigenvectors
X.T,
trans=int(wide),
)
m = len(gram)
_, basis = eigh(
gram,
lower=False, # pragma: no mutate - eigh reads lower=None as False
subset_by_index=(m - rank, m - 1),
)
if wide:
return basis @ (basis.T @ X)
return (X @ basis) @ basis.T
def _opposes_convention(loadings: np.ndarray) -> np.ndarray:
"""Return, per column, whether its largest-magnitude entry is negative.
Every fitting path flips such a component, loading and score together, so the
largest-magnitude element of each loading is positive (Wold, Esbensen and
Geladi, PCA, CILS, 1987, p 42) and the fitted signs do not depend on LAPACK or
on the NIPALS seed. A tie goes to the first entry.
Parameters
----------
loadings : np.ndarray of shape (n_features, n_components)
Unit-length loadings, one per column.
Returns
-------
np.ndarray of shape (n_components,)
True for each column to flip.
"""
# The largest entry of a unit vector is never 0, so its sign bit says whether it is negative.
return np.signbit(loadings[np.argmax(np.abs(loadings), axis=0), np.arange(loadings.shape[1])])
def _pca_ekf_press( # noqa: PLR0913, PLR0915
X: np.ndarray,
max_components: int,
*,
n_folds: int,
n_repeats: int,
ekf_max_iter: int,
ekf_tol: float,
scale_inside_folds: bool,
random_state: int | None,
) -> _EkfPress:
"""Element-wise k-fold (ekf) PCA cross-validation.
Partitions the elements of ``X`` into ``n_folds`` element-folds (each cell
is held out exactly once across folds); for each fold and each candidate
component count, the held-out cells are masked, initialised from the
in-fold column means, and refined by iterating SVD reconstruction
(Expectation-Maximisation style). The fitted model never sees the
held-out true values, so the squared error of the prediction is an honest
out-of-sample PRESS - the independence requirement of Bro, Kjeldahl,
Smilde & Kiers (2008, *Anal. Bioanal. Chem.* 390:1241-1251) that the
row-wise CV scheme violates.
Parameters
----------
X : np.ndarray of shape (n_samples, n_features)
Data matrix, of floats. With ``scale_inside_folds=True`` the raw
unscaled matrix may be passed; with ``False`` the caller is expected
to have mean-centred (and usually unit-variance scaled) it.
max_components : int
Maximum number of components to evaluate; PRESS is computed for
``1 .. max_components``.
n_folds : int
Number of element-folds. Bro 2008 uses 7 as a typical default; 5 is
a faster choice that still gives a stable curve.
n_repeats : int
Number of times to repeat the ekf pass with a fresh random fold
permutation. Each repeat covers every cell exactly once; ``n_repeats
> 1`` averages over different element-fold partitions, narrowing
the per-component PRESS standard error at extra runtime.
ekf_max_iter : int
Maximum number of EM iterations per fold and component count. Named
for the loop it bounds: this is the imputation EM inside one fold, a
different loop from the per-component NIPALS loop that ``PCA.max_iter``
caps (#588).
ekf_tol : float
Relative change in the held-out cell predictions below which EM
stops early. Likewise distinct from ``PCA.tol``.
scale_inside_folds : bool
If True, fit per-column mean and unit-variance constants on each
fold's in-fold cells and apply them to the whole matrix before
running EM. This removes the centring/scaling leakage of the
previous default (which used a single set of constants iteratively
recomputed from the imputed matrix). PRESS is then measured in that
same space, so a variable's contribution reflects its correlation
structure rather than its units; ``press_input_units`` carries the
other scale for callers who need it. If False, the scheme reverts to
the prior behaviour: the caller is responsible for scaling, the
in-loop column-mean drifts with the EM imputation, and PRESS is in
the units of whatever matrix was passed.
random_state : int or None
Seed for the element-fold permutation, for reproducibility across
repeats.
Returns
-------
_EkfPress
A named tuple with five fields. Under ``scale_inside_folds=True`` the
first three live in the space each fold was fitted in, so every column
contributes to the total in proportion to its correlation structure
rather than its units; under ``False`` there is no in-fold scale and
they are in the input units, unchanged from earlier releases.
press : np.ndarray of shape (max_components,)
PRESS per component count, averaged over ``n_repeats`` passes so
the scale is comparable to a single-pass run.
per_fold_press : np.ndarray of shape (max_components, n_folds * n_repeats)
Per-fold PRESS contributions across every fold of every repeat;
drives the 1-SE rule's standard error.
per_column_press : np.ndarray of shape (max_components, n_features)
The same total, split by variable, for the per-variable Q2.
null_model_ss : float
What the null model ("predict every held-out cell by the mean of
the cells that were not held out") got wrong, measured the same
way ``press`` is: the reference ``press`` is compared against.
Computed here rather than by the caller because it needs each
fold's own centring and scaling constants.
per_column_null_ss : np.ndarray of shape (n_features,)
That reference, split by variable.
press_input_units : np.ndarray of shape (max_components,)
PRESS in the units of the matrix that was passed in, for callers
comparing prediction error against instrument error. Always in
input units, whatever ``scale_inside_folds`` is.
References
----------
Bro, R., Kjeldahl, K., Smilde, A. K., & Kiers, H. A. L. (2008).
Cross-validation of component models: a critical look at current
methods. *Anal. Bioanal. Chem.*, 390(5), 1241-1251. PMID 18214448.
Camacho, J., & Ferrer, A. (2012). Cross-validation in PCA models with
the element-wise k-fold (ekf) algorithm: theoretical aspects.
*J. Chemometrics*, 26(7), 361-373. DOI 10.1002/cem.2440.
"""
n, p = X.shape
rng = np.random.default_rng(random_state)
# A cell the caller never measured has no true value to predict, so it is
# never held out and never scored. It still has to be filled for the SVD,
# which cannot see a NaN, so EM imputes it alongside the held-out cells.
observed = ~np.isnan(X)
total_folds = n_folds * n_repeats
# NaN, not zero: a fold that holds out no cells has no PRESS, and leaving
# a structural zero there let the caller's standard error treat it as a
# real fold that happened to predict perfectly. On a matrix with fewer
# cells than folds every cell lands in the last fold, so the empty folds
# fabricated a standard error out of zeros. np.nansum below keeps the
# PRESS total itself unchanged, since an empty fold contributes nothing.
per_fold_press = np.full((max_components, total_folds), np.nan)
# PRESS is what the model got wrong; the null reference is what predicting
# the in-fold column mean would have got wrong on the same cells. Both are
# accumulated here, in the same space, because only this loop knows each
# fold's centring and scaling constants. The reference does not depend on
# the component count, hence one value per fold rather than a grid.
per_fold_null = np.zeros(total_folds)
per_column_press = np.zeros((max_components, p))
per_column_null = np.zeros(p)
press_input_units = np.zeros(max_components)
fold_counter = 0
observed_positions = np.flatnonzero(observed)
n_cells = observed_positions.size
fold_size = n_cells // n_folds
for _ in range(n_repeats):
# Deal the observed cells, in a random order, into folds of fold_size
# cells each, the last fold taking the rest (all of them when there are
# fewer cells than folds): every fold has ~n_cells/n_folds cells, not all
# the same column, not the same row. Unmeasured cells are in no fold.
dealt = observed_positions[rng.permutation(n_cells)]
sizes = np.full(n_folds, fold_size)
sizes[-1] = n_cells - fold_size * (n_folds - 1)
fold_of_dealt = np.repeat(np.arange(n_folds), sizes)
for k in range(n_folds):
# (n, p) bool, True where the cell is held out
mask = np.isin(np.arange(n * p), dealt[fold_of_dealt == k]).reshape(n, p)
if not mask.any():
fold_counter += 1
continue
# Fit per-column centring/scaling on the in-fold cells. The centre
# is needed either way, as the null model's prediction for every
# held-out cell in this fold; only the True path also centres and
# scales the matrix it fits, and recomputes nothing inside EM.
# Held out this fold, or never measured: both are unknown to the
# fit and both are filled by EM below. Only ``mask`` is scored.
impute = mask | ~observed
# MCUVScaler's, in any units: a column whose in-fold values are all
# equal (or one, or none) keeps a scale of 1, and its centre is the
# value (0 for none).
in_fold = MCUVScaler().fit(np.where(impute, np.nan, X)) # pragma: no mutate - a None reads as NaN too
col_centre = in_fold.center_.to_numpy()
col_scale = in_fold.scale_.to_numpy() if scale_inside_folds else np.ones(p)
# Initial imputation in original space: held-out cells take the
# in-fold column mean (which is ``col_centre`` under
# ``scale_inside_folds=True``, and the same value under False).
Xtr = X.copy()
for j in range(p):
m_j = impute[:, j]
if m_j.any():
Xtr[m_j, j] = col_centre[j]
# The cells' own scale and centre for this fold, one value per
# held-out cell. Under scale_inside_folds=False the scale is 1, so
# every measurement below stays in the input units, as it was.
cell_scale = np.broadcast_to(col_scale, X.shape)[mask]
cell_centre = np.broadcast_to(col_centre, X.shape)[mask]
column_of_cell = np.broadcast_to(np.arange(p), X.shape)[mask]
null_residual = (X[mask] - cell_centre) / cell_scale
per_fold_null[fold_counter] = float(np.sum(null_residual**2))
per_column_null += np.bincount(column_of_cell, weights=null_residual**2, minlength=p)
held_scale = np.broadcast_to(col_scale, X.shape)[impute]
for a in range(1, max_components + 1):
# Every component count starts EM again from the column means.
# Past the data's rank the EM has no unique fixed point, so the
# PRESS it reaches depends on the starting point and on
# ekf_max_iter. Starting from the cells imputed for A - 1 is
# faster but changes the curve, and with it the selection.
Xa = Xtr.copy()
# Convergence is judged on the held-out cells in the space the
# fold is fitted in, not in the input units. Otherwise a column
# re-expressed in different units dominates the norm and changes
# how many EM iterations are taken (#546). Under
# scale_inside_folds=False the divisor is 1 throughout.
prev_held = Xa[impute] / held_scale
for _iteration in range(ekf_max_iter):
if scale_inside_folds:
# Centre and scale by the FIXED in-fold constants; the
# in-fold cells now sit on the analysis scale and the
# held-out cells move under EM.
Xs = (Xa - col_centre) / col_scale
recon = _low_rank_reconstruction(Xs, a) * col_scale + col_centre
else:
# Prior behaviour: recompute the column mean each
# iteration. The mean drifts slightly as the held-out
# cells are updated.
col_mean = Xa.mean(axis=0)
recon = _low_rank_reconstruction(Xa - col_mean, a) + col_mean
Xa[impute] = recon[impute]
held = Xa[impute] / held_scale
delta = np.linalg.norm(held - prev_held)
# Relative to the held-out cells' own size, in any units. A floor of 1 stopped EM
# after one pass on data in small units, which changed the PRESS curve (#1070).
threshold = ekf_tol * float(np.linalg.norm(prev_held))
if delta < threshold: # pragma: no mutate - < and <= differ only at a delta of exactly threshold
break
prev_held = held
error_input_units = X[mask] - Xa[mask]
error = error_input_units / cell_scale
per_fold_press[a - 1, fold_counter] = float(np.sum(error**2))
per_column_press[a - 1] += np.bincount(column_of_cell, weights=error**2, minlength=p)
press_input_units[a - 1] += float(np.sum(error_input_units**2))
fold_counter += 1
# Average over repeats so PRESS stays on the per-cell scale. nansum so the
# empty-fold NaNs above do not propagate into the total.
return _EkfPress(
press=np.nansum(per_fold_press, axis=1) / n_repeats,
per_fold_press=per_fold_press,
per_column_press=per_column_press / n_repeats,
null_model_ss=float(np.sum(per_fold_null)) / n_repeats,
per_column_null_ss=per_column_null / n_repeats,
press_input_units=press_input_units / n_repeats,
)
def _ekf_null_reference(
ekf: _EkfPress,
X: np.ndarray,
*,
scale_inside_folds: bool,
) -> tuple[float, np.ndarray]:
"""Return the sum of squares :math:`Q^2` measures the ekf PRESS against.
The reference has to be measured the same way PRESS was, or the ratio of the
two is not a fraction of anything. Under in-fold scaling it therefore comes
from the fold loop, which is the only place that knows each fold's constants;
a reference computed out here could only be in the input units, and then a
single high-variance column would set the whole curve however the folds were
actually fitted (#546). Under ``scale_inside_folds=False`` there is no
in-fold scale to match, so the centred total sum of squares of the matrix the
caller passed stays the reference, unchanged from earlier releases.
Parameters
----------
ekf : _EkfPress
What the fold loop measured.
X : np.ndarray of shape (n_samples, n_features)
The matrix as it was passed to the cross-validation.
scale_inside_folds : bool
Whether each fold standardised the matrix before fitting it.
Returns
-------
tuple[float, np.ndarray]
The pooled reference, and the same split by variable (shape
``(n_features,)``).
"""
if scale_inside_folds:
return ekf.null_model_ss, ekf.per_column_null_ss
centred_ss = np.nansum((X - np.nanmean(X, axis=0)) ** 2, axis=0)
return float(centred_ss.sum()), centred_ss
def _preprocess_for_cell_schemes(X: np.ndarray, scheme: str) -> np.ndarray:
"""Mean-centre and unit-variance scale the whole matrix, once.
The schemes below are the classical ones, and all of them are defined on a
single preprocessed matrix rather than on per-fold constants: there are no
folds to fit constants on in ``"sacv"`` and ``"gcv"``, and the two SVD
families of ``"ek"`` are both taken from the same preprocessed matrix. That
is what ``FactoMineR`` and ``pcaMethods`` do, so a number computed here is
comparable with theirs.
Parameters
----------
X : np.ndarray of shape (n_samples, n_features)
Raw data matrix. Missing cells are not supported by these schemes.
scheme : str
Name of the calling scheme, used only in the error message.
Returns
-------
np.ndarray
The centred, unit-variance scaled matrix.
Raises
------
ValueError
If ``X`` holds missing cells. These schemes factorise ``X`` directly,
and an SVD cannot see a NaN.
"""
per_column = np.isnan(X).sum(axis=0)
n_missing = int(per_column.sum())
if n_missing:
# Where the misses sit decides what to do about them, so say. A block whose
# gaps are one column's is repaired by dropping that column, keeping every
# row; one whose gaps are scattered is not.
worst = int(np.argmax(per_column))
share = int(per_column[worst])
where = (
f"{share} of them in column {worst}"
if share > n_missing / 2
else f"spread over {int((per_column > 0).sum())} columns"
)
raise ValueError(
f"cv_scheme={scheme!r} cannot take a block with missing cells, because it "
f"factorises the matrix directly. This block has {n_missing} of {X.size}, "
f"{where}. Use cv_scheme='ekf', which imputes an unmeasured cell alongside "
"the held-out ones and never scores it, or drop what is missing first."
)
return MCUVScaler().fit(X).transform(X).to_numpy() # in any units, as for every other scheme
def _pca_ckf_press(
Z: np.ndarray,
max_components: int,
*,
n_folds: int,
random_state: int | None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, np.ndarray]:
"""Column-wise k-fold (ckf) PCA cross-validation (#374).
The columns are split into ``n_folds`` groups. For each group ``g`` and each candidate
component count ``a``, the loadings ``P`` come from an SVD of the whole block, but the
scores used to predict the held-out columns are computed from the *retained* columns
only::
T_g = Z[:, -g] P[-g, :] (P[-g, :]' P[-g, :])^-1
Zhat[:, g] = T_g P[g, :]'
so no held-out value ever appears in the score that predicts it. The reference is the
column mean, which is zero in this centred space, so ``Q2 = 1 - PRESS / SS(Z)``.
**ckf is not as clean as ekf, and the difference is measurable.** The loadings are
still fitted on a matrix that includes the held-out columns, so information leaks
through ``P`` even though it does not leak through ``T``. Camacho and Ferrer (2012)
make exactly this point, and it is why ``cv_scheme="ekf"`` remains the default here:
ekf holds out individual cells and imputes them, so nothing held out reaches the model
at all.
On a 60-by-12 block of pure noise, ekf returns a negative Q2 at every component count,
as it must when there is nothing to predict; ckf on the same block returns **+0.017**
at one component. That is the leakage, and a test pins it.
ckf is offered anyway, for two reasons: a great deal of published chemometrics uses
it, so a number that has to line up with a paper may need it, and it costs one
decomposition rather than ``n_folds * n_repeats * max_components``, which matters on a
block with far more columns than samples.
Parameters
----------
Z : np.ndarray of shape (n_samples, n_features)
Centred and scaled data, as produced by :func:`_preprocess_for_cell_schemes`.
max_components : int
PRESS is computed for ``1 .. max_components``.
n_folds : int
Number of column groups. Capped at the number of columns, since a group cannot be
empty; ``n_folds == n_features`` is leave-one-variable-out.
random_state : int or None
Seed for the column permutation that forms the groups.
Returns
-------
press : np.ndarray of shape (max_components,)
Total held-out sum of squared errors per component count. NaN for a count that a
column group cannot evaluate, because it leaves fewer columns in than there are
components: a total over the other groups predicts fewer columns, and would read as
a better fit than any count every group evaluates (#675).
per_fold_press : np.ndarray of shape (max_components, n_folds)
The same total, split by column group, so the 1-SE rule has a spread to work with;
NaN for a group that cannot evaluate the count.
per_column_press : np.ndarray of shape (max_components, n_features)
The same total, split by variable; NaN for the variables of such a group.
null_model_ss : float
``SS(Z)``: what predicting every held-out value by its column mean gets wrong.
per_column_null_ss : np.ndarray of shape (n_features,)
That reference, split by variable.
References
----------
Camacho, J., & Ferrer, A. (2012). Cross-validation in PCA models with the
element-wise k-fold (ekf) algorithm: theoretical aspects. *J. Chemometrics*,
26(7), 361-373. DOI 10.1002/cem.2440.
"""
n_features = Z.shape[1]
# No more folds than columns, so np.array_split never hands back an empty fold.
n_folds = max(1, min(int(n_folds), n_features))
rng = np.random.default_rng(random_state)
groups = np.array_split(rng.permutation(n_features), n_folds)
# One SVD serves every component count: the rank-a loadings are the first a columns of
# V, so the whole curve costs one decomposition rather than max_components of them.
# The only caller, select_n_components, caps max_components at min(N - 1, K), within
# the min(N, K) components this SVD returns.
_, _, vt = np.linalg.svd(Z, full_matrices=False) # pragma: no mutate - rows past min(N, K) are never read
press = np.full(max_components, np.nan)
per_fold_press = np.full((max_components, n_folds), np.nan)
per_column_press = np.full((max_components, n_features), np.nan)
for a in range(1, max_components + 1):
loadings = vt[:a, :].T # (n_features, a)
for fold, held_out in enumerate(groups):
keep = np.setdiff1d(np.arange(n_features), held_out)
if keep.size < a:
# Fewer retained columns than components: the least-squares score is not
# determined. Leave this fold's contribution as NaN rather than inventing
# one, which is what `per_fold_press` being NaN-initialised is for.
continue
scores = np.linalg.lstsq(loadings[keep, :], Z[:, keep].T)[0].T
residual = Z[:, held_out] - scores @ loadings[held_out, :].T
squared = residual**2
per_column_press[a - 1, held_out] = squared.sum(axis=0)
per_fold_press[a - 1, fold] = float(squared.sum())
if not np.isnan(per_fold_press[a - 1, :]).any(): # every column predicted, so the total is comparable
press[a - 1] = float(per_fold_press[a - 1, :].sum())
per_column_null_ss = (Z**2).sum(axis=0)
return press, per_fold_press, per_column_press, float(per_column_null_ss.sum()), per_column_null_ss
def _leverage_corrected_press(
X: np.ndarray,
max_components: int,
*,
method: typing.Literal["sacv", "gcv"],
) -> tuple[np.ndarray, np.ndarray, float, np.ndarray]:
r"""Leave-one-cell-out prediction error, approximated without refitting anything.
Holding out one cell at a time and refitting costs :math:`n \times p` fits per
component count. Josse and Husson (2012) avoid all of them: fit once, then
inflate the ordinary residual by the leverage of the cell that produced it,
the same device that turns a regression residual into its
leave-one-out counterpart via :math:`e_i / (1 - h_{ii})`.
With :math:`\mathbf{U}`, :math:`\mathbf{V}` the first ``a`` left and right
singular vectors, the row and column leverages are
:math:`a_i = \sum_k u_{ik}^2` and :math:`b_j = \sum_k v_{jk}^2`, and
- ``"sacv"`` divides each residual by :math:`(1 - 1/n - a_i)(1 - b_j)`,
cell by cell. This is the *smoothing approximation* of the
cross-validation criterion.
- ``"gcv"`` replaces both leverages by one averaged constant,
:math:`np / [(n - 1)p - a(n + p - a - 1)]`, where the subtracted term
counts the free parameters of a centred rank-``a`` bilinear model. This
is *generalised* cross-validation, and it is the cheaper, blunter cousin.
Neither holds any data out. They correct the circularity analytically
instead, which is why they cost one SVD in total rather than one per fold.
Both are first-order approximations and both degrade as ``a`` approaches
the number of variables. A column's leverage is the share of its variance
inside the retained subspace, so it reaches one when the components reach
the variables, and the divisor reaches zero with it. Read these criteria
over a component count well below the number of variables; ``FactoMineR``
defaults to five for the same reason.
Parameters
----------
X : np.ndarray of shape (n_samples, n_features)
Data matrix, already centred and scaled.
max_components : int
Largest component count to evaluate; the criterion is returned for
``1 .. max_components``.
method : {"sacv", "gcv"}
Which of the two corrections to apply.
Returns
-------
press : np.ndarray of shape (max_components,)
The corrected sum of squared residuals per component count.
per_column_press : np.ndarray of shape (max_components, n_features)
The same total, split by variable.
null_model_ss : float
What predicting every cell by its column mean gets wrong. Since ``X``
arrives centred, that is its total sum of squares.
per_column_null_ss : np.ndarray of shape (n_features,)
That reference, split by variable.
References
----------
Josse, J., & Husson, F. (2012). Selecting the number of components in
principal component analysis using cross-validation approximations.
*Computational Statistics & Data Analysis*, 56(6), 1869-1879.
"""
n, p = X.shape
U, S, Vt = np.linalg.svd(X, full_matrices=False) # pragma: no mutate - vectors past min(N, K) are never read
press = np.full(max_components, np.nan)
per_column_press = np.full((max_components, p), np.nan)
null_ss = np.sum(X**2, axis=0)
# The only caller caps max_components at min(N - 1, K), so every count has its
# singular vectors.
for a in range(1, max_components + 1):
residual = X - (U[:, :a] * S[:a]) @ Vt[:a]
if method == "sacv":
row_leverage = np.sum(U[:, :a] ** 2, axis=1)
column_leverage = np.sum(Vt[:a] ** 2, axis=0)
denominator = np.outer(1.0 - 1.0 / n - row_leverage, 1.0 - column_leverage)
# A cell whose leverage reaches one is reconstructed entirely by
# itself, so its leave-one-out residual is not defined. Drop it
# rather than let a division by ~0 dominate the total.
scored = denominator > epsqrt # pragma: no mutate - > and >= differ only at a denominator of exactly epsqrt
if not scored.any():
# Nothing is left to measure. Summing an empty total would
# report zero error, which reads as a perfect fit and would win
# the selection outright. Report nothing instead, as the
# generalised criterion does when its own guard fires. Both
# leverages only grow with a, so no later count scores a cell
# either, and stopping here would report the same.
continue # pragma: no mutate - see above: break gives the same result
per_column = np.sum((residual / np.where(scored, denominator, np.inf)) ** 2, axis=0)
# Both totals must cover the same cells. Where some were dropped, a
# column's error is rescaled onto its whole sum of squares, so that
# dividing it by the block's own gives the criterion over the cells
# that were actually scored.
scored_null = np.sum(np.where(scored, X**2, 0.0), axis=0)
rescaled = scored_null > epsqrt # pragma: no mutate - > and >= differ only at exactly epsqrt
per_column[rescaled] *= null_ss[rescaled] / scored_null[rescaled]
else:
remaining = (n - 1) * p - a * (n + p - a - 1)
if remaining <= 0:
# More free parameters than degrees of freedom: the criterion
# has nothing left to measure against. For a <= min(N - 1, K)
# the count of free parameters only grows with a, so no later
# count has any left either.
continue # pragma: no mutate - see above: break gives the same result
per_column = np.sum((residual * (n * p / remaining)) ** 2, axis=0)
per_column_press[a - 1] = per_column
press[a - 1] = float(np.sum(per_column))
return press, per_column_press, float(np.sum(X**2)), null_ss
def _eastment_krzanowski_press(
X: np.ndarray,
max_components: int,
*,
n_folds: int,
random_state: int | None,
) -> tuple[np.ndarray, np.ndarray, float, np.ndarray]:
r"""Cross-validated prediction error by the two-model scheme of Eastment and Krzanowski.
An element is predicted by a score taken from a model that never saw its
**column** and a loading taken from a model that never saw its **row**, so
:math:`x_{ij}` enters neither decomposition:
.. math::
\hat{x}_{ij} = \sum_{k=1}^{a} u^{(-j)}_{ik}\sqrt{d^{(-j)}_k}\,
\sqrt{d^{(-i)}_k}\, v^{(-i)}_{jk}
This is what Simca-P reports as its PCA :math:`Q^2`, and what
``pcaMethods::Q2(type = "krzanowski")`` computes in R, so it is the scheme
to use when a number has to line up with either of those.
Rows and columns are each split into ``n_folds`` groups rather than deleted
one at a time, which costs ``2 * n_folds`` decompositions instead of
``n + p``. The sign of every fold model's singular vectors is aligned to the
full-data fit, because an SVD fixes each vector only up to sign and two fold
models must agree before their product means anything.
One caveat, which the classical scheme and the reference implementations
share: centring and scaling are fitted once on the whole matrix, 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
:math:`1/n` and vanishes as rows are added, but it is not nothing, and it is
why this scheme is not offered as the default. ``"ekf"`` fits its constants
inside each fold and has no such path.
Parameters
----------
X : np.ndarray of shape (n_samples, n_features)
Data matrix, already centred and scaled, with no missing cells.
max_components : int
Largest component count to evaluate.
n_folds : int
Number of row groups, and of column groups.
random_state : int or None
Seed for the row and column permutations.
Returns
-------
press : np.ndarray of shape (max_components,)
Prediction error sum of squares per component count.
per_column_press : np.ndarray of shape (max_components, n_features)
The same total, split by variable.
null_model_ss : float
What predicting every cell by its column mean gets wrong.
per_column_null_ss : np.ndarray of shape (n_features,)
That reference, split by variable.
References
----------
Eastment, H. T., & Krzanowski, W. J. (1982). Cross-validatory choice of the
number of components from a principal component analysis. *Technometrics*,
24(1), 73-77.
"""
n, p = X.shape
rng = np.random.default_rng(random_state)
row_group = np.array_split(rng.permutation(n), n_folds)
column_group = np.array_split(rng.permutation(p), n_folds)
U_full, _, Vt_full = np.linalg.svd(X, full_matrices=False) # pragma: no mutate - full_matrices=None reads as False
def aligned(vectors: np.ndarray, reference: np.ndarray) -> np.ndarray:
"""Flip each column whose direction opposes the full-data fit.
A fold model is fitted on fewer rows or columns than the whole, so it has at
most as many vectors as the full-data fit: each is matched to the one of the
same rank.
"""
opposed = np.sum(vectors * reference[:, : vectors.shape[1]], axis=0) < 0
return np.where(opposed, -vectors, vectors)
# Loadings from models that never saw the rows of their group.
row_models: list[tuple[np.ndarray, np.ndarray, np.ndarray]] = []
for rows in row_group:
kept = np.setdiff1d(np.arange(n), rows)
_, S_r, Vt_r = np.linalg.svd(X[kept], full_matrices=False) # pragma: no mutate - None reads as False
row_models.append((rows, aligned(Vt_r.T, Vt_full.T), S_r))
# Scores from models that never saw the columns of their group.
column_models: list[tuple[np.ndarray, np.ndarray, np.ndarray]] = []
for columns in column_group:
kept = np.setdiff1d(np.arange(p), columns)
U_c, S_c, _ = np.linalg.svd(X[:, kept], full_matrices=False) # pragma: no mutate - None reads as False
column_models.append((columns, aligned(U_c, U_full), S_c))
# A model has as many vectors as singular values.
available = min(S.size for _, _, S in [*row_models, *column_models])
if available < max_components:
warnings.warn(
f"cv_scheme='ek' can evaluate at most {available} component(s) with "
f"{n_folds} folds on a {n} by {p} block, not {max_components}: each fold "
"model is fitted on fewer rows or columns than the whole. Component counts "
"above that are reported as NaN. Use fewer folds to evaluate more.",
SpecificationWarning,
stacklevel=user_stacklevel(),
)
press = np.full(max_components, np.nan)
per_column_press = np.full((max_components, p), np.nan)
for a in range(1, min(available, max_components) + 1):
prediction = np.zeros_like(X)
for rows, V_r, S_r in row_models:
for columns, U_c, S_c in column_models:
weight = np.sqrt(S_c[:a]) * np.sqrt(S_r[:a])
prediction[np.ix_(rows, columns)] = (U_c[rows, :a] * weight) @ V_r[columns, :a].T
squared = (X - prediction) ** 2
press[a - 1] = float(np.sum(squared))
per_column_press[a - 1] = np.sum(squared, axis=0)
return press, per_column_press, float(np.sum(X**2)), np.sum(X**2, axis=0)
[docs]
class PCA(_LatentVariableModel, TransformerMixin, BaseEstimator):
"""Principal Component Analysis with support for missing data.
Parameters
----------
n_components : int or None
Number of principal components to extract. ``None`` asks for as many
components as the data supports: it is resolved at fit time to
``min(n_samples, n_features)``, which is also the ceiling an explicit
request is clamped to (with a ``SpecificationWarning``). The resolved
count is available on the fitted attribute ``n_components_``.
algorithm : str, default="auto"
Algorithm to use for fitting the model.
- ``"auto"``: Uses SVD when data is complete, NIPALS when data has missing values.
- ``"svd"``: Singular Value Decomposition. Requires complete data.
- ``"nipals"``: Non-linear Iterative Partial Least Squares. Handles missing data.
- ``"tsr"``: Trimmed Score Regression. Handles missing data.
tol : float, default=``epsqrt`` (about 1.5e-8)
Relative convergence tolerance for the iterative algorithms: the loop
stops once the norm of the difference between two successive score
vectors, relative to the norm of the score vector, falls below this.
See :func:`terminate_check`. Ignored by ``algorithm="svd"``, which is
direct rather than iterative, but checked by every ``fit``: it must lie
between machine precision and 10.
max_iter : int, default=500
Maximum number of iterations per component for the iterative
algorithms. A component that reaches the cap without converging emits
a :class:`SpecificationWarning`. Ignored by ``algorithm="svd"``, but
checked by every ``fit``: it must be at least 1. The same default as
every other iterative estimator in the package (#588).
missing_data_settings : dict or None, default=None
PCA reads no keys from it; it is kept for the same signature as the
estimators that name their imputation method here. The keys
``md_tol`` and ``md_max_iter`` were removed in 2.0 and raise a
:class:`ValueError`; set ``tol`` and ``max_iter`` instead.
Attributes (after fitting)
--------------------------
n_components_ : int
The resolved number of components actually fitted (the constructor
parameter clamped to ``min(n_samples, n_features)``; the parameter
itself is left as the user set it, including ``None``).
scores_ : pd.DataFrame of shape (n_samples, n_components)
The score matrix (T).
loadings_ : pd.DataFrame of shape (n_features, n_components)
The loading matrix (P).
r2_per_component_ : pd.Series of length n_components
Fractional R² (on X) explained by each component.
r2_cumulative_ : pd.Series of length n_components
Cumulative R² after each component.
r2_per_variable_ : pd.DataFrame of shape (n_features, n_components)
Per-variable cumulative R² after each component.
spe_ : pd.DataFrame of shape (n_samples, n_components)
Per-row SPE diagnostic; stored as the square root of the row
sum-of-squared X-residuals (so it is on the residual scale, not the
squared scale). **One column per component, not one value per row**:
column ``a`` is the SPE of the model truncated at ``a`` components, and
the last column is the value at the full fitted model. Reach for a
single number per observation with ``model.spe_.iloc[:, -1]``, not with
``np.asarray(model.spe_).ravel()``: ravel happens to give the right
answer at one component and silently gives ``n_samples * n_components``
values above it. A row with nothing observed has no residual, and an
SPE of NaN, as in ``diagnose``.
hotellings_t2_ : pd.DataFrame of shape (n_samples, n_components)
Cumulative Hotelling's T² statistic. Per-component, exactly as ``spe_``
above: column ``a`` uses the first ``a`` components and the last column
is the value at the full fitted model.
explained_variance_ : np.ndarray of shape (n_components,)
Variance explained by each component.
scaling_factor_for_scores_ : pd.Series of length n_components
Standard deviation per score (sqrt of explained variance).
has_missing_data_ : bool
Whether the training data contained missing values.
fitting_info_ : dict
Timing and iteration info from the fitting algorithm.
"""
_valid_algorithms: typing.ClassVar[list[str]] = ["auto", "svd", "nipals", "tsr"]
_parameter_constraints: typing.ClassVar = {
"n_components": [numbers.Integral, None],
"algorithm": [str],
"tol": [numbers.Real],
"max_iter": [numbers.Integral],
"missing_data_settings": [dict, None],
}
def __init__(
self,
n_components: int | None,
*,
algorithm: str = "auto",
tol: float = epsqrt,
max_iter: int = 500,
missing_data_settings: dict | None = None,
):
self.n_components = n_components
self.algorithm = algorithm
self.tol = tol
self.max_iter = max_iter
self.missing_data_settings = missing_data_settings
# ENG-17: the convenience methods (score_plot, vip, spe_limit, ...),
# hotellings_t2_limit, ellipse_coordinates and the rename __getattr__ are
# inherited from _LatentVariableModel. PCA supplies only its rename map.
_ATTRIBUTE_RENAMES: typing.ClassVar[dict[str, str]] = {
"x_scores": "scores_",
"loadings": "loadings_",
"x_loadings": "loadings_",
"squared_prediction_error": "spe_",
"R2": "r2_per_component_",
"R2cum": "r2_cumulative_",
"R2X_cum": "r2_per_variable_",
"hotellings_t2": "hotellings_t2_",
"scaling_factor_for_scores": "scaling_factor_for_scores_",
"N": "n_samples_",
"K": "n_features_in_",
"A": "n_components",
"extra_info": "fitting_info_",
}
_RENAME_CONTEXT: typing.ClassVar[str] = "PCA"
# Fitted diagnostics: per-component arrays (SVD/NIPALS) or scalar totals (TSR).
fitting_info_: dict[str, np.ndarray | int | float]
# ENG-18: public DataFrame views built lazily from the private ndarrays.
scores_ = _LazyFrame("_scores", index="_sample_index", columns="_component_names")
loadings_ = _LazyFrame("_loadings", index="_feature_names", columns="_component_names")
spe_ = _LazyFrame("_spe", index="_sample_index", columns="_component_names")
def __sklearn_tags__(self):
"""Declare sklearn capability tags (sklearn 1.6+).
``allow_nan=True`` because the NIPALS / TSR algorithm paths
thread missing data through; the SVD path raises explicitly
when NaN is supplied. The default for any one ``PCA(...)``
instance therefore allows NaN at the Pipeline / ``check_array``
boundary and lets the fitting algorithm reject it later if
needed (the user-facing error message is more informative than
the sklearn check-array one).
"""
tags = super().__sklearn_tags__()
tags.input_tags.allow_nan = True
return tags
[docs]
def get_feature_names_out(self, input_features=None) -> np.ndarray: # noqa: ANN001, ARG002
"""Return the output column names of :meth:`transform`.
:class:`PCA`'s ``transform`` produces scores, one column per
component, named ``["PC1", "PC2", ..., "PC{n_components}"]``.
The ``input_features`` argument is accepted (Pipeline
introspection passes it through) but unused: the output column
count is the fitted ``n_components``, not the input feature
count.
Used by :meth:`set_output` (sklearn 1.2+) to label the
:class:`~pandas.DataFrame` view of the scores when
``set_output(transform="pandas")`` is on, and by Pipeline
introspection.
"""
check_is_fitted(self, "loadings_")
return np.asarray([f"PC{a}" for a in self._component_names])
[docs]
@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X: DataMatrix, y: DataMatrix | None = None) -> PCA: # noqa: ARG002, PLR0915
"""Fit a principal component analysis (PCA) model to the data.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training data. May contain NaN values for missing data (the
NIPALS / TSR algorithms thread them through; the SVD path
rejects them).
y : ignored
Returns
-------
self : PCA
Raises
------
NotEnoughVarianceError
If every column of X is constant or has no observed value, as across the
family (#1144): such an X tells no row from another. A fit that another
routine of this package runs on a matrix it derived, such as ASCA's effect
matrix of a term with no effect, is not refused.
Warns
-----
UserWarning
If a column has no observed value. It is left out of the model, as
across the family (#1098): its loadings are 0, and a value a new row
has in it adds nothing to that row's scores or SPE.
"""
# Capture the original DataFrame's index/columns before validate_data
# converts X to an ndarray, then rebuild the DataFrame view downstream
# code expects. validate_data also sets n_features_in_ / feature_
# names_in_ and runs the sklearn input rejections (sparse, complex,
# empty, dtype-object) with the standard error messages.
# Reject sparse before validate_data does, so the message names the
# ColumnTransformer knob rather than `.toarray()` (#399).
_reject_sparse(X, "PCA")
refuse_empty_block(X, "X")
# Copies, so that a name the caller gives these axes later does not rename the model's frames
sample_index = X.index.copy() if isinstance(X, pd.DataFrame) else None
feature_columns = X.columns.copy() if isinstance(X, pd.DataFrame) else None
X_arr = validate_data(
self,
X,
reset=True,
accept_sparse=False,
ensure_min_samples=2,
ensure_min_features=1,
dtype="numeric",
ensure_all_finite="allow-nan",
)
if feature_columns is None:
feature_columns = pd.RangeIndex(X_arr.shape[1]) # type: ignore[assignment]
if sample_index is None:
sample_index = pd.RangeIndex(X_arr.shape[0]) # type: ignore[assignment]
X = pd.DataFrame(X_arr, index=sample_index, columns=feature_columns)
N, K = X.shape
self.n_samples_ = N
# n_features_in_ already set by validate_data; reassert for clarity.
self.n_features_in_ = K
self._feature_names = X.columns
self._sample_index = X.index
# A column with nothing observed gets a zero loading from NIPALS and TSR; it is named, as
# across the family, and a value a new row has in it is ignored (#1098).
self._unobserved_columns = unobserved_columns(X, "X")
# An X in which nothing varies tells no row apart (#1144). A matrix another routine of this
# package derived and fits, as ASCA's effect matrix of a term with no effect, is fitted as it was:
# that routine reads its zero scores as the answer.
if _fit_asked_for_by_the_caller(self):
check_x_varies([X], explains_y=False)
constant_columns(X, "X", self)
A = components_to_fit(self.n_components, N, K)
# The resolved (possibly clamped) component count is a fitted attribute;
# the constructor parameter n_components is left exactly as the user set
# it, per the sklearn clone/get_params contract (#505).
self.n_components_ = A
# Detect missing data and resolve algorithm
self.has_missing_data_ = bool(np.any(X.isna()))
algo = self.algorithm.lower()
if algo not in self._valid_algorithms:
raise ValueError(
f"Algorithm '{self.algorithm}' is not recognized. Must be one of {self._valid_algorithms}."
)
if algo == "auto":
algo = "nipals" if self.has_missing_data_ else "svd"
if algo == "svd" and self.has_missing_data_:
raise ValueError("SVD algorithm cannot handle missing data. Use 'nipals', 'tsr', or 'auto'.")
self.algorithm_ = algo
# One resolver for every iterative estimator, so the loop settings and the
# refusal of the removed ``md_*`` keys are defined in a single place (#588). The
# bounds are checked under ``"svd"`` too, where the loop settings govern nothing:
# a value out of range is refused whichever algorithm would read it, as across
# the family (#1071).
settings = resolve_loop_settings(
tol=self.tol,
max_iter=self.max_iter,
missing_data_settings=self.missing_data_settings,
estimator_name="PCA",
)
# Storage for numpy results (set by _fit_* methods)
X_values = np.asarray(X.copy())
# Dispatch. "auto" was resolved and anything unknown refused above, so these three
# branches are every value algo can hold.
if algo == "svd":
self._fit_svd(X_values, N, K, A)
elif algo == "nipals":
self._fit_nipals(X_values, N, K, A, settings)
else: # "tsr"
self._fit_tsr(X_values, N, K, A, settings)
# --- Common post-fit path ---
# ENG-18: scores_ / loadings_ / spe_ are stored as private ndarrays (the
# source of truth); the public DataFrame views are built lazily by the
# _LazyFrame descriptors from these arrays plus the index/column metadata
# (self._sample_index, self._feature_names, self._component_names).
self._component_names = list(range(1, A + 1))
component_names = self._component_names
self._loadings = self._loadings_np
self._scores = self._scores_np
# A row with nothing observed has no SPE, as diagnose gives it none (#1148); fit places it at the centre.
self._spe = np.where(_rows_with_data(X.to_numpy())[:, np.newaxis], self._spe_np, np.nan)
# The training residual block, kept for the trimmed-score-regression
# estimator in ``project`` / ``projection_matrix``: TSR regresses on
# the covariance of the observed columns, which is the model plane
# plus this spread around it. Kept as E (N x K) rather than its K x K
# Gram, because the estimator contracts it to A columns straight away.
# A missing training cell has no residual, so it enters as zero.
self._x_residuals = np.nan_to_num(X_values - self._scores_np @ self._loadings_np.T)
self.r2_per_component_ = pd.Series(
self._r2_np,
index=component_names,
name="R² per component",
)
self.r2_cumulative_ = pd.Series(
self._r2cum_np,
index=component_names,
name="Cumulative R²",
)
self.r2_per_variable_ = pd.DataFrame(
self._r2_per_var_np,
index=self._feature_names,
columns=component_names,
)
self.scaling_factor_for_scores_ = pd.Series(
np.sqrt(self.explained_variance_),
index=component_names,
name="Standard deviation per score",
)
# Hotelling's T² (cumulative across components). A component whose
# score standard deviation is ~0 (a rank-deficient fit, e.g. asking
# for min(N, K) components from mean-centred data, whose rank is at
# most min(N-1, K)) carries no T² information: its standardized score
# is 0/0. _t2_score_sd leaves it out, here and in diagnose and project,
# rather than dividing by ~zero, which poisoned T² with inf/NaN.
self.hotellings_t2_ = pd.DataFrame(
np.cumsum((self._scores / _t2_score_sd(self.scaling_factor_for_scores_.to_numpy())) ** 2, axis=1),
columns=component_names,
index=self._sample_index,
)
# Clean up temporary numpy staging names (the kept ndarrays above alias
# the same arrays, so they survive these del statements).
del self._loadings_np, self._scores_np, self._r2_np, self._r2cum_np, self._r2_per_var_np, self._spe_np
return self
def _fit_svd(self, X_values: np.ndarray, N: int, K: int, A: int) -> None:
"""Fit PCA using SVD decomposition (complete data only)."""
U, S, Vt = np.linalg.svd(X_values, full_matrices=False) # pragma: no mutate - vectors past A are never read
# Loadings are the first A right singular vectors (transposed to K x A),
# scores U * S for the first A components, signed as in every path.
flip = _opposes_convention(Vt[:A].T)
self._loadings_np = np.where(flip, -Vt[:A].T, Vt[:A].T)
self._scores_np = np.where(flip, -U[:, :A] * S[:A], U[:, :A] * S[:A])
# Explained variance; fit() refuses fewer than two rows, so N - 1 >= 1.
self.explained_variance_ = np.diag(self._scores_np.T @ self._scores_np) / (N - 1)
# Compute R2 and SPE via deflation
self._r2cum_np = np.zeros(A)
self._r2_per_var_np = np.zeros((K, A))
self._spe_np = np.zeros((N, A))
Xd = X_values
prior_ssx_col = ssq(Xd, axis=0)
base_variance = np.sum(prior_ssx_col)
for a in range(A):
Xd = Xd - self._scores_np[:, [a]] @ self._loadings_np[:, [a]].T
row_ssx = ssq(Xd, axis=1)
self._spe_np[:, a] = np.sqrt(row_ssx)
# Per-variable R^2 is undefined for a column with no variance to
# explain: dividing by NaN gives NaN there, without the RuntimeWarning
# a division by 0 raises. SEC-21 (#270) sub-item 4.
self._r2_per_var_np[:, a] = 1 - ssq(Xd, axis=0) / np.where(prior_ssx_col > 0, prior_ssx_col, np.nan)
self._r2cum_np[a] = 1 - np.sum(row_ssx) / base_variance if base_variance > 0 else np.nan
self._r2_np = np.diff(self._r2cum_np, prepend=0.0)
self.fitting_info_ = {"timing": np.full(A, np.nan), "iterations": np.full(A, np.nan)}
def _fit_nipals(self, X_values: np.ndarray, N: int, K: int, A: int, settings: dict) -> None:
"""Fit PCA using the NIPALS algorithm (handles missing data)."""
Xd = X_values
base_variance = ssq(Xd)
base_ss_col = ssq(Xd, axis=0) # the denominator of the cumulative per-variable R2, fixed before deflation
self._loadings_np = np.zeros((K, A))
self._scores_np = np.zeros((N, A))
self._r2cum_np = np.zeros(A)
self._r2_per_var_np = np.zeros((K, A))
self._spe_np = np.zeros((N, A))
timing = np.full(A, np.nan)
iterations = np.full(A, np.nan)
for a in range(A):
start_time = time.time()
start_ss_col = ssq(Xd, axis=0)
# Judged against the block's size before the first component, as OPLS judges it, so data in
# small units are not taken for none (#1070). A block of zeros has no first component.
if np.sqrt(np.sum(start_ss_col)) <= epsqrt * np.sqrt(base_variance):
emsg = (
"There is no variance left in the data array: cannot "
f"compute any more components beyond component {a}."
)
raise NotEnoughVarianceError(emsg)
# Seed the score from the column of X with the greatest
# sum-of-squares (variance, for mean-centred data) rather than the
# arbitrary first column (#195). NIPALS converges to the same
# component for any non-degenerate seed, but the highest-variance
# column is closest to the leading component, so it needs fewer
# iterations and is far more robust when the first column happens to
# be near-orthogonal to it. The deterministic sign convention applied
# below makes the fitted sign independent of this seed. A missing
# cell of the seed starts at 0.
t_a = np.nan_to_num(Xd[:, [int(np.argmax(start_ss_col))]]) + 1.0
# The first pass always runs: the ``+ 1.0`` offset that primes the
# loop can be negligible relative to a large-magnitude seed column,
# and the relative criterion (#504) would then report convergence at
# entry, silently returning the unrefined seed as the score.
# terminate_check ends the loop by md_max_iter at the latest.
for itern in itertools.count(1):
t_a_guess = t_a
# Regress X onto t_a to get loadings p_a
p_a = quick_regress(Xd, t_a, per="column")
# Floor the norm: a collapsed ``p_a`` would make this 0/0 -> NaN and
# poison every later component. `_mbpca` already guards the identical
# expression this way; see `_nz` (#513).
p_a = p_a / _nz(float(np.sqrt(ssq(p_a))))
# Regress X onto p_a to get scores t_a
t_a = quick_regress(Xd, p_a, per="row")
if terminate_check(t_a_guess, t_a, iterations=itern, settings=settings):
break
timing[a] = time.time() - start_time
iterations[a] = itern
# Whether the last pass met the tolerance, not whether it was the last one allowed: a
# component can converge on that pass (#1080).
converged = scores_converged(t_a_guess, t_a, settings["md_tol"])
if not converged:
warnings.warn(
f"PCA NIPALS: component {a + 1} reached the maximum number of "
f"iterations ({settings['md_max_iter']}) without converging.",
SpecificationWarning,
stacklevel=user_stacklevel(),
)
logger.debug(
"PCA NIPALS: component %d %s in %d iterations (tol=%g)",
a + 1,
"converged" if converged else "did not converge",
itern,
settings["md_tol"],
)
# Deflate
Xd = Xd - np.dot(t_a, p_a.T)
row_ssx = ssq(Xd, axis=1)
self._spe_np[:, a] = np.sqrt(row_ssx)
# Per-variable R^2 is undefined for a column with no variance to
# explain: dividing by NaN gives NaN there. SEC-21 (#270) sub-item 4.
self._r2_per_var_np[:, a] = 1 - ssq(Xd, axis=0) / np.where(base_ss_col > 0, base_ss_col, np.nan)
# base_variance is above 0: the first pass raised otherwise.
self._r2cum_np[a] = 1 - np.sum(row_ssx) / base_variance
flip = _opposes_convention(p_a)
self._loadings_np[:, a] = np.where(flip, -p_a, p_a).ravel()
self._scores_np[:, a] = np.where(flip, -t_a, t_a).ravel()
self._r2_np = np.diff(self._r2cum_np, prepend=0.0)
self.fitting_info_ = {"timing": timing, "iterations": iterations}
# Explained variance; fit() refuses fewer than two rows, so N - 1 >= 1.
self.explained_variance_ = np.diag(self._scores_np.T @ self._scores_np) / (N - 1)
def _fit_tsr(self, X_values: np.ndarray, N: int, K: int, A: int, settings: dict) -> None:
"""Fit PCA using the Trimmed Score Regression algorithm (handles missing data).
See papers by Abel Folch-Fortuny and also DOI: 10.1002/cem.750
"""
start_time = time.time()
Xd = X_values.copy()
base_variance = ssq(Xd)
mmap = np.isnan(Xd)
Xd[mmap] = 0.0
# No missing cells means nothing to impute: the EM loop does not run
# (``np.mean`` over the empty ``Xd[mmap]`` used to emit a RuntimeWarning
# every iteration and produce a NaN delta).
itern = 0
for itern in range(1, (settings["md_max_iter"] if mmap.any() else 0) + 1): # noqa: B007 - read after the loop
missing_X = Xd[mmap]
mean_X = np.mean(Xd, axis=0)
# The regression below is unchanged by the scale of S.
S = np.cov(Xd, rowvar=False, ddof=1) # pragma: no mutate - np.cov reads rowvar=None and ddof=None as given
Xc = Xd - mean_X
# svd(Xc) holds the loadings in the ROWS of Vt, whatever the shape
# of Xc, so V is (K, A): rows indexed by feature, columns by component.
_, _, Vt = np.linalg.svd(Xc, full_matrices=False) # pragma: no mutate - rows past A are never read
V = Vt[:A].T
for n in range(N):
row_mis = mmap[n, :]
if np.any(row_mis):
row_obs = ~row_mis
L = V[row_obs, 0 : min(A, sum(row_obs))]
S11 = S[row_obs, :][:, row_obs]
S21 = S[row_mis, :][:, row_obs]
z2 = (S21 @ L) @ np.linalg.pinv(L.T @ S11 @ L) @ L.T
Xc[n, row_mis] = z2 @ Xc[n, row_obs]
Xd = Xc + mean_X
delta = np.mean((Xd[mmap] - missing_X) ** 2)
if delta <= settings["md_tol"]: # pragma: no mutate - < and <= differ only at a change of exactly md_tol
break
# Final decomposition. The singular vectors of the covariance do not depend on its scale.
S = np.cov(Xd, rowvar=False, ddof=1) # pragma: no mutate - np.cov reads rowvar=None and ddof=None as given
_, _, V = np.linalg.svd(S) # S is square, so its reduced and full SVD are one
# Sign convention: flip so the largest-magnitude element in each
# loading is positive, matching _fit_svd / _fit_nipals (previously
# the TSR loading signs were LAPACK-dependent).
self._loadings_np = np.where(_opposes_convention(V[:A].T), -V[:A].T, V[:A].T) # K x A
# Project the imputed matrix as-is (no extra centring), exactly as
# transform() will project new data: the package convention is that X
# is preprocessed (e.g. MCUVScaler) before fit. The previous code
# centred here, so scores_ disagreed with transform(X) on the very
# same data, and the residuals below mixed centred reconstructions
# with uncentred data (inflating SPE and deflating R2).
self._scores_np = Xd @ self._loadings_np
# R2 and SPE
self._r2cum_np = np.zeros(A)
self._r2_per_var_np = np.zeros((K, A))
self._spe_np = np.zeros((N, A))
base_ss_col = ssq(X_values, axis=0)
for a in range(A):
residuals = self._scores_np[:, : a + 1] @ self._loadings_np[:, : a + 1].T - X_values
self._r2cum_np[a] = 1 - ssq(residuals) / base_variance
self._r2_per_var_np[:, a] = 1 - ssq(residuals, axis=0) / np.where(base_ss_col > 0, base_ss_col, np.nan)
self._spe_np[:, a] = np.sqrt(ssq(residuals, axis=1))
self._r2_np = np.diff(self._r2cum_np, prepend=0.0)
self.fitting_info_ = {"iterations": itern, "timing": time.time() - start_time}
# Explained variance; fit() refuses fewer than two rows, so N - 1 >= 1.
self.explained_variance_ = np.diag(self._scores_np.T @ self._scores_np) / (N - 1)
def _validated_rows(self, X: DataMatrix) -> tuple[np.ndarray, pd.Index | None]:
"""Return ``X`` checked against the fit, its columns in the fitted order, and its row labels.
Every method that takes new rows reads them through here, so a DataFrame's
columns are matched by name and a pandas missing value is a NaN in all of them.
Parameters
----------
X : array-like of shape (n_samples, n_features)
New rows; NaN marks a missing cell.
Returns
-------
X_arr : np.ndarray of shape (n_samples, n_features)
The validated rows, columns in the order of the fit.
sample_index : pd.Index or None
The DataFrame's row labels, or None for an array, which pandas numbers from 0.
"""
check_is_fitted(self, "loadings_")
# sklearn's validate_data(reset=False) checks feature-name *order*
# strictly, while _align_to_fit_features only needs set equality.
# Realign reordered DataFrame columns first so validate_data sees a
# name-ordered view; ndarrays / un-named DataFrames pass straight
# through.
if isinstance(X, pd.DataFrame):
X = _align_to_fit_features(X, self._feature_names)
X_arr = validate_data(
self,
X,
reset=False, # pragma: no mutate - validate_data reads reset=None as False
ensure_all_finite="allow-nan",
)
# A column the fit left out (#1098) is read at 0, so it adds nothing to a score or an SPE, and
# a NaN there does not make every score NaN. A model pickled before then left none out.
X_arr = np.where(
getattr(
self,
"_unobserved_columns",
False, # pragma: no mutate - np.where reads a condition of None as False
),
0.0,
X_arr,
)
# pandas numbers the rows of an array from 0 when given no index.
return X_arr, X.index if isinstance(X, pd.DataFrame) else None
def _scaled_rows(self, X: DataMatrix) -> pd.DataFrame:
"""Return :meth:`_validated_rows` as a frame labelled as fitted, the name :meth:`PLS._scaled_rows` has.
PCA neither centres nor scales, so the rows are in the space it was fitted in as they come.
The contribution functions read every model's new rows through this one name.
"""
X_arr, sample_index = self._validated_rows(X)
return pd.DataFrame(X_arr, index=sample_index, columns=self._feature_names)
[docs]
def diagnose(self, X: DataMatrix) -> Bunch:
"""Project new data and compute diagnostics (scores, Hotelling's T², SPE).
This logic was called ``predict`` until 1.38.1 (#396). The rename
matches :meth:`PLS.diagnose`, and the deprecated ``predict`` was
removed in 2.0: PCA is not a regressor, so it has no ``predict``.
Parameters
----------
X : array-like of shape (n_samples, n_features)
New rows. A row with a missing cell (NaN) is scored from the cells
it has, by trimmed score regression, the estimator :meth:`project`
uses by default, and its SPE is over those cells; until 2.4.0 its
scores, T² and SPE were NaN (#1144). A row with nothing observed
stays NaN. So does a row with too few observed cells for its
scores, typically fewer than ``n_components``, with a
:class:`SpecificationWarning`. :meth:`project` gives the other
estimators, a ``ridge``, and each row's conditioning.
Returns
-------
result : sklearn.utils.Bunch
With keys:
- ``scores`` (pd.DataFrame, ``(n_samples, n_components)``).
- ``cumulative_hotellings_t2`` (pd.DataFrame, ``(n_samples, n_components)``):
Hotelling's T², cumulative, so column ``a`` holds the T² over the first ``a``
components and the last column is the full model's, as in ``hotellings_t2_``.
Every ``diagnose`` of the family returns it in this shape (#983).
- ``spe`` (pd.Series, ``(n_samples,)``): SPE of the full model.
- ``hotellings_t2``: deprecated since 2.4.0, to be removed in 3.0;
the same DataFrame as ``cumulative_hotellings_t2``. Reading it warns.
"""
X_arr, sample_index = self._validated_rows(X)
scores = X_arr @ self._loadings
# A row with a missing cell is scored from the cells it has, as project() scores it by default (#1144)
rows, estimates = self._estimate_rows_with_gaps(X_arr)
scores[rows] = estimates.scores
# Hotelling's T² (cumulative), leaving out a component fit leaves out (#1026)
score_sd = _t2_score_sd(self.scaling_factor_for_scores_.to_numpy())
t2 = pd.DataFrame(
np.cumsum((scores / score_sd) ** 2, axis=1), index=sample_index, columns=self._component_names
)
# SPE: residual after reconstruction, over the cells a row has
residuals = X_arr - scores @ self._loadings.T
spe = np.sqrt(np.sum(residuals**2, axis=1))
spe[rows] = estimates.spe
return DiagnoseResult(
f"multivariate.{type(self).__name__}",
scores=pd.DataFrame(scores, index=sample_index, columns=self._component_names),
cumulative_hotellings_t2=t2,
hotellings_t2=t2,
spe=pd.Series(spe, index=sample_index, name="SPE"),
)
[docs]
def project(
self,
X: DataMatrix,
*,
method: str = "tsr", # pragma: no mutate - _validate_method lowercases it, so "TSR" is the same default
ridge: float = 0.0,
) -> Bunch:
"""Estimate scores and diagnostics for rows that may contain missing values.
This method estimates the scores of partially-observed rows from the
observed columns only, using the missing-data estimators of Arteaga
and Ferrer (2002): trimmed score regression (``"tsr"``, the default
and the statistically strongest), single-component projection
(``"scp"``), or projection to the model plane (``"pmp"``). Rows with
no missing values take the standard complete-data path, so their
scores are bitwise identical to :meth:`transform`. :meth:`transform`
and :meth:`diagnose` score a row with a missing cell by this method's
default, through the same code; this method adds the choice of
estimator, the ``ridge``, and each row's ``condition_number``, and
raises where they give NaN for a row the estimator cannot be inverted
for (#1144).
This is the "batch so far" primitive of online batch monitoring: the
future part of an unfolded batch row is missing by construction, and
the score estimate at each time sample is this projection (see
Garcia-Munoz, Kourti and MacGregor, 2004).
Parameters
----------
X : array-like of shape (n_samples, n_features)
New data in the model's (centred and scaled) space; NaN marks a
missing entry. A row with nothing observed sits at the model
centre, as :meth:`fit` places one: scores, T2 and SPE of 0, no
observed feature, and a ``condition_number`` of inf.
method : {"tsr", "scp", "pmp"}, default="tsr"
The score estimator; see
:mod:`process_improve.multivariate._projection`.
ridge : float, default=0.0
Non-negative regularisation added to the matrix inverted by the
``"tsr"`` and ``"pmp"`` estimators. Raise it above zero when
``condition_number`` reports near-singularity (typically very
early in a batch, when few columns are observed).
Returns
-------
result : sklearn.utils.Bunch
With keys ``scores`` (DataFrame, n_samples x n_components),
``hotellings_t2`` (Series; total over all components, computed
with the training score variances), ``spe`` (Series; the square
root of the residual sum of squares over the *observed* columns
only), ``condition_number`` (Series; the conditioning diagnostic
of each row's estimator, 1.0 when nothing is missing) and
``n_observed`` (Series; observed features per row). SPE and T2
for a partially-observed row must be compared against limits
built from the same missingness pattern, not against the
full-observation limits; see
:class:`process_improve.batch.BatchMonitor`.
"""
X_arr, sample_index = self._validated_rows(X)
raw = project_rows(x_scaled=X_arr, method=method, ridge=ridge, **self._projection_inputs())
scores = pd.DataFrame(raw.scores, index=sample_index, columns=self._component_names)
# A component fit leaves out of T² is left out here too (#1026).
s = _t2_score_sd(self.scaling_factor_for_scores_.to_numpy())
t2 = pd.Series(np.sum((raw.scores / s) ** 2, axis=1), index=sample_index, name="Hotelling's T2")
return Bunch(
scores=scores,
hotellings_t2=t2,
spe=pd.Series(raw.spe, index=sample_index, name="SPE"),
condition_number=pd.Series(raw.condition_number, index=sample_index, name="condition_number"),
n_observed=pd.Series(raw.n_observed, index=sample_index, name="n_observed"),
)
[docs]
def projection_matrix(
self,
observed: object,
*,
method: str = "tsr", # pragma: no mutate - _validate_method lowercases it, so "TSR" is the same default
ridge: float = 0.0,
) -> Bunch:
"""Build the fixed linear operator mapping observed columns to score estimates.
For a fixed missingness pattern, every estimator in :meth:`project`
is a fixed linear map ``t_hat = M @ z_observed``. This method exposes
that matrix so callers that reuse one pattern many times (an online
monitor at time sample ``k``, or an optimiser treating candidate
columns as observed) can precompute it once.
Parameters
----------
observed : array-like
Either a boolean mask of length ``n_features_in_`` (True =
observed), or a list of feature labels to treat as observed.
method : {"tsr", "scp", "pmp"}, default="tsr"
ridge : float, default=0.0
Returns
-------
result : sklearn.utils.Bunch
With keys ``matrix`` (DataFrame, n_components x n_observed,
columns labelled by the observed features), ``condition_number``
(float) and ``method``.
"""
check_is_fitted(self, "loadings_")
mask = coerce_observed_mask(observed, self._feature_names)
op = operator_for_pattern(observed=mask, method=method, ridge=ridge, **self._projection_inputs())
matrix = pd.DataFrame(
op.matrix,
index=self._component_names,
columns=pd.Index(self._feature_names)[mask],
)
return Bunch(matrix=matrix, condition_number=op.condition_number, method=op.method)
def _projection_inputs(self) -> dict[str, typing.Any]:
"""Return the fitted matrices the missing-data estimators of :meth:`project` read, by their argument names.
For PCA the loadings both generate the scores and reconstruct the rows, and are the
weights single-component projection projects onto, the default of the estimators.
"""
return {
"x_loadings": self._loadings,
"guide": self._loadings,
"score_variances": self.explained_variance_,
"x_residuals": getattr(self, "_x_residuals", None),
}
[docs]
def score(self, X: DataMatrix, y: DataMatrix | None = None) -> float: # noqa: ARG002
"""Negative mean squared reconstruction error (higher is better).
Follows the sklearn convention where higher scores indicate better
model fit. This makes PCA compatible with ``cross_val_score``,
``GridSearchCV``, and other sklearn model-selection utilities.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Test data to score.
y : ignored
Returns
-------
score : float
Negative mean squared reconstruction error.
Examples
--------
>>> from sklearn.model_selection import cross_val_score
>>> scores = cross_val_score(PCA(n_components=2), X_scaled, cv=5)
>>> print(f"Mean CV score: {scores.mean():.4f}")
"""
# The rows as transform() reads them: a DataFrame's columns by name, not by position.
X_arr, _ = self._validated_rows(X)
residuals = X_arr - (X_arr @ self._loadings) @ self._loadings.T
return -float(np.mean(residuals**2))
[docs]
@classmethod
def minka_mle(cls, X: DataMatrix) -> int:
"""Minka (2000) automatic-dimensionality estimate for PCA.
Closed-form Bayesian model selection on the PPCA evidence (Minka,
T. P. 2000. *Automatic Choice of Dimensionality for PCA*. NIPS 13,
pp. 598-604). Operates only on the covariance eigenvalues of ``X``
and is therefore very cheap; in the simulations Minka reports it
beats cross-validation. Use it alongside the ekf-CV recommendation
from :meth:`select_n_components` as a fast cross-check.
Internally ``X`` is mean-centred before estimation (a PPCA
assumption); it is **not** unit-variance scaled, because dividing
each column by its standard deviation compresses the noise
eigenvalues to near-zero values the MLE misreads as additional
latent signal. If your columns are on wildly different scales,
pass the analysis-scale ``X`` produced by your own preprocessing
(e.g. SNV for spectral data) and accept the centring this method
applies.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Data matrix.
Returns
-------
n_components : int
The MLE estimate of the effective dimensionality.
References
----------
Minka, T. P. (2000). Automatic Choice of Dimensionality for PCA.
Advances in Neural Information Processing Systems, 13, 598-604.
See Also
--------
parallel_analysis : Horn (1965) eigenvalue-vs-null retention.
select_n_components : ekf cross-validation; pass
``return_consensus=True`` to report all three side by side.
"""
if not isinstance(X, pd.DataFrame):
X = pd.DataFrame(X) # a 1-D X is one column
# sklearn centres X itself, and converts it to floats. ``svd_solver="full"``
# is the only solver that supports ``n_components="mle"`` in sklearn.
sk = _SkPCA(
n_components="mle",
svd_solver="full", # pragma: no mutate - "auto" picks a full decomposition for "mle" too
)
sk.fit(X)
return int(sk.n_components_)
[docs]
@classmethod
def parallel_analysis(
cls,
X: DataMatrix,
*,
n_simulations: int = 200,
quantile: float = 0.95,
surrogate: typing.Literal["normal", "permutation"] = "normal",
scale: bool = True,
random_state: int | None = None,
) -> Bunch:
"""Horn (1965) parallel analysis component-count estimate.
Generates ``n_simulations`` random matrices of the same shape as
``X``, computes their eigenvalues, and retains every observed
component whose eigenvalue exceeds the ``quantile`` of the null
distribution at the same rank. Widely regarded in psychometrics
as the best simple retention rule for PCA.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Data matrix.
n_simulations : int, default 200
Number of random matrices drawn to build the null
eigenvalue distribution.
surrogate : {"normal", "permutation"}, default "normal"
How the null matrices are built.
- ``"normal"``: independent standard-normal entries, which is
Horn's original proposal. Fast, and exactly right when the
columns really are Gaussian.
- ``"permutation"``: each column of the real data is permuted
independently (Buja and Eyuboglu, 1992). This breaks the
correlation between columns, which is what parallel analysis
is testing for, while leaving every column's own distribution
untouched. Prefer it on process data, where a tag may be
skewed, heavy-tailed, bounded at zero or quantised by its
instrument: a Gaussian null then answers a question about
Gaussian data rather than about this block.
quantile : float, default 0.95
Quantile of the null eigenvalues used as the retention
threshold. Horn's original proposal was the mean (0.5);
the more conservative 95th-percentile threshold is the
modern recommendation.
scale : bool, default True
Mean-centre and unit-variance scale ``X`` before estimation.
Unlike :meth:`minka_mle`, which is mean-centred but never
unit-variance scaled, parallel analysis defaults to autoscaling
here so wildly different column scales do not dominate the null
comparison.
random_state : int, optional
Seed for the null-matrix simulations.
Returns
-------
result : sklearn.utils.Bunch
With keys:
- ``n_components`` - number of components retained (can be 0
on pure noise).
- ``observed_eigenvalues`` - eigenvalues of ``X`` after
centring/scaling (np.ndarray of length ``min(n, p)``).
- ``null_threshold`` - per-rank ``quantile`` of the null
eigenvalue distribution (same length as
``observed_eigenvalues``). Plot it against
``observed_eigenvalues`` for the scree-versus-null picture the
method is usually read from.
- ``surrogate`` - which null was used, echoed back so a stored
result says how it was produced.
References
----------
Horn, J. L. (1965). A rationale and test for the number of
factors in factor analysis. *Psychometrika*, 30(2), 179-185.
Buja, A., & Eyuboglu, N. (1992). Remarks on parallel analysis.
*Multivariate Behavioral Research*, 27(4), 509-540.
See Also
--------
minka_mle : closed-form PPCA evidence rule.
select_n_components : ekf cross-validation; pass
``return_consensus=True`` to report all three side by side.
"""
if not isinstance(X, pd.DataFrame):
X = pd.DataFrame(X) # a 1-D X is one column
if scale:
X = MCUVScaler().fit_transform(X)
X_arr = np.asarray(X, dtype=float)
n, p = X_arr.shape
# Observed eigenvalues from the centred X. Centring removes one DoF
# so the smallest singular value is at-or-near zero; that's expected.
observed = np.linalg.svdvals(X_arr - X_arr.mean(axis=0)) ** 2 / max(1, n - 1)
known_surrogates = ("normal", "permutation")
if surrogate not in known_surrogates:
msg = f"surrogate must be one of {known_surrogates}; got {surrogate!r}."
raise ValueError(msg)
rng = np.random.default_rng(random_state)
null_eigs = np.zeros((n_simulations, observed.size))
for i in range(n_simulations):
if surrogate == "normal":
R = rng.standard_normal((n, p))
else:
# Permute each column on its own. Shuffling whole rows would leave the
# correlation structure exactly as it is and produce no null at all.
R = np.column_stack([rng.permutation(X_arr[:, j]) for j in range(p)])
null_eigs[i] = np.linalg.svdvals(R - R.mean(axis=0)) ** 2 / max(1, n - 1)
null_threshold = np.quantile(null_eigs, quantile, axis=0)
# Standard PA: retain consecutive components from the top while
# observed > null, so the count is the length of that leading run.
# Components past the first failure are not retained even if their
# observed eigenvalue happens to exceed the null (a rare numerical
# coincidence on real data).
above = observed > null_threshold # pragma: no mutate - > and >= differ only at an exact tie
n_retained = int(np.cumprod(above).sum())
return Bunch(
n_components=n_retained,
observed_eigenvalues=observed,
null_threshold=null_threshold,
surrogate=surrogate,
)
[docs]
@classmethod
def select_n_components( # noqa: PLR0913, PLR0915
cls,
X: DataMatrix,
*,
max_components: int | None = None,
cv: int | BaseCrossValidator = 5,
cv_scheme: typing.Literal["ekf", "ckf", "ek", "sacv", "gcv"] = "ekf",
n_repeats: int = 1,
selection_rule: SelectionRule = "min",
min_q2_increase: float = Q2_MIN_INCREMENT,
scale_inside_folds: bool = True,
ekf_max_iter: int = 50,
ekf_tol: float = 1e-6,
random_state: int | None = None,
return_consensus: bool = False,
) -> Bunch:
"""Select the number of PCA components via cross-validation.
Evaluates every component count ``1, 2, ..., max_components`` and
recommends one via the configured ``selection_rule``. The default
``cv_scheme="ekf"`` is the element-wise k-fold algorithm of Bro,
Kjeldahl, Smilde & Kiers (2008, *Anal. Bioanal. Chem.* 390:1241-1251),
which holds out individual cells of ``X`` and predicts them via
EM-style imputation from a model that never sees their true values,
so PRESS does not shrink monotonically with the component count the
way it does when a held-out row is scored through its own values.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training data. With the default ``scale_inside_folds=True`` pass
the raw, unscaled X; mean-centring and unit-variance scaling are
fit on each fold's in-fold cells. Pre-scale it yourself only
with ``scale_inside_folds=False``.
max_components : int, optional
Maximum number of components to evaluate. Default is
``min(n_samples - 1, n_features)``.
cv : int or sklearn CV splitter, default 5
For ``cv_scheme="ekf"``: the integer number of element-folds
(splitter objects are ignored). For ``cv_scheme="ek"``: the number
of row groups, and of column groups. Ignored by ``"sacv"`` and
``"gcv"``, which hold nothing out.
cv_scheme : {"ekf", "ckf", "ek", "sacv", "gcv"}, default "ekf"
How a held-out value is produced. Every one of these keeps the
prediction independent of the value being predicted; they differ in
what they hold out and what they cost. ``"row_wise"`` was removed
in 2.0 and raises a :class:`ValueError`: it scored a held-out row
through its own values, so it measured compression rather than
prediction and could not pick a component count.
- ``"ekf"``, the default: element-wise k-fold with EM imputation.
Scattered cells are held out and each is imputed from a model
that never saw it. Fits its centring and scaling inside every
fold, and is the only scheme here that takes a block with
missing cells. Cost: ``n_folds * n_repeats * max_components``
decompositions.
- ``"ckf"``: column-wise k-fold. Groups of columns are held out,
and their values are predicted from scores computed on the
retained columns only, so no held-out value appears in the score
that predicts it. Information still reaches the model through the
loadings, which are fitted on the whole block, and that is the
respect in which ``"ekf"`` is stricter. Offered because a great
deal of published chemometrics uses it, so a number that has to
line up with a paper may need it, and because on a block with far
more columns than samples it costs one decomposition rather than
``n_folds * n_repeats * max_components``. A component count that
leaves some group with fewer retained columns than components
cannot be scored there, so its PRESS and Q2 are NaN and it is
not recommended.
- ``"ek"``: the two-model scheme of 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. This is
what Simca-P reports and what ``pcaMethods::Q2`` computes by
default, so use it when a number has to line up with either.
Cost: ``2 * n_folds`` decompositions.
- ``"sacv"``: leave-one-cell-out, approximated by inflating each
residual by the leverage of the cell that produced it, after
Josse and Husson (2012). The cheap version of holding cells out
one at a time. Cost: one decomposition.
- ``"gcv"``: the same idea with a single averaged leverage instead
of one per cell. Blunter, and it tends to keep more components.
Cost: one decomposition. This and ``"sacv"`` are the defaults in
``FactoMineR`` and ``missMDA``.
``"ek"``, ``"sacv"`` and ``"gcv"`` factorise the matrix directly,
so they raise on a block with missing cells and they ignore
``scale_inside_folds``, ``n_repeats``, ``ekf_max_iter`` and ``ekf_tol``.
They also report no per-fold spread, so ``selection_rule="1se"``
has nothing to work with; ``"min"`` takes the global optimum,
where ``FactoMineR`` instead stops at the first local worsening,
which can return a smaller count on the same data.
n_repeats : int, default 1
Repeat the ekf pass with a fresh random fold permutation this
many times. Each repeat covers every cell exactly once;
``n_repeats > 1`` narrows the per-component PRESS standard
error (helpful when the 1-SE rule sits on a borderline) at
roughly linear extra runtime.
selection_rule : {"min", "1se", "q2_increment"}, default "min"
How the recommended component count is chosen. ``"min"`` is the
GlobalMin criterion Bro 2008 pairs with ekf - the component
count with the lowest pooled PRESS. ``"1se"`` is the one-
standard-error rule (needs the per-fold spread that ``"ekf"`` and
``"ckf"`` report). ``"q2_increment"`` is the Wold's-R-style
cumulative-:math:`Q^2` threshold from PR #371; ``min_q2_increase``
sets the threshold.
min_q2_increase : float, default 0.01
Threshold used only when ``selection_rule="q2_increment"``.
scale_inside_folds : bool, default True
With the default, mean-centring and unit-variance scaling
constants are fit on each fold's in-fold cells and applied to
the whole matrix before EM, removing the centring/scaling
leakage of the prior implementation. Set to ``False`` to
reproduce the previous behaviour (column mean recomputed each
EM iteration from the imputed matrix, no scaling); this is
useful only when ``X`` is already pre-scaled, and a
:class:`SpecificationWarning` is emitted because scaling
constants fit on the full matrix leak into every element-fold.
Pass the **raw, unscaled** X under the default. In-fold
re-standardisation overwrites whatever scaling the caller
applied, so two deliberately different strategies (autoscale
versus Pareto, say) become the same model and report the same
PRESS: a comparison between them shows no difference for
reasons that have nothing to do with the data. A
:class:`SpecificationWarning` is emitted when ``X`` arrives
already centred and unit-variance scaled, which is the
detectable half of that case; a block scaled some other way
cannot be recognised, so the rule is the caller's to keep.
Same contract as :meth:`PLS.select_n_components`.
ekf_max_iter, ekf_tol : int and float, default 50 and 1e-6
EM iteration cap and convergence tolerance for the ekf imputation
step. Named for the loop they bound: that EM loop is not the
per-component NIPALS loop, which ``PCA(tol=..., max_iter=...)``
governs (#588). Their old names ``n_iter`` and ``tol`` were removed
in 2.0.
random_state : int, optional
Seed for the ekf element-fold permutation.
return_consensus : bool, default False
When ``True``, also cross-check the CV recommendation against
two cheap alternative selectors: Minka's PPCA MLE
(:meth:`minka_mle`) and Horn's parallel analysis
(:meth:`parallel_analysis`). The result Bunch then gains the
``minka_n_components``, ``parallel_analysis_n_components``,
``consensus``, and ``consensus_counts`` keys (see Returns). The
parallel analysis runs with its default ``scale=True`` whatever
``scale_inside_folds`` is.
Returns
-------
result : sklearn.utils.Bunch
With keys:
- ``n_components`` - recommended number of components (int).
- ``press`` - pooled PRESS per component count (pd.Series,
indexed ``1..A_max``). Under ``cv_scheme="ekf"`` with
``scale_inside_folds=True`` this is measured in the space each
fold was fitted in, so every variable weighs the same; see
``press_input_units`` for the other scale.
- ``press_input_units`` - the same curve in the units of the
matrix that was passed in, for comparing prediction error
against instrument error (pd.Series, indexed ``1..A_max``).
- ``per_fold_press`` - per-fold PRESS contributions
(pd.DataFrame, ``A_max`` rows x ``n_folds * n_repeats`` columns
under ekf).
- ``se_press`` - standard error of the per-fold PRESS curve
(pd.Series, indexed ``1..A_max``). Drives the 1-SE rule.
- ``q2_se`` - the same standard error rescaled onto the Q2 scale
(``se_press / null_model_ss``, pd.Series indexed ``1..A_max``),
i.e. the half-width of a +/-1 SE band around ``q2``.
- ``press_ratio`` - ``PRESS_a / PRESS_{a-1}`` for inspection
(pd.Series, indexed ``2..A_max``).
- ``q2`` - cross-validated :math:`R^2_X` per component count
(pd.Series, indexed ``1..A_max``). Computed as
``1 - press / null_model_ss``, where the null model predicts
each held-out cell by its in-fold column mean, measured the
same way ``press`` is. Directly comparable to
``r2_cumulative_`` and to PLS's ``r2y_validated``.
- ``q2_per_variable`` - that same quantity split by variable
(pd.DataFrame, ``A_max`` rows x ``K`` columns), which is what
shows whether one column is carrying the pooled figure.
- ``cv_scores`` - alias of ``per_fold_press``.
- ``cv_scheme`` - the scheme used, echoed back (one of
``"ekf"``, ``"ckf"``, ``"ek"``, ``"sacv"`` or ``"gcv"``).
- ``selection_rule`` - the rule used to pick ``n_components``.
When ``return_consensus=True``, the Bunch additionally carries:
- ``minka_n_components`` - the Minka PPCA MLE estimate (int).
- ``parallel_analysis_n_components`` - Horn's parallel-analysis
estimate (int).
- ``consensus`` - ``"agree"`` if the three integer estimates
(CV recommendation, Minka, parallel analysis) span at most
1, otherwise ``"disagree"``.
- ``consensus_counts`` - the tuple
``(recommended, minka_n, parallel_analysis_n)``.
References
----------
Bro, R., Kjeldahl, K., Smilde, A. K., & Kiers, H. A. L. (2008).
Cross-validation of component models: a critical look at current
methods. *Anal. Bioanal. Chem.*, 390(5), 1241-1251.
Camacho, J., & Ferrer, A. (2012). Cross-validation in PCA models
with the element-wise k-fold (ekf) algorithm: theoretical aspects.
*J. Chemometrics*, 26(7), 361-373.
Eastment, H. T., & Krzanowski, W. J. (1982). Cross-validatory choice
of the number of components from a principal component analysis.
*Technometrics*, 24(1), 73-77.
Josse, J., & Husson, F. (2012). Selecting the number of components in
principal component analysis using cross-validation approximations.
*Computational Statistics & Data Analysis*, 56(6), 1869-1879.
"""
X = _float_frame(X) # an array is numbered from 0; a nullable column's pd.NA is a missing cell, as in fit
check_at_least_one("max_components", max_components)
N, K = X.shape
max_components = min(N - 1, K) if max_components is None else min(int(max_components), N - 1, K)
if max_components < 1:
raise ValueError("No components can be evaluated; the data is too small.")
component_index = pd.Index(range(1, max_components + 1), name="n_components")
X_arr = np.asarray(
X,
dtype=float, # pragma: no mutate - _float_frame has made X a frame of floats
)
# The number of folds of "ekf", of column groups of "ckf", and of row and column
# groups of "ek". A splitter object is not read: the default is used.
n_folds = cv if isinstance(cv, int) else 5
# ``press`` under ekf is the TOTAL PRESS per repeat (n_folds times the
# mean per-fold PRESS); the per-fold SE computed further below must be
# rescaled by this factor so the 1-SE rule compares like with like.
press_scale_multiplier = 1
if cv_scheme == "ekf":
press_scale_multiplier = n_folds
check_at_least_one("n_repeats", n_repeats)
# Same two traps as PLS.select_n_components, checked only here
# because the other schemes ignore the flag.
_warn_scaling_traps(X, scale_inside_folds=scale_inside_folds, fold="element-fold", metric="PRESS")
ekf = _pca_ekf_press(
X_arr,
max_components,
n_folds=n_folds,
n_repeats=n_repeats,
ekf_max_iter=ekf_max_iter,
ekf_tol=ekf_tol,
scale_inside_folds=scale_inside_folds,
random_state=random_state,
)
raw_press, raw_per_fold, raw_per_column = ekf.press, ekf.per_fold_press, ekf.per_column_press
raw_press_input_units = ekf.press_input_units
null_model_ss, per_column_null = _ekf_null_reference(ekf, X_arr, scale_inside_folds=scale_inside_folds)
elif cv_scheme == "ckf":
raw_press, raw_per_fold, raw_per_column, null_model_ss, per_column_null = _pca_ckf_press(
_preprocess_for_cell_schemes(X_arr, cv_scheme),
max_components,
n_folds=n_folds,
random_state=random_state,
)
# ckf scales the whole block once, like its cell-scheme siblings, so the PRESS
# it reports is already in the scaled space and there is no separate one to
# report in the caller's units.
raw_press_input_units = raw_press
elif cv_scheme in {"ek", "sacv", "gcv"}:
Z = _preprocess_for_cell_schemes(X_arr, cv_scheme)
if cv_scheme == "ek":
raw_press, raw_per_column, null_model_ss, per_column_null = _eastment_krzanowski_press(
Z, max_components, n_folds=n_folds, random_state=random_state
)
else:
raw_press, raw_per_column, null_model_ss, per_column_null = _leverage_corrected_press(
Z, max_components, method=cv_scheme
)
# None of these three splits the data into folds that could disagree:
# "ek" pools every cell into one total, and the other two hold nothing
# out at all. A per-fold spread would be fabricated, so it is absent
# rather than zero, and the 1-SE rule has nothing to work with.
raw_per_fold = np.full((max_components, 1), np.nan)
raw_press_input_units = raw_press
else:
removed = (
" 'row_wise' was removed in 2.0: it scored a held-out row through its own values, so it "
"could not pick a component count."
if cv_scheme == "row_wise"
else ""
)
raise ValueError(
f"Unknown cv_scheme {cv_scheme!r}; expected one of 'ekf', 'ckf', 'ek', 'sacv' or 'gcv'.{removed}"
)
press = pd.Series(raw_press, index=component_index, name="PRESS")
per_fold_press = pd.DataFrame(
raw_per_fold,
index=component_index,
columns=[f"fold_{i + 1}" for i in range(raw_per_fold.shape[1])],
)
press_input_units = pd.Series(raw_press_input_units, index=component_index, name="PRESS (input units)")
# Q2 = 1 - PRESS / null_model_ss. A block with nothing to predict has an
# undefined Q2 rather than a perfect one, and so does a column with no
# spread across the held-out cells. A column has nothing to predict when its
# observed values are all equal, in any units, as MCUVScaler reads one: its
# null sum of squares is then 0, or rounding where it was centred at a mean.
# A floor of epsqrt on that sum took every column in small units for one (#1070).
varies = np.fmax.reduce(X_arr) > np.fmin.reduce(X_arr) # one per column: reduce runs over axis 0; both skip NaN
null_ss = null_model_ss if varies.any() else np.nan
q2 = (1.0 - press / null_ss).rename("Q2")
column_null_ss = np.where(varies, per_column_null, np.nan)
q2_per_variable = pd.DataFrame(
1.0 - raw_per_column / column_null_ss, index=component_index, columns=list(X.columns)
)
# PRESS ratio: computable under every scheme, kept for inspection.
# A block with no variation anywhere gives PRESS 0 at every component
# count, and 0/0 should reach the caller as a NaN ratio rather than as a
# numpy warning pointing into this function. The standard error below
# already guards the same case.
with np.errstate(divide="ignore", invalid="ignore"):
ratio_values = {a: press[a] / press[a - 1] for a in range(2, max_components + 1)}
press_ratio = pd.Series(ratio_values, name="PRESS ratio")
press_ratio.index.name = "n_components"
# Per-component standard error across the folds that evaluated the count;
# one fold has no spread to report.
per_fold_arr = per_fold_press.to_numpy()
n_folds_per_a = np.sum(~np.isnan(per_fold_arr), axis=1)
se_values = np.full(max_components, np.nan)
spread = n_folds_per_a > 1
se_values[spread] = np.nanstd(per_fold_arr[spread], axis=1, ddof=1) / np.sqrt(n_folds_per_a[spread])
# Rescale the per-fold-mean SE onto the same scale as ``press`` (see
# press_scale_multiplier above). Previously the ekf band was ~n_folds
# times too narrow, silently degenerating selection_rule="1se" to
# "min".
se_values = se_values * press_scale_multiplier
se_press = pd.Series(se_values, index=component_index, name="SE(PRESS)")
# The same constant null-model sum-of-squares that normalises Q2
# (Q2 = 1 - PRESS / null_model_ss) also rescales the PRESS standard
# error onto the Q2 scale, giving a directly usable +/-1 SE band around
# the Q2 curve without callers having to re-derive the normalisation.
q2_se = pd.Series(se_values / null_ss, index=component_index, name="SE(Q2)")
press_arr_final = press.to_numpy()
if np.all(np.isnan(press_arr_final)):
raise RuntimeError(
"Cross-validation produced NaN PRESS for every component count; no recommendation can be made."
)
recommended = _select_n_components(
selection_rule,
mean_error=press_arr_final,
se_error=se_values,
q2_cumulative=q2.to_numpy(),
min_q2_increase=min_q2_increase,
)
# A criterion that never turns over has not found an optimum: it has run
# out of components to evaluate. The count returned is then the largest
# one tried rather than an answer, and saying so is the difference
# between a recommendation and a number. The leverage approximations
# can fall into it on data whose R2 approaches one, where the residual
# they inflate has almost nothing left in it.
evaluated = q2.to_numpy()[~np.isnan(q2.to_numpy())]
if evaluated.size > 1 and np.all(np.diff(evaluated) > 0) and int(recommended) == max_components:
# Name the schemes that hold data out, less whichever is running:
# pointing a caller back at the scheme being warned about is no
# remedy. One of the two always survives, since a scheme cannot be
# both.
alternatives = " or ".join(f"cv_scheme={name!r}" for name in ("ekf", "ek") if name != cv_scheme)
warnings.warn(
f"cv_scheme={cv_scheme!r} did not turn over: its Q2 rises at every one of the "
f"{evaluated.size} component counts it could evaluate, so {recommended} is the "
f"largest count tried rather than an optimum. Evaluate more components, or use "
f"a scheme that holds data out ({alternatives}), before reading this as an answer.",
SpecificationWarning,
stacklevel=user_stacklevel(),
)
consensus_fields: dict[str, object] = {}
if return_consensus:
# Two cheap cross-checks: Minka's PPCA MLE (mean-centred input,
# no unit-variance scaling) and Horn's parallel analysis. The latter
# always autoscales: its null is unit-variance noise, and rescaling a
# block that is already scaled changes nothing. ``scale_inside_folds``
# is about the CV folds and must not reach it, or a raw block passed
# with False is compared in its own units against that null (#513).
minka_n = cls.minka_mle(X)
pa_result = cls.parallel_analysis(X, random_state=random_state)
counts = (int(recommended), int(minka_n), int(pa_result.n_components))
consensus = "agree" if max(counts) - min(counts) <= 1 else "disagree"
consensus_fields = {
"minka_n_components": int(minka_n),
"parallel_analysis_n_components": int(pa_result.n_components),
"consensus": consensus,
"consensus_counts": counts,
}
return Bunch(
n_components=recommended,
press=press,
per_fold_press=per_fold_press,
se_press=se_press,
q2_se=q2_se,
press_ratio=press_ratio,
q2=q2,
q2_per_variable=q2_per_variable,
press_input_units=press_input_units,
cv_scores=per_fold_press,
cv_scheme=cv_scheme,
selection_rule=selection_rule,
**consensus_fields,
)
[docs]
def detect_outliers(self, conf_level: float = 0.95, *, judged_rows: JudgedRows | None = None) -> list[dict]:
"""Detect outlier observations using SPE and Hotelling's T² diagnostics.
Combines two approaches:
1. **Statistical limits** - observations exceeding the SPE or T² limit
at ``conf_level`` are flagged.
2. **Robust ESD test** - the generalized ESD test identifies
observations that are unusual *relative to the rest of the data*,
even if they fall below the statistical limit. The mean/std variant
is used here; the underlying ``detect_outliers_esd`` also offers an
opt-in robust median/MAD variant.
An observation can be flagged for one or both reasons.
Parameters
----------
conf_level : float, default 0.95
Confidence level in [0.8, 0.999]. Controls both the statistical
limits and the ESD test's significance level (alpha = 1 - conf_level).
judged_rows : {"fitted", "new"} or None, optional
Which rows the T² limit judges (see :func:`hotellings_t2_limit`). This
method judges the rows the model was fitted to, so ``"fitted"`` gives
their limit. ``None``, the default, gives the ``"new"`` limit until 3.0
and the ``"fitted"`` one from 3.0, with a :class:`FutureWarning`.
Returns
-------
outliers : list of dict
Sorted from most severe to least. Each dict contains:
- ``observation`` - index label of the observation
- ``outlier_types`` - list of ``"spe"`` and/or ``"hotellings_t2"``
- ``spe`` - SPE value for this observation
- ``hotellings_t2`` - T² value for this observation
- ``spe_limit`` - SPE limit at the given confidence level
- ``hotellings_t2_limit`` - T² limit at the given confidence level
- ``severity`` - max(spe/spe_limit, t2/t2_limit), rounded to 4
decimals. A ratio whose denominator is 0 (perfect-fit SPE limit)
or non-finite (T2 limit when A == N) is treated as 0 and does
not contribute to the ranking.
- ``hotellings_t2_limit_note`` - only when ``judged_rows`` is
``None`` and the limit changes in 3.0: the text of the warning.
Examples
--------
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> outliers = pca.detect_outliers(conf_level=0.95, judged_rows="fitted")
>>> for o in outliers:
... print(f"{o['observation']}: {o['outlier_types']} (severity={o['severity']})")
"""
check_is_fitted(self, "spe_")
if not (0.8 <= conf_level <= 0.999):
raise ValueError(f"conf_level must be between 0.8 and 0.999, got {conf_level}.")
N = self.n_samples_
# Full-model SPE and cumulative T² (last column)
spe_values = self.spe_.iloc[:, -1]
t2_values = self.hotellings_t2_.iloc[:, -1]
# Statistical limits
spe_lim = self.spe_limit(conf_level=conf_level)
t2_lim = self.hotellings_t2_limit(conf_level=conf_level, judged_rows=judged_rows)
# What the warning said, for a caller that reads the records rather than the warnings.
note = judged_rows_note(conf_level, self.n_components_, N, judged_rows)
# Robust ESD outlier detection on each series
max_outliers = max(1, N // 5)
alpha = 1 - conf_level
spe_outlier_idx, _ = detect_outliers_esd(spe_values.to_numpy(), max_outliers_detected=max_outliers, alpha=alpha)
t2_outlier_idx, _ = detect_outliers_esd(t2_values.to_numpy(), max_outliers_detected=max_outliers, alpha=alpha)
# Collect all flagged observations: ESD outliers, and any observation
# above the statistical limit
spe_flagged = set(spe_outlier_idx) | {int(i) for i in np.flatnonzero(spe_values.to_numpy() > spe_lim)}
above_t2_limit = t2_values.to_numpy() > t2_lim # pragma: no mutate - > and >= tie only at the limit itself
t2_flagged = set(t2_outlier_idx) | {int(i) for i in np.flatnonzero(above_t2_limit)}
# Merge into result dicts
all_flagged = spe_flagged | t2_flagged
results = []
for i in all_flagged:
types = []
if i in spe_flagged:
types.append("spe")
if i in t2_flagged:
types.append("hotellings_t2")
spe_val = float(spe_values.iloc[i])
t2_val = float(t2_values.iloc[i])
# A degenerate limit carries no information about how severe an
# observation is, so it must not contribute to the ranking. A
# perfect fit gives spe_lim == 0, which used to raise
# ZeroDivisionError outright, and a limit at machine-epsilon scale
# produced impressive-looking severities that were ratios of
# floating-point noise. The ratio itself is scale-invariant (value
# and limit share units), so the guard is only against a limit of
# zero. Deliberately NOT an absolute epsilon, which would make
# severity depend on the units of the data. The T2 limit is
# positive, and infinite when A == N, where its ratio is 0 by itself.
severity = max(spe_val / spe_lim if spe_lim > 0.0 else 0.0, t2_val / t2_lim)
results.append(
{
"observation": spe_values.index[i],
"outlier_types": types,
"spe": spe_val,
"hotellings_t2": t2_val,
"spe_limit": spe_lim,
"hotellings_t2_limit": t2_lim,
"severity": round(severity, 4),
}
| ({} if note is None else {"hotellings_t2_limit_note": note})
)
results.sort(key=lambda d: d["severity"], reverse=True)
return results
