当遗传算法遇到 Bayesian Optimization

前文很多篇文章详细介绍了如何将遗传算法应用到结构设计的优化中。遗传算法有很多优点，比如能够处理离散非线性问题，能够跳过局部最优收敛于全局最优等。但是同时它有一个致命的缺点，就是对样本空间的大量搜索会耗费非常多的时间，特别对于实际的结构设计优化问题，算一个普通的模型就得大概两分钟。假如给出一个8个变量的优化问题，每一维变量假设可以取10个不同的值，假设算一个模型花费2分钟，将每一种可能都算一遍，就得花费380余年。

如果算模型不需要时间就太好了，那是不可能的。由于在找到全局最优之前，搜索优化的过程中并不需要知道模型精准的观察值，只要趋势没错就行，因此我们可以想办法用不花时间的近似模型来代替真实的模型。近似的方法在数学上就是逼近，逼近是指什么，古人又是如何逼近的呢？

我们知道已知三个点，（x1， y1），（x2，y2），（x3， y3），并且已知未知曲线为二次多项式的时候，

我们可以通过解方程组得到这个曲线：

但是拉格朗日别具一格地通过三根曲线叠加得到了想要的解：

其中y1f1经过（x1，y1），y2f2经过（x2，y2），y3f3经过（x3，y3）。

这也就是著名的拉格朗日插值法。

但是拉格朗日插值每增加一个点，就要重新算一次，为了解决这个问题，就出现了牛顿插值法。

牛顿插值问题可以总结为：已知n+1个点相对多项式函数f(x)的值为：

求此函数f(x)。

如果将此问题稍微改造一下，要求函数f（x），且它在点x0，x0+△x，x0+2△x，x0+3△x，x0+4△x，x0+5△x…的值已知，就求出了著名的泰勒公式。

翻开干干净净的高数书惊喜地发现：

哦买噶，原来这些曲线都是可以用多项式来逼近的呀=>=

如果我们将x换成未知的f（x），最古老的可能是傅里叶级数，傅里叶发现任何周期函数都可以展开为：

我们已经大概知道了古人是如何逼近的，但是非常遗憾的是，现实中的优化问题的变量大多不止一个，当变量空间的维数大于1时，以上的插值逼近方法突然就完全失效了，无解或无唯一解。怎么办呢？

以地质学中的插值为例：地质空间中要在平面上确定一点，至少得要两个变量（x， y）。当我们在某个区域内勘探了n个点，得知了n个点每个点有多少石油，如果我们想知道其它点有多少石油，怎么办，对该区域每一个点都进行勘探肯定是不可能的，严格来说一片区域存在无穷多个点，最起码对于离已经勘探过的点很近的点就没有必要勘探了。那么多近才能算很近呢？

这个问题背后实际上是要我们对未知点和一个已知点之间的相关关系进行一个估计。显然石油一般都是连续的，如果一个已知点的石油含量有8000，我们就可以认为，离这个点越近的未知点的含量越接近8000。人们自然地提出了反距离插值方法：

其中zi是已经勘探过的点。

这样插值显然很粗糙，平面内各点的相关程度不可能永远都是与距离成反比。那么如何真实地反应各点的相关性呢？地质学家Kriging重新思考了这个问题。

我们知道，虽然现在都是确定的，但是在很多很多年前，每一个点在今天的石油存储量都是一系列随机事件发生后的实现。有没有石油，有多少石油，取决于此处数万年前有无森林，是否发生了适合的构造运动，是否有合适的温度使其从有机质转化为石油，覆土是否不是砂性的等等。

Kriging别具一格地将视角置身于万年之前，假定这一切都还没有发生，未勘探的点的石油存储量都是一个随机函数的实现，那么要求这些随机函数的值就变成了用估计的方法，在已知一些点的实现值（实际上是勘察得到的）的条件下，估计每一个点的随机函数的实现值，也就是求以下的条件数学期望：

每一个点是相互相关的。

这个解是一个最小方差无偏估计，然而求解这个几乎是不可能的，因为首先要求随机函数F（x）在n个点上的联合分布，就算知道了联合分布求条件期望也很复杂，需要对n个点求n重积分。

为了解决这个问题，Kriging假定随机函数符合以下线性模型：

在此条件下求未知点的最小方差线性无偏估计。推导过程变得容易且清晰。

以下小编绕开地质统计学原始的推导，别具一格地从线性回归模型出发，结合我们实际的结构设计优化问题，为大家梳理一下推导过程，以及我们实际应用时的近似逼近思路。

回归分析结构的小伙伴一定都很熟悉，比如如果我们想知道混凝土轴心抗压强度与立方体抗压强度的关系，我们做很多很多的试验，得到下面的结果：

然后我们用如下回归模型就可以得到它们之间的相关关系：

《规范》以大量试验为基础，给出了b的估计值以及取值方法。

所以回归分析就是研究相关关系的数学工具，而Krijing要研究的就是已勘探点和未勘探点的相关关系。

线性回归模型可以写成如下形式：

其中，fh（x）是关于x的线性或非线性模型，前面的是未知系数，最后一项是服从N（0，v）分布且独立的误差项。

在梳理推导之前，我们再想一想我们实际面临的结构设计优化问题，

我们知道当我们做以上混凝土抗压强度的试验，由于各种试验误差的存在，每次相同的试验条件下，试验结果也不可能完全相同，结果误差独立且服从正态分布，我们可以通过大量的试验来消除误差的影响。但是当我们跑结构模型的时候，只要变量相同，结果肯定一样，此时就不存在独立随机的误差，此时的误差是我们逼近产生的误差，这种误差不再满足相互独立的假定，可以写成关于变量的函数。

