机器学习 07：PCA 及其 sklearn 源码解读_sklearn源码分析

概述

---

Principal Component Analysis 即主成分分析的主要思想是将 n 维特征映射到 k 维上，这 k 维是全新的正交特征也被称为主成分，是在原有 n 维特征的基础上重新构造出来的 k 维特征。PCA 的工作就是从原始的空间中顺序地找一组相互正交的坐标轴，其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第 1,2 个轴正交的平面中方差最大的，依次类推，得到 n 个这样的坐标轴。通过这种方式获得的新的坐标轴，大部分方差都包含在前面 k 个坐标轴中，后面的坐标轴所含的方差几乎为 0 。于是，我们可以只保留前 k 维，从而实现对数据特征的降维处理。

PCA 的计算有两种方式：

1. 分解协方差矩阵：具体实现过程主要包括三步：
   1. 计算数据矩阵的协方差矩阵
   2. 计算协方差矩阵的特征值特征向量
   3. 选择特征值最大（即方差最大）的 k 个特征所对应的特征向量组成的矩阵
2. 分解内积矩阵

分解协方差矩阵

---

协方差矩阵

协方差用来刻画两个随机变量之间的相关性。简单比较一下方差与协方差：

均值：$$\bar{x}=\frac{1}{n}\sum _{i=1}^N{\vec{x}}_i$$

样本无偏方差：$$S^2=\frac{1}{n-1}\sum _{i=1}^N(\bar{x}-{\vec{x}}_i{)}^2$$

样本 $X$ 和样本 $Y$ 的协方差：$$Cov(X,Y)=\frac{1}{n-1}\sum _{i=1}^N(\bar{x}-{\vec{x}}_i)(\bar{y}-{\vec{y}}_i)$$

我们可以发现：方差的计算是针对一维特征，而协方差则要求至少二维。方差是协方差的特殊情况。协方差为正时，说明 $X,Y$ 是正相关关系，为负时，则说明是负相关关系，为 0 则说明两个样本相互独立。当有多个样本时，两两之间的协方差就构成了协方差矩阵。

协方差矩阵与散度矩阵关系密切，由散度矩阵的定义：$$S=\sum _{i=1}^N({\vec{x}}_i-\bar{x})({\vec{x}}_i-\bar{x}{)}^T$$我们发现，对于一组数据 $X$： $S=(n-1)Cov(X)$，因此它们有相同的特征值和特征向量。

特征值分解和 SVD

详细过程可以参考我的另一篇文章：特征分解和奇异值分解
计算出 n 个特征向量后，我们对其排序，并取前 k 个即可。

分解内积矩阵

---

在 sklearn 库中的 PCA 就是使用此方法计算的。

假设 $X$ 矩阵是一个样本集，它已经中心化过（PCA 要求数据中心化），那么协方差矩阵就是：$\frac{1}{n}X{X}^T$
我们对 $X$ 进行 SVD 分解：$X=U\Sigma V^T$，则 $X^T=V\Sigma U^T$
因此就有：$\frac{1}{n}X{X}^T=U{\Sigma }^2{U}^T$

因此，我们可以对 $X$ 进行 SVD 分解，然后求 $\Sigma$ 的内积来计算特征向量。

关于源码

---

花一晚上阅读了 sklearn 的官方代码，大概记录一下。源码篇幅过长，所以仅记录重点部分。

计算 PCA

通过 svd_solver 参数，我们可以选择不同方式对数据进行求解，但是基本思路相同，这里仅以 svd_solver == 'full' 也就是使用所有数据进行特征分解这一情况为例。

fit() 方法的实现：
其中，explained_variance_ 为方差值，explained_variance_ratio_ 为每个方程的贡献率，也就是占比，可以通过这个参数来判断 n 具体取值。

注意，n_components 参数可以取好几种值：

