Multivariate Analysis

Principal component analysis and maximum likelihood factor analysis with varimax and promax rotations.

Multivariate analysis module.

Provides principal component analysis and factor analysis, matching R’s implementations.

Public API:

pca(X) - Principal Component Analysis (matches prcomp) factor_analysis(X) - Maximum Likelihood Factor Analysis (matches factanal)

pystatistics.multivariate.pca(X, *, center=True, scale=False, n_components=None, names=None, backend=None, solver=None, force=False, device_resident=False, oversample=10, n_iter=4, seed=0)[source]

Principal Component Analysis via SVD.

Matches R’s stats::prcomp().

Algorithm:

  1. Center X (subtract column means).

  2. Optionally scale (divide by column standard deviations).

  3. Thin SVD: X_centered = U @ diag(S) @ V’.

  4. sdev = S / sqrt(n - 1).

  5. rotation = V (right singular vectors = loadings).

  6. x = U @ diag(S) (scores = X_centered @ V).

Parameters:

  • X (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Data matrix (n x p), n observations, p variables.

  • center (bool) – Whether to center columns (subtract mean). Default True.

  • scale (bool) – Whether to scale columns (divide by SD). Default False. Equivalent to R’s prcomp(scale. = TRUE).

  • n_components (int | None) – Number of components to retain. Default: min(n, p).

  • names (list[str] | None) – Variable names for the p columns.

  • backend (str | None) –

    Compute backend. The default depends on the input type:

    • numpy array / CPU DataSource → 'cpu' (R-reference path; matches R prcomp() to rtol = 1e-10). This is the regulated-industry default: unspecified means the validated-against-R path.

    • torch.Tensor on GPU (e.g. ds.to('cuda')['X']) → 'gpu'. Creating a GPU DataSource is already an opt-in to the GPU path, so we don’t demand the user repeat the choice with an explicit backend='gpu'.

    Explicit values:

    • 'cpu': force the numpy reference path (float64). Raises if given a GPU tensor (no silent device migration — Rule 1).

    • 'gpu': force GPU, float32. Raises if no GPU is available. Results match the CPU path at the GPU_FP32 tolerance tier (rtol ≈ 1e-4) — statistically equivalent, not bitwise-equivalent.

    • 'gpu_fp64': force GPU, float64. CUDA only (MPS lacks float64; raises on Apple Silicon). CPU-matching precision at a performance cost — consumer NVIDIA parts run FP64 at ~1/64× FP32, so this is usually slower than CPU LAPACK.

    • 'auto': prefer GPU (float32) when CUDA is available, else CPU.

  • solver (str | None) –

    Numerical algorithm to use. Only matters when actually running on GPU (the CPU path always uses SVD). The default (None) resolves by device: 'svd' on CUDA, and 'randomized' on Apple Silicon (MPS), where it is the only genuinely on-device path. 'svd': SVD of X. Always safe. Moderate (~3–4×) speedup

    over CPU LAPACK. CUDA only — on MPS torch.linalg.svd silently falls back to the CPU, so an explicit solver='svd' on Apple Silicon raises (no silent CPU substitution).

    'gram': Eigendecomposition of X’X — turns PCA into a

    big GEMM + small symmetric eigendecomp, both GPU sweet spots. For tall-skinny well-conditioned data (n ≫ p) this is typically 30–100×+ faster than the SVD path. Squares the condition number, so raises NumericalError on ill-conditioned data unless force=True. Precision gates: cond(X) ≲ 1e6 for FP64, cond(X) ≲ 1e3 for FP32. CUDA only — on MPS torch.linalg.eigh has no Metal kernel, so an explicit solver='gram' on Apple Silicon raises.

    'randomized': Halko–Martinsson–Tropp randomized truncated

    SVD with CholeskyQR2 orthonormalization. Every X-sized operation runs on-device via matmul + cholesky (only a tiny l×l eigendecomposition touches the host), so this is the path that makes PCA run — and beat the CPU 3–8× — on Apple Silicon. Device-agnostic, also selectable on CUDA. Approximate (top-k only) and seeded; see oversample / n_iter / seed. Squares the condition number inside CholeskyQR, so it gates on cond(X) ≲ 1e3 (FP32) just like the Gram path; force=True bypasses.

    'auto': On CUDA, uses 'gram' when n > 2p AND the

    condition check passes; falls back to 'svd' otherwise. On MPS, routes to 'randomized'.

  • force (bool) – Bypass the condition-number gate (Gram and randomized paths). Numerical results will be unreliable on truly ill-conditioned inputs — use only when you understand the data is well-conditioned despite the automated estimator disagreeing.

  • oversample (int) – Randomized solver only — extra sketch dimensions drawn beyond n_components for accuracy (default 10).

  • n_iter (int) – Randomized solver only — number of power iterations (default 4). Higher is more accurate on slowly-decaying spectra at proportional cost.

  • seed (int) – Randomized solver only — RNG seed for the Gaussian sketch. The default is fixed, so a given input is reproducible; pass a different value to draw an independent sketch. GPU RNG is not bit-identical to the CPU RNG (statistically equivalent, not bitwise-equal).

  • device_resident (bool) – When True and the GPU backend runs, the returned PCASolution holds its numeric fields (sdev, rotation, center, scale, x) as torch.Tensor instances on the fit’s device rather than materialising them as numpy arrays. This saves the ~150 ms D2H copy of the scores matrix on 1M × 100 FP32 data, which otherwise dominates any multi-step GPU pipeline that consumes PCA output. Call result.to_numpy() or result.to('cpu') to materialise a numpy-backed copy. Ignored on the CPU path (result is always numpy-backed there). Default False preserves 1.8.0 behavior.

Returns:

PCASolution with sdev, rotation (loadings), scores, etc.

Raises:

ValidationError – If inputs are invalid.

Return type:

PCASolution

Validates against: R stats::prcomp().

pystatistics.multivariate.factor_analysis(X, *, n_factors, rotation='varimax', method='ml', names=None, lower=0.005, tol=1e-08, max_iter=1000)[source]

Maximum likelihood factor analysis.

Matches R’s stats::factanal().

The model is: Sigma = Lambda Lambda’ + Psi, where Lambda is the loadings matrix (p x m) and Psi = diag(psi_1, …, psi_p) is the uniqueness matrix.

Algorithm:

  1. Minimise the ML objective over uniquenesses psi using L-BFGS-B, parameterised in log(psi) to ensure psi > 0.

  2. Extract loadings from the eigendecomposition of Psi^{-1/2} S Psi^{-1/2}.

  3. Apply rotation (varimax, promax, or none).

  4. Compute chi-squared goodness-of-fit test.

Parameters:

  • X (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Data matrix (n x p).

  • n_factors (int) – Number of factors to extract.

  • rotation (str) – ‘varimax’ (orthogonal), ‘promax’ (oblique), or ‘none’.

  • method (str) – ‘ml’ (only ML supported).

  • names (list[str] | None) – Variable names.

  • lower (float) – Lower bound on the uniquenesses, enforced during the optimisation. Matches R’s factanal(lower=) (default 0.005). Guards against degenerate Heywood solutions where a uniqueness collapses to ~0. Must satisfy 0 < lower < 1.

  • tol (float) – Convergence tolerance.

  • max_iter (int) – Maximum iterations.

Returns:

FactorSolution with loadings, uniquenesses, test statistics, etc.

Raises:
Return type:

FactorSolution

Validates against: R stats::factanal().

class pystatistics.multivariate.PCASolution(_result)[source]

Bases: SolutionReprMixin

Result from principal component analysis.

Matches the structure of R’s prcomp() return value. Wraps a Result [PCAParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Device residency

The numeric fields (sdev, rotation, center, scale, x) may be NDArray or torch.Tensor. When the GPU backend is invoked with a device-resident input tensor and device_resident=True is passed to pca(), the result’s fields stay on the original device – the x scores matrix especially, which is typically the largest payload. Downstream computations can chain straight off result.x without paying the host<->device round-trip that otherwise dominates multi-step GPU pipelines. Call to_numpy() (or to() with a CPU device) to materialise a numpy-backed copy.

property sdev: ndarray[tuple[Any, ...], dtype[_ScalarT]]
property rotation: ndarray[tuple[Any, ...], dtype[_ScalarT]]
property center: ndarray[tuple[Any, ...], dtype[_ScalarT]]
property scale: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None
property x: ndarray[tuple[Any, ...], dtype[_ScalarT]]
property n_obs: int
property n_vars: int
property var_names: tuple[str, ...] | None
property info: dict
property timing: dict[str, float] | None
property backend_name: str
property warnings: tuple[str, ...]
property device: str

'cpu', 'cuda', or 'mps'.

A numpy-backed result reports 'cpu'. A tensor-backed result reports the underlying tensor’s device type.

Type:

Device hosting the numeric fields

to_numpy()[source]

Return a new PCASolution with numpy-backed numeric fields.

Idempotent on numpy-backed results (returns self). On a tensor-backed result this is the explicit “I want CPU numpy” escape hatch, matching the rest of the library’s no-silent-migration contract.

Return type:

PCASolution

to(device)[source]

Move the result to device ('cpu', 'cuda', 'mps').

device='cpu' returns a numpy-backed PCASolution (same as to_numpy()). Any other device requires torch and materialises the fields as torch.Tensor instances on that device.

Parameters:

device (str)

Return type:

PCASolution

property explained_variance_ratio

Proportion of variance explained by each component.

Returned as the same array type as sdev – numpy for a numpy-backed result, torch.Tensor on the same device for a tensor-backed result. Both numpy ndarrays and torch tensors support .sum() and elementwise arithmetic so the impl is type-polymorphic.

property cumulative_variance_ratio

Cumulative proportion of variance explained.

summary()[source]

R-style summary matching summary(prcomp(...)).

Returns:

Formatted string with standard deviations, proportion of variance, and cumulative proportion for each component.

Return type:

str

Parameters:

_result (Result[PCAParams])

class pystatistics.multivariate.FactorSolution(_result)[source]

Bases: SolutionReprMixin

Result from factor analysis.

Matches the structure of R’s factanal() return value. Wraps a Result [FactorParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Parameters:

_result (Result[FactorParams])

property loadings: ndarray[tuple[Any, ...], dtype[_ScalarT]]
property uniquenesses: ndarray[tuple[Any, ...], dtype[_ScalarT]]
property communalities: ndarray[tuple[Any, ...], dtype[_ScalarT]]
property rotation_matrix: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None
property chi_sq: float | None
property p_value: float | None
property dof: int
property n_factors: int
property n_obs: int
property n_vars: int
property var_names: tuple[str, ...] | None
property method: str
property rotation_method: str
property converged: bool
property n_iter: int
property objective: float
property info: dict
property timing: dict[str, float] | None
property backend_name: str
property warnings: tuple[str, ...]
summary()[source]

R-style summary matching print(factanal(...)).

Returns:

Formatted string with loadings, uniquenesses, SS loadings, proportion and cumulative variance, and the chi-squared test.

Return type:

str

class pystatistics.multivariate.PCAParams(sdev, rotation, center, scale, x, n_obs, n_vars, var_names)[source]

Bases: object

Immutable parameter payload for principal component analysis.

The numeric fields may be numpy arrays (CPU path) or torch.Tensor instances (device-resident GPU path); see PCASolution for the device-residency contract.

Parameters:
sdev

Standard deviations of principal components (length min(n, p)).

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]]

rotation

Loadings matrix (p x n_components) – columns are eigenvectors.

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]]

