主成分分析原理的数学证明

  问题为:假设有$m$个$n$维实数域的样本$x_1,x_2,...,x_m$,其中$x_i=(x_{i1},x_{i2},...,x_{in})$.PCA希望将样本的维度从$n$维减少到$k(k<n)$维,当然承受一定程度(尽可能小的程度)的信息损失.首先给出主成分分析的结论,再进行证明.

1. 将$m$个样本沿着维度的方向0均值化后按行排成一个矩阵$X_{m,n}$并计算得到协方差矩阵$A=X^{T}X$;
2. 计算协方差矩阵$A$的$n$个特征值${\lambda }_1,{\lambda }_2,...,{\lambda }_n$(可能有重复值)和对应的特征向量$u_1,u_2,...,u_n$;
3. 将$n$个特征值按照从大到小的顺序进行排序,相应的将$n$个特征向量也排序.
4. 选前$k$个特征向量,按列排列得到一个变换矩阵$T_{n,k}$,每个样本经过变换矩阵后得到降维后的样本(由$n$维变为$k$维).

  后面的证明主要就集中在对这个样本矩阵的协方差矩阵$X^{T}X$来进行.这个矩阵有一些很好的性质提前写出来,下面的证明过程会用到.性质2的证明方法是构造法,性质1和3的证明都比较简单不再赘述,可以参考任意一本线性代数的教材.

性质1: $X^{T}X$是一个$n$维实对称矩阵,其所有特征值非负.

性质2: $X^{T}X$一定有$n$个线性无关的特征向量,即其一定能够对角化.

性质3: $X^{T}X$的属于不同特征值的特征向量两两正交,又其属于同一特征值的特征向量可以通过施密特正交化,再结合性质2可知,存在正交矩阵$Q$可以将$X^{T}X$对角化,即$Q^{-1}{X}^{T}XQ=diag(X^{T}X)$,其中$diag(X^{T}X)$为对角矩阵,对角线上的元素分别为与特征向量组成的矩阵$Q$对应的特征值.

  下面给出PCA步骤的证明.

  主成分分析的目的在于放弃原来的$n$维样本特征空间,另外找一个$k$维的样本特征空间,并将原来的每个$n$维的样本转换到这个维度更小的空间中来.转换的原则是:

1. 新的$k$维空间是一个正交空间,即$k$个基是两两正交的.
2. 在新的$k$维空间中,样本的新坐标具有最大的方差.这里稍微解释一下,根据方差的定义,最大的方差意味着样本的整体分布更加分散,而不是都集中在某个很小的区域里.样本在空间里分布的越分散,越容易体现他们的差异性和各自的特点.

  具体来说就是:寻找$k$个两两正交的单位方向向量$u$,使得样本在这些方向上的投影方差是所有方向向量里分别是第一大的,第二大的,…,第k大的.

  怎么寻找这些方向向量$u$呢?对于某一个任意的方向向量$u$,由于之前对样本做过零均值的归一化操作,又偏置值为0的线性变换（相当于对向量进行绕原点的旋转）不改变样本的均值,故样本在该方向上投影后的均值仍然为0.将样本在变换前的空间内记做$x_i$,在变换后的空间内记做$x_i^{\prime }$.根据方差的定义,所有样本在变换后的空间内的平均方差为:

  $Var(X^{\prime })=\frac{1}{n}{\sum }_{i=1}^n(x_i^{\prime }-E(x^{\prime }){)}^2=\frac{1}{n}{\sum }_{i=1}^n||x_i^{\prime }|{|}_2^2=\frac{1}{n}{\sum }_{i=1}^n(x_i,u{)}^2=\frac{1}{n}{u}^T{X}^{T}Xu,s.t.||u||=1$

  上式等价于,寻找方向向量$u$使得$f(u)=u^T{X}^{T}Xu,s.t.||u||=1$取得最大值.假设样本的各个维度的特征相互独立,我们将自变量$u=(u_1,u_2,...,u_n{)}^T$看做一个含有$n$个独立变量的多元变量,那么这个问题就显然是一个多元函数的条件极值问题,可以通过求解拉格朗日函数的极大值来解决:

  $L(u)=u^T{X}^{T}Xu+\lambda (1-u^{T}u)$

