机器学习（十一）——因子分析法

11.因子分析法

对于高斯混合模型而言，需要足够的样本才能对模型进行拟合。但是如果出现样本数量远远少于特征数量的情况下怎么办？

首先，我们需要了解为什么在这种情况下单一高斯模型不行？

在高斯模型中，我们是利用最大似然来估计高斯分布中的参数（平均值、协方差）：
ϕ = 1 m ∑ i = 1 m 1 { y ( i ) = 1 } μ 0 = ∑ i = 1 m 1 { y ( i ) = 0 } x ( i ) ∑ i = 1 m 1 { y ( i ) = 0 } μ 1 = ∑ i = 1 m 1 { y ( i ) = 1 } x ( i ) ∑ i = 1 m 1 { y ( i ) = 1 } Σ = 1 m ∑ i = 1 m ( x ( i ) − μ y ( i ) ) ( x ( i ) − μ y ( i ) ) T \begin{aligned} \phi & = \frac {1}{m} \sum^m_{i=1}1\{y^{(i)}=1\}\\ \mu_0 & = \frac{\sum^m_{i=1}1\{y^{(i)}=0\}x^{(i)}}{\sum^m_{i=1}1\{y^{(i)}=0\}}\\ \mu_1 & = \frac{\sum^m_{i=1}1\{y^{(i)}=1\}x^{(i)}}{\sum^m_{i=1}1\{y^{(i)}=1\}}\\ \Sigma & = \frac{1}{m}\sum^m_{i=1}(x^{(i)}-\mu_{y^{(i)}})(x^{(i)}-\mu_{y^{(i)}})^T\\ \end{aligned} ϕμ0​μ1​Σ​=m1​i=1∑m​1{y(i)=1}=∑i=1m​1{y(i)=0}∑i=1m​1{y(i)=0}x(i)​=∑i=1m​1{y(i)=1}∑i=1m​1{y(i)=1}x(i)​=m1​i=1∑m​(x(i)−μy(i)​)(x(i)−μy(i)​)T​
分布写出来的具体形式如下：
$$\begin{array}{rl}p(y)& ={\varphi }^y(1-\varphi {)}^{1-y}\\ p(x|y=0)& =\frac{1}{(2\pi {)}^{n/2}|\Sigma {|}^{1/2}}exp(-\frac{1}{2}(x-{\mu }_0{)}^T{\Sigma }^{-1}(x-{\mu }_0))\\ p(x|y=1)& =\frac{1}{(2\pi {)}^{n/2}|\Sigma {|}^{1/2}}exp(-\frac{1}{2}(x-{\mu }_1{)}^T{\Sigma }^{-1}(x-{\mu }_1))\end{array}$$
我们可以发现这里的 $\Sigma$ （将样本考虑为空间中的向量，向量的数量n远远小于向量的维数m（特征数），所以n个向量不能完全表示m维空间中的所有向量，因此由这n个向量所表示的协方差矩阵必然不是满秩的）是一个奇异矩阵，这也就意味着其逆矩阵 ${\Sigma }^{-1}$ 不存在，而 $1/|\Sigma {|}^{1/2}=1/0$。 在样本数远远小于特征数的情况下，高斯分布会拟合的非常差。

1.$\Sigma$ 的约束条件

如果我们没有充足的数据来拟合一个完整的协方差矩阵（covariance matrix），就可以对矩阵空间 $\Sigma$ 给出某些约束条件（restrictions）。例如，我们可以选择去拟合一个对角（diagonal）的协方差矩阵 $\Sigma$。这样，读者很容易就能验证这样的一个协方差矩阵的最大似然估计（maximum likelihood estimate）可以由对角矩阵（diagonal matrix） $\Sigma$ 满足：
$${\Sigma }_{jj}=\frac{1}{m}\sum _{i=1}^m(x_j^{(i)}-{\mu }_j{)}^2$$
因此，${\Sigma }_{jj}$ 就是对数据中第 $j$ 个坐标位置的方差值的经验估计（empirical estimate）。

回忆一下，高斯模型的密度的形状是椭圆形的。对角线矩阵 $\Sigma$ 对应的就是椭圆长轴（major axes）对齐（axis- aligned）的高斯模型。

