$$f(x)=\frac{1}{\sqrt{2\pi^k\left|\Sigma \right|}}\exp(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu))\tag{3}$$

其中

x
$x$是



k

$k$维向量，

μ
$\mu$是

k
$k$维向量，



Σ

$\Sigma$是

k∗k
$k*k$协方差矩阵。

很明显在多元高斯分布模型下，参数是

μzx,Σ
$\mu_{zx},\Sigma$——它们是由

x,z
$x,z$的联合分布生成的，所以我们可以用我们的原始参数

μ,Λ,ψ
$\mu,\Lambda,\psi$来表示

μzx,Σ
$\mu_{zx},\Sigma$，求得

x
$x$的边缘分布，再把相关参数带入式（3），这就得到了关于我们参数的概率分布，然后就可以通过最大似然估计来求取我们的参数。

求取$\mu_{zx}$

$\mu_{zx}$是$x,z$联合分布的期望值（期望的定义：所有结果*相应概率的总和）：

μzx=E[zx]=[E(z)E(x)](4)

$$\mu_{zx}=E\left[ \begin{matrix}z\\x \end{matrix}\right]=\left[ \begin{matrix}E(z)\\E(x) \end{matrix}\right]\tag{4}$$

由$z$~$N(0,I)$我们可以简单获得$E(z)=0$。
类似地由$\varepsilon$~$N(0,\psi)$，$x=\mu+\Lambda z+\varepsilon$，$\mu$是一个常数，我们有：

E[x]=E[μ+Λz+ε]=E[μ]+ΛE[z]+E[ε]=μ+0+0=μ(5)

$$\begin{align}E[x]&=E[\mu+\Lambda z+\varepsilon]\\&=E[\mu]+\Lambda E[z]+E[\varepsilon]\\&=\mu+0+0\\&=\mu\tag{5}\end{align}$$
所以：

μzx=[0⃗ μ](6)

$$\mu_{zx}=\left[ \begin{matrix}\vec 0\\\mu \end{matrix}\right]\tag{6}$$

求取$\Sigma$

$\Sigma$是$x,z$联合分布的协方差矩阵。
方差，度量随机变量与期望之间的偏离程度，定义如下：

Var(X)=E((X−E(X))2)=E(X2)−(E(X)2)(7)

$$Var(X)=E((X-E(X))^2)=E(X^2)-(E(X)^2)\tag{7}$$
协方差，两个变量总体误差的期望，定义如下：

Cov(X,Y)=E((X−E(X))(Y−E(Y)))(8)

$$Cov(X,Y)=E((X-E(X))(Y-E(Y)))\tag{8}$$
协方差、方差、期望之间的一些相互关系如下：

Cov(X,X)=Cov(X)=Var(X)=E(XXT)=σ2(9)

$$Cov(X,X)=Cov(X)=Var(X)=E(XX^T)=\sigma^2\tag{9}$$

下面开始求取$\Sigma$。

Σ=Cov[zx]=[ΣzzΣxzΣzxΣxx]=E[(z−E(z))(z−E(z))T(x−E(x))(z−E(z))T(z−E(z))(x−E(x))T(x−E(x))(x−E(x))T](10)

$$\begin{align}\Sigma&=Cov\left[ \begin{matrix}z\\x \end{matrix}\right]\\&=\left[ \begin{matrix}\Sigma_{zz}&\Sigma_{zx}\\\Sigma_{xz}&\Sigma_{xx} \end{matrix}\right]\\&=E\left[ \begin{matrix}(z-E(z))(z-E(z))^T& (z-E(z))(x-E(x))^T \\(x-E(x))(z-E(z))^T&(x-E(x))(x-E(x))^T \end{matrix}\right]\tag{10}\end{align}$$

由$z$~$N(0,I)$，可以简单得到：

Σzz=Cov(z)=σ2=I(11)

$$\Sigma_{zz}=Cov(z)=\sigma^2=I\tag{11}$$
由 $\varepsilon$~ $N(0,\psi)$， $x=\mu+\Lambda z+\varepsilon$， $E(x)=\mu$，并且 $z$和$\varepsilon$是相互独立，有：

Σzx=E[(z−E(z))(x−E(x))T]=E[(z−0)(μ+Λz+ε−μ)T]=E[zzT]ΛT+E[zεT]=IΛT+0=ΛT(12)

$$\begin{align}\Sigma_{zx}&=E[ (z-E(z))(x-E(x))^T]\\&=E[(z-0)(\mu+\Lambda z+\varepsilon-\mu)^T]\\&=E[zz^T]\Lambda^T+E[z\varepsilon^T]\\&=I\Lambda^T+0\\&=\Lambda^T\tag{12}\end{align}$$
类似地，我们可以得到：

Σxx=E[(x−E(x))(x−E(x))T]=E[(μ+Λz+ε−μ)(μ+Λz+ε−μ)T]=ΛE[zzT]ΛT+E[εεT]=ΛIΛT+ψ=ΛΛT+ψ(13)

$$\begin{align}\Sigma_{xx}&=E[ (x-E(x))(x-E(x))^T]\\&=E[(\mu+\Lambda z+\varepsilon-\mu)(\mu+\Lambda z+\varepsilon-\mu)^T]\\&=\Lambda E[zz^T]\Lambda^T+E[\varepsilon\varepsilon^T]\\&=\Lambda I\Lambda^T + \psi\\&=\Lambda \Lambda^T + \psi\tag{13}\end{align}$$