根据多元微分学的结论可知,上述拉格朗日函数在自变量$u_0=(u_{10},u_{20},...,u_{n0}{)}^T$取得极大值的充要条件是在该点一阶偏导数全部等于零且Hessian矩阵半负定.即:

  $\frac{\partial L}{\partial u_1}{|}_{u=u_0}=0,\frac{\partial L}{\partial u_2}{|}_{u=u_0}=0,...,\frac{\partial L}{\partial u_n}{|}_{u=u_0}=0$

  $Hessian(L)={\left|\begin{array}{cccc}\frac{{\partial }^{2}L}{\partial {u_1}^2}& \frac{{\partial }^{2}L}{\partial u_1\partial u_2}& ...& \frac{{\partial }^{2}L}{\partial u_1\partial u_n}\\ \frac{{\partial }^{2}L}{\partial u_2\partial u_1}& \frac{{\partial }^{2}L}{\partial {u_2}^2}& ...& \frac{{\partial }^{2}L}{\partial u_2\partial u_n}\\ ...& ...& ...& ...\\ \frac{{\partial }^{2}L}{\partial u_n\partial u_1}& \frac{{\partial }^{2}L}{\partial u_n\partial u_2}& ...& \frac{{\partial }^{2}L}{\partial {u_n}^2}\end{array}\right|}_{u=u_0}$负定

分别求$L(u)$的一阶导数和Hessian矩阵,结果为:

  $\frac{\partial L}{\partial u}=2(X^{T}Xu-\lambda u)$

  $Hessian(L)=2(X^{T}X-\lambda I)$

  上述两个式子看起来很简单,不过要求这个结果没点向量函数的导数知识还是比较难的,这里有两个路子可供大家选择.第一种办法是去看一下标量函数对向量变量的求导(算一阶偏导)和向量函数对向量变量的求导(算Hessian矩阵).第二种办法比较直观但是也比较麻烦,可以把$X^{T}X=\left|\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1n}\\ a_{21}& a_{22}& ...& a_{2n}\\ ...& ...& ...& ...\\ a_{n1}& a_{n2}& ...& a_{nn}\end{array}\right|$代入展开,再利用微积分的求偏导的方法去求就可以了.

  一阶导数等于0得到:$X^{T}Xu=\lambda u$,即$\lambda$取$X^{T}X$的特征值,$u$取其对应的特征向量$u_0$即可.这里我们假设先让$\lambda$取所有特征值里最大的那个,来看看此时得到的$u_0$是否能使得拉格朗日函数取得最大值.

根据$X^{T}X$的性质3, 此时Hessian矩阵:

$$\begin{gathered} Hessian(L)=X^{T}X-\lambda I=Q^{-1}diag(X^{T}X)Q-{\lambda }_{max}{Q}^{-1}IQ \\ =Q^{-1}(diag(X^{T}X)-{\lambda }_{max}I)Q=Q^{-1}\left|\begin{array}{cccc}{\lambda }_1-{\lambda }_{max}& 0& ...& 0\\ 0& {\lambda }_2-{\lambda }_{max}& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& {\lambda }_n-{\lambda }_{max}\end{array}\right|Q \end{gathered}$$

则Hessian矩阵一定相合于标准型$\left|\begin{array}{cccc}{\lambda }_1-{\lambda }_{max}& 0& ...& 0\\ 0& {\lambda }_2-{\lambda }_{max}& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& {\lambda }_n-{\lambda }_{max}\end{array}\right|$,由于${\lambda }_{max}$是所有特征值中最大的,则Hessian矩阵一定是半负定的,则此时拉格朗日函数一定取得极大值.至此就得到了一个合适的方向向量,即最大特征值对应的特征向量,也就是所谓的第一主成分.对于剩下的主成分也是类似来求,只不过按照特征值从大到小的顺序遍历就行了.

  这里还有一个小问题交代一下,通过上述分析我们知道,对于最大的特征值可以得到Hessian矩阵相合于一个半负定的矩阵.但是当取第二大的特征值时,相合矩阵的第一行第一列值就为正了,此时标准型是不定的.那怎么解释这个问题呢?我的想法是由于属于不同特征值的特征向量是相互正交的(属于统一特征值的特征向量也可以正交化),所以取第二大的特征值时,虽然标准型的第一个数字不为0或负数,但是此时特征向量向第一个特征值对应的特征向量的投影分量等于0,所以二次型的计算结果中,系数虽然不为0,但是相当于向量的等效模为0,所以可以忽略第一个特征值的部分.仅仅看第二大特征值之后的部分,该二次型还是半负定的.其余特征值情况同理.暂时没想到这个方法的不成立之处,也没想到其他更好的证明方法,希望能得到更多的讨论.