有时候，我们还要对这个协方差矩阵（covariance matrix）给出进一步的约束，不仅设为对角的（major axes），还要求所有对角元素（diagonal entries）都相等。这时候，就有 $\Sigma ={\sigma }^{2}I$，其中 ${\sigma }^2$ 是我们控制的参数。对这个 ${\sigma }^2$ 的最大似然估计则为：
$${\sigma }^2=\frac{1}{mn}\sum _{j=1}^n\sum _{i=1}^m(x_j^{(i)}-{\mu }_j{)}^2$$
这种模型对应的是密度函数为圆形轮廓的高斯模型（在二维空间也就是平面中是圆形，在更高维度当中就是球（spheres）或者超球体（hyperspheres））。

如果我们对数据要拟合一个完整的，不受约束的（unconstrained）协方差矩阵 $\Sigma$，就必须满足 $m\ge n+1$，这样才使得对 $\Sigma$ 的最大似然估计不是奇异矩阵（singular matrix）。在上面提到的两个约束条件之下，只要 $m\ge 2$，我们就能获得非奇异的（non-singular） $\Sigma$。

然而，将 $\Sigma$ 限定为对角矩阵，也就意味着对数据中不同坐标（coordinates）的 $x_i，x_j$建模都将是不相关的（uncorrelated），且互相独立（independent）。通常，还是从样本数据里面获得某些有趣的相关信息结构比较好。如果使用上面对 $\Sigma$ 的某一种约束，就可能没办法获取这些信息了。在本章讲义里面，我们会提到因子分析模型（factor analysis model），这个模型使用的参数比对角矩阵 $\Sigma$ 更多，而且能从数据中获得某些相关性信息（captures some correlations），但也不能对完整的协方差矩阵（full covariance matrix）进行拟合。

2.多重高斯模型的边缘分布和条件分布

在讲解因子分析之前，我们要先说一下一个联合多元高斯分布下的随机变量的条件和边缘分布。

假设我们有一个向量随机变量
$$x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$$
其中 $x_1\in R^r,x_2\in R^s$，因此 $x\in R^{r+s}$。设 $x\sim N(\mu ,\Sigma )$，则这两个参数为：
$$\mu =\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right]\Sigma =\left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right]$$
其中， ${\mu }_1\in R^r,{\mu }_2\in R^s,{\Sigma }_{11}\in R^{r\times r},{\Sigma }_{12}\in R^{r\times s}$，以此类推。由于协方差矩阵是对称的，所以有 ${\Sigma }_{12}={\Sigma }_{21}^T$。不难看出 $x_1$ 的期望 $E[x_1]={\mu }_1$

要用 $x_1$ 和 $x_2$的联合方差的概念：
$$\begin{array}{rl}Cov(x)& =\Sigma \\ & =\left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right]\\ & =E[(x-\mu )(x-\mu {)}^T]\\ & =E\left[\begin{array}{cc}\left(\begin{array}{c}{x}_1-{\mu }_1\\ x_2-{\mu }_2\end{array}\right)& {\left(\begin{array}{c}{x}_1-{\mu }_1\\ x_2-{\mu }_2\end{array}\right)}^T\end{array}\right]\\ & =\left[\begin{array}{cc}(x_1-{\mu }_1)(x_1-{\mu }_1{)}^T& (x_1-{\mu }_1)(x_2-{\mu }_2{)}^T\\ (x_2-{\mu }_2)(x_1-{\mu }_1{)}^T& (x_2-{\mu }_2)(x_2-{\mu }_2{)}^T\end{array}\right]\end{array}$$
所以我们可以发现：协方差 $Cov(x_1)=E[(x_1-{\mu }_1)(x_1-{\mu }_1)]={\Sigma }_{11}$。因此我们就可以给出一个正态分布 $x_1\sim N({\mu }_{,}{\Sigma }_{11})$ 来作为 $x_1$ 的边缘分布。