1. “mle”： 自动确定 n 值。
2. $0<\text{n\_components}<1$：表示希望保留的信息量的比例。PCA 会自动选使得 信息量 $\ge \text{n\_components}$ 的特征数量。
3. $\text{n\_components}\ge 1$ 且为整数：选前 n 个

def _fit_full(self, X, n_components):
	n_samples, n_features = X.shape
	
	... # 对传入的数据进行验证，暂且省略
	
	# Center data
	self.mean_ = np.mean(X, axis=0)
	X -= self.mean_
	
	# 使用分解内积矩阵方法求特征向量
	U, S, Vt = linalg.svd(X, full_matrices=False)
	U, Vt = svd_flip(U, Vt) # 由于有时候可能求得负值，所以翻转特征向量符号以强制执行确定性输出
	components_ = Vt
	
	explained_variance_ = (S**2) / (n_samples - 1)
	total_var = explained_variance_.sum()
	explained_variance_ratio_ = explained_variance_ / total_var
	singular_values_ = S.copy()  # Store the singular values.
	
	# 根据 n_components 确定取几个特征，也就是 n 的值。
	if n_components == "mle":
	    n_components = _infer_dimension(explained_variance_, n_samples)
	elif 0 < n_components < 1.0:
	    ratio_cumsum = stable_cumsum(explained_variance_ratio_)
	    n_components = np.searchsorted(ratio_cumsum, n_components, side="right") + 1
	# Compute noise covariance using Probabilistic PCA model
	# The sigma2 maximum likelihood (cf. eq. 12.46)
	if n_components < min(n_features, n_samples):
	    self.noise_variance_ = explained_variance_[n_components:].mean()
	else:
	    self.noise_variance_ = 0.0
	
	... # 对参数进行赋值，暂且省略
	
	return U, S, Vt
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

transform() 方法的实现：

这里运用了白化操作，即对数据进行去相关，以降低输入数据的冗余性，白化后的数据有两个特性：消除了特征之间的相关性，且所有特征的方差都为 1。

def transform(self, X):
        ... # 对程序和数据进行验证，暂且忽略
        
        if self.mean_ is not None:
            X = X - self.mean_
        X_transformed = np.dot(X, self.components_.T)
        # 白化
        if self.whiten:
            X_transformed /= np.sqrt(self.explained_variance_)
        return X_transformed
1
2
3
4
5
6
7
8
9
10

fit_transform() 方法的实现：

相关的公式都被作者贴心的写在注释里了。

def fit_transform(self, X, y=None):
        U, S, Vt = self._fit(X)
        U = U[:, : self.n_components_]

        if self.whiten:
            # X_new = X * V / S * sqrt(n_samples) = U * sqrt(n_samples)
            U *= sqrt(X.shape[0] - 1)
        else:
            # X_new = X * V = U * S * Vt * V = U * S
            U *= S[: self.n_components_]
        return U
1
2
3
4
5
6
7
8
9
10
11

自动选择 n 的值

通过计算对数似然的方法来判断 n 具体该取几，也就是保留几个维度。相关论文可见此链接

