Kernelized Principal Component Analysis.

主成分分析是一种线性降维方法，如果原始数据是线性不可分的，则降维后的数据仍然是线性不可分的。

核主成分分析(kernelized principal component analysis,KPCA)是一种非线性降维方法。其基本思想是先通过核方法把原始的线性不可分数据映射到高维空间，使其在高维空间线性可分；再通过主成分分析将其线性降维到低维空间，并在该低维空间中仍然线性可分。

KPCA的基本流程是对于维度为$d$的原始样本空间中的样本点$x \in \Bbb{R}^d$，引入映射$\phi$将其投影到维度为$d’$的高维特征空间$z=\phi(x) \in \Bbb{R}^{d’}$，再在特征空间中实施PCA方法，得到降维结果$y=Wz \in \Bbb{R}^{k}$。

1. 直接求解KPCA

KPCA中对样本点$x \in \Bbb{R}^d$引入映射$\phi$投影到高维特征空间$z=\phi(x) \in \Bbb{R}^{d’}$。映射$\phi$通常不显式给出，而是定义核函数计算高维特征空间中的内积：

常用的核函数包括：

• 线性核：$\kappa(x_i,x_j)=(ax_1^Tx_2+c)^d$
• 高斯核：$\kappa(x_i,x_j)=\exp(-\gamma ||x_1-x_2||^2)$
• Sigmoid核：$\kappa(x_i,x_j)=\tanh(\alpha x_1^Tx_2+\gamma)$

进一步引入矩阵$K \in \Bbb{R}^{n\times n}$，其中第$(i,j)$个元素$\kappa(x_i,x_j)$表示样本$x_i$和$x_j$映射到高维空间后的特征相似度。选定核函数后，矩阵$K$可以由原始空间的输入样本$X \in \Bbb{R}^{n\times d}$直接计算得到。

对于PCA，计算过程依赖于样本的协方差矩阵$S=X^TX \in \Bbb{R}^{d\times d}$，协方差矩阵衡量样本的不同特征之间的相似度。对于KPCA，样本映射到高维空间中的协方差矩阵为$S=\phi(X)^T\phi(X) \in \Bbb{R}^{d’\times d’}$，由于映射$\phi$通常是未知的，因此协方差矩阵无法直接计算。

另一方面，注意到矩阵$K=\phi(X)\phi(X)^T$已经得到，该矩阵衡量不同样本之间的相似度。下面寻找矩阵$S$和矩阵$K$对应的特征向量和特征值之间的关系。假设矩阵$K$具有特征向量$u\in \Bbb{R}^{n}$和特征值$\lambda$（矩阵$K$可以进行特征值分解），

上式两端左乘$\phi(X)^T$，

由上式可以看出，矩阵$S$具有特征向量$\phi(X)^Tu\in \Bbb{R}^{d’}$和特征值$\lambda$。

由于PCA要求矩阵$S$的特征向量是归一化的，而$\phi(X)^Tu$不满足；因此对其归一化：

对于KPCA，降维结果可以由矩阵$S$的前$k$个最大特征值对应的特征向量$v \in \Bbb{R}^{d’\times k}$表示为：

上式表明，求解KPCA的过程可以避开对协方差矩阵$S$进行特征值分解，只需要对矩阵$K$进行特征值分解即可。

值得一提的是，对于PCA通常对数据进行归一化处理，即减去所有样本数据的平均值。对于KPCA，无法显式地对高维映射后的数据$\phi(X)$进行处理，因此选择对矩阵$K=\phi(X)\phi(X)^T$进行预处理。

记单位矩阵$I_n \in \Bbb{R}^{n \times n}$和全$1$矩阵$\Bbb{1}_n \in \Bbb{R}^{n \times 1}$，则中心矩阵(centering matrix) $H=I_n - \frac{1}{N} \Bbb{1}_n \Bbb{1}_n^T$可以表示样本矩阵$\phi(X)$的中心化过程：

则中心化的矩阵$K$表示为：

综上所述，KPCA的一般步骤如下：