条件分布 $x_1|x_2\sim N({\mu }_{1|2},{\Sigma }_{1|2})$为（推导过程略）：
$$\begin{array}{rlr}& {\mu }_{1|2}={\mu }_1+{\Sigma }_{12}{\Sigma }_{22}^{-1}(x_2-{\mu }_2)& (1)\\ & {\Sigma }_{1|2}={\Sigma }_{11}-{\Sigma }_{12}{\Sigma }_{22}^{-1}{\Sigma }_{21}& (2)\end{array}$$

3.因子分析模型

我们可以假设，每个样本的背后都存在一个种类$z$，所以对于同一种类的样本数据，都是由种类$z$通过仿射变换映射到空间中的。我们制定在 $(x,z)$ 上的一个联合分布，如下所示，其中 $z\in R^k$ 是一个潜在随机变量（latent random variable）：
$$\begin{array}{rl}z& \sim N(0,I)\\ x|z& \sim N(\mu +\Lambda z,\Psi )\end{array}$$
上面的式子中，我们这个模型中的参数是向量 $\mu \in R^n$ ，矩阵 $\Lambda \in R^{n\times k}$，以及一个对角矩阵 $\Psi \in R^{n\times n}$。$k$ 的值通常都选择比 $n$ 小一点的。在 $\mu +\Lambda z^{(i)}$ 上加上协方差 $\Psi$ 作为噪音，就得到了 $x^{(i)}$。

所以得到使用下面的设定：
$$\begin{array}{rl}z& \sim N(0,I)\\ ϵ& \sim N(0,\Psi )\\ x& =\mu +\Lambda z+ϵ\end{array}$$
其中的 $ϵ$ 和 $z$ 是互相独立的。

然后咱们来确切地看看这个模型定义的分布。其中，随机变量 $z$ 和 $x$ 有一个联合高斯分布：
$$\left[\begin{array}{c}z\\ x\end{array}\right]\sim N({\mu }_{zx},\Sigma )$$
然后咱们要找到 ${\mu }_{zx}$ 和 $\Sigma$。

我们知道 $z$ 的期望 $E[z]=0$，这是因为 $z$ 服从的是均值为 $0$ 的正态分布 $z\sim N(0,I)$。 此外我们还知道：
$$\begin{array}{rl}E[x]& =E[\mu +\Lambda z+ϵ]\\ & =\mu +\Lambda E[z]+E[ϵ]\\ & =\mu \end{array}$$
综合以上这些条件，就得到了：
$${\mu }_{zx}=\left[\begin{array}{c}\vec{0}\\ \mu \end{array}\right]\Sigma =\left[\begin{array}{cc}{\Sigma }_{zz}& {\Sigma }_{zx}\\ {\Sigma }_{xz}& {\Sigma }_{xx}\end{array}\right]$$
所以：
$$\begin{array}{rl}{\Sigma }_{zz}& =I\\ & \\ {\Sigma }_{zx}& =E[(z-E[z])(x-E[x]{)}^T]\\ & =E[z(\mu +\Lambda z+ϵ-\mu {)}^T]\\ & =E[z{z}^T]{\Lambda }^T+E[z{ϵ}^T]\\ & ={\Lambda }^T\\ & \\ {\Sigma }_{zx}& ={\Sigma }_{zx}^T\\ & =\Lambda \\ & \\ {\Sigma }_{xx}& =E[(x-E[x])(x-E[x]{)}^T]\\ & =E[\mu +\Lambda z+ϵ-\mu )(\mu +\Lambda z+ϵ-\mu {)}^T]\\ & =E[\Lambda z{z}^T{\Lambda }^T+ϵz^T{\Lambda }^T+\Lambda z{ϵ}^T+ϵ{ϵ}^T]\\ & =\Lambda E[z{z}^T]{\Lambda }^T+E[ϵ{ϵ}^T]\\ & =\Lambda {\Lambda }^T+\Psi \end{array}$$
注意结论 $E[z{z}^T]=Cov(z)$（因为 $z$ 的均值为 $0$），而且 $E[z{ϵ}^T]=E[z]E[{ϵ}^T]=0$）（因为 $z$ 和 $ϵ$ 相互独立，因此乘积的期望等于期望的乘积）。

同时我们可以发现 $x$ 的边缘分布为 $x\sim N(\mu ,\Lambda {\Lambda }^T+\Psi )$。所以，我们下面使用EM算法求解最大似然(直接求最大似然过于复杂)

