2.3 Centered kernels and kernel PCA

Motivation

Ordinary PCA finds linear directions in the measured features. A curved pattern may become easier to describe after a nonlinear feature map: a scalar measurement can, for example, contribute both its value and its square. Applying PCA to those features gives coordinates that are nonlinear functions of the original measurement. Kernel PCA performs this calculation through inner products, so the transformed features need not be explicitly stored.

The important changes are mathematical rather than procedural. Centering now takes place in feature space. An eigenvector of the observation-by-observation Gram matrix is not itself a unit principal direction in that space. Finally, a new observation must be measured relative to the same feature mean as the training observations. Deriving these three facts gives a complete projection rule, including its normalization and every mean term.

Definitions

Let $x_1,\ldots,x_n\in\mathbb R^d$, with $n\ge2$, be training observations. As in ordinary PCA, observations occupy rows. A feature map $\phi$ sends each observation into an inner-product space. When this space has finite dimension $D$, write $F\in\mathbb R^{n\times D}$ with row $i$ equal to $\phi(x_i)^T$. The formulas using finite sums of feature vectors also apply when the feature space has infinitely many dimensions.

Definition

Kernel and Gram matrix

A kernel computes feature inner products:

$$k(x,y)=\langle\phi(x),\phi(y)\rangle,\qquad K_{ij}=k(x_i,x_j).$$

Thus $K=FF^T\in\mathbb R^{n\times n}$ in the finite-dimensional case. A positive semidefinite kernel is a symmetric function for which every finite Gram matrix satisfies $c^TKc\ge0$ for all real coefficient vectors $c$.

The matrix $K$ compares observations; $F^TF$ acts on feature directions. These are different matrices, although their positive eigenvalues agree. Working with $K$ is attractive when the feature space is much larger than the number of observations. Its entries describe similarities through the chosen feature geometry, rather than arbitrary resemblance scores.

Two standard choices, with parameters that preserve positive semidefiniteness, are the Gaussian radial basis kernel and the polynomial kernel:

$$k_{\mathrm{RBF}}(x,y)=\exp\!\left(-\gamma\|x-y\|_2^2\right),\quad\gamma>0, \qquad k_{\mathrm{poly}}(x,y)=(x^Ty+c)^q,\quad c\ge0,\quad q\in\{1,2,\ldots\}.$$

For the polynomial kernel, expansion expresses the function as a sum of inner products of monomial features with nonnegative weights. The Gaussian kernel can be written as a positive factor depending on $x$, a positive factor depending on $y$, and $\exp(2\gamma x^Ty)$. The exponential series is a sum of nonnegative multiples of polynomial inner-product kernels. This explains why the squared distance and the parameter sign matter. Positive semidefiniteness concerns quadratic forms; individual kernel entries can be negative, as they can for an ordinary linear inner product.

Definition

Centered features and centered kernel

Define the training-feature mean and centering matrix by

$$\bar\phi=\frac1n\sum_{i=1}^n\phi(x_i),\qquad \widetilde\phi(x)=\phi(x)-\bar\phi,\qquad H=I_n-\frac1n\mathbf1\mathbf1^T.$$

The centered Gram matrix has entries $\widetilde K_{ij}=\langle\widetilde\phi(x_i),\widetilde\phi(x_j)\rangle$. Here $\widetilde\phi(x)$ always subtracts the training mean, including when $x$ is a new observation.

For a nonlinear feature map, centering the raw observations does not generally center the features. The mean of a square, for example, is not the square of the mean. Kernel centering is therefore required even if each original feature column already has mean zero.

Theorem and proposition

Theorem

Feature centering gives a positive semidefinite centered Gram matrix

The centered Gram matrix is

$$\widetilde K=HKH =K-\frac1n\mathbf1\mathbf1^TK-\frac1nK\mathbf1\mathbf1^T +\frac1{n^2}\mathbf1\mathbf1^TK\mathbf1\mathbf1^T.$$