用最大似然估计法求解参数

经过上面的步骤，我们就把$\mu_{zx},\Sigma$用我们的参数$\mu,\Lambda,\psi$表示出来了：

[zx]
$\left[ \begin{matrix}z\\x \end{matrix}\right]$~

N(μzx,Σ)
$N(\mu_{zx},\Sigma)$~

N([0⃗ μ],[IΛΛTΛΛT+ψ])
$N(\left[ \begin{matrix}\vec 0\\\mu \end{matrix}\right],\left[ \begin{matrix}I&\Lambda^T\\\Lambda&\Lambda \Lambda^T + \psi \end{matrix}\right])$

然后我们可以求得 $x$的边缘分布：


$x$~ $N(\mu,\Lambda \Lambda^T + \psi)$
因此，给定一个训练集 $\begin{Bmatrix} x^{(i)};i=1,2,\cdots,m \end{Bmatrix}$，把参数带入式（3），我们可以写出下面的似然函数：

l(μ,Λ,ψ)=log∏i=1m12πn∣∣ΛΛT+ψ∣∣−−−−−−−−−−−√exp(−12(x(i)−μ)T(ΛΛT+ψ)−1(x(i)−μ))(14)

$$l(\mu,\Lambda,\psi)=\log\prod_{i=1}^m\frac{1}{\sqrt{2\pi^n\left|\Lambda \Lambda^T + \psi \right|}}\exp(-\frac{1}{2}(x^{(i)}-\mu)^T (\Lambda \Lambda^T + \psi)^{-1}(x^{(i)}-\mu))\tag{14}$$
对此似然函数做最大似然估计，就能求得我们的参数。

1.4、因子分析的EM求解

可以感受到，直接对这个似然函数求解是很困难的，在这个情况下，用EM算法就登场了——当一个似然函数难以直接求解其最大值的时候，可以通过EM算法不断建立下界、最大化下界的方式不断逼近该似然函数真实的最大值，当EM算法收敛，我们就认为已经求得了此最大值。

E-step

对于EM算法的E-step，我们有：

Qi(z(i)):=p(z(i)∣x(i);μ,Λ,ψ)(15)

$$Q_i(z^{(i)}):=p(z^{(i)} \mid x^{(i)};\mu,\Lambda,\psi)\tag{15}$$
进一步地：

Qi(z(i))=12πk∣∣Σz(i)∣x(i)∣∣−−−−−−−−−−√exp(−12(z(i)−μz(i)∣x(i))TΣ−1z(i)∣x(i)(z(i)−μz(i)∣x(i)))(16)

$$Q_i(z^{(i)})=\frac{1}{\sqrt{2\pi^k\left|\Sigma_{z^{(i)}\mid x^{(i)}} \right|}}\exp(-\frac{1}{2}(z^{(i)}-\mu_{z^{(i)}\mid x^{(i)}})^T\Sigma_{z^{(i)}\mid x^{(i)}}^{-1}(z^{(i)}-\mu_{z^{(i)}\mid x^{(i)}}))\tag{16}$$
其中：

μz(i)∣x(i)Σz(i)∣x(i)=ΛT(ΛΛT+ψ)−1(x(i)−μ)=I−ΛT(ΛΛT+ψ)−1Λ(17)

$$\begin{align}\mu_{z^{(i)}\mid x^{(i)}}&=\Lambda^T(\Lambda\Lambda^T+\psi)^{-1}(x^{(i)}-\mu)\\\Sigma_{z^{(i)}\mid x^{(i)}}&=I-\Lambda^T(\Lambda\Lambda^T+\psi)^{-1}\Lambda\tag{17}\end{align}$$
$\mu_{z^{(i)}\mid x^{(i)}},\Sigma_{z^{(i)}\mid x^{(i)}}$是讲义与课上直接给出的，这里也不进行推导。

M-step

在M-step中，我们需要最大化如下公式来求取参数$\mu,\Lambda,\psi$：

∑i=1m∫z(i)Qi(z(i))logp(x(i),z(i);μ,Λ,ψ)Qi(z(i))dz(i)(18)

$$\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)})}dz^{(i)}\tag{18}$$

视为期望，打开log

在这里，因为$z$是连续的，所以使用积分；如果是离散的，则使用累加。
并且，积分部分可以当成$z$服从$Q$分布时，函数$\log \frac{p(x^{(i)},z^{(i)};\mu,\Lambda,\psi)}{Q_i(z^{(i)})}$的期望，这里将会用E表示，省略$z^{(i)}$~$Q_i$的下标；对于函数中$x,z$的联合分布，我们可以用贝叶斯公式把它打开$p(x,z)=p(x\mid x)p(z)$；为了方便计算我们还要把log函数打开——经过这些分析，我们有如下推导：

∑i=1m∫z(i)Qi(z(i))logp(x(i),z(i);μ,Λ,ψ)Qi(z(i))dz(i)=∑i=1mE[logp(x(i),z(i);μ,Λ,ψ)Qi(z(i))]=∑i=1mE[logp(x(i)∣z(i);μ,Λ,ψ)p(z(i))Qi(z(i))]=∑i=1mE[logp(x(i)∣z(i);μ,Λ,ψ)+logp(z(i))−logQi(z(i))](19)

