PCA的计算方法

看到网上有一堆“博客”，明显是抄袭的，前后矛盾，自己摸索着写了一个PCA的计算过程。
假设有5个学生的6门功课：语文、数学、地理、化学、英语、历史，成绩如下：

X = np.array(
    [[84,65,61,72,79,81],
     [64,77,77,76,55,70],
     [65,67,63,49,57,67],
     [74,80,69,75,63,74],
     [84,74,70,80,74,82]])

注意，行是样本（表示一个学生），列是特征（表示一门课）。

首先要搞明白什么是协方差。定义：（下面的n是样本数）

均值（假设权重概率都为1）：

\[\mu = \frac{1}{n}\sum_{i=1}^n x_i \]

标准差（除以n-1表示无偏估计）：

\[std = \sqrt{\frac{1}{n-1}\sum_{i=1}^n(x_i-\mu)^2} \]

方差：

\[var = std^2 = \frac{1}{n-1}\sum_{i=1}^n(x_i-\mu)^2 \]

协方差：两个特征之间的方差

\[cov(X,Y) = \frac{1}{n-1}\sum_{i=1}^n (x_i - \mu_x)(y_i-\mu_y) \]

也就是计算所有样本的语文成绩与数学成绩之间的方差，或者化学成绩与英语成绩之间的方差。

用python实现方差协方差计算

def my_mean(data):
    return np.sum(data) / len(data)

def cov(a,b):
    assert(len(a) == len(b))
    mean_a = my_mean(a)
    mean_b = my_mean(b)
    p = (a - mean_a)
    q = (b - mean_b)
    r = np.dot(p,q.T)
    return r/(len(a)-1)

协方差矩阵：多个特征之间的方差的矩阵

\[c = \begin{pmatrix} cov(x,x) & cov(x,y) & cov(x,z) \\ cov(y,x) & cov(y,y) & cov(y,z) \\ cov(z,x) & cov(z,y) & cov(z,z) \\ \end{pmatrix} \]

可以用上面的函数来计算，当然numpy也有现成的函数：np.cov()。使用这个函数时，注意矩阵的形状。以上面的X的数据为例，一共5个样本6个特征，按照协方差矩阵的定义，应该是一个6x6的方阵。
如果用np.cov(X)计算出来的是一个5x5的方阵，不满足条件。仔细阅读cov()函数的文档，X的例子是列为特征，所以应该用：

cc = np.cov(X, rowvar=False)
#或者 
dd = np.cov(X.T)

cc= 
[[ 95.2  -13.9  -23.75  62.15 100.35  63.05]
 [-13.9   41.3   32.75  44.95 -26.95  -5.1 ]
 [-23.75  32.75  40.    42.5  -33.    -8.5 ]
 [ 62.15  44.95  42.5  151.3   53.7   53.85]
 [100.35 -26.95 -33.    53.7  110.8   65.9 ]
 [ 63.05  -5.1   -8.5   53.85  65.9   43.7 ]]

dd= 
[[ 95.2  -13.9  -23.75  62.15 100.35  63.05]
 [-13.9   41.3   32.75  44.95 -26.95  -5.1 ]
 [-23.75  32.75  40.    42.5  -33.    -8.5 ]
 [ 62.15  44.95  42.5  151.3   53.7   53.85]
 [100.35 -26.95 -33.    53.7  110.8   65.9 ]
 [ 63.05  -5.1   -8.5   53.85  65.9   43.7 ]]

当然，也可以用协方差矩阵的定义来求：

