2.11. Kernel PCA

KernelPCA runs PCA in an implicit feature space. It does not decompose the covariance of your data. Instead, it decomposes a centered kernel (Gram) matrix. This finds nonlinear structure that a linear projection cannot see.

Kernel PCA works with pairwise kernel values between samples. Its working object is an n x n matrix, not a d x d matrix. This fact alone drives the memory cost, the solver choice, and the main limitation: Kernel PCA has no inverse_transform.

This page is based on src/machine_learning/decomposition/kernel_pca.rs, the shared kernel types in src/machine_learning/types.rs, and the integration tests in tests/machine_learning/kernel_pca.rs.

The kernel machinery, the KernelType enum and the Gamma coefficient, is the same code that drives Support Vector Machines. Kernel PCA reuses this code unchanged. It does not define its own kernel types.

2.11.1. The kernel trick, and why plain PCA cannot see rings

PCA finds the directions of maximum variance in the input space. This works well when the structure is linear. It fails when the structure is not linear.

The classic failure case is 2 concentric rings. The inner ring has radius 0.5. The outer ring has radius 3.0. Radius alone separates the 2 classes perfectly. No straight line separates them, so no linear projection can either.

Feed these points to PCA, and the leading components capture only the angular spread of the outer ring. The radial information that distinguishes the classes never appears.

The kernel trick maps each point x through a nonlinear feature map phi(x). This maps x into a much higher-dimensional space. Ordinary PCA then runs in that space. In this lifted space, the rings can become linearly separable.

Kernel PCA never builds phi(x) directly. This is the trick. PCA in the lifted space needs only inner products of the form phi(xi) * phi(xj). A kernel function K(xi, xj) computes this inner product directly, without building phi.

For the RBF kernel, K(x, y) = exp(-gamma * ||x - y||^2). Its implicit feature space has infinite dimensions. Each kernel value is still a single scalar, and you can compute it directly.

The RBF value depends only on the distance between 2 points. This encodes the radial structure that plain PCA discards. This is why the RBF kernel splits the 2 rings. 2.11.7 has a runnable example.

2.11.2. Double centering: the non-obvious core

PCA in the lifted space needs centered features: phi_centered(xi) = phi(xi) - (1/n) * sum_k phi(xk). You cannot subtract this mean directly, because you never have phi in hand.

The object Kernel PCA actually decomposes is the matrix of inner products of the centered features. You can write these inner products entirely in terms of the raw kernel matrix K. Expanding phi_centered(xi) * phi_centered(xj) gives the double-centering identity:

Kc[i, j] = K[i, j] - row_means[i] - row_means[j] + overall_mean

row_means[i] is the mean of row i of K. overall_mean is the mean of the whole matrix. In matrix form, this is Kc = H * K * H. H is the centering matrix H = I - (1/n) * J. I is the identity matrix, and J is the n x n all-ones matrix.

This is called double centering because H applies on both sides of K. It subtracts the row mean and the column mean, then adds the grand mean back, so the formula does not remove it twice. A custom Kernel PCA implementation often misses that + overall_mean term. Missing it biases every projection.

fit implements this. It computes the per-row means and the overall mean of the training kernel matrix (kernel_means). It then rewrites each entry in place as K[i,j] - row_mean[i] - row_mean[j] + overall_mean (center_kernel_matrix).

H * K * H has a direct consequence. Every row of Kc sums to zero, so every column of the resulting projection has a mean of zero. The test test_centering_training_output_has_near_zero_column_means checks that each projected component averages within 1e-9 of zero.

New points need a different, asymmetric formula. Many naive implementations give up here and refuse out-of-sample transforms.

When you project a new sample, its cross-kernel row against the training set needs the training statistics, not its own. center_cross_kernel_matrix subtracts the training row means and the mean of the new row itself, then adds the training overall mean. RustyML implements this, so transform on unseen data works correctly. See 2.11.6.

2.11.3. Constructing the estimator

The constructor takes the kernel and the number of components. It validates both and returns a Result:

pub fn new(kernel: KernelType, n_components: usize) -> Result<Self, Error>
pub fn with_eigen_solver(self, eigen_solver: EigenSolver) -> Self

