SVD奇异值分解与PCA主成分分析

SVD奇异值分解

SVD的思想

Eigen::svd和 np.linalg.svd的不同之处

工程里，用一些比较简单的东西去表达一个复杂的东西，是很有用的。对应到数学，比如一个复杂的函数，我们可以用泰勒级数去模拟他，这样泰勒级数前几项，就可以近似当做这个函数的值。更自然的是三角函数的级数去模拟一个函数， 工具是傅里叶分析。在矩阵来说，我们想用一个简单的矩阵，或者说几个数字，去代表一个复杂的矩阵。奇异值就是很好的一个选择。

SVD的直观理解

奇异值分解（singular value decomposition, SVD），将矩阵分解为奇异向量（singular vector）和奇异值（singular value）。通过奇异值分解，我们会得到一些与特征分解相同类型的信息。每个实数矩阵都有一个奇异值分解。

首先先给出奇异值的物理意义：奇异值往往对应着矩阵中隐含的重要信息，且重要性和奇异值大小正相关。每个矩阵$A$都可以表示为一系列秩为1的“小矩阵”之和，而奇异值则衡量了这些“小矩阵”对于$A$的权重。

我们都知道，图片实际上对应着一个矩阵，矩阵的大小就是像素大小，比如这张图对应的矩阵阶数就是450×333，矩阵上每个元素的数值对应着像素值。我们记这个像素矩阵为A。

现在我们对矩阵A进行奇异值分解。直观上，奇异值分解将矩阵分解成若干个秩一矩阵之和，用公式表示就是：
$$A={\sigma }_1{u}_1{v}_1^T+{\sigma }_2{u}_2{v}_2^T+...+{\sigma }_r{u}_r{v}_r^T(1)$$
其中(1)式右边每一项前的系数$\sigma$就是奇异值，u和v分别表示列向量，秩一矩阵的意思是矩阵秩为1。注意到每一项$u{v}^T$都是秩为1的矩阵。我们假定奇异值满足
$${\sigma }_1\ge {\sigma }_2\ge ...\ge {\sigma }_r>0$$
（奇异值大于0是个重要的性质，但这里先别在意），如果不满足的话重新排列顺序即可，这无非是编号顺序的问题。

既然奇异值有从大到小排列的顺序，我们自然要问，如果只保留大的奇异值，舍去较小的奇异值，这样(1)式里的等式自然不再成立，那会得到怎样的矩阵——也就是图像？

当k从1不断增大时，Ak不断的逼近A。让我们回到公式
$$\begin{array}{r}A=\sigma 1u1v{1}^T+\sigma 2u2v{2}^T+...+\sigma rur{v}_r^T(1)\end{array}$$
矩阵A表示一个450×333的矩阵，需要保存450×333=149850个元素的值。等式右边u和v分别是450×1和333×1的向量，每一项有1+450+333=784个元素。如果我们要存储很多高清的图片，而又受限于存储空间的限制，在尽可能保证图像可被识别的精度的前提下，我们可以保留奇异值较大的若干项，舍去奇异值较小的项即可。例如在上面的例子中，如果我们只保留奇异值分解的前50项，则需要存储的元素为784×50=39200，和存储原始矩阵A相比，存储量仅为后者的26%。

在图像处理领域，奇异值不仅可以应用在数据压缩上，还可以对图像去噪。如果一副图像包含噪声，我们有理由相信那些较小的奇异值就是由于噪声引起的。当我们强行令这些较小的奇异值为0时，就可以去除图片中的噪声。如下是一张25×15的图像。

但往往我们只能得到如下带有噪声的图像（和无噪声图像相比，下图的部分白格子中带有灰色）：

通过奇异值分解，我们发现矩阵的奇异值从大到小分别为：14.15，4.67，3.00，0.21，……，0.05。除了前3个奇异值较大以外，其余奇异值相比之下都很小。强行令这些小奇异值为0，然后只用前3个奇异值构造新的矩阵，得到

可以明显看出噪声减少了（白格子上灰白相间的图案减少了）。

奇异值分解还广泛的用于主成分分析（Principle Component Analysis，简称PCA）和推荐系统（如Netflex的电影推荐系统）等。在这些应用领域，奇异值也有相应的意义。

SVD的几何含义

**奇异值分解的几何含义为：对于任何的一个矩阵，我们要找到一组两两正交单位向量序列，使得矩阵作用在此向量序列上后得到新的向量序列保持两两正交。**下面我们要说明的是，奇异值的几何含义为：这组变换后的新的向量序列的长度。

奇异值分解的几何本质为：我们总能发现一个正交的网格，能够被变换为另一个正交的网格。

接下来我们使用向量来解释：已知经过恰当选择的正交单位向量$v_1$和$v_2$，向量$M{v}_1$和$M{v}_2$是正交的。(注意：并不是同轴的，这一点和特征值是不一样的)

$\ast M\to$

将向量$M{v}_1$和$M{v}_2$同轴的单位向量记为$u_1$和$u_2$，将向量$M{v}_1$和$M{v}_2$的长度记为${\sigma }_1$和${\sigma }_2$，用来描述网格在这些特定的方向上被拉伸的倍数。这些值就被称为矩阵$M$的奇异值。

因此我们有
$$\begin{array}{r}M{v}_1={\sigma }_1{u}_1\\ M{v}_2={\sigma }_2{u}_2\end{array}$$
我们现在简单描述下矩阵M如何对一般的向量$x$做处理。由于向量$v_1$和$v_2$是正交单位向量，所以有
$$x=(v_1\cdot x)v_1+(v_2\cdot x)v_2$$
这意味着
$$\begin{array}{rl}Mx& =(v_1\cdot x)M{v}_1+(v_2\cdot x)M{v}_2\\ & =(v_1\cdot x){\sigma }_1{u}_1+(v_2\cdot x){\sigma }_2{u}_2\end{array}$$
我们知道，向量$v$和$x$的点乘可以用向量$v$转置乘以向量$x$表示，即
$$v\cdot x=v^{T}x$$
则
$$\begin{array}{r}Mx=u_1{\sigma }_1{v}_1^{T}x+u_2{\sigma }_2{v}_2^{T}x\\ M=u_1{\sigma }_1{v}_1^T+u_2{\sigma }_2{v}_2^T\end{array}$$
通常将其表示为
$$M=U\Sigma V^T$$
其中，矩阵$U$的列向量分别为$u_1$和$u_2$，$\Sigma$是对角矩阵，其对角元为${\sigma }_1$和${\sigma }_2$，并且矩阵$V$的列向量分别为$v_1$和$v_2$，T表示转置。

上式表明了如何将矩阵分解为三个矩阵相乘：V描述了一组正交基，U描述了值域中的一组正交基，$\Sigma$描述了V中的向量在U中被拉伸了多少。

SVD更严格的数学推导

