主成分分析|机器学习推导系列(五)

一、简介

1. 为什么需要降维

数据的维度过高容易造成维数灾难(Curse of Dimensionality)。.

维数灾难:通常是指在涉及到向量的计算的问题中,随着维数的增加,计算量呈指数倍增长的一种现象。

这里可以举两个几何的例子来看一下维数过高的影响:

上图表示一个多维空间(以二维为例),则其中图形的体积有如下关系:

V_{超立方体}=1\\ V_{超球体}=K\cdot 0.5^{D}\\ \lim_{D\rightarrow \infty }V_{超球体}=0

上式也就表明当数据的维度过高时,数据主要存在于空间的边边角角的地方,这也就造成了数据的稀疏性。

上图也表示一个多维空间(以二维为例),则其中图形的体积有如下关系:

V_{外}=K\cdot 1^{D}=K\\ V_{环形带}=V_{外}-V_{内}=K-K\cdot (1-\varepsilon )^{D}\\ \frac{V_{环形带}}{V_{外}}=\frac{K-K\cdot (1-\varepsilon )^{D}}{K}=1-(1-\varepsilon )^{D}\\ \lim_{D\rightarrow \infty }\frac{V_{环形带}}{V_{外}}=\lim_{D\rightarrow \infty }1-(1-\varepsilon )^{D}=1

可以看到当数据的维度过高时,数据主要存在于球壳上,类似于人的大脑皮层。

2. 降维的方法

降维可以作为一种防止过拟合的方式,其具体的方法包含下列几种:

