2.3 模型优化概率方法

机器学习领域中有很多模型优化方法通过概率工具实现问题求解。从直观上看，算法应该追求尽可能地稳定，引入概率方法有时会破坏算法的稳定性。但机器学习领域的很多模型优化问题都需要处理巨大的求解空间或可能存在的不确定性信息，此时难以通过基本优化方法实现求解，概率方法则可以通过适当的概率工具有效地解决这些问题。因此，模型优化的概率方法是机器学习理论研究和应用开发领域非常重要的基础知识。本节主要介绍三种常用的模型优化概率方法，即随机梯度法、最大期望法和蒙特卡洛法。

2.3.1 随机梯度法

随机梯度法是对梯度下降法的一种改进。梯度下降方法在机器学习模型进行优化时，其搜索方向的每次更新都需要计算所有训练样本。例如，对于模型f（X；β），其中β为模型参数向量，即β={β₁，β₂，…，β_t}，若用训练样本集合S={（X₁，y₁），（X₂，y₂），…，（X_n，y_n）}对该模型进行训练并将经验风险作为优化目标，则优化的目标函数为

使用梯度下降法对上述目标函数进行优化时，第k次迭代的搜索方向计算公式为

其中，β_k为第k次迭代的起始点。

由于训练样本集S包含n个样本，故在计算P_k时需分别计算这n个样本的损失函数对参数向量β的梯度并求出它们的均值。当n较大时，计算梯度需要耗费大量时间。由于提高模型泛化性能需要尽可能多的训练样本，故不能通过减少训练样本的方式来简化计算。

事实上，当训练样本集S较大时，可用随机梯度法实现对目标函数的优化。随机梯度法的基本思想是随机选择少量训练样本来计算梯度，并将该梯度作为在全部训练样本上梯度的近似代替值用于梯度下降的迭代计算。随机梯度法有很多具体的实现算法，随机梯度下降法是其中最基本也是最常用的一种。下面具体介绍随机梯度下降法。

随机梯度下降法每次迭代的方向计算只与单个样本有关。对于机器学习模型f（X；β）及其训练样本集S={（X₁，y₁），（X₂，y₂），…，（X_n，y_n）}，该方法每次随机选取S中单个样本代替全部样本对模型参数进行一次更新，然后换另一个未参与训练的样本进行下次更新，当S中的全部样本都参与更新计算之后，随机调整S中所有样本排列次序后重复以上过程。

具体地说，假设随机梯度下降法在进行第k次迭代时随机选择的样本为（X_i，y_i），则此次参数更新的搜索方向可通过该样本的损失函数计算得到。令β_k为模型的当前参数向量，则用函数F′（β）=L（f_βk（X_i），y_i）实现此次迭代搜索方向的更新。需要注意的是，函数F′（β）为随机选择到的样本（X_i，y_i）的损失函数，仅用于计算参数搜索的更新方向，模型优化的目标函数则保持不变。方向更新的具体计算公式为

显然，使用上述公式仅需计算一次梯度便可实现方向更新。用代替P_k，并进行第k次参数更新的结果为

由以上分析可知：梯度下降法在n元训练样本集S上对模型进行优化时，每次参数更新需要执行n次梯度计算。使用随机梯度下降法只需执行一次梯度计算，可以有效减小计算量，使得模型优化过程在大样本场合变得可行。随机梯度下降法的具体过程如下：

（1）初始化t=0，k=0；

（2）对S中的全部训练样本进行随机排序，得到新的样本集合S^t；

（3）若S^t中的样本均参与过训练过程，则令t=t+1，并重复步骤（2），否则随机选择S^t中未参与训练的单个样本，并根据该样本数据计算参数更新方向P_k；

（4）根据式（2-36）更新参数向量得第k次参数更新的结果β_k+1；

（5）判断迭代是否满足算法终止条件，若满足则终止算法；否则令k=k+1，重复步骤（3）、（4）。

事实上，目前的计算机大多采用多核架构，具有一定的并行计算能力，若每次只使用单个样本确定参数搜索方向的更新，则会造成计算能力的浪费。因此，小批量随机梯度下降法也是一种常用的随机梯度法。