n_components == 0 gets rejected right away, with Error::InvalidParameter naming the field. The relationship n_components <= n_samples cannot be checked at construction time, because the sample count is not known yet. fit enforces it instead, again with Error::InvalidParameter. A failed fit leaves the model untouched. The test test_fit_n_components_greater_than_n_samples_returns_invalid_parameter confirms the fitted-state getters stay None after a failed fit. You never end up with a half-mutated estimator.

The eigen solver defaults to EigenSolver::Dense. Set it with the chaining builder with_eigen_solver. The Default implementation gives you an RBF kernel with gamma = 0.1, n_components = 2, and the dense solver:

use rustyml::machine_learning::decomposition::kernel_pca::{EigenSolver, KernelPCA};
use rustyml::machine_learning::{Gamma, KernelType};
use ndarray::array;

fn main() {
    // 6 points in 2-D. RBF kernel, keep 2 components.
    let x = array![
        [1.0, 0.0],
        [0.0, 1.0],
        [-1.0, 0.0],
        [0.0, -1.0],
        [2.0, 0.5],
        [-0.5, 2.0],
    ];

    let mut kpca = KernelPCA::new(KernelType::RBF { gamma: Gamma::Value(0.5) }, 2)
        .unwrap()
        .with_eigen_solver(EigenSolver::Dense);

    let projected = kpca.fit_transform(&x).unwrap();
    assert_eq!(projected.nrows(), 6);
    assert_eq!(projected.ncols(), 2);

    // Fitted state is exposed through getters.
    println!("kept {} components", kpca.get_n_components());
    println!("training samples: {:?}", kpca.get_n_samples()); // Some(6)
    let eigenvalues = kpca.get_eigenvalues().unwrap();
    println!("leading eigenvalue: {}", eigenvalues[0]);
}

The getters mirror the internal state. get_kernel, get_n_components, and get_eigen_solver return by value. get_n_samples and get_n_features return Option<usize> (None before fitting). get_eigenvalues and get_eigenvectors return Option<&Array1<f64>> and Option<&Array2<f64>>.

The stored eigenvectors are the columns from the centered kernel matrix’s eigendecomposition. These are the per-sample coefficients, traditionally written alpha, with shape n_samples x n_components. They are not the input-space directions that PCA gives you. Kernel PCA has no meaningful loading vector to inspect. This is the flip side of working in an implicit space.

2.11.4. Kernels and choosing gamma

KernelType has 5 variants, shared with SVC:

Linear reduces Kernel PCA back to ordinary PCA, up to the centering convention. Use it only as a baseline. RBF is the default. Use RBF when you suspect nonlinear, distance-based structure, like the rings.

Poly captures polynomial interactions. Cosine normalizes away magnitude and keeps only direction. This helps with high-dimensional sparse data.

Sigmoid is not a true (Mercer) kernel. Its centered Gram matrix can be indefinite. This interacts with the eigenvalue handling in 2.11.6.

Constructor validation is strict, and specific to each kernel. Poly requires degree > 0, a positive finite gamma, and a finite coef0. RBF requires a positive finite gamma. Sigmoid only requires its parameters to be finite. gamma = 0 is accepted for Sigmoid (test test_new_sigmoid_gamma_zero_accepted), because a zero coefficient is a legitimate, if degenerate, sigmoid. Any violation returns Error::InvalidParameter, naming the field.

The gamma coefficient has type Gamma. It is either an explicit value or a data-dependent rule resolved at fit time:

fit resolves Scale and Auto once, using the training data’s variance and feature count. It stores the resolved value, so the training matrix and every later transform call use the same coefficient. Gamma::Scale fails with Error::InvalidInput if the data has zero variance (all-constant features), because the formula divides by it.

For the RBF kernel, gamma is an inverse squared bandwidth: gamma = 1 / (2 * sigma^2). A large gamma (small bandwidth) makes the kernel see only near neighbors. The Gram matrix approaches the identity matrix, every point looks maximally distinct, and the projection overfits noise. A small gamma (large bandwidth) makes every pair look similar. The Gram matrix approaches a constant matrix, and the leading component captures nothing.

