本文主要基于在2021年发表于 npj Computational Materials 上的文章 Bayesian optimization with adaptive surrogate models for automated experimental design [10]。

引言

无论是在日常生活之中，还是在学校课堂上，或者在学术研究之中，实验一直以来都帮助着人们理解各种形形色色现象之间的因果关系。通常来说，实验就是在一定条件下，通过控制变量产生不同的结果，来验证某种假设。比如物理可以通过实验可以探索物体运动的规律，材料学则可以分析某个化合物的性质。但早期的研究更注重通过观察来探究自然规律，比如古希腊科学家亚里士多德就通过大量的观察来进行研究，属于论证科学。再比如，阿基米德，他成功地结合了逻辑推理和观察实验这两种方法，但这只是为了检验逻辑推论 [18]。
随着时间的推移，人们逐渐开始不满足于已有的研究模式，于是培根等人开始推崇实验科学，强调“从实验中获得知识”。之后，伽利略也呼吁“科学的真理不应该在古代圣人的蒙着灰尘的书上去找，而应该在实验中和以实验为基础的理论中去找”。他初步开创了“观察--假设--推理--实验--得到结论”的研究方法，促进了科学研究的发展。如今，做实验已经成为无数物理，生物，化学，材料等学科学者经常要进行的工作。

Humen-centric vs Data-centric

但随着时代的发展，科技的进步，人们尝试解决的问题越来越困难，维度越来越高。许多时候，每次实验可能都对应着大量的时间和金钱的消耗，传统的最直接的枚举法，及尝试所有的可能情况，开始变得不切实际。所以人们开始对实验有了更高的要求，实验设计开始变得越来越重要。比如如果能通过更少的实验达到同样的效果，就能节约许多时间和金钱，可能一些之前预算内无法完成的研究变得可能。拿今天我们主要要聊的材料学来说，一个常见的问题是“在众多可能的材料结构中找到性质最优的那个”。如图 1 所示，我们想要在 和 中找到具有最优性质的结构。其中 ， 和 分别可以在红色，蓝色和黑色区域选择一种元素，基于材料学知识排除一些不稳定的组合后，还剩下403种可行的结构。于是问题就变成了如何在这403种材料里找到我们需要的那个。通常情况下，预算是无法支持我们尝试完全部的403种可能的结构。
除了枚举法，另一种方式是人们根据自己的经验 (human-centric) 来安排进行实验的化合物结构的先后顺序。这种方式优于简单粗暴的枚举法，尤其是对于经验丰富的学者来说。但这可能对进行实验的人会有较高的要求，同时当面对一些新的问题，以往的经验也不一定奏效，于是可能还是需要不少的尝试才能找到目标的结构。
不管是枚举法还是根据经验的实验设计，都有一个问题，就是没有充分利用已有的信息。假设我们已经进行了 次实验，并且记录好了相关的数据。那么当我们准备进行 次实验时，前 次的结果是可以对我们这次的实验选择进行指导的。比如从直觉上，如果我们发现了某种结构的化合物性质比较好，那么和它类似的化合物可能这种性质也会比较好。
这时候问题就变成了如何才能更充分利用已知信息。随着统计学、机器学习算法的发展，在对复杂函数进行建模的问题中，计算机相对于人脑的优势越来越明显。人们已经可以利用计算机和设计好的模型算法，更好地从数据中发掘出有价值的信息。“Taking human out of the loop" [14] 和 "data-centric" 的呼声越来越高，材料发现问题中的自动化实验设计的概念 [11] 也迅速成为一个活跃的研究领域 [1, 7, 9]。贝叶斯优化 (Bayesian optimization，BO) 就是其中一种让数据来指挥实验设计（自动化实验设计）的方法。
图1：MAX phases: Ternary (and higher order) layered carbides and nitrides with properties intermediate to those ofmetals and ceramics [15]。

浅谈贝叶斯优化

