从生成模型到GDA再到GMM和EM算法

在学习生成模型之前，先学习了解下密度估计和高斯混合模型。为什么呢？因为后面的VAE\GANs模型都需要把训练样本，也就是输入的图像样本看作是一个复杂的、多维的分布。

1. 知乎上关于图像频率的解释

作者：耳东陈
链接：https://www.zhihu.com/question/20099543/answer/13971906
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

首先说说图像频率的物理意义。图像可以看做是一个定义为二维平面上的信号，该信号的幅值对应于像素的灰度（对于彩色图像则是RGB三个分量），如果我们仅仅考虑图像上某一行像素，则可以将之视为一个定义在一维空间上信号，这个信号在形式上与传统的信号处理领域的时变信号是相似的。不过是一个是定义在空间域上的，而另一个是定义在时间域上的。所以图像的频率又称为空间频率，它反映了图像的像素灰度在空间中变化的情况。例如，一面墙壁的图像，由于灰度值分布平坦，其低频成分就较强，而高频成分较弱；而对于国际象棋棋盘或者沟壑纵横的卫星图片这类具有快速空间变化的图像来说，其高频成分会相对较强，低频则较弱（注意，是相对而言）。再来谈一谈如何定量的测量图像的空间频率，最为常用的方法就是二维傅里叶变换。图像经过二维傅里叶变换后会形成与图像等大的复数矩阵，取其幅值形成幅度谱，取其相位形成相位谱。图像的频率能量分布主要体现在幅度谱中。通常习惯将低频成分放在幅度谱的中央，而将高频成分放在幅度谱边缘。大多数自然图像的幅度谱在统计上呈现1/f^2分布，也就是频率成分的能量与频率的平方成反比。所以从绝对数值上看，低频能量通常是要高于高频能量的，这一规则也称为power law。power law并非是上帝的无心之作，事实上power law的出现时源于自然图像的尺度不变性（scale invariance）。这一点在很多文献中被解释为从不同的距离观察同样的自然场景，获得的图像的幅度谱是基本相同的。相关内容可以搜索关键字power law & natural image statistics。除了傅里叶变换外，正弦变换、余弦变换、Gabor变换、小波变换、WH变换也可以用来对图像频率分布进行定量测量。目前小波变换是研究的热点，因为小波变换不但能够反映频率能量的分布，同时还保留了图像特征的空间分布特性。

2. 生成学习算法(cs229,Ng)

2.1 生成算法和判别算法的区别

举个栗子：

我们要区分elephants(y=1)和dogs(y=0)

对判别模型（discriminative），以logistic回归为例：
- logistic回归模型：p(y|x;θ),hθ=g(θTx)其中g是sigmoid函数。通过logistic回归，我们找到一条决策边界decision boundary，能够区分elephants和dogs.
而对于生成模型（generative）：
- 模型：p(x|y)，在给定了样本所属的类的条件下，对样本特征建立概率模型。
- p(x|y=1)是elephants的分类特征模型
- p(x|y=0)是dogs的分类特征模型
已知p(y)是先验类别的分布，生成模型p(x|y)可通过学习得到，则根据Bayes公式，可以推导出后验概率分布,即通过特征数据进行预测
$p (y = 1 | x) = p (x | y = 1) p (x) p (x)$
其中分母可以这么计算得到p(x)=p(y=0|x)p(x)+p(y=1|x)p(x)，
但p(y)和p(x|y)我们已经知道了，那么可以不用计算p(x),计算方法如下：
$a r g max y p (y | x) = a r g max y p (x | y) p (y) p (x) = a r g max y p (x | y) p (y)$

总结下就是：
- 生成模型：一般是学习一个代表目标的模型，然后通过它去搜索图像区域，然后最小化重构误差。类似于生成模型描述一个目标，然后就是模式匹配了，在图像中找到和这个模型最匹配的区域，就是目标了。

判别模型：以分类问题为例，然后找到目标和背景的决策边界。它不管目标是怎么描述的，那只要知道目标和背景的差别在哪，然后你给一个图像，它看它处于边界的那一边，就归为哪一类。
由生成模型可以得到判别模型，但由判别模型得不到生成模型。

