三维点云课程—PCA介绍

三维点云课程---PCA介绍

三维点云课程---PCA介绍
- 1. 什么是PCA
- 2.知识铺垫
- - 2.1 SVD分解（奇异值分解）
  - 2.1 谱定理
  - 2.2 Rayleigh商
- 3.PCA的数学推导过程
- 4.PCA的总结
- 5.PCA注意点

1. 什么是PCA

PCA中文名是主成分分析，物理意义是将这些数据点投影到一个非常有特征性的方向上，每个数据点在这个向量上的投影就是主成分

2.知识铺垫

2.1 SVD分解（奇异值分解）

其中U和VTV^TVT 都是正交矩阵（旋转矩阵 U∗UT=UT∗U=IU*U^T=U^T*U=IU∗UT=UT∗U=I），∑为非负对角阵，里面的特征值由大到小进行排列，这些特征值所对应的特征向量就是描述这个矩阵变化的。通过特征值分解（m×mm \times mm×m）可以得到特征值（表示这个特征有多重要），特征向量 （表示这个特征是什么）。
但是对于$ m \times n$时，特征值分解就不行了，这时需要SVD分解了。
Am×n=Um×m×Σm×n×Vn×nTA_{m \times n}=U_{m \times m} \times \Sigma_{m \times n} \times V^T_{n \times n} Am×n=Um×m×Σm×n×Vn×nT

SVD具有以下性质

V就是ATAA^TAATA的特征向量，同理U就是AATAA^TAAT的特征向量
证明如下：

A=UΣVTAT=VΣTUTATA=VΣTUTUΣVT=VΣ2VT=VΛVTA = U\Sigma {V^T} \\ {A^T} = V{\Sigma ^T}{U^T} \\ {A^T}A = V{\Sigma ^T}{U^T}U\Sigma {V^T} = V{\Sigma ^2}{V^T} = V\Lambda {V^T} A=UΣVTAT=VΣTUTATA=VΣTUTUΣVT=VΣ2VT=VΛVT

任意的矩阵A都可以分解成三个矩阵
V表示原始域的标准正交基
U表示经过M变换后的新标准正交基
Σ\SigmaΣ表示了V中的向量与U中相对应向量之间的比例（伸缩）关系，每个奇异值会按从大到小排好顺序，值越大代表该维度重要性越高。

SVD简单的实例
分解矩阵A,其定义如下:
A=(102101)A = \begin{pmatrix} 1&0\\ 2&1\\ 0&1 \end{pmatrix} A=⎝⎛120011⎠⎞
首先计算ATAA^{T} AATA和AATA A^TAAT
ATA=(120011)(102101)(5222){A^T}A = \begin{pmatrix} 1&2&0\\ 0&1&1 \end{pmatrix}\begin{pmatrix} 1&0\\ 2&1\\ 0&1 \end{pmatrix}\begin{pmatrix} 5&2\\ 2&2 \end{pmatrix} ATA=(102101)⎝⎛120011⎠⎞(5222)
AAT=(102101)(120011)=(120251011)A{A^T} =\begin{pmatrix} 1&0\\ 2&1\\ 0&1 \end{pmatrix} \begin{pmatrix} 1&2&0\\ 0&1&1 \end{pmatrix} = \begin{pmatrix} 1&2&0\\ 2&5&1\\ 0&1&1 \end{pmatrix} AAT=⎝⎛120011⎠⎞(102101)=⎝⎛120251011⎠⎞

计算ATAA^TAATA的特征值和特征向量
λ1=6;v1=(2/51/5)λ2=1;v2=(−1/52/5){\lambda _1} = 6;{\rm{ }}{v_1} = \left( \begin{array}{l} 2/\sqrt 5 \\ 1/\sqrt 5 \end{array} \right)\\ {\lambda _2} = 1;{\rm{ }}{v_2} = \left( \begin{array}{l} -1/\sqrt 5 \\ 2/\sqrt 5 \end{array} \right)\\ λ1=6;v1=(2/51/5)λ2=1;v2=(−1/52/5)