这里先介绍一下贝叶斯优化 (BO)，这里举一个简单的例子，比如我们想要找未知函数 的最大值，于是我们基于已经观测了的函数值对 进行建模，然后根据训练的模型找到最可能是最大值点的 然后进行采样，接着将得到的新的样本点加入已观测数据之中。不断重复这个过程，直到循环终止条件得到满足。如图 2 所示，蓝色曲线是真实的未知函数，橘黄色曲线和区域是根据模型拟合出来的曲线及其置信区间，黑色圆点代表已观测的数据，红色曲线对应计算的采样函数 expected improvement，而灰色三角形则为下一个采样点。可以看到，每一次灰色三角形都在红色曲线的最高点，然后将新的数据加入数据集中并更新橘黄色曲线和区域，然后得到新的红色曲线，然后不断循环这个过程。
图2：Schematic illustration of Bayesian optimization (BO).
BO 是一个常用的优化目标函数的流程。当我们想基于最少的信息，来知道某个函数的最大值点（或者最小值点），BO 是一个很好的选择。近年来，BO 在实验设计和超参数的调试中受到了极大的欢迎。具体来说，BO 一般包括以下步骤：
图3：贝叶斯优化 (BO)
从算法中我们可以发现，使用 BO 进行自动化实验设计，有两个非常重要的点需要提前设定：（1）基于数据进行训练的概率模型；（2）决定下一个采样点的采集函数。概率模型是用来建立自变量（实验前可知的特征）和因变量（感兴趣的材料性质）的关系。而采集函数则是我们选择下一个样本的标准，代表了我们如何定义出现最大值的可能性。
关于采集函数，常用的包括但不限于 expected improvement，probability of improvement, upper confidence bound，以及 Thompson sampling。使用者可以基于自己的需求选择其中之一。这里我们以 expected improvement (EI) 为例进行后续分析。EI 会从期望的角度选择对 函数值提高最多的样本点，其公式如下：
其中， 是已观测的函数值中的最大值， 是对基于观测数据得到的后验分布取期望，。

跳出高斯过程，也许别有洞天

前面已经提到，除了采集函数，BO 另一个关键点是其使用的概率模型。假设我们有自变量 以及因变量 ，我们想要找到两者之间的关系，并且当新的 到来时能较好地预测其对应的 。通常我们会假设因变量是关于自变量的函数再加上一个额外的噪音项 。通常一个非常常用的选择是高斯过程回归 (Gaussian Process Regression, GPR)，不管是 python 还是 R 中都有成熟的包，所以使用起来比较方便，同时效果也不错。

并非万能的高斯过程回归

在 GPR 中，我们会给未知函数 加一个高斯过程的先验（ 是均值函数， 是所选的核函数）：
当新的输入 ，其预测后验分布为：
常用的均值函数是常数函数，而对于核函数，我们有很多选择，包括但不限于：
  • Radial-basis function (RBF) kernel：可用于拟合平稳和各向同性的特征。
  • RBF kernels with automatic relevance detection~(ARD)：在 RBF kernel 的基础上为每个特征分配不同的尺度参数，从而进行变量选择。
  • Dot-product kernels：可用于捕捉非平稳的特征。
  • Deep kernels：是一种更加灵活的核函数，但是参数较多，一般需要更多数据进行训练来达到较好的效果。
虽然 GPR 非常流行并且在许多场景下有不错的效果，但其并不是万能的，也有一些缺点。比如平稳以及各向同性的 GPR 面对高维问题时表现会下降，尤其当许多不重要的自变量混在了真实变量之中但建模者并不知道哪些是关键变量的时候。可是这种情况在材料学中却十分常见。所以如何在建模的时候，同时进行变量选择就变得十分重要了。近年来有一些基于 GP 的模型在尝试变量选择，比如前面提到过的 automatic relevance detection [2]。还有基于混合高斯模型和贝叶斯模型平均的方法 [15]。但在所有的这些基于 GP 的模型中，他们使用的核函数通常给拟合的函数引入了平滑性和连续性的假设。但在许多材料学的问题中，真实函数可能并不是平滑的或者函数值会有突然的跳跃。除了上述的研究，还有许多关于灵活的非平稳协方差核函数的文献 [12]，比如 Deep kernels [17]，但由于参数比较多，其在数据量较少的情况下不一定能获得好的表现，而在许多材料学研究中，昂贵的实验开销导致数据量不会太大。所以也许跳出并非万能的 GPR，会有不一样的风景。

贝叶斯多元自适应回归样条

贝叶斯多元自适应回归样条 (Bayesian multivariate adaptive regression splines, BMARS) [5] 就是一个不错的替代 GPR 的选择。BMARS 是一个非常灵活的非参数回归模型，通过产生乘积样条基函数对未知函数进行建模并自动识别变量之间的非线性相互作用。模型的具体形式如下：
表示基函数 对应的系数，而 的构造如下：
其中，， 代表变量的索引，对于每个基函数，变量集 是不重复的。
其模型训练和参数估计主要基于 reversible jump Metropolis-Hastings algorithm [6]。其中模型通过随机使用下列三个动作之一进行更新，包括 (a) 分割节点位置的变化；(b) 创建新的基函数；(c) 去除一个已有的基函数，然后再根据新的基函数，计算相关的接受概率。
除了更灵活的拟合能力，另一方面，我们还可以在 BMARS 构建的基函数中得到关于变量重要性的信息。具体来说，我们可以通过每个变量被包含在基函数中的次数来判断其重要性，次数越多则相对应的特征越重要。

集成学习 & 贝叶斯累加回归树

