主要思想

PCA将$D$维特征$\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N]\in {\mathbb{R}}^{D\times N}$（${\mathbf{x}}_i\in {\mathbb{R}}^D$）映射到$d(d\ll D)$维空间中（$\mathbf{Y}=[{\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_N]={\mathbf{W}}^T\mathbf{X}\in {\mathbb{R}}^{d\times N}$, $\mathbf{W}=[{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d]\in {\mathbb{R}}^{D\times d}$，${\mathbf{w}}_i\in {\mathbb{R}}^D$），使得在降维后的空间中，特征的方差最大，即保留主成分。

推导方法

根据方差最大，可以确定优化目标为
arg max ⁡ w 1 , w 2 , ⋯   , w d 1 N ∑ j = 1 d ∑ i = 1 N ( w j T x i − w j T x ‾ ) 2

w1​,w2​,⋯,wd​argmax​N1​j=1∑d​i=1∑N​(wjT​xi​−wjT​x)2​​
其中
$$\overline{\mathbf{x}}=\frac{1}{N}\sum _{i=1}^N{\mathbf{x}}_i$$
则，可进一步将优化目标化简为矩阵形式：
arg max ⁡ w 1 , w 2 , ⋯   , w d 1 N ∑ j = 1 d ∑ i = 1 N ( w j T x i − w j T x ‾ ) 2 arg max ⁡ w 1 , w 2 , ⋯   , w d ∑ j = 1 d w j T ( 1 N ∑ i = 1 N ( x i − x ‾ ) ( x i − x ‾ ) T ) w j

\argmaxw1,w2,⋯,wd1N∑j=1d∑i=1N(wjTxi−wjTx¯)2\argmaxw1,w2,⋯,wd∑j=1dwjT(1N∑i=1N(xi−x¯)(xi−x¯)T)wj" role="presentation" style="position: relative;">\argmaxw1,w2,⋯,wd\argmaxw1,w2,⋯,wd1N∑j=1d∑i=1N(wTjxi−wTjx¯¯¯)2∑j=1dwTj(1N∑i=1N(xi−x¯¯¯)(xi−x¯¯¯)T)wj$$\begin{array}{rl}{\text{argmax}}_{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}& \frac{1}{N}\sum _{j=1}^d\sum _{i=1}^N({\mathbf{w}}_j^T{\mathbf{x}}_i-{\mathbf{w}}_j^T\overline{\mathbf{x}}{)}^2\\ {\text{argmax}}_{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}& \sum _{j=1}^d{\mathbf{w}}_j^T\left(\frac{1}{N}\sum _{i=1}^N({\mathbf{x}}_i-\overline{\mathbf{x}})({\mathbf{x}}_i-\overline{\mathbf{x}}{)}^T\right){\mathbf{w}}_j\end{array}$$

w1​,w2​,⋯,wd​argmax​w1​,w2​,⋯,wd​argmax​​N1​j=1∑d​i=1∑N​(wjT​xi​−wjT​x)2j=1∑d​wjT​(N1​i=1∑N​(xi​−x)(xi​−x)T)wj​​​​
令 X ~ = [ x 1 − x ‾ , x 2 − x ‾ , ⋯   , x N − x ‾ ] ∈ R D × N \widetilde{\mathbf{X}}=

[x1−x¯,x2−x¯,⋯,xN−x¯]" role="presentation" style="position: relative;">[x1−x¯¯¯,x2−x¯¯¯,⋯,xN−x¯¯¯]$$\left[\begin{array}{c}{\mathbf{x}}_1-\overline{\mathbf{x}},{\mathbf{x}}_2-\overline{\mathbf{x}},\cdots ,{\mathbf{x}}_N-\overline{\mathbf{x}}\end{array}\right]$$

\in\mathbb{R}^{D\times N} X =[x1​−x,x2​−x,⋯,xN​−x​]∈RD×N，则
arg max ⁡ w 1 , w 2 , ⋯   , w d ∑ j = 1 d w j T ( 1 N ∑ i = 1 N ( x i − x ‾ ) ( x i − x ‾ ) T ) w j arg max ⁡ w 1 , w 2 , ⋯   , w d ∑ j = 1 d w j T X ~ X ~ T N w j arg max ⁡ w 1 , w 2 , ⋯   , w d t r a c e ( W T X ~ X ~ T N W )

\argmaxw1,w2,⋯,wd∑j=1dwjT(1N∑i=1N(xi−x¯)(xi−x¯)T)wj\argmaxw1,w2,⋯,wd∑j=1dwjTX~X~TNwj\argmaxw1,w2,⋯,wdtrace(WTX~X~TNW)" role="presentation" style="position: relative;">\argmaxw1,w2,⋯,wd\argmaxw1,w2,⋯,wd\argmaxw1,w2,⋯,wd∑j=1dwTj(1N∑i=1N(xi−x¯¯¯)(xi−x¯¯¯)T)wj∑j=1dwTjX˜X˜TNwjtrace(WTX˜X˜TNW)$$\begin{array}{rl}{\text{argmax}}_{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}& \sum _{j=1}^d{\mathbf{w}}_j^T\left(\frac{1}{N}\sum _{i=1}^N({\mathbf{x}}_i-\overline{\mathbf{x}})({\mathbf{x}}_i-\overline{\mathbf{x}}{)}^T\right){\mathbf{w}}_j\\ {\text{argmax}}_{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}& \sum _{j=1}^d{\mathbf{w}}_j^T\frac{\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T}{N}{\mathbf{w}}_j\\ {\text{argmax}}_{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}& trace({\mathbf{W}}^T\frac{\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T}{N}\mathbf{W})\end{array}$$

w1​,w2​,⋯,wd​argmax​w1​,w2​,⋯,wd​argmax​w1​,w2​,⋯,wd​argmax​​j=1∑d​wjT​(N1​i=1∑N​(xi​−x)(xi​−x)T)wj​j=1∑d​wjT​NX X T​wj​trace(WTNX X T​W)​​​
一般来说引入${\mathbf{W}}^T\mathbf{W}=\mathbb{I}$的约束，这样能确保只关注${\mathbf{w}}_i$的方向而不是尺度大小，并且对任意$i\ne j$${\mathbf{w}}_i$和${\mathbf{w}}_j$都是不相关的，从而能够保留更多的原数据信息。新的优化目标为
arg max ⁡ W    t r a c e ( W T X ~ X ~ T N W ) s . t . W T W = I

\argmaxW  trace(WTX~X~TNW)s.t.WTW=I" role="presentation" style="position: relative;">\argmaxW  s.t.trace(WTX˜X˜TNW)WTW=I$$\begin{array}{rl}{\text{argmax}}_{\mathbf{W}}& trace({\mathbf{W}}^T\frac{\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T}{N}\mathbf{W})\\ s.t.& {\mathbf{W}}^T\mathbf{W}=\mathbb{I}\end{array}$$

Wargmax​  s.t.​trace(WTNX X T​W)WTW=I​​​
采用拉格朗日乘子的方法求解以上优化问题，引入$\mathbf{\Lambda }\in {\mathbb{R}}^{d\times d}$，由于约束条件是对称的，因此$\mathbf{\Lambda }$也满足${\mathbf{\Lambda }}^T=\mathbf{\Lambda }$，则
$$L\left(\mathbf{W},\mathbf{\Lambda }\right)=-trace\left({\mathbf{W}}^T\frac{{\tilde{\mathbf{X}}\tilde{\mathbf{X}}}^T}{N}\mathbf{W}\right)+trace\left(\mathbf{\Lambda }\left({\mathbf{W}}^T\mathbf{W}-\mathbb{I}\right)\right)$$
则
∂ L ( W , Λ ) ∂ W = − 2 X ~ X ~ T N W + 2 W Λ = 0 X ~ X ~ T N W = W Λ [ X ~ X ~ T N w 1 , ⋯   , X ~ X ~ T N w d ] = [ w 1 , ⋯   , w d ] [ λ 11 λ 12 ⋯ λ 1 d λ 21 λ 22 ⋯ λ 2 d ⋮ ⋮ ⋱ ⋮ λ d 1 λ d 2 ⋯ λ d d ] \frac{\partial L\left( \mathbf{W},\mathbf{\Lambda } \right)}{\partial \mathbf{W}}=-2\frac{\mathbf{\tilde{X}\tilde{X}}^T}{N}\mathbf{W}+2\mathbf{ W\Lambda}=0 \\ \frac{\mathbf{\tilde{X}\tilde{X}}^T}{N}\mathbf{W}=\mathbf{W\Lambda } \\ \left[ \frac{\mathbf{\tilde{X}\tilde{X}}^T}{N}\mathbf{w}_1,\cdots ,\frac{\mathbf{\tilde{X}\tilde{X}}^T}{N}\mathbf{w}_d \right] =\left[ \mathbf{w}_1,\cdots ,\mathbf{w}_d \right] \left[

λ11λ12⋯λ1dλ21λ22⋯λ2d⋮⋮⋱⋮λd1λd2⋯λdd" role="presentation" style="position: relative;">λ11λ21⋮λd1λ12λ22⋮λd2⋯⋯⋱⋯λ1dλ2d⋮λdd$$\begin{array}{cccc}{\lambda }_{11}& {\lambda }_{12}& \cdots & {\lambda }_{1d}\\ {\lambda }_{21}& {\lambda }_{22}& \cdots & {\lambda }_{2d}\\ ⋮& ⋮& \ddots & ⋮\\ {\lambda }_{d1}& {\lambda }_{d2}& \cdots & {\lambda }_{dd}\end{array}$$

\right] \\ ∂W∂L(W,Λ)​=−2NX~X~T​W+2WΛ=0NX~X~T​W=WΛ[NX~X~T​w1​,⋯,NX~X~T​wd​]=[w1​,⋯,wd​]⎣ ⎡​λ11​λ21​⋮λd1​​λ12​λ22​⋮λd2​​⋯⋯⋱⋯​λ1d​λ2d​⋮λdd​​⎦ ⎤​
X ~ X ~ T N w i = [ w 1 , ⋯   , w d ] [ λ 1 i ⋮ λ d i ] w i T X ~ X ~ T N w i = λ i i \frac{\mathbf{\tilde{X}\tilde{X}}^T}{N}\mathbf{w}_i=\left[ \mathbf{w}_1,\cdots ,\mathbf{w}_d \right] \left[

λ1i⋮λdi" role="presentation" style="position: relative;">λ1i⋮λdi$$\begin{array}{c}{\lambda }_{1i}\\ ⋮\\ {\lambda }_{di}\end{array}$$

\right] \\ \mathbf{w}_i^T\frac{\mathbf{\tilde{X}\tilde{X}}^T}{N}\mathbf{w}_i=\lambda _{ii} NX~X~T​wi​=[w1​,⋯,wd​]⎣ ⎡​λ1i​⋮λdi​​⎦ ⎤​wiT​NX~X~T​wi​=λii​
以上优化问题的最终结果，实质上就是$\frac{\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T}{N}$的最大的$d$特征值（${\lambda }_{11},{\lambda }_{22},\cdots ,{\lambda }_{dd}$）所对应的特征向量。

import matplotlib.pyplot as plt
import numpy as np

from sklearn import datasets
from sklearn.decomposition import PCA

def plot_2d(X_r, y, target_names, name):
    plt.subplot(1,2,1 if 'Sklearn' in name else 2)
    colors = ["navy", "turquoise", "darkorange"]
    lw = 2

    for color, i, target_name in zip(colors, [0, 1, 2], target_names):
        plt.scatter(
            X_r[y == i, 0], X_r[y == i, 1], color=color, alpha=0.8, lw=lw, label=target_name
        )
    plt.legend(loc="best", shadow=False, scatterpoints=1)
    plt.xlabel(f"{name} of IRIS dataset")

def sklearn_pca(X, y, target_names):
    pca = PCA(n_components=2)
    X_r = pca.fit(X).transform(X)

    print(pca.explained_variance_)

    plot_2d(X_r, y, target_names, 'Sklearn PCA')

    return X_r

def my_pca(X, y, target_names):
    n_components=2
    # 去中心化
    N = X.shape[0]
    X_mean = np.mean(X, axis=0)
    X_ = X - X_mean
    # 构建协方差矩阵
    XX = 1. / N * np.matmul(X_.T, X_)
    # 求特征向量
    values, vectors = np.linalg.eig(XX)
    idxs = np.argsort(values)
    
    W = vectors[:,[idxs[-i] for i in range(1,n_components+1)]]
    
    X_r = np.matmul(X_, W)

    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致

    plot_2d(X_r, y, target_names, 'My Imple. PCA')

    return X_r



if __name__ == '__main__':
    iris = datasets.load_iris()

    X = iris.data
    y = iris.target
    target_names = iris.target_names

    plt.figure()

    X_r1 = sklearn_pca(X, y, target_names)

    X_r2 = my_pca(X, y, target_names)

    print(np.sum(np.abs(X_r1 - X_r2)))

    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

说明自己实现的方法与sklearn内核一致。

核技巧

利用$\varphi$将${\mathbf{x}}_i\in {\mathbb{R}}^D$映射到高维度空间，得到${\mathbf{X}}_{\varphi }=[\varphi ({\mathbf{x}}_1),\varphi ({\mathbf{x}}_2),\cdots ,\varphi ({\mathbf{x}}_N)]$，在该高维空间中进行PCA操作。首先，去中心化
$$\widetilde{{\mathbf{X}}_{\varphi }}=\left[\varphi ({\mathbf{x}}_1)-\overline{\varphi (\mathbf{x})},\varphi ({\mathbf{x}}_2)-\overline{\varphi (\mathbf{x})},\cdots ,\varphi ({\mathbf{x}}_N)-\overline{\varphi (\mathbf{x})}\right]$$
其中
$$\overline{\varphi (\mathbf{x})}=\frac{1}{N}\sum _{i=1}^N\varphi ({\mathbf{x}}_i)$$
则
X ϕ ~ X ϕ ~ T N = 1 N ∑ i = 1 N ( ϕ ( x i ) − 1 N ∑ j = 1 N ϕ ( x j ) ) ( ϕ ( x i ) − 1 N ∑ j = 1 N ϕ ( x j ) ) T

NXϕ​ ​Xϕ​ ​T​​=N1​i=1∑N​(ϕ(xi​)−N1​j=1∑N​ϕ(xj​))(ϕ(xi​)−N1​j=1∑N​ϕ(xj​))T​
需计算以上矩阵的特征向量${\mathbf{w}}_k$（对应第$k$大的特征值为${\lambda }_k$所对应的特征向量）
1 N ∑ i = 1 N ( ϕ ( x i ) − 1 N ∑ j = 1 N ϕ ( x j ) ) ( ϕ ( x i ) − 1 N ∑ j = 1 N ϕ ( x j ) ) T w k = λ k w k

1N∑i=1N(ϕ(xi)−1N∑j=1Nϕ(xj))(ϕ(xi)−1N∑j=1Nϕ(xj))Twk=λkwk" role="presentation" style="position: relative;">1N∑i=1N(ϕ(xi)−1N∑j=1Nϕ(xj))(ϕ(xi)−1N∑j=1Nϕ(xj))Twk=λkwk$$\begin{array}{rl}\frac{1}{N}\sum _{i=1}^N(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j))(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j){)}^T{\mathbf{w}}_k& ={\lambda }_k{\mathbf{w}}_k\end{array}$$

N1​i=1∑N​(ϕ(xi​)−N1​j=1∑N​ϕ(xj​))(ϕ(xi​)−N1​j=1∑N​ϕ(xj​))Twk​​=λk​wk​​​​
所以
w k = [ ∑ i = 1 N ϕ ( x i ) T w k − 1 N ∑ j = 1 N ϕ ( x j ) T w k N λ k ϕ ( x i ) ] = ∑ i = 1 N α k i ϕ ( x i )

\begin{aligned} \mathbf{w}_k&=\left[ \sum_{i=1}^N{ {\color{red} \frac{\phi (\mathbf{x}_i)^T\mathbf{w}_k-\frac{1}{N}\sum_{j=1}^N{\phi (\mathbf{x}_j)^T\mathbf{w}_k}}{N\lambda _k}}\phi (\mathbf{x}_i)} \right]\\ &=\sum_{i=1}^N{ {\color{red} \alpha _{ki}}\phi (\mathbf{x}_i)}\\ \end{aligned}" role="presentation" style="position: relative;">\begin{aligned} \mathbf{w}_k&=\left[ \sum_{i=1}^N{ {\color{red} \frac{\phi (\mathbf{x}_i)^T\mathbf{w}_k-\frac{1}{N}\sum_{j=1}^N{\phi (\mathbf{x}_j)^T\mathbf{w}_k}}{N\lambda _k}}\phi (\mathbf{x}_i)} \right]\\ &=\sum_{i=1}^N{ {\color{red} \alpha _{ki}}\phi (\mathbf{x}_i)}\\ \end{aligned}$$\text{begin{aligned} mathbf{w}\_k\&=left[ sum\_{i=1}^N{ {color{red} frac{phi (mathbf{x}\_i)^Tmathbf{w}\_k-frac{1}{N}sum\_{j=1}^N{phi (mathbf{x}\_j)^Tmathbf{w}\_k}}{Nlambda \_k}}phi (mathbf{x}\_i)} right] \&=sum\_{i=1}^N{ {color{red} alpha \_{ki}}phi (mathbf{x}\_i)} end{aligned}}$$

wk​​=[i=1∑N​Nλk​ϕ(xi​)Twk​−N1​∑j=1N​ϕ(xj​)Twk​​ϕ(xi​)]=i=1∑N​αki​ϕ(xi​)​
其中${\alpha }_{ki}=\left(\varphi ({\mathbf{x}}_i{)}^T{\mathbf{w}}_k-\frac{1}{N}{\sum }_{j=1}^N\varphi ({\mathbf{x}}_j{)}^T{\mathbf{w}}_k\right)/(N{\lambda }_k)$，则下面将求出相应的${\alpha }_{ki}$，将上式带入公式（5），可得
1 N ∑ i = 1 N ( ϕ ( x i ) − 1 N ∑ j = 1 N ϕ ( x j ) ) ( ϕ ( x i ) − 1 N ∑ j = 1 N ϕ ( x j ) ) T ∑ i = 1 N α k i ϕ ( x i ) = λ k ∑ i = 1 N α k i ϕ ( x i ) 1 N ( X ϕ X ϕ T − 1 N X ϕ 1 N × 1 1 N × 1 T X ϕ T ) X ϕ α k = λ k X ϕ α k 1 N ( X ϕ T X ϕ X ϕ T X ϕ α k − X ϕ T X ϕ 1 N × 1 1 N × 1 T N X ϕ T X ϕ α k ) = λ k X ϕ T X ϕ α k

1N∑i=1N(ϕ(xi)−1N∑j=1Nϕ(xj))(ϕ(xi)−1N∑j=1Nϕ(xj))T∑i=1Nαkiϕ(xi)=λk∑i=1Nαkiϕ(xi)1N(XϕXϕT−1NXϕ1N×11N×1TXϕT)Xϕαk=λkXϕαk1N(XϕTXϕXϕTXϕαk−XϕTXϕ1N×11N×1TNXϕTXϕαk)=λkXϕTXϕαk" role="presentation" style="position: relative;">1N∑i=1N(ϕ(xi)−1N∑j=1Nϕ(xj))(ϕ(xi)−1N∑j=1Nϕ(xj))T∑i=1Nαkiϕ(xi)1N(XϕXTϕ−1NXϕ1N×11TN×1XTϕ)Xϕαk1N(XTϕXϕXTϕXϕαk−XTϕXϕ1N×11TN×1NXTϕXϕαk)=λk∑i=1Nαkiϕ(xi)=λkXϕαk=λkXTϕXϕαk$$\begin{array}{rl}\frac{1}{N}\sum _{i=1}^N(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j))(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j){)}^T\sum _{i=1}^N{\alpha }_{ki}\varphi ({\mathbf{x}}_i)& ={\lambda }_k\sum _{i=1}^N{\alpha }_{ki}\varphi ({\mathbf{x}}_i)\\ \frac{1}{N}({\mathbf{X}}_{\varphi }{\mathbf{X}}_{\varphi }^T-\frac{1}{N}{\mathbf{X}}_{\varphi }{1}_{N\times 1}{1}_{N\times 1}^T{\mathbf{X}}_{\varphi }^T){\mathbf{X}}_{\varphi }{\alpha }_k& ={\lambda }_k{\mathbf{X}}_{\varphi }{\alpha }_k\\ \frac{1}{N}({\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k-{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }\frac{1_{N\times 1}{1}_{N\times 1}^T}{N}{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k)& ={\lambda }_k{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k\end{array}$$

N1​i=1∑N​(ϕ(xi​)−N1​j=1∑N​ϕ(xj​))(ϕ(xi​)−N1​j=1∑N​ϕ(xj​))Ti=1∑N​αki​ϕ(xi​)N1​(Xϕ​XϕT​−N1​Xϕ​1N×1​1N×1T​XϕT​)Xϕ​αk​N1​(XϕT​Xϕ​XϕT​Xϕ​αk​−XϕT​Xϕ​N1N×1​1N×1T​​XϕT​Xϕ​αk​)​=λk​i=1∑N​αki​ϕ(xi​)=λk​Xϕ​αk​=λk​XϕT​Xϕ​αk​​​​
而核矩阵$\mathbf{K}={\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }$，所以
1 N ( K K α k − K 1 N × 1 1 N × 1 T N K α k ) = λ k K α k 1 N ( I − 1 N × N N ) K α k = λ k α k K ~ α k = λ k α k

1N(KKαk−K1N×11N×1TNKαk)=λkKαk1N(I−1N×NN)Kαk=λkαkK~αk=λkαk" role="presentation" style="position: relative;">1N(KKαk−K1N×11TN×1NKαk)1N(I−1N×NN)KαkK˜αk=λkKαk=λkαk=λkαk$$\begin{array}{rl}\frac{1}{N}(\mathbf{K}\mathbf{K}{\alpha }_k-\mathbf{K}\frac{1_{N\times 1}{1}_{N\times 1}^T}{N}\mathbf{K}{\alpha }_k)& ={\lambda }_k\mathbf{K}{\alpha }_k\\ \frac{1}{N}\left(\mathbf{I}-\frac{1_{N\times N}}{N}\right)\mathbf{K}{\alpha }_k& ={\lambda }_k{\alpha }_k\\ \tilde{\mathbf{K}}{\alpha }_k& ={\lambda }_k{\alpha }_k\end{array}$$

N1​(KKαk​−KN1N×1​1N×1T​​Kαk​)N1​(I−N1N×N​​)Kαk​K αk​​=λk​Kαk​=λk​αk​=λk​αk​​
那么可以求出$\tilde{\mathbf{K}}$的第$k$大的特征值${\lambda }_k$所对应的特征向量${\alpha }_k$，从而得到
$${\mathbf{w}}_k={\mathbf{X}}_{\varphi }{\alpha }_k$$
对一个新的数据${\mathbf{x}}_{new}$，降维到${\mathbf{w}}_k$的维度
w k T ( ϕ ( x n e w ) − 1 N X ϕ 1 N × 1 ) = α k T X ϕ T ϕ ( x n e w ) − 1 N α k T X ϕ T X ϕ 1 N × 1 = α k T K ( ⋅ , x n e w ) − 1 N α k T K 1 N × 1

\begin{aligned} {\color{red} \mathbf{w}_{k}^{T}}\left( \phi \left( \mathbf{x}_{new} \right) -\frac{1}{N}\mathbf{X}_{\phi}1_{N\times 1} \right) &={\color{red} \alpha _{k}^{T}\mathbf{X}_{\phi}^{T}}\phi \left( \mathbf{x}_{new} \right) -\frac{1}{N}{\color{red} \alpha _{k}^{T}\mathbf{X}_{\phi}^{T}}\mathbf{X}_{\phi}1_{N\times 1} \\ &=\alpha _{k}^{T}\mathbf{K}\left( \cdot ,\mathbf{x}_{new} \right) -\frac{1}{N}\alpha _{k}^{T}\mathbf{K}1_{N\times 1} \end{aligned}" role="presentation" style="position: relative;">\begin{aligned} {\color{red} \mathbf{w}_{k}^{T}}\left( \phi \left( \mathbf{x}_{new} \right) -\frac{1}{N}\mathbf{X}_{\phi}1_{N\times 1} \right) &={\color{red} \alpha _{k}^{T}\mathbf{X}_{\phi}^{T}}\phi \left( \mathbf{x}_{new} \right) -\frac{1}{N}{\color{red} \alpha _{k}^{T}\mathbf{X}_{\phi}^{T}}\mathbf{X}_{\phi}1_{N\times 1} \\ &=\alpha _{k}^{T}\mathbf{K}\left( \cdot ,\mathbf{x}_{new} \right) -\frac{1}{N}\alpha _{k}^{T}\mathbf{K}1_{N\times 1} \end{aligned}$$\text{begin{aligned} {color{red} mathbf{w}\_k^T}left( phi left( mathbf{x}\_{new} right) -frac{1}{N}mathbf{X}\_{phi}1\_{Ntimes 1} right) \&={color{red} alpha \_k^{T}mathbf{X}\_{phi}^T}phi left( mathbf{x}\_{new} right) -frac{1}{N}{color{red} alpha \_k^{T}mathbf{X}\_{phi}^T}mathbf{X}\_{phi}1\_{Ntimes 1} \&=alpha \_k^{T}mathbf{K}left( cdot ,mathbf{x}\_{new} right) -frac{1}{N}alpha \_k^{T}mathbf{K}1\_{Ntimes 1} end{aligned}}$$

wkT​(ϕ(xnew​)−N1​Xϕ​1N×1​)​=αkT​XϕT​ϕ(xnew​)−N1​αkT​XϕT​Xϕ​1N×1​=αkT​K(⋅,xnew​)−N1​αkT​K1N×1​​

import matplotlib.pyplot as plt
import numpy as np

from sklearn import datasets
from sklearn.decomposition import PCA, KernelPCA
from sklearn.datasets import make_circles
from sklearn.model_selection import train_test_split

def plot_2d(X, y, X_pca, X_kpca, X_my_kpca):
    fig, (orig_data_ax, pca_proj_ax, kernel_pca_proj_ax, kernel_pca_proj_ax2) = plt.subplots(
    ncols=4, figsize=(17, 4)
    )

    orig_data_ax.scatter(X[:, 0], X[:, 1], c=y)
    orig_data_ax.set_ylabel("Feature #1")
    orig_data_ax.set_xlabel("Feature #0")
    orig_data_ax.set_title("Testing data")

    pca_proj_ax.scatter(X_pca[:, 0], X_pca[:, 1], c=y)
    pca_proj_ax.set_ylabel("Principal component #1")
    pca_proj_ax.set_xlabel("Principal component #0")
    pca_proj_ax.set_title("my Imple. PCA")

    kernel_pca_proj_ax.scatter(X_kpca[:, 0], X_kpca[:, 1], c=y)
    kernel_pca_proj_ax.set_ylabel("Principal component #1")
    kernel_pca_proj_ax.set_xlabel("Principal component #0")
    _ = kernel_pca_proj_ax.set_title("Sklearn KernelPCA")

    kernel_pca_proj_ax2.scatter(X_my_kpca[:, 0], X_my_kpca[:, 1], c=y)
    kernel_pca_proj_ax2.set_ylabel("Principal component #1")
    kernel_pca_proj_ax2.set_xlabel("Principal component #0")
    _ = kernel_pca_proj_ax2.set_title("my Imple. KernelPCA")

    plt.show()

def my_pca(X_tr, X_te):
    n_components=2
    # 去中心化
    N = X_tr.shape[0]
    X_mean = np.mean(X_tr, axis=0)
    X_ = X_tr - X_mean
    # 构建协方差矩阵
    XX = 1. / N * np.matmul(X_.T, X_)
    # 求特征向量
    values, vectors = np.linalg.eig(XX)
    idxs = np.argsort(values)
    
    W = vectors[:,[idxs[-i] for i in range(1,n_components+1)]]
    
    X_r = np.matmul(X_te - X_mean, W)

    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致

    return X_r

def my_kpca(X_tr, X_te, gamma=10):
    # exp(-||x_i-x_j||^2*gamma)
    N = X_tr.shape[0]
    K_tr0 = np.sum(X_tr**2, axis=1, keepdims=True)
    K_tr0 = np.exp(-(K_tr0 + K_tr0.T - 2 * np.matmul(X_tr, X_tr.T))*gamma)
    K_ = 1. / N * (K_tr0 - np.ones([N, N]) @ K_tr0 /N)
    # 求特征向量
    values, vectors = np.linalg.eig(K_)
    idxs = np.argsort(values)
    
    alphas = vectors[:,[idxs[-i] for i in range(1,3)]] # N_tr x 2
    K_new = np.sum(X_tr**2, axis=1, keepdims=True) + \
            np.sum(X_te**2, axis=1, keepdims=True).T - \
            2 * X_tr @ X_te.T
    K_new = np.exp(-K_new * gamma) # N_tr x N_te
    X_r = K_new.T @ alphas - 1. / N * np.ones([1, N]) @ K_tr0 @ alphas

    X_r[:,0] = - X_r[:,0] #确保与sklearn得到的结果一致
    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致

    return X_r

def sklearn_kpca(X_tr, X_te, gamma=10):
    kernel_pca = KernelPCA(
        n_components=None, kernel="rbf", gamma=gamma, fit_inverse_transform=False
    )
    X_test_kernel_pca = kernel_pca.fit(X_train).transform(X_test)
    return X_test_kernel_pca



if __name__ == '__main__':
    X, y = make_circles(n_samples=1000, factor=0.3, noise=0.05, random_state=0)
    X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=0)

    X_pca = my_pca(X_train, X_test)
    X_kpca = sklearn_kpca(X_train, X_test, gamma=10.0)
    X_my_kpca = my_kpca(X_train, X_test, gamma=10.0)

    plot_2d(X_test, y_test, X_pca, X_kpca, X_my_kpca)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

虽然与sklearn的KPCA在倍数上有点出入，但是不想深究了，就这样吧！