计算AATAA^TAAT的特征值和特征向量
λ1=6;u1=(2/305/301/30)λ2=1;u2=(1/50−2/5)λ3=0;u3=(2/61/6−1/6){\lambda _1} = 6;{\rm{ }}{u_1} = \begin{pmatrix} 2/\sqrt 30 \\ 5/\sqrt 30 \\ 1/\sqrt 30 \end{pmatrix}\\ {\lambda _2} = 1;{\rm{ }}{u_2} = \begin{pmatrix} 1/\sqrt 5 \\ {\quad 0}\\ -2/\sqrt 5 \end{pmatrix} \\ {\lambda _3} = 0;{\rm{ }}{u_3} = \begin{pmatrix} 2/\sqrt 6 \\ 1/\sqrt 6\\ -1/\sqrt 6 \end{pmatrix} λ1=6;u1=⎝⎛2/305/301/30⎠⎞λ2=1;u2=⎝⎛1/50−2/5⎠⎞λ3=0;u3=⎝⎛2/61/6−1/6⎠⎞

最终得到A的奇异值分解
A=UΣVT=(2/301/52/65/3001/61/30−2/5−1/6)(600100)(2/5−1/51/52/5)A = U\Sigma {V^T} = \begin{pmatrix} {2/\sqrt {30} }&{1/\sqrt 5 }&{2/\sqrt 6 }\\ {5/\sqrt {30} }&0&{1/\sqrt 6 }\\ {1/\sqrt {30} }&{ - 2/\sqrt 5 }&{ - 1/\sqrt 6 } \end{pmatrix} \begin{pmatrix} {\sqrt 6 }&0\\ 0&1\\ 0&0 \end{pmatrix}\begin{pmatrix} {2/\sqrt 5 }&{ - 1/\sqrt 5 }\\ {1/\sqrt 5 }&{2/\sqrt 5 } \end{pmatrix} A=UΣVT=⎝⎛2/305/301/301/50−2/52/61/6−1/6⎠⎞⎝⎛600010⎠⎞(2/51/5−1/52/5)

参考奇异值分解(SVD)原理及实例分析

2.1 谱定理

定义：设A是n×nn \times nn×n的对称矩阵，即A∈Rn×nA \in R^{n \times n}A∈Rn×n 。并且λi∈R\lambda_i \in Rλi∈R,i=1,2,…,n为R的特征值。那么存在一系列的正交向量ui∈Rnu_i \in R_nui∈Rn,i=1,2,…,n,使得Aui=λiuiAu_i=\lambda_i u_iAui=λiui。等价地，存在一个正交矩阵U=[u1,...,un]U=\left[\begin{matrix} u_1, & ..., & u_n\end{matrix} \right]U=[u1,...,un] ,使得
A=UΛUT=∑i=1nλiuiuiT,Λ=diag(λ1,...,λn)A = U\Lambda {U^T} = \sum\limits_{i = 1}^n {{\lambda _i}{u_i}{u_i}^T{\rm{ }},\Lambda {\rm{ = diag(}}{\lambda _1}{\rm{,}}...{\rm{,}}{\lambda _n}{\rm{)}}} A=UΛUT=i=1∑nλiuiuiT,Λ=diag(λ1,...,λn)

2.2 Rayleigh商

对于一个对称矩阵A，那么
λmin(A)≤xTAxxTx≤λmax(A){\lambda _{min}}(A) \le \frac{{{x^T}Ax}}{{{x^T}x}} \le {\lambda _{max}}(A) λmin(A)≤xTxxTAx≤λmax(A)

其中λmin(A){\lambda _{min}}(A)λmin(A) ，λmax(A){\lambda _{max}}(A)λmax(A)分别表示A矩阵的最小和最大的特征值。其实通过A的SVD分解可以得知，A其实是由两个旋转矩阵U,V 和一个对角矩阵Σ\SigmaΣ构成，那么决定A的最大值和最小值就和旋转矩阵没有关系了，只与对角矩阵有关，拉伸的大小取决于A的最大和最小特征值大小。