def _assess_dimension(spectrum, rank, n_samples):
    """Compute the log-likelihood of a rank ``rank`` dataset.
    The dataset is assumed to be embedded in gaussian noise of shape(n,
    dimf) having spectrum ``spectrum``.
    Parameters
    ----------
    spectrum : ndarray of shape (n_features,)
        Data spectrum.
    rank : int
        Tested rank value. It should be strictly lower than n_features,
        otherwise the method isn't specified (division by zero in equation
        (31) from the paper).
    n_samples : int
        Number of samples.
    Returns
    -------
    ll : float
        The log-likelihood.
    """

    n_features = spectrum.shape[0]
    if not 1 <= rank < n_features:
        raise ValueError("the tested rank should be in [1, n_features - 1]")

    eps = 1e-15

    if spectrum[rank - 1] < eps:
        # When the tested rank is associated with a small eigenvalue, there's
        # no point in computing the log-likelihood: it's going to be very
        # small and won't be the max anyway. Also, it can lead to numerical
        # issues below when computing pa, in particular in log((spectrum[i] -
        # spectrum[j]) because this will take the log of something very small.
        return -np.inf

    pu = -rank * log(2.0)
    for i in range(1, rank + 1):
        pu += (
            gammaln((n_features - i + 1) / 2.0)
            - log(np.pi) * (n_features - i + 1) / 2.0
        )

    pl = np.sum(np.log(spectrum[:rank]))
    pl = -pl * n_samples / 2.0

    v = max(eps, np.sum(spectrum[rank:]) / (n_features - rank))
    pv = -np.log(v) * n_samples * (n_features - rank) / 2.0

    m = n_features * rank - rank * (rank + 1.0) / 2.0
    pp = log(2.0 * np.pi) * (m + rank) / 2.0

    pa = 0.0
    spectrum_ = spectrum.copy()
    spectrum_[rank:n_features] = v
    for i in range(rank):
        for j in range(i + 1, len(spectrum)):
            pa += log(
                (spectrum[i] - spectrum[j]) * (1.0 / spectrum_[j] - 1.0 / spectrum_[i])
            ) + log(n_samples)

    ll = pu + pl + pv + pp - pa / 2.0 - rank * log(n_samples) / 2.0
    return ll


def _infer_dimension(spectrum, n_samples):
    """Infers the dimension of a dataset with a given spectrum.
    The returned value will be in [1, n_features - 1].
    """
    ll = np.empty_like(spectrum)
    ll[0] = -np.inf  # we don't want to return n_components = 0
    for rank in range(1, spectrum.shape[0]):
        ll[rank] = _assess_dimension(spectrum, rank, n_samples)
    return ll.argmax()
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

样本得分

通过计算每个样本的对数似然来反应性能。相关论文可见此链接

def get_covariance(self):
    """Compute data covariance with the generative model.
    ``cov = components_.T * S**2 * components_ + sigma2 * eye(n_features)``
    where S**2 contains the explained variances, and sigma2 contains the
    noise variances.
    Returns
    -------
    cov : array of shape=(n_features, n_features)
        Estimated covariance of data.
    """
    components_ = self.components_
    exp_var = self.explained_variance_
    if self.whiten:
        components_ = components_ * np.sqrt(exp_var[:, np.newaxis])
    exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.0)
    cov = np.dot(components_.T * exp_var_diff, components_)
    cov.flat[:: len(cov) + 1] += self.noise_variance_  # modify diag inplace
    return cov

def get_precision(self):
    """Compute data precision matrix with the generative model.
    Equals the inverse of the covariance but computed with
    the matrix inversion lemma for efficiency.
    Returns
    -------
    precision : array, shape=(n_features, n_features)
        Estimated precision of data.
    """
    n_features = self.components_.shape[1]

    # handle corner cases first
    if self.n_components_ == 0:
        return np.eye(n_features) / self.noise_variance_

    if np.isclose(self.noise_variance_, 0.0, atol=0.0):
        return linalg.inv(self.get_covariance())

    # Get precision using matrix inversion lemma
    components_ = self.components_
    exp_var = self.explained_variance_
    if self.whiten:
        components_ = components_ * np.sqrt(exp_var[:, np.newaxis])
    exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.0)
    precision = np.dot(components_, components_.T) / self.noise_variance_
    precision.flat[:: len(precision) + 1] += 1.0 / exp_var_diff
    precision = np.dot(components_.T, np.dot(linalg.inv(precision), components_))
    precision /= -(self.noise_variance_**2)
    precision.flat[:: len(precision) + 1] += 1.0 / self.noise_variance_
    return precision

