\(\newcommand{\ci}{\perp\!\!\!\perp}\) \(\newcommand{\cA}{\mathcal{A}}\) \(\newcommand{\cB}{\mathcal{B}}\) \(\newcommand{\cC}{\mathcal{C}}\) \(\newcommand{\cD}{\mathcal{D}}\) \(\newcommand{\cE}{\mathcal{E}}\) \(\newcommand{\cF}{\mathcal{F}}\) \(\newcommand{\cG}{\mathcal{G}}\) \(\newcommand{\cH}{\mathcal{H}}\) \(\newcommand{\cI}{\mathcal{I}}\) \(\newcommand{\cJ}{\mathcal{J}}\) \(\newcommand{\cK}{\mathcal{K}}\) \(\newcommand{\cL}{\mathcal{L}}\) \(\newcommand{\cM}{\mathcal{M}}\) \(\newcommand{\cN}{\mathcal{N}}\) \(\newcommand{\cO}{\mathcal{O}}\) \(\newcommand{\cP}{\mathcal{P}}\) \(\newcommand{\cQ}{\mathcal{Q}}\) \(\newcommand{\cR}{\mathcal{R}}\) \(\newcommand{\cS}{\mathcal{S}}\) \(\newcommand{\cT}{\mathcal{T}}\) \(\newcommand{\cU}{\mathcal{U}}\) \(\newcommand{\cV}{\mathcal{V}}\) \(\newcommand{\cW}{\mathcal{W}}\) \(\newcommand{\cX}{\mathcal{X}}\) \(\newcommand{\cY}{\mathcal{Y}}\) \(\newcommand{\cZ}{\mathcal{Z}}\) \(\newcommand{\bA}{\mathbf{A}}\) \(\newcommand{\bB}{\mathbf{B}}\) \(\newcommand{\bC}{\mathbf{C}}\) \(\newcommand{\bD}{\mathbf{D}}\) \(\newcommand{\bE}{\mathbf{E}}\) \(\newcommand{\bF}{\mathbf{F}}\) \(\newcommand{\bG}{\mathbf{G}}\) \(\newcommand{\bH}{\mathbf{H}}\) \(\newcommand{\bI}{\mathbf{I}}\) \(\newcommand{\bJ}{\mathbf{J}}\) \(\newcommand{\bK}{\mathbf{K}}\) \(\newcommand{\bL}{\mathbf{L}}\) \(\newcommand{\bM}{\mathbf{M}}\) \(\newcommand{\bN}{\mathbf{N}}\) \(\newcommand{\bO}{\mathbf{O}}\) \(\newcommand{\bP}{\mathbf{P}}\) \(\newcommand{\bQ}{\mathbf{Q}}\) \(\newcommand{\bR}{\mathbf{R}}\) \(\newcommand{\bS}{\mathbf{S}}\) \(\newcommand{\bT}{\mathbf{T}}\) \(\newcommand{\bU}{\mathbf{U}}\) \(\newcommand{\bV}{\mathbf{V}}\) \(\newcommand{\bW}{\mathbf{W}}\) \(\newcommand{\bX}{\mathbf{X}}\) \(\newcommand{\bY}{\mathbf{Y}}\) \(\newcommand{\bZ}{\mathbf{Z}}\) \(\newcommand{\ba}{\mathbf{a}}\) \(\newcommand{\bb}{\mathbf{b}}\) \(\newcommand{\bc}{\mathbf{c}}\) \(\newcommand{\bd}{\mathbf{d}}\) \(\newcommand{\be}{\mathbf{e}}\) \(\newcommand{\bg}{\mathbf{g}}\) \(\newcommand{\bh}{\mathbf{h}}\) \(\newcommand{\bi}{\mathbf{i}}\) \(\newcommand{\bj}{\mathbf{j}}\) \(\newcommand{\bk}{\mathbf{k}}\) \(\newcommand{\bl}{\mathbf{l}}\) \(\newcommand{\bm}{\mathbf{m}}\) \(\newcommand{\bn}{\mathbf{n}}\) \(\newcommand{\bo}{\mathbf{o}}\) \(\newcommand{\bp}{\mathbf{p}}\) \(\newcommand{\bq}{\mathbf{q}}\) \(\newcommand{\br}{\mathbf{r}}\) \(\newcommand{\bs}{\mathbf{s}}\) \(\newcommand{\bt}{\mathbf{t}}\) \(\newcommand{\bu}{\mathbf{u}}\) \(\newcommand{\bv}{\mathbf{v}}\) \(\newcommand{\bw}{\mathbf{w}}\) \(\newcommand{\bx}{\mathbf{x}}\) \(\newcommand{\by}{\mathbf{y}}\) \(\newcommand{\bz}{\mathbf{z}}\) \(\newcommand{\RR}{\mathbb{R}}\) \(\newcommand{\NN}{\mathbb{N}}\) \(\newcommand{\balpha}{\boldsymbol{\alpha}}\) \(\newcommand{\bbeta}{\boldsymbol{\beta}}\) \(\newcommand{\btheta}{\boldsymbol{\theta}}\) \(\newcommand{\hpi}{\widehat{\pi}}\) \(\newcommand{\bpi}{\boldsymbol{\pi}}\) \(\newcommand{\hbpi}{\widehat{\boldsymbol{\pi}}}\) \(\newcommand{\bxi}{\boldsymbol{\xi}}\) \(\newcommand{\bmu}{\boldsymbol{\mu}}\) \(\newcommand{\bepsilon}{\boldsymbol{\epsilon}}\) \(\newcommand{\bzero}{\mathbf{0}}\) \(\newcommand{\T}{\text{T}}\) \(\newcommand{\Trace}{\text{Trace}}\) \(\newcommand{\Cov}{\text{Cov}}\) \(\newcommand{\Var}{\text{Var}}\) \(\newcommand{\E}{\mathbb{E}}\) \(\newcommand{\Pr}{\text{Pr}}\) \(\newcommand{\pr}{\text{pr}}\) \(\newcommand{\pdf}{\text{pdf}}\) \(\newcommand{\P}{\text{P}}\) \(\newcommand{\p}{\text{p}}\) \(\newcommand{\One}{\mathbf{1}}\) \(\newcommand{\argmin}{\operatorname*{arg\,min}}\) \(\newcommand{\argmax}{\operatorname*{arg\,max}}\) \(\newcommand{\dtheta}{\frac{\partial}{\partial\theta} }\) \(\newcommand{\ptheta}{\nabla_\theta}\) \(\newcommand{\alert}[1]{\color{darkorange}{#1}}\) \(\newcommand{\alertr}[1]{\color{red}{#1}}\) \(\newcommand{\alertb}[1]{\color{blue}{#1}}\)

1 Motivation

Principal Component Analysis (PCA) is a widely used technique for dimensionality reduction and feature extraction. However, traditional PCA is limited to linear relationships in the data. Kernel Principal Component Analysis (KPCA) extends PCA to capture nonlinear relationships by using kernel methods. This allows KPCA to find more complex structures in the data, making it a powerful tool.

2 The Two-Moon Example

To illustrate the effectiveness of KPCA, we will use the two-moon dataset, which consists of two interleaving half circles. This dataset is a classic example where traditional PCA fails to capture the underlying structure due to its nonlinear nature.

  library(ggplot2)
  library(dplyr)
  set.seed(546)
  
  #generate two-moon data
  n_per_moon <- 200
  theta <- runif(n_per_moon, 0, pi)
  x1 <- cbind(cos(theta), sin(theta)) + matrix(rnorm(2*n_per_moon, sd=0.05), n_per_moon, 2)
  x2 <- cbind(1 - cos(theta), 1 - sin(theta) - 0.5) + matrix(rnorm(2*n_per_moon, sd=0.05), n_per_moon, 2)
  
  data <- rbind(x1, x2) %>% as.data.frame()
  colnames(data) <- c("x1", "x2")
  data$label <- factor(rep(c("Class 1", "Class 2"), each = n_per_moon))
  
  # plot the data
  ggplot(data, aes(x = x1, y = x2, color = label)) +  
    geom_point(size = 2) + 
    theme_minimal() + 
    labs(title = "Two-Moon Dataset", x = "X1", y = "X2") +   
    theme(legend.title = element_blank())

A standard PCA on this dataset cannot unfold the two curved structures: any linear projection loses part of their geometry. We should note that many successful tools for such a dataset require a nonlinear transformation, such as spectral clustering.

Let’s first quickly review the idea of PCA. Suppose a data matrix \(\bX \in \RR^{n \times p}\) is centered, i.e., each column has mean zero. Let \(\bx_i\in\RR^p\) denote its \(i\)th row as a column vector. PCA seeks to find the directions that maximize the variance of the projected data. This is equivalent to solving the following optimization problem for the first principal direction:

\[ \begin{aligned} & \max_{\bv \in \RR^p, \lVert \bv \rVert = 1} \frac{1}{n}\bv^\T \bX^\T \bX \bv \\ =& \max_{\bv \in \RR^p, \lVert \bv \rVert = 1} \frac{1}{n} \sum_{i=1}^n (\bx_i^\T \bv)^2 \\ =& \max_{\bv \in \RR^p, \lVert \bv \rVert = 1} \frac{1}{n} \sum_{i=1}^n \bv^\T (\bx_i \bx_i^\T) \bv \\ =& \max_{\bv \in \RR^p, \lVert \bv \rVert = 1} \bv^\T \left( \frac{1}{n} \sum_{i=1}^n \bx_i \bx_i^\T \right) \bv \\ =& \max_{\bv \in \RR^p, \lVert \bv \rVert = 1} \bv^\T \boldsymbol{\Sigma} \bv \end{aligned} \]

where \(\boldsymbol{\Sigma} = \frac{1}{n} \bX^\T \bX = \frac{1}{n} \sum_{i=1}^n \bx_i \bx_i^\T\) is the sample covariance matrix. The solution is the eigenvector corresponding to the largest eigenvalue of \(\boldsymbol{\Sigma}\), i.e., \(\boldsymbol{\Sigma}\bv=\lambda\bv\). Subsequent principal directions use the remaining eigenvectors.

3 A New Perspective of PCA

In the original PCA, we are living in the Euclidean space \(\RR^p\). The concepts of distance and inner product are defined in the standard way. However, if we want to perform nonlinear PCA, we can first lift each input \(x_i\) to a feature representation in a higher-dimensional space. In our case, this space can be an infinite-dimensional RKHS, where \(x_i\) is represented by the kernel section \(K(\cdot,x_i)\). Let’s think about each concept in the PCA setting and find its counterpart in the RKHS setting.

For this conceptual comparison, write \(f_i\) for the centered feature representation of \(x_i\).

  • Data points: In PCA, each data point is a vector in \(\RR^p\). In the RKHS setting, each data point represents a function \(f_i\), an infinite-dimensional vector in the RKHS.
  • A direction: In PCA, a direction is a vector in \(\RR^p\), the same space as the data points. In the RKHS setting, a direction is also a function in the RKHS.
  • Projecting data onto a direction: In PCA, projecting a data point \(\bx_i\) onto a direction \(\bv\) is done by computing the inner product \(\bx_i^\T\bv\). In the RKHS setting, projecting a function \(f_i\) onto a direction \(g\) is done by computing the inner product \(\langle f_i,g\rangle_\cH\), which is also called the score.
  • Variance of the projected data: In PCA, the variance of the projected data onto a direction \(\bv\) is given by \(\frac{1}{n}\sum_{i=1}^n(\bx_i^\T\bv)^2\). In the RKHS, the variance is given by \(\frac{1}{n}\sum_{i=1}^n\langle f_i,g\rangle_\cH^2\).

Here comes a tricky part: how do we compute the covariance in the RKHS? In PCA, we have the covariance matrix \(\boldsymbol{\Sigma} = \frac{1}{n} \sum_{i=1}^n \bx_i \bx_i^\T\), which is a \(p \times p\) matrix. In the RKHS, we define a covariance operator that maps functions to functions:

\[ \Sigma_\cH = \frac{1}{n} \sum_{i=1}^n f_i \otimes f_i. \]

Here, \(f_i \otimes f_i\) is the outer-product operator defined by \((f_i \otimes f_i)(g) = \langle f_i, g \rangle_{\cH} f_i\) for any \(g \in \cH\). The covariance operator \(\Sigma_\cH\) is linear and maps functions in the RKHS to functions in the RKHS. What does this even mean? Let’s use the feature-map representation to understand this idea.

4 The Feature Mapping View

Based on our previous understanding of RKHS, we can create a feature map \(\Phi\) that maps the data points to a higher-dimensional space where linear PCA can be applied. For the moment, let’s use a finite-dimensional representation; the RKHS view below extends the same calculation to an infinite-dimensional space:

\[ \Phi(x) = (\phi_1(x), \phi_2(x), \ldots, \phi_m(x))^\T \]

We will search for a direction in the span of the observed feature representations, i.e.,

\[ \bv = \sum_{i=1}^n \alpha_i \Phi(x_i). \]

Why can we do this? Any component of a direction that is orthogonal to every \(\Phi(x_i)\) contributes zero to every projected score. Therefore, any principal direction with positive variance can be chosen in their span. In this high-dimensional space, say \(\cG\), we can define the inner product and variance in the same way as before. Following our previous discussion, the PCA problem can be formulated as:

\[ \max_{\bv \in \cG, \lVert \bv \rVert_\cG = 1} \bv^\T \left( \frac{1}{n} \sum_{i=1}^n \Phi(x_i) \Phi(x_i)^\T \right) \bv. \]

where \(\frac{1}{n} \sum_{i=1}^n \Phi(x_i) \Phi(x_i)^\T\) represents the covariance matrix when the feature representations are centered. We will return to this centering step shortly. Now we can plug in our representation of \(\bv\) into the above optimization problem. The objective function is

\[ \begin{aligned} & \left( \sum_{i=1}^n \alpha_i \Phi(x_i) \right)^\T \left( \frac{1}{n} \sum_{j=1}^n \Phi(x_j) \Phi(x_j)^\T \right) \left( \sum_{k=1}^n \alpha_k \Phi(x_k) \right) \\ =& \frac{1}{n} \sum_{i=1}^n \sum_{j=1}^n \sum_{k=1}^n \alpha_i \alpha_k \Phi(x_i)^\T \Phi(x_j) \Phi(x_j)^\T \Phi(x_k) \\ =& \frac{1}{n} \sum_{j=1}^n \left( \sum_{i=1}^n \alpha_i \Phi(x_i)^\T \Phi(x_j) \right) \left( \sum_{k=1}^n \alpha_k \Phi(x_j)^\T \Phi(x_k) \right) \\ =& \frac{1}{n} \sum_{j=1}^n \left( \sum_{i=1}^n \alpha_i K(x_i,x_j) \right)^2 \\ =& \frac{1}{n} \sum_{j=1}^n \left( \sum_{i=1}^n \alpha_i K_{ij} \right) \left( \sum_{k=1}^n \alpha_k K_{kj} \right) \\ =& \frac{1}{n} \sum_{i,k} \alpha_i \alpha_k \left( \sum_j K_{ij} K_{kj} \right) \\ =& \frac{1}{n} \sum_{i,k} \alpha_i \alpha_k (\bK^2)_{ik} \\ =& \frac{1}{n} \balpha^\T \bK^2 \balpha \end{aligned} \]

where \(\bK\in\RR^{n\times n}\) is the kernel matrix with \(K_{ij}=K(x_i,x_j)=\Phi(x_i)^\T\Phi(x_j)\). The constraint can be simplified as:

\[ \left\lVert \sum_{i=1}^n \alpha_i \Phi(x_i) \right\rVert_\cG^2 = \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j K(x_i,x_j) = \balpha^\T \bK \balpha = 1. \]

Therefore, the optimization problem can be rewritten as:

\[ \begin{aligned} \max_{\balpha \in \RR^n} \quad & \frac{1}{n} \balpha^\T \bK^2 \balpha\\ \text{ subj. to } \quad & \balpha^\T \bK \balpha = 1. \end{aligned} \]

This is a generalized Rayleigh quotient problem. We can also view it through the Lagrangian

\[ \mathcal{L}(\balpha, \lambda) = \frac{1}{n} \balpha^\T \bK^2 \balpha - \lambda (\balpha^\T \bK \balpha - 1). \]

Taking the derivative with respect to \(\balpha\) and setting it to zero gives

\[ \bK^2\balpha=n\lambda\bK\balpha. \]

Because \(\bK\) may be singular, we cannot simply cancel it from both sides. Instead, let

\[ \bK\bu_j=\rho_j\bu_j, \qquad \|\bu_j\|=1, \qquad \rho_j>0. \]

The coefficient vector and the corresponding unit-norm principal direction are

\[ \balpha_j=\frac{\bu_j}{\sqrt{\rho_j}}, \qquad \bv_j=\sum_{i=1}^n\frac{u_{ij}}{\sqrt{\rho_j}}\Phi(x_i). \]

This normalization is important: \(\balpha_j^\T\bK\balpha_j=1\). The variance captured by this direction is \(\rho_j/n\), while the vector of training scores is

\[ \bK\balpha_j=\sqrt{\rho_j}\bu_j. \]

For a new data point \(x\), the corresponding score is

\[ s_j(x)=\langle \Phi(x),\bv_j\rangle_\cG =\frac{1}{\sqrt{\rho_j}}\sum_{i=1}^n u_{ij}K(x,x_i). \]

There is still one important issue: the feature representations used by PCA must be centered.

5 Double-Centering of the Kernel Matrix

Specifically, we need to make sure the data in the feature space are centered. Otherwise, PCA can find a direction that primarily points toward the feature-space mean. In the feature-map view, we center each observation as

\[ \Phi_c(x_i) = \Phi(x_i) - \frac{1}{n} \sum_{\ell=1}^n \Phi(x_\ell) = \Phi(x_i) - \overline{\Phi}. \]

This leads to the following centering procedure on the kernel matrix:

\[ \begin{aligned} K_c(i,j) &= \langle \Phi(x_i)-\overline{\Phi},\; \Phi(x_j)-\overline{\Phi}\rangle \\ &= K(x_i,x_j) - \left\langle \Phi(x_i), \frac{1}{n}\sum_{k=1}^n \Phi(x_k) \right\rangle - \left\langle \frac{1}{n}\sum_{k=1}^n \Phi(x_k), \Phi(x_j) \right\rangle\\ &\quad + \left\langle \frac{1}{n}\sum_{k=1}^n \Phi(x_k), \frac{1}{n}\sum_{m=1}^n \Phi(x_m) \right\rangle \\ &= K(x_i,x_j) - \frac{1}{n}\sum_{k=1}^n K(x_i,x_k) - \frac{1}{n}\sum_{k=1}^n K(x_k,x_j) + \frac{1}{n^2}\sum_{k=1}^n\sum_{m=1}^n K(x_k,x_m) \end{aligned} \]

This takes out the mean of each row and each column and adds back the overall mean, a procedure known as double centering. In matrix form, we can write it as

\[ \begin{aligned} \bK_c &= \bH \bK \bH \\ &= \big(\bI - \tfrac{1}{n}\One \One^\T\big)\, \bK \,\big(\bI - \tfrac{1}{n}\One \One^\T\big) \\ &= \bK \;-\; \tfrac{1}{n}\,\One \One^\T \bK \;-\; \tfrac{1}{n}\,\bK \One \One^\T \;+\; \tfrac{1}{n^2}\,\One \One^\T \bK \One \One^\T, \end{aligned} \]

All formulas above now apply with \(\bK\) replaced by \(\bK_c\). Therefore, the final step of KPCA is to solve

\[ \bK_c\bu_j=\rho_j\bu_j, \qquad \|\bu_j\|=1, \qquad \rho_j>0. \]

The normalized direction uses \(\balpha_j=\bu_j/\sqrt{\rho_j}\), and the training-score vector is \(\sqrt{\rho_j}\bu_j\).

6 The Two-Moon Example Revisited

Now we can apply KPCA to the two-moon dataset. Let’s use a Gaussian kernel. We will compute the kernel matrix, center it, and then solve the eigenvalue problem to obtain the principal components.

  # compute Gaussian kernel using the dist() function
  dist_matrix <- as.matrix(dist(data[, c("x1", "x2")]))
  sigma <- 0.2
  K <- exp(-dist_matrix^2 / (2 * sigma^2))
  
  # center the kernel matrix
  n_total <- nrow(K)
  H <- diag(n_total) - matrix(1, n_total, n_total) / n_total
  Kc <- H %*% K %*% H
  Kc <- (Kc + t(Kc)) / 2
  
  # eigen decomposition
  eig <- eigen(Kc, symmetric = TRUE)
  positive <- eig$values > max(eig$values) * 1e-10
  eig_values <- eig$values[positive]
  eig_vectors <- eig$vectors[, positive, drop = FALSE]
  alpha <- sweep(eig_vectors, 2, sqrt(eig_values), "/")
  scores <- Kc %*% alpha

  # plot the first two principal component scores
  kpca_data <- data.frame(PC1 = scores[, 1], PC2 = scores[, 2], label = data$label)
  ggplot(kpca_data, aes(x = PC1, y = PC2, color = label)) +  
    geom_point(size = 2) + 
    theme_minimal() + 
    labs(title = "Kernel PCA on Two-Moon Dataset", x = "KPCA Score 1", y = "KPCA Score 2") +
    theme(legend.title = element_blank())

7 Predicting New Data Points

To predict the principal-component score for a new data point \(x_{\mathrm{new}}\), first define its kernel vector relative to the training data:

\[ \bk(x_{\mathrm{new}}) =\big(K(x_1,x_{\mathrm{new}}),\ldots,K(x_n,x_{\mathrm{new}})\big)^\T. \]

Because the training feature representations were centered, this new kernel vector must be centered using the training-sample means:

\[ \begin{aligned} \bk_c(x_{\mathrm{new}}) &=\bk(x_{\mathrm{new}})-\frac{1}{n}\bK\One -\frac{1}{n}\One\One^\T\bk(x_{\mathrm{new}})\\ &\quad+\frac{1}{n^2}\One\One^\T\bK\One. \end{aligned} \]

Here, \(\One\in\RR^n\) is the vector of ones. The direction for component \(j\) is

\[ \bv_j=\sum_{i=1}^n\alpha_{ij}\Phi_c(x_i), \qquad \balpha_j=\frac{\bu_j}{\sqrt{\rho_j}}. \]

The score for the new data point is then

[ \[\begin{aligned} s_j(x_{\mathrm{new}}) &=\langle \Phi(x_{\mathrm{new}})-\overline{\Phi},\bv_j\rangle_\cG\\ &=\balpha_j^\T\bk_c(x_{\mathrm{new}}). \end{aligned}\]

8 The RKHS View

In fact, we can also view KPCA in the RKHS without explicitly constructing a feature map. Recall that each data point \(x_i\) is mapped to the kernel section \(h_i=K(\cdot,x_i)\). Define the centered kernel sections by

\[ h_i=K(\cdot,x_i), \qquad \overline h=\frac{1}{n}\sum_{r=1}^n h_r, \qquad f_i=h_i-\overline h. \]

The covariance operator in the RKHS is

\[ \Sigma_\cH=\frac{1}{n}\sum_{i=1}^n f_i\otimes f_i. \]

This corresponds to the covariance matrix in feature space, but it maps functions to functions. For any \(f,h,g\in\cH\), the outer-product operator works as

\[ (f\otimes h)(g)=\langle h,g\rangle_\cH f. \]

The PCA problem in the RKHS can be formulated as:

\[ \begin{aligned} &\max_{g \in \cH, \lVert g \rVert_\cH = 1} \langle g, \Sigma_\cH g \rangle_\cH \\ =& \max_{g \in \cH, \lVert g \rVert_\cH = 1} \left\langle g, \left\{\frac{1}{n} \sum_{i=1}^n f_i \otimes f_i\right\}(g) \right\rangle_\cH \\ =& \max_{g \in \cH, \lVert g \rVert_\cH = 1} \frac{1}{n} \sum_{i=1}^n \langle f_i, g \rangle_\cH^2 \end{aligned} \]

Since the functions \(f_i\) are centered, we have \(\frac{1}{n}\sum_{i=1}^n f_i=0\), which means \(\frac{1}{n}\sum_{i=1}^n\langle f_i,g\rangle_\cH=0\). The projected data therefore have mean zero. More generally, their variance is

\[ \Var\big(\langle f_i,g\rangle_\cH\big)=\frac{1}{n}\sum_{i=1}^n\langle f_i,g\rangle_\cH^2-\left(\frac{1}{n}\sum_{i=1}^n\langle f_i,g\rangle_\cH\right)^2. \]

This is exactly the variance of the projected data onto the direction \(g\). Now we can optimize \(g\) to maximize this variance. If a principal direction has covariance-operator eigenvalue \(\nu_j=\rho_j/n>0\), then

\[ \Sigma_\cH g=\nu_j g \quad\Longrightarrow\quad g=\frac{1}{n\nu_j}\sum_{i=1}^n\langle f_i,g\rangle_\cH f_i. \]

Therefore, the eigenvalue equation directly shows that every positive-variance direction lies in the span of the centered kernel sections. We can write

\[ g=\sum_{i=1}^n\alpha_i f_i. \]

Plugging this representation into the PCA problem, we have: \[ \begin{aligned} \langle g, \Sigma_\cH g \rangle_\cH &= \frac{1}{n} \sum_{i=1}^n \left\langle f_i, \sum_{j=1}^n \alpha_j f_j \right\rangle_\cH^2 \\ &= \frac{1}{n} \sum_{i=1}^n \left( \sum_{j=1}^n \alpha_j \langle f_i, f_j \rangle_\cH \right)^2 \\ &= \frac{1}{n} \sum_{i=1}^n \left( \sum_{j=1}^n \alpha_j K_c(i,j) \right)^2 \\ &= \frac{1}{n} \balpha^\T \bK_c^2 \balpha \end{aligned} \]

The unit-norm constraint is

\[ \|g\|_\cH^2=\balpha^\T\bK_c\balpha=1. \]

And we get exactly the same optimization problem as before. Hence the feature-map view can be skipped if we directly work in the RKHS. For a new data point \(x_{\mathrm{new}}\), define its centered kernel section as

\[ f_{\mathrm{new}}=K(\cdot,x_{\mathrm{new}})-\overline h. \]

Its score along \(g=\sum_{i=1}^n\alpha_i f_i\) is

\[ \begin{aligned} s(x_{\mathrm{new}}) &=\langle f_{\mathrm{new}},g\rangle_\cH\\ &=\sum_{i=1}^n\alpha_i\langle f_{\mathrm{new}},f_i\rangle_\cH\\ &=\balpha^\T\bk_c(x_{\mathrm{new}}). \end{aligned} \]

We can visualize the first KPCA score function in our two-moon example:

  # Create a grid of points
  x_seq <- seq(-1.5, 2.5, length.out = 100)
  y_seq <- seq(-1, 1.5, length.out = 100)
  grid <- expand.grid(x = x_seq, y = y_seq)
  
  # Compute the centered score function on the grid
  train_x <- as.matrix(data[, c("x1", "x2")])
  alpha1 <- alpha[, 1]
  train_kernel_mean <- rowMeans(K)
  overall_kernel_mean <- mean(K)

  compute_score <- function(point, train_x, alpha, train_kernel_mean,
                            overall_kernel_mean, sigma) {
    squared_distance <- rowSums(sweep(train_x, 2, point, "-")^2)
    k_new <- exp(-squared_distance / (2 * sigma^2))
    k_new_centered <- k_new - train_kernel_mean - mean(k_new) +
      overall_kernel_mean
    sum(alpha * k_new_centered)
  }
  
  grid$z <- apply(
    grid, 1, compute_score,
    train_x = train_x,
    alpha = alpha1,
    train_kernel_mean = train_kernel_mean,
    overall_kernel_mean = overall_kernel_mean,
    sigma = sigma
  )
  
  # Plot the function
  ggplot(grid, aes(x = x, y = y)) +
    geom_raster(aes(fill = z), interpolate = TRUE) +
    scale_fill_viridis_c() +
    geom_contour(aes(z = z), color = "white", alpha = 0.5) +
    geom_point(data = data, aes(x = x1, y = x2, color = label), size = 1) +
    theme_minimal() +
    labs(title = "First KPCA Score Function", x = "X1", y = "X2") +
    theme(legend.title = element_blank())

9 Other Unsupervised Learning Methods in RKHS

Many unsupervised learning methods that can be formulated through inner products or squared distances admit kernel versions, including kernel \(K\)-means and several spectral methods. However, each method requires its own derivation. In KPCA, centering introduced the empirical feature-space mean

\[ \widehat\mu=\frac{1}{n}\sum_{i=1}^n\Phi(x_i). \]

The next lecture studies this mean directly through kernel mean embeddings and compares two such means using MMD.

For the original formulation of KPCA, see Schölkopf, Smola, and Müller (1998), Nonlinear Component Analysis as a Kernel Eigenvalue Problem.