除了模型混合 (Model Mixing)，集成学习 (Ensemble Learning) [13] 提供了另一种组合模型的方法，并且取得了很好的效果。其主要想法是构造多个弱学习器然后将它们聚合成一个更强的学习器 [8]。在某些情况下，单个模型本身很难准确地勾画整个未知的复杂机制。因此，在集成学习框架中使用分而治之的方法是更好的策略，它允许每个模型拟合函数的一小部分。集成学习处理复杂数据的强大性能使其成为 BO 流程中概率模型的绝佳候选，然后其在 BO 的问题中并没有被充分研究。一个非常合适的贝叶斯集成学习模型便是贝叶斯累加回归树 (Bayesian Additive Regression Trees, BART) [4]，同时具有集成学习和贝叶斯框架的优势。它是一种基于树的模型，没有固有的关于平滑性的假设， 所以在对材料学中经常遇到的非平滑目标函数进行建模时，是一种更灵活的模型。
BART 也属于非参数回归模型，其具体公式如下：
其中 代表二叉回归树， 代表二叉回归树 的叶节点 的均值向量， 则是将 分配给 的函数。
BART 通过给各个二叉回归树添加限制其复杂度的先验，使得每棵树都比较小，只拟合未知函数的一小部分，从而实现集成学习的效果。模型的训练主要基于 [4] 提出的一个 backfitting MCMC algorithm，其抽样的效果很好所以已经被广泛使用。
同时，BART 也可以进行自动变量选择。每个自变量的重要性取决于其作为拆分规则的频率，频率越高则相应的因素越重要。[3] 进一步提出了一种基于排列的变量选择方法，这也是确定因子显著性的一个很好的选择。

当实验设计遇到 BO

通过在 BO 中使用 BMARS 或 BART 作为拟合黑盒函数的概率模型，[10] 提出了一个高效的自动实验设计的平台，旨在减少所需的试验次数以及寻找和发现最佳材料结构者的总费用。框架如图 4 所示。
图4：自动化实验设计框架的工作流程[10]
在这个工作流程中，实验者从一个初始已知数据集开始，样本大小可以很小（比如两个）。然后，实验者基于观察到的数据集训练 BMARS 或 BART 并从得到的后验分布中进行抽样。基于这些样本，计算每个可能的最佳材料结构（还未进行实验的材料结构）对应的采集函数。接着选择得分最高的候选者并进行相应的实验。得到新的结果后，将新数据加入已知数据集。如果此时循环停止的准则被满足，则结束循环，否则则基于扩充后的数据更新概率模型并指导下一轮的实验。
在这个完全自动化的框架中，我们需要提供是初始样本和循环停止标准。起始数据集可以是该项目之前的一些已经获得的相关数据。如果我们没有这种信息，我们可以随机进行少量实验来填充数据库并初始化框架中使用的概率模型。对于循环停止标准，材料学研究中常用的是找到满足要求的材料或者达到了最多能承受的实验次数（实验预算）[15]。

第一回合：测试函数上的初步对决

为了对比 BMARS，BART 和 GPR 的表现，[10] 选择了2个常用的测试函数（Rosenbrock 函数 和 Rastrigin 函数）作为第一回合的考题，目标是找到这两个函数的全局最小值点。如图 5 (a) 所示，Rosenbrock 函数有一个狭长的抛物线形平坦山谷，定位到这个山谷并不难，但是找到全局最优点就并不简单了。而图 5 (b) 则对应着Rastrigin 函数，我们可以看到函数曲面存在着众多的局部最小值点，这使得这个任务变得更加具有挑战性，因为算法很容易陷入这些局部最优值。
具体参与对比的 GPR 模型包括如下核函数：Radial-basis function (GP RBF)，RBK kernel with automatic relevance detection (GP RBK ARD) [2]，非平稳的 dot-product kernel (GP Dot) [16]，更灵活的 deep kernel (GP DKNet) [17]. 同时也和贝叶斯模型平均进行了比较 (BMA1 & BMA2) [15]。
图 6 展示了各模型具体的结果，曲线代表了各个每次迭代后以观测值中的最小值，所以曲线下降越快代表着模型表现越好。其中蓝色实线和红色实线分别代表 BMARS 和 BART 的结果。可以发现在第一回合，BMARS 优势明显，同时 BART 也展现出和许多基于 GPR 的模型旗鼓相当的表现。
图5：Plots of black-box functions [10].
图6：The average minimum observed based on each model in each iteration [10].

第二回合：MAX phases 中再次碰撞

为了进一步比较各个模型的效果，第二回合选择在了一个实际的材料发现的问题，如图 1 所示。这是一个前面简单提到过的任务，我们想要在 和 中找到性质最优的结构。其中 ， 和 分别可以在图 1 中红色，蓝色和黑色区域选择一种元素，剔除一些不稳定的组合后，可能具有最优性质的结构共计403种。图 7 展示了找到最优结构所需的平均试验次数，次数越少说明表现越好。我们可以看到，BART 和 BMARS 一直是表现的最好的几个模型。
图7
  • 更多的讨论和数据分析结果详见 [10]。

