如何计算协方差矩阵的特征根R语言 协方差矩阵特征值

协方差矩阵的特征值

由于求协方差矩阵的特征值具有非常重要的地位，为此，我们专讨论它。

（0）随机变量的方差、两随机变量的协方差

对于随机变量$x$的观察值（样本）集$\{x_j\}_{j=1}^m$，有均值和方差（许多同学总是把方差（估值）中的系数错误地记成了$\frac{1}{m}$）：
$\begin{align} \text{（样本）均值：}\mu & =\frac{1}{m}\sum_{j=1}^mx_j \tag{D1} \\ \text{（样本）方差：}\sigma ^2 & =\frac{1}{m-1}\sum_{j=1}^m(x_j-\mu)^2 \tag{D2} \end{align}$
注：若定义随机变量$x$的数学期望$\mathbb{E}\, (x)$和方差$\mathrm{Var}(x)$分别为$\mu ,\sigma ^2$，则这里的样本均值和样本方差分别为其对应的估值，这时为区别，这里就应改用$\hat{\mu},\hat{\sigma}$，后续情况类似。

若数据集已“中心化” ，则
$\begin{align} \mu =0,\quad \sigma ^2 =\frac{1}{m-1}\sum_{j=1}^mx_j^2 \tag{D3} \end{align}$

设随机变量$x$与$y$的样本集分别为$\{x_j\}_{j=1}^m$、$\{y_j\}_{j=1}^m$，则可定义
$\begin{align} \text{协方差：}\sigma_{xy} & =\frac{1}{m-1}\sum_{j=1}^m(x_j-\mu_x)(y_j-\mu_y) \tag{D4} \\ \text{相关系数：}\rho _{xy} & =\frac{\sigma_{xy} }{\sigma_{x} \sigma_{y} } \tag{D5} \end{align}$
其中，$\mu_x$对应于式(D1)，$\sigma_{x}^2$对应于式(D2)，其余类推。

（1）随机向量$\boldsymbol{x}$样本集$\mathbf{X}$的协方差矩阵

现在把随机变量扩展成随机向量$\boldsymbol{x}=(x_1;x_2;\cdots;x_d)$，样本集为$\{\boldsymbol{x}_j\}_{j=1}^m$，每个分量$x_i$均为随机变量，将样本集作为矩阵
$\begin{align*} \mathbf{X} & =(\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_m)\notag \\ & = \left(\begin{matrix} x_{11} & x_{12} & \cdots & x_{1m} \\ x_{21} & x_{22} & \cdots & x_{2m} \\ \vdots & \vdots & \vdots & \vdots \\ x_{d1} & x_{d2} & \cdots & x_{dm} \\ \end{matrix}\right) \end{align*}$
其中，矩阵每列是一个向量样本（共$m$个样本），每行是一个分量的样本集，分量$x_i$的样本集为$\{x_{ij}\}_{j=1}^m$，
依上述定义，每个分量$x_i$有上述样本均值和样本方差，两分量$x_{i_1},x_{i_2}$间有协方差和相关系数。