It is symmetric positive semidefinite, $\widetilde K\mathbf1=0$, and $\operatorname{rank}(\widetilde K)\le n-1$.

If $\widetilde F=HF$, then $\widetilde K=\widetilde F\widetilde F^T$. The corresponding sample covariance in feature space is $\widetilde F^T\widetilde F/(n-1)$. If a positive Gram eigenvalue is $\mu_j$, the associated covariance eigenvalue is $\mu_j/(n-1)$. Consequently, explained variance fractions are ratios of the positive $\mu_j$ themselves. These eigenvalues already equal squared singular values; squaring them again would change the variance criterion.

Theorem

Unit principal directions and projection scores

Suppose $\widetilde Ka_j=\mu_ja_j$, with $\mu_j>0$ and $a_j^Ta_j=1$. Choose the $a_j$ orthonormally. Then

$$\alpha_j=\frac{a_j}{\sqrt{\mu_j}},\qquad v_j=\sum_{i=1}^n\alpha_{ji}\widetilde\phi(x_i)$$

define orthonormal principal directions in feature space. For a training observation, its score on component $j$ is

$$z_{ij}=\langle\widetilde\phi(x_i),v_j\rangle =\sqrt{\mu_j}\,a_{ji}.$$

For any new observation $x_*$, form $\widetilde k_*=(\langle\widetilde\phi(x_i),\widetilde\phi(x_*)\rangle)_i$. Its component score is $z_{*j}=\alpha_j^T\widetilde k_*$.

The loading direction $v_j$, Gram eigenvector $a_j$, and coefficient vector $\alpha_j$ have different meanings. The last two both have $n$ entries, but their normalizations differ. A unit Euclidean norm for $a_j$ becomes a unit feature-space norm for $v_j$ only after division by $\sqrt{\mu_j}$.

A finite eigenproblem for nonlinear features

The principal direction solves the same variance problem as before, now in feature space:

$$\max_{\|v\|=1}\frac1{n-1} \sum_{i=1}^n\langle\widetilde\phi(x_i),v\rangle^2.$$

Every direction perpendicular to the span of the centered training features has zero score on every training observation. It therefore contributes zero training variance. All positive-variance directions can be found inside this span, whose dimension is at most $n-1$. This is the reason a finite Gram matrix suffices even for an infinite-dimensional feature map: the eigenproblem describes the part of feature space supported by the training observations. The nonlinear map changes that space's geometry, while the PCA step remains an orthogonal linear projection inside it.

For a retained number of components $s$, the total squared reconstruction error of the centered training features is

$$\sum_i\left\|\widetilde\phi(x_i)- \sum_{j=1}^s z_{ij}v_j\right\|^2 =\operatorname{tr}(\widetilde K)-\sum_{j=1}^s\mu_j.$$

The trace adds the squared lengths of the centered feature vectors, and the retained eigenvalues add the energy captured by their orthogonal projections. Dividing this energy by $n-1$ converts it to sample variance. Ratios do not need that division because the same factor appears in numerator and denominator. Retaining all positive components reproduces the training features exactly, even though the feature map may have many more coordinates.

This reconstruction is a feature-space statement. A projected feature vector need not equal $\phi(x)$ for any original input $x$, since a nonlinear map's image need not be a linear subspace. Kernel PCA supplies coordinates and an optimal feature-space projection; converting that projection back to an original measurement is a further problem. Its reconstruction error should therefore be interpreted in the kernel's geometry, rather than automatically as squared error in the original measured variables.

Kernel PCA also determines a distance geometry. The distance between two feature vectors follows from the two diagonal Gram entries and their pairwise inner product. Subtracting a shared mean leaves pairwise differences unchanged, so centering changes the coordinate origin without changing those feature-space distances. The choice of kernel introduces the nonlinear geometry; centering then positions that geometry around its training mean before the projection selects its dominant directions.