数学严谨证明：
通过谱定理可以得知，对称矩阵A可以分解为A=UΛUTA = U\Lambda {U^T}A=UΛUT,那么
xTAx=xTUΛUTx=x‾TΛx‾=∑i=1nλix‾i2{x^T}Ax = {x^T}U\Lambda {U^T}x = {\overline x ^T}\Lambda \overline x = \sum\limits_{i = 1}^n {{\lambda _i}\overline x _i^2} xTAx=xTUΛUTx=xTΛx=i=1∑nλixi2
其中x‾=UTx{\overline x=U^Tx}x=UTx,可以得到
λmin∑i=1nx‾i2≤∑i=1nλix‾i2≤λmax∑i=1nx‾i2{\lambda _{min}}\sum\limits_{i = 1}^n {\overline x _i^2} \le \sum\limits_{i = 1}^n {{\lambda _i}\overline x _i^2} \le {\lambda _{max}}\sum\limits_{i = 1}^n {\overline x _i^2} λmini=1∑nxi2≤i=1∑nλixi2≤λmaxi=1∑nxi2
其实，对称矩阵U并不改变任何向量的范数
∑i=1nxi2=xTx=xTUUTx=(UTx)T(UTx)=x‾Tx‾=∑i=1nx‾i2\sum\limits_{i = 1}^n {x_i^2} = {x^T}x = {x^T}U{U^T}x = {({U^T}x)^T}({U^T}x) = {\overline x ^T}\overline x = \sum\limits_{i = 1}^n {\overline x _i^2} i=1∑nxi2=xTx=xTUUTx=(UTx)T(UTx)=xTx=i=1∑nxi2
那么联合上式，可以得到
(λminxTx=λmin∑i=1nxˉi2=λmin∑i=1nxi2)≤(∑i=1nλixˉi2=xTAx)≤(λmax∑i=1nxˉi2=λmax∑i=1nxˉi2=λmax∑i=1nxi2=λmaxxTx)进而化简得到：λminxTx≤xTAx≤λmaxxTx({\lambda _{min}}{x^T}x = {\lambda _{min}}\sum\limits_{i = 1}^n {\bar x_i^2} = {\lambda _{min}}\sum\limits_{i = 1}^n {x_i^2} ) \le (\sum\limits_{i = 1}^n {{\lambda _i}\bar x_i^2 = {x^T}Ax} ) \\ \le ({\lambda _{max}}\sum\limits_{i = 1}^n {\bar x_i^2 = } {\lambda _{max}}\sum\limits_{i = 1}^n {\bar x_i^2} = {\lambda _{max}}\sum\limits_{i = 1}^n {x_i^2} = {\lambda _{max}}{x^T}x) \\ 进而化简得到：{\lambda _{min}}{x^T}x \le {x^T}Ax \le {\lambda _{max}}{x^T}x (λminxTx=λmini=1∑nxˉi2=λmini=1∑nxi2)≤(i=1∑nλixˉi2=xTAx)≤(λmaxi=1∑nxˉi2=λmaxi=1∑nxˉi2=λmaxi=1∑nxi2=λmaxxTx)进而化简得到：λminxTx≤xTAx≤λmaxxTx

3.PCA的数学推导过程