现在假定对向量样本集$\mathbf{X}$的各分量作中心化（即对行作中心化）得到向量样本集为$\mathbf{X$，
则由式(D1)、式(D2)有
$\begin{align} \mathbf{X$
其中，右侧除了系数$(m-1)$外，矩阵的主对角线为分量的方差，其他元素为两分量的协方差，因此，定义$\mathbf{X$为样本集$\mathbf{X}$的协方差矩阵。

（2）求协方差矩阵不必中心化

设$\mathbf{X}$由$\mathbf{X}$中心化所得，即
$\begin{align} \boldsymbol{x}$
其中，$\boldsymbol{1}$为元素全为1的（列）向量，定义向量$\boldsymbol{\mu}\mathop{=} \limits^{\mathrm{def}} \frac{1}{m}\mathbf{X} \boldsymbol{1}$，定义各列全为向量$\boldsymbol{\mu}$的矩阵$\mathbf{A}$
$\begin{align} \mathbf{A} & \mathop{=} \limits^{\mathrm{def}}(\boldsymbol{\mu},\boldsymbol{\mu},\cdots,\boldsymbol{\mu})\notag \\ & =\boldsymbol{\mu}(1,1,\cdots,1)\notag \\ & =\boldsymbol{\mu}\boldsymbol{1}^\mathrm{T}\notag \\ & =\frac{1}{m}\mathbf{X} \boldsymbol{1}\boldsymbol{1}^\mathrm{T} \tag{D9} \end{align}$
$\begin{align*} \mathbf{X}$
其中，$\mathbf{H}=\mathbf{I}-\frac{1}{m}\boldsymbol{1}\boldsymbol{1}^\mathrm{T}$。

由此，协方差矩阵
$\begin{align} \mathbf{X}$
这说明不需要显式的中心化过程，即可通过式(D10)右边求出样本集$\mathbf{X}$的协方差矩阵$\mathbf{X}$。

（3）通过$\mathbf{Y}$的奇异值分解得到$\mathbf{Y}\mathbf{Y}^{\mathrm{T}}$的特征分解

我们讨论如何通过矩阵$\mathbf{Y}$的奇异值分解$\boldsymbol{\Sigma }$得到$\mathbf{Y}\mathbf{Y}^{\mathrm{T}}$的特征分解$\boldsymbol{\Lambda }$。

由【西瓜书附录式(A.33)】，$\mathbf{Y}$的奇异值分解表述为如下：
$\begin{align} \begin{cases} \mathbf{Y}=\mathbf{U}\boldsymbol{\Sigma }\mathbf{V}^{\mathrm{T}} \\ \mathbf{Y}\in \mathbb{R} ^{d \times m},\mathbf{U}\in \mathbb{R} ^{d \times d},\boldsymbol{\Sigma }\in \mathbb{R} ^{d \times m},\mathbf{V}\in \mathbb{R} ^{m \times m} \\ \mathbf{U }^{\mathrm{T}}\mathbf{U}=\mathbf{I}_{d \times d},\ \mathbf{V}^{\mathrm{T}}\mathbf{V}=\mathbf{I}_{m \times m} \\ \boldsymbol{\Sigma }= \begin{cases} (\mathrm{diag}({\sigma }_1,{\sigma }_2,\cdots,{\sigma }_d)\ ,\ \mathbf{0 }),&\quad (d\leqslant m) \\ (\mathrm{diag}({\sigma }_1,{\sigma }_2,\cdots,{\sigma }_m)\ ;\ \mathbf{0 }),&\quad (d> m) \\ \end{cases} \\ {\sigma }_1\geqslant {\sigma }_2\geqslant \cdots\geqslant 0 \\ \end{cases} \tag{D11} \end{align}$
其中，$\mathrm{diag}({\sigma }_1,{\sigma }_2,\cdots,{\sigma }_d)$为对角矩阵。

则
$\begin{align} \mathbf{Y}\mathbf{Y}^{\mathrm{T}} & =(\mathbf{U}\boldsymbol{\Sigma }\mathbf{V}^{\mathrm{T}})(\mathbf{U}\boldsymbol{\Sigma }\mathbf{V}^{\mathrm{T}})^{\mathrm{T}}\notag \\ & =\mathbf{U}\boldsymbol{\Sigma }\mathbf{V}^{\mathrm{T}}\mathbf{V}\boldsymbol{\Sigma }\mathbf{U}^{\mathrm{T}}\notag \\ & =\mathbf{U}\boldsymbol{\Sigma }\boldsymbol{\Sigma }\mathbf{U}^{\mathrm{T}}\notag \\ & =\mathbf{U}\boldsymbol{\Lambda }\mathbf{U}^{\mathrm{T}} \tag{D12} \end{align}$
式(D12)即为$\mathbf{Y}\mathbf{Y}^{\mathrm{T}}$的特征分解，其中，对角阵
$\begin{align} \boldsymbol{\Lambda }=\boldsymbol{\Sigma }\boldsymbol{\Sigma }=\boldsymbol{\Sigma }^2,\quad \text{即： }(\boldsymbol{\Lambda })_{ii}={\lambda }_i={\sigma}_i^2 \tag{D13} \end{align}$
因为，由式(D12)有
$\begin{align} \mathbf{Y}\mathbf{Y}^{\mathrm{T}}\mathbf{U} & =\mathbf{U}\boldsymbol{\Lambda }\mathbf{U}^{\mathrm{T}}\mathbf{U}\notag \\ & =\mathbf{U}\boldsymbol{\Lambda } \tag{D14} \end{align}$
即
$\begin{align} \mathbf{Y}\mathbf{Y}^{\mathrm{T}}\boldsymbol{u}_i & =\boldsymbol{u}_i{\lambda }_i \tag{D15} \end{align}$
由式(D15)及特征向量的定义可知，$\boldsymbol{u}_i$为特征向量，${\lambda }_i$为对应的特征值，这就通过$\mathbf{Y}$的奇异值分解$\boldsymbol{\Sigma }$实现了$\mathbf{Y}\mathbf{Y}^{\mathrm{T}}$的特征分解$\boldsymbol{\Lambda }=\boldsymbol{\Sigma }^2$，且特征值为奇异值的平方${\lambda }_i={\sigma}_i^2$，次序为从大到小排序。

实际问题（例如，10.2 低维嵌入（立交桥就是嵌入三维空间中的二维）特征分解进行降维的技巧的式(10.17)中）中，往往最终不是求协方差矩阵$\mathbf{X}$（求它只是过程中的过渡），而是要求出协方差矩阵的特征值（特征矩阵），由（1）知，只需求出$\mathbf{X}\mathbf{H}\mathbf{H}^\mathrm{T}\mathbf{X}^\mathrm{T}$的特征分解。
由（2）知，只需通过$\mathbf{X}\mathbf{H}$（视作$\mathbf{Y}$）的奇异值分解来实现$\mathbf{X}\mathbf{H}\mathbf{H}^\mathrm{T}\mathbf{X}^\mathrm{T}$的特征分解。 即协方差矩阵的特征分解归结为$\mathbf{X}\mathbf{H}$的奇异值分解，而奇异值分解又有现成的计算工具。

综上，求样本集$\mathbf{X}$（不必中心化）协方差矩阵的特征矩阵的方法：对$\mathbf{X}\mathbf{H}$进行奇异值分解，得$\boldsymbol{\Sigma }$，再由${\lambda }_i={\sigma}_i^2$组成特征矩阵$\boldsymbol{\Lambda }$。