2.10. Principal Component Analysis

Principal Component Analysis (PCA) finds an orthonormal set of directions in feature space. These directions are the principal axes. The first axis captures the most variance in the data. The second axis is orthogonal to the first and captures the most of the remaining variance. Each later axis follows the same pattern.

Projecting the data onto the top k axes gives the k-dimensional linear subspace that keeps the most variance. This is the same subspace that gives the smallest squared reconstruction error. The subspace that keeps the most variance also discards the least information. This dual property is why PCA is a common first step to compress, denoise, or plot high-dimensional data.

RustyML’s PCA mirrors sklearn.decomposition.PCA. Fit it on a feature matrix. Use transform to get scores. Use inverse_transform to map scores back to feature space.

RustyML’s PCA differs from scikit-learn’s PCA in 3 ways. First, it centers the data but never scales it. The values you feed it matter more than in a pipeline that standardizes for you.

Second, the solver setting is an enum named SVDSolver, with 3 concrete strategies. The right choice depends on the number of features and the number of components you want.

Third, every factorization runs on hand-written pure Rust. There is no LAPACK and no BLAS. This makes the tradeoffs between solvers concrete rather than theoretical.

2.10.1. The surface you work with

PCA and SVDSolver live under rustyml::machine_learning::decomposition. Both types are also available through the prelude. The model is unsupervised. It takes an f64 feature matrix with samples as rows and features as columns. It uses no labels.

// Construction. n_components must be > 0, checked here.
PCA::new(n_components: usize) -> Result<PCA, Error>
PCA::default()                                         // n_components = 2, Full solver
pca.with_svd_solver(solver: SVDSolver) -> PCA          // builder, consumes and returns self

// Learning and projecting. fit takes &mut self.
pca.fit(&x)                -> Result<&mut PCA, Error>
pca.transform(&x)          -> Result<Array2<f64>, Error>   // (n_samples, n_components)
pca.fit_transform(&x)      -> Result<Array2<f64>, Error>
pca.inverse_transform(&scores) -> Result<Array2<f64>, Error> // (n_samples, n_features)

// Fitted state. mean, components, variances, and the sample/feature counts are None before
// fit, Some(&...) after. n_components and svd_solver are always available.
pca.get_mean()                     -> Option<&Array1<f64>>   // per-feature centering mean
pca.get_components()               -> Option<&Array2<f64>>   // (n_components, n_features), rows are axes
pca.get_explained_variance()       -> Option<&Array1<f64>>
pca.get_explained_variance_ratio() -> Option<&Array1<f64>>
pca.get_singular_values()          -> Option<&Array1<f64>>
pca.get_n_components() -> usize
pca.get_svd_solver()   -> SVDSolver
pca.get_n_samples()    -> Option<usize>
pca.get_n_features()   -> Option<usize>

The constructor checks only that n_components > 0. The harder bound, n_components <= min(n_samples, n_features), is checked at fit time, because it depends on the data. Use the getters to check whether the model is fitted. Each one returns None until a successful fit call sets it.

get_components returns the axes as rows of an (n_components, n_features) matrix. This is the same layout as scikit-learn’s components_. So transform is centered data times components.T, and inverse_transform is scores times components.

This example runs the full loop on a small 2-D data set and keeps 1 component:

use ndarray::array;
use rustyml::machine_learning::decomposition::PCA;

fn main() {
    let x = array![
        [2.5, 2.4],
        [0.5, 0.7],
        [2.2, 2.9],
        [1.9, 2.2],
        [3.1, 3.0],
        [2.3, 2.7],
    ];

    // fit_transform needs &mut self. It runs fit, then transform, on the same data.
    let mut pca = PCA::new(1).unwrap();
    let scores = pca.fit_transform(&x).unwrap();

    println!("scores shape: {:?}", scores.shape()); // [6, 1]
    println!("mean:         {:?}", pca.get_mean().unwrap());
    println!("component:    {:?}", pca.get_components().unwrap());
    println!("variance:     {:?}", pca.get_explained_variance().unwrap());
    println!("ratio:        {:?}", pca.get_explained_variance_ratio().unwrap());
    println!("singular:     {:?}", pca.get_singular_values().unwrap());
}

fit_transform(&x) is not a fused shortcut that gives a different result. It calls fit, then transform, on the same matrix. This is identical to running the two steps by hand. Use the two-step form to project a different matrix, for example new samples, through an already-fitted model. fit returns &mut Self, so you can chain a getter call directly onto it. The normal path, though, is to read state back with the get_* methods.