小批量随机梯度下降法的基本步骤与随机梯度下降法类似，基本思想是首先从S中随机选择一小批训练样本代替全部样本以实现对模型参数搜索方向的更新，然后在下次迭代搜索中再换另一小批不同样本进行搜索方向的更新计算，当训练样本集S中的全部样本都参与过一次参数的更新之后，将S中的全部样本进行随机排序后重新划分成不同小批次再去更新搜索方向，直至目标函数的取值接近最小值。显然，随机梯度下降法可以看作是小批量随机梯度下降法在每个小批量样本数均为1时的一种特殊情况。

现结合实例说明小批量随机梯度下降法的具体过程。假设训练集S共有m=100000个样本，将其中全部样本随机排序后得到集合S′，即有

首先，将S′划分为若干小批量训练样本集，例如，将S′等分为1000个规模为100的小批量训练样本集S₁，S₂，…，S₁₀₀₀，即有

第1个小批量训练集，…，，

第2个小批量训练集，…，，……

第1000个小批量训练集，…，。

然后，在每次迭代过程中随机选择一个当前未参与训练的小批次训练集进行搜索方向的更新计算。假设当前模型参数向量为β_k并选择第i个批次进行方向更新，则有

据此计算更新方向为

依据上述更新方向对当前参数向量β_k进行更新

当所有划分的批次均参与了训练过程而未达到算法终止条件时，对S′进行随机排序后重新划分成若干小批次样本集合，并重复上述搜索方向和参数更新过程，直至达到算法终止条件时结束算法。

随机梯度法每次仅选择单个或小批量样本来确定搜索方向和实现对模型参数的更新，可有效减少花费在梯度计算上的时间。该方法的计算结果存在一定的随机性，即不能保证每次更新的方向都是目标函数减小的方向，参数更新偶尔会使得目标函数取值增大。然而，取自相同总体的样本服从相同的概率分布，并且所有训练样本均参与迭代计算，故在参数更新的整体过程中，随机梯度下降法可以保证目标函数的取值不断下降。随机梯度法的参数更新过程如图2-6所示。

图2-6 随机梯度法参数更新示意图

另一方面，随机梯度方法的参数更新方向存在着一定的随机性，这也为避免陷入局部最优提供了可能。正常情况下梯度下降算法很难跳出局部最优，但是随机梯度法打破了参数更新的稳定性，这有可能会使迭代过程中出现跳出当前局部最小值范围的情况。

目前，很多机器学习任务的训练样本集规模都很大，由以上分析可知随机梯度方法能够在保证优化效果的前提下有效减少在大规模训练样本集上对模型进行优化的计算量。因此，随机梯度法是目前最常用的模型优化方法之一。

2.3.2 最大期望法

有些机器学习模型含有隐含变量，通常难以对这些具有隐含变量的模型直接进行参数求解或估计，此时可通过迭代求解隐含变量（或其函数）数学期望最大值的方法实现对带隐含变量模型的参数优化求解。这类优化方法通常称为最大期望法或EM算法，主要包括两个计算步骤，即计算数学期望E的步骤和计算函数最大值优化M的步骤。可用EM算法解决带隐含变量模型参数的最大似然估计或者最大后验估计问题。这里主要讨论如何通过EM算法求解参数的最大似然估计，最大后验估计的EM求解方法与此类似。

假设某学校男生和女生身高分别服从两个参数不同的高斯（正态）分布，现从该学校所有学生当中随机抽取100名学生，分别测量他们的身高数据但不记录性别信息。根据这些数据对男生身高和女生身高所服从正态分布进行参数估计时，则存在一个隐含的性别数据。由于不知道性别信息，无法直接得知性别数据的分布，故无法直接求得已知样本数据所满足的似然函数。然而，在给定参数情况下，可以结合已知样本计算出某个同学性别的概率分布情况，并且可以知道在给定参数情况下性别与身高的联合概率分布。