https://zhuanlan.zhihu.com/p/26306568

求奇异值

• pytorch

import torch

A = torch.randn(5, 3)
U, S, Vh = torch.linalg.svd(A, full_matrices=False)
U.shape, S.shape, Vh.shape
torch.dist(A, U @ torch.diag(S) @ Vh)

U, S, Vh = torch.linalg.svd(A)
U.shape, S.shape, Vh.shape
torch.dist(A, U[:, :3] @ torch.diag(S) @ Vh)

A = torch.randn(7, 5, 3)
U, S, Vh = torch.linalg.svd(A, full_matrices=False)
torch.dist(A, U @ torch.diag_embed(S) @ Vh)

• numpy

import numpy as np

center = np.mean(data, axis=0)
centered_data = data - center
U, S, Vt = np.linalg.svd(centered_data)
k = 3
U_reduced = U[:, :k]
S_reduced = np.diag(S[:k])
Vt_reduced = Vt[:k, :]

# 重构原始数据
reconstructed_data = np.dot(U_reduced, np.dot(S_reduced, Vt_reduced)) + np.mean(data, axis=0)

print(f"原始数据:\n{data[:5]}")
print(f"重构的数据:\n{reconstructed_data[:5]}")
print(f"U:\n{U_reduced}")
print(f"S:\n{S}")
print(f"Vt:\n{Vt_reduced}")

SVD和矩阵特征分解的关系

花书的第二十八页：事实上，我们可以用与A相关的特征分解去解释A的奇异值分解。A的左奇异向量是$A{A}^T$的特征向量。A的右奇异向量是$A^{T}A$的特征向量。A的非零奇异值是$A^{T}A$特征值的平方根，同时也是$A{A}^T$特征值的平方根。

应用实例

数据压缩

奇异值分解可用于有效的表达数据。假设我们希望传送下面的图片，这个图片由$15\times 25$的黑白像素矩阵组成。

我们把这个图像用15×25矩阵表示，0代表黑色像素，1代表白色像素。这样的话，矩阵一共有375个元素。

如果我们对矩阵M进行奇异值分解，会发现只有三个非零奇异值。
$$\begin{array}{r}{\sigma }_1=14.72\\ {\sigma }_2=5.22\\ {\sigma }_3=3.31\end{array}$$
因此，矩阵可以表示为
$$M=u_1{\sigma }_1{v}_1^T+u_2{\sigma }_2{v}_2^T+u_3{\sigma }_3{v}_3^T$$
其中，$v_i$有三个，每一个有15个元素，$u_i$有三个，每一个有25个元素，奇异值有三个。我们仅用了3×(1+15+25)=123个数就能表示375个像素的图片了。这说明，奇异值分解能够发现矩阵中的冗余元素，并且能够消除它。

为什么只有三个非零奇异值？因为非零奇异值的数目等于矩阵的秩，在本例中，矩阵中有三个线性无关的列，即秩为3。

降噪

上一个例子说明了我们如何处理许多奇异值为零的情况。通常大的奇异值指向感兴趣的信息。举个例子，当我们把一幅图用扫描仪扫进电脑中时，会在图像中引入一些瑕疵点（通常称为噪声）。

我们用和上例相同的方式来处理：将数据用15×25的矩阵来表示，并进行奇异值分解。奇异值如下所示：
$$\begin{array}{rl}& {\sigma }_1=14.15\\ & {\sigma }_2=4.67\\ & {\sigma }_3=3.00\\ & {\sigma }_4=0.21\\ & {\sigma }_5=0.19\\ & ...\\ & {\sigma }_{15}=0.05\end{array}$$
很明显，前三个奇异值是最重要的，所以我们假设其他的都是图像中的噪声引起的，所以只保留前三项做近似
$$M\approx u_1{\sigma }_1{v}_1^T+u_2{\sigma }_2{v}_2^T+u_2{\sigma }_2{v}_2^T$$
这样就可以改善图像质量。

数据分析

噪声在手机数据的时候也会毫无例外的出现：不论设备有多好，测量值总是有一些误差在里面。如果我们还记得在矩阵中，大的奇异值指向重要的特征，则当数据收集完后，很自然的就会想到用奇异值分解来研究数据。

作为一个例子，假设我们收集到了如下数据。

我们把这些数据放到矩阵中：
$$M=\left[\begin{array}{cccccccccc}-1.03& 0.74& -0.02& 0.51& -1.31& 0.99& 0.69& -0.12& -0.72& 1.11\\ -2.23& 1.61& -0.02& 0.88& -2.39& 2.02& 1.62& -0.35& -1.67& 2.46\end{array}\right]$$

并进行奇异值分解，得到奇异值为
$$\begin{array}{r}{\sigma }_1=6.04\\ {\sigma }_2=0.22\end{array}$$
由于第一个奇异值远远要比第二个奇异值大，我们有理由认为，较小的${\sigma }_2$是数据中的噪声引起的，所以理想情况下${\sigma }_2$等于零。在这种情况下，矩阵的秩为1($y_i=k{x}_i$)意味着所有的数据都是都在一条由$u_i$定义的线上。

这个简单的例子其实是给 PCA (主成分分析：principal component analysis)开了个小头，PCA使用奇异值分解来检测数据的相关性和冗余程度。

类似的，奇异值分解还可用于检测数据中的群组类别，这就解释了为什么奇异值分解能改善Netflix的电影推荐系统。你观看过的电影的评分可让程序把你排进一个群组，这个群组中的其他人观看过的电影评分和你类似。这样就可以选择一个你所在的群组中的其他人看过的电影评分较高的电影推荐给你。

坐标系变换

奇异值分解（Singular Value Decomposition，SVD）是一种将一个矩阵分解为三个矩阵乘积的方法，其应用非常广泛。在计算机视觉和几何处理中，SVD 可以用于估计点云的旋转矩阵。

当我们使用SVD分解来获得旋转矩阵时，它实际上是将原始的变换分解为三个基本的变换：旋转、缩放和再次旋转。

在点云处理中，假设我们有两个点云，分别表示为 A 和 B。我们希望找到一个旋转矩阵 R，将点云 A 旋转到与点云 B 对齐。这个问题可以转化为求解最佳的旋转矩阵 R，使得点云 A 经过旋转 R 后与点云 B 的误差最小。

SVD 可以帮助我们找到这个最佳的旋转矩阵 R。具体步骤如下：

1. 首先，我们将点云 A 和点云 B 分别转换为 3 × N 的矩阵，其中 N 是点的数量，每一列代表一个点的坐标。

2. 然后，我们计算两个点云的质心（centroid）并将它们移到原点，使得它们的质心重合。

3. 接下来，我们计算点云 A 和点云 B 的协方差矩阵 H，并对其进行 SVD 分解。其中，SVD 分解可以将矩阵 H 分解为三个矩阵 U, S, V^T，满足 H = U · S · V^T。

