PCA

推导

矩阵乘法的几何意义

假定有一个行向量组成的矩阵

$$P=\begin{pmatrix}p_1\\p_2\\\vdots\\p_n\end{pmatrix}$$

和列向量组成的矩阵

$$A=\begin{pmatrix}a_1,a_2,\cdots,a_m\end{pmatrix}$$

, 则有

$$PA=\begin{bmatrix} {p_1a_1}&{p_1a_2}&{\cdots}&{p_1a_m}\\ {p_2a_1}&{p_2a_2}&{\cdots}&{p_2a_m}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {p_na_1}&{p_na_2}&{\cdots}&{p_na_m}\\ \end{bmatrix}$$

考虑向量点积的几何意义 , 当 时， 即为 在 方向的上投影。当 由一组正交基组成时， 可以将每一个列向量 投影到由原来的特征线性组合而成的新的坐标系中，其中每个行向量 都可以被视为一个组合出的新特征，而每个列向量 则是一个样本。

特征分解

上一节介绍了矩阵乘法的几何意义。PCA 的核心思想即选取原来特征的线性组合组成新特征，然后将原来的样本映射到新的空间中，如果选取得当，这些新特征承载了更多的数据，那么就可以用更少的特征来描述样本，从而进行数据的降维。被选取的特征即为主成分 principal component。

现在的主要问题在于，如何选取主成分。假定有列向量组成的数据矩阵 和行向量组成的投影矩阵 ，我们希望被投影的结果尽可能保留原来数据中的信息，常见的方法是使投影后的数据在主成分的方向上方差最大。此外，当我们选取多个主成分时，我们希望不同的主成分是线性无关的，从而可以使得不同的主成分之间没有冗余的信息。

因此我们对投影所得的新数据集 提出两个目标：

• 在选取的主成分方向上，数据的方差尽可能大
• 不同的主成分之间相关系数为0

幸运的是，线性代数中已经有了强有力的工具来分析这两个目标。考虑协方差矩阵