def score_samples(self, X):
	"""Return the log-likelihood of each sample.
	Parameters
	----------
	X : array-like of shape (n_samples, n_features)
	    The data.
	"""
	... # 对程序和数据的验证，暂且忽略
	
	Xr = X - self.mean_
	n_features = X.shape[1]
	precision = self.get_precision()
	log_like = -0.5 * (Xr * (np.dot(Xr, precision))).sum(axis=1)
	log_like -= 0.5 * (n_features * log(2.0 * np.pi) - fast_logdet(precision))
	return log_like

def score(self, X, y=None):
    """Return the average log-likelihood of all samples."""
    return np.mean(self.score_samples(X))
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

---

2023-8-2 更新

---

PCA 公式推导

---

首先需要两个前置知识：

谱分解定理 Spectral Theorem

Let $A\in R^n$ be symmetric, and ${\lambda }_i\in R,i\in [1,n]$ be the eigenvalues of $A$. There exists a set of orthonormal vectors ${\vec{u}}_i\in R^n,i\in [1,n]$, such that $A{\vec{u}}_i={\lambda }_i{\vec{u}}_i$. Equivalently, there exists an orthogonal matrix $U=[{\vec{u}}_1,\cdots ,{\vec{u}}_n]$ (i.e. $U{U}^T=U^{T}U=I_n$) such that: A = U Λ U T = ∑ i = 1 n λ i u ⃗ i u ⃗ i T , Λ = [ λ 1 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ λ n ] , u ⃗ i = [ u i 1 ⋮ u i n ] A=U\Lambda U^T=\sum^n_{i=1}\lambda_i\vec u_i\vec u_i^T, \Lambda=

$$\begin{bmatrix}\lambda_1&\cdots&0\\\vdots&\ddots&\vdots \\0&\cdots&\lambda_n\end{bmatrix}$$ ,\vec u_i= $$\begin{bmatrix}u_{i1}\\\vdots\\u_{in}\end{bmatrix}$$ A=UΛUT=i=1∑n​λi​u i​u iT​,Λ= ​λ1​⋮0​⋯⋱⋯​0⋮λn​​ ​,u i​= ​ui1​⋮uin​​ ​

瑞利商 Rayleigh Quotients

Given a symmertric matrix $A\in S^n$,$$\begin{gathered} {\lambda }_{\min }(A)\le \frac{{\vec{x}}^{T}A\vec{x}}{{\vec{x}}^T\vec{x}}\le {\lambda }_{\max }(A),\forall \vec{x}\ne 0 \\ {\lambda }_{\max }(A)=\max {\vec{x}}^{T}A\vec{x} \\ {\lambda }_{\max }(A)=\max {\vec{x}}^{T}A\vec{x} \end{gathered}$$The maximum and minimum are attained for $\vec{x}={\vec{u}}_1$ and for $\vec{x}={\vec{u}}_n$, where ${\vec{u}}_1$ and ${\vec{u}}_n$ are the largest and smallest eigenvector of $A$. 其中 $\vec{x}$ 要做标准化处理，即 $\Vert \vec{x}{\Vert }_2=1$

Proof
对 ${\vec{x}}^{T}A\vec{x}$ 使用谱定理： x ⃗ T A x ⃗ = x ⃗ T U Λ U T x ⃗ = x ⃗ ′   T Λ x ⃗ ′ = ∑ i = 1 n λ i x i ′ 2 ,   x ⃗ ′ = U T x ⃗ = [ x 1 ′ ⋮ x n ′ ] \vec x^TA\vec x=\vec x^TU\Lambda U^T\vec x=\vec x^{\prime~T}\Lambda\vec x^\prime=\sum^n_{i=1}\lambda_ix_i^{\prime2},~\vec x^\prime=U^T\vec x=