降维\left\{\begin{matrix} 直接降维(特征选择)\\ 线性降维(PCA、MDS)\\ 非线性降维(流形) \end{matrix}\right.

特征选择是一种直接剔除主观认为不重要的特征的过程。

本文接下来的部分主要介绍主成分分析(PCA)。

二、样本均值与样本方差

1. 概述

假设有以下数据:

x_{i}\in \mathbb{R}^{p},i=1,2,\cdots ,N\\ X=(x_{1},x_{1},\cdots ,x_{N})^{T}=\begin{pmatrix} x_{1}^{T}\\ x_{2}^{T}\\ \vdots \\ x_{N}^{T} \end{pmatrix}=\begin{pmatrix} x_{11} & x_{12} & \cdots &x_{1p} \\ x_{21} & x_{22}& \cdots &x_{2p} \\ \vdots & \vdots & \ddots &\vdots \\ x_{N1}& x_{N2} & \cdots & x_{Np} \end{pmatrix}_{N \times p}

2. 样本均值与样本方差

以下定义了数据的样本均值与样本方差:

Sample\: Mean:\bar{x}_{p\times 1}=\frac{1}{N}\sum_{i=1}^{N}x_{i}\\ Sample\: Covariance:S_{p\times p}=\frac{1}{N}\sum_{i=1}^{N}(x_{i}-\bar{x})(x_{i}-\bar{x})^{T}

接下来需要对样本均值与样本方差进行一些变换来获得其另一种表示形式:

规定向量1_{N}=\begin{pmatrix} 1\\ 1\\ \vdots \\ 1 \end{pmatrix}_{N\times 1}\\ \bar{x} =\frac{1}{N}\sum_{i=1}^{N}x_{i}=\frac{1}{N}\underset{X^{T}}{\underbrace{\begin{pmatrix} x_{1} & x_{2} & \cdots & x_{N} \end{pmatrix}}}\begin{pmatrix} 1\\ 1\\ \vdots \\ 1 \end{pmatrix}=\frac{1}{N}X^{T}1_{N}\\ S=\frac{1}{N}\sum_{i=1}^{N}(x_{i}-\bar{x})(x_{i}-\bar{x})^{T}\\ =\frac{1}{N}\begin{pmatrix} x_{1}-\bar{x} & x_{2}-\bar{x} & \cdots & x_{N}-\bar{x} \end{pmatrix}\begin{pmatrix} (x_{1}-\bar{x})^{T}\\ (x_{2}-\bar{x})^{T}\\ \vdots \\ (x_{N}-\bar{x})^{T} \end{pmatrix}\\ 上式中\begin{pmatrix} x_{1}-\bar{x} & x_{2}-\bar{x} & \cdots & x_{N}-\bar{x} \end{pmatrix}\\ =\begin{pmatrix} x_{1} & x_{2} & \cdots & x_{N} \end{pmatrix}-\begin{pmatrix} \bar{x} & \bar{x} & \cdots & \bar{x} \end{pmatrix}\\ =X^{T}-\bar{x}\begin{pmatrix} 1 & 1 & \cdots & 1 \end{pmatrix}\\ =X^{T}-\bar{x}1_{N}^{T}\\ =X^{T}-\frac{1}{N}X^{T}1_{N}1_{N}^{T}\\ =X^{T}(I_{N}-\frac{1}{N}1_{N}1_{N}^{T})\\ 则S=\frac{1}{N}X^{T}\underset{H}{\underbrace{(I_{N}-\frac{1}{N}1_{N}1_{N}^{T})}}(I_{N}-\frac{1}{N}1_{N}1_{N}^{T})^{T}X\\ (H称为中心矩阵,centering\; matrix)\\ =\frac{1}{N}X^{T}HH^{T}X

中心矩阵H具备以下性质:

①\; H^{T}=H\\ H^{T}=(I_{N}-\frac{1}{N}1_{N}1_{N}^{T})^{T}=I_{N}-\frac{1}{N}1_{N}1_{N}^{T}=H\\ ②\; H^{n}=H\\ H^{2}=H\cdot H=(I_{N}-\frac{1}{N}1_{N}1_{N}^{T})(I_{N}-\frac{1}{N}1_{N}1_{N}^{T})\\ =I_{N}-\frac{2}{N}1_{N}1_{N}^{T}+\frac{1}{N^{2}}1_{N}1_{N}^{T}1_{N}1_{N}^{T}\\ =I_{N}-\frac{2}{N}\begin{pmatrix} 1\\ 1\\ \vdots \\ 1 \end{pmatrix}\begin{pmatrix} 1 & 1 & \cdots & 1 \end{pmatrix}+\frac{1}{N^{2}}\begin{pmatrix} 1\\ 1\\ \vdots \\ 1 \end{pmatrix}\begin{pmatrix} 1 & 1 & \cdots & 1 \end{pmatrix}\begin{pmatrix} 1\\ 1\\ \vdots \\ 1 \end{pmatrix}\begin{pmatrix} 1 & 1 & \cdots & 1 \end{pmatrix}\\ =I_{N}-\frac{2}{N}\begin{bmatrix} 1 & 1 & \cdots & 1 \\ 1 & 1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 1 \end{bmatrix}_{N\times N}+\frac{1}{N^{2}}\begin{bmatrix} 1 & 1 & \cdots & 1 \\ 1 & 1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 1 \end{bmatrix}_{N\times N}\begin{bmatrix} 1 & 1 & \cdots & 1 \\ 1 & 1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 1 \end{bmatrix}_{N\times N}\\ =I_{N}-\frac{2}{N}\begin{bmatrix} 1 & 1 & \cdots & 1 \\ 1 & 1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 1 \end{bmatrix}_{N\times N}+\frac{1}{N^{2}}\begin{bmatrix} N & N & \cdots & N \\ N & N & \cdots & N \\ \vdots & \vdots & \ddots & \vdots \\ N & N & \cdots & N \end{bmatrix}_{N\times N}\\ =I_{N}-\frac{2}{N}\begin{bmatrix} 1 & 1 & \cdots & 1 \\ 1 & 1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 1 \end{bmatrix}_{N\times N}+\frac{1}{N}\begin{bmatrix} 1 & 1 & \cdots & 1 \\ 1 & 1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 1 \end{bmatrix}_{N\times N}\\ =I_{N}-\frac{1}{N}1_{N}1_{N}^{T}\\ =H\\ 则H^{n}=H

因此最终可以得到

\bar{x} =\frac{1}{N}X^{T}1_{N}\\ S=\frac{1}{N}X^{T}HX

三、主成分分析的思想

总结起来就是:

一个中心:PCA是对原始特征空间的重构,将原来的线性相关的向量转换成线性无关的向量;
两个基本点:最大投影方差和最小重构距离,这是本质相同的两种方法,在接下来的部分将具体介绍。

PCA首先要将数据中心化(即减去均值)然后投影到一个新的方向上,这个新的方向即为重构的特征空间的坐标轴,同时也要保证投影以后得到的数据的方差最大,即最大投影方差,这样也保证了数据的重构距离最小。

四、最大投影方差

假设投影方向为u,由于我们只关注投影的方向,因此将u的模设置为1,即u^{T}u=1,则中心化后的数据在u方向上的投影为(x_{i}-\bar{x} )^{T}u,是一个标量。按照最大投影方差的思想,我们定义损失函数如下:

J(u)=\frac{1}{N}\sum_{i=1}^{N}((x_{i}- \bar{x})^{T}u)^{2}\\ =\sum_{i=1}^{N}\frac{1}{N}u^{T}(x_{i}- \bar{x})(x_{i}- \bar{x} )^{T}u\\ =u^{T}\underset{S}{\underbrace{[\frac{1}{N}\sum_{i=1}^{N}(x_{i}- \bar{x})(x_{i}- \bar{x})^{T}]}}u\\ =u^{T}Su

因此该问题就转换为以下最优化问题:

\left\{\begin{matrix} \hat{u}=\underset{u}{argmax}\;u^{T}Su\\ s.t.\; u^{T}u=1 \end{matrix}\right.

然后使用拉格朗日乘子法进行求解:

L(u,\lambda )=u^{T}Su+\lambda (1-u^{T}u)\\ \frac{\partial L}{\partial u}=2Su-2\lambda u=0\\ S\underset{特征向量}{u}=\underset{特征值}{\lambda} u

最后解得符合条件的向量是协方差矩阵S的特征向量。如果想要降到q维(q<p),则只需要将对应特征值最大的前q个特征向量取出来作为投影方向然后获得数据在这些方向上的投影即为重构的坐标,即:

\begin{pmatrix} x_{1}^{T}\\ x_{2}^{T}\\ \vdots \\ x_{N}^{T} \end{pmatrix}_{N\times p}\begin{pmatrix} u_{1} & u_{2} & \cdots & u_{q} \end{pmatrix}_{p\times q}=\begin{bmatrix} x_{1}^{T}u_{1}& x_{1}^{T}u_{2}& \cdots & x_{1}^{T}u_{q}\\ x_{2}^{T}u_{1}& x_{2}^{T}u_{2}& \cdots & x_{2}^{T}u_{q}\\ \vdots & \vdots & \ddots & \vdots \\ x_{N}^{T}u_{1}& x_{N}^{T}u_{2}& \cdots & x_{N}^{T}u_{q} \end{bmatrix}_{N\times q}

特征向量表示投影变换的方向,特征值表示投影变换的强度。通过降维,我们希望减少冗余信息,提高识别的精度,或者希望通过降维算法来寻找数据内部的本质结构特征。找最大的特征值是因为 ,在降维之后要最大化保留数据的内在信息,并期望在所投影的维度上的离散最大。

五、最小重构距离

最小重构距离是另一种求解的方法,其本质上和最大投影方差是相同的。

我们知道有p个投影方向符合条件,因此原来的数据可以表示为以下形式,降维的数据也就是舍弃掉第q+1到第p这几个方向上的信息。

原来的中心化了的数据x_{i}-\bar{x}=\sum_{k=1}^{p}((x_{i}-\bar{x})^{T}u_{k})u_{k}\\降维的数据\hat{x}_{i}=\sum_{k=1}^{q}((x_{i}-\bar{x})^{T}u_{k})u_{k}\\(u_{1}到u_{p}分别对应从大到小的特征值)

因此重构距离也就是指x_{i}-\hat{x}_{i},本着最小化重构距离的思想我们可以设置新的损失函数如下:

J=\frac{1}{N}\sum_{i=1}^{N}\left \| (x_{i}- \bar{x})-\hat{x}_{i}\right \|^{2}\\ =\frac{1}{N}\sum_{i=1}^{N}\left \| \sum_{k=q+1}^{p}((x_{i}-\bar{x})^{T}u_{k})u_{k}\right \|^{2}\\ =\frac{1}{N}\sum_{i=1}^{N}\sum_{k=q+1}^{p}((x_{i}-\bar{x})^{T}u_{k})^{2}\\ =\sum_{k=q+1}^{p}\underset{u_{k}^{T}Su_{k}}{\underbrace{\frac{1}{N}\sum_{i=1}^{N}((x_{i}-\bar{x})^{T}u_{k})^{2}}}\\ =\sum_{k=q+1}^{p}u_{k}^{T}Su_{k}\\ s.t.\; u_{k}^{T}u_{k}=1

然后就可以转化为以下最优化问题:

\left\{\begin{matrix} \hat{u}=argmin\sum_{k=q+1}^{p}u_{k}^{T}Su_{k}\\ s.t.\; u_{k}^{T}u_{k}=1 \end{matrix}\right.

显然这里的每个u_{k}是可以单独求解的,最终也可以解得u_{k}是协方差矩阵S的特征向量,只不过这里的u_{k}是对应特征值较小的几个特征向量。

六、SVD角度看PCA和PCoA

协方差矩阵S的特征分解:

S=GKG^{T},其中G^{T}G=I,K=\begin{bmatrix} k_{1} & & & \\ & k_{2} & & \\ & & \ddots & \\ & & & k_{p} \end{bmatrix},k_{1}\geq k_{2}\geq \cdots \geq k_{p}

X中心化的结果HX做奇异值分解:

HX=U\Sigma V^{T},其中\left\{\begin{matrix} U_{N\times N}是正交矩阵\\ V_{p\times p}是正交矩阵\\ \Sigma _{N\times p}是对角矩阵 \end{matrix}\right.

接下里可以做以下变换:

S_{p\times p}=X^{T}HX=X^{T}H^{T}HX=V\Sigma^{T} U^{T}U\Sigma V^{T}=V\Sigma^{T}\Sigma V^{T}\\ (V\Sigma^{T}\Sigma V^{T}是S的特征值分解,\Sigma^{T}\Sigma 即为K。)

接下来我们构造矩阵T_{N\times N}

T_{N\times N}=HXX^{T}H^{T}=U\Sigma V^{T}V\Sigma^{T} U^{T}=U\Sigma \Sigma^{T} U^{T}\\ (U\Sigma \Sigma^{T} U^{T}是T的特征值分解,\Sigma \Sigma^{T} 为特征值矩阵。)

对比S_{p\times p}T_{N\times N},我们可以发现:
①将S进行特征分解然后得到投影的方向,也就是主成分,然后矩阵HXV即为重构坐标系的坐标矩阵;
②将T进行特征分解可以直接获得坐标矩阵U\Sigma
(注意应保证ST特征分解得到的特征向量是单位向量。)

关于为什么将T进行特征分解可以直接获得坐标矩阵,现做以下解释:

坐标矩阵HXV=U\Sigma V^{T}V=U\Sigma \\ 也就是说U\Sigma 即为坐标矩阵\\ 接着T{\color{Red} {U\Sigma}} =U\Sigma \Sigma^{T} U^{T}U\Sigma ={\color{Red} {U\Sigma}} (\Sigma^{T} \Sigma )\\ 也就是说U\Sigma是T的特征向量组成的矩阵

使用T进行特征分解的方法叫做主坐标分析(Principal Co-ordinates Analysis,PCoA)。

这两种⽅法都可以得到主成分,但是由于⽅差矩阵是p\times p的,⽽TN\times N的,所以对样本量较少的时候可以采⽤ PCoA的⽅法。

七、概率PCA(p-PCA)

1. 概述

假设有以下数据:

x\in \mathbb{R}^{p},z\in \mathbb{R}^{q},q< p

其中x是原始数据,z是降维后的数据,可以将z看做隐变量(latent variable),x看做观测变量(observed variable),则p-PCA就可以看做生成模型。

xz满足以下关系:

\left\{\begin{matrix} z\sim N(0_{q\times 1},I_{q\times q})\\ x=Wz+\mu +\varepsilon \\ \varepsilon \sim N(0_{p\times 1},\sigma ^{2}I_{p\times p}) \end{matrix}\right.

这是一个线性高斯模型,其中\varepsilon是噪声,\varepsilonz是独立的。求解这个模型要经过两个阶段:
①inference:求P(z|x)
②learning:使用EM算法求解参数W、\mu,\sigma ^{2}

x的生成过程如下:

生成过程

上图中数据空间为⼆维,潜在空间为⼀维。⼀个观测数据点x的⽣成⽅式为:⾸先从潜在变量的先验分布p(z)中抽取⼀个潜在变量的值\hat{z},然后从⼀个各向同性的⾼斯分布(⽤红⾊圆圈表示)中抽取⼀个x的值,这个各向同性的⾼斯分布的均值为W\hat{z}+\mu,协⽅差为σ^{2}I。绿⾊椭圆画出了边缘概率分布p(x)的密度轮廓线。

2. 推断(inference)

求解P(z|x)的过程如下:

P(z)\rightarrow P(x|z)\rightarrow P(x)\rightarrow P(z|x)

  • P(x|z)

E\left [x|z\right ]=E\left [Wz+\mu +\varepsilon \right ]=Wz+\mu +0=Wz+\mu \\ Var\left [x|z\right ]=Var\left [Wz+\mu +\varepsilon \right ]=\sigma ^{2}I\\ \Rightarrow x|z\sim N(Wz+\mu ,\sigma ^{2}I)

  • P(x)

E[x]=E[Wz+\mu +\varepsilon ]=E[Wz+\mu ]+E[\varepsilon ]=\mu \\ Var[x]=Var[wz+\mu +\varepsilon ]=Var[Wz]+Var[\varepsilon ]=WW^{T}+\sigma ^{2}I\\ \Rightarrow x\sim N(\mu ,WW^{T}+\sigma ^{2}I)

  • P(z|x)

该问题和高斯分布|机器学习推导系列(二)中第六部分的问题是类似的。

\begin{pmatrix} x\\ z \end{pmatrix}\sim N\left (\begin{bmatrix} \mu \\ 0 \end{bmatrix} ,\begin{bmatrix} WW^{T}+\sigma ^{2}I & \Delta \\ \Delta ^{T}& I \end{bmatrix}\right )\\ \Delta =Cov(x,z)\\ =E[(x-\mu )(z-0)^{T}]\\ =E[(Wz+\mu +\varepsilon -\mu )z^{T}]\\ =E[(Wz+\varepsilon )z^{T}]\\ =E[Wzz^{T}+\varepsilon z^{T}]\\ =E[Wzz^{T}]+E[\varepsilon z^{T}]\\ =WE[zz^{T}]+E[\varepsilon ]E[z^{T}]\\ =WE[(z-0 )(z-0 )^{T}]+0\\ =WVar[z]\\ =W\\ 因此\begin{pmatrix} x\\ z \end{pmatrix}\sim N\left (\begin{bmatrix} \mu \\ 0 \end{bmatrix} ,\begin{bmatrix} WW^{T}+\sigma ^{2}I & W \\ W^{T}& I \end{bmatrix}\right )

利用高斯分布|机器学习推导系列(二)中第五部分的公式可以求解P(z|x)

P(z|x)\sim N(W^{T}(WW^{T}+\sigma ^{2}I)^{-1}(x-\mu ),I-W^{T}(WW^{T}+\sigma ^{2}I)^{-1}W)

3. 学习(learning)

使用EM算法求解,这里不做展示。

参考资料

ref:降维时为什么找最大的特征值对应的特征向量
ref:《模式识别与机器学习》