更一般地，在一个包含隐含数据Z的模型参数估计问题中，X为可直接观测数据，即已知样本，β为模型参数向量。由于隐含数据Z的存在，通常无法直接得知已知样本X取值下的似然函数L（β|X）=p（X|β），但可知参数给定情况下X和Z的联合概率分布p（X，Z|β），以及在参数向量和可观测数据给定情况下Z取值状态的条件概率分布p（Z|X，β）。

现根据以上信息对模型参数向量β进行最大似然估计。最大似然估计通过最大化似然函数实现，然而已知样本X取值状态的似然函数L（β|X）=p（X|β）却由于隐含数据Z的存在而难以直接求得，故难以直接使用似然函数对模型参数向量β进行最大似然估计。由于已知参数β给定条件下X和Z的联合概率分布为p（X，Z|β），故将上述似然函数转化为

考虑其对数似然，有

根据最大似然估计思想，只需求得对数似然lnL（β|X）的最大值即可得到模型参数β的最大似然估计值。但是由于上式存在隐含数据Z和积分的对数，直接求解其最大值较为困难，故用迭代逼近方法实现对数似然最大值的估计。为此，将对数似然做如下变形

其中，p（Z）为隐含数据Z的某一分布。

由于对数函数为凸函数，故如下不等式成立

由此可得对数似然的下界函数

取p（Z）=p（Z|X，β_t），则可得

略去下界函数B（β，β_t）中与待求参数向量β无关的项，得到如下Q函数

上式表示对隐含变量Z的函数L（β|X，Z）在概率分布p（Z|X，β_t）下的数学期望。由于B（β，β_t）≤lnL（β|X），故可通过迭代选取不同下界函数B（β，β_t）最大值的方式逐步逼近对数似然lnL（β|X）的最大值，迭代逼近的具体过程如图2-7所示。

图2-7 EM算法逼近似然函数最大值

由于Q（β，β_t）是通过省略B（β，β_t）中与待求参数向量β无关项得到的，用Q（β，β_t）迭代求解对数似然最大值与用B（β，β_t）求解对数似然最大值具有相同的效果，故EM算法通常直接使用Q函数进行优化计算。由以上分析可得EM算法的基本步骤如下。

（1）设置初始参数β₀和迭代停止条件；

（2）E步（期望步）：根据可直接观测数据X和当前参数向量取值β_t计算Q（β，β_t）：

Q（β，β_t）=∫lnL（β|Z）p（Z|X，β_t）dZ

（3）M步（最大化步）：最大化Q（β，β_t）并根据Q（β，β_t）的最大值更新参数β_t的取值：

（4）判断是否满足迭代停止条件，若满足则停止迭代；否则令t=t+1并返回步骤（2）。

EM算法可随机选择初始参数β₀，但应注意EM算法对初始参数β₀具有一定的敏感性。如果初值β₀选取不当，则可能会使迭代结果陷入局部最优。

通常将迭代停止条件设为相邻两次参数值变化不大或前后两次参数更新使得Q（β，β_t）值变化很小时停止更新，即有

|β_t+1-β_t|<ε或|Q（β_t+1，β_t）-Q（β_t，β_t）|<ε

【例题2.9】假设X₁，X₂取自指数分布Y（θ）=θe^-xθ，且X₁，X₂不相关。若X₁=10而X₂缺失，试用EM算法实现参数θ的最大似然估计。

【解】由于X₁，X₂为取自指数分布Y（θ）=θe^-xθ的离散值，则有

由于X₁，X₂同分布，则有，从而有

为求出使得Q（θ，θ_t）最大化的参数θ，可将使得Q（θ，θ_t）最大化的参数作为第t+1次的估计值θ_t+1，即得迭代公式

当迭代收敛时，有θ_*=2θ-/（10θ_*+1），解得θ_*=0.1，即迭代算法求得参数估计值收敛于0.1，即参数θ的最大似然估计值为θ_*=0.1。□

2.3.3 蒙特卡洛法

蒙特卡洛法是一类以概率统计理论为基础的随机型数值算法。一般来说，当一个随机算法满足采样次数越多，其输出结果越近似最优解这一特性时，便可称为蒙特卡洛法。该方法通常将待求解问题与某一概率模型联系起来，利用从大量样本中获得的信息来完成对所求参数的概率估计，由此实现对实际问题的求解。

