<统计>主分量分析 / 主成分分析

主成分分析(PCA)是获取高维数据的一系列技术之一，它使用变量之间的依赖关系将其表示为更易处理、更低维的形式，而不会丢失太多信息。它已被广泛用于数据压缩和去噪。然而，它的整个数学过程有时对用户来说是模棱两可的。

在这篇文章中，我想从数学上讨论PCA的整个过程，包括PCA的投影和重建，并提供了大部分的推导和证明。在文章的最后，我从头开始实现了PCA投影和重建。读完这篇文章后，PCA中应该不再有黑匣子了。

在线性代数中，正交矩阵是列和行都是正交单位向量(正交向量)的实方矩阵。

根据可逆矩阵的定义，这意味着矩阵$A$是可逆的，并且$A^{-1}=A^{\top}$。

由于$\det(A)\neq 0$，矩阵$A$是可逆的。我们在正交矩阵定义的两侧将$A^{-1}$相乘。

类似地，如果一个复方阵$A$的转置共轭$A^{\dagger}$也是它的逆，则它是酉阵。

由于共轭性质和$\overline{A}=A$因为$A$是实值矩阵，

设$\λ\in\mathbb{C}$是对称矩阵$A$的特征值。$Av=\lambda v$和$v\neq 0$。我们将$v^{\dagger}$($v^{\dagger}=\overline{v}^{\top}$)乘到两边，并且由于$A^{\top}=A$和我们刚刚导出的属性$A\overline{v}=\overline{\lambda}\overline{v}$，

如果$x^{\dagger}A x\geq 0$对$x\in{C}^n$，则$n\x\n$对称矩阵$A$被定义为半正定的，如果$x^{\dagger}A x\geq 0$对$x\in\mathbb{C}^n$.。

因为$x^{\dagger}$x\geq$对于$x\in\mathbb{C}^n$，假设$x$是$A$的特征向量，并且$Ax=\lambda x$其中$x\neq 0$，

因为$x^{\dagger}x$必须是实数，并且$x^{\dagger}x>；0$，所以我们有$\lambda\geq 0$。

协方差矩阵是半正定的。这意味着协方差矩阵的特征值是非负的。

协方差必须是半正定的证明可以在我上一篇关于多元高斯和协方差矩阵的文章中找到。

$m\x n$矩阵$A$的奇异值$\sigma_1$，$\sigma_2$，$\cdots$，$\sigma_r$是相关Gram矩阵$K=A^{\dagger}A$的非负特征值的平方根$\sigma_i=\sqrt{\lambda_i}$。$K$对应的特征向量称为$A$的奇异向量。

请注意，关联的Gram矩阵$K=A^{\dagger}A$是实数和对称的，因此$K$的特征值都是实数。

因此，$K=A^{\dagger}A$的所有特征值都是非负的，并且它们都有对应的奇异值$A$。

在线性代数中，奇异值分解(SVD)是实矩阵或复矩阵的因式分解，它通过极分解的推广将平方正规矩阵的特征分解推广到任意$m×n$矩阵。

具体地说，$m\x n$实矩阵或复矩阵$M$的奇异值分解是$U\Sigma V^{\dagger}$形式的因式分解，其中$U$是$m\x m$实矩阵或复酉矩阵，$\Sigma$是对角线上具有非负实数的$m\x n$矩形对角矩阵，$V$是$n\x n$实矩阵或复酉矩阵。

$\sigma$的对角条目$\sigma_{i}=\sigma_{ii}$称为$M$的奇异值。非零奇异值的个数等于$M$的秩。

我们将跳过为什么每个矩阵都有奇异值分解的证明以及奇异值分解的算法。

请注意，$\Sigma^{\dagger}\Sigma$是正方形对角线矩阵。根据特征值和特征向量的定义，$Sigma^{\dagger}\Sigma$的对角值(包括零点)是$A^{\dagger}A$的特征值。$V$的所有列都是$A^{\dagger}A$的相应特征向量。

请注意，$\Sigma\Sigma^{\dagger}$是正方形对角线矩阵。根据特征值和特征向量的定义，$\Sigma\Sigma^{\dagger}$的对角值(包括零点)是$A^{\dagger}$的特征值。$U$的所有列都是$A^{\dagger}$的相应特征向量。

我们从$p$维向量开始，并希望通过向下投影到$q$维子空间(其中$q\leq p$)来总结它们。我们的总结是将原始向量投影到横跨该子空间的主轴$q$方向上。