然鹅，生成模型p(x|y)怎么得到呢？不慌，我们先了解下多维正态分布～

2.2 多维正态分布(the multivariate nirmal distribution)

从生成模型到GDA再到GMM和EM算法

关于一维正态分布怎么推导出多维正态分布的概率密度函数，可参考知乎:多维高斯分布是如何由一维发展而来的？

首先一维正态分布:

p(x)=12π‾‾‾√exp(−x22)

二维标准正态分布，就是两个独立的一维标准正态分布随机变量的联合分布：

p(x,y)=p(x)p(y)=12πexp(−x2+y22)

把两个随机变量组合成一个随机向量：v=[xy]T

p(v)=12πexp(−12vTv) 显然x,y相互独立的话，就是上面的二维标准正态分布公式～

然后从标准正态分布推广到一般正态分布，通过一个线性变化：v=A(x−μ)

p(x)=|A|2πexp[−12(x−μ)TATA(x−μ)]

注意前面的系数多了一个|A|（A的行列式）。

可以证明这个分布的均值为μ，协方差为(ATA)−1。记Σ=(ATA)−1，那就有

p (x) = 1 2 π | Σ | 1 / 2 exp [- 12 (x - μ) T Σ - 1 (x - μ)]

推广到n维：

p (x) = 1 (2 π) n / 2 | Σ | 1 / 2 exp [- 12 (x - μ) T Σ - 1 (x - μ)]

需要注意的是：这里的二维、n维到底指的是什么？
- 以飞机检测的数据点为例，假设它由heat和time决定，那么这就是个二维正态分布，数据点的生成所处的位置由其概率决定，也就是p(x)
- 如果这个数据有n个特征，那么其分布就是n维正态分布。
- 之前一直理解的是，n维正态分布是两个向量巴拉巴拉。。好像一直没搞懂。。

再顺便了解下协方差矩阵吧～

2.3 关于协方差矩阵，参考blog

对多维随机变量X=[X1,X2,…,Xn]T，我们往往需要计算各维度之间的协方差，这样协方差就组成了一个n×n的矩阵，称为协方差矩阵。协方差矩阵是一个对角矩阵，对角线上的元素是各维度上随机变量的方差,非对角线元素是维度之间的协方差。我们定义协方差为Σ, 矩阵内的元素Σij为:

Σ i j = c o v (X i, X j) = E [(X i - E (X i)) (X j - E (X j))]

则协方差矩阵为:

Σ = E [(X - E (X)) (X - E (X)) T]

= ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ c o v (X 1, X 1) c o v (X 2, X 1) ⋮ c o v (X n, X 1) c o v (X 1, X 2) c o v (X 2, X 2) ⋮ c o v (X n, X 2,) \dots \dots ⋮ \dots c o v (X 1, X n) c o v (X 2, X n) ⋮ c o v (X n, X n) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

如果X~N(μ,Σ),则Cov(X)=Σ

了解了多维正态分布和协方差，我们再回到生成模型p(x|y)。。其实我们就是假设对于n维特征，p(x|y)是n维正态分布～怎么理解呢，下面就说！

2.4 高斯判别分析模型The Gaussian Discriminant Analysis model

高斯判别模型就是：假设p(x|y)是一个多维正态分布

则生成模型：

y ~ Bernoulli(ϕ) 伯努利分布，又称两点分布，0-1分布

x|y=0 ~ N(u0,Σ)

x|y=1 ~ N(u1,Σ)

这里可以看作是一个二分类，y=0和y=1,可以看作是伯努利分布，则p(y)=ϕy(1−ϕ)1−y，要学的参数之一: ϕ=p(y=1)，试想如果是多分类呢，那么要学习的参数就有ϕ1,ϕ2,....ϕk

其中类别对应的特征x|y=0,x|y=1服从正态分布。怎么理解呢？就是既然你们都是一类人，那么你们的身高啊，体重啊等等应该满足正态分布。。有几维特征就满足几维正态分布

这里x是n维特征，身高，体重，颜值…balabala，所以x|y=0满足n维正态分布～x|y=1也是啦，只不过对于不同的类，对应n维特征的均值不一样，奇怪为什么协方差矩阵是一样的？？？？？？？？？？