4. 最后，我们将矩阵 V^T · U^T 作为旋转矩阵 R。

这样，我们就得到了一个将点云 A 旋转到与点云 B 对齐的最佳旋转矩阵 R。通过应用这个旋转矩阵，我们可以将点云 A 进行旋转，使其与点云 B 对齐，从而实现点云的配准（registration）。

协方差矩阵描述了数据集中两个或多个变量之间的线性关系和各个变量的方差。在计算机视觉和机器人学中，协方差矩阵常用于描述点云数据的分布和结构。
假设有一个点云，其中的点可以表示为
P = { p 1 , p 2 , … , p n } \mathbf{P} = \{ \mathbf{p}_1, \mathbf{p}_2, \ldots, \mathbf{p}_n \} P={p1​,p2​,…,pn​}，其中${\mathbf{p}}_i$是第 ( i ) 个点的坐标。对于这样一个点云，其协方差矩阵 $\Sigma$可以表示为：
$$\Sigma =\frac{1}{n}\sum _{i=1}^n({\mathbf{p}}_i-\bar{\mathbf{p}})({\mathbf{p}}_i-\bar{\mathbf{p}}{)}^T$$其中，$\bar{\mathbf{p}}$是点云的平均位置，定义为：$$\bar{\mathbf{p}}=\frac{1}{n}\sum _{i=1}^n{\mathbf{p}}_i$$接下来，考虑一个点云经过旋转变换，使用旋转矩阵 $\mathbf{R}$进行旋转。新的点云可以表示为 P ′ = { R p 1 , R p 2 , … , R p n } \mathbf{P}' = \{ \mathbf{R}\mathbf{p}_1, \mathbf{R}\mathbf{p}_2, \ldots, \mathbf{R}\mathbf{p}_n \} P′={Rp1​,Rp2​,…,Rpn​}。对于旋转后的点云，其协方差矩阵 ${\Sigma }^{\prime }$可以表示为：${\Sigma }^{\prime }=\frac{1}{n}{\sum }_{i=1}^n(\mathbf{R}{\mathbf{p}}_i-\overline{\mathbf{Rp}})(\mathbf{R}{\mathbf{p}}_i-\overline{\mathbf{Rp}}{)}^T$这里，$\overline{\mathbf{Rp}}$是旋转后点云的平均位置，定义为：$$\overline{\mathbf{Rp}}=\frac{1}{n}\sum _{i=1}^n\mathbf{R}{\mathbf{p}}_i$$从上述公式可以看出，旋转矩阵 $\mathbf{R}$会影响点云的协方差矩阵。但是，旋转不会改变点云的主要方向（主成分）。因此，通过计算点云的主成分，可以得到旋转后的点云的主要方向，从而实现旋转矩阵的估计。

当我们考虑点云数据时，我们通常想要了解数据的分布以及数据之间的关系。协方差矩阵是用来描述这种数据之间关系的一种工具。它告诉我们不同维度之间的相关性以及每个维度的方差。

假设我们有一个二维点云，其中每个点都有两个坐标：$(x,y)$。我们可以将每个点看作一个二维向量，表示为 ${\mathbf{p}}_i=\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]$。点云中共有 $n$ 个这样的点。

点云的协方差矩阵描述了点在各个维度上的变化程度以及各个维度之间的相关性。对于二维情况，协方差矩阵如下所示：

$$\Sigma =\left[\begin{array}{cc}\text{Var}(x)& \text{Cov}(x,y)\\ \text{Cov}(x,y)& \text{Var}(y)\end{array}\right]$$

其中，$\text{Var}(x)$和$\text{Var}(y)$分别是 ( x ) 和 ( y ) 维度上的方差，$\text{Cov}(x,y)$是 $x$和$y$之间的协方差。

现在假设我们对这个点云进行了旋转变换，旋转矩阵为 $\mathbf{R}$，那么新的点云可以表示为 P ′ = { R p 1 , R p 2 , … , R p n } \mathbf{P}' = \{ \mathbf{R}\mathbf{p}_1, \mathbf{R}\mathbf{p}_2, \ldots, \mathbf{R}\mathbf{p}_n \} P′={Rp1​,Rp2​,…,Rpn​}

对于旋转后的点云，其协方差矩阵为：

$${\Sigma }^{\prime }=\left[\begin{array}{cc}\text{Var}(x^{\prime })& \text{Cov}(x^{\prime },y^{\prime })\\ \text{Cov}(x^{\prime },y^{\prime })& \text{Var}(y^{\prime })\end{array}\right]$$

其中，( (x’, y’) ) 是旋转后点的新坐标。由于 ${\mathbf{P}}^{\prime }=\mathbf{RP}$，我们可以将$(x^{\prime },y^{\prime })$表示为 $\mathbf{R}$与原始坐标$(x,y)$的乘积。

在实际应用中，我们可以通过计算旋转后的点云的协方差矩阵来推断旋转矩阵$\mathbf{R}$。这样做的一个常见方法是使用主成分分析（PCA），它可以找到协方差矩阵的特征向量，这些特征向量代表了数据中的主要方向，而特征向量对应的特征值表示了数据在这些方向上的方差。通过比较旋转后的点云的主成分与原始点云的主成分，我们可以推断出旋转矩阵的近似值。

理论上，通过协方差计算变换矩阵的过程是可能的，但它可能会涉及到一些复杂的数学推导和优化算法。一种常见的方法是使用主成分分析（PCA），它可以从协方差矩阵中提取主要的特征向量，从而推断出数据的主要方向，从而得到变换矩阵。
下面是一个简化的步骤：

1. 计算原始点云数据的协方差矩阵。
2. 对协方差矩阵进行特征值分解，得到特征向量和对应的特征值
3. 选择最大的几个特征值对应的特征向量，它们表示数据的主要方向。
4. 构建变换矩阵，使得这些特征向量成为变换后数据的新坐标系的基向量。
   这种方法可以使得变换后的数据在主要方向上具有最大的方差，从而保留了原始数据的主要结构。具体实现时，你可以使用数学库（如NumPy或SciPy）来计算协方差矩阵和进行特征值分解。然后，根据选定的主要特征向量构建变换矩阵。需要注意的是，实际应用中可能需要考虑一些额外的因素，例如数据的噪声和误差，以及如何选择主要的特征向量。