$$\begin{align}&\quad\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)})}dz^{(i)}\\&=\sum_{i=1}^m E [\log \frac{p(x^{(i)},z^{(i)};\mu,\Lambda,\psi)}{Q_i(z^{(i)})}]\\&=\sum_{i=1}^m E [\log \frac{p(x^{(i)}\mid z^{(i)};\mu,\Lambda,\psi)p(z^{(i)})}{Q_i(z^{(i)})}]\\&=\sum_{i=1}^m E [\log p(x^{(i)}\mid z^{(i)};\mu,\Lambda,\psi)+\log p(z^{(i)})-\log Q_i(z^{(i)})]\tag{19}\end{align}$$

去掉无关项后带入具体分布

这就比较清爽了，然后，记住我们的目标是求得参数$\mu,\Lambda,\psi$，但是它们不能一起求解，所以下面以参数$\Lambda$为例，对公式进行求解——在式（19）中，对参数$\Lambda$求偏导。另外式（19）中的$p(z^{(i)}$与$Q_i(z^{(i)})$与$\Lambda$无关，可以忽略掉，所以实际上就是对下式求偏导：

∑i=1mE[logp(x(i)∣z(i);μ,Λ,ψ)](20)

$$\begin{align}\sum_{i=1}^m E [\log p(x^{(i)}\mid z^{(i)};\mu,\Lambda,\psi)]\tag{20}\end{align}$$

在对式（20）求偏导之前，还可以对其进行一些处理——由式（2），并且$x\mid z$服从多元高斯分布，所以有：

∑i=1mE[logp(x(i)∣z(i);μ,Λ,ψ)]=∑i=1mE[log12πn|ψ|−−−−−−√exp(−12(x(i)−(μ+Λz(i)))Tψ−1(x(i)−(μ+Λz(i))))]=∑i=1mE[−12log|ψ|−n2log(2π)−12(x(i)−μ−Λz(i))Tψ−1(x(i)−μ−Λz(i)))](21)

$$\begin{align}&\quad \sum_{i=1}^m E [\log p(x^{(i)}\mid z^{(i)};\mu,\Lambda,\psi)]\\&=\sum_{i=1}^m E [\log \frac{1}{\sqrt{2\pi^n\left|\psi \right|}}\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\left|\psi \right|-\frac{n}{2}\log(2\pi)-\frac{1}{2}(x^{(i)}-\mu-\Lambda z^{(i)})^T\psi^{-1}(x^{(i)}-\mu-\Lambda z^{(i)}))]\tag{21}\end{align}$$

去掉无关项后求偏导

同样地，我们的目标是与$\Lambda$有关的项，所以忽略掉前面的无关项之后，我们实际上是对下式求偏导并求解：

∇Λ∑i=1mE[−12(x(i)−μ−Λz(i))Tψ−1(x(i)−μ−Λz(i)))]=∑i=1m∇Λ−E[12((x(i)T−μT)ψ−1A−z(i)TΛTψ−1B)((x(i)−μ)C−Λz(i)D)]=∑i=1m∇Λ−E[12(AC−AD−BC+BD)](22)

$$\begin{align}&\quad \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[\frac{1}{2}(\underbrace{(x^{(i)^T}-\mu^T)\psi^{-1}}_A-\underbrace{z^{(i)^T}\Lambda^T\psi^{-1}}_B)(\underbrace{(x^{(i)}-\mu)}_C-\underbrace{\Lambda z^{(i)}}_D)]\\&=\sum_{i=1}^m \nabla_\Lambda -E[\frac{1}{2}(AC-AD-BC+BD)]\tag{22}\end{align}$$
打开后我们可以发现， $AC$这一项是与 $\Lambda$无关的，把这一项忽略掉，所以由式（22）继续推导有：

∑i=1m∇Λ−E[12(−AD−BC+BD)]=∑i=1m∇Λ−E[12(−(x(i)T−μT)ψ−1Λz(i)E−z(i)TΛTψ−1(x(i)−μ)F+z(i)TΛTψ−1Λz(i))](23)

$$\begin{align}&\quad\sum_{i=1}^m \nabla_\Lambda -E[\frac{1}{2}(-AD-BC+BD)]\\&=\sum_{i=1}^m \nabla_\Lambda -E[\frac{1}{2}(-\underbrace{(x^{(i)^T}-\mu^T)\psi^{-1} \Lambda z^{(i)}}_E-\underbrace{z^{(i)^T}\Lambda^T\psi^{-1}(x^{(i)}-\mu)}_F+z^{(i)^T}\Lambda^T\psi^{-1}\Lambda z^{(i)})]\tag{23}\end{align}$$
因为期望是一个常数，又因为 $a=tr(a)$，所以可以直接对上式求迹；
因为 $trA=trA^T$，可以对E求转置，又因为对角矩阵的转置是它本身—— $(\psi^{-1})^T=\psi^{-1}$，所以有 $trE=trE^T=trF$，对式（23）继续推导有：

∑i=1m∇Λ−E[12(−(x(i)T−μT)ψ−1Λz(i)E−z(i)TΛTψ−1(x(i)−μ)F+z(i)TΛTψ−1Λz(i))]=∑i=1m∇Λ−E[tr12(−z(i)TΛTψ−1(x(i)−μ)ET−z(i)TΛTψ−1(x(i)−μ)F+z(i)TΛTψ−1Λz(i))]=∑i=1m∇ΛE[−tr12z(i)TΛTψ−1Λz(i)+trz(i)TΛTψ−1(x(i)−μ)](24)

