一种基于序贯蒙特卡罗方法的贝叶斯优化方法与流程

本发明属于全局优化技术领域，具体涉及一种基于序贯蒙特卡罗方法的贝叶斯优化方法。

背景技术：

在科学和工程领域的许多问题可以被描述为找到一个很难去估计的未知函数的最小值或最大值。贝叶斯优化是一种广泛使用的概率方法，用于解决这一问题。为了解释未知目标函数的特征，高斯过程特别适用于模型预测的解释，并为学习和模型选择提供合理的框架。由于其在调制和学习方面的优势，它已成为机器学习中最流行的基于非参数核的全局优化方法。它广泛应用于许多应用或研究领域，例如，在强化学习、大数据、无线传感器网络以及许多其他工业应用领域，如天气预报、计算生物学、医疗应用、生存分析等。

为了找到函数的全局最优值，贝叶斯优化引入了一个采集函数，用于评估并平衡下一数据点的搜索和利用。几十年来，人们已经提出了许多采集函数的变体。最常用的是后验概率改进(pi)、期望改进(ei)、上界置信区间(ucb)、熵搜索(es)。目前，我们在贝叶斯优化上做了大量的工作，以提高高斯过程模型对大数据的计算效率的有效性，对非高斯似然的后验近似和协方差函数的超参数优化。尽管贝叶斯优化效率很高，但实践者仍然面临着几个问题。在实际实现中，贝叶斯优化通常需要在每次迭代中使用辅助全局优化器来优化采集函数，例如分割矩形(direct)算法和协方差矩阵自适应进化策略(cma-es)。引入的辅助优化器带来了其他问题：如何理解这种辅助确实找到了采集函数的最大值，以及辅助的使用可能是不必要的代价，因为它必须在每次迭代中执行。为了避免辅助优化器，munos提出了同步乐观优化(soo)，它能够在不了解其平滑性的情况下全局优化目标函数。随后作者提出了其他类似的方法如随机soo、并行乐观优化(poo)、贝叶斯多尺度乐观优化(bamsoo)等。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种通过序贯蒙特卡罗方法计算权重分布；然后通过学生t过程回归而不是高斯过程回归将改进的序贯蒙特卡罗方法扩展到全局优化中，可以在较少样本数的情况下更有效地得到最大值的分布，具有更强的探索能力和异常值适应性能的基于序贯蒙特卡罗方法的贝叶斯优化方法。

本发明的目的是通过以下技术方案来实现的：一种基于序贯蒙特卡罗方法的贝叶斯优化方法，包括以下步骤：

s1、建立目标函数，得到目标函数的后验分布；

s2、加入噪声信号，计算目标函数的边缘student-t分布；

s3、通过序贯蒙特卡罗近似的方法对目标函数进行优化。

进一步地，所述步骤s1具体实现方法为：

对于d维随机变量x，多变量student-t分布由下式给出：

其中，md(x，μ，∑)＝(x-μ)^t∑^-1(x-μ)；md(x，μ，∑)表示相对于协方差矩阵σ，从x到均值μ的马哈拉诺比斯距离平方；马哈拉诺比斯距离测量当前参数(μ,σ)对观测x的可解释性；附加参数ν>2，描述分布的自由度，并在变化时充当平滑因子；当ν→∞时，student-t分布收敛于具有相同均值和协方差矩阵的高斯分布；

通过均值函数m(x)和协方差函数k(x,x′)的student-t过程来构建目标函数f(x)：

yi＝f(xi)，i＝1，2，...，p

来自student-t过程先验的样本也是student-t分布：

令均值μ为0，然后给定观测数据集{x，y}＝{(xi，yi)，i＝1，2，...，p}，得到后验分布：

其中θ是协方差矩阵k中的超参数和自由度ν的集合；

基于历史数据集{x，y}，给定新的输入数据潜在函数f*的后验预测为：

其中，

其中，k*是p维向量，其中第i个条目是k(x*，xi)，具体定义为：