center

Column means used for centering (length p).

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]]

scale

Column SDs used for scaling (None if scale=False).

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]] | None

x

Scores matrix (n x n_components).

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]]

n_obs

Number of observations.

Type:

int

n_vars

Number of variables.

Type:

int

var_names

Variable names, or None.

Type:

tuple[str, …] | None

sdev: ndarray[tuple[Any, ...], dtype[_ScalarT]]
rotation: ndarray[tuple[Any, ...], dtype[_ScalarT]]
center: ndarray[tuple[Any, ...], dtype[_ScalarT]]
scale: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None
x: ndarray[tuple[Any, ...], dtype[_ScalarT]]
n_obs: int
n_vars: int
var_names: tuple[str, ...] | None
class pystatistics.multivariate.FactorParams(loadings, uniquenesses, communalities, rotation_matrix, chi_sq, p_value, dof, n_factors, n_obs, n_vars, var_names, method, rotation_method, converged, n_iter, objective)[source]

Bases: object

Immutable parameter payload for factor analysis.

Parameters:
loadings

Rotated loadings matrix (p x n_factors).

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]]

uniquenesses

Uniqueness for each variable (length p).

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]]

communalities

1 - uniquenesses.

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]]

rotation_matrix

Rotation matrix, or None if no rotation.

