perf: compute only top-rank principal directions in PCAEM

[jamgochiana]

Closes #18

Problem

PCAEM.fit called pca_fit(PCA, x; maxoutdim=rank). maxoutdim only truncates an already-computed full decomposition; it does not make MultivariateStats compute fewer components. For d < n PCA does a full eigen of the d x d covariance, and for d >= n a full svd of the centered data. When rank is far smaller than min(d, n), almost all of that work is wasted.

Approach

PCAEM.fit only needs two things from PCA: the loadings P (d x rank) and the mean mu. Everything downstream (reduced data, reconstruction, and the diagonal noise D) derives from those.

The new _truncated_pca helper eigendecomposes only whichever Gram matrix is smaller:

• d <= n branch: eigendecompose the d x d covariance Symmetric((Z Z') / (n-1)) and take the top-rank eigenvectors as P.
• d > n branch: eigendecompose the smaller n x n Gram matrix Symmetric((Z' Z) / (n-1)), then recover feature-space loadings via U = Z W and normalize columns.

Choosing the smaller matrix by the d <= n branch is essential: eigendecomposing the larger Gram matrix is catastrophically slow (e.g. eigs(ZZ') on d=20000, n=1000 is ~12x slower than baseline).

Small-dimension dense fallback

When the smaller dimension is small, Arpack's iterative eigs has high overhead (and can fail to converge), and eigs requires nev < size - 1. The helper falls back to a dense eigen + truncation when min(d, n) <= 100 or rank >= min(d, n) - 1. The existing PCAEM tests (d=20, n=1000) and edge cases hit this fallback, which is correct and expected. The existing rank <= n_features ArgumentError check is preserved.

New dependency

Adds Arpack (eigs) to [deps] and [compat].

Equivalence tests

Added a "Truncated PCA equivalence" testset in test/test_fitting.jl verifying, against pca_fit(PCA, x; maxoutdim=rank), that:

• The mean mu matches (~1e-10).
• The reconstruction subspace projector P*P' matches P_full*P_full' (sign/rotation-invariant) to ~1e-8.
• The reconstruction-error variance diagonal D matches to ~1e-8.

These cover both branches and both the dense-fallback path (d=20/n=1000, d=200/n=30) and the iterative eigs path (d=150/n=400, d=400/n=150). Loadings are only defined up to per-column sign / rotation, so only sign/rotation-invariant quantities are asserted. The existing PCAEM testset (LRDMvNormal components, correct rank, finite logpdf, valid weights) and edge-case tests remain unchanged and pass.

Results

Full suite passes: LRDMvNormal 53/53, Fitting Methods 104/104, Prediction and Marginal Functions 24/24. (Unrelated GaussianMixtures "Too low occupancy" warnings are expected.) Sanity benchmark on d=5000, n=800, rank=5 showed ~9x speedup vs the full PCA.

🤖 Generated with Claude Code