2.10.2. It centers, it does not scale

Before it decomposes the data, fit subtracts the per-feature mean and stores it in get_mean(). This is the only preprocessing step fit runs. It does not divide by the standard deviation. It does not whiten the data. It does not touch the units of the columns.

transform reuses the stored training mean to center new samples. inverse_transform adds the mean back. The mean is part of the fitted model, not a value computed fresh on each call.

Skipping this step is the most common cause of a misleading PCA result. Variance is not scale-invariant. Converting a feature from meters to millimeters multiplies its values by 1000, so its variance grows by 1000^2, a million times. PCA then assigns its first component to that one column.

The fix is to standardize the features before fitting. Standardizing puts every column on unit variance, so PCA compares them on equal footing. Use standardize with StandardizationAxis::Column. The next example shows the difference:

use ndarray::array;
use rustyml::machine_learning::decomposition::PCA;
use rustyml::utils::standardize::{standardize, StandardizationAxis};

fn main() {
    // Feature 0 uses tiny units. Feature 1 uses huge units. They carry the same
    // information, but very different variance.
    let x = array![
        [0.1, 1000.0],
        [0.2, 1005.0],
        [0.3, 1002.0],
        [0.4,  998.0],
        [0.5, 1010.0],
    ];

    let mut raw = PCA::new(2).unwrap();
    raw.fit(&x).unwrap();
    println!("raw ratios:          {:?}", raw.get_explained_variance_ratio().unwrap());

    let xs = standardize(&x, StandardizationAxis::Column).unwrap();
    let mut scaled = PCA::new(2).unwrap();
    scaled.fit(&xs).unwrap();
    println!("standardized ratios: {:?}", scaled.get_explained_variance_ratio().unwrap());
}

On the raw data, the first ratio is close to 1.0. Column 1’s numeric spread swamps the signal from column 0, so PC1 points almost exactly along column 1. After column standardization, the two features contribute about equally. The ratios then reflect the real correlation structure between them.

Whether to standardize is a modeling decision. If the features already share the same meaningful unit, for example pixel intensities, centering alone is enough. Make this decision on purpose, because rustyml does not make it for you. Statisticians describe this choice as PCA on the covariance matrix versus PCA on the correlation matrix. Standardizing first is the same as running PCA on the correlation matrix.

2.10.3. Choosing a solver, and the pure-Rust machinery behind it

SVDSolver selects how PCA computes the components. The name is not fully accurate. Only 1 of the 3 variants literally forms an SVD. The interface stays the same across all 3 variants. Whichever solver you choose, you get the same component rows, explained variances, and singular values. Each solver computes them through a different route, with a different cost and a different accuracy.

pub enum SVDSolver {
    Full,             // default: exact eigendecomposition of the covariance matrix
    Randomized(u64),  // randomized range finder, seeded by the u64
    PowerIteration,   // deflated power iteration for the top-k eigenpairs
}

The cost of each solver follows a clear pattern. Full and PowerIteration both build the dense d x d covariance. Both pay O(n * d^2) to form it and O(d^2) memory to hold it. The difference is what happens next. Full runs a complete O(d^3) eigendecomposition.

PowerIteration extracts k eigenpairs one at a time, using a bounded number of matrix-vector products for each, and deflates each one before finding the next. This wins on compute time when k is far smaller than d.

Randomized never builds a d x d matrix at all. Its working set is the n x k sketch and the k x d reduced projection. This makes it the memory-frugal choice when d is large. On small problems, all 3 solvers agree to several digits. The approximate solvers only show their speed advantage, or their numerical drift, once the matrices grow large.

use ndarray::array;
use rustyml::machine_learning::decomposition::{PCA, SVDSolver};

fn main() {
    let x = array![
        [2.5, 2.4],
        [0.5, 0.7],
        [2.2, 2.9],
        [1.9, 2.2],
        [3.1, 3.0],
        [2.3, 2.7],
        [2.0, 1.6],
        [1.0, 1.1],
    ];

    for solver in [SVDSolver::Full, SVDSolver::Randomized(42), SVDSolver::PowerIteration] {
        let mut pca = PCA::new(1).unwrap().with_svd_solver(solver);
        pca.fit(&x).unwrap();
        let sigma = pca.get_singular_values().unwrap();
        println!("{:?}: sigma_1 = {:.6}", solver, sigma[0]);
    }
}