每个分类对应的n维特征的分布显然不是独立的，比如体重和颜值还是有关系的吧～他们的协方差，方差就统统都在Σ协方差矩阵里面了

p (x | y = 0) = 1 (2 π) n / 2 | Σ | 1 / 2 exp [- 12 (x - μ 0) T Σ - 1 (x - μ 0)]

p (x | y = 0) = 1 (2 π) n / 2 | Σ | 1 / 2 exp [- 12 (x - μ 1) T Σ - 1 (x - μ 1)]

这样，模型中我们要学习的参数有ϕ,Σ,μ0,μ1，对于训练数据，就是观测到的数据x,y，既然他们出现了，那么他们的联合概率，也就是似然函数∏mi=1p(x,y)就要最大～其对数似然log-likelihood：

\begin{equation}\begin{aligned}L(\phi,\Sigma, \mu_0,\mu_1) &= log\prod_{i=1}^mp(x^{(i)},y^{(i)};\phi,\Sigma, \mu_0,\mu_1$) \\  &= log\prod_{i=1}^mp(x^{(i)}|y^{(i)};\phi,\Sigma, \mu_0,\mu_1$) p(y^{(i)};\phi)\\\end{aligned}\end{equation}\label{eq2}

其中p(y(i);ϕ)是已知的，也就是先验概率(class priors)，p(x(i)|y(i))就是上面推导的～代入后，分别对参数求导即可：

从生成模型到GDA再到GMM和EM算法

3. 高斯混合模型GMM

其实和GDA非常相似，不过这里没有了类标签，只有一堆样本特征，{x(1),x(2),...,x(m)},
我们不知道这些样本属于几个类别，也不知道有哪些类了。但虽然不知道，我们确定他们是存在的，只是看不见而已。我们可以假设存在k类，{z(1),z(2),...,z(k)},看不见的，我们就叫它们隐藏随机变量(latent random variable)，

这样一来，就训练样本就可以用这样的联合分概率模型表示了，p(x(i),z(i))=p(x(i)|z(i))p(z(i))

同GDA不一样的是，这里是多分类，可假定z(i)∼Multinomial(ϕ)，多项式分布（二项分布的拓展～），那么p(z(i))=ϕj

同GDA相同的是，对于每一个类别，其对应的样本满足n维正态分布，也就是：x(i)|z(i)=j∼N(μj,Σj),但注意哦，这里每个高斯分布使用了不同的协方差矩阵Σj

p (x | z (1)) = 1 (2 π) n / 2 | Σ 0 | 1 / 2 exp [- 12 (x - μ 0) T Σ - 1 0 (x - μ 0)]

p (x | z (2)) = 1 (2 π) n / 2 | Σ 1 | 1 / 2 exp [- 12 (x - μ 1) T Σ - 1 1 (x - μ 1)]

. . . .

p (x | z (k)) = 1 (2 π) n / 2 | Σ k | 1 / 2 exp [- 12 (x - μ k) T Σ - 1 k (x - μ k)]

然后带入到训练样本的对数似然（log-likelihood）:

\begin{equation}\begin{aligned}L(\phi,\mu,\Sigma)&=\sum_{i=1}^{m}logp(x^{(i)};\phi,\mu,\Sigma) \\  &= \sum_{i=1}^{m}log\sum_{z^{(i)}=1}^{k}p(x^{(i)}|z^{(i)};\mu,\Sigma) p(z^{(i)};\phi)\\\end{aligned}\end{equation}\label{eq2}

这里需要注意下标：对于类别有k类，第一个求和符号是对第i个样本在k个类别上的联合概率，第二个求和符号是m个样本的联合概率。

我们可以注意到，如果我们知道z(i),那么这个似然函数求极大值就很容易了，类似于高斯判别分析，这里的z(i)相当于标签，分别对参数求导可得：
GMM2.png

解释一下参数:1{z(i)=j}表示第i个样本为j类时，这个值就为１，那么ϕj=1m∑mi=11{z(i)=j}表示样本中类别为j的概率

其中p(z(i);ϕ)是根据伯努利分布得到的，跟GDA不同的是，p(y|ϕ)是已知的频率概率。So…z(i)到底有多少个分类？每个类别的概率是多少？