k*＝[k(x*，x1)k(x*，x2)...k(x*，xn)]

k是由成对核k(xi，xj)构成p×p的协方差矩阵，具体定义为：

y是p维观测值向量；我们令均值μ和μ*为0和0。

进一步地，所述步骤s2具体实现方法为：

将噪声纳入核心函数：

其中是噪声标度，δij是克罗内克函数；

x为具有p自变量的student-t分布，其具有自由度v，均值向量μ和协方差矩阵σ以及以下分布：

其中xi具有自由度v的pi变量student-t分布，均值向量μi＝μ(xi)，协方差矩阵

∑ii＝∑(xi，xi)，并且xj具有类似参数(p＝pi+pj)的pj变量student-t分布；

f(xi)的边缘student-t分布由下式给出：

进一步地，所述步骤s3具体实现方法为：

s31、进行最大值分布估计：对于固定粒子数n，将它们随机分成n个区间，n≤n；对于区间xi中的每个粒子i，随机选择另一个区间xj，从以下联合分布中取样：

对于每个样本如果则将粒子从区间xi移动到xj；重复上述步骤，直到满足预设的停止标准，最大值的分布密度可由n个区间中粒子数的直方图构成；

s32、序贯蒙特卡罗的方法对样本粒子点进行抽取；

定义重要性权重为：

其中，q(xn)是在输入空间中展开的粒子的分布：

其中m表示成分指数，m表示混合成分的总数，αm表示满足αm>0的第m个成分的概率质量，每个m∈{1,2,…,m}和由分布q(·)的所有可调参数组成；

q′(xn)为建议分布：

q′(xn)＝αpm(xn)+(1-α)q(xn)

其中，0≤α≤1，为预设参数，用于调整新形成空间和旧形成空间之间的重采样率；

推导出重要性权重的最大值的分布为：

使用核密度估计来适应连续环境以及平滑绘图后得到：

当重新采样旧粒子时，从估计的最大值分布中取样α比例的样本，并且从student-t分布中取样(1-α)的新样本，得到权重更新公式为：

当新粒子的输出值大于旧粒子的输出值时，旧的粒子被新例子取代，并且利用上述权重更新公式对权重进行更新。

本发明的有益效果是：本发明使用序贯蒙特卡罗方法来解决黑盒函数的最大化问题，通过序贯蒙特卡罗方法计算权重分布；然后通过学生t过程回归而不是高斯过程回归将改进的序贯蒙特卡罗方法扩展到全局优化中，可以在较少样本数的情况下更有效地得到最大值分布，具有更强的探索能力和异常值适应性能。

附图说明

图1为本发明的基于序贯蒙特卡罗方法的贝叶斯优化方法的流程图；

图2为采样高斯过程求解目标函数最大值分布的示意图；

图3为蒙特卡罗方法、准蒙特卡罗近似和二分法近似的最大值分布比较图。

具体实施方式

下面结合附图进一步说明本发明的技术方案。

如图1所示，一种基于序贯蒙特卡罗方法的贝叶斯优化方法，包括以下步骤：

s1、建立目标函数，得到目标函数的后验分布；具体实现方法为：

对于d维随机变量x，多变量student-t分布由下式给出：

其中，md(x，μ，∑)＝(x-μ)^t∑^-1(x-μ)；md(x，μ，∑)表示相对于协方差矩阵σ，从x到均值μ的马哈拉诺比斯距离平方；马哈拉诺比斯距离测量当前参数(μ,σ)对观测x的可解释性；附加参数ν>2，描述分布的自由度，并在变化时充当平滑因子；当ν→∞时，student-t分布收敛于具有相同均值和协方差矩阵的高斯分布；

通过均值函数m(x)和协方差函数k(x,x′)的student-t过程来构建目标函数f(x)：

yi＝f(xi)，i＝1，2，...，p

来自student-t过程先验的样本也是student-t分布：

