主成分分析PCA(Principal Component Analysis)

PCA 是一种常见的基于线性变换的数据降维方法，能够将原始数据变换为一组各维度线性无关的表示。

算法步骤

1. 构造n行m列的矩阵X；
2. 按行对矩阵进行0均值化；
3. 求出协方差矩阵 $C=\frac{1}{m}X{X}^T$ 的特征值以及对应的特征向量；
4. 按特征值大小对特征向量进行排序，之后取出前k行(期望的维度)组成矩阵P；
5. 返回$Y=PX$。

R语言代码实现PCA函数

# @param X: matrix 要求内容为可计算的数字
# @param k: 降维后矩阵的尺寸（小于X的列数）
pca <- function(X,k){
  # 行0均值化
  for(i in 1:nrow(X)) X[i,] <- X[i,]-mean(X[i,])
  # 按特征值大小对行排序后的协方差矩阵的特征向量矩阵(data.frame)，返回前n列与X的积
  with(eigen((X%*%t(X))/ncol(X)),vectors[order(values),])[1:k,]%*%X
}

算法的数学原理

原始数据

数据的行和列分别代表什么含义？

# 随机产生两组数据各100个
A <- matrix(runif(300)*20,nrow=3)
B <- matrix(runif(300)*20+15,nrow=3)

使用随机数产生两组使用矩阵表示的数据。这两个矩阵各有100列，每列表示一条记录（Record），每行表示记录中的一个字段（Attribute）同时也代表数据的维度。3行就代表三个维度。以下是部分数据的样例。

          [,1]      [,2]       [,3]      [,4]      [,5]      [,6]      [,7]      [,8]
[1,] 18.960956  1.748467  6.2212104 18.322197  9.317954  2.825602 11.268678  3.739493
[2,]  2.551108 15.809941  0.6862373  6.144451  9.101383 16.783569  0.281998  1.353761
[3,] 15.463714 14.475244 10.3760060  7.748462 11.792671  5.040896  3.564827 15.413125

使用plot3d函数可以直观的看到数据的分布状态

plot3d(c(A[1,],B[1,]),c(A[2,],B[2,]),c(A[3,],B[3,]),col = rep(c("red","blue"),each=100),xlab="x",ylab="y",zlab="z")

方差

为什么要求降维后数据的方差最大？

如果一组数据之间离得很近，代表这组数据的相似度很高，也就是说明了这组数据所包含的信息量更少。在机器学习中，数据中的信息量越大，学习效果也就相应的越好。那么应当如何去衡量一组数据的相似度呢？方差就是一种简单有效的方式。计算每个字段方差的函数如下：
$$Var(A)=\frac{1}{m}\sum _{i=1}^m(a_i-\bar{A}{)}^2$$
然而这个式子里有一个均值$\bar{A}$的存在，导致计算复杂度提升，那么有没有什么方法能够去除均值呢？

均值归零

均值归零的意义是什么？

如果我们能够让每个字段的均值为0，那么字段方差函数可以简化为如下形式：
$$Var(A)=\frac{1}{m}\sum _{i=1}^m{a}_i^2$$
均值归零的方法表示如下：
$$a_i^{\prime }=a_i-\bar{A},i=1,2,...,m$$

协方差

为什么要求协方差为0？

当我们仅使用方差最大的标准来降维时，可能会出现相同字段间数据分散，但不同字段间数据聚合的现象，也就是说不同字段间会存在线性相关性。那么字段间就会存在信息重复的情况。我们可以使用协方差函数计算两组数据之间相关性系数：
$$Cov(A,B)=\frac{1}{m}\sum _{i=1}^m(a_i-\bar{A})(b_i-\bar{B})$$
当然在行均值归零化之后，协方差函数简化为下面的形式：
$$Cov(A,B)=\frac{1}{m}\sum _{i=1}^m{a}_i{b}_i$$
相关性系数的绝对值越小，说明两组数据的相关程度越低。

协方差矩阵

如何使方差最大同时让协方差为0？

假设我们有矩阵$Y$：
$$Y=\left[\begin{array}{cccc}{a}_1& a_2& ...& a_m\\ b_1& b_2& ...& b_m\end{array}\right]$$
那么
$$D=\frac{1}{m}Y{Y}^T=\left(\begin{array}{cc}\frac{1}{m}{\sum }_{i=1}^m{a}_i^2& \frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i\\ \frac{1}{m}{\sum }_{i=1}^m{a}_i{b}_i& \frac{1}{m}{\sum }_{i=1}^m{b}_i^2\end{array}\right)=\left(\begin{array}{cc}Var(A)& Cov(A,B)\\ Cov(A,B)& Var(B)\end{array}\right)$$

Y矩阵是什么？

若D矩阵除主对角线外的元素全为0，就表示矩阵Y各字段的协方差均为0。显然，Y矩阵就是我们期望X降维后的目标矩阵。$Y=f(X)$

如何变换X才能得到Y？

我们假设Y和X之间存在一种线性变换关系f，那么
$$Y=f(X)=PX$$
其中P是一个未知的矩阵。于是Y的协方差矩阵就被转换为了
$$\begin{gathered} D=\frac{1}{m}Y{Y}^T \\ =\frac{1}{m}(PX)(PX{)}^T \\ =\frac{1}{m}PX{X}^T{P}^T \\ =P(\frac{1}{m}X{X}^T)P^T \\ =PC{P}^T \end{gathered}$$
其中C是X的协方差矩阵，显然C是已知矩阵。那么问题就从“如何转换X到Y”变成了“对矩阵C对角化”。对角化过程请参照“线性代数”相关书籍。

假设我们得到了n组特征值与特征向量，我们将特征值按大小顺序取前k个对应的特征向量组成k*m的矩阵$P=[{\xi }_1,{\xi }_2,...,{\xi }_k{]}^T$，那么
$$PC{P}^T=\Lambda =\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_k\end{array}\right)$$

得到矩阵P之后P左乘X就得到了Y。

示例

以下是三维降至二维的结果。

以下是三维直接降至一维的结果

以下是全部示例代码

# ----------------------PCA函数--------------------------
# @param X: matrix 要求内容为可计算的数字
# @param k: 降维后矩阵的尺寸（小于X的列数）
pca <- function(X,k){
  # 行0均值化
  for(i in 1:nrow(X)) X[i,] <- X[i,]-mean(X[i,])
  # 按特征值大小对行排序后的协方差矩阵的特征向量矩阵(data.frame)，返回前n列与X的积
  with(eigen((X%*%t(X))/ncol(X)),vectors[order(values),])[1:k,]%*%X
}

--------------------模拟数据--------------------
# 随机产生两组数据各100个
A <- matrix(runif(300)*20,nrow=3)
B <- matrix(runif(300)*20+15,nrow=3)

# 绘制出可拖动的3D立体分布图
if(!require(rgl)){
  install.packages("rgl")
  library(rgl)
}
plot3d(c(A[1,],B[1,]),c(A[2,],B[2,]),c(A[3,],B[3,]),col = rep(c("red","blue"),each=100),xlab="x",ylab="y",zlab="z")

# -----------------对三维数据进行降维到二维--------------

# two_dim <- pca(cbind(A,B),2)
# plot(two_dim[1,],two_dim[2,], pch = 16,col = rep(c('red','blue'),each = 100),xlab="x", ylab="y") 

# ---------------三维降至一维---------------------
# plot(pca(cbind(A,B),1),rep.int(0,200), pch = 16,col = rep(c('red','blue'),each = 100),xlab="x", ylab="y")