譬如上式中∑mi=11{z(i)=j}这个没法求对吧～

specific bootstrap procedure

这个时候EM算法就登场了～～～

3.1 EM算法

上面也提到了，如果z(i)是已知的话，那么ϕj=1m∑mi=11{z(i)=j}表示类别j的概率p(z(i)=j)也就已知了，但是呢？我们不知道。。所以我们要猜测p(z(i)=j)这个值，也就是EM算法的第一步：

Repeat until convergence:{

(E-step):for each i,j,set:

$w (i) j : = p (z (i) = j | x (i); ϕ, μ, Σ)$

w(i)j什么意思呢?就是对于i样本,它是j类的后验概率。在GDA里面，x_i的类别是确定的，在GMM里面呢？不知道它的类别，所以只能假设k类都有可能，它是j类别的概率就是w(i)j，在GDA里面，它仅仅取决于ϕj,而在GMM里面，它取决于ϕj,μj,Σj

The values wj calculated in the E-step represent our “soft” guesses for
the values of z(i) .

The term “soft” refers to our guesses being probabilities and taking values in [0, 1]; in
contrast, a “hard” guess is one that represents a single best guess (such as taking values
in {0, 1} or {1, … , k}).

这样一来，参数更新就可以这样写了，也就是EM算法的第二步：

(M-step) Updata the parameters:

这一步怎么推导的，后面会详细介绍

$ϕ j : = 1 m \sum i = 1 m w (i) j$

$μ j : = \sum m i = 1 w (i) j x (i) \sum m i = 1 w (i) j$

$Σ j : = \sum m i = 1 w (i) j (x (i) - μ j) (x (i) - μ j) T \sum m i = 1 w (i) j$

｝

w(i)j表示第i个样本为ｊ类别的概率，而ϕj表示m个样本中j类别的概率，μj,Σj分别表示j类别对应的n维高斯分布的期望和协方差矩阵

所以，求出w(i)j，一切就都解决了吧？对于后验概率p(z(i)=j|x(i))可以根据Bayes公式：

p (z (i) = j | x (i); ϕ, μ, Σ) = p (x (i) | z (i) = j; μ, Σ) p (z (i) = j; ϕ) \sum k l = 1 p (x (i) | z (i) = l; μ, Σ) p (z (i) = l; ϕ)

其中先验概率p(x(i)|z(i)=j;μ,Σ)可以根据高斯分布的密度函数来求，z(i)=j∼N(μj,Σj)
类先验(class priors)p(z(i)=j;ϕ)可以取决于多项式分布中j类的概率ϕj

The EM-algorithm is also reminiscent of the K-means clustering algorithm, except that instead of the “hard” cluster assignments c(i), we instead have the “soft” assignments wj . Similar to K-means, it is also susceptible to local optima, so reinitializing at several different initial parameters may be a good idea.

EM算法使我们联想起了k-means,区别在于k-means的聚类是通过欧氏距离c(i)来定义的，而EM是通过wjprobabilities来分类的。同k-means一样，这里的EM算法也是局部优化，因此最好采用不同的方式初始化～

3.2 convergence?

我们知道k-means一定是收敛的，虽然结果不一定是全局最优解，但它总能达到一个最优解。但是EM算法呢，一定收敛吗？

4. The EM algorithm

前面我们讲的是基于高斯混合模型的EM算法，但一定所有的类别都是高斯分布吗？还有卡方分布，泊松分布等等呢，接下来我们就将讨论EM算法的一般性。

在学习一般性的EM算法前，先了解下Jensen’s inequality

4.2 Jensen’s inequality

如果函数f，其二阶导恒大与等于0 (f″≥0)，则它是凸函数f(convec function)。
如果凸函数的输入是向量vector-valued inputs，那么它的海森矩阵(hessian)H是半正定的。Jensen’s 不等式：

Let f be a convex function, and let X be a random variable.
Then:
$E [f (X)] \geq f (E X) .$
Moreover, if f is strictly convex, then E[f(X)]=f(EX) holds true if and
only if X=E[X] with probability 1 (i.e., if X is a constant).

举个栗子来解释jensen不等式：
从生成模型到GDA再到GMM和EM算法

假设输入随机变量X是一维的哈，然后Ｘ取a,b的概率都是0.5,那么

E X = (a + b) / 2, f (E X) = f (a + b 2)

E [f (X)] = f (a) + f (b) 2

因为是凸函数，所以 f(EX)≤E[f(X)]

同理，如果是凹函数(concave function),那么不等式方向相反f(EX)≥E[f(X)]。后面EM算法里面就要用到log(X)，log(x)就是个典型的凹函数～

4.2 The EM algorithm

首先，问题是：我们要基于给定的m个训练样本{x(1),x(2),...,x(m)}来进行密度估计～

像前面一样，创建一个参数模型p(x,z)来最大化训练样本的对数似然：

\begin{equation}\begin{aligned}L(\theta)&=\sum_{i=1}^mlogp(x;\theta)\\&=\sum_{i=1}^mlog\sum_zp(x,z;\theta)\\\end{aligned}\end{equation}\label{eq2}

可以看到，z(i)是隐藏的随机变量(latent random variable),关于参数\theta的最大似然估计就很难计算了。

其实同GMM一样，我们也可以把z(i)看做是类别，每个类别其分布函数是Qi，在GMM里面，Qi是多项式分布概率z(i)∼multinomial(ϕ)，但一般性的EM算法，我们不知道它是啥，只知道∑kzQi(z)=1,即z(i)∼Qi

解释下公式中的推导：
- 这里是针对样本i来说，对于样本i，它可能是z1,z2,...,zk都有可能，但他们的probability之和为１，也就是∑zQi(z)=1

(2)到(3)的推导：可以将[p(x(i),z(i);θ)Qi(z(i))]看做随机变量Ｘ取z^{(i)})},随机变量Ｘ有k个取值，那么（２）式中的后半部分 log∑z(i))[p(x(i),z(i);θ)Qi(z(i))]就是log(EX)了，logx是一个凹函数，则其大于E[log(x)]
log(θ)无法直接求导找到其最大值，但是可以先找到其下界凹函数，然后取凹函数的最大值，一次迭代下去。。。下界凹函数就是(3)式