Type:

numpy.ndarray[tuple[Any, …], numpy.dtype[numpy._typing._array_like._ScalarT]] | None

chi_sq

Goodness-of-fit chi-squared statistic, or None.

Type:

float | None

p_value

p-value for chi-sq test, or None.

Type:

float | None

dof

Degrees of freedom.

Type:

int

n_factors

Number of factors extracted.

Type:

int

n_obs

Number of observations.

Type:

int

n_vars

Number of variables.

Type:

int

var_names

Variable names, or None.

Type:

tuple[str, …] | None

method

Estimation method (e.g. ‘ml’).

Type:

str

rotation_method

Rotation method (e.g. ‘varimax’, ‘promax’, ‘none’).

Type:

str

converged

Whether the optimisation converged.

Type:

bool

n_iter

Number of iterations used.

Type:

int

objective

Final objective function value.

Type:

float

loadings: ndarray[tuple[Any, ...], dtype[_ScalarT]]
uniquenesses: ndarray[tuple[Any, ...], dtype[_ScalarT]]
communalities: ndarray[tuple[Any, ...], dtype[_ScalarT]]
rotation_matrix: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None
chi_sq: float | None
p_value: float | None
dof: int
n_factors: int
n_obs: int
n_vars: int
var_names: tuple[str, ...] | None
method: str
rotation_method: str
converged: bool
n_iter: int
objective: float