The useful range sits where the typical value of gamma * ||x - y||^2, for nearby points, is close to 1. A practical starting point is Gamma::Scale. From there, decrease gamma if the projection looks like noise. Increase gamma if distinct clusters merge together. Kernel PCA is unsupervised, so it has no built-in cross-validation for this choice. The concentric-rings separability metric in the test suite (class_separability) is the kind of downstream signal you can tune against.

2.11.5. Eigen solvers: exact versus iterative

EigenSolver selects how RustyML extracts the top n_components eigenpairs of the centered kernel matrix. All 3 solvers are pure-Rust, in-house implementations. The crate does not depend on a second linear-algebra library for this step (nalgebra remains only a dev-dependency, for test cross-checks).

The solver choice affects speed and the numerical path, not the result. Dense, Lanczos, and PowerIteration agree on the leading eigenvalues. They produce projections that are identical up to a per-column sign flip. The tests confirm this directly: test_eigensolver_dense_vs_lanczos_agree and test_eigensolver_dense_vs_power_iteration_agree compare column norms (sign-agnostic) and the top eigenvalue, to tolerances of 1e-5 and 1e-4. The sign ambiguity comes from the eigenvectors themselves, not from a bug. Fix the sign yourself downstream if you need a canonical one.

The iterative solvers do not save memory. All 3 solvers operate on the same n x n centered kernel matrix, which must exist in full before any decomposition starts. Lanczos and PowerIteration skip the O(n^3) cost of a full decomposition, when n_components is far smaller than n_samples. Neither one avoids the O(n^2) matrix itself. 2.11.8 covers that limit.

// A few components from a large kernel matrix: skip the full O(n^3) decomposition.
let kpca = KernelPCA::new(KernelType::RBF { gamma: Gamma::Scale }, 3)
    .unwrap()
    .with_eigen_solver(EigenSolver::Lanczos);

Kernel PCA carries no randomness in any solver. There is no seed to set. 2 runs on the same data and machine produce bit-identical output. The test test_determinism_dense_solver asserts exact equality (assert_allclose(..., 0.0)). This differs from t-SNE, whose stochastic initialization needs seeding.

2.11.6. Fitting, transforming, and out-of-sample projection

The 3 entry points are inherent methods. KernelPCA also implements the crate’s Fit, Transform, and FitTransform traits. These traits forward to the inherent methods, so generic code can treat KernelPCA like PCA:

pub fn fit<S>(&mut self, x: &ArrayBase<S, Ix2>) -> Result<&mut Self, Error>
pub fn transform<S>(&self, x: &ArrayBase<S, Ix2>) -> Result<Array2<f64>, Error>
pub fn fit_transform<S>(&mut self, x: &ArrayBase<S, Ix2>) -> Result<Array2<f64>, Error>

fit needs at least 2 samples. 1 row returns Error::InvalidInput. 0 rows return Error::EmptyInput. Non-finite input returns Error::NonFinite.

fit resolves gamma, builds and centers the training kernel matrix, extracts the eigenpairs, and stores everything a later transform needs. This includes a copy of the full training matrix, which every transform call reuses.

transform projects any matrix with the same feature count as the training data, including data the model has never seen. It builds the cross-kernel matrix between the new points and the stored training samples. It centers this matrix with the training statistics (see 2.11.2), then projects it onto the stored eigenvectors.

The projected coordinate on component k is (Kc * v_k) / sqrt(lambda_k). This 1/sqrt(lambda) scaling turns raw eigenvectors into properly normalized principal components. Calling transform before fit returns Error::NotFitted. A feature-count mismatch returns Error::DimensionMismatch.

use rustyml::machine_learning::decomposition::kernel_pca::KernelPCA;
use rustyml::machine_learning::{Gamma, KernelType};
use ndarray::array;

