Nonlinear Component Analysis as a Kernel Eigenvalue Problem

Scholkopf B, Smola A J, Muller K, et al. Nonlinear component analysis as a kernel eigenvalue problem[J]. Neural Computation, 1998, 10(5): 1299-1319.

引

普通的PCA將下式進(jìn)行特征分解（用論文的話講就是對角化）：
$C = \frac{1}{M} \sum \limits_{j=1}^M x_j x_j^T$
其中 $x_j \in \mathbb{R}^{N}, j = 1, \ldots, M$ 葵萎，且 $\sum \limits_{j=1}^M x_j = 0$ (中心化)导犹。

而kernel PCA試圖通過一個非線性函數(shù)：
$\Phi:\mathbb{R}^N \rightarrow F, x \rightarrow X$
其中 $F$ 是一個高維空間（甚至是無限維)。
所以我們要解決這么一個問題：
$\bar{C} = \frac{1}{M} \sum_{j=1}^M \Phi (x_j) \Phi(x_j)^T$

其實我們面對的第一個問題不是維度的問題而是 $\Phi$ 的選擇或者說構(gòu)造羡忘。我們?yōu)槭裁匆褦?shù)據(jù)映射到高維的空間谎痢？因為當(dāng)前數(shù)據(jù)的結(jié)構(gòu)（或者說分布）并不理想。

比如滿足 $(x-1)^2+(y-1)^2=4$ 的點卷雕，我們可以擴充到高維空間 $(x^2, x, y, y^2)$ 节猿，在高維空間是線性的（雖然這個例子用在kernel SVM 比較好）。

因為 $\Phi(\cdot)$ 的構(gòu)造蠻麻煩的漫雕，即便有一些先驗知識滨嘱。我們來看一種比較簡單的泛用的映射：
$(x_1, x_2, x_3) \rightarrow (x_1^3, x_2^3, x_3^3, x_1^2x_2,x_1^2x_3,x_1x_2^2,x_1x_3^2,x_2^2x_3,x_2x_3^2,x_1x_2x_3)$
這種樣子的映射，很容易把維度擴充到很大很大蝎亚，這個時候求解特征問題會變得很麻煩九孩。

kernel PCA

假設(shè) $\sum \limits_{i=1}^M \Phi(x_i)=0$ (如何保證這個性質(zhì)的成立在最后講，注意即便 $\sum \limits_{i=1}^M x_i = 0$ 发框， $\sum \limits_{i=1}^M \Phi(x_i)=0$ 也不一定成立)躺彬。

假設(shè)我們找到了 $\bar{C}$ 的特征向量 $V \ne 0$ :
$\bar{C}V = \lambda V$
因為 $V$ 是 $\Phi(x_i),i=1,\ldots, M$ 的線性組合（這個容易證明），所以梅惯， $V$ 可以由下式表示：
$V = \sum \limits_{i=1}^M \alpha_i \Phi(x_i)$

所以：
$\lambda V^T \Phi(x_j) = V^T\bar{C} \Phi(x_j), \quad for \: all \: j=1,\ldots, M$
等價于(記 $\Phi = [\Phi(x_1), \ldots, \Phi(x_M)]$ )：
$\begin{array}{ll} \lambda \sum \limits_{i=1}^M \alpha_i (\Phi^T(x_i)\Phi(x_j)) &= \lambda \{ \Phi^T \Phi(x_j)\} ^T \alpha \\ & =\frac{1}{M} \sum \limits_{i=1}^M \alpha_i \Phi^T(x_i) \Phi \Phi^T \Phi(x_j) \\ & = \frac{1}{M} \{\Phi^T \Phi \Phi^T \Phi(x_j)\}^T \alpha \end{array}$
對于 $j=1,\ldots, M$ 均成立宪拥，其中 $\alpha = [\alpha_1, \ldots, \alpha_M]^T$ 。

等價于：
$M \lambda \Phi^T \Phi \alpha = \Phi^T \Phi \Phi^T \Phi \alpha$
令 $K = \Phi^T \Phi$ 铣减，那么可寫作：
$M \lambda K \alpha = K^2\alpha$
其中 $K_{ij} = \Phi^T(x_i) \Phi(x_j)$