在推导之前，先来看一下输入和输出是什么。在三维点云中，输入的时n个点。输出的k个主向量，数学记法：
Input: xi∈Rn,i=1,2,...,nx_i \in R^n ,i=1,2,...,nxi∈Rn,i=1,2,...,n
**Output:**主向量z1,z2,...,zk∈Rn,k≤nz_1,z_2,...,z_k \in R^n,k \le nz1,z2,...,zk∈Rn,k≤n
目前存在一个问题，点云的第一主方向怎么确定，第二主方向怎么确定，第三，第四呢
第一主方向：使投影数据在该方向上的方差为最大的方向，即投影的点在那个方向上分布的比较散
第二主方向：从数据点删除最重要的部分，即数据点减去投影。查找剩下数据的最重要组成部分
第三主方向：依此类推
PCA的证明：
由上面的主方向的确定可知，需要用到数据的方差，那么
1.将数据标准化为零平均值，注意这里可以理解x~1\widetilde x_1x1为数
X~=[x~1,...,x~m],x~i=xi−x‾,i=1,...,mx‾=1m∑i=1mxi\widetilde X = [{\widetilde x_1},...,{\widetilde x_m}],\widetilde x_i=x_i-\overline x,i=1,...,m \quad \quad \overline x=\frac{1}{m}\sum\nolimits_{i = 1}^m {{x_i}} X=[x1,...,xm],xi=xi−x,i=1,...,mx=m1∑i=1mxi
2.假设存在一个主方向z,使得PCA在此方向上的投影方差最大,αi\alpha_iαi为投影
αi=x~iTz,i=1,2,...,m∣∣z∣∣2=1{\alpha _i} = {\widetilde x_i}^Tz ,i=1,2,...,m \quad \quad ||z|{|_2} = 1 αi=xiTz,i=1,2,...,m∣∣z∣∣2=1
3.协方差为
1m∑i=1mαi2=1m∑i=1mzTx~ix~iTz=1mzTX~X~Tz\frac{1}{m}\sum\limits_{i = 1}^m {\alpha _i^2} = \frac{1}{m}\sum\limits_{i = 1}^m {{z^T}{\widetilde x_i}{\widetilde x_i}^Tz} = \frac{1}{m}{z^T} \widetilde X{\widetilde X^T}z m1i=1∑mαi2=m1i=1∑mzTxixiTz=m1zTXXTz
4.因此最大化上式
max⁡z∈RnzT(X~X~T)z,s.t:∣∣z∣∣2=1\mathop {\max }\limits_{z \in {R^n}} {z^T}(\widetilde X{\widetilde X^T})z,s.t:||z|{|_2} = 1 z∈RnmaxzT(XXT)z,s.t:∣∣z∣∣2=1
5.根据Rayleigh商，可以得到
λmin(X~X~T)≤(zTX~X~TzzTz=zTX~X~Tz)≤λmax(X~X~T){\lambda _{min}}(\widetilde X{\widetilde X^T}) \le (\frac{{{z^T}\widetilde X{\widetilde X^T}z}}{{{z^T}z}} = {z^T}\widetilde X{\widetilde X^T}z) \le {\lambda _{max}}(\widetilde X{\widetilde X^T}) λmin(XXT)≤(zTzzTXXTz=zTXXTz)≤λmax(XXT)
令H=X~X~TH=\widetilde X{\widetilde X^T}H=XXT,易知H为对称矩阵，根据谱定理可知，X~=UrΣVrT=∑i=1σiuiviT\widetilde X={U_r}\Sigma {V_r}^T=\sum\limits_{i = 1} {{\sigma _i}{u_i}v_i^T}X=UrΣVrT=i=1∑σiuiviT关于为什么成立，详见附录证明
H=X~X~T=UrΣVrT×VrΣTUrT=UrΣ2UrH =\widetilde X{ \widetilde X^T} = {U_r}\Sigma {V_r}^T \times {V_r}{\Sigma ^T}U_r^T = {U_r}{\Sigma ^2}{U_r} H=XXT=UrΣVrT×VrΣTUrT=UrΣ2Ur
综上所述，第一主方向为$ z_1=u_1,,,u_1是是是U_r$的第一列