4.因子分析的EM算法

步骤E:

我们只需要计算$Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\mu ,\Lambda ,\Psi )$即可。我们利用上面求得的结果可知：
$$\begin{array}{rl}{\mu }_{z^{(i)}|x^{(i)}}& ={\Lambda }^T(\Lambda {\Lambda }^T+\Psi {)}^{-1}(x^{(i)}-\mu )\\ {\Sigma }_{z^{(i)}|x^{(i)}}& =I-{\Lambda }^T(\Lambda {\Lambda }^T+\Psi {)}^{-1}\Lambda \end{array}$$
所以：
$$Q_i(z^{(i)})=\frac{1}{(2\pi {)}^{k/2}|{\Sigma }_{z^{(i)}|x^{(i)}}{|}^{1/2}}exp(-\frac{1}{2}(z^{(i)}-{\mu }_{z^{(i)}|x^{(i)}}{)}^T{\Sigma }_{z^{(i)}|x^{(i)}}^{-1}(z^{(i)}-{\mu }_{z^{(i)}|x^{(i)}}))$$
步骤M:

我们这里需要去最大化下面这个关于参数 $\mu ,\Lambda$, $\Psi$ 的函数值：
$$\begin{array}{rl}& \sum _{i=1}^m{\int }_{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\mu ,\Lambda ,\Psi )}{Q_i(z^{(i)})}d{z}^{(i)}\\ & =\sum _{i=1}^m{\int }_{z^{(i)}}{Q}_i(z^{(i)})[logp(x^{(i)}|z^{(i)};\mu ,\Lambda ,\Psi )+logp(z^{(i)})-log{Q}_i(z^{(i)})]d{z}^{(i)}\\ & =\sum _{i=1}^m{E}_{z^{(i)}\sim Q_i}[logp(x^{(i)}|z^{(i)};\mu ,\Lambda ,\Psi )+logp(z^{(i)})-log{Q}_i(z^{(i)})]\end{array}$$
上面的等式中，$“z^{(i)}\sim Q_i”$ 这个下标（subscript），表示的意思是这个期望是关于从 $Q_i$ 中取得的 $z^{(i)}$ 的。在去除无关项之后我们可以发现，只需要最大化：
$$\begin{array}{rl}& \sum _{i=1}^{m}E[logp(x^{(i)}|z^{(i)};\mu ,\Lambda ,\Psi )]\\ & =\sum _{i=1}^{m}E[log\frac{1}{(2\pi {)}^{n/2}|\Psi {|}^{1/2}}exp(-\frac{1}{2}(x^{(i)}-\mu -\Lambda z^{(i)}{)}^T{\Psi }^{-1}(x^{(i)}-\mu -\Lambda z^{(i)}))]\\ & =\sum _{i=1}^{m}E[-\frac{1}{2}log|\Psi |-\frac{n}{2}log(2\pi )-\frac{1}{2}(x^{(i)}-\mu -\Lambda z^{(i)}{)}^T{\Psi }^{-1}(x^{(i)}-\mu -\Lambda z^{(i)})]\end{array}$$
我们首先计算关于 $\Lambda$ 的最大化，可以发现只有最后一项与$\Lambda$有关，我们取导数(梯度)：
$$\begin{array}{rl}{\nabla }_{\Lambda }& \sum _{i=1}^m-E[\frac{1}{2}(x^{(i)}-\mu -\Lambda z^{(i)}{)}^T{\Psi }^{-1}(x^{(i)}-\mu -\Lambda z^{(i)})]\\ & =\sum _{i=1}^m{\nabla }_{\Lambda }E[-tr\frac{1}{2}{z}^{(i)T}{\Lambda }^T{\Psi }^{-1}\Lambda z^{(i)}+tr{z}^{(i)T}{\Lambda }^T{\Psi }^{-1}(x^{(i)}-\mu )]\\ & =\sum _{i=1}^m{\nabla }_{\Lambda }E[-tr\frac{1}{2}{\Lambda }^T{\Psi }^{-1}\Lambda z^{(i)}{z}^{(i)T}+tr{\Lambda }^T{\Psi }^{-1}(x^{(i)}-\mu )z^{(i)T}]\\ & =\sum _{i=1}^{m}E[-{\Psi }^{-1}\Lambda z^{(i)}{z}^{(i)T}+{\Psi }^{-1}(x^{(i)}-\mu )z^{(i)T}]\end{array}$$
设置导数为 $0$，然后简化，就能得到：
$$\sum _{i=1}^m\Lambda E_{z^{(i)}\sim Q_i}[z^{(i)}{z}^{(i)T}]=\sum _{i=1}^m(x^{(i)}-\mu )E_{z^{(i)}\sim Q_i}[z^{(i)T}]$$
接下来，求解 $\Lambda$，就能得到：
$$\Lambda =(\sum _{i=1}^m(x^{(i)}-\mu )E_{z^{(i)}\sim Q_i}[z^{(i)T}])(\sum _{i=1}^m{E}_{z^{(i)}\sim Q_i}[z^{(i)}{z}^{(i)T}]{)}^{-1}$$