看下Ng画的图：
从生成模型到GDA再到GMM和EM算法

怎么确定下界函数在哪里呢？找到下界函数和l(θ)交接的点，这就是EM算法的第一步：
我们要让不等式相等,即Jensen’s inequality中的随机变量取值是一个常量，看(2)式：

p (x (i), z (i); θ) Q i (z (i)) = c

根据∑zQi(z)=1,联立上式求解可得：

\begin{equation}\begin{aligned}Q_i(z^{(i)})&=\dfrac{p(x^{(i)},z^{(i)};\theta)}{\sum_zp(x^{(i)},z^{(i)};\theta)}\\&=\dfrac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)}\\&=p(z^{(i)}|x^{(i)};\theta)\\\end{aligned}\end{equation}\label{eq2}

对比GMM来看会发现：
- 在高斯混合模型中，我们做了两个假设，一个是隐藏随机变量z(i)有k类，且满足ｋ多项式分布z(i)∼multinomial(ϕi)，另一个就是对于每一个分类，其对应的n维样本特征满足n维正态分布z(j)i∼N(μi,Σi)。
- 而在EM算法中，将他们一般化了，也就是泛化了。

EM总结下来：

Repeat until convergence {

(E-step)
For each i,找到下界函数, set:
$Q i (z (i)) : = p (z (i) | x (i); θ)$

(M-step)找到下界凹函数的最大值,也就是(3)式 Set:
$θ : = a r g max θ \sum i m \sum z (i) k Q i (z (i)) l o g p (x (i), z (i); θ) Q i (z (i))$
}

4.3 training and will converge?

首先说是否收敛，答案是肯定收敛的。。懒得输公式了。。直接贴图吧，这个比较好理解：
这里写图片描述

理论推导模模糊糊似懂非懂，感觉是那么回事儿。但是真让去推导或是代码实现，感觉又一团糟了。所以，EM算法是迭代算法，每一步迭代到底做了些啥？参数有多少？到底是怎么更新参数的？

我们知道了GMM是EM算法的一种特殊形式，那么就以高斯混合模型为例，进行训练：

4.4 Mixture of Gaussians revisited

从生成模型到GDA再到GMM和EM算法