1. 给定输入样本$X \in \Bbb{R}^{n\times d}$和降维维度$k$;
2. 选定核函数，计算核矩阵$K \in \Bbb{R}^{n\times n}$；
3. 对矩阵$K$进行中心化：$\tilde{K}=(I_n - \frac{1}{N} \Bbb{1}_n \Bbb{1}_n^T)K (I_n - \frac{1}{N} \Bbb{1}_n \Bbb{1}_n^T)^T$;
4. 对矩阵$\tilde{K}$进行特征值分解；
5. 选择前$k$个最大特征值$\lambda$对应的特征向量$u \in \Bbb{R}^{n\times k}$;
6. 降维：$Y = \frac{1}{\sqrt{λ}}Ku$。

2. 根据表示定理求解KPCA

KPCA可以表示为高维特征空间的线性模型$y=Wz \in \Bbb{R}^{k}$。根据表示定理，线性模型的参数$W \in \Bbb{R}^{k \times d’}$的最优解可以表示为所有样本的线性组合：

则其第$j$维的投影向量$w_j \in \Bbb{R}^{d’}$可以表示成：

其中$\alpha^j=(\alpha_1^j,\alpha_2^j,…,\alpha_N^j) \in \Bbb{R}^N$。

根据主成分分析的结果，投影向量$w_j$应满足：

由于$\sum_{i=1}^{N}\phi(x_i)\phi(x_i)^T=\sum_{i=1}^{N}z_iz_i^T=ZZ^T$及$w_j =Z\alpha^j$，上式可以表示为：

一般情况下，我们不清楚映射$\phi$的具体形式，因此引入核函数：

则原式进一步表示为：

若下式满足，则上式自动满足：

即$\alpha^j$是$K$的特征向量。取$K$最大的$k$个特征向量${\alpha^1,\alpha^2,…,\alpha^k}$，则原始样本$x$投影后的第$j$维坐标$y_j$为：

3. 实现KPCA

① KPCA from scratch

由上述介绍，KPCA的一般步骤如下：

1. 给定输入样本$X \in \Bbb{R}^{n\times d}$和降维维度$k$;
2. 选定核函数，计算核矩阵$K \in \Bbb{R}^{n\times n}$；
3. 对矩阵$K$进行中心化：$\tilde{K}=(I_n - \frac{1}{N} \Bbb{1}_n \Bbb{1}_n^T)K (I_n - \frac{1}{N} \Bbb{1}_n \Bbb{1}_n^T)^T$;
4. 对矩阵$\tilde{K}$进行特征值分解；
5. 选择前$k$个最大特征值$\lambda$对应的特征向量$u \in \Bbb{R}^{n\times k}$;
6. 降维：$Y = \frac{1}{\sqrt{λ}}Ku$。

def kernel(xi, xj, gamma=1):
    return np.exp(-gamma**2*np.sum((xi-xj)**2))

def KPCA(data, k):
    n, d = data.shape
    K = np.zeros([n, n])
    # 计算核函数
    for i in range(n):
        for j in range(n):
            K[i, j] = kernel(data[i], data[j])
    # K矩阵归一化
    H = np.eye(n)-np.ones([n,n])/n # (n, n)
    K = H.dot(K).dot(H.T) # (n, n)
    # 计算特征值和特征向量
    eig_values, eig_vectors = np.linalg.eig(K)
    # 选择前k个最大的特征值标号
    index = np.argsort(-eig_values)[:k]
    # 选择对应的特征向量(主方向)
    KPCA_vectors = eig_vectors[index, :] # (k, n)
    # 降维
    KPCA_vectors /= np.expand_dims(np.sqrt(eig_values[index]),axis=-1)
    reduced_data = np.dot(K, KPCA_vectors.T) # (n, k)
    return reduced_data

② KPCA from sklearn

KPCA也可以通过sklearn库快速实现：

from sklearn.datasets import load_digits
from sklearn.decomposition import KernelPCA
import matplotlib.pyplot as plt

digits = load_digits()
X_kpca = KernelPCA(n_components=2, kernel="rbf", gamma=1).fit_transform(digits.data)

plt.figure(figsize=(10, 5))
plt.scatter(X_kpca[:, 0], X_kpca[:, 1], c=digits.target,label="KernelPCA")
plt.legend()
plt.show()