机器学习领域中有很多模型优化方法通过概率工具实现问题求解。从直观上看,算法应该追求尽可能地稳定,引入概率方法有时会破坏算法的稳定性。但机器学习领域的很多模型优化问题都需要处理巨大的求解空间或可能存在的不确定性信息,此时难以通过基本优化方法实现求解,概率方法则可以通过适当的概率工具有效地解决这些问题。因此,模型优化的概率方法是机器学习理论研究和应用开发领域非常重要的基础知识。本节主要介绍三种常用的模型优化概率方法,即随机梯度法、最大期望法和蒙特卡洛法。
随机梯度法是对梯度下降法的一种改进。梯度下降方法在机器学习模型进行优化时,其搜索方向的每次更新都需要计算所有训练样本。例如,对于模型 f ( X ; β ),其中 β 为模型参数向量,即 β ={ β 1 , β 2 ,…, β t },若用训练样本集合 S ={( X 1 , y 1 ),( X 2 , y 2 ),…,( X n , y n )}对该模型进行训练并将经验风险作为优化目标,则优化的目标函数为
使用梯度下降法对上述目标函数进行优化时,第 k 次迭代的搜索方向计算公式为
其中, β k 为第 k 次迭代的起始点。
由于训练样本集 S 包含 n 个样本,故在计算 P k 时需分别计算这 n 个样本的损失函数对参数向量 β 的梯度并求出它们的均值。当 n 较大时,计算梯度需要耗费大量时间。由于提高模型泛化性能需要尽可能多的训练样本,故不能通过减少训练样本的方式来简化计算。
事实上,当训练样本集 S 较大时,可用随机梯度法实现对目标函数的优化。随机梯度法的基本思想是随机选择少量训练样本来计算梯度,并将该梯度作为在全部训练样本上梯度的近似代替值用于梯度下降的迭代计算。随机梯度法有很多具体的实现算法,随机梯度下降法是其中最基本也是最常用的一种。下面具体介绍随机梯度下降法。
随机梯度下降法每次迭代的方向计算只与单个样本有关。对于机器学习模型 f ( X ; β )及其训练样本集 S ={( X 1 , y 1 ),( X 2 , y 2 ),…,( 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 1 , S 2 ,…, S 1000 ,即有
第1个小批量训练集
,…,
,
第2个小批量训练集
,…,
,……
第1000个小批量训练集
,…,
。
然后,在每次迭代过程中随机选择一个当前未参与训练的小批次训练集进行搜索方向的更新计算。假设当前模型参数向量为 β k 并选择第 i 个批次进行方向更新,则有
据此计算更新方向为
依据上述更新方向对当前参数向量 β k 进行更新
当所有划分的批次均参与了训练过程而未达到算法终止条件时,对 S′ 进行随机排序后重新划分成若干小批次样本集合,并重复上述搜索方向和参数更新过程,直至达到算法终止条件时结束算法。
随机梯度法每次仅选择单个或小批量样本来确定搜索方向和实现对模型参数的更新,可有效减少花费在梯度计算上的时间。该方法的计算结果存在一定的随机性,即不能保证每次更新的方向都是目标函数减小的方向,参数更新偶尔会使得目标函数取值增大。然而,取自相同总体的样本服从相同的概率分布,并且所有训练样本均参与迭代计算,故在参数更新的整体过程中,随机梯度下降法可以保证目标函数的取值不断下降。随机梯度法的参数更新过程如图2-6所示。
图2-6 随机梯度法参数更新示意图
另一方面,随机梯度方法的参数更新方向存在着一定的随机性,这也为避免陷入局部最优提供了可能。正常情况下梯度下降算法很难跳出局部最优,但是随机梯度法打破了参数更新的稳定性,这有可能会使迭代过程中出现跳出当前局部最小值范围的情况。
目前,很多机器学习任务的训练样本集规模都很大,由以上分析可知随机梯度方法能够在保证优化效果的前提下有效减少在大规模训练样本集上对模型进行优化的计算量。因此,随机梯度法是目前最常用的模型优化方法之一。
有些机器学习模型含有隐含变量,通常难以对这些具有隐含变量的模型直接进行参数求解或估计,此时可通过迭代求解隐含变量(或其函数)数学期望最大值的方法实现对带隐含变量模型的参数优化求解。这类优化方法通常称为最大期望法或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 | β ),故将上述似然函数转化为
考虑其对数似然,有
根据最大似然估计思想,只需求得对数似然ln L ( β | X )的最大值即可得到模型参数 β 的最大似然估计值。但是由于上式存在隐含数据 Z 和积分的对数,直接求解其最大值较为困难,故用迭代逼近方法实现对数似然最大值的估计。为此,将对数似然做如下变形
其中, p ( Z )为隐含数据 Z 的某一分布。
由于对数函数为凸函数,故如下不等式成立
由此可得对数似然的下界函数
取 p ( Z )= p ( Z | X , β t ),则可得
略去下界函数 B ( β , β t )中与待求参数向量 β 无关的项,得到如下 Q 函数
上式表示对隐含变量 Z 的函数 L ( β | X , Z )在概率分布 p ( Z | X , β t )下的数学期望。由于 B ( β , β t )≤ln L ( β | X ),故可通过迭代选取不同下界函数 B ( β , β t )最大值的方式逐步逼近对数似然ln L ( β | X )的最大值,迭代逼近的具体过程如图2-7所示。
图2-7 EM算法逼近似然函数最大值
由于 Q ( β , β t )是通过省略 B ( β , β t )中与待求参数向量 β 无关项得到的,用 Q ( β , β t )迭代求解对数似然最大值与用 B ( β , β t )求解对数似然最大值具有相同的效果,故EM算法通常直接使用 Q 函数进行优化计算。由以上分析可得EM算法的基本步骤如下。
(1)设置初始参数 β 0 和迭代停止条件;
(2)E步(期望步):根据可直接观测数据 X 和当前参数向量取值 β t 计算 Q ( β , β t ):
(3)M步(最大化步):最大化 Q ( β , β t )并根据 Q ( β , β t )的最大值更新参数 β t 的取值:
(4)判断是否满足迭代停止条件,若满足则停止迭代;否则令 t = t +1并返回步骤(2)。
EM算法可随机选择初始参数 β 0 ,但应注意EM算法对初始参数 β 0 具有一定的敏感性。如果初值 β 0 选取不当,则可能会使迭代结果陷入局部最优。
通常将迭代停止条件设为相邻两次参数值变化不大或前后两次参数更新使得 Q ( β , β t )值变化很小时停止更新,即有
【例题2.9】 假设 X 1 , X 2 取自指数分布 Y ( θ )= θ e - xθ ,且 X 1 , X 2 不相关。若 X 1 =10而 X 2 缺失,试用EM算法实现参数 θ 的最大似然估计。
【解】 由于 X 1 , X 2 为取自指数分布 Y ( θ )= θ e - xθ 的离散值,则有
由于
X
1
,
X
2
同分布,则有
,从而有
为求出使得
Q
(
θ
,
θ
t
)最大化的参数
θ
,可将使得
Q
(
θ
,
θ
t
)最大化的参数
作为第
t
+1次的估计值
θ
t
+1
,即得迭代公式
当迭代收敛时,有 θ * =2 θ -/(10 θ * +1),解得 θ * =0.1,即迭代算法求得参数估计值收敛于0.1,即参数 θ 的最大似然估计值为 θ * =0.1。
蒙特卡洛法是一类以概率统计理论为基础的随机型数值算法。一般来说,当一个随机算法满足采样次数越多,其输出结果越近似最优解这一特性时,便可称为蒙特卡洛法。该方法通常将待求解问题与某一概率模型联系起来,利用从大量样本中获得的信息来完成对所求参数的概率估计,由此实现对实际问题的求解。
例如,可用蒙特卡洛法近似计算一个不规则湖面的面积。假设围住湖面的长方形面积为 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 1 , x 2 ,…, 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 2 ( x ) p 2 ( x ) /q 2 ( 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 1 , X 2 ,…, 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 0 中随机选择多条马尔可夫链的初始状态, p 0 的具体形式为
其中,
表示在时序为0的初始状态分布下
X
=
X
i
的概率。
由于这些马尔可夫链的状态空间是一个
n
维空间,故可将概率分布
p
0
表示为一个
n
维向量
p
0
=(
,
,…,
)
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 0 ;
(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采样的改进,这里不再赘述。