fn main() {
    let x_train = array![
        [1.0, 0.0],
        [0.0, 1.0],
        [-1.0, 0.0],
        [0.0, -1.0],
        [2.0, 0.5],
        [-0.5, 2.0],
        [1.5, -1.5],
        [-2.0, 1.0],
    ];

    // Default solver is Dense. No builder call needed.
    let mut kpca = KernelPCA::new(KernelType::RBF { gamma: Gamma::Value(0.5) }, 2).unwrap();
    kpca.fit(&x_train).unwrap();

    // Points the model has never seen, same feature count.
    let x_new = array![[0.3, 0.3], [1.8, -1.2]];
    let projected_new = kpca.transform(&x_new).unwrap();
    assert_eq!(projected_new.nrows(), 2);
    assert_eq!(projected_new.ncols(), 2);
    println!("{projected_new:?}");
}

fit_transform is a shortcut. It fits the model, then transforms the same data. The test test_fit_transform_equals_fit_then_transform confirms it matches the 2-call path to 1e-10.

fit_transform is a convenience, not an optimization. Internally, it just calls fit, then transform, on the same matrix. transform rebuilds the kernel matrix from scratch either way, so the cost is the same as calling both separately. Use the 2-step form when you need to project a different matrix through an already-fitted model.

For a proper Mercer kernel on distinct points, the centered Gram matrix is positive semidefinite. Every kept eigenvalue is strictly positive (test test_eigenvalues_are_positive_after_fit also confirms they come out sorted in descending order).

Kernel PCA does not hard-fail on non-positive eigenvalues, though. A centered Gram matrix is only positive semidefinite up to round-off error. Non-Mercer kernels, like Sigmoid, can produce genuinely negative trailing eigenvalues. fit only rejects non-finite eigenvalues (NaN or Inf maps to Error::Computation), instead of rejecting the whole fit. Any component whose eigenvalue is not meaningfully positive, below a relative 1e-12 * lambda_max threshold, gets a projection scale of 0.0. This zeroes that column instead of producing Inf or NaN.

This keeps the requested n_components dimensionality, while quietly discarding degenerate directions. The test test_fit_indefinite_kernel_negative_eigenvalue_is_tolerated drives this path with a Sigmoid kernel. It checks that the offending column comes out all zeros, and the rest stays finite.

2.11.7. Worked example: separating concentric rings

This is the case plain PCA cannot handle. There are 2 rings, radially separable but linearly entangled. Run the RBF kernel, and the 2 classes land in distinguishable regions of component space:

use rustyml::machine_learning::decomposition::kernel_pca::KernelPCA;
use rustyml::machine_learning::{Gamma, KernelType};
use ndarray::Array2;
use std::f64::consts::PI;

fn main() {
    // Inner ring r = 0.5, outer ring r = 3.0. No line separates them.
    let n = 12;
    let mut data: Vec<f64> = Vec::new();
    for i in 0..n {
        let a = 2.0 * PI * i as f64 / n as f64;
        data.push(0.5 * a.cos());
        data.push(0.5 * a.sin());
    }
    for i in 0..n {
        let a = 2.0 * PI * i as f64 / n as f64;
        data.push(3.0 * a.cos());
        data.push(3.0 * a.sin());
    }
    let x = Array2::from_shape_vec((2 * n, 2), data).unwrap();

    let mut kpca = KernelPCA::new(KernelType::RBF { gamma: Gamma::Value(0.5) }, 2).unwrap();
    let proj = kpca.fit_transform(&x).unwrap();

    // The RBF kernel encodes radial distance, so the rings separate along a component.
    let inner_mean: f64 = (0..n).map(|i| proj[[i, 0]]).sum::<f64>() / n as f64;
    let outer_mean: f64 = (n..2 * n).map(|i| proj[[i, 0]]).sum::<f64>() / n as f64;
    println!("inner-ring mean of component 0: {inner_mean:.4}");
    println!("outer-ring mean of component 0: {outer_mean:.4}");
    println!("gap between ring means: {:.4}", (inner_mean - outer_mean).abs());
}

Swap KernelType::RBF { .. } for KernelType::Linear, and the gap between the ring means collapses. The linear projection is dominated by the outer ring’s angular variation. It never encodes the radius. The test test_rbf_separates_radial_clusters_better_than_linear makes this quantitative, with a Fisher-style separability score. It asserts that the RBF projection beats the linear one by a comfortable margin.

2.11.8. What Kernel PCA cannot do