c1 = np.array(
    [[cov(x[:,0],x[:,0]),cov(x[:,0],x[:,1]),cov(x[:,0],x[:,2]),cov(x[:,0],x[:,3]),cov(x[:,0],x[:,4]),cov(x[:,0],x[:,5])],
     [cov(x[:,1],x[:,0]),cov(x[:,1],x[:,1]),cov(x[:,1],x[:,2]),cov(x[:,1],x[:,3]),cov(x[:,1],x[:,4]),cov(x[:,1],x[:,5])],
     [cov(x[:,2],x[:,0]),cov(x[:,2],x[:,1]),cov(x[:,2],x[:,2]),cov(x[:,2],x[:,3]),cov(x[:,2],x[:,4]),cov(x[:,2],x[:,5])],
     [cov(x[:,3],x[:,0]),cov(x[:,3],x[:,1]),cov(x[:,3],x[:,2]),cov(x[:,3],x[:,3]),cov(x[:,3],x[:,4]),cov(x[:,3],x[:,5])],
     [cov(x[:,4],x[:,0]),cov(x[:,4],x[:,1]),cov(x[:,4],x[:,2]),cov(x[:,4],x[:,3]),cov(x[:,4],x[:,4]),cov(x[:,4],x[:,5])],
     [cov(x[:,5],x[:,0]),cov(x[:,5],x[:,1]),cov(x[:,5],x[:,2]),cov(x[:,5],x[:,3]),cov(x[:,5],x[:,4]),cov(x[:,5],x[:,5])]])
print("c1=",c1)

c1= 
[[ 95.2  -13.9  -23.75  62.15 100.35  63.05]
 [-13.9   41.3   32.75  44.95 -26.95  -5.1 ]
 [-23.75  32.75  40.    42.5  -33.    -8.5 ]
 [ 62.15  44.95  42.5  151.3   53.7   53.85]
 [100.35 -26.95 -33.    53.7  110.8   65.9 ]
 [ 63.05  -5.1   -8.5   53.85  65.9   43.7 ]]

可以看到cc,dd,c1三者是一样的。

得到协方差矩阵后，使用函数求其特征值和特征向量：

eig_value, eig_vector = np.linalg.eig(cc)

eig_value: 
[3.06293191e+02 1.63510310e+02 9.89302953e+00 2.60347035e+00 3.77194607e-15 1.37974973e-14]

eig_vector: 
[[-0.53264253  0.20279107 -0.34433806  0.39437042 -0.61869481 -0.55543331]
 [ 0.00876193 -0.46059524 -0.81597078  0.02185232  0.25842516  0.34848844]
 [ 0.04593605 -0.47328385  0.37877077  0.70892582 -0.03144886  0.21014772]
 [-0.51955599 -0.64238594  0.24891406 -0.45230979 -0.15412561 -0.22434743]
 [-0.55131936  0.32775478  0.09651389 -0.13044526  0.29446728  0.67491022]
 [-0.37445103  0.05145202  0.0297077   0.34614812  0.66255449  0.14160509]]

因为要降维，原来是6维，我们降为2维，所以取eig_value中两个最大的值：

eig_sort_index = np.argsort(eig_value)
# 得到结果 [4,5,3,2,1,0] 是个从小到大的排序数组下标，表示第0个值最大，第1个值其次，第4个值最小
eig_big_index = eig_sort_index[:-3:-1] # 两个最大值，[0,1]
eig_big_vector = eig_vector[:,eig_big_index]

得到eig_big_vector的结果：

array([[-0.53264253,  0.20279107],
       [ 0.00876193, -0.46059524],
       [ 0.04593605, -0.47328385],
       [-0.51955599, -0.64238594],
       [-0.55131936,  0.32775478],
       [-0.37445103,  0.05145202]])

这实际上是两个特征向量的组合，每一列是一组特征向量。

然后得到X的mean均值矩阵：

mean = np.mean(x, axis=0)
print(mean)
x_mean = x - mean
print("x_mean=",x_mean)

mean=[74.2 72.6 68.  70.4 65.6 74.8]
x_mean= 
[[  9.8  -7.6  -7.    1.6  13.4   6.2]
 [-10.2   4.4   9.    5.6 -10.6  -4.8]
 [ -9.2  -5.6  -5.  -21.4  -8.6  -7.8]
 [ -0.2   7.4   1.    4.6  -2.6  -0.8]
 [  9.8   1.4   2.    9.6   8.4   7.2]]

最后用特征向量与x_mean相乘，得到PCA的降维结果：

r = np.dot(x_mean, eig_big_vector)
print("r=", r)

r= 
[[-16.14860528  12.48396235]
 [ 10.61676743 -15.67317428]
 [ 23.40212697  13.607117  ]
 [ -0.43966353  -7.77054621]
 [-17.43062559  -2.64735885]]