PowerIteration and Randomized expose only a small part of their internal settings. PowerIteration has no public setting for iteration count or tolerance. Internally, it runs up to 1000 iterations per component, to a 1e-6 eigenvalue tolerance, with a fixed seed. None of these values is tunable through the public API.

PowerIteration can also fail to converge. This can happen when the data has less effective rank than the number of components you request. For example, duplicate columns or an exact linear dependency between features can cause this. When the deflation step runs out of variance to extract, fit returns Error::NotConverged. Full and Randomized never fail this way, because neither depends on a per-component convergence check.

Randomized exposes only its seed. The oversampling amount and the number of subspace iterations are fixed. So the full configuration surface for PCA is n_components plus with_svd_solver. This is small on purpose.

All 3 solvers run on the crate’s own machine_learning::linalg module. This module reimplements symmetric eigendecomposition, SVD (one-sided Jacobi), and QR (modified Gram-Schmidt) directly on ndarray arrays. There is no LAPACK dependency to install, link, or vendor. This matches the pure-Rust numerics approach the rest of the crate takes. The covariance and projection multiplies use the parallel gemmkit backend.

2.10.4. Reading the variance and choosing n_components

After a fit, 3 vectors describe how the variance spreads across the kept axes. get_singular_values returns the singular values, in descending order, of the centered data. get_explained_variance returns each singular value squared, divided by (n - 1). This is the actual variance along each component, in the units of the original data.

get_explained_variance_ratio divides each of those values by the total variance of the centered data. That total is the full trace, summed over all features, not just the ones you kept. So each entry in the ratio is the fraction of the overall variance that component explains.

This denominator choice matters. If you keep fewer components than features, the ratios sum to less than 1. The shortfall is exactly the variance you discarded.

This property makes the ratio the right tool for choosing n_components. The standard method uses a scree plot and a cumulative rule. Fit the full-rank decomposition once. Read the cumulative ratio. Keep the smallest number of axes that clears a variance target, for example 90%, 95%, or 99%.

use ndarray::array;
use rustyml::machine_learning::decomposition::PCA;

fn main() {
    // 8 samples, 3 features. Feature 2 is almost a linear echo of feature 0,
    // so the data is effectively low-rank.
    let x = array![
        [1.0, 0.2, 1.19],
        [2.0, 0.1, 2.09],
        [3.0, 0.4, 3.38],
        [4.0, 0.2, 4.19],
        [5.0, 0.5, 5.48],
        [6.0, 0.3, 6.29],
        [7.0, 0.1, 7.09],
        [8.0, 0.6, 8.58],
    ];

    // Fit the full decomposition (n_components = min(n_samples, n_features) = 3),
    // then decide how many axes to actually keep.
    let mut pca = PCA::new(3).unwrap();
    pca.fit(&x).unwrap();

    let ratio = pca.get_explained_variance_ratio().unwrap();
    let mut cumulative = 0.0;
    for (i, r) in ratio.iter().enumerate() {
        cumulative += r;
        println!("PC{}: individual = {:.4}, cumulative = {:.4}", i + 1, r, cumulative);
    }

    // Smallest k whose cumulative ratio clears 95%.
    let mut acc = 0.0;
    let mut k = ratio.len();
    for (i, r) in ratio.iter().enumerate() {
        acc += r;
        if acc >= 0.95 {
            k = i + 1;
            break;
        }
    }
    println!("keep {} component(s) to retain >= 95% variance", k);
}

One detail matters here. You fit once at full rank to read the spectrum. Then refit with the chosen k if you want the smaller projection. Refitting is cheap, and it gives you a model whose transform outputs exactly k columns.

If refitting is expensive for the data, use a shortcut instead. The first k component rows of the full fit are the same axes as a k-component fit would produce. The only difference is the sign convention, covered in the next section. So you can slice the full fit rather than refit it. Do not pick k from a single ratio in isolation. Read the cumulative curve instead, to see where the returns flatten out.

2.10.5. Reconstruction and inverse_transform

inverse_transform maps scores back to feature space, using the formula reconstructed = scores * components + mean. When you keep all min(n_samples, n_features) components, the round-trip is lossless, up to floating-point error. When you keep fewer components, the reconstruction is the orthogonal projection of the data onto the retained subspace. The residual is the variance in the axes you dropped. Measuring that residual is a direct way to put a number on the cost of compression.

use ndarray::array;
use rustyml::machine_learning::decomposition::PCA;