结语

随着数据科学的不断发展，越来越多之前 human-centric 的任务可以通过 data-centric 的算法得到更好的解决，从而促进相应学科的发展。这个过程其实也同时反过来促进了数据科学的发展。此文就是一个这样的例子，通过更好地从历史数据中提取有效信息，来设计更合理高效的实验计划。相信不少学科中还有许多这种可以和数据科学很好地结合同时相互促进的任务，它们都在等待着一双发现问题的眼睛。

参考文献

[1]  M. Aldeghi, F. H ̈ase, R. J. Hickman, I. Tamblyn, and A. Aspuru-Guzik.  Golem:  An algorithmfor robust experiment and process optimization, 2021.
[2]  S. A. Aye and P. Heyns.  An integrated Gaussian process regression for prediction of remaininguseful  life  of  slow  speed  bearings  based  on  acoustic  emission.Mechanical  Systems  and  SignalProcessing, 84:485–498, 2017.
[3]  J. Bleich, A. Kapelner, E. I. George, and S. T. Jensen. Variable selection for bart:  an applicationto gene regulation.The Annals of Applied Statistics, pages 1750–1781, 2014.
[4]  H. A. Chipman, E. I. George, R. E. McCulloch, et al.  Bart:  Bayesian additive regression trees.The Annals of Applied Statistics, 4(1):266–298, 2010. [5]  D.  G.  Denison,  B.  K.  Mallick,  and  A.  F.  Smith.   Bayesian  mars.Statistics  and  Computing,8(4):337–346, 1998.
[6]  P. J. Green.  Reversible jump markov chain monte carlo computation and Bayesian model deter-mination.Biometrika, 82(4):711–732, 1995.
[7]  F. Hase, M. Aldeghi, R. Hickman, L. Roch, E. Liles, M. Christensen, J. Hein, and A. Aspuru-Guzik.  Olympus:  a benchmarking framework for noisy optimization and experiment planning.Machine Learning:  Science and Technology, 2021.
[8]  B. Krawczyk, L. L. Minku, J. Gama, J. Stefanowski, and M. Wo ́zniak. Ensemble learning for datastream analysis:  A survey.Information Fusion, 37:132–156, 2017.
[9]  A.  G.  Kusne,  H.  Yu,  C.  Wu,  H.  Zhang,  J.  Hattrick-Simpers,  B.  DeCost,  S.  Sarker,  C.  Oses,C.  Toher,  S.  Curtarolo,  et  al.   On-the-fly  closed-loop  materials  discovery  via  bayesian  activelearning.Nature communications, 11(1):1–11, 2020.
[10]  B. Lei, T. Q. Kirk, A. Bhattacharya, D. Pati, X. Qian, R. Arroyave, and B. K. Mallick. Bayesianoptimization with adaptive surrogate models for automated experimental design.npj  Computa-tional Materials, 7(194), 2021.7
[11]  P.  Nikolaev,  D.  Hooper,  F.  Webber,  R.  Rao,  K.  Decker,  M.  Krein,  J.  Poleski,  R.  Barto,  andB. Maruyama.  Autonomy in materials research:  a case study in carbon nanotube growth.npjComputational Materials, 2(1):1–6, 2016.
[12]  C.  J.  Paciorek  and  M.  J.  Schervish.   Nonstationary  covariance  functions  for  gaussian  processregression.  InNIPS, pages 273–280. Citeseer, 2003.
[13]  O.  Sagi  and  L.  Rokach.   Ensemble  learning:  A  survey.Wiley  Interdisciplinary  Reviews:   DataMining and Knowledge Discovery, 8(4):e1249, 2018.
[14]  B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas.  Taking the human out ofthe loop:  A review of bayesian optimization.Proceedings of the IEEE, 104(1):148–175, 2015.
[15]  A.  Talapatra,  S.  Boluki,  T.  Duong,  X.  Qian,  E.  Dougherty,  and  R.  Arr ́oyave.   Autonomousefficient experiment design for materials discovery with bayesian model averaging.Physical ReviewMaterials, 2(11):113803, 2018.
[16]  C. K. Williams and C. E. Rasmussen.Gaussian processes for machine learning, volume 2.  MITpress Cambridge, MA, 2006.
[17]  A.  G.  Wilson,  Z.  Hu,  R.  Salakhutdinov,  and  E.  P.  Xing.   Deep  kernel  learning.   InArtificialintelligence and statistics, pages 370–378. PMLR, 2016.
[18]  Zengjian  Lv.Galileo-the  pioneer  of  experimental  science.http://html.rhhz.net/kjdb/20170818.htm, 2017.
继续阅读
阅读原文