理解为什么特征向量（eigenvectors）可以表示旋转矩阵需要一些线性代数的知识。让我们从简单的情况开始。
考虑一个二维空间中的旋转。对于一个二维空间中的向量，旋转矩阵会将这个向量旋转到新的位置。而特征向量是在变换（比如旋转）后方向不变的向量。这意味着，如果我们有一个二维空间中的向量，它是一个特征向量，那么当我们对这个向量应用旋转矩阵时，它仍然保持在同一条线上，只是可能会被拉伸或收缩。这就是特征向量的一个重要性质：它们在变换后的方向保持不变。现在考虑一个二维空间的旋转矩阵。对于这样的矩阵，它有两个特征向量，分别对应于旋转的主轴。这两个特征向量描述了旋转的方向。当我们对这些特征向量应用旋转矩阵时，它们保持不变，因为它们是特征向量。这就是为什么特征向量可以表示旋转矩阵的原因。在更高维度的情况下，旋转矩阵可能有更多的特征向量，每个特征向量对应于旋转的一个主轴。这些特征向量描述了旋转的方向，而特征值则描述了在每个主轴上的缩放因子（即旋转的量）。因此，特征向量和特征值可以完整描述一个旋转矩阵的变换效果。总而言之，特征向量可以表示旋转矩阵，因为它们描述了旋转矩阵的旋转方向，而特征值描述了旋转的量。
在三维空间中，旋转矩阵可以通过特征向量和特征值来表示。一个三维空间中的旋转矩阵将三维向量旋转到新的位置，而特征向量是在旋转后方向不变的向量。这意味着，如果我们有一个三维空间中的向量，它是一个特征向量，那么当我们对这个向量应用旋转矩阵时，它仍然保持在同一条线上，只是可能会被拉伸或收缩。
对于一个三维旋转矩阵，它有三个特征向量，每个特征向量对应于旋转的一个主轴。这三个特征向量描述了旋转的方向。当我们对这些特征向量应用旋转矩阵时，它们保持不变，因为它们是特征向量。而特征值则描述了在每个主轴上的缩放因子（即旋转的量）。因此，在三维空间中，通过特征向量和特征值，我们可以完整描述一个旋转矩阵的变换效果。特征向量表示旋转的方向，而特征值描述了旋转的量。

在计算协方差矩阵的特征值分解或者奇异值分解时，通常会对数据进行中心化处理。这是因为中心化后的数据更容易处理，并且可以更好地捕捉数据的变化情况，从而更准确地找到数据集的主要方向。

中心化数据的过程是通过减去数据的均值来实现的，这样做可以将数据的中心移到原点附近，使得数据的平均值为零。这对于协方差矩阵的计算是很重要的，因为协方差矩阵的计算涉及到数据点与均值之间的差异。

具体来说，对点云数据进行中心化处理的原因包括：

1. 简化计算：中心化后的数据具有零均值，这样可以简化协方差矩阵的计算，因为协方差矩阵的元素可以通过简单的内积来计算。
2. 减小数值误差：通过减去均值，可以减小数据中的数值范围，从而有助于提高数值稳定性和减少数值误差。
3. 保留数据结构：中心化后的数据保留了原始数据的结构，但更容易处理。这有助于确保计算出的特征向量和特征值更好地反映数据的本质。

因此，在进行特征值分解或奇异值分解之前，通常会对数据进行中心化处理，以确保更好的结果。

PCA主成分分析

在实际生产生活中，我们所获得的数据集在特征上往往具有很高的维度，对高维度的数据进行处理时消耗的时间很大，并且过多的特征变量也会妨碍查找规律的建立。如何在最大程度上保留数据集的信息量的前提下进行数据维度的降低，是我们需要解决的问题。

对数据进行降维有以下优点：

• 使得数据集更易使用
• 降低很多算法的计算开销
• 去除噪声
• 使得结果易懂

降维技术作为数据预处理的一部分，即可使用在监督学习中也能够使用在非监督学习中。其中主成分分析PCA应用最为广泛，本文也将详细介绍PCA。

主成分分析（Principal Component Analysis）：一种统计方法，它对多变量表示数据点集合寻找尽可能少的正交矢量表征数据信息特征。

PCA既是一种压缩数据的方式，也是一种学习数据表示的无监督学习方法。

直觉分析

直观上而言，思路应该聚焦在这两方面：

一是，要考虑去除掉特征之间的相关性，想法是创造另一组新的特征来描述样本，并且新的特征必须彼此之间不相关。

二是，在新的彼此无关的特征集中，舍弃掉不重要的特征，保留较少的特征，实现数据的特征维度降维，保持尽量少的信息损失。

期望与方差

我们的最终目标是分析如何提取数据的主成分，如何对手头的数据进行降维，以便后续的进一步分析。往往问题的切入点就是数据各个维度之间的关系以及数据的整体分布。因此，我们有必要先花点功夫，来梳理一下如何对数据的整体分布情况进行描述。

期望衡量的是一组变量$X$取值分布的平均值，一般记作：$E[X]$，反映的是不同数据集的整体水平。比如，在一次期末考试中，一班的平均成绩是90分，二班的平均成绩是85分，那么从这两个班级成绩的均值来看，就反映出一班的成绩在总体上要优于二班。

方差这个概念大家也不陌生，方差的定义是$V[X]=E[(X-\mu )]$，其中，$\mu =E[X]$，它反映的是一组数据的离散程度。通俗的说就是：对于一组数据，其方差越大，数据的分布就越发散，方差越小，数据的分布就越集中。在一组样本集的方差计算中，我们采用$\frac{1}{n-1}{\sum }_1^n(x_i-\mu {)}^2$作为样本$X$的方差估计。

期望（mean）和方差（variance）是描述随机变量分布的重要统计量。

1. 期望：随机变量的期望是对所有可能取值的加权平均，其表示了随机变量的中心位置。对于离散型随机变量 (X)，其期望 (E[X]) 计算公式为：

$$E[X]=\sum _i{x}_i\cdot P(X=x_i)$$

其中，(x_i) 是随机变量 (X) 可能取到的值，(P(X = x_i)) 是随机变量 (X) 取到 (x_i) 的概率。

对于连续型随机变量 (X)，其期望 (E[X]) 计算公式为：

$$E[X]={\int }_{-\infty }^{\infty }x\cdot f(x)dx$$

其中，(f(x)) 是随机变量 (X) 的概率密度函数。

2. 方差：随机变量的方差是随机变量与其期望之差的平方的期望，其表示了随机变量的离散程度。方差越大，随机变量的取值越分散。方差的计算公式为：
   $$Var(X)=E[(X-E[X]{)}^2]$$
   对于离散型随机变量 (X)，方差 (Var(X)) 的计算公式为：

$$Var(X)=\sum _i(x_i-E[X]{)}^2\cdot P(X=x_i)$$

对于连续型随机变量 (X)，方差 (Var(X)) 的计算公式为：

$$Var(X)={\int }_{-\infty }^{\infty }(x-E[X]{)}^2\cdot f(x)dx$$

其中，(f(x)) 是随机变量 (X) 的概率密度函数。

协方差与协方差矩阵

大多数人都会有这样一种感觉，一般而言一个人如果数学成绩好，他的物理成绩大概率也会不错，反之亦然。但是数学成绩与英语成绩的却没有那么相关，至少说不像是数学成绩和物理成绩关联那样密切。
协方差是用来衡量两个随机变量之间的关系的统计量，它描述了两个随机变量一起变化的趋势。协方差可以通过以下公式计算：

