PCA算法,英文全称principal component analysis,翻译为主成分分析。
在机器学习中,特征向量的维度有时候成百上千,甚至更多。面对这么多的影响因素,人们自然而然的就想到了抓住主要矛盾,只找出对结果影响最大的因素进行分析。pca给人们提供了这样的思路。

基本思想

举个例子,给各门功课的成绩分配不同的权重,依次来区分学生的总成绩。
区分学生的成绩,也就是意味着使总成绩的方差更大。
x1,x2,......xp 表示 p 门功课成绩,用 s 表示总成绩,权重分别为 c1,c2,c3,c4......cp

s=c1x1+c2x2+c3x3+.....+cpxp
一共有 n 名学生。用 X1,X2,X3,.....,Xp 表示样本值为 x1,x2,......xp 的随机变量。
这时候,问题就更加清晰了:

find C=[c1,c2,c3,c4......cp]

to maximize Var(c1X1+c2X2+c3X3+.....+cpXp)

但是,这样的解有无限个,一般加上限制条件:
c21+c22+c23+c24......c2p=1

每得到一个C就得到一个主成分值Z, 每个C向量间正交。因为,作为主成分,主成分间应该没有关系。

另一个角度

假设我们要得到 k 个主成分,矩阵 Zn,k 作为降维后的数据,矩阵 Xn,d 表示原始的 n d 维数据, 矩阵 Pd,k k 列主成分。那么有:
Z=XP
此时, P k 个列向量正交。
Z作为要处理的数据,Z的各个feature之间应该保持正交。
所以 ZTZ=PTXTXP 是对角矩阵。
所以问题简化为:求矩阵P使得 XTX 对角化。由于其为实对称矩阵那么 P 一定存在。

从统计学来看,矩阵 ZTZ
Var(Z1)=1nni=1(zi1μ1)2
如果所有数据平均值为0.
Var(Z1)=1nni=1(zi1)2
Cov(Z1,Z2)=1nni=1(zi1zi2)

所以,

1nZTZ=σ11σ21.σk1σ12σ22.σk2σ13σ23.σk3.......σk4σ1kσ2k.σk5

实现

主成分回归分析的实现。

序号 x1 x2 x3 x4 y
1 7 26 6 60 78.5
2 1 29 15 52 74.3
3 11 56 8 20 104.3
4 11 31 8 47 87.6
5 7 52 6 33 95.9
6 11 55 9 22 109.2
7 3 71 17 6 102.7
8 1 31 22 44 72.5
9 2 54 18 22 93.1
10 21 47 4 26 115.9
11 1 40 23 34 83.8
12 11 66 9 12 113.3
13 10 68 8 12 109.4
#!usr/bon/env python
#-*- coding:utf-8 -*-

import numpy as np


data = np.array([[7, 26, 6, 60],\
                 [1, 29, 15, 52],\
                 [11, 56, 8, 20],\
                 [11, 31, 8, 47],\
                 [7, 52, 6, 33],\
                 [11, 55, 9, 22],\
                 [3, 71, 17, 6],\
                 [1, 31, 22, 44],\
                 [2, 54, 18, 22],\
                 [21, 47, 4, 26],\
                 [1, 40, 23, 34],\
                 [11, 66, 9, 12],\
                 [10, 68, 8, 12]])
y = np.array([78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7,\
              72.5, 93.1, 115.9, 83.8, 113.3, 109.4])
#regularation
c = (y - np.mean(y))/np.std(y)
sd = np.std(data, axis=0) 
avg = np.mean(data, axis=0)
data = data - avg
data /= sd
x = np.dot(data.T, data)
x /= 13
# x = 
#array([[ 1. , 0.22857947, -0.82413376, -0.24544511],
# [ 0.22857947, 1. , -0.13924238, -0.972955 ],
# [-0.82413376, -0.13924238, 1. , 0.029537 ],
# [-0.24544511, -0.972955 , 0.029537 , 1. ]])
#
eigval, eigvec = np.linalg.eig(x)


z = np.dot(data, eigvec[:,:3])
a = np.dot(z.T,z)
b = np.dot(z.T,c)
np.linalg.solve(a,b)
# array([ 0.65695805, -0.00830863, 0.30277024])

eigval 为 array([ 2.23570403e+00, 1.57606607e+00, 1.86606149e-01, 1.62374573e-03])
可以看出最后一个特征值很小,选取从大到小的三个特征值分析即可。
这三个特征值形成的累积贡献率为99.96%
特征向量为:
array([[ 0.47595517, 0.50897938, 0.67550019, 0.24105218],
[ 0.56387024, -0.41393149, -0.31442044, 0.64175607],
[-0.39406653, -0.60496908, 0.63769109, 0.26846611],
[-0.54793119, 0.45123511, -0.19542096, 0.67673402]])
所以需要保留的特征向量为:
array([[ 0.47595517, 0.50897938, 0.67550019],
[ 0.56387024, -0.41393149, -0.31442044],
[-0.39406653, -0.60496908, 0.63769109],
[-0.54793119, 0.45123511, -0.19542096,])

z1=0.476x1+0.564x20.394x30.548x4
z2=0.509x10.414x20.605x3+0.4512x4
z3=0.676x10.314x2+0.638x30.195x4
拟合得到:
y=0.657z10.0083z2+0.303z3
x1x2,x3,x4 表示为:

y=0.513x1+0.279x20.061x30.423x4