$$\frac{1}{m}DD^T = \begin{bmatrix} {\frac{1}{m}\sum_{i=1}^md_{i1}^2}&{\frac{1}{m}\sum_{i=1}^md_{i1}d_{i2}}&{\cdots}&{\frac{1}{m}\sum_{i=1}^md_{i1}d_{in}}\\ {\frac{1}{m}\sum_{i=1}^md_{i2}d_{i1}}&{\frac{1}{m}\sum_{i=1}^md_{i2}^2}&{\cdots}&{\frac{1}{m}\sum_{i=1}^md_{i2}d_{in}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {\frac{1}{m}\sum_{i=1}^md_{in}d_{i1}}&{\frac{1}{m}\sum_{i=1}^md_{in}d_{i2}}&{\cdots}&{\frac{1}{m}\sum_{i=1}^md_{in}^2}\\ \end{bmatrix}$$

，当 时，即在任意特征上均值都为0时，有

$$\frac{1}{m}DD^T=\begin{bmatrix} {Var(D_1)}&{Cov(D_1,D_2)}&{\cdots}&{Cov(D_1,D_n)}\\ {Cov(D_2,D_1)}&{Var(D_2)}&{\cdots}&{Cov(D_2,D_n)}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {Cov(D_n,D_1)}&{Cov(D_n,D_2)}&{\cdots}&{Var(D_n)}\\ \end{bmatrix}$$

，其中 为样本特征数， 为样本容量， 表示样本 的第 个特征， 表示数据集中第 个特征的随机变量。

要满足上述的两个条件，即希望我们能得到一个对角矩阵，只有主对角线上的元素有值，其他位置均为0。换言之，我们希望 。

考虑 ，则有

$$\Lambda=\frac{1}{m}DD^T=P\cdot \frac{1}{m}AA^T\cdot P^T$$

，即

$$\frac{1}{m}AA^T=P^{-1}\Lambda (P^T)^{-1}$$

如果我们对 进行特征值分解，由于该矩阵为对称矩阵(显然具有对称性)，分解可得一组标准正交的特征向量，即分解可得

$$\frac{1}{m}AA^T=E\Lambda E^T$$

，其中 为正交的，有 。令 ，可得

$$\frac{1}{m}AA^T=P^T\Lambda P$$

因此对 进行特征分解，转置即可得到所需的投影矩阵 。由此可得 。

特征选取

虽然得到了投影矩阵，但是从这种特征中选取有价值的部分仍然是一个问题。一个自然的结论是，更大的特征值对应的特征向量将原数据投影到的特征上后方差更大，因此应该选择特征值更大的特征向量。以下给出这一结论的证明。

考虑一个特征向量 和其对应的特征值 ，假设数据矩阵 的协方差矩阵为 ，则用 对数据 投影所得数据之方差为 ，我们要做的就是最大化这个方差，即

$$\arg\max_u u^T\Sigma u$$

分解协方差矩阵可以得到 ，则

$$u^T\Sigma u=\sum_{i=1}^{n}\lambda_i(u^Tp_i)^2$$

，其中 为第 个特征向量。因为 是正交的，因此

$$\sum_{i=1}^{n}(u^Tp_i)^2=1$$

，因此可以将问题转化为如下优化问题

$$\max\sum_{i=1}^{n}\lambda_iz_i^2\quad s.t.\quad \sum_{i=1}^{n}z_i^2=1$$

不妨假设最大特征值为 ，则我们希望 ，且 。由于 为正交标准矩阵，当 时取得最大值，且最大值为 。

综上所示，选取最大特征值对应特征向量可以得到最大的方差，且方差为该特征向量。

实现

算法实现

与推导不同，实现时我们往往使用行向量来处理样本；而在 sklearn 等常见的科学计算包内，特征值分解所得的特征向量矩阵往往以列向量的形式给出，因此在实际书写时，顺序会略有偏差。

考虑数据矩阵 ，首先将数据每一特征平移至均值为0的位置

n_sample, n_feature = X.shape 
self.mean_ = np.array([np.mean(X[:, i]) for i in range(n_feature)]) 
X = X - self.mean_

然后计算协方差矩阵

cov = np.dot(X.T, X) / (n_sample)

使用 sklearn 计算矩阵特征值，并按照特征值大小进行排序，重新组装成矩阵并求转置，得到按照特征值大小排序的投影矩阵

D, Q = np.linalg.eig(cov) # cov = QDQ^T
eig_pairs = [(np.abs(D[i]), Q[:,i]) for i in range(len(D))]
eig_pairs.sort(key=lambda x: x[0], reverse=True)
self.feature_ = np.array([eig_pairs[i][1] for i in range(n_feature)]).T

最后选取前 个特征向量作为新的特征对数据集进行变换

def transform_k(self, X, k):
    X = X - self.mean_
    return np.dot(X, self.feature_[:,:k])

需要注意的是被变换的数据集也要对均值进行处理，并且因为如前文所述，两个矩阵均为推导过程中的转置，所以计算式也进行相应的变化。

SVD 解法

考虑矩阵 ，使用特征值分解解决 PCA 存在一个问题，即协方差矩阵 可能维度较高，特征值分解成本巨大。因此考虑更简单高效的替代方案。

不难发现，SVD 奇异值分解和特征值分解求 PCA 有着相似的问题形式。任意形状的矩阵均可得到如下分解 ，其中 为正交标准基， 主对角线为奇异值。

要求解 ，考虑

$$AA^T=U\Sigma V^T(U\Sigma V^T)^T=U\Sigma V^TV\Sigma U^T=U\Sigma^2 U^T$$

，对 进行特征值分解，即可得到 和 。不难发现，该式形式与特征值分解协方差的表达式一致，因此有 ，即 ；且 即对应的特征向量。（同理，要求解 只需要对 做特征值分解即可。）因此，只要能求解 SVD ，得到左奇异向量 ，即可得到 PCA 所需的投影矩阵，并按照奇异值大小进行排序进行特征提取。

事实上，SVD 同样也是一种降维算法。对于数据矩阵 ，不妨假设其为 条 维的列向量，表示 条特征数为 的数据。对 做奇异值分解得到 ，其中左奇异向量 ，奇异值 ，右奇异向量 。我们可以使用部分左奇异向量 来提取特征，压缩维度；同理，也可以使用部分右奇异向量 来提取数据，去除冗余样本。

因此，从某种角度来说，SVD 的功能是 PCA 的超集，使用 SVD 同样也可以求解 PCA 。但问题是，使用特征值分解来求解 SVD 的运算量甚至更大，因此我们考虑采用更高效准确的迭代解法，避免了分解特征值。

关于 SVD 的迭代解法，可以参考 SVD 的牛顿迭代法解法，本文中不多赘述。使用迭代解法求解 SVD，可以高效准确的得到 PCA 的解，因此科学计算包中的 PCA 往往使用 SVD 求解。

简单实现

import numpy as np

class PCA:
    def fit(self, X):
        '''
        fit PCA to X
        '''
        n_sample, n_feature = X.shape 
        self.mean_ = np.array([np.mean(X[:, i]) for i in range(n_feature)]) 
        X = X - self.mean_

        # covariance matrix
        cov = np.dot(X.T, X) / (n_sample)

        # eigenvalue and eigenvector
        D, Q = np.linalg.eig(cov) # cov = QDQ^T
        eig_pairs = [(np.abs(D[i]), Q[:,i]) for i in range(len(D))]
        eig_pairs.sort(key=lambda x: x[0], reverse=True)

        self.feature_ = np.array([eig_pairs[i][1] for i in range(n_feature)]).T

    def transform(self, X):
        X = X - self.mean_
        return np.dot(X, self.feature_)

    def transform_k(self, X, k):
        '''
        transform X to k-dimensional space
        '''
        X = X - self.mean_
        return np.dot(X, self.feature_[:,:k])

class PCA_SVD(PCA):
    def fit(self, X):
        '''
        fit PCA to X
        '''
        n_sample, n_feature = X.shape 
        self.mean_ = np.array([np.mean(X[:, i]) for i in range(n_feature)]) 
        X = X - self.mean_

        U,S,V_T = np.linalg.svd(X)
        D = (S**2)/(n_sample)

        eig_pairs = [(np.abs(D[i]), V_T[:,i]) for i in range(len(D))]
        eig_pairs.sort(key=lambda x: x[0], reverse=True)

        self.feature_ = np.array([eig_pairs[i][1] for i in range(n_feature)])


if __name__ == '__main__':
    # X = np.array([[-1,1],[-2,-1],[-3,-2],[1,1],[2,1],[3,2]])
    # X = np.array([[-1,-2],[-1,0],[0,0],[2,1],[0,1]])
    X = np.array([[-1,1,0],[-4,3,0],[1,0,2]])

    pca = PCA()
    pca.fit(X)
    print(pca.transform_k(X,1))

    pca2 = PCA_SVD()
    pca2.fit(X)
    print(pca2.transform_k(X,1))

    from sklearn.decomposition import PCA as pcca
    pca3 = pcca(n_components=1)
    pca3.fit(X)
    print(pca3.transform(X))

    print(pca.feature_)
    print(pca2.feature_)