Kernel PCA

September 2026 · 17 min read

Kernel PCA performs PCA in a feature space defined by a kernel. The map into that space can be nonlinear, while the PCA operation there remains linear. A kernel supplies the pairwise dot products needed for the calculation without requiring us to build every mapped feature.

We will use five scalar inputs and a feature map small enough to draw. That makes it possible to check the kernel calculation against ordinary PCA on explicitly constructed features before moving to more general kernels.

First change the representation

Take inputs $x\in\{-2,-1,0,1,2\}$. Map each into two coordinates, $\phi(x)=(x,x^2)$. The mapped points lie on a parabola. Their mean is (0, 2), and ordinary PCA on these mapped points chooses the vertical direction first:

This example expands one input into two features and keeps one score. It is designed to expose the mechanism, not to claim useful compression of an already one-dimensional input. In a larger problem, the map can contain many features or even be infinite-dimensional.

Replace explicit features with their dot products

For two scalar inputs $a$ and $b$, the dot product of their mapped features can be computed directly from the inputs:

For example, $k(-2,1)=(-2)(1)+4(1)=2$. Collect one such value for each pair of training inputs into the Gram matrix $K$, where $K_{ij}=k(x_i,x_j)$. Its rows and columns index observations, not feature dimensions.

The kernel changes which geometry PCA sees. Here it accounts for both a signed value and its square. A positive semidefinite kernel produces a Gram matrix with nonnegative eigenvalues, up to numerical error, for any finite dataset. That is the property expected of genuine dot products; an arbitrary similarity table need not have it.

Center in the mapped space

The original inputs already have mean zero, but the mapped features have mean (0, 2). Centering $x$ before squaring therefore does not center $\phi(x)$. Kernel PCA must remove the mapped mean.

We can do this using matrix averages. Let $r_i$ be row $i$'s mean in the raw training matrix, $c_j$ column $j$'s mean, and $g$ the mean of all entries. The centered Gram entry is

For this symmetric matrix, the row and column mean lists are both [8, 2, 0, 2, 8], and the grand mean is 4. Click a cell to follow its calculation:

Turn the centered matrix into scores

Find the positive eigenvalues of $K_c$ and their unit-length eigenvectors. Write one pair as $(\eta_j,v_j)$. Each entry of $v_j$ belongs to a training observation. Scaling these entries by $\sqrt{\eta_j}$ gives the training scores:

Our two positive Gram eigenvalues are 14 and 10. With the first direction oriented toward increasing $x^2$, its five scores are [2, −1, −2, −1, 2], exactly the explicit feature-space PCA result. The second component recovers the centered linear feature, up to sign.

Gram eigenvalues are not yet sample covariance eigenvalues. Divide by $N-1$ for the convention used in our PCA chapter: with five rows, the feature-space variances are 3.5 and 2.5. The first retains $14/(14+10)=7/12$, or about 58.3%, of the variance in this mapped space.

Project a new point using the training means

Save projection weights $\alpha=v/\sqrt{\eta}$ for the component we want to keep. Here $v$ is its unit Gram eigenvector and $\eta$ its positive eigenvalue. The weights tell us how to combine a new point’s comparisons with the training inputs.

For a new input $x$, compute its kernel values against the saved training inputs. Subtract the mean of that new row and the stored training column means, then add the stored training grand mean. Finally take the dot product with the fitted weights $\alpha$.

These operations do not require explicitly constructing $\phi(x)$. That is the point of the kernel formulation. The full fit below uses the explicit map only as an independent check.

import numpy as np
from sklearn.decomposition import KernelPCA

x = np.array([-2, -1, 0, 1, 2], dtype=float)
def kernel(a, b):
    product = float(a @ b)
    return product + product ** 2

# The explicit map is available only to check this small example.
Phi = np.column_stack([x, x*x])
K = np.array([[kernel(np.array([a]), np.array([b])) for b in x] for a in x])
np.testing.assert_allclose(K, Phi @ Phi.T)
column_mean = K.mean(axis=0)
grand_mean = K.mean()
Kc = K - K.mean(axis=1)[:, None] - column_mean[None, :] + grand_mean
Phi_centered = Phi - Phi.mean(axis=0)
np.testing.assert_allclose(Kc, Phi_centered @ Phi_centered.T)

values, vectors = np.linalg.eigh(Kc)
order = np.argsort(values)[::-1]
values, vectors = values[order], vectors[:, order]
positive = values > 1e-10  # Numerical cutoff for this example's scale.
values, vectors = values[positive], vectors[:, positive]
# Orient the first component toward increasing x^2, purely for display.
if vectors[:, 0] @ (x*x - 2) < 0:
    vectors[:, 0] *= -1