6.对于第二主方向证明如下
x~i(1)=x~i−u1(u1Tx~i),i=1,...,mX~(1)=[x~i(1),...,x~m(1)]=(In−u1u1T)X~{\widetilde x_i}^{(1)} = {\widetilde x_i} - {u_1}({u_1}^T{\widetilde x_i}),i = 1,...,m\\ {\widetilde X^{(1)}} = [{\widetilde x_i}^{(1)},...,\widetilde x_m^{(1)}] = ({I_n} - {u_1}u_1^T)\widetilde X xi(1)=xi−u1(u1Txi),i=1,...,mX(1)=[xi(1),...,xm(1)]=(In−u1u1T)X
那么由于u是一个正交矩阵，u1T∗u1=E,u1T∗u2=0u^T_1 *u_1=E,u^T_1 *u_2=0u1T∗u1=E,u1T∗u2=0
X(1)=∑i=1rσiuiviT−(u1u1T)∑i=1rσiuiviT=∑i=1rσiuiviT−∑i=1rσiu1u1TuiviT=∑i=1rσiuiviT−σ1u1v1T=∑i=2rσiuiviT\begin{array}{l} {X^{(1)}} = \sum\limits_{i = 1}^r {{\sigma _i}{u_i}v_i^T - ({u_1}u_1^T)\sum\limits_{i = 1}^r {{\sigma _i}{u_i}v_i^T} } \\ \quad \quad {\rm{ = }}\sum\limits_{i = 1}^r {{\sigma _i}{u_i}v_i^T} - \sum\limits_{i = 1}^r {{\sigma _i}{u_1}u_1^T{u_i}v_i^T} \\ \quad \quad {\rm{ = }}\sum\limits_{i = 1}^r {{\sigma _i}{u_i}v_i^T} - {\sigma _1}{u_1}{v_1}^T\\ \quad \quad {\rm{ = }}\sum\limits_{i = 2}^r {{\sigma _i}{u_i}v_i^T} \end{array} X(1)=i=1∑rσiuiviT−(u1u1T)i=1∑rσiuiviT=i=1∑rσiuiviT−i=1∑rσiu1u1TuiviT=i=1∑rσiuiviT−σ1u1v1T=i=2∑rσiuiviT
故第二主方向为$ z_2=u_2,,,u_2是是是U_r$的第一列

那么第三、第四主方向依次类推

4.PCA的总结

在三维点云中，算法输入为$ X_{m \times n}$
1.按列计算数据集X的均值XmeanX_{mean}Xmean,然后令Xnew=X−XmeanX_{new}=X-X_{mean}Xnew=X−Xmean
2.求解矩阵XnewX_{new}Xnew的协方差矩阵，并将其记为Cov
3.计算协方差矩阵Cov的特征值和特征向量
4.将特征值按照从大到小进行排序，选择其中最大的k个，然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵WnxkW_{nxk}Wnxk
5.计算XnewWX_{newW}XnewW,即将数据集XnewX_{new}Xnew投影到选取的特征向量上，这样就得到了我们需要的已经降维的数据集XnewWX_{newW}XnewW

5.PCA注意点

需要强调一点的是协方差矩阵求解的方式,因为矩阵中的数据按行排列与按列排列求出的协方差矩阵是不同的,上述关于PCA的推导都是基于矩阵H是按行进行排列的,因为X~=[x~1,...,x~m]\widetilde X = [{\widetilde x_1},...,{\widetilde x_m}]X=[x1,...,xm]，以这种方式进行存储数据，数据的维度为3×m3 \times m3×m 。但是在写程序的过程中，在numpy数据是以列的形式进行排列的，数据的维度为m×3m \times 3m×3 ,那么相应的协方差的求法就和推导有些不同。
补充说明，在numpy中通过covX=np.cov(X,rowvar=0)计算协方差矩阵，通过rowvar进行列与行的切换。如果rowvar为True（默认值），则每行代表一个变量X，另一个行为变量Y。否则，每列代表一个变量X，另一个列为变量Y,显然在程序中使用rowvar=0

协方差矩阵的一些计算例子，
X=[123311]=[c1c2c3]X=\begin{bmatrix} 1 & 2 &3 \\ 3 & 1 & 1 \end{bmatrix}=\begin{bmatrix} c1 & c2 &c3 \end{bmatrix} X=[132131]=[c1c2c3]