1. 离散型随机变量：

$$Cov(X,Y)=E[(X-E[X])(Y-E[Y])]$$

其中，(E[X]) 和 (E[Y]) 分别是随机变量 (X) 和 (Y) 的期望。

2. 连续型随机变量：

$$Cov(X,Y)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }(x-E[X])(y-E[Y])\cdot f(x,y)dxdy$$

其中，(f(x, y)) 是随机变量 (X) 和 (Y) 的联合概率密度函数。

协方差的正负号表示了两个变量之间的线性关系方向，正值表示正相关，负值表示负相关。而协方差的绝对值大小表示了两个变量之间的线性关系的强度，绝对值越大表示两个变量之间的关系越强。
我们为了具体量化这种现象，引入协方差这个概念，对于随机变量$X$和$Y$，二者的协方差定义如下：
$$Cov[X,Y]=E[(X-\mu )(Y-\nu )]$$
其中，$\mu$和$\nu$分别是$X$和$Y$的期望。其实，我们可以对比观察一下上面讲过的方差的公式：
$$V[X]=E[(X-\mu {)}^2]=E[(X-\mu )(X-\mu )]$$
不难发现可以将方差看做是协方差的一种特殊形式。

通过上面协方差的定义式，我们可以用通俗的话来概况这里面的含义：

当$X$和$Y$的协方差为正的时候，表示当$X$增大的时候，$Y$也倾向于随之增大；协方差为负的时候，表示当$X$增大，$Y$倾向于减小；当协方差为0的时候，表示当$X$增大时，$Y$没有明显的增大或减小的倾向，二者是相互独立的。

我们可以引入一个协方差矩阵，将一组变量$X_1,X_2,X_3$两两之间的协方差用矩阵的形式统一进行表达：
$$\left[\begin{array}{ccc}V[X_1]& Cov[X_1,X_2]& Cov[X_1,X_3]\\ & & \\ Cov[X_2,X_1]& V[X_2]& Cov[X_2,X_3]\\ & & \\ Cov[X_3,X_1]& Cov[X_3,X_2]& V[X_3]\end{array}\right]$$
协方差矩阵的第$i$行，第$j$列表示$X_i$和$X_j$之间的协方差，而对角线上表示随机变量自身的方差。我们从定义中可以看出，$Cov[X_i,X_j]$和$Cov[X_j,X_i]$显然是相等的。因此，协方差矩阵是一个对称矩阵，结合上一讲关于对称矩阵的介绍，我们知道协方差矩阵拥有很多优良的性质，后面的应用中会派上大用场。

通过协方差矩阵可以看出变量之间的关系，主要有以下几个方面：

1. 对角线元素（方差）：
   对角线上的元素表示各个变量的方差，它们描述了各个变量自身的变化程度。方差越大，表示变量的取值波动越大。

2. 非对角线元素（协方差）：
   非对角线上的元素表示各个变量之间的协方差，它们描述了不同变量之间的线性关系。协方差的正负号表示了变量之间的关系方向，正值表示正相关，负值表示负相关，而协方差的绝对值大小表示了变量之间的线性关系的强度。

3. 协方差的大小与相关性：
   协方差的绝对值越大，表示变量之间的线性关系越强，即一个变量的值的变化会引起另一个变量值的变化。而协方差接近于零，则表示两个变量之间基本上没有线性关系。

4. 协方差矩阵的特征值和特征向量：
   协方差矩阵的特征值和特征向量可以提供更多关于变量之间关系的信息。特征向量表示了协方差矩阵的主要方向，而特征值表示了在这些方向上的变量之间的方差大小。通过分析特征值和特征向量，可以确定变量之间的主要关系方向和重要性。

综合来看，通过分析协方差矩阵的各个元素以及其特征值和特征向量，可以全面地了解变量之间的关系，从而进行数据分析、模型建立和决策制定。

主成分分析法降维的思路

直接去掉特征？太盲目

能不能直接挑选一个特征属性，然后将其去掉，通过这种方法简单的实现数据降维？直觉告诉我们，似乎是不行的。这种方法最为关键的一点就是没有考虑到特征属性之间彼此紧密耦合的相关关系。这样独立的对各个特征属性进行分析，是很盲目的。

那么具体该如何做？我们还是举一个非常直观的例子来说明：

假设我们研究的对象有两个特征属性$X$和$Y$，对5个样本进行数据采样的结果如下：

取值	样本1	样本2	样本3	样本4	样本5
$X$	2	2	4	8	4
$Y$	2	6	6	8	8

我们的目标是对其降维，只用一维特征来表示每个样本。我们首先将其绘制在二维平面图中进行整体观察，如下图所示：

在5个样本的特征数据分布图中，我们能挖掘出几个重要的规律：

首先，我们来探究这两个特征属性的数据分布和数据相关情况，该数据的协方差矩阵为：
$$\left[\begin{array}{cc}6& 4\\ 4& 6\end{array}\right]$$
结合特征变量$X$和$Y$的图像分布以及我们的观察发现，5个样本的特征$X$和特征$Y$呈现出正相关性，数据彼此之间存在着影响。

其次，我们如果真要硬来，直接去掉特征$X$或特征$Y$来实现特征降维，可行吗？显然是不行的，如图3所示，如果我们直接去掉特征$X$，那么样本2和样本3，以及样本4和样本5，就变得无法区分了。但是实质上，在降维前的原始数据中，他们显然是有明显区别的。直接去掉特征$Y$也是一样，这种简单粗暴的降维方法显然是不可取的，他忽视了数据中的内在结构关系，并且带来了非常明显的信息损失。

从解除特征之间的相关性

我们应该如何调整我们的降维策略？基于上面的现象，直观上而言，思路应该聚焦在这两方面：

一是，要考虑去除掉特征之间的相关性，想法是创造另一组新的特征来描述样本，并且新的特征必须彼此之间不相关。

二是，在新的彼此无关的特征集中，舍弃掉不重要的特征，保留较少的特征，实现数据的特征维度降维，保持尽量少的信息损失。

剖析主成分分析的细节

构造彼此无关的新特征

首先我们来看第一点，我们尝试用2个新的特征来对样本进行描述，那么为了让这两个新特征满足彼此无关的要求，就得让这两个新特征的协方差为0，构成的协方差矩阵是一个对角矩阵。

而目前我们所使用的原始特征$X$和$Y$的协方差不为0，其协方差矩阵是一个普通的对称矩阵。

看到这里，你是不是突然觉得眼前一亮，我们的工作不就是要将当前的对称矩阵变换到对角矩阵么？大的方向的确如此，我们就沿着这个方向，一步一步从头到尾来实施一遍主成分分析(PCA)的具体方法。