给定一个数据集$X\in\mathbb{R}^{n\Times p}$，其行是$0\leq i\leq n$($\sum_{i=0}^{n-1}x_{i}=0$)的居中数据向量$x_i\in\mathbb{R}^p$($|w|=1$)，如果我们有一个单位向量$w\in\mathbb{R}^p$($|w|=1$

数据向量$x_i$在$w$上的投影向量是$\langx_i，w\rangw$。

残差是从数据向量$x_i$到$w$的距离，是向量$x_i-\langx_i，w\rangw$的长度。

让我们检查一下剩余平方$|x_i-\langle x_i，w\rangle w|^2$是什么。

投影的优化目标是最小化均方误差$\text{mse}(W)$，它是残差平方和的平均值。

记住方差和期望值之间的关系，$\mathbb{V}(X)=\mathbb{E}(X^2)-\mathbb{E}(X)^2$，我们有。

因此，$\frac{1}{n}\sum_{i=0}^{n-1}\langx_i，w\rang^2$实际上是数据向量与投影单位向量$\mathbb{V}\BIG(\langx，w\Rangle\Big)$之间的内积的方差。

由于$\mathbb{E}\BIG(\langx，x\rangelBIG)$不依赖于$w$，最小化剩余平方和$\text{mse}(W)$的均值等价于最大化数据向量与投影单位向量之间的内积的方差。

数据向量和投影单位向量之间的内积的方差也可以使用数据集矩阵$X$来表示。

如果我们将数据向量投影到$q$正交单位向量$w_0$，$w_1$，$\cdots$，$w_{q-1}$。这些单位向量以列的形式组合成矩阵$W\in\mathbb{R}^{p\x q}$。然后我们的优化目标是最大化所有这些投影向量的残差平方和。

假设我们将数据投影到$p$正交单位向量$w\in\mathbb{R}^{p}$。这些单位向量按列组合成矩阵$W\in\mathbb{R}^{p\x p}$，$W$是正交矩阵。

让我们先尝试最大化一个$\mathscr{L}(W)=w^{\top}X^{\top}X w$，约束是$w^{\top}w=1$。我们将使用拉格朗日乘子$\lambda$来帮助我们解决优化问题。

根据定义，$λ$和$w$是矩阵$X^{\top}X$的特征值和特征向量。

由于$X^{TOP}X$是实对称矩阵且是半正定的，因此$X^{TOP}X$的所有特征值都是非负的。

$X^{\top}X$具有精确的$p$特征值$\lambda_0，\lambda_1，\cdots，\lambda_{p-1}$及其对应的特征向量$w_0，w_1，\cdots，w_{p-1}$。矩阵$W$由特征向量$w_0，w_1，\cdots，w_{p-1}$组成，恰好是$\mathscr{L}(W)$的全局最大解。$\mathscr{L}(W)$的全局最大值是非负特征值的和。

因为这些特征值是非负的，并且它们对全局最大值的贡献可以根据它们的值来排序。

我们把与最大特征值对应的特征向量称为第一主轴，将与第二大特征值对应的特征向量称为第二主轴，依此类推。

假设我们将数据投影到$q$正交单位向量$w\in\mathbb{R}^{p}$和$q\leq p$。这些单位向量以列的形式组合成矩阵$W\in\mathbb{R}^{p\x q}$。最优化问题变成了最大化下面的方程。

然而，解决方案是显而易见的。我们只需挑选矩阵$X^{\top}X$的$Q$最大特征值及其对应的特征向量。

如果我们已经计算了矩阵$X^{\top}X$的特征值和特征向量，我们就可以对任意数量的主分量进行投影，而不需要做大量的计算。

简而言之，PCA没有什么特别之处，只是计算矩阵$X^{\top}X$的特征值和特征向量。要将数据矩阵$X\in\mathbb{R}^{n\Times p}$投影到$q$正交单位向量，只需取与$q$最大特征值对应的$q$对应的特征向量，将它们组合为矩阵$W\in\mathbb{R}^{p\Times q}$的列，$XW$是投影降维数据集。给定任何新的数据矩阵$A\in\mathbb{R}^{m\x p}$，类似地，投影降维数据集将是$AW$。

使用计算机直接从矩阵$X^{\TOP}X$计算特征值和特征向量有时会导致数值不稳定，由于计算机的精度有限，我不在这里举例说明这一点。

人们发现用奇异值分解(SVD)计算矩阵$X^{TOP}X$的特征值和特征向量具有更好的数值稳定性。

正如我在“前提条件”的“范数矩阵的奇异值分解”一节中所证明的那样，特征向量就是矩阵$X$的奇异值分解中矩阵$V$的列，相应的特征值就是矩阵$\Sigma^{\top}\Sigma$的对角值。

特征值是$\Sigma_{n\Times p}^{\top}\Sigma_{n\Times p}$的对角值，特征向量是$V_{p\Times p}$的列。

主成分分析是将数据矩阵从高维投影到低维的过程。然而，PCA重建是一个将数据从低维投影到高维的过程，这是一个与PCA相反的过程。

假设我们有一个数据矩阵$X\in\mathbb{R}^{n\次p}$，PCA后的投影数据矩阵是$X^{\Prime}\in\mathbb{R}^{n\次q}$，其中$X^{\Prime}=XW$，投影矩阵$W\in\mathbb{R}^{p\次q}$，我们的优化目标是找到投影矩阵$Z\in\mathbb{R}^{q\×p}$，使得重构矩阵$X^{\素数}=X^{\素数}Z=X^{\素数}Z=XWZ$与原始矩阵$X$，$\lVert X-X^{\素数}\rVert_2^{2}$之间的重构误差最小。

结果表明，当$Z=W^{\top}=V{：，：q}$时，重构误差最小。

既然我们已经讨论了这么多关于PCA的基本原理，让我们使用Numpy和Scipy从头开始实现PCA，并将结果与使用Scikit-Learning的PCA进行比较。

将numpy导入为np导入sklearn.datets，skLearning。从scipy import linalg导入时间类ManualPCA(Object)分解：def__init__(self，n_component=None，Backend="；Svd"；，Verbose=false)：self。N_Components=n_Components Self。BACKEND=BACKEND IF BACKEND！="；SvD"；AND BACKEND！="；EIG"；：引发异常("；不支持的PCA后端{}！"；。格式化(自身。后端))自我。详细=详细定义SvD_FIT(SELF，X)：START_TIME=时间。时间()自我。MU=NP。平均值(X，轴=0)X=X自身。Mu#self.V的列是X^T，X，X的本征向量。你，赛尔夫。S，Self。VH=线性。SVD(X)自身。V=自我。VH。Transspose()end_time=时间。如果是SELF，则时间()。Verbose==True：打印("；SVD拟合时间：{：.6f}秒"；。FORMAT(END_TIME-START_TIME))def EIG_FIT(SELF，X)：START_TIME=TIME。时间()自我。MU=NP。平均值(X，轴=0)X=X自身。µ协方差=Np。点(X.。Transspose()，X)#可以使用EIGH而不是EIG作为对称矩阵，以获得更快的速度。W，赛尔夫。V=linalg。EIG(协方差)IDX_ORDER=NP。ArSort(自我分类)。W)[：：-1]自我。W=自我。W[IDX_ORDER]自身。V=自我。V[：，IDX_ORDER]SELF。VH=SELF。五.。Transspose()end_time=时间。如果是SELF，则时间()。Verbose==True：打印("；特征匹配时间：{：.6f}秒"；。Format(End_Time-Start_Time))def Fit(Self，X)：If Self。Backend=="；Svd"；：Self。SVD_FIT(X)ELIF SELF。Backend=="；EIG"；：Self。EIG_FIT(X)ELSE：引发异常("；不支持的PCA后端{}！"；。格式化(自身。后端))def变换(self，X)：X=X-self。如果是自己的话。N_Components==NONE：返回NP。点(X，SELF.。V)否则：返回NP。点(X，SELF.。V[：，：Self.。N_Components])def Inverse_Transform(self，X)：If self。N_Components==NONE：返回NP。点(X，SELF.。VH)+SELF。MUELSE：返回NP。点(X，SELF.。VH[：Self.。N_组件])+自身。Mu def main()：x=skLearning。数据集。Load_iris()。DATA n_Components=None#n_Components=2 SCRICKIT_PCA=skLearning。腐烂。PCA(n_Components=n_Components)SCRICKIT_PCA。拟合(X)X_Projected=SCISKIT_PCA。Transform(X)X_Projected_Inversed=SCISKIT_PCA。Inverse_Transform(X_Projected)断言NP。Allclose(X，X_Projected_Inversed)==True，"；Scikit PCA有问题！"；SvD_PCA=ManualPCA(n_Components=n_Components，Backend="；SvD"；，Verbose=True)SvD_PCA。拟合(X)X_Projected=SvD_PCA。Transform(X)X_Projected_Inversed=SvD_PCA。Inverse_Transform(X_Projected)断言NP。Allclose(X，X_Projected_Inversed)==True，"；手动SVD PCA有问题！"；EIG_PCA=ManualPCA(n_Components=n_Components，Backend="；EIG"；，Verbose=True)EIG_PCA。拟合(X)X_Projected=EIG_PCA。Transform(X)X_Projected_Inversed=EIG_PCA。Inverse_Transform(X_Projected)断言NP。Allclose(X，X_PROPECTED_INVERSED)==True，"；手动特征PCA有问题！"；if__name__=="；__main__"；：main()。

请注意，对于任何特征，数据集矩阵$X$必须始终居中，否则PCA将无效。