所以她君，我們可以通過下式來求解 $\alpha$ :
$M\lambda \alpha = K \alpha$
即 $\alpha$ 是 $K$ 的特征向量(注意，當(dāng) $\alpha$ 為特征向量的時候是一定符合 $M \lambda K \alpha = K^2\alpha$ 的葫哗，反之也是的缔刹，即二者是等價的)球涛。

假設(shè) $\lambda_1 \ge \lambda_2 \ge \ldots \ge \lambda_M$ 對應(yīng) $\alpha^1, \ldots, \alpha^M$ ，那么相應(yīng)的 $V$ 也算是求出來了校镐。

需要注意的是亿扁， $\|\alpha\|$ 往往不為1，因為我們希望 $\|V\|=1$ 鸟廓，所以：
$V^TV = \alpha^T K \alpha = \lambda \|\alpha\|^2 = 1$
所以 $\|\alpha\| = \frac{1}{\sqrt{\lambda}}$

PCA當(dāng)然需要求主成分从祝，假設(shè)有一個新的樣本 $x$ ，我們需要求：
$\Phi(x)^TV = \Phi^T(x) \Phi \alpha = \sum \limits_{i=1}^M \alpha_i \Phi^T(x_i) \Phi(x)$

注意引谜，我們只需要計算 $\Phi^T(x_i) \Phi(x)$ 牍陌。

現(xiàn)在回到kernel PCA 上的關(guān)鍵kernel上。注意到员咽，無論是K毒涧，還是最后計算主成分，我們都只需要計算 $\Phi^T(x)\Phi(y)$ 就可以了骏融，所以如果我們能夠找到一個函數(shù) $k(x,y)$ 來替代就不必顯示將 $x$ 映射到 $\Phi(x)$ 了链嘀，這就能夠避免了 $\Phi(\cdot)$ 的選擇問題和計算問題萌狂。

kernel 的選擇

顯然档玻，PCA的 $\lambda \ge 0$ ，所以我們也必須保證 $K$ 為半正定矩陣茫藏，相應(yīng)的核函數(shù) $k$ 稱為正定核误趴，Mercer定理有相應(yīng)的構(gòu)建。

也有現(xiàn)成的正定核：

多項式核

$k(x, y) = (x^Ty + 1)^d$
論文中是 $(x^Ty)^d$

高斯核函數(shù)

$k(x, y) = \exp \{\ -\frac{\|x-y\|^2}{2\sigma^2}\}$

性質(zhì)

在這里插入圖片描述

論文用上面的一個例子來說明务傲，kernel PCA可能更準(zhǔn)確地抓住數(shù)據(jù)的結(jié)構(gòu)凉当。

kernel PCA具有普通PCA的性質(zhì)，良好的逼近（從方差角度）售葡，以及擁有最多的互信息等等看杭。并且，如果 $k(x, y) = k(x^Hy)$ 挟伙，那么kernel PCA還具有酉不變性楼雹。

因為普通的PCA處理的是一個 $N \times N$ 的協(xié)方差矩陣，所以尖阔，至多獲得 $N$ 個載荷向量贮缅，而kernel PCA至多獲得 $M$ 個載荷向量(特征值非零)。所以介却，kernel PCA有望比普通PCA更加精準(zhǔn)谴供。

一些問題

中心化