为计算方便，令均值μ为0，然后给定观测数据集{x，y}＝{(xi，yi)，i＝1，2，...，p}，得到后验分布：

其中θ是协方差矩阵k中的超参数和自由度ν的集合；

基于历史数据集{x，y}，给定新的输入数据潜在函数f*的后验预测为：

其中，

其中，k*是p维向量，其中第i个条目是k(x*，xi)，具体定义为：

k*＝[k(x*，x1)k(x*，x2)...k(x*，xn)]

k是由成对核k(xi，xj)构成p×p的协方差矩阵，具体定义为：

y是p维观测值向量；我们令均值μ和μ*假定为0和0。

s2、加入噪声信号，计算目标函数的边缘student-t分布；具体实现方法为：

将噪声纳入核函数：

其中是噪声标度，δij是克罗内克函数；

x为具有p自变量的student-t分布，其具有自由度ν，均值向量μ和协方差矩阵σ以及以下分布：

其中xi具有自由度ν的pi变量student-t分布，均值向量μi＝μ(xi)，协方差矩阵∑ii＝∑(xi，xi)，并且xj具有类似参数(p＝pi+pj)的pj变量student-t分布；

f(xi)的边缘student-t分布由下式给出：

s3、通过序贯蒙特卡罗近似的方法对目标函数进行优化；具体实现方法为：

s31、进行最大值分布估计：对于固定粒子数n，将它们随机分成n个区间，n≤n；对于区间xi中的每个粒子i，随机选择另一个区间xj，从以下联合分布中取样：

假设每次采样到的样本为如果则将粒子从区间xi移动到xj；重复上述步骤，直到满足预设的停止标准，最大值分布可由n个区间中粒子数的直方图构成；

s32、采用基于序贯蒙特卡罗的方法对样本粒子点进行抽取；

定义重要性权重为：

其中，q(xn)是在输入空间中展开的粒子的分布：

其中m表示成分指数，m表示混合成分的总数，αm表示满足αm>0的第m个成分的概率质量，每个m∈{1,2,…,m}和由分布q(·)的所有可调参数组成；

q′(xn)为建议分布：

q′(xn)＝αpm(xn)+(1-α)q(xn)

其中，0≤α≤1，为预设参数，用于调整新形成空间和旧形成空间之间的重采样率；

推导出重要性权重的最大值分布为：

使用核密度估计来适应连续环境以及平滑绘图后得到：

当重新采样旧粒子时，从估计的最大值分布中取样α比例的样本，并且从student-t分布中取样(1-α)的新样本，得到权重更新公式为：

当新粒子的输出值大于旧粒子的输出值时，旧的粒子被新例子取代，并且利用上述权重更新公式对权重进行更新。

下面进一步对本发明的贝叶斯优化方法的推导进行说明。

(一)对于高斯过程(gp)的框架，用后验均值函数m(x)和后协方差函数k(x,x′)来假设未知函数f(x)，即函数f(x)遵循高斯过程，定义为：

高斯过程回归的目标是基于观测数据x，找到待观测输入数据x*的后验概率p(f*|x,y)；观测数据x对应的结果y和一些试验输入x*分别记为f＝f(x)和f*＝f(x*)；

当数据x被观察到输出y时，函数f被认为是高斯函数：

利用函数f的高斯过程属性得到函数f和f*是联合高斯分布：

f和f*分别是n维观测输入数据x和n*维未知输入数据x*对应的函数值，x和x*分别是n维观测输入数据和n*维待观测输入数据(测试数据)；k(x，x*)表示观测数据x和测试数据x*之间的协方差矩阵(n×n*)，k(x，x)、k(x*，x*)、k(x*，x)则分别是观测数据x、测试数据x*以及测试数据x*、观测数据x之间协方差矩阵(n×n、n*×n*、n*×n)；为简洁起见，设k＝k(x，x)、k**＝k(x*，x*)、k*＝k(x，x*)、

协方差矩阵k，k*，k**分别定义为：

在公式(3)中，f为隐变量，不能直接被观测到，只有对y噪声的观察；

