机器学习-无监督学习-主成分分析(PCA)-数学原理、应用以及Python和R实现

数学原理简要说明

基本原理

PCA所做的事情实际上是将多维的数据降维到一个更低的维度，但同时保证数据之间尽可能的分开，而衡量这个分散程度的标准就是方差的大小。比如假设我们有X, Y, Z三个坐标轴，构成的空间内部有一系列数据点，我们想要降维到2维，则需要找PC1、PC2这两个轴，使得当前数据点投影到这两个轴后，两个轴方向的方差最大（注意每个轴是独立查找的), 即最大投影方差。

首先，我们需要对数据进行中心化，设A为一个中心化后的数据矩阵，考虑一个轴投影，s为投影后的方差总和，v为投影方向的单位向量，长度为1, Si对应为每个A中列向量向v投影后得到的平方，则有:

$s = \sum\limits_{i=1}^{n}\frac{S_i}{n-1} = \frac{v^TAA^Tv}{n-1} = v^T\frac{AA^T}{n-1}v = v^TCv$

那么接下来我们实际上是要求S的最大值，并且其对应的v即投影方向，以及有着限制条件 $v^Tv = 1$ ，可以发现，这实际上就是一个拉格朗日乘数法的应用：

$F(v, \lambda) = v^TCv + \lambda(1-v^Tv)$

$F'_v = 2Cv - \lambda2v = 0$

可得:

$Cv = λv$

可以发现，我们要找的v，实际上就是协方差矩阵C的特征向量的方向，接着我们再来看特征值，

$s = v^TCv = v^T\lambda v = \lambda$

可以发现，向某一特征向量投影后，所对应的方差就是特征向量对应的特征值，顺便，我们知道C是实对称矩阵，所以特征向量都是相互正交的。

接下来我们也许想要知道之前每个维度对投影到现有特征向量的贡献，你可以想象现在的特征向量，实际是原来轴的Linear combination，将特征向量长度定为1，每个原始维度通过线性组合得到这个特征向量，则每个原始维度前的系数就是Loading Score。

用每个特征向量对应的特征值，即方差，除以方差也即特征值总和，可以量化投影到每个特征向量使得数据分散的程度，即可量化每个特征向量的贡献大小，即方差占比。

求解方法

根据上述，我们知道我们所要做的事情就是得到协方差矩阵的特征值、特征向量，有了特征向量我们就可以将原数据向其投影，得到降维的新数据。有两种求解方法，分别为EVD和SVD，即特征值分解和奇异值分解。

EVD

若A为一个N x N方阵，且有N个线性独立的特征向量qi(i = 1,..., N)。则有：

$A = Q\Lambda Q^{-1}$

Q为NxN，列向量为A的特征向量， $\Lambda$ 是对角阵，对角线上为A的对应特征向量的特征值。

因此我们可以对协方差矩阵C做EVD，得到相应的特征向量和特征值。

SVD

即 $A_{m × n} = U_{m × m}\Sigma_{m × n} V^T_{n × n}$

$U$ 和 $V^T$ 为正交矩阵，分别称为左奇异矩阵和右奇异矩阵，其中，∑为对角阵，对角的元素σi称为矩阵A的奇异值。

$A^TA = (U\Sigma V^T)^TU\Sigma V^T = V\Sigma^TU^TU\Sigma V^T = V\Sigma^T\Sigma V^T = V\Sigma^2 V^T$

则可以发现， $V$ 是 $A^TA$ 的特征向量按列组成的正交阵， $\Sigma^2$ 是 $A^TA$ 特征值组成的对角阵，因此可发现，A的奇异值，实际上是 $A^TA$ 特征值的平方根。

而 $协方差C = \frac{1}{n-1}A^TA$ ，因此我们对A进行SVD，Σ可以拿到协方差矩阵的特征值，V可以拿到协方差矩阵的特征向量。

SVD和EVD比较

由于一般来说，数据集的行数，即维度，都很高， $A^TA$ 这样直接做EVD的计算量是非常大的，而SVD避免了直接计算协方差矩阵，SVD 比典型的特征值分解过程更稳定，尤其是对于机器学习而言。在算法上，SVD通常应用分而治之算法，而EVD应用QR算法，分治法在一些研究中被认为比QR法更加稳定。

代码实现

为了Python和R的方便比较，统一使用iris数据集。

Python实现