我们首先回顾一下协方差的定义公式：
$$Cov(X,Y)=\frac{1}{n-1}\sum _1^n(x_i-\mu )(y_i-\nu )$$
感觉公式还不够简练，我们发现如果变量$X$和变量$Y$的均值为0，那么协方差的式子就会变得更加简单。

这个其实好办，我们将$X$中的每一个变量都减去他们的均值$\mu$，同样将$Y$中的每一个变量都减去他们的均值$\nu$，这样经过零均值化处理后，特征$X$和特征$Y$的平均值都变为了0。很显然，这样的处理并不会改变方差与协方差的值，因为数据的离散程度并未发生改变，同时，从公式的表达上来看也非常清楚。

经过零均值处理后，$X$和$Y$的协方差矩阵就可以写作：
$$\begin{gathered} \left[\begin{array}{cc}V[X]& Cov[X,Y]\\ Cov[Y,X]& V[Y]\end{array}\right] \\ =\left[\begin{array}{cc}\frac{1}{n-1}{\sum }_1^n{x}_i^2& \frac{1}{n-1}{\sum }_1^n{x}_i{y}_i\\ \frac{1}{n-1}{\sum }_1^n{y}_i{x}_i& \frac{1}{n-1}{\sum }_1^n{y}_i^2\end{array}\right] \\ =\frac{1}{n-1}\left[\begin{array}{cc}{\sum }_1^n{x}_i^2& {\sum }_1^n{x}_i{y}_i\\ {\sum }_1^n{y}_i{x}_i& {\sum }_1^n{y}_i^2\end{array}\right] \end{gathered}$$
我们看一下
$$\left[\begin{array}{cc}{\sum }_1^n{x}_i^2& {\sum }_1^n{x}_i{y}_i\\ {\sum }_1^n{y}_i{x}_i& {\sum }_1^n{y}_i^2\end{array}\right]$$
这个矩阵，不难发现它就是零均值化后的样本矩阵
$$A=\left[\begin{array}{ccccc}{x}_1& x_2& x_3& ...& x_n\\ y_1& y_2& y_3& ...& y_n\end{array}\right]$$
与自身转置矩阵$A^T$相乘的结果，即：协方差矩阵
$$\begin{gathered} C=A{A}^T=\left[\begin{array}{cc}{\sum }_1^n{x}_i^2& {\sum }_1^n{x}_i{y}_i\\ {\sum }_1^n{y}_i{x}_i& {\sum }_1^n{y}_i^2\end{array}\right] \\ =\left[\begin{array}{ccccc}{x}_1& x_2& x_3& ...& x_n\\ y_1& y_2& y_3& ...& y_n\end{array}\right]\left[\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\\ x_3& y_3\\ .& .\\ x_n& y_n\end{array}\right] \end{gathered}$$
一般而言，$X$和$Y$满足线性无关，所以依据我们之前介绍的知识，协方差矩阵是对称的、正定的满秩方阵。

接下来，我们就想把各样本的特征$X$和特征$Y$的取值用另外一组基来表示，由于要求在新的基下表示的新特征彼此无关（即新特征的协方差矩阵是对角化的，非对角线的协方差元素为零），那么必然的，新选择的两个基必须彼此正交。

我们选择一组新的基，$p_1$和$p_2$，他们的模长为1，彼此之间满足标准正交，前面我们学过，在此条件下，向量$a=\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]$与向量$p_i$的点积就表示向量$a$在$p_i$方向上的投影，同时由于$p_i$是单位向量，那么点积的结果$p_i^{T}a$就代表了向量$a$在基向量$p_i$上的投影坐标值。

那么向量$a$在由标准正交向量$p_1$和$p_2$构成的新基底上的坐标即为$\left[\begin{array}{c}{p}_1^{T}a\\ p_2^{T}a\end{array}\right]$，如果我们令$P=\left[\begin{array}{cc}{p}_1& p_2\end{array}\right]$，那么可以进一步将其表示为：$P^{T}a$。

梳理一下这里所做的工作：
$$A=\left[\begin{array}{ccccc}{x}_1& x_2& x_3& ...& x_n\\ y_1& y_2& y_3& ...& y_n\end{array}\right]$$
，就是$n$个样本的特征$X$和特征$Y$的原始采样值。而
$$P^{T}A=\left[\begin{array}{ccccc}{p}_1^T{a}_1& p_1^T{a}_2& p_1^T{a}_3& ...& p_1^T{a}_n\\ p_2^T{a}_1& p_2^T{a}_2& p_2^T{a}_3& ...& p_2^T{a}_n\end{array}\right]$$
（其中$a_i=\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]$）。描述的就是这$n$个样本在新构建的特征下的取值。

此时，我们来实现第一个目标：

即，试图让第一个新特征$\left[p_1^T{a}_1,p_1^T{a}_2,...,p_1^T{a}_n\right]$和第二个新特征$\left[p_2^T{a}_1,p_2^T{a}_2,...,p_2^T{a}_n\right]$满足线性无关。换句话说，就是要求这两个特征的协方差为0，即让协方差矩阵为一个对角矩阵。
$$\begin{array}{rl}D& =\frac{1}{n-1}{P}^{T}A(P^{T}A{)}^T\\ & =\frac{1}{n-1}{P}^{T}A{A}^{T}P\\ & =\frac{1}{n-1}{P}^{T}CP\\ & =\frac{1}{n-1}\Lambda \end{array}$$
，这就转化成了我们再熟悉不过的问题了：去寻找让矩阵$C$对角化的$P$矩阵，而$D=\frac{1}{n-1}\Lambda =\frac{1}{n-1}{P}^{T}CP$，所以当$P$是$C$矩阵的特征向量时，采用新特征值的样本矩阵$P^{T}A$的协方差矩阵$D$为单位阵，即不同新特征之间的协方差为零，新特征之间线性无关。这里还有一个好的前提条件：$C$是对称矩阵、正定矩阵，它保证了矩阵一定能够对角化，且特征值全都为正。这里需要理解的是，让协方差为0是为了投射后的空间的基之间线性无关，方差保留前几个最大的值是为了保留原始数据所在空间的最大信息。

因为线性代数中有这么一个定理：对称矩阵$C$一定可以得到由一组标准正交特征向量构成的特征矩阵$P$，即$\Lambda =P^{T}CP$（其中，$P=[p_1,p_2,...,p_n]$，$q_1$到$q_n$就是协方差矩阵$C$的$n$个标准正交特征向量）。由此，我们再通过$P^{T}A$就能计算得出彼此线性无关的新特征。

PCA的严谨证明

上述的过程可能会让人有点莫名其妙，凭什么正交矩阵的基向量上的方差就一定最大了？没道理呀，那现在就证明一下。

PCA的推导有两种出发点：

• 基于最小投影距离（具体证明见PCA两种推导方式）