$$\begin{align}&\quad\sum_{i=1}^m \nabla_\Lambda -E[\frac{1}{2}(-\underbrace{(x^{(i)^T}-\mu^T)\psi^{-1} \Lambda z^{(i)}}_E-\underbrace{z^{(i)^T}\Lambda^T\psi^{-1}(x^{(i)}-\mu)}_F+z^{(i)^T}\Lambda^T\psi^{-1}\Lambda z^{(i)})]\\&=\sum_{i=1}^m \nabla_\Lambda -E[tr\frac{1}{2}(-\underbrace{z^{(i)^T}\Lambda^T\psi^{-1}(x^{(i)}-\mu)}_{E^T}-\underbrace{z^{(i)^T}\Lambda^T\psi^{-1}(x^{(i)}-\mu)}_F+z^{(i)^T}\Lambda^T\psi^{-1}\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)]\tag{24}\end{align}$$
然后利用 $trAB=trBA$，把式（25）中的 $z^{(i)^T}$放到它们自己的后面，再把求导切换到期望里面——求导是针对 $\Lambda$，期望是针对 $z^{(i)}$，所以是可以切换的：

∑i=1m∇ΛE[−tr12z(i)TΛTψ−1Λz(i)+trz(i)TΛTψ−1(x(i)−μ)]=∑i=1m∇ΛE[−tr12ΛTψ−1Λz(i)z(i)T+trΛTψ−1(x(i)−μ)z(i)T]=∑i=1m(E[−∇Λtr12ΛTψ−1Λz(i)z(i)T]+E[∇ΛtrΛTψ−1(x(i)−μ)z(i)T])(25)

$$\begin{align} &\quad\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\left( E[-\nabla_\Lambda tr\frac{1}{2}\Lambda^T\psi^{-1}\Lambda z^{(i)}z^{(i)^T}]+E[\nabla_\Lambda tr\Lambda^T\psi^{-1}(x^{(i)}-\mu)z^{(i)^T}]\right)\tag{25}\end{align}$$
对于式（25），先用矩阵的迹的性质 $\nabla_{A^T}f(A)=(\nabla_A f(A))^T$处理一下：

∑i=1m(E[−(∇ΛTtr12ΛTψ−1Λz(i)z(i)T)T]+E[(∇ΛTtrΛTψ−1(x(i)−μ)z(i)T)T])(26)

$$\sum_{i=1}^m\left( E\left[-\left(\nabla_{\Lambda^T} tr\frac{1}{2}\Lambda^T\psi^{-1}\Lambda z^{(i)}z^{(i)^T}\right)^T\right]+E\left[\left(\nabla_{\Lambda^T} tr\Lambda^T\psi^{-1}(x^{(i)}-\mu)z^{(i)^T}\right)^T\right]\right)\tag{26}$$
对式（26）的第一项使用 $\nabla_A trABA^TC=CAB+C^TAB^T$的性质，对第二项使用 $\nabla_A trAB=B^T$ 的性质：

∑i=1m⎛⎝⎜⎜E⎡⎣⎢⎢−⎛⎝⎜∇ΛTAtr12ΛTAψ−1BΛATz(i)z(i)TC⎞⎠⎟T⎤⎦⎥⎥+E⎡⎣⎢⎢⎛⎝⎜∇ΛTAtrΛTAψ−1(x(i)−μ)z(i)TB⎞⎠⎟T⎤⎦⎥⎥⎞⎠⎟⎟=∑i=1mE((−12z(i)z(i)TΛTψ−1−12(z(i)z(i)T)TΛT(ψ−1)T)T+((ψ−1(x(i)−μ)z(i)T)T)T)=∑i=1mE((−z(i)z(i)TΛTψ−1)T+ψ−1(x(i)−μ)z(i)T)=∑i=1mE(−ψ−1Λz(i)z(i)T+ψ−1(x(i)−μ)z(i)T)(27)

$$\begin{align} &\quad\sum_{i=1}^m\left( E\left[-\left(\nabla_{\underbrace{\Lambda^T}_A} tr\frac{1}{2}\underbrace{\Lambda^T}_A\underbrace{\psi^{-1}}_B\underbrace{\Lambda}_{A^T} \underbrace{z^{(i)}z^{(i)^T}}_C\right)^T\right]+E\left[\left(\nabla_{\underbrace{\Lambda^T}_A} tr\underbrace{\Lambda^T}_A\underbrace{\psi^{-1}(x^{(i)}-\mu)z^{(i)^T}}_B\right)^T\right]\right)\\&=\sum_{i=1}^mE\left(\left(-\frac{1}{2}z^{(i)}z^{(i)^T}\Lambda^T\psi^{-1}-\frac{1}{2}(z^{(i)}z^{(i)^T})^T\Lambda^T(\psi^{-1})^T\right)^T+\left((\psi^{-1}(x^{(i)}-\mu)z^{(i)^T})^T\right)^T\right)\\&=\sum_{i=1}^mE\left(\left(-z^{(i)}z^{(i)^T}\Lambda^T\psi^{-1}\right)^T+\psi^{-1}(x^{(i)}-\mu)z^{(i)^T}\right)\\&=\sum_{i=1}^mE\left(-\psi^{-1}\Lambda z^{(i)}z^{(i)^T}+\psi^{-1}(x^{(i)}-\mu)z^{(i)^T}\right)\tag{27}\end{align}$$
令式（27）=0，并化简，就可以求得参数 $\Lambda$：

