揭秘PCA与协方差的关系

这两天学习了PCA，学完后一头雾水，怎么求出协方差矩阵的特征值和特征向量就可以降维了，它的原理到底是什么？经过一番查找，终于明显这其中的弯弯绕绕了。

1 为什么需要PCA

PCA（主成分分析），他是特征提取的一种方法，主要目的是对数据进行降维，减小计算量。

特征可分为三类：

• 相关特征：对于学习任务（例如分类问题）有帮助，可以提升学习算法的效果；例如猫狗分类中，嘴巴这一维特征。
• 无关特征：对于我们的算法没有任何帮助，不会给算法的效果带来任何提升；例如猫狗分类中，颜色这一特征，我们并不能根据颜色来判断是猫还是狗，因此颜色是无关特征。
• 冗余特征：不会对我们的算法带来新的信息，或者这种特征的信息可以由其他的特征推断出。例如猫狗分类中，其中的特征包括体重、头部宽度、身长、腿的数量和尾巴长度。如果我们已经有了体重和身长这两个特征，那么尾巴长度就可能是一个冗余特征，因为它与体重和身长之间可能存在高度相关性。

当特征维度过多时，学习任务的计算量和复杂性会急剧增加。为了应对这个问题，可以使用主成分分析（Principal Component Analysis，PCA）对特征进行降维。注意，PCA是特征提取的一种方法，它得到的结果是原特征的线性组合，是新的特征，在是原特征中不存在的。

2 PCA与协方差

为什么是协方差矩阵？

首先，方差的意义是什么？
对于一组数据，它的方差越大，说明坐标点越分散，那么改属性就能够比较好的反映源数据。
方差的计算公式是：

协方差就是一种用来度量两个随机变量关系的统计量。
同一元素的协方差就表示该元素的方差，不同元素之间的协方差就表示它们的相关性。
协方差的计算公式是：

由公式可以看出，方差就是协方差的一种特殊情况，同一个元素的协方差就是它的方差。

协方差矩阵

比如，三维(x,y,z)的协方差矩阵：

PCA算法优化的目标是：① 降维后同一维度的方差最大； ② 不同维度之间的相关性为0

仔细一看，按照这个要求的话，就是说降维后的特征，它的协方差矩阵需要是一个除了对角线元素外，其他地方全为0，且对角线元素值从上到下按从大到小进行排列，这…这不就是矩阵的对角化嘛！接下来我们看一下原矩阵与变换后的协方差矩阵的关系，此处引用了以为大佬的博客。

http://blog.codinglabs.org/articles/pca-tutorial.html

现在问题全部都聚焦在矩阵对角化中，这在线性代数中是非常常见的，这里原理不再赘述，总之，P就是协方差矩阵的特征向量单位化后按行排列出的矩阵。终于真相大白了！

3 代码实现

3.1 不调库pca代码实现

import numpy as np
 