• 基于最大投影方差？

但这两种的推导最后都是殊途同归。为了方便起见，选第二种思路最大投影方差来证明：

样本在新坐标系中的投影为$P^{T}A$，其协方差矩阵为
$$\begin{array}{rl}D& =\frac{1}{n-1}{P}^{T}A(P^{T}A{)}^T\\ & =\frac{1}{n-1}{P}^{T}A{A}^{T}P\\ & =\frac{1}{n-1}{P}^{T}CP\end{array}$$
使投影方差和最大，就是最大化协方差矩阵的迹：
$$\begin{array}{rl}{\mathrm{argmax}}_P& tr(P^{T}CP)\\ s.t.& P^{T}P=I\end{array}$$
利用拉格朗日乘子法:
$$J(W)=tr(P^{T}CP)+\lambda (I-P^{T}P)$$
对$P$求导，结果即为
$$CP+C^{T}P-2\lambda P=0(C=C^T)$$
可得
$$CP=\lambda P$$
则矩阵$P$为矩阵$C$的特征向量。

新得到的特征该如何取舍

在上面的一系列操作中，我们构造了两个新特征。原始特征的方向分别是$x$轴正方向和$y$轴的正方向，两个原始特征彼此相关。而新构造的两个特征，方向垂直，在这两个方向上构造的新特征，协方差为0，线性无关。

那么，我们将原始特征的采样值变换到两个新选取的投影基上，就得到新的一组特征取值。

我们在此基础上再进行特征维度的降维，让二维数据变成一维。由于这两个新的特征彼此无关，我们可以放心大胆的保留一个，去掉一个。

具体保留哪一个，我们的判定标准是方差，方差越大，表明这个特征里数据分布的离散程度就越大，特征所包含的信息量就越大；反之，如果特征里数据的方差小，分布集中，则表明其包含的信息量就小。那么，我们自然选择保留信息量大的那个特征了。

这里，主特征值对应的特征向量方向上的特征对应的方差最大（即${\sigma }_1^2=\frac{1}{n-1}{\lambda }_1$），而次特征值对应的特征向量方向上的特征对应的方差是${\sigma }_2^2=\frac{1}{n-1}{\lambda }_2$。我们决定保留主特征值对应的特征向量方向上的特征取值。注意，这里方差和特征值的关系是比例关系，即${\sigma }^2=\frac{1}{n-1}\lambda$，这是因为
$$D=\frac{1}{n-1}\left[\begin{array}{cc}{\sum }_1^n{x}_i^2& {\sum }_1^n{x}_i{y}_i\\ {\sum }_1^n{y}_i{x}_i& {\sum }_1^n{y}_i^2\end{array}\right]=\frac{1}{n-1}{P}^{T}CP=\frac{1}{n-1}\Lambda$$
此时，我们完成了主成分的提取，样本如果用一维数据来描述，最终的取值就是$p_1^{T}A$，从而实现了数据的降维。

最终所保留的就是以$p_1$为基向量的坐标值，如下图所示：

衡量信息的损失

那么，我们怎么衡量数据降维过程中的信息损失？或者反过来说，我们保留下了多少信息？这里介绍主成分贡献率的概念，上面的例子里，原本一共有两个特征，我们保留第一个作为主成分，我们用方差来衡量主成分贡献率：

主成分贡献率为： $\frac{{\lambda }_1}{{\lambda }_1+{\lambda }_2}$。

曾经在项目中做过一个尝试，对于512维度进行不同维度压缩后保留的信息：

512->256: 0.9998875105696291
512->128: 0.8964985919468982
512->64 : 0.6758438520261455

由此可见，我们在数据降维的过程中，如果将512维度变为256维度，数据压缩率为50%，但保留了99.99%的信息，效果还是不错的。

推广到n个特征的降维

这里我们推广到一般情况，如果我们有$n$个特征，如何进行数据降维。

这里的思路是完全一样的。

第一步：针对采样得到的$m$个样本，我们得到了一个$n\times m$的样本数据矩阵$A$；

第二步：对每一个特征的取值，进行零均值化处理，如果这些特征不在一个数量级上，我们还应该将其除以标准差$\sigma$。

第三步：利用预处理后的矩阵$A$，求$n$个特征的协方差矩阵：$C=\frac{1}{n-1}A{A}^T$

第四步：对协方差矩阵进行对角化处理，求得矩阵$C$的$n$个标准正交特征向量，并按照对应的特征值的大小依次排列。

第五步：按照事先规定的主成分贡献度，提取满足该数值要求的前$k$个新构造的特征作为主成分，构造成数据压缩矩阵：$P=[p_1,p_2,...,p_k]$。

第六步：通过矩阵相乘$P^{T}A$实现前$k$个主成分的提取，将$n$维特征压缩到$k$维。

实例分析

原始数据零均值化，求协方差

原始数据为
$$\begin{array}{rl}X& =\left[\begin{array}{ccccc}2& 2& 4& 8& 4\end{array}\right]\\ Y& =\left[\begin{array}{ccccc}2& 6& 6& 8& 8\end{array}\right]\end{array}$$
即
$$A=\left[\begin{array}{ccccc}2& 2& 4& 8& 4\\ 2& 6& 6& 8& 8\end{array}\right]$$
零均值化为
$$A=\left[\begin{array}{ccccc}-2& -2& 0& 4& 0\\ -4& 0& 0& 2& 2\end{array}\right]$$
其协方差矩阵为（注：$A$零均值化了）
$$\begin{array}{rl}C& =\frac{1}{n-1}\left[\begin{array}{cc}{\sum }_1^n{x}_i^2& {\sum }_1^n{x}_i{y}_i\\ {\sum }_1^n{y}_i{x}_i& {\sum }_1^n{y}_i^2\end{array}\right]\\ & =\frac{1}{n-1}\left[\begin{array}{ccccc}{x}_1& x_2& x_3& ...& x_n\\ y_1& y_2& y_3& ...& y_n\end{array}\right]\left[\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\\ x_3& y_3\\ .& .\\ x_n& y_n\end{array}\right]\\ & =\frac{1}{n-1}A{A}^T\\ & =\left[\begin{array}{cc}6& 4\\ 4& 6\end{array}\right]\end{array}$$
用代码表示为

import numpy as np

# 原始数据
x = [2, 2, 4, 8, 4]
y = [2, 6, 6, 8, 8]

a = x - np.mean(x)
b = y - np.mean(y)

# 零均值化
S = np.vstack((x, y))
print(S)
# [[-2. -2.  0.  4.  0.]
#  [-4.  0.  0.  2.  2.]]

# 协方差
print(np.cov(S))
# [[ 6.  4.]
#  [ 4.  6.]]
# 原始数据的协方差
print(np.cov(np.vstack((a, b))))
# [[6. 4.]
#  [4. 6.]]