⟹⟹⟹∑i=1mE(−ψ−1Λz(i)z(i)T+ψ−1(x(i)−μ)z(i)T)=0∑i=1m−ψ−1ΛE[z(i)z(i)T]+∑i=1m−ψ−1(x(i)−μ)E[z(i)T]=0∑i=1mΛE[z(i)z(i)T]=∑i=1m(x(i)−μ)E[z(i)T]Λ=(∑i=1m(x(i)−μ)E[z(i)T])(∑i=1mE[z(i)z(i)T])−1(28)

$$\begin{align} & \sum_{i=1}^mE\left(-\psi^{-1}\Lambda z^{(i)}z^{(i)^T}+\psi^{-1}(x^{(i)}-\mu)z^{(i)^T}\right)=0\\\Longrightarrow &\sum_{i=1}^m-\psi^{-1}\Lambda E[z^{(i)}z^{(i)^T}]+\sum_{i=1}^m -\psi^{-1} (x^{(i)}-\mu)E[z^{(i)^T}]=0\\\Longrightarrow &\sum_{i=1}^m\Lambda E[z^{(i)}z^{(i)^T}]= \sum_{i=1}^m(x^{(i)}-\mu)E[z^{(i)^T}]\\\Longrightarrow & \Lambda=\left(\sum_{i=1}^m(x^{(i)}-\mu)E[z^{(i)^T}]\right)\left(\sum_{i=1}^m E[z^{(i)}z^{(i)^T}]\right)^{-1}\tag{28}\end{align}$$
我们发现，这里的公式与线性回归中最小二乘法的矩阵形式相似。
相似原因：在因子分析中， $x$是$z$的线性函数，在E-step中给出 $z$的$Q$分布之后，在M-tep中寻找 $x$与$z$的映射关系 $\Lambda$；在线性回归的最小二乘中，也是寻找 $x$与$y$的线性关系。
不同之处：最小二乘只用到了 $z$的最优估计，因子分析还用到了$z^{(i)}z^{(i)^T}$的估计。

对于参数$\Lambda$，这里还有未知数$E[z^{(i)^T}]$与$E[z^{(i)}z^{(i)^T}]$，并且此处的期望是在$z^{(i)}$服从$Q_i$ 前提下计算的，所以对于前者，通过式（17）我们有：

E[z(i)T]=μTz(i)∣x(i)(29)

$$E[z^{(i)^T}]=\mu_{z^{(i)}\mid x^{(i)}}^T\tag{29}$$
对于后者，由式（7）~式（9）方差与协方差的性质，我们有：

Cov(z)⟹E[z(i)z(i)T]=E(zzT)−E(z)E(zT)=E(z(i))E(z(i)T)+Cov(z)=μz(i)∣x(i)μTz(i)∣x(i)+Σz(i)∣x(i)(30)

$$\begin{align} Cov(z)&=E(zz^T)-E(z)E(z^T)\\\Longrightarrow E[z^{(i)}z^{(i)^T}]&=E(z^{(i)})E(z^{(i)^T})+Cov(z)\\&=\mu_{z^{(i)}\mid x^{(i)}}\mu_{z^{(i)}\mid x^{(i)}}^T+\Sigma_{z^{(i)}\mid x^{(i)}}\tag{30}\end{align}$$
注意这里的 $E[z^{(i)}z^{(i)^T}]$ 不仅仅等于 $E(z^{(i)})E(z^{(i)^T})$，后面还有加上一个后验概率 $p(z\mid x)$协方差，要特别注意。

到这里，我们就可以把参数$\Lambda$的最终形式给出来了：

Λ=(∑i=1m(x(i)−μ)μTz(i)∣x(i))(∑i=1mμz(i)∣x(i)μTz(i)∣x(i)+Σz(i)∣x(i))−1(31)

$$\Lambda=\left(\sum_{i=1}^m(x^{(i)}-\mu)\mu_{z^{(i)}\mid x^{(i)}}^T\right)\left(\sum_{i=1}^m \mu_{z^{(i)}\mid x^{(i)}}\mu_{z^{(i)}\mid x^{(i)}}^T+\Sigma_{z^{(i)}\mid x^{(i)}}\right)^{-1}\tag{31}$$
另外对于其他两个参数 $\Lambda,\psi$，使用相同的方法可以求得，这里直接给出结果
：

μ=1m∑i=1mx(i)(32)

$$\mu=\frac{1}{m}\sum_{i=1}^{m}x^{(i)}\tag{32}$$

ϕ=1m∑i=1mx(i)x(i)T−x(i)μTz(i)∣x(i)ΛT−Λμz(i)∣x(i)x(i)T+Λμz(i)∣x(i)μTz(i)∣x(i)+Σz(i)∣x(i)ΛT(33)

$$\phi=\frac{1}{m}\sum_{i=1}{m}x^{(i)}x^{(i)^T}-x^{(i)}\mu_{z^{(i)}\mid x^{(i)}}^T\Lambda^T-\Lambda\mu_{z^{(i)}\mid x^{(i)}}x^{(i)^T}+\Lambda\mu_{z^{(i)}\mid x^{(i)}}\mu_{z^{(i)}\mid x^{(i)}}^T+\Sigma_{z^{(i)}\mid x^{(i)}}\Lambda^T\tag{33}$$