有一个很有意思的地方需要注意，上面这个等式和用最小二乘线性回归推出的正则方程有密切关系：
$$“{\theta }^T=(y^{T}X)(X^{T}X{)}^{-1}”$$

由于我们定义 $Q_i$ 是均值（mean）为 ${\mu }_{z^{(i)}|x^{(i)}}$，协方差（covariance）为 ${\Sigma }_{z^{(i)}|x^{(i)}}$ 的一个高斯分布，所以很容易能得到：
$$\begin{array}{rl}{E}_{z^{(i)}\sim Q_i}[z^{(i)T}]& ={\mu }_{z^{(i)}|x^{(i)}}^T\\ E_{z^{(i)}\sim Q_i}[z^{(i)}{z}^{(i)T}]& ={\mu }_{z^{(i)}|x^{(i)}}{\mu }_{z^{(i)}|x^{(i)}}^T+{\Sigma }_{z^{(i)}|x^{(i)}}\end{array}$$
上面第二个等式的推导依赖于下面这个事实：对于一个随机变量 $Y$，协方差 $Cov(Y)=E[Y{Y}^T]-E[Y]E[Y{]}^T$ ，所以 $E[Y{Y}^T]=E[Y]E[Y{]}^T+Cov(Y)$。把这个代入到等式，就得到了 $M$ 步骤中 $\Lambda$ 的更新规则：
$$\Lambda =(\sum _{i=1}^m(x^{(i)}-\mu ){\mu }_{z^{(i)}|x^{(i)}}^T)(\sum _{i=1}^m{\mu }_{z^{(i)}|x^{(i)}}{\mu }_{z^{(i)}|x^{(i)}}^T+{\Sigma }_{z^{(i)}|x^{(i)}}{)}^{-1}$$
最后，我们还可以发现，在 $M$ 步骤对参数 $\mu$ 和 $\Psi$ 的优化。不难发现其中的 $\mu$ 为：
$$\mu =\frac{1}{m}\sum _{i=1}^m{x}^{(i)}$$
由于这个值不随着参数的变换而改变（也就是说，和 $\Lambda$ 的更新不同，这里等式右侧不依赖 $Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\mu ,\Lambda ,\Psi )$，这个 $Qi(z^{(i)})$ 是依赖参数的），这个只需要计算一次就可以，在算法运行过程中，也不需要进一步更新。类似地，对角矩阵 $\Psi$ 也可以通过计算下面这个式子来获得：
$$\Phi =\frac{1}{m}\sum _{i=1}^m{x}^{(i)}{x}^{(i)T}-x^{(i)}{\mu }_{z^{(i)}|x^{(i)}}^T{\Lambda }^T-\Lambda {\mu }_{z^{(i)}|x^{(i)}}{x}^{(i)T}+\Lambda ({\mu }_{z^{(i)}|x^{(i)}}{\mu }_{z^{(i)}|x^{(i)}}^T+{\Sigma }_{z^{(i)}|x^{(i)}}){\Lambda }^T$$
然后只需要设 ${\Psi }_{ii}={\Phi }_{ii}$（也就是说，设 $\Psi$ 为一个仅仅包含矩阵 $\Phi$ 中对角线元素的对角矩阵）。