另外我们无法每次都知道像以上混凝土抗压强度间a+bx那样简单纯粹的关系，也就是线性回归模型的第一项的fh（x）我们没办法知道，况且如果我们知道的话就不存在如何逼近的问题了。

对于相互相关的逼近误差，我们假定其和已勘探点的距离直接相关，但不不用欧拉距离，因为欧拉距离默认地将每一维变量的敏感程度视为一致，我们假定：

于是θh就可以反应每个变量的敏感程度，而ph反应了平滑程度。两点之间逼近误差的相关关系就为：

那对于回归模型的第一项我们怎么办？人们发现，在定义了以上的相关的误差项之后，误差项表现出了极强的逼近能力，将fh简单的选为常数就能得到很好的结果。也就是我们不用再纠结回归模型的第一项，也就是说像前面混凝土强度的例子用到的模型，a+bx，我们甚至可以不要bx，只保留a。

于是，我们的回归模型变成了这样：

我们的回归分析任务变成了这样：

对一片区域，我们已经勘探了如上图的小黑点，知道了其石油存储量，现在我们要对这片区域的石油存储量的平均值进行估计，同时根据已知点，求出未知点偏离均值程度的最小方差线性无偏估计，而每个点相互之间对均值的偏差的相关关系我们是可以知道的，如上所述和距离直接相关。

至于最佳线性无偏估计，如下图例子：

第一个射击结果的均值偏离了中心均值10分，期望不一致就是有偏，第二个期望一致就是无偏，而第三个方差更小就更优秀，自然方差最小就是最佳。

于是我们推得了y（x*）的最佳线性无偏估计为：

其中R是相关系数矩阵，就是已知点两两之间误差的相关关系，

r为x*和所有已知点误差的相关关系

到此为止，貌似我们已经完成了我们的逼近任务，只要输入x*，好像就可以得到我们想要的估计值y（x*）。但是前面在假定误差相关关系的时候，我们定义了敏感程度相关的参数θh和平滑程度相关的ph，在完成逼近之前，我们要知道这些参数的值。

为了找到这些参数的值，我们采用最大似然估计，就是将视角置于万年之前，对如今已经勘探了的点，先求出其实现值同时发生的事件的概率。既然它已经发生了，我们就认为这种事件发生的概率是很大的，那么问题就转换为求使这个函数值最大的参数的值。

它们同时发生的概率为：

写成对数函数的形式：

然后对它对参数求导等于0，就求出了：

代回上式就得到同时发生的概率为：

这个函数仅跟未知参数θh，ph有关，到这里就没办法直接解了，但是我们可以用遗传算法在找到θh，ph的最优值。

找出参数值后，求得μ，我们的逼近任务就算完成了。

我们用此方法对Branin函数进行逼近，如下图，左边为真实函数值，右边为逼近值：

十个样本时

二十个样本时

四十个样本时

可以发现随着样本点逐渐增多，逼近模型也越来越精确，四十个样本时，逼近模型几乎和原函数值一模一样。

现在考虑将此方法应用到结构设计优化问题中去，我们知道优化过程并不需要模型特别的精确，而我们的目标就是算尽量少的结构模型。如果用以上10个样本点时的逼近模型，遗传算法也能很快找到全局最优的大概位置（左边那块），然后再在这块加入新的样本点，模型在这里就会变得更精确，慢慢的就找到了全局最优。这是一个不错的思路，这时问题就变成了在什么位置加入新的样本点。

前面推导逼近模型的时候我们知道，我们最终的逼近值是建立在对参数θh，ph的最大似然估计上的，也就是说我们将y（x*）的最大似然值作为我们的逼近值，或者说我们对y（x*）似然值最大时最有信心，这个毋庸置疑，如果我们换一个角度，我们将y（x*）的最大似然函数曲线画出来：

对于这两种情况，我们对哪个更有信心呢？

显然我们都取y(x*)的最大似然值作为我们最终的逼近值，而左边这种相对平缓的情况，我们虽然确定了最大值的位置，但是其他位置和最大值也没差多少，说不定就比最大值位置的更接近真实值呢。而右边的情况，稍微挪一点位置，似然值快速下降，因此对比这两种情况，不难得出我们对于右边的逼近更有信心的结论。

这种信心的比较很有帮助，因为越是没信心的地方，就意味着它的逼近误差也就应该越大，信心越强的地方，逼近就越精确，也就是说我们可以对样本空间任意一点的逼近误差进行一个估计。

可以发现y（x*）的似然函数的二阶导数为：

我们用其倒数来评价对逼近的误差估计，如果我们直接推导回归模型的均方差，得到的结果：

和上式的倒数除去最后一项一模一样，而最后一项是由μ的不确定性导致，往往值非常小，可以忽略不计。我们将s（x*）称为点x*处逼近的标准误差。

上图是对某函数逼近的过程，绿色的线为标准误差，可以看出样本点越稀疏的地方，估计误差就越大。此时我们有了该在何处加点的思路，当然是估计误差最大的地方。但是我们何时停止这种加点呢，因为在遍历完样本空间之前，理论上永远存在估计误差。

我们进一步思考，由于最终的目的是找到全局最优，所以希望加入的点的实现值最好能够比已知逼近模型的最优值更优，标准误差给我们提供了评价任意点实现值比目前的逼近值更优的可能性大小的机会：

之前我们曾假定每个点实现值的误差服从N(0,v）分布，此时我们将v赋值标准误差s，如上图，那么该点实现值比目前的逼近值更优的可能性大小就为图中阴影部分的面积：

甚至我们可以求出该点能够提升值的期望：

然后我们在每一代此期望值最大的地方加入新的样本点，如下对一个简单函数的优化过程：