PCA處理的是協(xié)方差矩陣，正如我們最開始所假設(shè)的齿坷， $\sum \limits_{i=1}^{M} \Phi(x_i)=0$ ,即中心化桂肌。因為 $\Phi(\cdot)$ 并不是線性函數(shù)数焊，所以，即便 $\sum \limits_{i=1}^M x_i = 0$ 也不能保證 $\sum \limits_{i=1}^{M} \Phi(x_i)=0$ 崎场，不過有別的方法處理昌跌。
令
$\tilde{\Phi}(x_i) = \Phi(x_i) - \frac{1}{M}\sum \limits_{k=1}^M \Phi(x_k) \\ \tilde{K}_{ij} = \tilde{\Phi}^T(x_i) \Phi(x_j) \\ 1_{M} = \{1\}_{ij}^{M \times M}$
可以得到：
$\begin{array}{ll} \tilde{K}_{ij} &= \tilde{\Phi}^T(x_i) \Phi(x_j) \\ &= \big(\Phi(x_i) - \frac{1}{M}\sum \limits_{k=1}^M \Phi(x_k)\big)^T \big(\Phi(x_j) - \frac{1}{M}\sum \limits_{k=1}^M \Phi(x_k)\big) \\ &= K_{ij} - \frac{1}{M} \sum \limits_{k=1}^M K_{kj} - \frac{1}{M} \sum \limits_{k=1}^M K_{ik} + \frac{1}{M^2} \sum \limits \limits_{m,n=1}^M K_{mn} \\ &= (K - 1_MK - K1_M + 1_MK1_M)_{ij} \end{array}$
于是，我們通過 $K$ 可以構(gòu)造出 $\tilde{K}$ 照雁。只需再求解 $\tilde{K}\tilde{\alpha} = M \lambda \tilde{\alpha}$ 即可蚕愤。

恢復(fù)

我們知道，根據(jù)PCA選出的載荷向量以及主成分饺蚊，我們能夠恢復(fù)出原數(shù)據(jù)（或者近似萍诱，如果我們只選取了部分載荷向量）。對于kernel PCA污呼，比較困難裕坊，因為我們并沒有顯式構(gòu)造 $\Phi(\cdot)$ ，也就沒法顯式找到 $V$ 燕酷，更何況籍凝，有時候我們高維空間找到 $V$ 在原空間中并不存在原像。
或許苗缩，我們可以通過：
$\min \limits_{x} \quad \|\Phi(x) - \Phi(\hat{x})\|^2$
來求解饵蒂，注意到，上式也只和核函數(shù) $k(x,y)$ 有關(guān)酱讶。

代碼


import numpy as np

class KernelPCA:

    def __init__(self, data, kernel=1, pra=3):
        self.__n, self.__d = data.shape
        self.__data = np.array(data, dtype=float)
        self.kernel = self.kernels(kernel, pra)
        self.__K = self.center()

    @property
    def shape(self):
        return self.__n, self.__d

    @property
    def data(self):
        return self.data

    @property
    def K(self):
        return self.__K

    @property
    def alpha(self):
        return self.__alpha

    @property
    def eigenvalue(self):
        return self.__value

    def kernels(self, label, pra):
        """
        數(shù)據(jù)是一維的時候可能有Bug
        :param label: 1:多項式;2:exp
        :param pra:1: d; 2: sigma
        :return: 函數(shù) 或報錯
        """
        if label is 1:
            return lambda x, y: (x @ y) ** pra
        elif label is 2:
            return lambda x, y: \
                np.exp(-(x-y) @ (x-y) / (2 * pra ** 2))
        else:
            raise TypeError("No such kernel...")

    def center(self):
        """中心化"""
        oldK = np.zeros((self.__n, self.__n), dtype=float)
        one_n = np.ones((self.__n, self.__n), dtype=float)
        for i in range(self.__n):
            for j in range(i, self.__n):
                x = self.__data[i]
                y = self.__data[j]
                oldK[i, j] = oldK[j, i] = self.kernel(x, y)
        return oldK - 2 * one_n @ oldK + one_n @ oldK @ one_n

    def processing(self):
        """實際上就是K的特征分解,再對alpha的大小進(jìn)行一下調(diào)整"""
        value, alpha = np.linalg.eig(self.__K)
        index = value > 0
        value = value[index]
        alpha = alpha[:, index] * (1 / np.sqrt(value))
        self.__alpha = alpha
        self.__value = value / self.__n

    def score(self, x):
        """來了一個新的樣本退盯，我們進(jìn)行得分"""
        k = np.zeros(self.__n)
        for i in range(self.__n):
            y = self.__data[i]
            k[i] = self.kernel(x, y)
        return k @ self.__alpha





"""

import matplotlib.pyplot as plt

x = np.linspace(-1, 1, 100)
y = x ** 2 + [np.random.randn() * 0.1 for i in range(100)]
data = np.array([x, y]).T

test = KernelPCA(data, pra=1)
test.processing()
print(test.alpha.shape)
print(test.alpha[:, 0])

"""