考虑到在公式(2)中显示的高斯似然，将公式(3)变成：

通过高斯条件作用得到：

其中，

为简单起见，在公式(1)m(x)≡0中设置了先前的均值函数；因此，公式(7)的均值函数都是零向量；假设我们有一组数据点x1:n和它们对应的近似概率密度f1:n；一个新的数据点x*也遵循f*＝f(x*)，类似于公式(5)，f1:n和f*都是联合高斯函数：

其中

k*＝[k(x*，x1)k(x*，x2)...k(x*，xn)](9)

测试输入x*的后验分布p(f*|x,y,x*)根据公式(6)来预测：

其中期望均值和方差为：

从公式(11)可以看出，未知输入数据x*的后验概率分布p(f*|x,y,x*)的均值和方差均独立于f*的，我们只需利用观测数据x的信息、对应的输出y和测试点x*来计算f*。然而，这些有利的特性给我们带来方便时，它也产生了一些不足。

(二)将找到函数f(x)的最大值的问题转化为建立概率密度函数p(x*)，将这种概率称为最大值分布，表示pm(x*)；采样高斯过程求解目标函数最大值分布如图2所示；获得最大值的分布与全局优化类似，除了全局最优和其位置都是变量。依靠基于序贯蒙特卡罗的方法依次选择最大值点，可以得到其最大值的分布密度；而不是依赖于最大值点选择的采集函数标准。首先考虑直接和直观的准蒙特卡罗采样方式，然后我们转向二分法抽样，最后提出了基于序贯蒙特卡罗的方法。

具体包括以下子步骤：

(1)将蒙特卡罗的思想应用于蒙特卡罗抽样方法。整体算法总结如下：对于固定粒子数n，将它们随机分成n个区间，n≤n；对于区间xi中的每个粒子i，随机选择另一个区间xj，从以下联合分布中取样：

对于从中抽取的每个样本如果则将粒子从区间xi移动到xj；重复上述步骤，直到满足预设的停止标准，最大值分布则由n个区间中粒子数的直方图构成。

通常，区间数n等于训练点的数量，同时不大于粒子数n。对于较大的粒子数，算法在有限的输入执行中收敛。图3显示了蒙特卡罗近似方法、理论极限近似方法和采用10000个样本函数的暴力近似方法的最大值分布。蒙特卡罗方法的采样数为10000，执行10轮。很明显，蒙特卡罗方法对最大值分布的置信度最低。尽管分布会收敛到近似最大值分布，但这种原始算法在计算上是密集的，因为它会产生样本成本(n³)。这种算法不仅收敛速度慢，而且结果的置信度低于其他近似值。

(2)为了提高算法的效率和最大值分布的置信度，考虑序贯蒙特卡罗方法。首先简要讨论序贯蒙特卡罗的方法。

序贯蒙特卡罗(smc)或粒子滤波方法基于迭代重要性采样，应用于利用一系列时变目标分布过滤问题。得到后验分布的经验点质量近似模型为：

其中是描述训练点xn和δ(·)特征的样本或粒子是狄拉克-德尔塔函数；

由于我们不知道真正的后验，我们通常定义一个易于实现的分布，称为建议分布，并绘制i.i.d样本从中估计函数的均值为：

其中，是归一化的重要性权重：

重要性权重表示具有最大值的粒子的概率；在粒子滤波方法中，它被定义为后验和建议分布的比率。

为了找到最大值的概率，或简单地找到最大值的概率位置，将重要性权重定义为：

其中，q(xn)是在输入空间中展开的粒子的分布，通常使用简单的平坦分布q(xn)＝c，c>0；

建议分布的选择q′(xn)在性能中起着至关重要的作用。自然选择只是估计的最大值分布。然而，估计的最大值分布不能是真正的最大值分布，并且可能具有很大的偏差，这将完全导致错误的分布。为避免这种情况，使用混合重要性抽样和防御重要性抽样的想法来权衡先前估计分布与原始平面分布之间的重新取样，因此得到建议分布q′(xn)的计算方法为：