在上面的过程中，我们对原始数据两个特征的采样值进行了零均值化处理，同时也验证了：零均值化的处理后，协方差矩阵确实不变。

对样本的协方差矩阵求特征值和特征向量

对于协方差矩阵$C$，其特征值$\Lambda$和特征向量$P$

特征值为$[10,2]$，对应的特征向量分别为$\left[\begin{array}{c}0.707\\ 0.707\end{array}\right]$和$\left[\begin{array}{c}-0.707\\ 0.707\end{array}\right]$。

用代码表示为

import numpy as np

C = np.array([[6, 4],
              [4, 6]])

evalue, evector = np.linalg.eig(C)
print(evalue) # 特征值
# [ 10.+0.j   2.+0.j]
print(evector)  # 特征向量
#[[ 0.70710678 -0.70710678]
# [ 0.70710678  0.70710678]]

在新坐标系上降维

在上面的一系列操作中，我们构造了两个新特征。原始特征的方向分别是$x$轴正方向和$y$轴的正方向，两个原始特征彼此相关。而新构造的两个特征，方向分别是$\left[\begin{array}{c}0.707\\ 0.707\end{array}\right]$和$\left[\begin{array}{c}-0.707\\ 0.707\end{array}\right]$，在这两个方向上构造的新特征，协方差为0，线性无关。 那么，我们将原始特征的采样值变换到两个新选取的投影基上，就得到新的一组特征取值$P^{T}A$：
$$\begin{array}{rl}{P}^{T}A& ={\left[\begin{array}{cc}0.707& -0.707\\ 0.707& 0.707\end{array}\right]}^T\left[\begin{array}{ccccc}-2& -2& 0& 4& 0\\ -4& 0& 0& 2& 2\end{array}\right]\\ & =\left[\begin{array}{ccccc}-4.242& -1.414& 0& 4.242& 1.414\\ -1.414& 1.414& 0& -1.414& 1.414\end{array}\right]\end{array}$$
我们得到的两个结果，就是在原始特征和新构建特征上分别得到的值。 我们在此基础上再进行特征维度的降维，让二维数据变成一维。由于这两个新的特征彼此无关，我们可以放心大胆的保留一个，去掉一个。

具体保留哪一个，我们的判定标准是方差，方差越大，表明这个特征里数据分布的离散程度就越大，特征所包含的信息量就越大；反之，如果特征里数据的方差小，分布集中，则表明其包含的信息量就小。那么，我们自然选择保留信息量大的那个特征了。

这里，$\left[\begin{array}{c}0.707\\ 0.707\end{array}\right]$方向上的特征对应的方差是2.5，而$\left[\begin{array}{c}-0.707\\ 0.707\end{array}\right]$方向上的特征对应的方差是0.5。我们决定保留$\left[\begin{array}{c}0.707\\ 0.707\end{array}\right]$方向上的特征取值。

用代码表示为

import numpy as np

x = [2,2,4,8,4]
y = [2,6,6,8,8]

x = x - np.mean(x)
y = y - np.mean(y)

A = np.vstack((x,y))

p_1 = [0.707, 0.707]
p_2 = [-0.707, 0.707]

P = np.vstack((p_1,p_2))

# 零均值化后的A
print(A)
# [[-2. -2.  0.  4.  0.]
#  [-4.  0.  0.  2.  2.]]
# 新坐标系P^T下的数据特征P^TA
print(np.dot(P,A))
# [[-4.242 -1.414  0.     4.242  1.414]
#  [-1.414  1.414  0.    -1.414  1.414]]

主成分贡献率

上面的例子里，原本一共有两个特征，我们保留第一个作为主成分，我们用方差来衡量主成分贡献率： 主成分贡献率为：
$$\frac{{\lambda }_1}{{\lambda }_1+{\lambda }_2}=\frac{2.5}{2.5+0.5}=83.3$$
由此可见，我们在数据降维的过程中，数据压缩率为50%，但保留了83.3%的信息，感觉效果还是不错的。

PCA的优缺点

作为一个非监督学习的降维方法，它只需要特征值分解，就可以对数据进行压缩，去噪。因此在实际场景应用很广泛。

PCA算法的主要优点有：

• 仅仅需要以方差衡量信息量，不受数据集以外的因素影响。

• 各主成分之间正交，可消除原始数据成分间的相互影响的因素。

• 计算方法简单，主要运算是特征值分解，易于实现。

PCA算法的主要缺点有：

• 主成分各个特征维度的含义具有一定的模糊性，不如原始样本特征的解释性强。
• 方差小的非主成分也可能含有对样本差异的重要信息，因降维丢弃可能对后续数据处理有影响。
• PCA只能对数据进行线性变换，这对于一些线性不可分的数据是不利的。为了解决PCA只能进行线性变换的问题，Schölkopf, Bernhard在1998年提出了Kernel PCA，能够对数据进行非线性变换。
• PCA的结果容易受到每一维数据的大小的影响，如果我们对每一维数据乘以一个不同的权重因子之后再进行PCA降维，得到的结果可能与直接进行PCA降维得到的结果相差比较大。
• PCA要求数据每一维的均值都是0，在将原始数据的每一维的均值都变成0时可能会丢失掉一些信息。

PCA的python代码

#!/usr/bin/bash
# -*-coding:utf-8 -*-
import numpy as np

def pca(X, n_components):
    # 1. 数据预处理：标准化
    X_std = (X - np.mean(X, axis=0)) / np.std(X, axis=0)
    
    # 2. 计算协方差矩阵
    cov_matrix = np.cov(X_std, rowvar=False)
    
    # 3. 计算特征值和特征向量
    eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
    
    # 4. 选择主成分
    idx = eigenvalues.argsort()[::-1]
    top_eigenvectors = eigenvectors[:, idx[:n_components]]
    
    # 5. 投影数据
    projected_data = np.dot(X_std, top_eigenvectors)
    
    return projected_data

# 示例数据
data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

# PCA 分析
n_components = 2
projected_data = pca(data, n_components)
print("降维后的数据：")
print(projected_data)
# [[-2.12132034e+00  1.68306266e-16]
#  [ 0.00000000e+00  0.00000000e+00]
#  [ 2.12132034e+00 -1.68306266e-16]]

参考资料

• 如何通俗易懂地讲解什么是 PCA 主成分分析？石溪的回答

• PCA两种推导方式

• PCA主成分分析学习总结

• The Matrix Cookbook page13

• 理解主成分分析 (PCA)

• hinton对pca的理解

• Dimensionality Reduction——PCA原理篇

• Dimensionality Reduction——PCA实现篇

• A Tutorial on Principal Component Analysis

• https://zhuanlan.zhihu.com/p/26306568

• 奇异值的物理意义是什么？-知乎

• [We Recommend a Singular Value Decomposition（Feature Column from the AMS）

• A Singularly Valuable Decomposition: The SVD of a Matrix