Kernel PCA has no inverse_transform. Plain PCA has one. You can map a low-dimensional code back to the input space, because the projection is a linear map with a clean transpose.

Kernel PCA cannot do this. This is not an oversight. It is the pre-image problem. A projected point lives in the implicit feature space. To invert it, you need an input x whose feature map phi(x) lands at that location.

For most kernels, RBF above all, the feature map is nonlinear, infinite-dimensional, and not surjective. An arbitrary point in feature space usually has no exact pre-image. It has only approximate ones, found through a separate nonlinear optimization.

RustyML does not ship that approximation. Kernel PCA is strictly a forward, one-way projection. Use it for visualization, for denoising by projection, or as a nonlinear feature stage that feeds a downstream classifier. Do not use it for reconstruction.

The Gram matrix is O(n^2). That is the real ceiling. fit builds an n x n matrix of f64 values. Memory grows as 8 * n^2 bytes, regardless of feature count. This is roughly 800 MB at n = 10,000, and 3.2 GB at n = 20,000.

Time is worse. Building the matrix is O(n^2 * d), through a single parallel GEMM. The dense eigendecomposition is O(n^3). Switching to Lanczos or PowerIteration trims the decomposition cost, when you need only a handful of components. Nothing removes the O(n^2) matrix itself.

In practice, Kernel PCA is comfortable into the low thousands of samples. It starts to hurt in the tens of thousands. Past that point, subsample a representative set to fit on, then transform the rest (each new batch pays O(m * n * d) for this). Or use a method that never forms the full kernel matrix.

Every transform call also carries the training set with it. The projection is defined relative to the stored training samples, so transform rebuilds an m x n cross-kernel matrix for m new points. PCA’s transform cost is independent of the training size. Kernel PCA’s transform cost scales with n forever. Budget for it.

The parallel machinery starts automatically above internal size gates, keyed on the element count of the kernel matrix. Roughly, the centering scans parallelize once n^2 clears a few hundred thousand elements. The elementwise centering parallelizes once it clears a few million elements. The kernel GEMM has its own FLOPs gate. You do not configure any of this per call. See 7.3. Performance Tuning and Parallelism for the tunable gates.

2.11.9. Persistence

KernelPCA derives Serialize and Deserialize, and exposes the standard pair:

pub fn save_to_path(&self, path: &str) -> Result<(), Error>
pub fn load_from_path(path: &str) -> Result<Self, Error>

Serialization uses the compact postcard binary format. The .bin, .dat, or any other extension in the path is just a filename. The bytes are always binary.

A round-tripped model reproduces transform output exactly. The test test_save_load_round_trip asserts equality to 1e-12.

Check what gets serialized. A fitted Kernel PCA stores the entire training matrix, along with the eigenvectors and the centering statistics, because transform needs all of it. The saved file grows with your training set. This is another consequence of the same O(n^2)/stored-samples design. Remember this before you persist a model fit on a large corpus.

use rustyml::machine_learning::decomposition::kernel_pca::KernelPCA;
use rustyml::machine_learning::{Gamma, KernelType};
use ndarray::array;
use std::fs;

fn main() {
    let x = array![
        [1.0, 0.0],
        [0.0, 1.0],
        [-1.0, 0.0],
        [0.0, -1.0],
        [2.0, 0.5],
        [-0.5, 2.0],
        [1.5, -1.5],
        [-2.0, 1.0],
    ];

    let mut kpca = KernelPCA::new(KernelType::RBF { gamma: Gamma::Value(0.5) }, 2).unwrap();
    kpca.fit(&x).unwrap();
    let before = kpca.transform(&x).unwrap();

    let path = "kpca_model.bin";
    kpca.save_to_path(path).unwrap();
    let loaded = KernelPCA::load_from_path(path).unwrap();
    let after = loaded.transform(&x).unwrap();

    assert_eq!(before.shape(), after.shape());
    fs::remove_file(path).unwrap();
}

A missing file surfaces as Error::Io (test test_load_from_nonexistent_path_returns_io_error). For the general error taxonomy, see 1.6. Error Handling. For persistence patterns across the crate, see 7.2. Model Persistence in Depth.