q′(xn)＝αpm(xn)+(1-α)q(xn)(17)

其中，0≤α≤1，是手动设置参数，用于调整新形成空间和旧形成空间之间的重采样率；

推导出具有重要性权重的最大值分布密度为：

使用核密度估计来适应连续环境以及平滑绘图得到：

基于公式(17)和(18)可以计算公式(16)中的新的重采样粒子的权重，但是该权重计算的代价较高，为了降低计算资源，当重新采样旧粒子时，我们采用从核函数中抽取一个α的比例和(1-α)的比例来自q(xn)中的新样本使得权重计算更有效，如下所述：

当新的重采样粒子通过公式(12)获得比旧粒子更高估计时，递归地更新权重；即，当新粒子的输出值大于旧粒子的输出值时，粒子被取代，并且权重被更新。从理论上讲，新粒子与旧粒子竞争的次数越多，我们获得的性能就越好。但是，鉴于效率和性能之间的平衡，只有一个新样本更好。

(三)尽管高斯过程在贝叶斯优化中取得了巨大成功，但仍存在一些缺点。首先，由于高斯分布的概率，目标函数的异常值的后验通常具有低概率，而不管采样数据如何。第二，从公式(11)中，我们注意到后验方差不依赖于返回的客观值，而是仅取决于输入评估位置，这意味着如果样本集中的客观值与在高斯先验的预期下差别很大，则后验方差不会更高。因此，我们考虑一个更通用的模型，即student-t过程(tp)。

(1)student-t过程(tp)是贝叶斯优化的另一种选择。对应于多元高斯分布表示gp的边际分布，多变量student-t分布也表示stp。对于d维随机变量x，多变量student-t分布由下式给出：

其中，

md(x，μ，∑)＝(x-μ)^t∑^-1(x-μ)

md(x，μ，∑)表示相对于协方差矩阵σ，从x到均值μ的马哈拉诺比斯距离平方；马哈拉诺比斯距离测量当前参数(μ,σ)对观测x的可解释性；附加参数ν>2，描述分布的自由度，并在变化时充当平滑因子；当ν→∞时，student-t分布收敛于具有相同均值和协方差矩阵的高斯分布；

类似于公式(1)，考虑目标函数f(x)，实现具有均值函数m(x)和协方差函数k(x,x′)的student-t过程：

yi＝f(xi)，i＝1，2，...，p

与高斯过程回归类似，来自student-t过程先验的样本也是student-t分布：

为简单起见，假设均值μ为0，然后给定观测数据集{x，y}＝{(xi，yi)，i＝1，2，...，p}，得到后验分布为：

其中，θ是协方差矩阵k中的超参数和自由度ν的集合；为简洁起见，除必要外，一般将省略此参数。

基于历史数据集{x，y}，给定新的输入数据潜在函数f*的后验预测为：

其中k*是p维向量，其中第i个条目是k(x*，xi)，如公式(9)中所定义；k是由成对核k(xi,xj)构成p×p的协方差矩阵，并且在公式(4)中定义；y是p维观测值。均值μ和μ*假定为0和0。

用公式(10)比较公式(25)和公式(26)，显然tp和gp具有相同的后验均值，而后验协方差不同。与gp不同，tp的后验协方差具有额外的主导项，导致与输出相关的后验协方差，这意味着更好地考虑不确定性。该主导项对后验协方差进行缩放，并且马哈拉诺比斯距离y^tk^-1y确定与gp相比的尺度水平。后协方差的尺度随着马哈拉诺比斯距离的减小而减小。

随着公式(24)并且越来越多的数据被采样，可以估计样本的最大值，最后通过直方图获取最大概率。