class CPCA(object):
    '''用PCA求样本矩阵X的K阶降维矩阵Z
    Note:请保证输入的样本矩阵X shape=(m, n)，m行样例，n个特征
    '''
    def __init__(self, X, K):
        '''
        :param X,训练样本矩阵X
        :param K,X的降维矩阵的阶数，即X要特征降维成k阶
        '''
        self.X = X       #样本矩阵X
        self.K = K       #K阶降维矩阵的K值
        self.centrX = [] #矩阵X的中心化
        self.C = []      #样本集的协方差矩阵C
        self.U = []      #样本矩阵X的降维转换矩阵
        self.Z = []      #样本矩阵X的降维矩阵Z
        
        self.centrX = self._centralized()
        self.C = self._cov()
        self.U = self._U()
        self.Z = self._Z() #Z=XU求得
        
    def _centralized(self):
        '''矩阵X的中心化'''
        print('样本矩阵X:\n', self.X)
        centrX = []
        mean = np.array([np.mean(attr) for attr in self.X.T]) #样本集的特征均值
        print('样本集的特征均值:\n', mean)
        centrX = self.X - mean ##样本集的中心化
        print('样本矩阵X的中心化centrX:\n', centrX)
        return centrX
        
    def _cov(self):
        '''求样本矩阵X的协方差矩阵C'''
        #样本集的样例总数
        ns = np.shape(self.centrX)[0]
        #样本矩阵的协方差矩阵C
        C = np.dot(self.centrX.T, self.centrX)/(ns - 1)
        print('样本矩阵X的协方差矩阵C:\n', C)
        return C
        
    def _U(self):
        '''求X的降维转换矩阵U, shape=(n,k), n是X的特征维度总数，k是降维矩阵的特征维度'''
        #先求X的协方差矩阵C的特征值和特征向量
        a, b = np.linalg.eig(self.C) #特征值赋值给a，对应特征向量赋值给b。函数doc：https://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.linalg.eig.html
        print('样本集的协方差矩阵C的特征值:\n', a)
        print('样本集的协方差矩阵C的特征向量:\n', b)
        #给出特征值降序的topK的索引序列
        ind = np.argsort(-1*a)
        #构建K阶降维的降维转换矩阵U
        UT = [b[:, ind[i]] for i in range(self.K)] #  UT 是一个包含 K 个列向量的列表。
        U = np.transpose(UT) # U 是一个 K x M 的降维转换矩阵，其中每一行代表一个降维后的特征向量。
        print('%d阶降维转换矩阵U:\n'%self.K, U)
        return U
        
    def _Z(self):
        '''按照Z=XU求降维矩阵Z, shape=(m,k), n是样本总数，k是降维矩阵中特征维度总数'''

        # 降维矩阵 U 的形状为 (M, K)，原始数据矩阵 X 的形状为 (N, M)，降维后的矩阵 Z 的形状为 (N, K)。

        Z = np.dot(self.X, self.U)
        print('X shape:', np.shape(self.X))
        print('U shape:', np.shape(self.U))
        print('Z shape:', np.shape(Z))
        print('样本矩阵X的降维矩阵Z:\n', Z)
        return Z
        
if __name__=='__main__':
    '10样本3特征的样本集, 行为样例，列为特征维度'
    X = np.array([[10, 15, 29],
                  [15, 46, 13],
                  [23, 25, 30],
                  [11, 20,  35],
                  [42, 45, 11],
                  [9,  38, 5],
                  [11, 21, 14],
                  [8,  5,  15],
                  [11, 12, 21],
                  [21, 20, 25]])
    K = np.shape(X)[1] - 1
    print('样本集(10行3列，10个样例，每个样例3个特征):\n', X)
    pca = CPCA(X,K)

3.2 简化版本的pca代码

import numpy as np
class PCA():
    def __init__(self,n_components):
        self.n_components = n_components
    
    def fit_transform(self,X):
        self.n_features_ = X.shape[1]
        # 求协方差矩阵
        X = X - X.mean(axis=0)
        self.covariance = np.dot(X.T,X)/X.shape[0]
        # 求协方差矩阵的特征值和特征向量
        eig_vals, eig_vectors = np.linalg.eig(self.covariance)
        # 获得降序排列特征值的序号
        idx = np.argsort(-eig_vals)
        # 降维矩阵
        self.components_ = eig_vectors[:,idx[:self.n_components]]
        # 对X进行降维
        return np.dot(X,self.components_)
 
# 调用
pca = PCA(n_components=2)
X = np.array([[-1, 2, 66, -1], [-2,6,58,-1], [-3,8,45,-2], [1,9,36,1], [2,10,62,1], [3,5,83,2]])  #导入数据，维度为4
newX=pca.fit_transform(X)
print(newX)                  #输出降维后的数据

3.3 sklearn库中的PCA

import numpy as np
from sklearn.decomposition import PCA
X = np.array([[-1,2,66,-1], [-2,6,58,-1], [-3,8,45,-2], [1,9,36,1], [2,10,62,1], [3,5,83,2]])  #导入数据，维度为4
pca = PCA(n_components=2)   #降到2维
pca.fit(X)                  #训练
newX=pca.fit_transform(X)   #降维后的数据
# PCA(copy=True, n_components=2, whiten=False)
print(pca.explained_variance_ratio_)  #输出贡献率
print(newX)                  #输出降维后的数据