Proof sketch or proof idea

Proof

Deriving centering and normalization

Expand $\langle\phi(x_i)-\bar\phi,\phi(x_j)-\bar\phi\rangle$. The four terms are the raw kernel entry, two training averages, and the training grand average. Collecting them in matrix form gives $HKH$. Alternatively, $\widetilde F=HF$ immediately gives $\widetilde F\widetilde F^T=HFF^TH$. For any $c$, $c^T\widetilde Kc=\|\sum_i c_i\widetilde\phi(x_i)\|^2\ge0$. Since $H\mathbf1=0$ and $H$ has rank $n-1$, the null-vector and rank claims follow.

Next,

$$\langle v_j,v_\ell\rangle =\frac{a_j^T\widetilde Ka_\ell}{\sqrt{\mu_j\mu_\ell}} =\frac{\mu_\ell a_j^Ta_\ell}{\sqrt{\mu_j\mu_\ell}}.$$

This equals one when $j=\ell$ and zero otherwise. Multiplying $\widetilde K\alpha_j$ gives $\sqrt{\mu_j}a_j$, proving the training-score formula. The covariance operator applied to $v_j$ gives $\mu_jv_j/(n-1)$, so these are precisely feature-space PCA directions.

For the new-observation rule, define the fixed training row means, grand mean, and raw kernel vector by

$$m_i=\frac1n\sum_{\ell=1}^nK_{i\ell},\qquad m=\frac1{n^2}\sum_{i,\ell=1}^nK_{i\ell},\qquad k_{*i}=k(x_i,x_*).$$

Expanding the centered feature inner product gives

$$\widetilde k_{*i} =k_{*i}-m_i-\frac1n\sum_{\ell=1}^n k_{*\ell}+m.$$

The new vector's average appears because it computes $\langle\bar\phi,\phi(x_*)\rangle$ using the existing training mean. It does not add the new observation to the training set or recompute that mean. Saving the training observations, $m_i$, $m$, and the normalized coefficient vectors is enough to apply the same projection later.

Positive-eigenvalue eigenvectors have zero sum. Indeed, $\mathbf1^T\widetilde Ka_j=\mu_j\mathbf1^Ta_j=0$, so $\mathbf1^Ta_j=0$. This can simplify the score calculation: terms constant over the training index disappear when dotted with $\alpha_j$. The projection may thus also be written as

$$z_{*j}=\alpha_j^T(k_*-m_{\mathrm{row}}),\qquad m_{\mathrm{row}}=(m_1,\ldots,m_n)^T.$$

The stored row means remain essential. Their contribution is not generally constant across observations, and zero-sum coefficients do not remove it. The full four-term formula remains useful for constructing and interpreting the centered kernel vector itself.

Worked examples

Worked example

1. Raw centering does not center nonlinear features

Take scalar observations $x_1=-1$, $x_2=0$, $x_3=1$, and $\phi(x)=(x,x^2)^T$. Their raw mean is zero, but

$$F=\begin{bmatrix}-1&1\\0&0\\1&1\end{bmatrix},\qquad \bar\phi=\begin{bmatrix}0\\2/3\end{bmatrix},\qquad \widetilde F=\begin{bmatrix}-1&1/3\\0&-2/3\\1&1/3\end{bmatrix}.$$

The kernel is $k(x,y)=xy+(xy)^2$. Computing it without storing features gives

$$K=\begin{bmatrix}2&0&0\\0&0&0\\0&0&2\end{bmatrix},\qquad \widetilde K=\frac19 \begin{bmatrix}10&-2&-8\\-2&4&-2\\-8&-2&10\end{bmatrix}.$$

Each row of $\widetilde K$ sums to zero. Direct multiplication of $\widetilde F\widetilde F^T$ gives the same entries. Leaving the raw scalar data unchanged because its mean is zero would leave the square feature with mean $2/3$; it would solve an uncentered feature-space problem.