值得对自由度ν进行更多的认识，它是用来控制重尾过程行为的自由度。ν越低，尾部越重，反之亦然。对于ν→∞，tp退化为gp，这意味着gp回归可以被解释为tp回归的特例。直觉上，功能的预测方差会受到自身观察的影响是合理的。但对于gp来说，情况并非如此，因为它只与输入点有关。虽然在某些应用中，在获得函数观察之前知道预测方差是有利的，但是当呈现有限的数据集并且需要准确的不确定性量化时，gp是低效的。tp填补了这一空白，并且优于披露更多的不确定性。

(2)直到现在，似乎tp的预测与gp的预测一样简单。然而，在实际应用中，不能直接观察到函数值f(xi)，而是仅获得具有加性噪声的yi。由于信号被噪声纠缠，导致观测结果不再是独立信号和噪声之和，因此对tp的噪声观测更为复杂。为了近似具有独立学生t噪声的后验，qingtangtao等提出了一种具有学生t似然(tprt)的student-t过程回归方法。由于其复杂性，本发明将噪声纳入核函数用以简化计算：

其中是噪声标度，δij是克罗内克函数；这里的噪声与隐函数不相关，但不是独立的，因为缩放参数ν受平方指数核以及噪声核的影响。利用这种噪声模型，两个独立的多元student-t分布的总和是可分析处理的。

假设x具有p自变量student-t分布，其自由度为v，均值μ和协方差矩阵σ以及以下分布：

其中xi是具有自由度v的pi维student-t分布，均值μi＝μ(xi)，协方差矩阵∑ii＝∑(xi，xi)，并且xj具有类似参数(p＝pi+pj)的pj变量student-t分布。然后f(xi)的边缘student-t分布由下式给出：

f(xj)的分布类似于f(xi)的分布。

(3)现在考虑序贯蒙特卡罗近似。它非常类似于基于高斯过程的蒙特卡罗方法。然而，我们对建议分布做出了不同的选择，因为它对粒子滤波至关重要。建议分布与后验的相似性越高，算法的性能越好。

由于假设目标函数有一些噪声观测值，必须利用这些噪声观测值来进行预测。嘈杂的测量通常具有更多的波动，这些波动不能通过高斯分布描绘，沿着平坦分布。student-t分布具有比高斯分布更高的峰度，能更好地表征波动的特征。此外，student-t分布使得去除异常值更有可能，并且当观察到的样本变化时，后验方差增大或减小。因此认为student-t分布作为建议分布的选择使我们更有利于算法更有效地探索输入空间。

注意我们从student-t建议分布中抽取样本发生在迭代阶段，而不是在初始化。在初始化时，仍然从均值分布中抽取样本，因为所有权重都设置为相同的值。

在初始化之后，考虑权重来选择随机粒子(尽管在初始化步骤中权重是相同的，但是它们在迭代时变化)。然后，基于(20)得到粒子的新样本的权重，其中采用具有混合和防御重要性抽样的想法。也就是说，从估计的最大值分布产生具有α比例的新样本，部分(1-α)比例的样本从student-t分布产生。通常情况下，为了获得更好的性能，也可以使用student-t分布的混合。也就是公式(17)中的分布q(x)由如下公式给出：

其中m表示成分指数，m表示混合成分的总数，αm表示满足αm>0的第m个成分的概率质量，每个m∈{1,2,…,m}和由建议分布q(·)的所有可调参数组成。由于计算成本，我们只采用公式(21)中所示的student-t分布。

为了判断当前样本是否更好地表达其对应粒子基于公式(28)来评估输出概率。并比较两个输出以确定是否更换旧粒子。最后，权重由公式(20)更新。

本发明阐述的是序贯蒙特卡罗方法在高斯过程回归中的应用。我们首先改进了迭代顺序蒙特卡罗方法时用于计算权重的分布。然后，我们通过student-t过程回归而不是高斯过程回归将序贯蒙特卡罗方法扩展到全局优化中。所提出的算法都具有与传统ucb，pi和ei更优异的性能。我们发现，由于全局优化的不确定性，student-t过程因其更强的探索能力和异常值适应性而发挥更强大的优化工具。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。