1.求每个维度的平均值
c‾=[21.52]=[c1‾c2‾c3‾]\overline c = \begin{bmatrix} 2&{1.5}&2 \end{bmatrix} = \begin{bmatrix} {\overline {c1} }&{\overline {c2} }&{\overline {c3} } \end{bmatrix} c=[21.52]=[c1c2c3]
2.将X的每一列减去平均值
X=[−10.511−0.5−1]X = \begin{bmatrix} { - 1}&{0.5}&1\\ 1&{ - 0.5}&{ - 1} \end{bmatrix} X=[−110.5−0.51−1]
3.计算协方差矩阵
Cov=1m−1XTX=12−1[2−1−2−10.51−212]Cov = \frac{1}{{m - 1}}{X^T}X = \frac{1}{{2 - 1}} \begin{bmatrix} 2&{ - 1}&{ - 2}\\ { - 1}&{0.5}&1\\ { - 2}&1&2 \end{bmatrix} Cov=m−11XTX=2−11⎣⎡2−1−2−10.51−212⎦⎤
注意这就和推导的协方差矩阵不同，但是不影响关于最终主方向的性质。

参考协方差矩阵计算方法

(一) 三维点云课程---PCA介绍相关推荐

三维点云：PCA（下）open3d
三维点云之PCA应用下在三维点云上的应用主成分方向降维求解点云的法向量拟合平面 AABB框点云匹配在三维点云上的应用 PCA 是有损的数据压缩方式,它常用于对高维数据进行降维,也就是把高 ...
(四) 三维点云课程---PointNet-Pytorch运行
三维点云课程-PointNet-Pytorch运行三维点云课程---PointNet-Pytorch运行三维点云课程---PointNet-Pytorch运行 1.分类---Classificat ...
(七) 三维点云课程---ICP应用
(七)三维点云课程-ICP应用三维点云课程---ICP应用 (七)三维点云课程---ICP应用 1.代码解读 1.1滤波 1.2 计算特征向量 1.3 特征向量的配准 1.3.1 get_poten ...
深蓝学院的三维点云课程：第一章
一.前言为什么现在点云应用这么广泛,就是因为他有深度信息. 像人脸识别用来解锁手机,比如Iphnoe手机在前边有一个深度摄像头,所以它产生的点云真的是一个三维点云:然后很多手机他可能就没有深度摄像头 ...
三维点云课程（七）——特征点描述
1. 什么是特征点 1.1 图像特征点 ORB slam 1.2 点云特征点点云配准 :ICP要求有足够好的初始平移旋转矩阵,且有一定的重合率 2. 怎么提取特征点 2.1 图像提取特征点 2.1. ...
三维点云课程第一章：应用
核PCA的核心思想就是我们可以把高维空间的运算转化成低维空间的一个核函数,这个思想叫做Kernel track.也就是说我们在学习SVM支持向量机的时候,也会经常的用到这个核的方法. 应用1:如何去 ...
三维点云课程（六）——三维目标检测
目录 1. 图像目标检测 1.1 评价检测好坏 1.2 物体检测的方法 1.2.1 Two-Stage 1.2.2 One-Stage 2. 点云目标检测 2.1 VoxelNet 2.2 P ...
课程笔记-三维点云处理01 ——Introduction and Basic Algorithms
课程笔记-三维点云处理01 --Introduction and Basic Algorithms 本系列笔记是对深蓝学院所开设的课程:<三维点云处理>的笔记课程每周更新,我也会努力将每 ...
基于三维点云数据的主成分分析方法（PCA）的python实现
主成分分析(PCA)获取三维点云的坐标轴方向和点云法向量 # 实现PCA分析和法向量计算,并加载数据集中的文件进行验证import open3d as o3d # import os import n ...

(一) 三维点云课程---PCA介绍

三维点云课程—PCA介绍