嘛，就不整合成函数了，感觉反而应该更简洁清晰明了的，首先是sklearn的实现（也是基于SVD)，接着是手推numpy基于EVD，然后是SVD。
注意，sklearn中的预处理中ddof都默认为0，它基于numpy, 而pandas则默认为1，但实际上不会有太大影响，也可以调用scipy.stats中的zscore来标准化，或者手动标准化。

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import seaborn as sns
from sklearn.decomposition import PCA
from sklearn.preprocessing import LabelEncoder, StandardScaler
from scipy.stats import zscore

# 借助sklearn实现
iris = sns.load_dataset("iris")
prop = iris.iloc[:, 0:4]
label = iris.iloc[:, -1]
# 转换成因子
label_factor = LabelEncoder()
label_factor.fit(label.unique())
f_label = label_factor.transform(label)
maptable = {f: l for f, l in zip(f_label, label)}
# scale
zscore = StandardScaler()
zscore.fit(prop)
z_prop = zscore.transform(prop)  # 或 z_prop = zscore(prop, ddof=1)
# PCA
pca = PCA(n_components=2)
pca.fit(z_prop)
expv = pca.explained_variance_ratio_
plt.subplot(221)
barplot = plt.bar(["PC"+str(i) for i in range(1, len(expv)+1)], expv)
for rect in barplot:
    y = rect.get_height()
    plt.text(rect.get_x()+rect.get_width()/2, y +
             0.05, str(round(y*100, 2)), ha='center')
plt.ylim(0, 1)
pca_transform = pca.transform(z_prop)
plt.subplot(222)
for f in maptable.keys():
    plt.scatter(pca_transform[label == maptable[f]][:, 0],
                pca_transform[label == maptable[f]][:, 1], label=maptable[f])
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.legend()
plt.title("sklearn")

## 自己实现
# ----------------
# EVD
z_prop = (prop-prop.mean())/prop.std()
cov = z_prop.cov()
eigvalue, eigvector = np.linalg.eig(cov)
evd_data = z_prop @ eigvector[:, 0:2]
plt.subplot(223)
for f in maptable.keys():
    plt.scatter(evd_data[label == maptable[f]][0],
                evd_data[label == maptable[f]][1], label=maptable[f])
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.legend()
plt.title("EVD")
# SVD
u, s, vt = np.linalg.svd(z_prop)
svd_data = z_prop @ vt[0:2, :].T
plt.subplot(224)
for f in maptable.keys():
    plt.scatter(svd_data[label == maptable[f]][0],
                svd_data[label == maptable[f]][1], label=maptable[f])
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.legend()
plt.title("SVD")
plt.tight_layout()

image.png

R实现

rm(list=ls())
library(FactoMineR)
library(factoextra)
library(ggplot2)
library(gridExtra)


prop <- iris[, 1:4]
pca_recs <- PCA(prop, graph=F)
pca_recs$eig
p1 <- fviz_eig(pca_recs, addlabels=TRUE, ylim=c(0, 100), 
               geom=c("bar", "line"), barfill="orange", barcolor="white", 
               linecolor = "blue")+
    ggtitle("Variance coverage")+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Principal Components")
p2 <- fviz_pca_ind(pca_recs,
             label = "none",
             habillage = iris$Species, # color by groups
)+
    ggtitle("Plotted by factoextra")+
    theme(plot.title = element_text(hjust=0.5))+
    xlab(paste("PC1", " (", round(pca_recs$eig[, 2][1], 2), "%)", sep=""))+
    ylab(paste("PC2", " (", round(pca_recs$eig[, 2][2], 2), "%)", sep=""))
# EVD
z_prop <- scale(prop)
z_prop.cov <- cov(z_prop)
z_prop.cov.eigen <- eigen(z_prop.cov)
evd.data <- z_prop %*% z_prop.cov.eigen$vectors[, 1:2]
evd.data <- data.frame(evd.data, as.vector(iris$Species))
names(evd.data) <- c("PC1", "PC2", "Species")
p3 <- ggplot(evd.data, aes(x=PC1, y=PC2))+
    geom_point(aes(col=Species))+
    ggtitle("EVD")+
    theme(plot.title = element_text(hjust = 0.5))
# SVD
z_prop.svd <- svd(z_prop)
svd.data <- z_prop %*% z_prop.svd$v[, 1:2]
svd.data <- data.frame(svd.data, as.vector(iris$Species))
names(evd.data) <- c("PC1", "PC2", "Species")
p4 <- ggplot(evd.data,aes(x=PC1, y=PC2))+
    geom_point(aes(col=Species))+
    ggtitle("SVD")+
    theme(plot.title = element_text(hjust = 0.5))
grid.arrange(p1, p2, p3, p4, nrow=2)