This simple map separates two nonlinear effects. The first feature preserves the sign of the scalar, while the second records its squared magnitude. The observations at minus one and one agree in the second feature even though they differ in the first. Kernel PCA can find variation along both of these feature directions, using only the pairwise values in the displayed Gram matrix.

Worked example

2. Eigenvalues, normalized coefficients, and training scores

For the centered matrix above, orthonormal eigenvectors can be chosen as

$$\mu_1=2,\quad a_1=\frac1{\sqrt2}(-1,0,1)^T,\qquad \mu_2=\frac23,\quad a_2=\frac1{\sqrt6}(1,-2,1)^T.$$

The remaining eigenvalue is zero, with eigenvector $\mathbf1/\sqrt3$. Dividing by the square roots of the positive eigenvalues gives

$$\alpha_1=\tfrac12(-1,0,1)^T,\qquad \alpha_2=\tfrac12(1,-2,1)^T,\qquad Z=\begin{bmatrix}-1&1/3\\0&-2/3\\1&1/3\end{bmatrix}.$$

The columns of $Z$ are $\sqrt{\mu_1}a_1$ and $\sqrt{\mu_2}a_2$. In the explicit feature space, $\widetilde F^T\alpha_1=(1,0)^T$ and $\widetilde F^T\alpha_2=(0,1)^T$, confirming unit directions. Sample score variances are $1$ and $1/3$ because $n-1=2$. The first component explains $2/(2+2/3)=3/4$ of feature-space variation. Its discarded squared reconstruction error in feature space is $2/3$.

The eigenvalues of the centered Gram matrix are already nonnegative energies. They correspond to $\sigma_j^2$ in the SVD of the centered feature table. This is why ordinary PCA's squared-singular-value formula becomes an unsquared Gram-eigenvalue formula. The two methods agree exactly when the same explicit features and the same centering are used.

Worked example

3. Projecting a new observation with the training mean

Keep the same three training observations and consider $x_*=2$. The raw kernel vector, stored row means, and stored grand mean are

$$k_*=(2,0,6)^T,\qquad (m_1,m_2,m_3)^T=(2/3,0,2/3)^T,\qquad m=4/9.$$

The new vector's average is $8/3$. Thus

$$\widetilde k_*=(-8/9,-20/9,28/9)^T,\qquad (z_{*1},z_{*2})=(2,10/3).$$

For example, $\alpha_1^T\widetilde k_*=(8/9+28/9)/2=2$. The explicit centered feature vector is also $\phi(2)-\bar\phi=(2,4-2/3)^T=(2,10/3)^T$. This independently checks both coordinates. Including the new point when recomputing the mean would instead change the fitted feature space coordinates for every observation.

Computing the projection

The following code uses only the kernel formula, then compares the magnitudes of the new scores with the exact example. An eigensolver may reverse the sign of an eigenvector; the corresponding coordinate reverses with it.

import numpy as np

def kernel(x, y):
    t = x * y
    return t + t**2

train = np.array([-1., 0., 1.])
K = kernel(train[:, None], train[None, :])
row_mean = K.mean(axis=1)
grand_mean = K.mean()
Kc = K - row_mean[:, None] - row_mean[None, :] + grand_mean
mu, a = np.linalg.eigh(Kc)
order = np.argsort(mu)[::-1]
mu, a = mu[order], a[:, order]
tol = 1e-12 * np.linalg.norm(Kc, ord=2)
assert mu.min() >= -tol
positive = mu > tol
mu, a = mu[positive], a[:, positive]
alpha = a / np.sqrt(mu)
training_scores = Kc @ alpha
raw_new = kernel(train, 2.)
centered_new = raw_new - row_mean - raw_new.mean() + grand_mean
new_scores = centered_new @ alpha
print(mu)                  # [2., 0.66666667]
print(np.abs(new_scores))  # [2., 3.33333333]