$$\begin{bmatrix}x^\prime_{1}\\\vdots\\x^\prime_{n}\end{bmatrix}$$ x TAx =x TUΛUTx =x ′ TΛx ′=i=1∑n​λi​xi′2​, x ′=UTx = ​x1′​⋮xn′​​ ​ 则有：$${\lambda }_{\min }\sum _{i=1}^n{x}_i^{\prime 2}\le \sum _{i=1}^n{\lambda }_i{x}_i^{\prime 2}\le {\lambda }_{\max }\sum _{i=1}^n{x}_i^{\prime 2}$$
由于 $U$ 是正交矩阵，所以：$$\sum _{i=1}^n{x}_i^{\prime 2}={\vec{x}}^{\prime T}{\vec{x}}^{\prime }=(U^T\vec{x}{)}^T(U^T\vec{x})={\vec{x}}^{T}U{U}^T\vec{x}={\vec{x}}^T\vec{x}=\sum _{i=1}^n{x}_i^2$$
所以：$${\lambda }_{\min }{\vec{x}}^T\vec{x}\le {\vec{x}}^{T}A\vec{x}\le {\lambda }_{\max }{\vec{x}}^T\vec{x}$$

接下来，我们可以开始进行 PCA 的推导，首先我们要解决如下三个问题：

1. 如何判断何为最重要的主成分：可以计算样本在该主成分方向上的方差，方差越大说明在该方向上样本分布越离散，特征也就越明显，方差最大的就是最重要的那个主成分。
2. 如何找次重要的主成分：把最重要的主成分去掉，也就是去掉每个样本在该方向上的分量，然后再选取最重要的主成分。
3. 如何找到第 n 重要的主成分：如果 $n\ne 1$，重复第二个问题的解决方法。

现在我们分别证明上面两个问题的解决方式是正确的：

• 对于 “如何判断何为最重要的主成分” 问题：
  • 首先对样本进行中心化：$$\tilde{X}=[{\tilde{x}}_1,\cdots ,{\tilde{x}}_m],{\tilde{x}}_i={\vec{x}}_i-\bar{x}$$
  • 然后我们假设最重要的主成分为 ${\vec{z}}_1,{\vec{z}}_1\in R^n,\Vert \vec{z}{\Vert }_2=1$，那么样本在 ${\vec{z}}_1$ 上的投影大小为：${\alpha }_i={\tilde{x}}_i^T{\vec{z}}_1$
  • 然后计算方差：$$\frac{1}{m}\sum _{i=1}^m{\alpha }_i^2=\frac{1}{m}{\vec{z}}_1^T\tilde{X}{\tilde{X}}^T{\vec{z}}_1$$
  • 我们使方差最大化，以此来计算 ${\vec{z}}_1$：$$\underset{{\vec{z}}_1\in R^n}{\max }{z}_1^T(\tilde{X}{\tilde{X}}^T){\vec{z}}_1$$
    • 对 $\tilde{X}$ 使用 SVD：$\tilde{X}=U\Sigma V^T={\sum }_{i=1}^m{\sigma }_i{\vec{u}}_i{\vec{v}}_i$
    • 因此 $\tilde{X}{\tilde{X}}^T=\tilde{X}=U\Sigma V^{T}V\Sigma U^T=U{\Sigma }^2{U}^T$（$\Sigma$ 是对角矩阵，因此这里直接用平方形式表示），此时问题转化为：$$\underset{{\vec{z}}_1\in R^n}{\max }{z}_1^{T}U{\Sigma }^2{U}^T{\vec{z}}_1$$
    • 对它使用瑞利商可知，当 ${\vec{z}}_1={\vec{u}}_1$ 时取最大值，也就是最重要的主成分就是第一个特征向量。