例如，可用蒙特卡洛法近似计算一个不规则湖面的面积。假设围住湖面的长方形面积为M，湖面面积S为未知数。由于湖面形状不规则无法直接求出S，为此向包含湖面的长方形区域内随机撒布n个点，假设其中有k个点落在湖面当中，则可得到湖面大小S的估计值

显然，撒布点数n越大，估计值就越接近于湖面的真实面积S。这是因为根据大数定律和中心极限定理，重复进行大量试验时，事件A出现的频率接近事件A发生的概率。

蒙特卡洛法最早用于近似求解难以精确计算的定积分，对于某函数φ（x）的积分

若能找到定义在（a，b）上的一个函数f（x）和概率密度函数p（x），满足φ（x）=f（x）p（x），则可将上述积分I转化为

如此一来，就将原积分I的求解问题转化为了求解函数f（x）在p（x）分布上的数学期望问题。假设在分布p（x）上做随机采样得到的样本集合为{x₁，x₂，…，x_n}，则可用这些样本来对I进行估计，即有

由于

故I_p是I的无偏估计，即对于任意分布p（x），积分I的蒙特卡洛法估计值是无偏的。

有时在分布p（x）上进行采样会出现一些问题，例如采样偏差较大等，甚至在某些时候根本无法对其进行采样。此时，可通过变换积分的分解形式找到另外一个易于采样的分布q（x），将积分I转化为如下形式

其中，q（x）为某一分布，可将f（x）[p（x）/q（x）]看作是某一函数。

令g（x）=f（x）p（x）/q（x），则有

此时可将积分I看作是函数g（x）在分布q（x）上的期望。假设在分布q（x）上进行随机采样获得的样本集合为{，，…，}，则可用该样本集合对I进行估计，即有

在新分布q（x）上进行采样的过程被称为重要性采样，p（x）/q（x）被称为重要性权重。由于积分I任意分解形式通过蒙特卡洛法得到的估计值都是无偏的，故应适当选择新分布q（x），使得估计值I_q的方差尽可能地达到最小。对于I_q的方差Var（I_q）

由于其后项与分布q（x）无关，故只需最小化期望E_q[f²（x）p²（x）/q²（x）]。

根据琴生不等式，有

故取q（x）为如下分布时可使得方差Var（I_q）达到最小

重要性采样可以改进原采样，如图2-8a所示，若在分布p（x）上进行采样，则很可能仅仅采样到函数f（x）的极端值；如图2-8b所示，若在分布q（x）上进行重要性采样，则函数f（x）的取值较为丰富。

在实际应用当中，重要性采样并非总能奏效。在很多情况下，重要性采样q（x）难以取得理论上的最佳分布q^*（x），而只能取到一个方差较大的可行分布，难以满足任务需求。此时用马尔可夫链蒙特卡洛法（Markov Chain Monte Carlo，MCMC）近似计算待求参数概率分布的方式解决。

假设离散型随机变量X的取值范围是{X₁，X₂，…，X_n}，则称该集合为X的状态空间。马尔可夫链是一个关于离散型随机变量取值状态的数列，从X的随机初始状态X_i开始，马尔可夫链依据一个只与前一时序状态相关的状态转移分布P（X^t+1|X^t）来确定下一时序的状态，其中X^t和X^t+1分别表示随机变量X在第t时序和第t+1时序的状态。

图2-8 不同采样方式对比示意图

a）在分布p（x）上进行采样 b）在分布q（x）上进行采样

MCMC法同时运行多条马尔可夫链，这些马尔可夫链在同一时序内具有相同的状态空间，即服从同一个状态概率分布。随着马尔可夫链的不断更新，其状态概率分布最终会收敛于某一分布。具体来说，MCMC法从一个初始概率分布p⁰中随机选择多条马尔可夫链的初始状态，p⁰的具体形式为

其中，表示在时序为0的初始状态分布下X=X_i的概率。

由于这些马尔可夫链的状态空间是一个n维空间，故可将概率分布p⁰表示为一个n维向量p⁰=（，，…，）^T。此时，马尔可夫链下一时序的分布可表示为