注意这里的参数是$\phi$不是$\psi$，得到$\phi$之后还需要将$\psi_{ii}=\phi_{ii}$，因为$\phi$不是对角矩阵，所以只需要取$\phi$对角线上的值即可。

2、主成分分析（Principal component analysis，简称PCA

PCA的意义

PCA技术的一大好处是对数据进行降维的处理。我们可以对新求出的“主元”向量的重要性进行排序，根据需要取前面最重要的部分，将后面的维数省去，可以达到降维从而简化模型或是对数据进行压缩的效果。同时最大程度的保持了原有数据的信息。

PCA将n个特征降维到k个，可以用来进行数据压缩，如果100维的向量最后可以用10维来表示，那么压缩率为90%。同样图像处理领域的KL变换使用PCA做图像压缩。但PCA要保证降维后，还要保证数据的特性损失最小。

预处理

运行PCA算法之前，数据一般要进行预处理。预处理步骤如下：
①令$\mu=\frac{1}{m}\sum_{i=1}^{m}x^{(i)}$
②用$(x^{(i)}-\mu)$代替$x^{(i)}$
③令$\sigma^2=\frac{1}{m}\sum_{i=1}^{m}(x^{(i)}_j)^2$
④用$x^{(i)}_j/\sigma^2$代替$x^{(i)}_j$
步骤①与②将数据的均值变为0；
步骤③与④将数据每个维度的方差变为1，使每个维度都在同一个维度下被度量。
如果已知数据均值为0，或者数据已在同样一个维度下，就无需进行上面的步骤。

最大方差理论解释PCA

Ng在课上说，PCA有9到10种解释方法，这里仅用其中的最大方差理论来解释PCA
在信号处理中认为信号具有较大的方差，噪声有较小的方差，信噪比就是信号与噪声的方差比，越大越好。如前面的图，样本在横轴上的投影方差较大，在纵轴上的投影方差较小，那么认为纵轴上的投影是由噪声引起的。
因此我们认为，最好的k维特征是将n维样本点转换为k维后，每一维上的样本方差都很大。

下面举例来说明如何寻找主方向。

图六

图六左图是经过预处理的样本点，均值为0，特征方差归一；
图六的中图和右图是将样本投影到某个维度上的表示， 该维度用一条过原点的直线表示，前处理的过程实质是将原点移到样本点的中心点。




图七

图七中，原点到样本在直线上的投影的距离 $x^{(i)^T}u$既是方差；样本点到直线上的距离是平方误差。
我们可以直观感受到，图六中间的图方差和是最大的，平方误差是最小的。下面给出用最大方差理论寻找该最佳方向的公式定义。

形式化

设$x^{(i)}$是数据集中的点，$u$是要求解的单位向量，那么方差最大化可以形式化为最大化下式：

1m∑i=1m(x(i)Tu)2=1m∑i=1muTx(i)x(i)Tu=uT(1m∑i=1mx(i)x(i)T)u(34)

$$\frac{1}{m}\sum_{i=1}^{m}(x^{(i)^T}u)^2=\frac{1}{m}\sum_{i=1}^{m}u^T x^{(i)}x^{(i)^T}u=u^T\left(\frac{1}{m}\sum_{i=1}^{m}x^{(i)}x^{(i)^T}\right)u\tag{34}$$
因为归一化后的数据，投影值的均值也为0，所以在方差计算中直接平方。
同时这个式子还有一个约束条件，即$\left \| u \right \|_2=1$。熟悉的最大化某个带约束的函数，引入拉格朗日乘子来求解：

l=uT(1m∑i=1mx(i)x(i)T)u−λ(∥u∥2−1)=uTΣu−λ(uTu−1)(35)

$$l=u^T\left(\frac{1}{m}\sum_{i=1}^{m}x^{(i)}x^{(i)^T}\right)u-\lambda(\left \| u \right \|_2-1)=u^T\Sigma u -\lambda(u^Tu-1)\tag{35}$$
对 $u$求导：

∇ul=∇u(uTΣu−λ(uTu−1))=∇utr(uTΣu)−λ∇utr(uTu)=(∇uTtr(uTΣu))T−(λ∇uTtr(uTu))T=((Σu)T)T−(λuT)T=Σu−λu(36)

$$\begin{align}\nabla_ul&=\nabla_u\left(u^T\Sigma u -\lambda(u^Tu-1)\right)\\&=\nabla_utr(u^T\Sigma u)-\lambda\nabla_utr(u^Tu)\\&=\left(\nabla_{u^T}tr(u^T\Sigma u)\right)^T-\left(\lambda\nabla_{u^T}tr(u^Tu)\right)^T\\&=\left((\Sigma u)^T\right)^T-\left(\lambda u^T\right)^T\\&=\Sigma u-\lambda u\tag{36}\end{align}$$
这里主要运用了 $\nabla_{A^T}f(A)=(\nabla_A f(A))^T$与 $\nabla_A trAB=B^T$的性质，处理方式与上面因子分析中式（26）与（27）的方式类似，不再赘述。
令导数为0，可知 $u$是$\Sigma$的特征向量—— $(A-\lambda)x=0$看起来会眼熟一些。
因为 $\Sigma=\frac{1}{m}\sum_{i=1}^{m}x^{(i)}x^{(i)^T}$是对称矩阵，因而可得到 $n$个相正交的特征向量：$\left[ \begin{matrix}u_1\\u_2\\\vdots\\u_n \end{matrix}\right]$
此时如何达到降维的效果？选取 最大的$k$个（$k<n$）特征值对应的特征向量即可。$u\in R^n$，$x^{(i)}\in R^n$，降维后的数据为：