alpha = vectors[:, 0] / np.sqrt(values[0])
print(values)             # [14, 10]: Gram eigenvalues.
print(Kc @ alpha)         # [2, -1, -2, -1, 2]: first scores.
print(values / (len(x)-1))  # [3.5, 2.5]: feature-space sample variances.

def transform(new):
    cross = np.array([[kernel(np.array([a]), np.array([b])) for b in x]
                      for a in np.asarray(new)])
    centered = cross - cross.mean(axis=1)[:, None] - column_mean[None, :] + grand_mean
    return centered @ alpha

print(transform([3]))  # [7]; use the stored training means.
reference = KernelPCA(n_components=2, kernel=kernel, eigen_solver="dense").fit(x[:, None])
reference_scores = reference.transform(x[:, None])[:, 0]
sign = np.sign(reference_scores @ (Kc @ alpha))
np.testing.assert_allclose(sign * reference_scores, Kc @ alpha, atol=1e-10)
np.testing.assert_allclose(sign * reference.transform([[3]])[:, 0], transform([3]))

The kernel and its scale are modeling choices

A linear kernel reproduces ordinary centered PCA scores, up to valid axis choices. Nonlinear kernels ask PCA to preserve variation in a different representation. They do not automatically discover the representation that best serves a prediction task.

Fit scaling and kernel PCA inside each training fold. Choose the kernel, its parameters and the retained component count on development data, and assess the complete downstream pipeline on independent evaluation data. For precomputed kernels, training uses a train-by-train matrix and new data uses a new-by-train matrix in the same training-column order.

The usual dense formulation stores an $N\times N$ Gram matrix. At 10,000 training rows, one matrix of 64-bit numbers occupies about 800 MB before other arrays. A full dense eigendecomposition also becomes expensive. Reducing the explicit feature count does not remove this dependence on the number of training observations.

A score does not necessarily identify an original input

Our first score is $x^2-2$, so it cannot distinguish $x$ from $-x$. No inverse can recover a unique sign from that score alone. More generally, a reconstructed point in feature space need not correspond to any exact original input.

Scikit-learn can learn an approximate inverse with fit_inverse_transform=True. That is a separately fitted pre-image model, not the exact transpose reconstruction used by ordinary PCA. Evaluate its reconstruction quality separately if recovering inputs matters.

Optional implementation: apply a learned kernel projection

The Python example above fits the eigenvectors. This C++ companion implements training-mean preparation and new-point transformation using the resulting fitted coefficients. Keeping those steps explicit makes it easy to verify the train-by-train and new-by-train distinction.

#include <cmath>
#include <stdexcept>
#include <vector>
struct Projection {
    std::vector<double> training, column_mean, alpha;
    double grand_mean;
};
double kernel(double a, double b) {
    double product = a*b, value = product + product*product;
    if (!std::isfinite(value)) throw std::overflow_error("kernel overflow");
    return value;
}

// alpha comes from an eigenvector of the centered training Gram matrix,
// divided by the square root of its positive eigenvalue (see Python fit).
Projection prepare_projection(const std::vector<double>& training,
                              const std::vector<double>& alpha) {
    if (training.empty() || alpha.size() != training.size())
        throw std::invalid_argument("invalid model shape");
    for (double x : training) if (!std::isfinite(x))
        throw std::invalid_argument("finite training inputs required");
    for (double a : alpha) if (!std::isfinite(a))
        throw std::invalid_argument("finite fitted weights required");
    const auto n = training.size();
    Projection model{training, std::vector<double>(n, 0), alpha, 0};
    for (std::size_t j = 0; j < n; ++j) {
        for (std::size_t i = 0; i < n; ++i)
            model.column_mean[j] += kernel(training[i], training[j]) / n;
        model.grand_mean += model.column_mean[j] / n;
    }
    return model;
}

double transform(const Projection& model, double x) {
    if (!std::isfinite(x)) throw std::invalid_argument("finite new input required");
    std::vector<double> row(model.training.size());
    double row_mean = 0, score = 0;
    for (std::size_t j = 0; j < row.size(); ++j) {
        row[j] = kernel(x, model.training[j]);
        row_mean += row[j] / row.size();
    }
    for (std::size_t j = 0; j < row.size(); ++j) {
        double centered = row[j] - row_mean - model.column_mean[j] + model.grand_mean;
        score += centered * model.alpha[j];
    }
    return score;
}

References and further reading

• Schölkopf, Smola and Müller: Nonlinear Component Analysis as a Kernel Eigenvalue Problem — feature-space PCA through kernels; Appendix B develops centering.
• Scikit-learn: KernelPCA — fitted eigenpairs, kernel choices, transform inputs and approximate inversion.
• Scikit-learn: RBF kernel — the distance-based kernel and its gamma parameter.