• 对于 “如何找次重要的主成分” 问题：
  • 我们记 ${\tilde{x}}_i^{(1)}$ 为去掉最重要的主成分方向上的分量之后的样本，计算方式如下：
    • 先计算 ${\tilde{x}}_i$ 在 $u_1$ 方向上的投影向量，再用 ${\tilde{x}}_i$ 减去这个向量：${\tilde{x}}_i^{(1)}={\tilde{x}}_i-{\vec{u}}_1({\vec{u}}_1^T{\tilde{x}}_i)$ $${\tilde{X}}^{(1)}=[{\tilde{x}}_1^{(1)},\cdots ,{\tilde{x}}_m^{(1)}]=(I_n-{\vec{u}}_1{\vec{u}}_1^T)\tilde{X}$$
    • 把 ${\tilde{X}}^{(1)}$ 展开计算，由于 $U$ 是正交矩阵，所以 ${\vec{u}}_1$ 点乘除它以外的其它特征向量都为 0：$$\begin{array}{rl}{\tilde{X}}^{(1)}& =\sum _{i=1}^m{\sigma }_i{\vec{u}}_i{\vec{v}}_i-({\vec{u}}_1{\vec{u}}_1^T)\sum _{i=1}^m{\sigma }_i{\vec{u}}_i{\vec{v}}_i^T\\ & =\sum _{i=1}^m{\sigma }_i{\vec{u}}_i{\vec{v}}_i-{\sigma }_1{\vec{u}}_1{\vec{v}}_1^T\\ & =\sum _{i=2}^m{\sigma }_i{\vec{u}}_i{\vec{v}}_i\end{array}$$
  • 因此我们计算下式就可以求出 ${\vec{z}}_2$ $$\underset{{\vec{z}}_2\in R^n}{\max }{z}_2^T({\tilde{X}}^{(1)}{\tilde{X}}^{(1)T}){\vec{z}}_2$$
  • 我们发现 ${\vec{z}}_2={\vec{u}}_2$，因此次重要的主成分其实就是第二个特征向量
• 对于 “如何找到第 n 重要的主成分” 问题：
  • 求出 $X^{(n-1)}$，按照上面的步骤计算
  • 显然，${\vec{z}}_n={\vec{u}}_n$，也就是说，特征向量就是主成分，而且越靠前的特征向量越重要

Kernel PCA

---

kernel PCA 通过使用核函数对数据进行非线性变换使得 PCA 具有了非线性分解的能力。

推导过程如下：