y(i)=⎡⎣⎢⎢⎢⎢⎢⎢uT1x(i)uT2x(i)⋮uTkx(i)⎤⎦⎥⎥⎥⎥⎥⎥∈Rk(37)

$$y^{(i)}=\left[ \begin{matrix}u_1^T x^{(i)}\\u_2^T x^{(i)}\\\vdots\\u_k^T x^{(i)} \end{matrix}\right] \in R^k\tag{37}$$
这样，对于$n$维的原始数据$x^{(i)}$，我们有了更低维的$k$维$y^{(i)}$去表示它。
以上就是PCA的数学表达式。

3、独立成分分析（Independent component analysis，简称ICA）

解决的问题

ICA解决的是原始数据分解的问题，更具体一些的例子就是将声音信号分离的问题。

经典的鸡尾酒宴会问题（cocktail party problem）。假设在party中有n个人，他们可以同时说话，我们也在房间中一些角落里共放置了个声音接收器用来记录声音。宴会过后，我们从$n$个麦克风中得到了一组数据$\begin{Bmatrix} x^{(i)};i=1,2,\cdots,m \end{Bmatrix}$，i表示采样的时间顺序，也就是说共得到了m组采样，每一组采样都是n维的。我们的目标是单单从这m组采样数据中分辨出每个人说话的信号。
同时我们假设$n$个信号源为$s=(s_1,s_2,\cdots,s_n)^T$，$s\in R^n$，每一维都是一个人的声音信号，每个人发出的声音信号独立。A是一个未知的混合矩阵（mixing matrix），用来组合叠加信号s，那么：

x=⎡⎣⎢|x(1)||x(2)|⋯|x(m)|⎤⎦⎥=⎡⎣⎢|As(1)||As(2)|⋯|As(m)|⎤⎦⎥=A⎡⎣⎢⎢⎢⎢⎢⎢s(1)1s(1)2⋮s(1)ns(2)1s(2)2⋮s(2)n⋯⋯⋱⋯s(m)1s(m)2⋮s(m)n⎤⎦⎥⎥⎥⎥⎥⎥=As(38)

$$x=\left[ \begin{matrix}|&|&&|\\x^{(1)}&x^{(2)}&\cdots&x^{(m)}\\|&|&&| \end{matrix}\right] =\left[ \begin{matrix}|&|&&|\\As^{(1)}&As^{(2)}&\cdots&As^{(m)}\\|&|&&| \end{matrix}\right]=A\left[ \begin{matrix}s_1^{(1)}&s_1^{(2)}&\cdots &s_1^{(m)}\\s_2^{(1)}&s_2^{(2)}&\cdots&s_2^{(m)}\\\vdots&\vdots&\ddots&\vdots\\s_n^{(1)}&s_n^{(2)}&\cdots&s_n^{(m)} \end{matrix}\right]=As\tag{38}$$
$x$不是一个向量，是一个矩阵。其中每个列向量是$x^{(i)}$，$x^{(i)}=As^{(i)}$。$A$和$s$都是未知的，$x$是观测收集到的声音信号，是已知的，我们要想办法从$x$推出$s$。这个过程也称为盲信号分离。
令$W=A^{-1}$，有$s^{(i)}=A^{-1}x^{(i)}=Wx^{(i)}$，所以$W$可以表示为：

W=⎡⎣⎢⎢⎢⎢⎢———WT1WT2⋮WTn———⎤⎦⎥⎥⎥⎥⎥(39)

$$W=\left[ \begin{matrix}—&W_1^T&—\\—&W_2^T&—\\&\vdots&\\—&W_ n^T&—\end{matrix}\right] \tag{39}$$
其中$W_i\in R^n$,其实就是将$W_i$写成行向量形式。那么得到：$s^{(i)}_j=Wx^{(i)}_j$。

ICA不确定性（Ambiguities）

下面是两个原信号不确定性：
第一个，由于$w$和$s$都不确定，那么在没有先验知识的情况下，无法同时确定这两个相关参数。比如上面的公式$s=wx$。当$w$扩大两倍时，$s$只需要同时扩大两倍即可，等式仍然满足，因此无法得到唯一的$s$。
第二个，如果将$s$的顺序打乱，变成另外一个顺序，那么只需要调换$A$的列向量顺序即可，因此也无法单独确定$s$。
另外，还有一种ICA不适用的情况，那就是信号不能是高斯分布的。当源信号是高斯分布的时候，可以由不同的混合矩阵$A$乘上$s$，得到相同分布的$x$，一样无法确定源信号。

密度函数和线性变换

假设我们的随机变量$s$有概率密度函数$p_s(s)$（连续值是概率密度函数，离散值是概率）。为了简单，我们再假设$s$是实数，还有一个随机变量$x=As$，$A$和$x$都是实数。令$p_x(x)$是x的概率密度，那么怎么求$p_x(x)$？

求解之前先看一个密度函数与分布函数之间的关系。
设一个概率密度函数（Probability density function,简称pdf）为$f(x)$，一个累积分布函数（Cumulative distribution function，简称cdf）为$F(x)$，它们之间的关系为：