An exact zero eigenvalue represents a zero linear combination of centered feature vectors, so it cannot supply a unit feature direction by division. It is omitted. In floating-point arithmetic, an eigenvalue very close to zero can have a small negative sign; classify it relative to the matrix scale. A substantially negative eigenvalue has a different meaning: the computed matrix is not positive semidefinite to that accuracy. Replacing it with its absolute value would change the geometry rather than normalize the same PCA.

Kernel choice determines which nonlinear variations receive large feature norms. For an RBF kernel, increasing $\gamma$ makes similarity decay over shorter original-space distances. For polynomial features, degree and input scale change the relative contributions of different powers. These parameters define the feature metric; maximizing its variance alone does not select them or optimize separation of labeled classes.

Units also interact directly with the Gaussian parameter. Multiplying every input coordinate by a scalar $b\ne0$ multiplies squared distances by $b^2$. Replacing $\gamma$ by $\gamma/b^2$ then gives exactly the same Gram matrix. Keeping $\gamma$ fixed instead selects a different similarity geometry. This simple identity helps distinguish a change of measurement units from an intentional change of the model's length scale.

Common mistakes

Common mistake

Confusing Gram eigenvectors with feature-space directions

Using $a_j$ directly as the coefficient vector creates a feature direction with squared norm $\mu_j$, not one. The correct coefficient is $a_j/\sqrt{\mu_j}$. Likewise, centering a new kernel vector by subtracting only its own average misses the stored training row means and grand mean. All four terms in the centered-kernel formula arise from one feature-space inner product.

Summary

A kernel supplies a positive semidefinite observation Gram matrix. Centering its features gives $HKH$. Positive eigenpairs provide normalized coefficients $a_j/\sqrt{\mu_j}$ and training scores $\sqrt{\mu_j}a_j$. A new observation uses kernel values centered with the existing training-feature mean. The resulting coordinates describe nonlinear features of the original data; their variance and reconstruction error are measured in that feature space.

Exercises

Checkpoint

1. For three observations with Gram matrix $K=I$, find the centered Gram matrix and its eigenvalues. What fraction of the feature-space variation does one principal component retain?

Use the centering matrix and its action on constant and zero-sum vectors.

Checkpoint

2. For six training observations, suppose a centered Gram matrix has positive eigenvalues $9$ and $4$. For a unit eigenvector $a$ associated with $9$, find its normalized coefficient vector, training score vector, and sample score variance.

Distinguish the Gram eigenvalue from the covariance eigenvalue.

Checkpoint

3. In the polynomial example, project the new observation $x=-2$ onto both displayed principal directions using the fixed training mean. Give its raw kernel vector, centered kernel vector, and scores.

The stored row means and grand mean remain unchanged.

Solutions

Solution · 1. A centered identity Gram matrix

Since $H^2=H$, $HKH=H=I_3-\mathbf1\mathbf1^T/3$. Its eigenvalues are $1,1,0$: it fixes every zero-sum vector and annihilates the constant vector. One component retains $1/(1+1)=1/2$ of the variation. The leading direction is not unique because the positive eigenvalue is repeated.

Solution · 2. The three normalizations

The coefficient vector is $\alpha=a/3$, the training score vector is $3a$, and its sample variance is $9/(6-1)=9/5$. The other component's variance is $4/5$. Gram eigenvalues must be divided by $n-1$ for sample variances, while the explained fraction for the first component is $9/13$.

Solution · 3. A new point on the opposite side

The raw vector is $(6,0,2)^T$, with average $8/3$. Subtracting the training row means and adding $4/9$ gives $\widetilde k_*=(28/9,-20/9,-8/9)^T$. Dotting with the two coefficient vectors gives scores $(-2,10/3)$. This agrees with $\phi(-2)-\bar\phi=(-2,4-2/3)^T$.

Related notes

• SVD and principal components
• Classical and nonmetric MDS
• Fisher discriminant and class scatter
• Course reading