• 对于输入数据 $x_i\in R^{n_0}$，使用一个非线性映射函数 $\varphi :R^{n_0}\to R^{n_1}$
• 假设 $\varphi (x_i)$ 已经中心化，即：$$\frac{1}{N}\sum _{i=1}^N\varphi (x_i)=0$$它的协方差矩阵为$$\tilde{H}=\frac{1}{N}\sum _{i=1}^N\varphi (x_i){\varphi }^T(x_i)$$我们要对 $\tilde{H}$ 进行特征值分解，即求解 $\tilde{H}\tilde{z}=\tilde{\lambda }\tilde{z}$
• 把特征值分解方程展开得：$$\begin{array}{rl}\tilde{H}\tilde{z}& =\tilde{\lambda }\tilde{z}\\ \frac{1}{N}\sum _{i=1}^N\varphi (x_i){\varphi }^T(x_i)\tilde{z}& =\tilde{\lambda }\tilde{z}\end{array}$$由于 ${\varphi }^T(x_i)\tilde{z}$ 的结果是一个标量，可以表示为一个数字 ${\alpha }_i$，则有 $$\tilde{z}=\sum _{j=1}^N{\alpha }_j\varphi (x_j)$$，即特征向量可以表示为所有输入的线性组合。我们将此等式带入方程得：$$\begin{array}{rl}\frac{1}{N}\sum _{i=1}^N\varphi (x_i){\varphi }^T(x_i)(\sum _{j=1}^N{\alpha }_j\varphi (x_j))& =\tilde{\lambda }\sum _{j=1}^N{\alpha }_j\varphi (x_j)\\ \frac{1}{N}\sum _{i=1}^N\varphi (x_i)(\sum _{j=1}^N{\alpha }_j{\varphi }^T(x_i)\varphi (x_j))& =\tilde{\lambda }\sum _{j=1}^N{\alpha }_j\varphi (x_j)\end{array}$$
• 我们定义核函数为 $k(x_i,x_j)={\varphi }^T(x_i)\varphi (x_j)$，该函数可以通过 $x_i,x_j$ 直接计算出一个标量，作为 $\varphi (x_i{)}^T\varphi (x_j)$ 的结果，这样可以节省计算 $\varphi$ 函数的时间，而且，由于 $\varphi$ 的计算方式以及维度很难确定，使用核函数也可以大大节省工作量。值得注意的是，并不是所有核函数都有对应的 $\varphi$ 函数。用核函数替换后的方程如下：$$\begin{array}{rl}\frac{1}{N}\sum _{i=1}^N\varphi (x_i)(\sum _{j=1}^N{\alpha }_{j}k(x_i,x_j))& =\tilde{\lambda }\sum _{j=1}^N{\alpha }_j\varphi (x_j)\\ \frac{1}{N}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_{j}k(x_k,x_i)k(x_i,x_j)& =\tilde{\lambda }\sum _{j=1}^N{\alpha }_{j}k(x_k,x_j)\end{array}$$
• 接着可以定义一个 Gram 矩阵 $K\in R^N,K(i,j)=k(x_i,x_j)$，简化上式：$$\begin{array}{rl}{K}^2\vec{\alpha }& =N\tilde{\lambda }K\vec{\alpha }\\ K\vec{\alpha }& =N\tilde{\lambda }\vec{\alpha }\\ K\vec{\alpha }& =\lambda \vec{\alpha }\end{array}$$现在我们可以把 $\vec{\alpha }$ 看作一个特征向量并进行求解，但这一步的前提是 1. $\tilde{z}$ 必须是单位向量 2. $\varphi (x_i)$ 是中心化的
  • $\tilde{z}$ 必须是单位向量：$$\begin{array}{rl}1& ={\tilde{z}}^T\tilde{z}\\ 1& =\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{\varphi }^T(x_i)\varphi (x_j)\\ 1& ={\alpha }^{T}K\alpha \end{array}$$结合 $K\alpha =\lambda \alpha$，可得 ${\alpha }^T\lambda \alpha =1$，因此我们要对 $\alpha$ 进行归一化，使其范数为 $\frac{1}{\lambda }$
  • 对 $\varphi (x_i)$ 进行中心化：$\tilde{\varphi }(x_i)=\varphi (x_i)-\frac{1}{N}{\sum }_{j=1}^N\varphi (x_j)$，则有 $$\begin{array}{rl}\tilde{k}(x_i,x_j)& =(\varphi (x_i)-\frac{1}{N}\sum _{k=1}^N\varphi (x_k){)}^T(\varphi (x_j-\frac{1}{N}\sum _{l=1}^N\varphi (x_l))\\ & =k(x_i,x_j)-\frac{1}{N}\sum _{l=1}^{N}k(x_i,x_l)-\frac{1}{N}\sum _{k=1}^{N}k(x_j,x_k)+\frac{1}{N^2}\sum _{k=1}^N\sum _{l=1}^{N}k(x_k,x_l)\end{array}$$因此有：$\tilde{K}=K-2I_{\frac{1}{N}}K+I\frac{1}{N}KI\frac{1}{N}$
• 通过以上步骤，我们可以求出 $\vec{\alpha }$，但是我们仍然不可以直接求解 $\tilde{z}$，因为涉及到 $\varphi$ 函数，而我们的目的是用核函数替换掉所有的 $\varphi$ 函数来简化计算，所以我们可以通过 $${\varphi }^T(x)\tilde{z}=\sum _{j=1}^N{\alpha }_{j}k(x,x_j)$$ 来计算出每个数据在 $\tilde{z}$ 方向上的投影，以此间接获取我们所需要的信息，该投影即为 $x_j$ 在 $\tilde{z}$ 方向上的坐标，而 $x_j$ 去掉该投影方向的分量后，又可以用来计算下一个主成分，所以即使不知道 $\tilde{z}$ 具体为多少也没关系。

同步更新于：SP-FA 的博客

---