F(x)f(x)=∫xf(x)dx=F′(x)(40)

$$\begin{align}F(x)&=\int_xf(x)dx\\f(x)&=F'(x)\tag{40}\end{align}$$
以正态分布为例，它的概率密度函数与累积分布函数的图像如下图：

图八

一个很直观的感受就是，后者是由前者积分得到的，这也是符合式（40）的公式表示的。

关于概率密度的直接定义是：

Fx(a)=P(X≤a)=∫a−∞f(x)dx(41)

$$F_x(a)=P(X\leq a)=\int_{-\infty}^af(x)dx\tag{41}$$
关于$p_x(x)$的推导如下：

Fx(x)px(x)=P(X≤x)=P(As≤x)=P(s≤Wx)=Fs(Wx)=F′x(x)=F′s(Wx)=ps(Wx)|W|(42)

$$\begin{align}F_x(x)&=P(X\leq x)=P(As\leq x)=P(s\leq Wx)=F_s(Wx)\\p_x(x)&=F'_x(x)=F'_s(Wx)=p_s(Wx)|W|\tag{42}\end{align}$$

ICA算法

这里使用最大似然估计来解释算法， 我们假定每个$s_i$有概率密度$p_s$，那么给定时刻原信号的联合分布就是：

p(s)=∏i=1nps(si)(43)

$$p(s)=\prod_{i=1}^n p_s(s_i)\tag{43}$$
这个式子有一个假设前提：每个人发出的声音信号各自独立。由式（42），我们有：

p(x)=ps(Wx)|W|=|W|∏i=1nps(WTix)(44)

$$p(x)=p_s(Wx)|W|=|W|\prod_{i=1}^n p_s(W_i^Tx)\tag{44}$$
前面提到过，如果没有先验知识，我们无法求得W和s。因此我们需要知道$p_s(s_j)$，我们打算选取一个概率密度函数赋给s，但是我们不能选取高斯分布的密度函数。在概率论里我们知道密度函数$p(x)$由累计分布函数（cdf）$F(x)$求导得到。$F(x)$要满足两个性质是：单调递增和在[0,1]。我们发现sigmoid函数很适合，定义域负无穷到正无穷，值域0到1，缓慢递增。我们假定s的累积分布函数符合sigmoid函数：

g(s)=11+e−s(45)

$$g(s)=\frac{1}{1+e^{-s}}\tag{45}$$
求导后有：

ps(s)=g′(s)=es(1+es)2(46)

$$p_s(s)=g'(s)=\frac{e^{s}}{(1+e^{s})^2}\tag{46}$$
知道$p_s(s)$之后，就剩下参数$W$了，给定训练样本$\begin{Bmatrix} x^{(i)};i=1,2,\cdots,m \end{Bmatrix}$后，样本对数似然估计如下：

l(W)=log∏i=1mp(x(i))=∑i=1m⎛⎝∑j=1nlogg′(WTix(i))+log|W|⎞⎠(47)

$$l(W)=\log\prod_{i=1}^mp(x^{(i)})=\sum_{i=1}^m\left(\sum_{j=1}^n\log g'(W_i^Tx^{(i)})+\log|W|\right)\tag{47}$$
接下来就是对W求导了，这里牵涉一个问题是对行列式|W|进行求导的方法，属于矩阵微积分，这里直接给出结果：

∇W|W|=|W|(W−1)T(48)

$$\nabla_W|W|=|W|(W^{-1})^T\tag{48}$$
$W$最后的迭代公式如下，$\log g'(s)$的导数为$1-2g(s)$，$\alpha$是梯度上升速率，人为指定。：

W:=W+α⎛⎝⎜⎜⎜⎜⎜⎡⎣⎢⎢⎢⎢⎢1−2g(WT1x(i))1−2g(WT2x(i))⋮1−2g(WTnx(i))⎤⎦⎥⎥⎥⎥⎥x(i)T+(WT)−1⎞⎠⎟⎟⎟⎟⎟(49)

$$W:=W+\alpha\left(\left[ \begin{matrix}1-2g(W_1^Tx^{(i)})\\1-2g(W_2^Tx^{(i)})\\\vdots\\1-2g(W_n^Tx^{(i)})\end{matrix}\right] x^{(i)^T+(W^T)^{-1}}\right)\tag{49}$$

4、后记

从11月初到3月底学习这份课程，再从4月到现在5月初把这6份总结写完，前后也是用了半年时间——这都是在业余时间完成的，鬼知道我经历了什么。最后写到这里，竟然是不知道说什么了。清明三天早10晚11把第一份总结写完，在假期最后一天9点多躺在床上的时候，没有什么快感，只有一种身体和脑子被掏空的感觉，现在也是。而且一直压制着想要直接做项目的冲动，我知道一旦直接去操作了，我应该就没有办法写下这样的东西了。
想要学习决策树，Boosting，神经网络，深度学习，还有这个课程中缺失的最后一部分关于强化学习的完整内容，还有其他的很多很多东西，到时候会边学边写，一个个来，像这次积累了这么多东西然后一次写完真是太难过了。但是回头一看，这些东西也还是太少了。接下来主要是做点项目，没学到又要用到的算法和知识到时候先求会用。
时间是最大的敌人，世界在变好，弱是原罪。
告一段落。