主成分分析
(原创文章，转载请注明出处！)
一、主成分分析的作用
主成分分析（PCA）是通过一些方法将高维的训练数据映射到低维，起到一个训练数据降维的作用。
这样一方面能消除训练数据中冗余特征，另一方面能提升训练算法的运行效率，还能减少对存储资源的占用；如果是降维到2D或者3D，还能用2D、3D的展示技术将数据显示出来，方便分析。
主成分分析这么多的好处，那是不是时时处处都要使用它？
答案是否定的。通过主成分分析来把数据降维时，原始数据中所包含的信息，是会有一定的丢失的（有损压缩）。
所以，当训练算法的时候，存储资源足够，算法的运行效率也符合要求，那么就不需要使用PCA来处理数据。
另外一点，PCA可以用来减少训练数据的维度，那么就相当于减少了算法模型中的特征数，相应的要学习的参数也就减少了，这样可以减少overfitting的可能性（不过，如果是想避免overfitting还是应该使用正则化技术，而不是PCA）。
在机器学习算法训练前，用训练样本集中的数据（不包含测试数据集）计算出，用于将数据从高维映射到低维的映射矩阵；在进行算法的测试的时候，先用已知的映射矩阵将测试数据集中的数据映射到低维，然后再进行计算。
二、主成分分析的过程
PCA的目标是寻找一个超平面，使训练数据点到这个超平面的距离之和（投影误差）最小。把训练数据点投影到以这个超平面空间。
这个超平面空间，由一组基向量所确定。主成分分析的目标就是寻找一组基向量，然后使用基向量来将训练数据点投影到低维的超平面空间。
要到达降维的作用，所确定的一组基向量的个数小于训练数据点的维数。
1. 数据预处理
- 使用均值正则化（使每个特征的均值为0）来处理训练数据点，比如：有5个3维的数据点：x1=(0.5, 11,100)T
x2=(0.4, 8,130)T
x3=(0.8, 9,120)T
x4=(0.2, 3,104)T
x5=(0.3, 5,114)T
第一个特征的均值： u(1)= (1/m){∑i=1m(xi(1))}, 其中m=5
= (1/5) * (0.5 + 0.4 + 0.8 + 0.2 + 0.3) = 0.44
那么，新的x1(1) = 0.5 - 0.44 = 0.06， x2(1) = 0.4 - 0.44 = -0.04
- 进行数量级的scaling
像上面例子中，三个feature的数量级不同，就需要使他们在数量级上一致，这样才能有可比较性。
可以使用如下方法：
d(1)= (1/(5-1))*{(0.5-0.44)2 + (0.4-0.44)2 + (0.8-0.44)2 + (0.2-0.44)2 + (0.3-0.44)2} = 0.053
σ(1)=sqrt(d(1)) = 0.23
那么，新的x1(1) = (0.5 - 0.44) / σ(1) = 0.26
u(2)= (1/5) * (11 + 8 +9 + 3 + 5) = 7.2 ,
d(2)= (1/(5-1))*{(11-7.2)2 + (8-7.2)2 + (9-7.2)2 + (3-7.2)2 + (5-7.2)2} = 10.2
σ(2)=sqrt(d(2)) = 3.19
那么，新的x1(2) = (11 - 7.2) / σ(2) = 3.8 / 3.19 = 1.19
u(3)= (1/5) * (100 + 130 + 120 + 104 + 114) = 113.6 ,
d(3)= (1/(5-1))*{(100-113.6)2 + (130-113.6)2 + (120-113.6)2 + (104-113.6)2 + (114-113.6)2} = 146.8
σ(3)=sqrt(d(3)) = 12.12
那么，新的x1(3) = (100 - 113.6) / σ(3) = -13.6 / 12.12= -1.22
通过以上的预处理，是数据样本点各个特征在相同的scale上。
2. PCA计算
计算样本点的协方差矩阵，然后对协方差矩阵进行Singlar Value Dicomposition, 得到超平面空间的基向量。
(1) 协方差矩阵 = (1/m) ∑i=1m(xi *xiT), m个样本点，每个样本点有n个特征，(xi *xiT)是一个n-by-n的方阵
m个方阵相加，然后除以m，最后的结果是n-by-n方阵。(xi要先经过预处理)
对上述的例子，经过预处理(均值正则化、scaling)后：
x1=(0.26, 1.19, -1.12)，x2=(-0.17, 0.25, 1.35)，x3=(1.57, 0.56, 0.53)，x4=(-1.04, -1.31, -0.79)，x5=(-0.61, -0.69, 0.03)
m=5, n=3, 计算协方差矩阵的第一行第二个元素c12：
x1 *x1T = 
= CX11， x2 *x2T = 
= CX22
x3 *x3T = 
= CX33, x4 *x4T = 
= CX44 ,
x5 *x5T = 
= CX55
c12等于这5个矩阵的第一行第二个元素之和，然后除以5 ：
c12 = (1/(5-1)) * [0.3094 + (-0.0425) + 0.8792 + 1.3624 + 0.4209] = 0.73 。
记第一个特征t1，第二个特征t2， {E[x1(1)-E(t1)] * E[x1(2)-E(t2)]}/{sqrt(D(t1)) * sqrt(D(t2))}, c12的计算式中的五个值都是有这个公式计算而得，而这个公式计算的是系数，那么c12就是第一、二个特征的相关系数的平均值。（总体的第一、二个特征的相关系数用样本的平均值代替）
如果不进行除以sd（标准方差）的scaling，则计算的是协方差。
对本例，最终的相关系数矩阵为：
(2) 对该协方差（相关系数）矩阵进行主成分分解
a. 特征值分解
经过标准化之后得到相关系数矩阵是对称正定矩阵，可以通过特征值分解得到主成分。
z=(z1,z2,z3,...,zn)T=(v1,v2,v3,...,vn)Tx=QTx, 其中：Q=(v1,v2,v3,...,vn), v1,v2,v3,...,vn是相关系数矩阵的特征向量，相互正交，对应的特征值为λ1≥λ2≥λ3≥...≥λn
zi=viTx，为第i个主成分。z1是第一主成分，z2是第二主成分，......
把n个原始变量x1,x2,x3,...,xn的总方差分解成n个不相关变量z1,z2,z3,...,zn的方差之和：∑i=1nVar(zi)=∑i=1nVar(xi)，z1,z2,z3,...,zn一起能完全解释原始变量x1,x2,x3,...,xn 。每个变量方差的大小描述了变量所携带的信息量，方差大说明变量携带的信息多。
λi/{∑i=1nλi}为第i个主成分zi的贡献率，z1的贡献率最大，表示它对原始变量x1,x2,x3,...,xn的解释能力最强。
主成分分析，希望用前k（p<n）个主成分来解释原始变量，且前k个主成分总贡献率要大于一定的阈值，比如：95%。这样数据的维度降低了，但降维所带来的数据信息丢失又是在可接收的范围之内。
可以看出使用特征值分解的主要工作就是求相关系数矩阵的特征值以及相应的特征向量，然后选择前k个（由大到小）特征值对应的特征向量，来进行将原始数据映射到新的子空间。
b. 奇异值分解（SVD）
奇异值分解将一个比较复杂的矩阵用更小更简单的几个子矩阵的相乘来表示，这些小矩阵描述了矩阵的重要特性。假设原始矩阵A是m-by-n，奇异值分解如下：
A=U∑VT，其中：U是m-by-m方阵，∑是m-by-n的对角阵，V是n-by-n方阵。U、V各自的向量正交。
ATAvi=λivi, 即V中的列向量是矩阵ATA的特征向量，λi是每个向量对应的特征值；
σi=sqrt(λi), 即对角阵∑对角线上的元素为相应的特征值λi开平方，称为奇异值；
ui=(1/σi)Avi, 为矩阵U的每列的向量。
∑中对角线上的元素：σ1≥σ2≥σ3≥...≥σn，即：λ1≥λ2≥λ3≥...≥λn，所以U和V中的列向量依相应特征值的顺序排列。
主成分计算：z=VTx，类似特征值分解，大的奇异值对应的变量所能反应的数据的特征就多，所以取前K大的奇异值对应的矩阵V中的列向量构成的n-by-k矩阵。
每个样本点x是 n-by-1的向量，VT是k-by-n，VTx 则得到k-by-1的向量，达到对数据降维的作用，当然降维会带来一定的信息损失。
在满足 {∑i=1kσi} / {∑i=1nσi} ≥ 一定阈值（比如：95%）条件下，选择尽量小的K。
(如果是UTA，选择U的前k个列向量，即：UT是k-by-m，而A是m-by-n，则得到k-by-n的矩阵，就是把样本点个数进行了压缩)
(3)从z得到x
当把x从n维映射到k维的z之后，通过z只能重构x的近似向量，而不能完全的回到x，即：
xapprox≈Vn-by-k*zk-by-1, 得到约等于x的n维向量。
计算主成分的贡献率，也可以通过公式： {∑i=1k||xioriginal-xiapprox||2} / {∑i=1n||xioriginal||2}
三、用R编程实现PCA
1、eigen
eigen函数用来计算一个矩阵的特征值和特征向量。通过对相关系数矩阵计算特征值和特征向量，来获得主成分，程序代码如下：
1 x1 <- c(0.5, 11,100) 2 x2 <- c(0.4, 8,130) 3 x3 <- c(0.8, 9,120) 4 x4 <- c(0.2, 3,104) 5 x5 <- c(0.3, 5,114) 6 original_x <- matrix(c(x1,x2,x3,x4,x5), nrow=5, ncol=3, byrow=TRUE) 7 cor_x <- cor(original_x) 8 eigen(cor_x) 9 10 # $values 11 # [1] 1.8136028 0.9475442 0.2388530 12 # 13 # $vectors 14 # [,1] [,2] [,3] 15 # [1,] 0.6932939 -0.07728855 0.7164985 16 # [2,] 0.6553525 -0.34594075 -0.6714448 17 # [3,] 0.2997610 0.93506763 -0.1891874
2、svd
svd函数用来对矩阵进行奇异值分解。能直接计算出U、∑、V。程序代码如下：
1 x1 <- c(0.5, 11,100) 2 x2 <- c(0.4, 8,130) 3 x3 <- c(0.8, 9,120) 4 x4 <- c(0.2, 3,104) 5 x5 <- c(0.3, 5,114) 6 original_x <- matrix(c(x1,x2,x3,x4,x5), nrow=5, ncol=3, byrow=TRUE) 7 col1 <- (original_x[,1] - mean(original_x[,1])) / sd(original_x[,1]) 8 col2 <- (original_x[,2] - mean(original_x[,2])) / sd(original_x[,2]) 9 col3 <- (original_x[,3] - mean(original_x[,3])) / sd(original_x[,3]) 10 std_x <- matrix(c(col1,col2,col3), nrow=5) 11 svd(std_x,nu=5,nv=3) 12 13 # $d 14 # [1] 2.6934014 1.9468377 0.9774517 15 # 16 # $u 17 # [,1] [,2] [,3] [,4] [,5] 18 # [1,] -0.2316664 0.76089582 -0.40903159 0.3601178 0.26528579 19 # [2,] -0.1668699 -0.61250817 -0.56141824 0.5302241 0.02569229 20 # [3,] -0.5984366 -0.09147745 0.65686800 0.3907370 0.22215896 21 # [4,] 0.6765052 0.10549224 0.29254870 0.6515667 -0.14521774 22 # [5,] 0.3204676 -0.16240244 0.02103313 -0.1093756 0.92656385 23 # 24 # $v 25 # [,1] [,2] [,3] 26 # [1,] -0.6932939 0.07728855 0.7164985 27 # [2,] -0.6553525 0.34594075 -0.6714448 28 # [3,] -0.2997610 -0.93506763 -0.1891874 29 30 31 eigen_stdx <- eigen(t(std_x) %*% std_x, symmetric=FALSE) 32 sqrt(eigen_stdx$values) 33 # [1] 2.6934014 1.9468377 0.9774517
代码中给出的计算结果中 $d 就是 奇异值，即：∑矩阵对角线上的元素。
代码的最后三行通过计算矩阵 ATA 特征值的方式来计算得到的奇异值，与直接使用svd函数是一致的。
从奇异值可以看出，第一、第二主成分的贡献率是相当高的，所以本例从3 维 降维到2维，信息损失很少。
3、princomp
princomp函数可以用来直接提取矩阵主成分。princomp函数使用的是上面提到的通过协方差矩阵，或相关系数矩阵来，计算主成分。
1 x1 <- c(0.5, 11,100) 2 x2 <- c(0.4, 8,130) 3 x3 <- c(0.8, 9,120) 4 x4 <- c(0.2, 3,104) 5 x5 <- c(0.3, 5,114) 6 original_x <- matrix(c(x1,x2,x3,x4,x5), nrow=5, ncol=3, byrow=TRUE) 7 prncp_obj <- princomp(original_x, cor=TRUE) 8 summary(prncp_obj, loadings=TRUE) 9 # Importance of components: 10 # Comp.1 Comp.2 Comp.3 11 # Standard deviation 1.3467007 0.9734188 0.48872587 12 # Proportion of Variance 0.6045343 0.3158481 0.07961766 13 # Cumulative Proportion 0.6045343 0.9203823 1.00000000 14 # 15 # Loadings: 16 # Comp.1 Comp.2 Comp.3 17 # [1,] 0.693 0.077 0.716 18 # [2,] 0.655 -0.346 -0.671 19 # [3,] 0.300 0.935 -0.189 20 21 predict(prncp_obj) 22 # Comp.1 Comp.2 Comp.3 23 # [1,] 0.6976202 -1.6561892 -0.44699965 24 # [2,] 0.5024977 1.3332041 -0.61353148 25 # [3,] 1.8020806 0.1991126 0.71784129 26 # [4,] -2.0371699 -0.2296176 0.31970432 27 # [5,] -0.9650286 0.3534901 0.02298552
上面的代码中，使用相关系数矩阵来计算主成分，“Cumulative Proportion” 这一行（第13行）表明前两个主成分的贡献率超过了92%。
predict给出了每个样本点主成分计算的结果，可以使用结果矩阵的前两列（第一、二主成分）作为学习算法的输入。