fn main() {
    let x = array![
        [1.0, 0.2, 1.19],
        [2.0, 0.1, 2.09],
        [3.0, 0.4, 3.38],
        [4.0, 0.2, 4.19],
        [5.0, 0.5, 5.48],
        [6.0, 0.3, 6.29],
    ];

    // Keep a single axis, then round-trip back into the original 3-D feature space.
    let mut pca = PCA::new(1).unwrap();
    pca.fit(&x).unwrap();

    let scores = pca.transform(&x).unwrap();                     // (6, 1)
    let reconstructed = pca.inverse_transform(&scores).unwrap(); // (6, 3)

    // Frobenius norm of the residual = the variance thrown away with PC2 and PC3.
    let err: f64 = (&reconstructed - &x).iter().map(|d| d * d).sum::<f64>().sqrt();
    println!("scores shape:         {:?}", scores.shape());
    println!("reconstruction shape: {:?}", reconstructed.shape());
    println!("reconstruction error: {:.6}", err);
}

Watch the shape contract, because it reverses transform. transform takes n_features columns and returns n_components columns. inverse_transform takes n_components columns and returns n_features columns. Pass it a matrix whose column count does not match n_components, and you get Error::DimensionMismatch, not a silent broadcast.

The sign fix, covered in the next section, flips a whole axis together with its scores. This makes reconstruction invariant to it. The product scores * components stays the same, whether or not an axis was negated. A sign flip never corrupts a round-trip.

2.10.6. Determinism: sign convention and seeds

Eigenvectors and singular vectors are only defined up to sign. -v spans the same axis as v. scikit-learn users know this as the reason PC signs sometimes flip between runs or between library versions.

RustyML fixes the sign after the decomposition, so the result stays deterministic. Each component row is negated, if needed, so its largest-magnitude loading becomes non-negative. As a result, all 3 solvers agree on the orientation of every axis on the same data. Repeated runs of any one solver also agree. The sign stays stable, and fit is reproducible with no extra step.

This has 2 consequences. First, the sign convention is specific to rustyml. It keys off the component vectors themselves, not off the U factor the way scikit-learn’s svd_flip does. So a component’s sign may differ from what scikit-learn prints for the same data. This difference is cosmetic. The axis, the variance it explains, and every reconstruction stay identical.

Second, the only real source of run-to-run variation is SVDSolver::Randomized. Its random sketch is seeded by the u64 you pass. The same seed gives bit-identical output. A different seed gives a slightly different approximate subspace. Full and PowerIteration carry no external randomness. PowerIteration seeds its own starting vector internally, with a fixed value.

For exact reproducibility across a pipeline, see Reproducibility and Random Seeds. The only lever here is the Randomized seed.

A fitted PCA serializes with save_to_path and load_from_path. These methods write and read the whole model (mean, components, variances, and singular values) as compact postcard binary. The file extension does not matter. The bytes stay postcard format regardless of the file name. A reloaded model transforms data identically to the original model:

pca.save_to_path("pca_model.bin")?;
let loaded = PCA::load_from_path("pca_model.bin")?;
let scores = loaded.transform(&x_new)?; // identical to the pre-save model

For the mechanics, including versioning concerns and when the binary format is safe for long-term storage, see Model Persistence in Depth.

2.10.7. Errors and edge cases

PCA validates its input and returns typed errors instead of panicking. The .unwrap() calls in these examples keep the code short. Do not use .unwrap() this way in production code. The table below lists the errors you can hit in practice, plus the one that is specific to the PowerIteration solver:

The NotConverged case is specific to PowerIteration. Full and Randomized never raise it, because neither depends on a strict per-component convergence check.

The 2-sample minimum is not arbitrary. Variance needs an n - 1 denominator, so a single row has nothing to decompose. The n_components ceiling of min(n_samples, n_features) is the rank bound. You cannot extract more orthogonal directions than the data spans.

Asking for more raises a parameter error, instead of a silently truncated result. rustyml checks this bound against the data, not against the model alone. So the same PCA::new(5) instance can succeed on a 100 x 20 matrix and fail on a 3 x 4 one.

PCA is linear by construction. It can only find directions that are linear combinations of the input features. Structure that lives on a curved manifold stays invisible to it. When a scree plot does not flatten, and reconstruction stays poor at every k, this is usually a signal to move to a nonlinear method.

Kernel PCA applies the same variance-maximizing idea in a kernel feature space. t-SNE is the tool for 2-D visualization of nonlinear neighborhood structure. If the goal is to separate labeled classes, rather than to capture raw variance, use the supervised counterpart, Linear Discriminant Analysis. It projects the data toward class separation instead of toward total spread.