上式说明马尔可夫链在第t+1时序的状态分布只与前一时序的状态分布p^t和状态转移分布P（X^t+1|X^t）相关。

接着考虑所有同时运行的马尔可夫链状态概率分布情况。若用矩阵A表示P（X^t+1|X^t），矩阵A中第i行第j列元素a_ij表示从状态X_j转移到状态X_i的概率P（X_i|X_j），则这些马尔可夫链在同一时序上状态概率分布的变化可表示为

状态概率分布会随着时序t不断变大而最终收敛，即有p^t+1=p^t。此时可用所求状态概率分布作为真实状态概率分布的一种近似分布。

事实上，由于贝叶斯学派认为模型的未知参数服从于某一概率分布，并将参数的概率分布看作一种状态概率分布，故可通过MCMC法求解其近似分布的方式实现模型优化。

MCMC法的基本思想如图2-8所示。对于参数的概率分布，MCMC法从随机初始点开始进行采样，随着采样的不断进行，其样本点的概率分布会逐步收敛。图中圆点表示概率分布还未收敛时采样到的点，叉点表示概率分布收敛后采样到的点。由于最终收敛的概率分布是对参数概率分布的一种近似，故在其中进行大量采样便可求出该近似概率分布，即求得参数概率分布的一种近似分布。

图2-9 MCMC法的基本思想

显然，MCMC法的关键在于如何确定状态转移概率分布P（X^t+1|X^t），如果P（X^t+1|X^t）的选择合适，则状态概率分布会收敛到待求参数的概率分布。从收敛性角度考虑，MCMC法的状态转移概率分布P（X^t+1|X^t）应保证多条马尔可夫链同时满足状态分布收敛的条件。

事实上，当状态概率分布p（X）与状态转移概率分布P（X^t+1|X^t）满足细致平稳条件时，状态转移分布P（X^t+1|X^t）就可以使状态概率分布收敛，此时的状态分布p（X）称为状态转移分布P（X^t+1|X^t）的平稳分布。

所谓细致平稳条件，是指在状态分布p（X）中，利用状态转移分布对其状态进行更新时，由X_i转移到X_j的概率应与由X_j转移到X_i的概率相同，即有p（X_i）P（X_j|X_i）=p（X_j）P（X_i|X_j）。

在实际应用中找到一个状态转移分布并使其满足细致平稳条件有时是一件比较困难的事情，通常需要使用某些采样方法构造一个满足条件的状态转移分布。其中最简单的采样方法为MCMC采样。具体地说，对于状态分布p（X），可随机选择一个状态转移分布Q（X^t+1|X^t），它们通常不满足细致平稳条件，即

但可选择一对转移算子α（X_i，X_j）和α（X_j，X_i）满足下列等式

可用转移算子α（X_i，X_j）和α（X_j，X_i）帮助确定是否接受采样结果。要想使得上述等式成立，最简单的方法是令

并取状态转移分布为P（X_i|X_j）=Q（X_j|X_i）α（X_i，X_j），则有状态分布p（X）即为P（X_i|X_j）的平稳分布，此时就可通过MCMC方法获得一个收敛的状态概率分布。

MCMC采样的基本步骤如下：

（1）随机选定一个状态转移分布Q（X_j|X_i）并记状态分布为p（X），设定收敛前状态转移次数上限τ和采样个数n；

（2）从任意分布中采样获得初始状态X₀；

（3）依据条件状态分布Q（X_j|X_i）进行采样，即采样值为X^t+1=Q（X|X^t）；

（4）从均匀分布U[0，1]中采样并记为u_t，若u_t<α（X^t，X^t+1）=p（X^t）Q（X^t+1|X^t），则接受X^t+1，并令t=t+1，否则，拒绝此次转移并保持t值不变；

（5）当t>τ+n时终止算法。此时共取到τ+n个样本，但前τ个样本在分布未收敛时获得，不能作为MCMC采样结果，故只将第τ+1到τ+n个样本作为采样结果。

除了MCMC采样之外，还有MH采样、Gibbs采样等获取样本方法。其中MH采样是对MCMC采样的改进，Gibbs采样是对MH采样的改进，这里不再赘述。