# 2.11. Kernel PCA
`KernelPCA` runs [PCA](./2.10._Principal_Component_Analysis.md) 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](./2.5._Support_Vector_Machines.md). 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](#2117-worked-example-separating-concentric-rings) 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:
```text
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](#2116-fitting-transforming-and-out-of-sample-projection).
## 2.11.3. Constructing the estimator
The constructor takes the kernel and the number of components. It validates both and returns a `Result`:
```rust,ignore
pub fn new(kernel: KernelType, n_components: usize) -> Result<Self, Error>
pub fn with_eigen_solver(self, eigen_solver: EigenSolver) -> Self
```
| Parameter | Type | Meaning |
| --- | --- | --- |
| `kernel` | `KernelType` | The kernel function and its parameters. RustyML validates it up front. See [2.11.4](#2114-kernels-and-choosing-gamma). |
| `n_components` | `usize` | How many leading components to keep. Must be `> 0`. At fit time, it must also be `<= n_samples`. |
`n_components == 0` gets rejected right away, with [`Error::InvalidParameter`](../Chapter-01/1.6._Error_Handling.md) 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:
```rust
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:
| Variant | Formula | Parameters |
| --- | --- | --- |
| `Linear` | `K(x, y) = x*y` | none |
| `Poly { degree, gamma, coef0 }` | `(gamma*x*y + coef0)^degree` | `degree: u32` (`> 0`), `gamma: Gamma`, `coef0: f64` |
| `RBF { gamma }` | `exp(-gamma*||x - y||^2)` | `gamma: Gamma` |
| `Sigmoid { gamma, coef0 }` | `tanh(gamma*x*y + coef0)` | `gamma: Gamma`, `coef0: f64` |
| `Cosine` | `(x*y) / (||x||*||y||)` | none |
`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](#2116-fitting-transforming-and-out-of-sample-projection).
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:
| `Gamma` variant | Resolves to | Use when |
| --- | --- | --- |
| `Gamma::Value(v)` | `v` | You have a specific bandwidth in mind. |
| `Gamma::Scale` | `1 / (n_features * Var(X))` | A default that adapts to feature spread (scikit-learn's `'scale'`). |
| `Gamma::Auto` | `1 / n_features` | A simpler `1/d` rule (scikit-learn's `'auto'`). |
`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).
| Variant | Strategy | Best for |
| --- | --- | --- |
| `Dense` (default) | Full symmetric eigendecomposition (Householder tridiagonalization, then implicit-shift QL), the classic EISPACK/JAMA algorithm pair. Takes the leading pairs after the full decomposition. | Small to mid-sized kernel matrices, where you can afford the full `O(n^3)` decomposition. |
| `Lanczos` | Krylov-subspace iteration with full reorthogonalization. Reduces the problem to a small tridiagonal problem and solves that exactly with the same dense solver. | A few leading components of a large kernel matrix. |
| `PowerIteration` | Power iteration with Hotelling deflation, one component at a time. | The simplest iterative option. Use it when Lanczos is more than you need. |
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](#2118-what-kernel-pca-cannot-do) covers that limit.
```rust,ignore
// 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](./2.12._t-SNE.md), whose stochastic initialization needs [seeding](../Chapter-07/7.1._Reproducibility_and_Random_Seeds.md).
## 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`:
```rust,ignore
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](#2112-double-centering-the-non-obvious-core)), 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`.
```rust
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:
```rust
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](./2.10._Principal_Component_Analysis.md) 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](../Chapter-07/7.3._Performance_Tuning_and_Parallelism.md) for the tunable gates.
## 2.11.9. Persistence
`KernelPCA` derives `Serialize` and `Deserialize`, and exposes the standard pair:
```rust,ignore
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.
```rust
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](../Chapter-01/1.6._Error_Handling.md). For persistence patterns across the crate, see [7.2. Model Persistence in Depth](../Chapter-07/7.2._Model_Persistence_in_Depth.md).