贝叶斯优化实战(一)

529 阅读1小时+

贝叶斯优化实战(一)

原文:zh.annas-archive.org/md5/5ea9246d4fedff957a1c887275e38f01

译者:飞龙

协议:CC BY-NC-SA 4.0

第一章:开头内容

前言

随着我们在机器学习和相关领域面临的问题复杂性不断增加,优化我们对资源的利用和有效地做出知情决策变得越来越重要。贝叶斯优化是一种强大的技术,用于找到昂贵评估的目标函数的最大值和最小值,已经成为解决这一挑战的非常有用的解决方案。其中一个原因是该函数可以被视为黑匣子,这使得研究人员和从业者可以用贝叶斯推断作为主要优化方法来处理非常复杂的函数。

由于其复杂性,贝叶斯优化对于初学者机器学习从业者而言比其他方法更难以掌握。然而,像贝叶斯优化这样的工具必须是任何想要取得最佳结果的机器学习从业者的工具包中的一部分。要精通这个主题,必须对微积分和概率有非常扎实的直觉。

这就是贝叶斯优化实战的用武之地。在这本书中,Quan 美妙而成功地揭开了这些复杂概念的神秘面纱。通过实践方法、清晰的图表、真实世界的例子和有用的代码示例,他从理论和实践两个角度剖析了这个主题的复杂性。

Quan 利用他作为数据科学家和教育工作者的丰富经验,给读者提供了这些技术的非常清晰的图景以及它们如何应用于解决实际问题。从贝叶斯推断的原理开始,本书逐渐建立了贝叶斯优化和高斯过程模型的概念。它教授了最先进的库,如 GPyTorch 和 BoTorch,并探讨了它们在几个领域的应用。

这本书是任何数据科学或机器学习从业者必读的书籍,他们想要利用贝叶斯优化真正解决实际问题的潜力。我强烈推荐给任何想通过贝叶斯推断掌握优化艺术的人。

—路易斯·塞拉诺,博士,人工智能科学家和推广者,

机器学习的领悟者的作者

工程师和科学家面临着一个共同的挑战,这是捕捉其研究和创造力价值的关键所在。他们需要优化。一个机器学习工程师找到使模型泛化的超参数。一群物理学家调整自由电子激光以获得最大脉冲能量。一个软件工程师配置 JVM 的垃圾收集器以最大化服务器的吞吐量。一个材料科学工程师选择使太阳能电池的光吸收最大化的微观结构形态。在每个例子中,都有不能基于第一原理做出的设计决策。它们取决于实验评估。

要通过实验评估某物,人们可能需要执行软件、运行硬件或构建新物体,同时测量其性能。要找到一个好的设计,就需要进行多次评估。这些评估需要时间、花费金钱,可能还会带来风险。因此,必须尽可能少地进行实验评估以找到最优设计。这就是贝叶斯优化的全部内容。

在过去的 20 年中,我在我的工作中使用了贝叶斯优化和相关的先驱方法。在这段时间里,学术研究和工业应用的报告改善了贝叶斯优化的性能和扩展了其适用性。现在已经存在高质量的软件工具和具有项目特定优化器的技术。

可以将目前的状态类比为使用线性模型进行预测的状态。一个想要构建线性模型的工程师将发现,像 sklearn 这样的软件工具使他们能够设计各种类型(例如连续或分类)和数量的输入和输出变量,执行自动变量选择,并测量一般化的质量。同样地,一个想要构建贝叶斯优化器的工程师将发现,构建在 GPyTorch、pyro 和 PyTorch 之上的 BoTorch 提供了优化不同变量类型、最大化多个目标、处理约束和更多问题的工具。

本书从最基本的组件——高斯过程回归和获取函数的数值优化——开始教授贝叶斯优化,直到最新的处理大量评估(也称为观测)和异类设计空间的方法。在这个过程中,它覆盖了你可能需要的所有特殊化工具:处理约束、多目标、并行化评估和通过成对比较进行评估。你将找到足够的技术深度来让你对工具和方法感到舒适,并且足够的真实代码可以让你快速地将这些工具和方法用于实际工作。

尽管贝叶斯优化取得了很多成就,但鲜有面向新手的文献。这本书很好地填补了这个空缺。

——David Sweet,叶史瓦大学兼职教授

《面向工程师的实验》作者,Cogneato.xyz

前言

2019 年秋季,我是一名大一博士生,不确定应该在研究中解决哪个问题。我知道我想要专注于人工智能(AI)——用计算机自动化思考过程很有吸引力——但 AI 是一个庞大的领域,我很难将我的研究范围缩小到一个具体的主题上。

当我参加了一门名为《机器学习的贝叶斯方法》的课程时,所有的不确定性都消失了。在此之前,我在本科阶段曾简要接触过贝叶斯定理,但正是在这门课的第一节讲座中,一切开始变得清晰起来!贝叶斯定理提供了一种直观的思考概率的方式,对我来说,这是一种人类信念的优雅模型:我们每个人都有一个先验信念(关于任何事情),我们从这个先验开始,当我们观察到支持或反对该先验的证据时,我们的信念会更新,结果是反映先验和数据的后验信念。贝叶斯定理将这种保持信念的优雅方式引入到人工智能中,并在许多问题上找到应用,这对我来说是一个强烈的信号,表明贝叶斯机器学习是一个值得追求的主题。

当我们到达关于贝叶斯优化(BayesOpt)的讲座时,我的决定已经做出:理论是直观的,应用是多样的,可以建立的可能性非常大。再次,我内心的某种东西(至今仍然如此)被自动化思维或更具体地说是决策吸引。BayesOpt 正是这个完美的吸引力。我加入了罗曼·加内特(Roman Garnett)教授教授的研究实验室,我的 BayesOpt 之旅开始了!

跳到 2021 年,我花了一些时间研究和实施 BayesOpt 解决方案,我对 BayesOpt 的欣赏只增加了。我会向朋友和同事推荐使用它来处理困难的优化问题,并承诺 BayesOpt 会表现良好。只有一个问题:我找不到一个好的资源可以指向。研究论文数学内容繁重,在线教程太短,无法提供实质性的见解,而 BayesOpt 软件的教程杂乱无章,缺乏良好的叙述。

然后,一个想法浮现在脑海中,以托尼·莫里森(Toni Morrison)的话说,“如果有一本你想读的书,但它还没有被写出来,那么你必须写它。”多么真实啊!这个前景让我兴奋,原因有两个:我可以写一本关于我心爱的事物的书,而写作无疑会帮助我获得更深层次的洞察。我拼凑了一个提案,并联系了曼宁,这是我最喜欢的书籍的出版商,风格正是我所设想的。

2021 年 11 月,我的收购编辑安迪·沃尔德龙(Andy Waldron)给我发了一封电子邮件,标志着曼宁(Manning)的第一次沟通。 2021 年 12 月,我签署了合同并开始写作,这比我最初想象的时间要长(我相信每本书都是如此)。 2023 年 4 月,在出版前的最后几个步骤之一,我写了这篇前言!

致谢

抚养一个孩子需要整个村庄的参与,写一本书也是如此。以下只是我自己村庄的一小部分,在写作过程中给予了我巨大帮助的人。

我首先要感谢我的父母 Bang 和 Lan,他们的持续支持使我能够毫无畏惧地探索未知:出国留学;攻读博士学位;当然,还有写书。我还要诚挚地感谢我的姐姐和知己 Nhu,她总是在我最艰难的时刻帮助我。

贝叶斯优化是我博士研究的重要组成部分,我要感谢那些在项目中真正让我的博士经历变得宝贵的人。特别感谢我的指导老师罗曼·加内特,他毫不费力地说服我去从事贝叶斯机器学习的研究。是你开启了这一切。我还要感谢来自主动学习实验室的朋友们:叶虎·陈、沙扬·莫纳德杰米和阿尔维塔·奥特利教授。他们说博士阶段的回报很少,而与你们一起工作正是构成了其中大部分回报。

接下来,我要由衷感谢曼宁公司的出色团队。我感谢我的开发编辑 Marina Michaels,她以最高水平的专业精神、关怀、支持和耐心领导着这艘船,从一开始就与我同舟共济。能够与你配对完成我们的项目,我感到非常幸运。感谢我的收购编辑 Andy Waldron,即使已经有一个更好的作者在写一本类似主题的书,他仍对这个想法充满信心,以及 Ivan Martinovic´,他帮助我解决了 AsciiDoc 的问题,并耐心修复了我的标记代码。

我要感谢审稿人们投入时间和精力,大大提高了写作质量:Allan Makura、Andrei Paleyes、Carlos Aya-Moreno、Claudiu Schiller、Cosimo Attanasi、Denis Lapchev、Gary Bake、George Onofrei、Howard Bandy、Ioannis Atsonios、Jesús Antonino Juárez Guerrero、Josh McAdams、Kweku Reginald Wade、Kyle Peterson、Lokesh Kumar、Lucian Mircea Sasu、Marc-Anthony Taylor、Marcio Nicolau、Max Dehaut、Maxim Volgin、Michele Di Pede、Mirerfan Gheibi、Nick Decroos、Nick Vazquez、Or Golan、Peter Henstock、Philip Weiss、Ravi Kiran Bamidi、Richard Tobias、Rohit Goswami、Sergio Govoni、Shabie Iqbal、Shreesha Jagadeesh、Simone Sguazza、Sriram Macharla、Szymon Harabasz、Thomas Forys 和 Vlad Navitski。

写书时不可避免地会有盲点,而审稿人们帮助填补了这些盲点,并让作者专注于真正重要的事情。衷心感谢 Kerry Koitzsch 提供的有见地的反馈和 James Byleckie 在代码和写作方面提出的优秀建议。

最后,我要感谢 GPyTorch 和 BoTorch 库背后的团队,这些库是为本书开发的代码的主要工作马力。我尝试过各种高斯过程和贝叶斯优化的库,但总是发现自己回到 GPyTorch 和 BoTorch。我希望本书能在这些库周围构建一个已经很棒的社区中发挥作用。

关于本书

过去要了解贝叶斯优化,人们需要搜索相关库的在线文章和教程,这些文章和教程零散分布,由于其性质,不会深入探讨具体细节。你也可以求助于技术教材,但它们通常太过密集和数学密集,这对于希望立即着手实践的从业者来说是一个挑战。

本书通过提供一系列实践讨论、对感兴趣的读者的更深入材料的引用以及可直接使用的代码示例来填补这一空白。它首先为贝叶斯优化的组成部分建立直觉,然后使用最先进的软件在 Python 中实现它们。

本书的精神是提供一个以数学和概率的高层次直觉为基础的贝叶斯优化易于理解的介绍。感兴趣的读者可以在全书中找到更多被引用的更深入的技术文本,以深入了解感兴趣的主题。

谁应该阅读本书?

对于对超参数调优、A/B 测试或实验以及更一般的决策制定感兴趣的数据科学家和 ML 从业者,本书将有所裨益。

化学、材料科学和物理等科学领域的研究人员面临着困难的优化问题,也会发现本书有所帮助。尽管大多数跟随内容所必需的背景知识将会被涵盖,但读者应该熟悉 ML 中的常见概念,例如训练数据、预测模型、多元正态分布等。

本书的组织方式:一条路线图

本书包括四个主要部分。每个部分都包含涵盖相应主题的几章:

  • 第一章介绍了使用真实用例的贝叶斯优化。它还包括了一个视觉示例,展示了贝叶斯优化如何加速寻找昂贵函数的全局最优解,而不涉及技术细节。

第一部分涵盖了高斯过程作为我们希望优化的函数的预测模型。其核心论点是,高斯过程提供了校准的不确定性量化,这在我们的贝叶斯优化框架中是必不可少的。本部分由两章组成:

  • 第二章显示高斯过程是从一些观察数据中学习回归模型的自然解决方案。高斯过程定义了函数的分布,并且可以根据一些观察数据来更新以反映我们对函数值的信念。

  • 第三章介绍了我们将先验信息合并到高斯过程中的两种主要方式:均值函数和协方差函数。均值函数指定了一般趋势,而协方差函数指定了函数的平滑程度。

第二部分列举了贝叶斯优化策略,这些策略是关于如何进行函数评估的决策过程,以便尽可能高效地确定全局最优值。虽然不同的策略由不同的目标驱动,但它们都旨在平衡探索和利用之间的权衡。这部分由三章组成:

  • 第四章讨论了一种自然的方法来决定哪个函数评估是最有利的:考虑从当前最佳函数值中获得的改进。由于基于高斯过程的对函数的信念,我们可以计算这些与改进相关的数量,并且可以在封闭形式下廉价地进行,从而实现了两种特定的贝叶斯优化策略:改进的概率和预期改进。

  • 第五章探讨了贝叶斯优化与另一种常见的问题类别之间的联系,称为多臂老虎机。我们学习如何将多臂老虎机策略转换为贝叶斯优化设置,并获得相应的策略:上置信界和汤普森采样。

  • 第六章考虑了一种减少我们对函数全局最优值信念中不确定性的策略。这构成了基于熵的策略,使用了一种称为信息论的数学子领域。

第三部分介绍了一些最常见的用例,这些用例不能很好地适应本书迄今为止开发的工作流程,并展示了如何修改贝叶斯优化来解决这些优化任务:

  • 第七章介绍了批量优化,在这种情况下,为了提高吞吐量,我们允许实验并行运行。例如,可以同时在一组 GPU 上并行训练大型神经网络的多个实例。这需要优化策略同时返回多个建议。

  • 第八章讨论了安全关键的用例,我们在这些情况下不能自由地探索搜索空间,因为一些函数评估可能会产生不利影响。这促使了这样一种设置,即对于所讨论的函数如何行为有一定的约束,并且我们需要在优化策略的设计中考虑这些约束。

  • 第九章表明,当我们可以以不同成本和精度级别观察函数值的多种方式时——通常称为多信度贝叶斯优化——考虑可变成本可以提高优化性能。

  • 第十章涵盖了成对比较,已经显示出比数字评估或评级更准确地反映了个人偏好,因为它们更简单,对标注者的认知负荷较轻。第十章将贝叶斯优化扩展到此设置,首先使用特殊的高斯过程模型,然后修改现有策略以适应这个成对比较的工作流程。

  • 一个人可能希望同时优化多个可能冲突的目标。第十一章研究了这个多目标优化问题,并展示了贝叶斯优化如何在这种情况下扩展。

第四部分涉及高斯过程模型的特殊变体,展示了它们在建模和提供不确定性校准预测方面的灵活性和效果,即使在贝叶斯优化环境之外也是如此:

  • 在第十二章中,我们了解到在某些情况下,无法获得经过训练的高斯过程的闭合形式解。然而,可以使用复杂的近似策略进行高保真度的近似。

  • 第十三章展示了由于 Torch 生态系统的存在,将 PyTorch 神经网络与 GPyTorch 高斯过程结合起来是一个无缝过程。这使得我们的高斯过程模型更加灵活和表达能力更强。

初学者将从前六章获益良多。有经验的从业者希望将贝叶斯优化应用于他们的案例中,可能会在第七章到第十一章中找到价值,这些章节可以独立阅读,任意顺序。长期使用高斯过程的用户很可能对最后两章感兴趣,我们在这些章节中开发了专门的高斯过程模型。

关于代码

您可以从本书的在线版本 livebook.manning.com/book/bayesian-optimization-in-action 获取可执行的代码片段。代码可以从 Manning 网站 www.manning.com/books/bayesian-optimization-in-action 和 GitHub github.com/KrisNguyen135/bayesian-optimization-in-action 下载;后者接受问题和拉取请求。

您将使用 Jupyter notebooks 来运行附带书籍的代码。Jupyter notebooks 提供了一种干净的方式来动态地使用代码,使我们能够探索每个对象的行为以及它与其他对象的交互。有关如何开始使用 Jupyter notebook 的更多信息,请查找它们的官方网站:jupyter.org。在我们的情况下,动态探索对象的能力特别有帮助,因为 Bayesian 优化工作流的许多组件是由 GPyTorch 和 BoTorch 实现的 Python 对象。

GPyTorch 和 BoTorch 是 Python 中用于高斯过程建模和贝叶斯优化的首选库。还有其他选择,例如 scikit-Learn 的 scikit-optimize 扩展或 GPflow 和 GPflowOpt,它们扩展了 TensorFlow 框架用于贝叶斯优化。然而,GPyTorch 和 BoTorch 的组合构成了最全面和灵活的代码库,其中包括许多来自贝叶斯优化研究的最新算法。根据我自己使用贝叶斯优化软件的经验,我发现 GPyTorch 和 BoTorch 在易于初学者使用和提供最新方法之间取得了良好的平衡。

有一件事需要注意:正是因为这些库正在积极地维护,书中展示的 API 可能在新版本中略有变化,因此重要的是您按照requirements.txt文件中指定的库版本来运行代码,以避免错误。你可以在官方 Python 文档中找到更多关于如何使用requirements.txt文件创建 Python 环境的说明,例如在packaging.python.org/en/latest/guides/installing-using-pip-and-virtual-environments。话虽如此,要使用新版本,您可能只需要对代码做些微小的修改。

当您阅读本书时,您会注意到文本往往只专注于代码的关键组件,省略了许多细节,比如库的导入和繁琐的代码。(当然,代码第一次使用时,它会在文本中被正确介绍。)简洁的讨论有助于我们专注于每一章中真正新颖的内容,避免重复。另一方面,Jupyter 笔记本中的代码是自包含的,可以单独运行每个笔记本,无需任何修改。

liveBook 讨论论坛

购买《Bayesian Optimization in Action》包括免费访问 liveBook,Manning 的在线阅读平台。使用 liveBook 独有的讨论功能,您可以将评论附加到全书或特定的部分或段落上。您可以轻松地为自己做笔记,提出和回答技术问题,并从作者和其他用户那里获得帮助。要访问论坛,请转到livebook.manning.com/book/bayesian-optimization-in-action/discussion。您还可以在livebook.manning.com/discussion了解更多关于 Manning 论坛和行为规则的信息。

Manning 对我们的读者的承诺是提供一个有意义的对话场所,既是个人读者之间的对话,也是读者与作者之间的对话。这并不是对作者参与的具体数量的承诺,他们对论坛的贡献仍然是自愿的(且未支付的)。我们建议您尝试向作者提出一些具有挑战性的问题,以免他们的兴趣消失!论坛和以前讨论的存档将在该书印刷期间都可以从出版商的网站上访问到。

关于作者

Quan Nguyen 是一位 Python 程序员和机器学习爱好者。他对解决涉及不确定性的决策问题感兴趣。Quan 撰写了几本关于 Python 编程和科学计算的书籍。他目前正在华盛顿大学圣路易斯分校攻读计算机科学博士学位,他在那里研究机器学习中的贝叶斯方法。

关于技术编辑

本书的技术编辑是 Kerry Koitzsch。Kerry 是一位作者,软件架构师,在企业应用程序和信息架构解决方案的实施方面拥有三十多年的丰富经验。Kerry 是一本关于分布式处理的书籍的作者,以及许多较短的技术出版物的作者,并拥有一项关于创新 OCR 技术的专利。他还是美国陆军成就奖的获得者。

关于封面插图

Bayesian Optimization in Action 封面上的图案标题为“波兰人”,摘自 Jacques Grasset de Saint-Sauveur 的一本 1797 年出版的作品集。每一幅插图都是手工精细绘制和上色的。

在那些日子里,通过人们的服装很容易辨别他们住在哪里,以及他们的行业或生活地位。Manning 凭借基于数个世纪前地区文化的丰富多样性的书籍封面,赞美了计算机业务的创造性和主动性,这些书籍封面是由这类集合中的图片重新唤起的。

第二章:贝叶斯优化简介

本章内容包括

  • 是什么促使了贝叶斯优化以及它是如何工作的

  • 贝叶斯优化问题的实际例子

  • 贝叶斯优化的一个玩具示例

你选择阅读本书是一个很棒的选择,我对你即将开始的旅程感到兴奋!从高层次来看,贝叶斯优化是一种优化技术,当我们试图优化的函数(或者一般情况下,当输入一个值时产生输出的过程)是一个黑盒且评估起来时间、金钱或其他资源成本很高时,可以应用此技术。这个设置涵盖了许多重要的任务,包括超参数调优,我们将很快定义它。使用贝叶斯优化可以加速搜索过程,并帮助我们尽快找到函数的最优解。

尽管贝叶斯优化在机器学习研究界一直受到持久的关注,但在实践中,它并不像其他机器学习话题那样常用或广为人知。但为什么呢?有些人可能会说贝叶斯优化具有陡峭的学习曲线:使用者需要理解微积分、使用概率,并且需要是一个经验丰富的机器学习研究者,才能在应用中使用贝叶斯优化。本书的目标是打破贝叶斯优化难以使用的观念,并展示该技术比想象的更直观、更易用。

在本书中,我们会遇到许多插图、图表和代码,旨在使讨论的主题更加简单明了和具体。你将了解贝叶斯优化的每个组成部分在高层次上是如何工作的,并学会如何使用 Python 中的最先进的库来实现它们。配套的代码还可帮助你快速上手你自己的项目,因为贝叶斯优化框架非常通用和“即插即用”。这些练习对此也非

总的来说,我希望这本书对你的机器学习需求有所帮助,并且是一本有趣的阅读。在我们深入讨论实际内容之前,让我们花点时间来讨论贝叶斯优化试图解决的问题。

1.1 寻找一个昂贵黑盒函数的最优解

如前所述,超参数调优是贝叶斯优化在机器学习中最常见的应用之一。我们在本节中探讨了这个问题以及其他一些问题,作为黑盒优化问题的一个例子。这将帮助我们理解为什么需要贝叶斯优化。

1.1.1 超参数调优作为昂贵黑盒优化问题的一个示例

假设我们想在一个大数据集上训练神经网络,但我们不确定这个神经网络应该有多少层。我们知道神经网络的架构是深度学习中的一个成功因素,因此我们进行了一些初步测试,并得到了表格 1.1 中显示的结果。

表格 1.1 超参数调优任务的示例

层数测试集准确率
50.72
100.81
200.75

我们的任务是决定神经网络在寻找最高准确率时应该有多少层。很难决定我们应该尝试下一个数字是多少。我们找到的最佳准确率为 81%,虽然不错,但我们认为通过不同数量的层,我们可以做得更好。不幸的是,老板已经设定了完成模型实施的截止日期。由于在我们的大型数据集上训练神经网络需要几天时间,我们只剩下几次试验的机会,然后就需要决定我们的网络应该有多少层。考虑到这一点,我们想知道我们应该尝试哪些其他值,以便找到提供最高可能准确率的层数。

这项任务旨在找到最佳设置(超参数值),以优化模型的某些性能指标,如预测准确率,在机器学习中通常被称为超参数调整。在我们的示例中,神经网络的超参数是其深度(层数)。如果我们使用决策树,常见的超参数包括最大深度、每个节点的最小数据点数和分裂标准。对于支持向量机,我们可以调整正则化项和核函数。由于模型的性能很大程度上取决于其超参数,超参数调整是任何机器学习流水线的重要组成部分。

如果这是一个典型的真实世界数据集,这个过程可能需要大量的时间和资源。来自 OpenAI 的图 1.1(openai.com/blog/ai-and-compute/)显示,随着神经网络变得越来越大和越来越深,所需的计算量(以 petaflop/s-days 为单位)呈指数增长。

图 1.1 训练大型神经网络的计算成本一直在稳步增长,使得超参数调整变得越来越困难。

这意味着在大型数据集上训练模型是相当复杂的,并且需要大量的工作。此外,我们希望确定能够提供最佳准确率的超参数值,因此需要进行多次训练。我们应该如何选择数值来对我们的模型进行参数化,以便尽快找到最佳组合?这是超参数调整的核心问题。

回到我们在第 1.1 节中的神经网络示例,我们应该尝试多少层才能找到高于 81%的准确度?在 10 层和 20 层之间的某个数值是有前途的,因为在 10 层和 20 层,我们的性能比在 5 层时更好。但我们应该检查哪个确切的数值仍然不明显,因为在 10 和 20 之间的数值仍然可能有很大变异性。当我们说变异性时,我们隐含地谈论了我们关于模型测试准确性如何随层数变化而变化的不确定性。即使我们知道 10 层导致 81%,20 层导致 75%,我们仍然不能确定例如 15 层会产生什么值。这就是说,当我们考虑 10 和 20 之间的这些值时,我们需要考虑我们的不确定性水平。

此外,如果某个大于 20 的数值为我们提供了最高可能的准确度怎么办?这对于许多大型数据集来说是一种情况,其中足够的深度对于神经网络学习任何有用的东西都是必要的。或者,尽管不太可能,少于 5 层的小层数实际上是我们需要的吗?

我们应该如何以有原则的方式探索这些不同的选择,以便在时间耗尽和我们必须向老板汇报时,我们可以充分自信地认为我们已经找到了我们模型的最佳层数?该问题是昂贵黑盒优化问题的一个例子,我们接下来会讨论这个问题。

1.1.2 昂贵黑盒优化问题

在这个子章节中,我们正式介绍了昂贵黑盒优化问题,这是贝叶斯优化的目标。理解为什么这个问题很难有助于我们理解,为什么贝叶斯优化比更简单的、更天真的方法更受欢迎,比如网格搜索(我们将搜索空间分为相等的段)或随机搜索(我们使用随机性来指导我们的搜索)。

在这个问题中,我们可以黑匣子方式访问函数(一些输入-输出机制),我们的任务是找到最大化此函数输出的输入。该函数通常称为目标函数,因为优化它是我们的目标,并且我们希望找到目标函数的最优解——产生最高函数值的输入。

目标函数的特点

术语黑盒意味着我们不知道目标的底层公式;我们唯一能够访问的是通过在某个输入处计算函数值进行观察时得到的函数输出。在我们的神经网络示例中,我们不知道如果我们逐层增加层数,我们的模型的准确性将如何变化(否则,我们将只选择最佳层)。

这个问题很昂贵,因为在许多情况下,做出观察(在某个位置评估目标)的成本非常高昂,使得像穷举搜索这样的天真方法难以处理。在机器学习和尤其是深度学习中,时间通常是主要的约束条件,正如我们之前讨论过的那样。

超参数调整属于这一类昂贵的黑盒优化问题,但不是唯一的!任何试图找到一些设置或参数来优化一个过程,而不知道不同设置如何影响和控制过程结果的程序都属于黑盒优化问题。此外,尝试特定设置并观察其对目标过程(目标函数)的结果是耗时的、昂贵的或在某种程度上成本高昂的。

定义 尝试特定设置的行为——即,在某个输入处评估目标函数的值——称为发出查询查询目标函数。整个过程总结如图 1.2。

图 1.2 黑盒优化问题的框架。我们反复查询不同位置的函数值以找到全局最优解。

1.1.3 其他昂贵的黑盒优化问题的实际例子

现在,让我们考虑一些属于昂贵的黑盒优化问题类别的实际例子。我们会发现这样的问题在这个领域中很常见;我们经常会遇到想要优化但只能评估少数次的函数。在这些情况下,我们希望找到一种方法来智能地选择在哪里评估函数。

第一个例子是药物发现——科学家和化学家识别具有理想化学特性的化合物,可以合成成药物的过程。正如你可以想象的那样,实验过程非常复杂并且成本很高。使这项药物发现任务令人生畏的另一个因素是近年来已经观察到的药物发现研发生产力下降趋势。这种现象被称为艾尔姆定律——摩尔定律的反转——它大致说明了每十亿美元批准的新药物数量在固定时间内减半。艾尔姆定律在杰克·W·斯坎内尔(Jack W. Scannell)、亚历克斯·布兰克利(Alex Blanckley)、海伦·波尔登(Helen Boldon)和布莱恩·沃灵顿(Brian Warrington)撰写的自然杂志论文“诊断药物研发效率下降”(www.nature.com/articles/nrd3681)的第 1 张图片中有可视化。(或者,你也可以简单地在谷歌上搜索“艾尔姆定律”的图像。)

赫尔姆斯定律显示,每十亿美元的药物研究与开发(R&D)投资所得到的药物发现能力,经过对数尺度上的线性下降。换句话说,对于固定数量的 R&D 投资,药物发现能力在最近几年呈指数级下降。虽然在不同年份存在局部趋势的起伏,但从 1950 年到 2020 年,指数下降趋势是明显的。

实际上,同样的问题适用于任何科学发现任务,其中科学家们通过使用需要顶级设备和可能需要几天或几周才能完成的实验,搜索罕见、新颖和有用的新化学品、材料或设计,用于某种度量标准。换句话说,他们试图针对极为昂贵的数据评估来优化它们各自的目标函数。

以表 1.2 为例,展示了真实数据集中的几个数据点。目标是找到具有最低混合温度的合金组成(来自四个母元素)。这是一个黑盒优化问题。在这里,材料科学家们研究了铅(Pb)、锡(Sn)、锗(Ge)和锰(Mn)的合金组成。每个给定的组合百分比对应于一个可以在实验室中合成和实验的潜在合金。

表 1.2 来自材料发现任务的数据

Pb 的百分比Sn 的百分比Ge 的百分比Mn 的百分比混合温度(°F)
0.500.500.000.00192.08
0.330.330.330.00258.30
0.000.500.500.00187.24
0.000.330.330.33188.54
来源: 作者的研究工作。

由于低的混合温度表示合金结构稳定、有价值,目标是找到混合温度尽可能低的组成。但是有一个瓶颈:通常需要数天时间才能确定给定合金的混合温度。我们要算法地解决的问题是类似的:给定我们看到的数据集,我们应该尝试下一个组合(在铅、锡、锗和锰的含量上如何)以找到最低的混合温度?

另一个例子是在采矿和石油钻探中,或者更具体地说,在一个大区域内找到具有最高价值矿物或石油产量的地区。这需要大量的规划、投资和劳动力,是一项昂贵的事业。由于挖掘作业对环境有重大负面影响,因此有相应的法规来减少采矿活动,在这个优化问题中对可以进行的函数评估数量设定了限制。

昂贵的黑盒优化中的中心问题是:如何决定在哪里评估这个目标函数,以便在搜索结束时找到其最优值?正如我们在后面的例子中看到的,简单的启发式方法——如随机或网格搜索,这些方法是流行的 Python 包(如 scikit-learn)实现的——可能会导致对目标函数的浪费性评估,从而导致整体优化性能较差。这就是贝叶斯优化发挥作用的地方。

1.2 介绍贝叶斯优化

考虑到昂贵的黑盒优化问题,现在我们介绍贝叶斯优化作为这个问题的解决方案。这给了我们一个高层次的关于贝叶斯优化是什么以及它如何使用概率机器学习来优化昂贵的黑盒函数的概念。

贝叶斯优化(BayesOpt)的定义是一种机器学习技术,它同时维护一个预测模型来学习关于目标函数的信息并且通过贝叶斯概率和决策理论决定如何获取新数据来完善我们对目标的知识。

通过数据,我们指的是输入输出对,每个对应一个输入值到该输入处的目标函数值的映射。在超参数调优的具体案例中,这些数据与我们旨在调整的机器学习模型的训练数据不同。

在贝叶斯优化过程中,我们根据贝叶斯优化算法的建议做出决定。一旦我们采取了贝叶斯优化建议的行动,贝叶斯优化模型将根据该行动的结果进行更新,并继续推荐下一步要采取的行动。这个过程重复进行,直到我们有信心找到了最优行动。

这个工作流程有两个主要组成部分:

  • 一个从我们所做的观察中学习并对未见数据点上的目标函数值进行预测的机器学习模型

  • 通过评估目标以定位最优值的优化策略

我们在以下小节中介绍每个组件。

1.2.1 用高斯过程建模

贝叶斯优化首先在我们试图优化的目标函数上拟合一个预测的机器学习模型——有时,这被称为替代模型,因为它充当我们从观察中相信的函数和函数本身之间的替代。这个预测模型的作用非常重要,因为它的预测结果会影响贝叶斯优化算法的决策,并且直接影响优化性能。

几乎在所有情况下,高斯过程(GP)被用于这种预测模型的角色,我们在本小节中对此进行了研究。在高层次上,与任何其他机器学习模型一样,高斯过程(GP)的运行原则是相似的数据点产生相似的预测。与岭回归、决策树、支持向量机或神经网络等其他模型相比,GPs 可能不是最受欢迎的模型类别。然而,正如我们在本书中一再看到的那样,GPs 带有一个独特且至关重要的特性:它们不像其他提到的模型那样产生点估计预测;相反,它们的预测是以概率分布的形式。以概率分布或概率预测的形式进行的预测在贝叶斯优化中是关键的,它使我们能够量化我们的预测不确定性,进而改善我们在做出决策时的风险-回报折衷。

首先让我们看看当我们在数据集上训练一个 GP 时它是什么样子的。比如说,我们有兴趣训练一个模型来从表 1.3 的数据集中学习,该数据集在图 1.3 中被可视化为黑色的 x

表 1.3 对应于图 1.3 的一个示例回归数据集

训练数据点标签
1.14701.8423
-4.07120.7354
0.96270.9627
1.24711.9859

图 1.3 非贝叶斯模型,比如岭回归器,做出的是点估计,而高斯过程则产生概率分布作为预测。因此,高斯过程提供了一个校准的不确定性量化,这是在做出高风险决策时的一个重要因素。

我们首先在这个数据集上拟合了一个岭回归模型,并在 -5 和 5 的范围内做出预测;图 1.3 的顶部面板显示了这些预测。岭回归模型是线性回归模型的改进版本,其中模型的权重被正则化,以便更偏爱较小的值,以防止过拟合。该模型在给定测试点时所做的每个预测都是一个单值数字,它并没有捕获我们对所学习函数行为的不确定性水平的认识。例如,给定 2 的测试点,该模型简单地预测为 2.2。

我们不需要过多地了解这个模型的内部工作原理。关键在于岭回归器产生没有不确定性度量的点估计,这也是许多机器学习模型的情况,比如支持向量机、决策树和神经网络。

那么,高斯过程是如何做出预测的呢?如图 1.3 的底部面板所示,高斯过程的预测是以概率分布的形式(具体来说,是正态分布)进行的。这意味着在每个测试点,我们有一个平均预测(实线)以及所谓的 95% 置信区间(CI)(阴影区域)。

需要注意的是,“CI”这个缩写词在统计学的频率主义中常用来缩写“置信区间”(confidence interval);在本书中,我只使用“CI”来指代“可信区间”(credible interval)。虽然这两个概念在技术上有许多不同之处,但从高层次上来看,我们仍然可以将本书中的 CI 视为一个区间,这个区间内有可能包含一个感兴趣的数量(在这种情况下,就是我们正在预测的函数的真实值)。

GP vs. 岭回归

有趣的是,当使用相同的协方差函数(也称为“核”)时,“高斯过程”(GP)和“岭回归模型”产生了相同的预测结果(对于 GP 来说是均值预测),如图 1.3 所示。我们会在第三章更深入地讨论协方差函数。这意味着,GP 具有岭回归模型所有的好处,同时还提供了额外的 CI 预测。

在实际测试位置上,这个 CI 有效地度量了我们对每个测试位置的价值的不确定性水平。如果一个位置的预测 CI 较大(比如图 1.3 中的-2 或 4),则这个值的可能值范围更广。换句话说,我们对这个值的确定性更低。如果一个位置的 CI 较小(图 1.3 中的 0 或 2),则我们对这个位置的值更有信心。GP 的一个很好的特点是,对于训练数据的每个点,预测 CI 接近于 0,这表示我们对其值没有任何不确定性。这是有道理的;毕竟,我们已经从训练集中知道了该值。

带噪音的函数评估

虽然在图 1.3 中不是这种情况,但是在我们的数据集中,数据点的标签可能是有噪声的。在实际世界中,观察数据的过程很可能会受到噪声的干扰。在这种情况下,我们可以使用 GP 进一步指定噪声水平,观察数据点的 CI 将不会降为 0,而是降至指定的噪声水平。这表明了使用 GP 建模所具有的灵活性。

能够将我们对不确定性的水平进行量化的能力(称为“不确定性量化”)在任何高风险的决策过程中都非常有用,比如贝叶斯优化。在 1.1 节中出现的情景再次设想一下,我们调整神经网络中的层数,并且只有时间尝试一个更多的模型。假设在那些数据点上训练之后,GP 预测 25 层的平均精度将为 0.85,相应的 95% CI 为 0.81 至 0.89。另一方面,对于 15 层,GP 预测我们的精度平均也是 0.85,但是 95% CI 为 0.84 至 0.86。在这种情况下,即使这两个数字具有相同的期望值,选择 15 层是相当合理的,因为我们更“确定”15 层将给我们带来好的结果。

清楚地说,GP 不会为我们做出任何决定,但它确实通过其概率预测为我们提供了一种方法。决策留给 BayesOpt 框架的第二部分:策略。

使用 BayesOpt 策略做决策

除了作为预测模型的 GP 之外,在 BayesOpt 中,我们还需要一个决策过程,我们将在本小节中探讨这个问题。这是 BayesOpt 的第二个组成部分,它接受 GP 模型所做的预测,并推理如何最好地评估目标函数,以便有效地找到最优解。

如前所述,95% CI 为 0.84 至 0.86 的预测要比 95% CI 为 0.81 至 0.89 的预测更好,特别是如果我们只有一次尝试的机会。这是因为前者更像是一件确定的事情,保证为我们带来一个好结果。在两个点的预测均值和预测不确定性可能不同的更一般情况下,我们应该如何做出这个决定?

这正是 BayesOpt 策略帮助我们做的事情:量化一个点的有用性,考虑到其预测概率分布。策略的工作是接受 GP 模型,该模型代表我们对目标函数的信念,并为每个数据点分配一个分数,表示该点在帮助我们识别全局最优解方面的帮助程度。这个分数有时被称为获取分数。我们的工作是选择最大化这个获取分数的点,并在该点评估目标函数。

我们在图 1.4 中看到与图 1.3 中相同的 GP,在底部面板中显示了一个名为Expected Improvement的特定 BayesOpt 策略如何在x-轴上的每个点(在我们的搜索空间内的-5 到 5 之间)得分。我们将在第四章学习这个名称的含义以及该策略如何对数据点进行评分。现在,让我们先记住,如果一个点具有较大的获取分数,这个点对于定位全局最优解是有价值的。

图 1.4 BayesOpt 策略通过其在定位全局最优解中的有用性对每个单独的数据点进行评分。该策略倾向于高预测值(其中回报更有可能)以及高不确定性(其中回报可能较大)。

在图 1.4 中,最佳点在 1.8 左右,这是有道理的,因为根据我们在顶部面板中的 GP,在那里我们也实现了最高的预测均值。这意味着我们将选择在 1.8 处评估我们的目标,希望从我们收集到的最高值中得到改进。

我们应该注意到,这不是一个一次性的过程,而是一个学习循环。在循环的每一次迭代中,我们根据我们从目标中观察到的数据训练一个高斯过程(GP),在这个高斯过程上运行贝叶斯优化策略,以得到一个希望帮助我们确定全局最优的建议,然后在推荐位置进行观察,将新点添加到我们的训练数据中,并重复整个过程,直到达到某个终止条件。事情可能变得有点混乱,所以是时候退后一步,看看贝叶斯优化的大局了。

与实验设计的联系

此时,贝叶斯优化的描述可能让你想起了统计学中的实验设计(DoE)的概念,它旨在通过调整可控设置来解决优化目标函数的问题。这两种技术之间存在许多联系,但是贝叶斯优化可以被看作是一种更一般的方法,它由机器学习模型高斯过程(GP)驱动。

1.2.3 组合高斯过程和优化策略形成优化循环

在本小节中,我们总结了我们迄今为止所描述的内容,并使过程更加具体。我们全面地看到了贝叶斯优化的工作流程,并更好地理解了各个组成部分是如何相互配合的。

我们从一个初始数据集开始,就像表 1.1、1.2 和 1.3 中的那样。然后,贝叶斯优化的工作流程在图 1.5 中进行了可视化,总结如下:

  1. 我们在这个数据集上训练了一个高斯过程(GP)模型,根据我们从训练数据中观察到的内容给出了对我们的目标在每个地方的信念。这种信念由实线和阴影区域表示,就像图 1.3 和 1.4 中的那样。

  2. 然后,贝叶斯优化策略接收这个高斯过程,并根据该点在域中的价值对每个点进行评分,这如图 1.4 中的下曲线所示。

  3. 最大化该分数的点是我们选择下一个要评估目标的点,然后将其添加到我们的训练数据集中。

  4. 这个过程会重复进行,直到我们无法再评估目标。

BayesOpt 循环

图 1.5 贝叶斯优化循环,结合了高斯过程(GP)建模和决策制定的策略。现在可以使用这个完整的工作流程来优化黑盒函数。

与监督学习任务不同,我们只需在训练数据集上拟合一个预测模型并在测试集上进行预测(只包括步骤 1 和 2),贝叶斯优化工作流程通常被称为主动学习。主动学习是机器学习中的一个子领域,我们可以决定我们的模型从哪些数据点中学习,而这个决策过程则由模型本身来决定。

正如我们所说,GP 和政策是这个 BayesOpt 过程的两个主要组成部分。如果 GP 没有很好地对客观函数进行建模,那么我们将无法很好地将训练数据中的信息通知给政策。另一方面,如果政策不能很好地给“好”点分配高分和给“坏”点分配低分(其中意味着有助于找到全局最优解),那么我们的后续决策将是错误的,并且很可能会取得糟糕的结果。

换句话说,如果没有一个好的预测模型,比如一个 GP,我们就无法通过校准的不确定性做出良好的预测。没有政策,我们可以做出良好的预测,但我们不会做出良好的决策

我们在本书中多次考虑的一个例子是天气预报。想象一下这样一个情景,你要决定在离开家去上班前是否带伞,并查看手机上的天气预报应用程序。

不用说,应用程序的预测需要准确可靠,这样你就可以自信地根据它们做出决定。一个总是预测晴天的应用程序是不够的。此外,你需要一种明智的方式来根据这些预测做出决策。无论多么可能下雨,永远不带伞都是一个糟糕的决策政策,当下雨时会让你陷入麻烦。另一方面,即使天气预报有 100% 的晴天可能性,也不是一个明智的决定。你希望根据天气预报自适应地决定是否带伞。

自适应地做出决策是 BayesOpt 的核心,为了有效地实现这一点,我们需要一个好的预测模型和一个好的决策政策。在这个框架的两个组成部分都需要注意;这就是为什么本章后面的两个主要部分分别涵盖了用 GP 进行建模和用 BayesOpt 政策进行决策。

1.2.4 BayesOpt 的实际应用

此时,你可能想知道所有这些复杂的机器真的是否有效—或者是否比一些简单的策略如随机抽样更有效。为了找出答案,让我们看一下 BayesOpt 在一个简单函数上的“演示”。这也是我们从抽象到具体的好方法,并揭示了我们在后续章节中能做什么。

假设我们试图优化的黑盒客观函数(特别是在这种情况下,最大化)是图 1.6 中的一维函数,从 -5 到 5 定义。同样,这个图片仅供参考;在黑盒优化中,我们实际上不知道客观函数的形状。我们看到客观函数在 -5(大约 -2.4 和 1.5)附近有几个局部极大值,但全局最大值在右侧大约是 4.3。此外,假设我们被允许最多评估客观函数 10 次。

图 1.6 要最大化的目标函数,在随机搜索中浪费资源于不利区域

在看到贝叶斯优化如何解决这个优化问题之前,让我们看看两种基准策略。第一种是随机搜索,在–5 到 5 之间均匀采样;我们得到的任何点都是我们将评估目标的位置。图 1.6 是这样一个方案的结果。这里找到的价值最高的点大约在x = 4 处,其值为f(x) = 3.38。

随机搜索的工作原理

随机搜索涉及在我们的目标函数域内均匀随机选择点。也就是说,我们最终到达域内某点的概率等于我们最终到达其他任何点的概率。如果我们认为搜索空间中有重要区域应该更加关注,我们可以从非均匀分布中抽取这些随机样本,而不是均匀采样。然而,这种非均匀策略需要在开始搜索之前知道哪些区域是重要的。

你可能会觉得不满意的是,这些随机抽样的点中有许多恰好落入 0 周围的区域。当然,许多随机样本聚集在 0 周围只是偶然的,而在另一个搜索实例中,我们可能会发现在另一个区域有许多样本。然而,我们仍然可能浪费宝贵的资源来检查函数的一个小区域,其中包含许多评估。直觉上,扩展我们的评估更有利于我们了解目标函数。

这种扩散评估的想法将我们带到了第二个基准:网格搜索。在这里,我们将搜索空间划分为均匀间隔的段,并在这些段的端点进行评估,就像图 1.7 所示。

图 1.7 网格搜索仍然无法有效缩小好的区域。

这次搜索中的最佳点是最右边的最后一个点,在 5 处进行评估,大约为 4.86。这比随机搜索更好,但仍然错过了实际的全局最优点。

现在,我们准备看贝叶斯优化如何运作!贝叶斯优化从一个随机抽样的点开始,就像随机搜索一样,如图 1.8 所示。

图 1.8 贝叶斯优化的开始与随机搜索相似。

图 1.8 的顶部面板表示对评估点进行训练的高斯过程,而底部面板显示了由期望改进策略计算的分数。请记住,这个分数告诉我们应该如何评价我们搜索空间中的每个位置,我们应该选择下一个评估分数最高的位置。有趣的是,在这一点上,我们的策略告诉我们,我们搜索的-5 到 5 之间的几乎整个范围都很有前景(除了围绕 1 的区域,我们已经进行了一次查询)。这应该是直观的,因为我们只看到了一个数据点,而且我们还不知道其他区域的目标函数是什么样子的。我们的策略告诉我们我们应该探索更多!现在让我们看看从第一个查询到第四个查询时我们模型的状态如何在图 1.9 中。

图 1.9

图 1.9 在四次查询之后,我们确定了第二个最佳点。

四次查询中有三次集中在点 1,这里有一个局部最优点,我们还看到我们的策略建议我们下一步查询这个区域的另一个点。此时,你可能担心我们会陷入这个局部最优区域无法摆脱,无法找到真正的最优点,但我们会看到情况并非如此。让我们快进到图 1.10 中的接下来两次迭代。

图 1.10

图 1.10 在充分探索局部最优点后,我们被鼓励看看其他区域。

在对这个局部最优区域进行五次查询后,我们的策略决定有其他更有前景的区域可供探索——即左边约为-2 和右边约为 4 的区域。这非常令人放心,因为它表明一旦我们足够探索一个区域,贝叶斯优化就不会陷入那个区域。现在让我们看看图 1.11 中进行八次查询后会发生什么。

图 1.11

图 1.11 贝叶斯优化成功地忽略了左边的大区域。

在这里,我们观察到了右边的另外两个点,这些点更新了我们的高斯过程模型和我们的策略。观察均值函数(实线,表示最可能的预测),我们看到它几乎与真实目标函数从 4 到 5 的情况完全匹配。此外,我们的策略(底部曲线)现在非常接近全局最优点,基本上没有其他区域。这很有趣,因为我们并没有彻底检查左边的区域(我们只有一个观察到左边的 0),但我们的模型认为与当前区域相比,无论那个区域的函数长什么样子,都不值得调查。在我们的情况下,这实际上是正确的。

最后,在进行了 10 次查询的搜索结束时,我们的工作流现在在图 1.12 中可视化。现在几乎没有疑问,我们已经确定了约为 4.3 的全局最优点。

图 1.12

图 1.12 贝叶斯优化在搜索结束时找到了全局最优点。

这个例子清楚地向我们展示了贝叶斯优化比随机搜索和网格搜索要好得多。这对我们来说应该是一个非常鼓舞人心的迹象,因为后两种策略是许多机器学习实践者在面临超参数调优问题时使用的方法。

例如,scikit-learn 是 Python 中最流行的机器学习库之一,它提供了model_selection模块用于各种模型选择任务,包括超参数调优。然而,随机搜索和网格搜索是该模块实现的唯一超参数调优方法。换句话说,如果我们确实使用随机或网格搜索调整超参数,我们有很大的提升空间。

总的来说,使用 BayesOpt 可能会导致优化性能显著提高。我们可以快速看几个真实世界的例子:

  • 一份名为“贝叶斯优化优于随机搜索的机器学习超参数调优”的 2020 年研究论文(arxiv.org/pdf/2104.10201.pdf)是 Facebook、Twitter、英特尔等公司的联合研究成果,发现 BayesOpt 在许多超参数调优任务中都非常成功。

  • 弗朗西斯·阿诺德(2018 年诺贝尔奖获得者,加州理工学院教授)在她的研究中使用 BayesOpt 来引导寻找高效催化理想化学反应的酶。

  • 一项名为“通过高通量虚拟筛选和实验方法设计高效的分子有机发光二极管”的研究(www.nature.com/articles/nmat4717)发表在自然上,将 BayesOpt 应用于分子有机发光二极管(一种重要类型的分子)的筛选问题,并观察到效率大幅提高。

还有很多类似的例子。

不适用 BayesOpt 的情况

同样重要的是要知道问题设置不合适的情况以及何时使用 BayesOpt。正如我们所说,当我们的有限资源阻止我们多次评估目标函数时,BayesOpt 是有用的。如果不是这种情况,评估目标函数是廉价的,我们没有理由在观察目标函数时吝啬。

如果我们能够彻底检查密集网格上的目标,那将确保找到全局最优解。否则,可能会使用其他策略,例如 DIRECT 算法或进化算法,这些算法在评估成本较低时通常擅长优化。此外,如果目标梯度的信息可用,梯度算法将更适合。

我希望这一章能激发你的兴趣,让你对即将发生的事情感到兴奋。在下一节中,我们将总结你将在本书中学到的关键技能。

1.3 你将在本书中学到什么?

本书深入理解 GP 模型和 BayesOpt 任务。您将学习如何使用最先进的工具和库在 Python 中实现 BayesOpt 流水线。您还将接触到一系列建模和优化策略,当处理 BayesOpt 任务时。到本书结束时,您将能够做到以下几点:

  • 使用 GPyTorch 实现高性能的 GP 模型,这是 Python 中的首选 GP 建模工具;可视化和评估它们的预测;为这些模型选择适当的参数;并实现扩展,例如变分 GP 和贝叶斯神经网络,以适应大数据

  • 使用最先进的 BayesOpt 库 BoTorch 实现各种 BayesOpt 策略,并与 GPyTorch 很好地集成,并检查以及理解它们的决策策略

  • 使用 BayesOpt 框架处理不同的专业化设置,例如批处理、约束和多目标优化

  • 将 BayesOpt 应用于真实任务,例如调整机器学习模型的超参数

进一步地,我们在练习中使用真实世界的例子和数据来巩固我们在每一章学到的知识。在整本书中,我们在许多不同的环境中运行我们的算法在相同的数据集上,以便我们可以比较和分析不同的方法。

摘要

  • 现实世界中的许多问题可以被看作是昂贵的黑盒优化问题。在这些问题中,我们只观察到函数值,没有任何额外的信息。此外,观察一个函数值是昂贵的,使得许多盲目的优化算法无法使用。

  • BayesOpt 是一种机器学习技术,通过设计目标函数的智能评估来解决这个黑盒优化问题,以便尽快找到最优解。

  • 在 BayesOpt 中,GP 充当预测模型,预测给定位置的目标函数的值。GP 不仅生成均值预测,还通过正态分布表示不确定性的 95% CI。

  • 要优化黑盒函数,BayesOpt 策略会迭代地决定在哪里评估目标函数。该策略通过量化每个数据点在优化方面的帮助程度来实现这一点。

  • 在 BayesOpt 中,GP 和策略是相辅相成的。前者用于进行良好的预测,后者用于做出良好的决策。

  • 通过以自适应的方式做出决策,BayesOpt 在优化方面比随机搜索或网格搜索更好,后者通常用作黑盒优化问题中的默认策略。

  • BayesOpt 在机器学习和其他科学应用中,如药物发现中的超参数调优中取得了显著的成功。

第一部分:使用高斯过程建模

预测模型在贝叶斯优化中起着至关重要的作用,通过准确的预测指导决策。正如我们在第 1.2.1 节中看到的,并且在这一部分中一再看到的,高斯过程提供了校准的不确定性量化,这是任何决策任务中的关键组成部分,也是许多机器学习模型缺乏的特性。

我们从第二章开始,该章解释了高斯过程作为函数分布的直观理解,以及作为无限维度中多元正态分布的一般化。我们通过贝叶斯定理探讨了如何更新高斯过程以反映我们对函数值的信念,考虑到新数据。

第三章展示了高斯过程的数学灵活性。这种灵活性使我们,用户,能够通过全局趋势和高斯过程预测的可变性将先验信息合并到预测中。通过结合高斯过程的不同组件,我们获得了数学建模广泛函数的能力。

我们在本部分的讨论中附带了代码实现,使用了最先进的 Python 高斯过程库 GPyTorch。当你阅读本部分的材料时,你将获得使用 GPyTorch 设计和训练高斯过程模型的实践经验。

第三章:高斯过程作为函数分布

本章内容包括

  • 对多元高斯分布及其属性的速成课程

  • 将 GPs 理解为无限维度中的多元高斯分布

  • 在 Python 中实现 GP

现在我们已经看到贝叶斯优化可以帮助我们做什么,我们已经准备好踏上掌握贝叶斯优化的旅程。正如我们在第一章中看到的,贝叶斯优化工作流程由两个主要部分组成:高斯过程(GP)作为预测模型或替代模型,以及用于决策的策略。使用 GP,我们不仅获得测试数据点的点估计作为预测,而且我们有一个完整的概率分布表示我们对预测的信念。

使用 GP,我们从相似的数据点产生相似的预测。例如,在天气预报中,当估计今天的温度时,GP 会查看与今天相似的几天的气候数据,即最近几天或一年前的这一天。另一个季节的天数不会在进行此预测时通知 GP。同样,当预测房屋价格时,GP 将会说预测目标所在地区的相似房屋比另一个州的房屋更具信息量。

数据点之间的相似程度是使用 GP 的协方差函数来编码的,此外,该函数还模拟了 GP 预测的不确定性。请记住,在第一章中我们对比了岭回归模型和 GP 的模型,再次显示在图 2.1 中。在这里,虽然岭回归器只产生单值预测,但 GP 在每个测试点输出一个正态分布。不确定性量化是将 GP 与其他 ML 模型区分开来的因素,特别是在不确定性决策的背景下。

图 2.1 岭回归和 GP 的预测。尽管 GP 的平均预测与岭回归的预测相同,但 GP 还提供了表示预测不确定性的 CI。

我们将看到如何通过高斯分布在数学上实现相关建模和不确定性量化,并学习如何在 GPyTorch 中实际实现 GP,这是 Python 中首选的 GP 建模工具。能够用 GP 对函数进行建模是迈向贝叶斯优化的第一步,我们将在本章中完成这一步。

为什么选择 GPyTorch?

在 Python 中还有其他的 GP 建模库,如 GPy 或 GPflow,但我们选择了 GPyTorch 作为本书的工具。基于 PyTorch 构建且处于积极维护状态,GPyTorch 提供了从数组操作到 GP 建模再到使用 BoTorch 进行贝叶斯优化的简化工作流程,我们将在第四章开始使用 BoTorch。

该库也在积极维护,并且已实现了许多最先进的方法。例如,第十二章介绍了使用 GPyTorch 对大型数据集进行缩放的方法,在第十三章中,我们学习将神经网络集成到 GP 模型中。

2.1 如何以贝叶斯方式出售您的房屋

在我们立即进入高斯过程之前,让我们考虑一个房价建模领域的示例场景,以及房子价格如何与其他房子相关确定的例子。这个讨论作为多元高斯分布中相关性如何工作的示例,是高斯过程的核心部分。

假设你是密苏里州的一位房主,正打算出售你的房子。你正在尝试确定一个合适的要价,并与朋友讨论如何做到这一点:

你:     我不确定该怎么办。我只是不知道我的房子值多少钱。

朋友: 你有个大概的估算吗?

你:     我猜大概在 15 万到 30 万之间。

朋友: 这个范围挺大的。

你:     是啊,我希望我认识已经卖掉房子的人。我需要一些参考。

朋友: 我听说爱丽丝卖了她的房子 25 万。

你:     在加利福尼亚的阿利克斯吗?这真让人吃惊!而且,我不认为加利福尼亚的房子会帮助我更好地估算自己的房子。它可能仍然在 15 万到 30 万之间。

朋友: 不,是住在你隔壁的爱丽丝。

你:     哦,我明白了。这实际上非常有用,因为她的房子和我的非常相似!现在,我猜我的房子估价在 23 万到 27 万之间。是时候和我的房地产经纪人谈谈了!

朋友: 很高兴我能帮上忙。

在这次对话中,你说使用你的邻居爱丽丝的房子作为参考是估算你自己价格的好策略。这是因为这两个房子在属性上相似,并且彼此物理上靠近,所以你期望它们卖出的价格相似。另一方面,阿利克斯的房子位于加利福尼亚,与我们的房子毫不相关,所以即使你知道她的房子卖了多少钱,你也无法获得任何关于你感兴趣的新信息:你自己的房子值多少钱。

我们刚刚进行的计算是关于我们对房子价格的信念的贝叶斯更新。你可能熟悉贝叶斯定理,如图 2.2 所示。有关贝叶斯定理和贝叶斯学习的优秀介绍,请参阅路易斯·塞拉诺的《精通机器学习》(Manning,2021)第八章。

贝叶斯定理给了我们一种更新我们对我们感兴趣的数量的信念的方法,这种数量在这种情况下是我们房子的合适价格。在应用贝叶斯定理时,我们从先验信念,即我们的第一个猜测,到关于所讨论数量的后验信念。这个后验信念结合了先验信念和我们观察到的任何数据的可能性。

图片

图 2.2 贝叶斯定理,它提供了一种更新对感兴趣的数量的信念的方法,表示为一个随机变量的概率分布。在观察到任何数据之前,我们对 X 有先验信念。在使用数据更新后,我们获得了关于 X 的后验信念。

在我们的例子中,我们首先有一个先验置信度,认为房价在 150k 到 300k 之间。正如你的朋友所说的那样,150k 到 300k 的范围很大,所以在这个初始先验置信度中没有太多信息,任何在 150k 到 300k 之间的价格都是可能的。当我们根据两个房子中任意一个的价格的新信息更新这个范围到后验置信度时,一件有趣的事情发生了。

首先,假设 Alix 在加利福尼亚的房子价值为 250k,我们对我们自己房子的后验置信度保持不变:从 150k 到 300k。同样,这是因为 Alix 的房子与我们的房子无关,她的房子的价格也无法告诉我们我们感兴趣的东西的数量。

其次,如果新的信息是 Alice 的房子,它就在我们的旁边,价值为 250k,那么我们的后验置信度就会从先验置信度大幅改变:变为 230k 到 270k 的范围。有了 Alice 的房子作为参考,我们已经根据观察值 250k 更新了我们的置信度,同时缩小了置信度的范围(从 150k 的差异缩小到 40k 的差异)。这是非常合理的事情,因为 Alice 的房子对我们房子的价格非常具有信息量。图 2.3 可视化了整个过程。

图 2.3 以贝叶斯方式更新我们房价的置信度。根据观察的房价与我们房子的相似程度,后验置信度要么保持不变,要么发生 drastīc 更新。

注意,示例中的数字不是精确的,只是为了使例子更具直观性。然而,我们将看到,使用多元高斯分布来建模我们的置信度,可以以可量化的方式实现这种直观的更新过程。此外,利用这样的高斯分布,我们可以确定一个变量(某人的房子)是否与我们感兴趣的变量(我们自己的房子)足够相似,以影响我们的后验置信度的程度。

2.2 用多元高斯分布和贝叶斯更新建模相关性

在本节中,我们学习多元高斯分布(或多元高斯分布,或简称高斯分布)以及它们如何促进我们之前看到的更新规则。这为我们后续讨论高斯过程奠定了基础。

2.2.1 使用多元高斯分布共同建模多个变量

在这里,我们首先介绍了多元高斯分布是什么以及它们可以模拟的内容。我们将看到,通过使用协方差矩阵,多元高斯分布描述了不仅是单个随机变量的行为,而且还描述了这些变量之间的相关性。

首先,让我们来考虑正态分布,也叫钟形曲线。正态分布在现实世界中非常常见,被用来模拟各种量,比如身高、智商、收入和出生体重。

当我们想要建模多于一个量时,我们将使用 MVN 分布。为此,我们将这些量聚合成一个随机变量向量,然后称此向量遵循 MVN 分布。这个聚合如图 2.4 所示。

图 2.4 MVN 分布将多个正态分布的随机变量组合在一起。虽然 MVN 的均值向量连接了均值,但协方差矩阵模拟了各个变量之间的相关性。

定义 考虑一个随机向量X = [X[1]X[2] ... X[n]],它遵循一个被标记为N(μ, Σ)的高斯分布,其中μ是长度为n的向量,Σ是一个n乘以n的矩阵。在这里,μ被称为均值向量,其各个元素表示X中相应随机变量的期望值,Σ是协方差矩阵,描述了各个变量的方差以及变量之间的相关性。

让我们花一点时间解析 MVN 分布的定义:

  • 首先,由于 MVN 的便利性质,向量X中的每个随机变量都遵循正态分布。具体来说,第i个变量X[i]的平均值为μ[i],这是 MVN 的均值向量μ的第i个元素。

  • 此外,X[i]的方差是协方差矩阵Σ的第i对角条目。

  • 如果我们有一个遵循 MVN 的随机变量向量,那么每个单独的变量对应于一个已知的正态分布。

如果协方差矩阵Σ中的对角线条目是各个变量的方差,那么非对角线条目呢?该矩阵中第i行和第j列的条目表示*X[i]X[j]*之间的协方差,这与两个随机变量之间的相关性有关。假设相关性为正,则以下结论成立:

  • 如果这种相关性很高,那么两个随机变量*X[i]X[j]*被认为是相关的。这意味着如果一个值增加,另一个值也倾向于增加,如果一个值减少,另一个值也会减少。你的邻居爱丽丝的房子和你自己的房子就是相关变量的例子。

  • 另一方面,如果这种相关性很低且接近零,则无论*X[i]的值是什么,我们关于X[j]*值的了解很可能不会发生太大变化。这是因为两个变量之间没有相关性。加利福尼亚州的阿利克斯的房子和我们的房子属于这个类别。

负相关性

前面的描述是针对正相关性的。相关性也可以是负的,表示变量朝相反的方向移动:如果一个变量增加,另一个变量就会减少,反之亦然。正相关性展示了我们在这里想要学习的重要概念,所以我们不会担心负相关性的情况。

为了使我们的讨论更具体,让我们定义一个 MVN 分布,同时对三个随机变量进行建模:我们房子的价格 A;邻居爱丽丝房子的价格 B;以及加利福尼亚的阿利克斯房子的价格 C。这个三维高斯分布的协方差矩阵也在图 2.4 中描述。

注意 通常方便假设这个高斯分布的均值向量归一化为零向量。这种归一化通常在实践中完成,以简化许多数学细节。

再次,对角线单元格告诉我们单个随机变量的方差。B 的方差(3)略大于 A(1),这意味着我们对 B 的值更不确定,因为我们对邻居的房子不了解所有情况,所以不能做出更准确的估计。另一方面,第三个变量 C 具有最大的方差,表示加利福尼亚的房屋价格范围更广。

注意 这里使用的值(1, 3, 10)是示例值,目的是说明随机变量方差越大,我们对该变量的值越不确定(在了解其值之前)。

我们家(A)和邻居家(B)之间的协方差为 0.9,这意味着两栋房子的价格存在着显著的相关性。这是有道理的,因为如果我们知道邻居房子的价格,我们就能更好地估算出我们自己房子的价格,因为它们位于同一条街上。还要注意的是,无论是 A 还是 B 与加利福尼亚房价都没有任何相关性,因为位置上来看,C 与 A 或 B 没有任何共同之处。另一种说法是,即使我们知道加利福尼亚房子的价格,我们也不会对我们自己房子的价格了解多少。现在让我们在图 2.5 中使用平行坐标图来可视化这个三维高斯分布。

图 2.5 平行坐标图可视化了来自房价示例的均值归一化 MVN。误差条表示相应正态分布的 95% CI,而淡化的线显示了从多元高斯中绘制的样本。

注意图中的粗体菱形和相应的误差条:

  • 粗体菱形代表我们高斯分布的均值向量,即零向量。

  • 误差条表示三个单独变量的 95%可信区间(CI)。从 A 到 B 到 C,我们观察到越来越大的 CI,对应着相应方差的增加。

可信区间

正态分布随机变量x的(1 – α) CI 是一个范围,其中x落入这个范围的概率恰好为(1 – α)。统计学家通常对正态分布使用 95% CI。这里并没有什么特别之处,只是因为 95% 是许多统计程序用来确定某事是否有意义的阈值。例如,一个t检验通常使用置信水平 1 – α = 0.95,对应着p值小于 α = 0.05 表示显著结果。关于正态分布的一个方便事实是μ ± 1.96σ是一个 95% CI(有些甚至使用μ ± 2σ),其中μ和σ是变量x的均值和标准差,这是一个容易计算的量。

图 2.5 表示我们关于三栋房子标准化价格的先验信念。从这个先验开始,我们猜测所有三栋房子的标准化价格都为零,并且对我们的猜测有不同程度的不确定性。此外,由于我们正在处理一个随机分布,我们可以从这个 MVN 中抽取样本。这些样本显示为相连的半透明菱形。

2.2.2 更新 MVN 分布

有了一个 MVN 分布,我们将看到如何在本小节中观察到一些数据后更新这个分布。具体地说,跟随本章开头的示例,我们想要根据观察到的 B 或 C 的值推导出关于这些价格的后验信念。这是一个重要的任务,因为这是 MVN 以及 GP 从数据中学习的方式。

定义 这个更新过程有时被称为条件设定:推导出一个变量的条件分布,在我们已知某个其他变量的值的情况下。更具体地说,我们正在将我们的信念——一个联合三元高斯——条件设定为 B 或 C 的值,获得这三个变量的联合后验分布。

在这里,利用图 2.2 中的贝叶斯定理,我们可以得出这个后验分布的闭式形式。然而,推导过程相当数学密集,所以我们不会在这里详细介绍。我们只需要知道,我们有一个公式,可以插入我们想要条件的 B 或 C 的值,然后这个公式会告诉我们 A、B 和 C 的后验分布是什么。令人惊讶的是,高斯的后验分布是根据同样是高斯的数据进行条件设定的,我们可以获得指定后验高斯的确切后验均值和方差。(在本章后面,我们会看到当我们在 Python 中实现 GP 时,GPyTorch 会为我们处理这个数学密集的更新。)

注意 对于感兴趣的读者,这个公式及其推导可以在 Carl Edward Rasmussen 和 Christopher K. I. Williams(MIT Press,2006)的书籍Gaussian Processes for Machine Learning的第二章第二部分中找到,这本书通常被认为是 GP 的圣经。

现在让我们重新生成平行坐标图,以 B = 4 作为 B 的示例值进行条件限制。结果如图 2.6 所示。

图 2.6 平行坐标图可视化了图 2.5 中 MVN 在 B = 4 条件下的情况。在这里,A 的分布被更新,所有绘制的样本都插值为 B = 4。

在用关于 B 的观察更新我们的信念后,我们的后验信念发生了一些变化:

  • A 的分布会发生变化,由于 A 和 B 之间的正相关关系,其均值会略微增加。此外,其误差范围现在变得更小。

  • B 的后验分布简单地变成了一个具有零方差的特殊正态分布,因为我们现在对后验中的 B 值完全确定。换句话说,对于 B 的值不再存在不确定性。

  • 与此同时,C 的分布在更新后保持不变,因为它与 B 没有相关性。

所有这些都是有道理的,并且与我们从房价示例中得到的直觉相符。具体来说,当我们得知邻居的房子价格时,关于我们自己房子的信念被更新为类似于观察到的价格,并且我们的不确定性也减少了。

当我们以 C 为条件时会发生什么?正如您可能猜到的那样,由于在 B 的值上进行条件限制后 C 保持不变的原因,当我们以 C 为条件时,A 和 B 的后验分布都保持不变。图 2.7 展示了 C = 4 的情况。

图 2.7 平行坐标图可视化了图 2.5 中 MVN 在 C = 4 条件下的情况。在这里,没有其他边缘分布改变。所有绘制的样本都插值为 C = 4。

当我们得知加利福尼亚州的一栋房子被卖掉时,我们发现这栋房子与我们在密苏里州的房子无关,因此对我们房子价格的信念保持不变。

图 2.6 和图 2.7 还有另一个有趣之处。请注意,在图 2.6 中,当我们以 B = 2 为条件时,我们绘制的所有后验 MVN 的样本都经过点 (B, 2)。这是因为在我们的后验信念中,我们对于 B 取值不再有任何不确定性,因此从后验分布中绘制的任何样本都需要满足此条件的约束。图 2.7 中的点 (C, 2) 也是同样道理。

从视觉上来看,您可以将这理解为当我们对一个变量进行条件限制时,我们将从先验分布(在图 2.5 中)绘制的样本在相应的变量条件处“绑定”成一个结,如图 2.8 所示。

图 2.8 在观察上对高斯进行条件限制类似于在该观察周围打结。后验分布中的所有样本都需要通过这个结,且在观察点没有不确定性。

最后,我们可以通过类似于图 2.3 的图表来描述我们刚刚进行过的贝叶斯条件过程。这在图 2.9 中显示。

图 2.9:以贝叶斯方式更新关于我们房子价格的信念。根据观察到的房价与我们房子的相似程度,后验信念要么保持不变,要么得到极大更新。

再次,如果我们将高斯分布条件设置为 C,则不相关变量的后验分布保持不变。然而,如果我们将其设置为 B,与之相关的变量 A 会得到更新。

2.2.3:用高维高斯分布建模多个变量

MVN 分布不仅需要包含三个随机变量;事实上,它可以同时模拟任意数量的变量。在本小节中,我们了解到高维高斯分布的工作方式与我们迄今所见的相同。所以,我们可以说,一个代表三座房子的 3 维高斯分布,我们可以将其替换为一个编码街道上一系列房屋信息的 20 维高斯分布。甚至更高维的高斯分布可以模拟城市或国家中的房屋。

此外,通过这些平行坐标图,我们可以同时可视化高维高斯分布的所有单个变量。这是因为每个变量对应一个单独的误差条,只占用x轴上的一个位置。

我们再次将其均值向量归一化为零向量,虽然显示其 20×20 的协方差矩阵不方便,但我们可以绘制一个热图来可视化这个矩阵,就像图 2.10 中所示的那样。

图 2.10:显示 20 维高斯分布的协方差矩阵的热图。相邻变量之间的相关性比远离的变量更大,这由较深的色调表示。

对角线条目,或者单个变量的方差,在这种情况下都是 1。此外,变量被排序,使得彼此靠近的变量是相关的;也就是说,它们的协方差取值较大。相反,彼此远离的变量则不太相关,它们的协方差接近于零。例如,在这个高斯分布中的任意一对连续变量(第一个和第二个,第二个和第三个等等)的协方差大约为 0.87。也就是说,任意两个相邻的房子的协方差为 0.87。如果我们考虑第 1 个和第 20 个变量——也就是,街道一端的房子和另一端的房子——它们的协方差实际上是零。

这非常直观,因为我们期望附近的房子价格相似,所以一旦我们知道一座房子的价格,我们就可以更多地了解该地区周围房屋的价格,而不是远处房屋的价格。

图 2.11 误差线和从先验(左)和后验(右)高斯分布中抽取的样本,以第 10 个变量的值为 2 条件。接近第 10 个变量的变量在后验中的不确定性减少,其均值更新为接近 2。

这在平行坐标图中如何体现?图 2.11 显示了我们左侧的先验高斯和右侧第 10 个变量的值为 2 的后验高斯。基本上,我们正在模拟以下事件,即我们发现第 10 座房子的价格为 2(其确切单位被省略):

  • 首先,我们再次看到这种现象,在后验分布中,误差线和样本围绕我们条件观察到的观察点打成结。

  • 其次,由于协方差矩阵施加的相关结构,接近第 10 个变量的变量其均值被“拉高”,以便均值向量现在平滑地插值到点 (10, 2)。这意味着我们更新了我们的信念,因为周围的房屋现在其价格增加以反映我们学到的信息。

  • 最后,在这一点 (10, 2) 条件之后,围绕这一点的不确定性(由误差线表示)减小。这是一个非常好的属性,因为直觉上,如果我们知道一个变量的值,我们应该对与我们所知变量相关的其他变量的值更加确定。也就是说,如果我们知道一座房子的价格,我们就会更加确定附近房屋的价格。这个属性是高斯过程提供的校准不确定性量化的基础,我们在本章的下一节中会看到。

2.3 从有限到无限的高斯分布

现在我们准备讨论什么是高斯过程。与前面三个变量 A、B 和 C,或者之前章节中的 20 个变量相同的方式,我们假设现在有无限多个变量,所有这些变量都属于多元正态分布。这个 无限维 高斯被称为 高斯过程

想象一下在一个非常大而密集的区域内预测房价。整个区域的规模如此之大,以至于如果我们离一座房子移动了一个非常小的距离,我们就会到达另一座房子。鉴于高斯分布中变量(房屋)的高密度,我们可以将整个区域视为有无限多个房屋;也就是说,高斯分布有无限多个变量。

图 2.12 使用不同数量的变量对加州房价建模。我们拥有的变量越多,我们的模型就越平滑,我们越接近无限维模型。

这在图 2.12 中使用一个包含加利福尼亚州 5,000 个房价的数据集进行了说明。在左上角的面板中,我们展示了散点图中的个别数据点。在其余的面板中,我们使用各种数量的变量来对数据进行建模,其中每个变量对应于加利福尼亚地图内的一个区域。随着变量数量的增加,我们的模型变得更加精细。当这个数量是无限的时候——也就是说,当我们可以在这张地图上的任何区域进行预测,无论多么小——我们的模型存在于一个无限维空间中。

这正是高斯过程的含义:在无限维空间中的高斯分布。在任何区域进行预测的能力帮助我们摆脱了有限维多元正态分布,并获得了一个机器学习模型。严格来说,当变量有无穷多个时,高斯分布的概念并不适用,因此技术上的定义如下。

定义 高斯过程是一组随机变量,使得这些变量的任意有限子集的联合分布是一个多元正态分布。

这个定义意味着,如果我们有一个高斯过程模型来描述一个函数 ƒ,那么在任何一组点处的函数值都由一个多元正态分布来建模。例如,变量向量 [ƒ(1) ƒ(2) ƒ(3)] 遵循一个三维高斯分布;[ƒ(1) ƒ(0) ƒ(10) ƒ(5) ƒ(3)] 遵循另一个不同的、五维的高斯分布;以及 [ƒ(0.1) ƒ(0.2) ... ƒ(9.9) ƒ(10)] 遵循另一个高斯分布。

图 2.13 不同高斯分布的平行坐标图。高斯过程的任何有限子集都是多元正态分布。随着变量数量趋于无穷,我们得到一个高斯过程,并且可以在域中的任何地方进行预测。

这在图 2.13 中有所说明。前三个面板显示了平行坐标图中的三元高斯分布,分别为 [ƒ(–2) ƒ(1) ƒ(4)],一个 11 元高斯分布,对应 [ƒ(–4.5) ƒ(–4) ... ƒ(4) ƒ(4.5)],以及一个在更密集的网格上的 101 元高斯分布。最后,在最后一个面板中,我们有无限多个变量,这给了我们一个高斯过程。

由于我们现在处于无限维空间,讨论均值向量和协方差矩阵已经没有意义了。相反,我们在高斯过程中有的是一个均值 函数 和一个协方差 函数,但这两个对象的角色仍然与多元正态分布相同:

  • 首先是均值函数,它接受一个输入 x,计算函数值 ƒ(x) 的期望。

  • 第二,协方差函数接受两个输入,x[1] 和 x[2],并计算两个变量 ƒ(x[1]) 和 ƒ(x[2]) 之间的协方差。如果 x[1] 等于 x[2],那么这个协方差值就是 ƒ(x) 的正态分布的方差。如果 x[1] 不同于 x[2],协方差表示两个变量之间的相关性。

由于均值和协方差是函数,我们不再受限于固定数量的变量——相反,我们有效地拥有无限数量的变量,并且可以在任何地方进行我们的预测,如图 2.13 所示。这就是为什么尽管高斯过程具有多元正态分布的所有特性,但高斯过程存在于无限维度的原因。

出于同样的原因,高斯过程可以被视为对函数的分布,就像本章的标题所暗示的那样。本章我们经历的从一维正态分布到高斯过程的过程,在表 2.1 中总结。

表 2.1 高斯分布对象及其模拟对象。使用高斯过程时,我们在无限维度下操作,模拟函数而不是数字或向量。

分布类型模拟变量数量描述
一维正态分布一个数字的分布
多元正态分布有限数量有限长度向量的分布
高斯过程无限数量函数的分布

要看高斯过程的实际应用,让我们重新审视本章开头图 2.1 中的曲线拟合过程,我们将我们的域限制在-5 到 5 之间。如图 2.14 所示。

图 2.14 高斯过程对零、一个、两个和四个观测条件下的预测

在每个面板中,以下内容为真:

  • 中间的实线是平均函数,类似于图 2.11 中连接菱形的实线。

  • 另一方面,阴影区域是域上的 95% CI,对应于图 2.11 中的误差条。

  • 各种曲线是从相应高斯过程中抽取的样本。

在观察任何数据之前,我们从左上角的先验高斯过程开始。就像先验多元正态分布一样,在没有训练数据的情况下,我们的先验高斯过程产生恒定的均值预测和不确定性。这是一个合理的行为。

当我们将高斯过程与各种数据点进行条件化时,有趣的部分就出现了。这在图 2.14 剩余的面板中进行了可视化。正如多元正态分布的离散情况一样,高斯过程在连续域中工作时,均值预测以及从后验分布中抽取的样本在训练集的数据点之间平滑插值,而关于函数值的不确定性,由置信区间(CI)量化,在这些观测点周围的区域平滑减少。这就是我们所说的校准不确定性量化,这是高斯过程的最大卖点之一。

高斯过程的平滑性

平滑性 特性是指要求相似点彼此相关的约束。换句话说,相似的点应该产生相似的函数值。这也是为什么当我们在顶部右侧面板的数据点处于 3 时,2.9 和 3.1 处的平均预测会更新为比它们的先验均值大的原因。这些点,2.9 和 3.1,与 3 相似,因为它们彼此接近。这种平滑性是通过 GP 的协方差函数来设置的,这是第三章的主题。虽然我们迄今为止看到的示例都是一维的,但当我们的搜索空间是更高维度时,这种平滑性仍然保持,正如我们后面看到的那样。

总的来说,我们已经看到,GP 在扩展到无限维度时是多元正态分布,由于高斯分布具有许多便利的数学特性,GP 不仅产生平均预测,而且通过其预测协方差以一种原则性的方式量化了我们对函数值的不确定性。平均预测正好通过训练数据点,并且不确定性在这些数据点处收敛。

建模非高斯数据

在现实生活中,并不是所有的数据都遵循高斯分布。例如,对于限制在数值范围内的值或不遵循钟形分布的变量,高斯分布是不合适的,可能会导致低质量的预测。

在这些情况下,我们可以应用各种数据处理技术来“转换”我们的数据点以遵循高斯分布。例如,Box-Muller 变换是一个从均匀分布的随机数生成一对正态分布的随机数的算法。有兴趣的读者可以在 Wolfram 的 MathWorld 上找到有关此算法的更多详细信息 (mathworld.wolfram.com/Box-MullerTransformation.xhtml)。

2.4 在 Python 中实现 GPs

在本章的最后一节中,我们迈出了实现 Python 中高斯过程(GPs)的第一步。我们的目标是熟悉我们将用于此任务的库的语法和 API,并学习如何重新创建我们迄今为止看到的可视化效果。这个动手实践部分还将帮助我们更深入地理解 GPs。

首先,确保您已经下载了本书的附带代码并安装了必要的库。关于如何做到这一点的详细说明包含在前言中。我们使用包含在 Jupyter 笔记本 CH02/01 - Gaussian processes.ipynb 中的代码。

2.4.1 设置训练数据

在我们开始实现 GP 模型的代码之前,让我们先花一些时间创建我们想要建模的目标函数和训练数据集。为此,我们需要导入 PyTorch 用于计算和操作张量,以及 Matplotlib 用于数据可视化:

import torch
import matplotlib.pyplot as plt

在本示例中,我们的目标函数是一维的 Forrester 函数。Forrester 函数是具有一个全局最大值和一个局部最大值的多模式函数(www.sfu.ca/~ssurjano/forretal08.xhtml),使得拟合和找到函数的最大值成为一个非平凡的任务。该函数具有以下公式:

这是如下实现的:

def forrester_1d(x):
    y = -((x + 1) ** 2) * torch.sin(2 * x + 2) / 5 + 1
    return y.squeeze(-1)

让我们快速在图表中绘制此函数。在这里,我们限制自己在-3 到 3 之间的域,并在此范围内的 100 个点的密集网格上计算此 Forrester 函数。我们还需要一些用于训练的样本点,我们通过torch.rand()进行随机采样并将其存储在train_x中;train_y包含这些训练点的标签,可以通过评估forrester_1d(train_x)获得。此图由以下代码生成,产生图 2.15:

xs = torch.linspace(-3, 3, 101).unsqueeze(1)
ys = forrester_1d(xs)

torch.manual_seed(0)
train_x = torch.rand(size=(3, 1)) * 6 - 3
train_y = forrester_1d(train_x)

plt.figure(figsize=(8, 6))

plt.plot(xs, ys, label="objective", c="r")
plt.scatter(train_x, train_y, marker="x", c="k", label="observations")

plt.legend(fontsize=15);

图 2.15 显示的是当前示例中使用的目标函数,如实线所示。标记表示训练数据集中的点。

我们看到的三个标记是我们随机选择包含在我们的训练数据集中的点。这些训练数据点的位置存储在train_x中,它们的标签(这些位置处的 Forrester 函数的值)存储在train_y中。这设置了我们的回归任务:在这三个数据点上实现和训练一个 GP,并在范围在-3 到 3 之间的点上可视化其预测。在这里,我们还创建了xs,这是在此范围内的密集网格。

2.4.2 实现 GP 类

在本小节中,我们将学习如何在 Python 中实现 GP 模型。我们使用 GPyTorch 库,这是一种现代 GP 建模的最先进工具。

重要提示:GPyTorch 的设计理念是遵循 DL 库 PyTorch,并使其所有模型类扩展基础模型类。如果您熟悉在 PyTorch 中实现神经网络,您可能会知道这个基类是torch.nn.Module。使用 GPyTorch,我们通常扩展gpytorch.models.ExactGP类。

要实现我们的模型类,我们使用以下结构:

import gpytorch

class BaseGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        ...

    def forward(self, x):
        ...

在这里,我们实现了一个名为BaseGPModel的类,它有两个特定方法:__init__()forward()。我们的 GP 模型的行为在很大程度上取决于我们如何编写这两个方法,无论我们想要实现什么样的 GP 模型,我们的模型类都需要具有这些方法。

让我们首先讨论__init__()方法。它的作用是接收由第一个和第二个参数train_xtrain_y定义的训练数据集,以及存储在变量likelihood中的似然函数,并初始化 GP 模型,即一个BaseGPModel对象。我们实现该方法如下:

def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ZeroMean()
        self.covar_module = gpytorch.kernels.RBFKernel()

在这里,我们简单地将三个输入参数传递给我们的超类的__init__()方法,而gpytorch.models.ExactGP的内置实现则为我们处理了繁重的工作。剩下的是均值和协方差函数的定义,正如我们所说的,这是高斯过程的两个主要组成部分。

在 GPyTorch 中,对于均值和协方差函数都有广泛的选择,我们将在第三章中进行探讨。现在,我们使用 GP 的最常见选项:

  • 对于均值函数,我们使用gpytorch.means.ZeroMean(),在先验模式下输出零均值预测

  • 对于协方差函数,我们使用gpytorch.kernels.RBFKernel(),它实现了径向基函数(RBF)核——这是高斯过程中最常用的协方差函数之一,它实现了数据点彼此接近时彼此相关的思想。

我们分别将这些对象存储在mean_modulecovar_module类属性中。这就是我们对__init__()方法需要做的全部。现在,让我们将注意力转向forward()方法。

forward()方法非常重要,因为它定义了模型如何处理其输入。如果您在 PyTorch 中使用过神经网络,那么您知道网络类的forward()方法会依次将其输入通过网络的层,并且最终层的输出就是神经网络产生的内容。在 PyTorch 中,每个层都被实现为一个模块,这是一个表示 PyTorch 中任何处理数据的对象的基本构建块的术语。

在 GPyTorch 中,高斯过程的forward()方法工作方式类似:高斯过程的均值和协方差函数被实现为模块,并且该方法的输入被传递给这些模块。我们不是将结果依次通过不同的模块传递,而是同时将输入传递给均值和协方差函数。这些模块的输出然后被组合以创建一个多变量正态分布。PyTorch 和 GPyTorch 之间的这种差异在图 2.16 中说明。

图 2.16 PyTorch 和 GPyTorch 如何在其各自的forward()方法中处理数据的过程。输入由不同的模块处理以产生最终输出,对于前馈神经网络来说,输出是一个数字,对于 GP 来说是一个多变量正态分布。

forward()方法在以下代码中实现:

def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

这里的逻辑非常简单:因为我们有一个均值函数和一个协方差函数,所以我们只需在输入x上调用它们来计算均值和协方差预测。最后,我们需要返回的是一个 MVN 分布,由gpytorch.distributions.MultivariateNormal类实现,带有相应的均值和协方差。换句话说,我们只是根据模型类的mean_ modulecovar_module属性计算出的均值向量和协方差矩阵创建了一个 MVN 分布。

就是这样!用 GPyTorch 实现 GP 模型是多么容易,令人惊讶。对我们来说最大的收获是,我们需要在 __init__() 方法中实现均值和协方差函数。在 forward() 方法中,当我们需要进行预测时,我们只需在输入上调用这两个函数即可。

2.4.3 用 GP 进行预测

有了 BaseGPModel 类,我们可以开始用 GP 进行预测了!回想一下,在 __init__() 方法中,我们需要传递一个似然函数 likelihood,以及我们的训练数据。在许多回归任务中,gpytorch.likelihoods.GaussianLikelihood 对象就足够了。我们可以这样创建这个对象:

likelihood = gpytorch.likelihoods.GaussianLikelihood()

现在,我们可以初始化我们的 BaseGPModel 对象了。但在我们将其初始化为我们的三项训练数据之前,我们首先想要用先验 GP 进行预测。

要初始化一个没有任何训练数据的 GP 对象,我们将 None 传递给训练特征 (train_x) 和标签 (train_y)。因此,我们的先验 GP 被创建如下:

model = BaseGPModel(None, None, likelihood)

最后,在我们进行任何预测之前,需要进行一些簿记工作。首先,我们设置 GP 的超参数:

lengthscale = 1
noise = 1e-4

model.covar_module.lengthscale = lengthscale
model.likelihood.noise = noise

model.eval()
likelihood.eval()

我们将在第三章讨论每个超参数控制的内容。目前,我们只使用我个人喜欢的默认值:长度尺度为 1,噪声方差为 0.0001。最后一个细节是通过调用相应对象的 eval() 方法在 GP 模型和其似然性中启用预测模式。

处理完这些簿记任务后,我们现在终于可以在我们的测试数据上调用这个 GP 模型进行预测了。我们可以这样做:

with torch.no_grad():
    predictive_distribution = likelihood(model(xs))

请记住,在模型类的 forward() 方法中,我们返回 MVN 分布,因此当我们通过 model(xs) 使用我们的模型传递一些测试数据时,这就是输出。(在 PyTorch 的语法中,调用 model(xs) 是对测试数据 xs 调用 forward() 方法的一种简写。)我们还将相同的输出通过似然函数 likelihood,将噪声方差合并到我们的预测中。简而言之,我们在 predictive_distribution 中存储的是代表测试点 xs 的预测的 MVN 分布。此外,我们在 torch.no_grad() 上下文中计算此对象,当我们不希望 PyTorch 跟踪这些计算的梯度时,这是一个好的做法。

注意 我们只想在优化模型参数时计算操作的梯度。但是当我们想要进行预测时,我们应该完全固定我们的模型,因此禁用梯度检查是适当的。

2.4.4 可视化 GP 的预测

有了这个预测的高斯分布,我们现在可以重新创建我们迄今为止见过的 GP 图。这些图中的每一个都包括一个均值函数 μ,我们可以从 MVN 中获取

predictive_mean = predictive_distribution.mean

另外,我们想显示 95% 置信区间。从数学上讲,这可以通过提取预测协方差矩阵Σ的对角元素来完成(请记住,这些元素表示单个方差σ²),取这些值的平方根以计算标准差σ,并计算 μ ± 1.96σ 的置信区间范围。

幸运的是,当我们使用 GP 时,计算 95% 置信区间是一个常见的操作,所以 GPyTorch 提供了一个方便的辅助方法,称为 confidence_ region(),我们可以直接从 MVN 分布对象中调用该方法:

predictive_lower, predictive_upper =
    ➥    predictive_distribution.confidence_region()

此方法返回两个 Torch 张量的元组,分别存储置信区间的下限和上限端点。

最后,我们可能希望从当前 GP 模型中为我们的图形绘制样本。我们可以通过直接从高斯对象 predictive_distribution 调用方法 sample() 来完成这个操作。如果我们不传入任何输入参数,该方法将返回一个样本。在这里,我们希望从我们的 GP 中抽取五次样本,操作如下:

torch.manual_seed(0)
    samples = predictive_distribution.sample(torch.Size([5]))

我们传递一个 torch.Size() 对象以表示我们希望返回五个样本。在抽样之前设置随机种子是一种确保代码可重现性的良好实践。有了这些,我们就可以开始绘制一些图形了!

第一件事是简单地绘制均值函数:

plt.plot(xs, predictive_mean.detach(), label="mean")

至于 95% 置信区间,我们通常使用像我们到目前为止看到的那样的阴影区域,这可以使用 Matplotlib 的 fill_between() 函数完成:

plt.fill_between(
    xs.flatten(),
    predictive_upper,
    predictive_lower,
    alpha=0.3,
    label="95% CI"
)

最后,我们绘制单独的样本:

for i in range(samples.shape[0]):
    plt.plot(xs, samples[i, :], alpha=0.5)

此代码将生成图 2.17 中的图形。

图 2.17 先验 GP 的预测,具有零均值和 RBF 核。虽然均值和置信区间是恒定的,但单个样本表现出复杂的非线性行为。

我们看到,在整个域上,我们的先验 GP 产生的均值函数在零处保持恒定,我们的 95% 置信区间也是恒定的。这是可以预期的,因为我们使用了一个 gpytorch .means.ZeroMean() 对象来实现均值函数,并且在没有任何训练数据的情况下,我们的先验预测默认为此 0 值。

话虽如此,均值和置信区间只是期望的测量值:它们表示我们预测的平均行为在许多不同的实现中的平均行为。然而,当我们绘制单个样本时,我们会看到每个样本都有一个非常复杂的形状,根本不是恒定的。所有这些都是说,虽然我们在任何点的预测期望值为零,但它可能取值范围很广。这表明 GP 可以以灵活的方式对复杂的非线性行为进行建模。

到目前为止,我们已经学会了在没有任何训练数据的情况下生成和可视化预测的先前 GP 的过程。现在,让我们实际上在我们随机生成的训练集上训练一个 GP 模型,并观察预测结果的变化。迄今为止,我们编码的好处就是这一切可以完全重复进行,只不过现在我们要用我们的训练数据来初始化 GP(请记住,之前我们在第一个和第二个参数中使用了 None):

model = BaseGPModel(train_x, train_y, likelihood)

重新运行代码将给我们 Figure 2.18。

图 2.18 由后验 GP 生成的预测。均值函数和随机绘制样本能很好地插值训练数据点,同时不确定性在这些数据点周围消失。

这正是我们想要看到的预测类型:均值线和样本完美地插值了我们观测到的数据点,并且我们的不确定性(由 CI 表示)在这些数据点周围也减小了。

我们已经可以看到这种不确定性量化在建模目标函数方面的作用。在仅观测到三个数据点之后,我们的 GP 对真实目标函数有一个相当好的近似。实际上,几乎所有的目标函数都在 95% CI 内,表明我们的 GP 成功地考虑了目标函数的行为方式,即使在我们尚未从函数中获得任何数据的区域也是如此。这种校准的量化在我们实际需要根据 GP 模型做出决策的情况下特别有益——也就是说,在我们决定在哪些点上观察函数值以找到最优值时——但是让我们将其保留到本书的下一部分。

2.4.5 超越一维目标函数

到目前为止,我们只看到了在一维目标函数上训练的 GP 的示例。然而,GP 并不局限于一个维度。实际上,只要我们的均值和协方差函数能够处理高维输入,GP 就可以在高维上高效地运算。在本小节中,我们将学习如何训练一个基于二维数据集的 GP。

我们按照前一节的步骤进行。首先,我们需要一个训练数据集。在这里,我将自行创建一个带有点(0, 0),(1, 2),和(–1, 1)以及相应标签为 0,–1 和 0.5 的虚拟集合。换句话说,我们从中学习的目标函数在(0, 0)处的值为 0,在(1, 2)处的值为–1,在(–1, 1)处的值为 0.5。我们希望在[–3, 3]-by-[–3, 3]的正方形内进行预测。

在 Python 中的设置如下:

# training data
train_x = torch.tensor(
    [
        [0., 0.],
        [1., 2.],
        [-1., 1.]
    ]
)

train_y = torch.tensor([0., -1., 0.5])

# test data
grid_x = torch.linspace(-3, 3, 101)               ❶

grid_x1, grid_x2 = torch.meshgrid(grid_x, grid_x,
➥indexing="ij")                                  ❷
xs = torch.vstack([grid_x1.flatten(), grid_x2.flatten()]).transpose(-1, -2)

❶ 一维点阵

❷ 二维点阵

变量 xs 是一个 10,201×2 的矩阵,其中包含了我们希望进行预测的正方形上的点阵。

重要提示:这里有 10,201 个点,因为我们在两个维度中各取了一个 101 点的网格。现在,我们只需重新运行之前用于训练 GP 并对这个二维数据集进行预测的 GP 代码。请注意,不需要修改我们的BaseGPModel类或任何预测代码,这非常了不起!

不过,我们确实需要更改的一件事是如何可视化我们的预测。由于我们正在操作二维空间,因此在单个图中绘制预测均值和 CI 变得更加困难。在这里,一个典型的解决方案是绘制预测均值的热力图和预测标准差的另一个热力图。虽然标准差不完全等于 95% CI,但这两个对象实质上量化了同样的内容:我们对函数值的不确定性。

因此,我们现在不再调用predictive_distribution.confidence_region(),而是这样提取预测标准差:

predictive_stddev = predictive_distribution.stddev

现在,要绘制热力图,我们使用 Matplotlib 中的imshow()函数。我们需要注意predictive_meanpredictive_stddev的形状。这里是一个长度为 10,000 的张量,因此在传递给imshow()函数之前,需要将其重塑为一个方阵。可以按如下方式完成:

fig, ax = plt.subplots(1, 2)

ax[0].imshow(
    predictive_mean.detach().reshape(101, 101).transpose(-1, -2),
    origin="lower",
    extent=[-3, 3, -3, 3]
)                           ❶

ax[1].imshow(
    predictive_stddev.detach().reshape(101, 101).transpose(-1, -2),
    origin="lower",
    extent=[-3, 3, -3, 3]
)                           ❷

❶ 预测均值的第一个热力图

❷ 预测标准差的第二个热力图

此代码生成了图 2.19 中的两个热力图。

图 2.19 由二维 GP 进行的预测。均值函数仍然与训练数据一致,并且在这些数据点周围的区域,不确定性再次消失。

我们看到,在一维情况下我们所拥有的东西在这个例子中也有所延续:

  • 左侧面板显示我们的均值预测与训练数据一致:左侧的亮斑点对应于 (-1, 1),其值为 0.5,而右侧的暗斑点对应于 (1, 2),其值为 -1(我们在 (0, 0) 处的观测值为 0,这也是先验均值,因此在左侧面板中并不像其他两个那样明显)。

  • 我们的不确定性(由预测标准差测量)在我们的训练数据的三个点周围接近零,如右面板所示。远离这些数据点,标准差平滑地增加到归一化的最大不确定性 1。

这意味着在高维情况下,GP 的所有良好属性,如平滑插值和不确定性量化,都得以保留。

这标志着第二章的结束。我们已经对 GP 的概念有了理解,并学会了如何使用 GPyTorch 在 Python 中实现基本的 GP 模型。如前所述,我们将在第三章深入探讨 GP 的均值和协方差函数,包括它们的超参数,以及每个组件如何控制我们的 GP 模型的行为。

2.5 练习

在这个练习中,我们在第一章中看到的一个真实数据集上训练了一个高斯过程(GP),该数据集在表 2.2 中再次显示。每个数据点(行)对应于通过在不同比例下混合铅(Pb)、锡(Sn)、锗(Ge)和锰(Mn)—这些称为父化合物—创建的合金(一种金属)。特征包含在前四列中,这些是父化合物的百分比。预测目标,混合温度,在最后一列中,表示合金形成的最低温度。任务是根据合金的组成百分比预测混合温度。

表 2.2 来自材料发现任务的数据。特征是材料结构以父化合物百分比表示,预测目标是混合温度。

Pb 的%Sn 的%Ge 的%Mn 的%混合温度(°F)
0.500.500.000.00192.08
0.330.330.330.00258.30
0.000.500.500.00187.24
0.000.330.330.33188.54

有多个步骤:

  1. 创建表 2.2 中包含的四维数据集。

  2. 将第五列标准化,方法是减去所有值的均值并将结果除以它们的标准差。

  3. 将前四列视为特征,第五列视为标签。在这些数据上训练一个高斯过程(GP)。您可以重复使用我们在该章中实现的 GP 模型类。

  4. 创建一个测试数据集,其中含有零百分比的锗和锰的组合。换句话说,测试集是一个在单位正方形上的网格,其轴是铅和锡的百分比。

    测试集应该如下 PyTorch 张量所示:

    tensor([[0.0000, 0.0000, 0.0000, 0.0000],
            [0.0000, 0.0100, 0.0000, 0.0000],
            [0.0000, 0.0200, 0.0000, 0.0000],
            ...,
            [1.0000, 0.9800, 0.0000, 0.0000],
            [1.0000, 0.9900, 0.0000, 0.0000],
            [1.0000, 1.0000, 0.0000, 0.0000]])
    

    注意第三列和第四列都是零。

  5. 对此测试集预测混合温度。也就是说,计算测试集中每个点的标准化混合温度的后验均值和标准差。

  6. 可视化预测。这涉及以与图 2.19 相同的方式将均值和标准差显示为热图。解决方案包含在 CH02/02 - Exercise.ipynb 中。

总结

  • 多元高斯分布(MVN)模型化了许多随机变量的联合分布。均值向量表示变量的期望值,而协方差矩阵模拟了这些变量的方差和它们之间的协方差。

  • 通过应用贝叶斯定理,我们可以计算 MVN 的后验分布。通过这种贝叶斯更新,与观察变量相似的变量被更新以反映这种相似性。总的来说,相似的变量产生相似的预测。

  • 高斯过程(GP)将多元正态分布(MVN)扩展到无限维,使其成为函数的分布。然而,GP 的行为仍然类似于 MVN 分布。

  • 即使没有任何训练数据,GP 仍然可以根据先验 GP 产生预测。

  • 训练完成后,高斯过程(GP)的平均预测平滑地插值了训练数据点。

  • 使用高斯过程的最大优势之一是模型提供的不确定性的校准量化:对于观察到的数据点周围的预测更加自信;另一方面,远离训练数据的预测则更不确定。

  • 使用多变量正态分布或高斯过程进行条件化在视觉上类似于在观察点处打结。这迫使模型确切地经过观察,并将不确定性减少到零。

  • 使用 GPyTorch 实现 GP 时,我们可以以模块化的方式编写一个扩展基类的模型类。具体来说,我们实现了两个特定的方法:__init__(),声明了 GP 的均值和协方差函数;forward(),为给定输入构造了一个多变量正态(MVN)分布。

第四章:使用均值和协方差函数自定义高斯过程

本章涵盖内容

  • 使用均值函数控制 GP 的预期行为

  • 使用协方差函数控制 GP 的平滑度

  • 使用梯度下降学习 GP 的最优超参数

在第二章中,我们了解到均值和协方差函数是高斯过程 (GP) 的两个核心组成部分。即使在实现我们的 GP 时使用了零均值和 RBF 协方差函数,你在这两个组成部分上可以选择很多选项。

通过选择均值或协方差函数的特定选择,我们实际上在为我们的 GP 模型指定先验知识。将先验知识纳入预测是我们在任何贝叶斯模型中都需要做的事情,包括 GP。虽然我说我们需要这样做,但将先验知识纳入模型中总是一件好事,特别是在数据获取昂贵的情况下,比如贝叶斯优化。

例如,在天气预报中,如果我们想要估计一月份在密苏里的典型天气的温度,我们不必进行任何复杂的计算就能猜到温度会相当低。在温度较高的加州的夏天,我们也可以猜到天气会相对炎热。这些粗略的估计可以用作贝叶斯模型中的初始猜测,实质上就是模型的先验知识。如果我们没有这些初始猜测,我们将需要进行更复杂的建模来进行预测。

正如我们在本章中学到的,将先验知识纳入 GP 中可以大大改变模型的行为,从而带来更好的预测性能(最终带来更有效的决策)。只有当我们对函数的行为没有任何好的猜测时,才应该不使用先验知识;否则,这等同于浪费信息。

在本章中,我们将讨论均值和协方差函数的不同选项,以及它们如何影响生成的 GP 模型。与第二章不同,我们在这里采用了一种实践的方法,并围绕 Python 中的代码实现展开讨论。到本章结束时,我们将开发一种选择适当的均值和协方差函数以及优化它们的超参数的流程。

3.1 先验概率对于贝叶斯模型的重要性

问题:为什么你似乎无法改变一些人的想法?答案:因为他们的先验概率。为了说明先验概率在贝叶斯模型中有多么重要,请考虑以下情景。

假设你和你的朋友鲍勃和爱丽丝在嘉年华上闲逛,你正在和一个声称自己是灵媒的人交谈。他们允许你通过以下程序来测试这一说法:你和你的朋友们每个人都想一个 0 到 9 之间的数字,而“灵媒”将告诉你们每个人在想什么数字。你可以重复这个过程任意次数。

现在,你们三个人都对这个所谓的灵媒感到好奇,你们决定进行这个测试 100 次。令人惊讶的是,在这 100 次测试之后,嘉年华上的这位自称的灵媒准确地猜出了你们每个人心里想的数字。然而,测试结束后,你们每个人的反应却各不相同,如图 3.1 所示。

图 3.1 显示了你们朋友中看到一个人连续猜对一个秘密数字 100 次后的反应。由于他们的先验信念不同,每个人得出了不同的结论。

关于贝叶斯定理的进一步阅读

如果你需要恢复记忆,请随时返回图 2.2,我们在那里研究了贝叶斯定理。在本书中,我们只是大致概述了这一过程,但我建议你阅读威尔·库尔特(Will Kurt)的《贝叶斯统计学的有趣方法》(No Starch Press,2019)的第 1 和第二章,如果你想深入研究这个过程。

你们三个人怎么可能观察同样的事件(嘉年华上的人连续猜对 100 次)却得出不同的结论呢?要回答这个问题,考虑使用贝叶斯定理更新信念的过程:

  1. 每个人都从特定的先验概率开始,认为这个人是个灵媒。

  2. 然后,你们每个人都观察到他们猜对你的数字一次的事件。

  3. 然后,你们每个人都计算可能性项。首先,鉴于他们确实是灵媒,他们的猜测正确的可能性是完全的 1,因为真正的灵媒总是能通过这个测试。其次,鉴于他们 不是 灵媒,他们的猜测正确的可能性是 10 个中的 1,因为每次,你们都是在 0 到 9 之间随机选择一个数字,所以这 10 个选项中的任何一个猜测都有相等的可能性:10 分之 1。

  4. 最后,通过将先验与这些可能性项相结合来计算这个人不是灵媒的后验概率,你们更新了自己的信念。具体来说,这个后验概率将与先验和第一个可能性项 相乘 成比例。

  5. 你们然后重复这个过程 100 次,每次使用前一次迭代的后验概率作为当前迭代的先验概率。

这里高层次的重要性在于,经过每次测试,你和你的朋友对这个人是灵媒的后验信念 从未减少,因为这个陈述与你观察到的数据不符。具体来说,图 3.2 显示了你们小组中每个人的渐进后验概率,作为“嘉年华中的灵媒”通过了多少次测试的函数。

图 3.2 进展后验概率,即在成功猜测次数的函数中,卡尼瓦尔上的女士是一个灵媒的概率。这个后验概率从不下降,但根据初始先验的不同而行为也不同。

如图所示,三条曲线中的每一条要么增加,要么保持不变——没有一条曲线实际下降,因为一个减少的成为灵媒的概率不符合 100 次连续成功猜测。但为什么这三条曲线看起来如此不同呢?你可能已经猜到了,曲线的起始位置——即每个人认为女士是灵媒的先验概率——是原因。

在鲍勃的案例中,在左侧面板上,他最初的先验相对较高,认为那个人是灵媒的概率为 1%。鲍勃是一个信徒。随着他观察到越来越多与这一信念一致的数据,他的后验概率也越来越高。

在你自己的情况下,在中间,作为怀疑论者,你的先验要低得多:10 的 14 次方中的 1。然而,由于你的观察确实表明女士是灵媒,随着更多数据的输入,你的后验概率也增加到最后的 1%。

另一方面,艾丽斯的情况则不同。从一开始,她不相信灵媒是真实存在的,因此她给了她的先验概率精确的零。现在,请记住,根据贝叶斯定理,后验概率与先验概率乘以似然性成比例。由于艾丽斯的先验概率恰好为零,贝叶斯更新中的这种乘法将总是产生另一个零。

由于艾丽斯最初的概率为零,在一次成功测试后,这个概率保持不变。一次正确的猜测后,艾丽斯的后验概率为零。两次猜测后,仍为零。所有 100 次正确的猜测后,这个数字仍然是零。所有的一切都符合贝叶斯更新规则,但由于艾丽斯的先验不允许灵媒存在的可能性,任何数据都无法说服她相反。

这突显了贝叶斯学习的一个重要方面——我们的先验确定了学习的方式(见图 3.3):

  • 鲍勃的先验相当高,因此在 100 次测试结束时,他完全相信那个人是灵媒。

  • 另一方面,你更加怀疑,你的初始先验比鲍勃低得多。这意味着你需要更多的证据才能得出高后验。

  • 艾丽斯对可能性的完全忽视,用她的零先验表示,使她的后验概率保持在零。

图 3.3 各人的先验信念如何被相同的数据更新。与鲍勃相比,你的先验更低,增长速度更慢。艾丽斯的先验为 0,始终保持为 0。

虽然我们的例子中的主张是关于某个人是否是通灵者的事件,但是同样的贝叶斯更新过程适用于所有情形,其中我们对某个事件有概率信念,并经常根据数据进行更新。事实上,这正是我们有时似乎无法改变某个人的想法的原因:因为他们的先验概率为零,没有任何东西可以将后验概率更新为非零。

从哲学角度来看,这个讨论非常有趣,因为它表明,为了能够说服某个人做某事,他们需要至少想到一种可能性,即指定事件的非零先验概率。更具体地说,这个例子说明了对贝叶斯模型拥有良好的先前知识的重要性。正如我们所说,我们通过均值和协方差函数来指定高斯过程的先验知识。每个选择都会在高斯过程的预测中产生不同的行为。

3.2 将先前已知内容并入高斯过程

在本节中,我们将确定在高斯过程中指定先前知识的重要性。这个讨论为我们在本章剩余部分的讨论下了动力。

先验高斯过程可能一开始具有恒定的平均值和 CI。如图 3.4 所示,该高斯过程然后被更新为平稳地内插观测到的数据点。也就是说,均值预测恰好穿过数据点,并且 95%的 CI 在那些区域消失。

图 3.4 先验高斯过程和后验高斯过程的比较。先验高斯过程包含有关目标函数的先前信息,而后验高斯过程将该信息与实际观测值相结合。

图 3.4 中的先验高斯过程不对我们正在建模的目标函数做出任何假设。这就是为什么这个高斯过程的平均预测值在任何地方都是零。但是,在许多情况下,即使我们不知道目标函数的确切形式,我们也了解目标的某些方面。

以以下为例:

  • 在超参数调整应用程序中模拟模型准确性时,我们知道目标函数的范围在 0 到 1 之间。

  • 在我们从第 2.1 节中的住房价格示例中,函数值(价格)严格为正,当房屋的理想属性(例如生活区)增加时,应该增加。

  • 在住房示例中,函数值对某些特征更加敏感。例如,与生活区面积的函数相比,房屋价格随着层数的增加更快——多一层楼会增加房屋价格,而多一个平方英尺的生活区则不会。

这种信息正是我们希望用 GP 表示的先验知识,使用 GP 的最大优势之一是我们有许多方法来融入先验知识。这样做有助于缩小 GP 代理与其建模的实际目标函数之间的差距,这也将更有效地引导后续的优化。

用 GP 结合先验知识

我们通过选择适当的平均和协方差函数以及设置它们的参数值来融入先验知识。特别是

  • 平均函数定义了目标函数的预期行为。

  • 协方差函数定义了目标的结构,或者更具体地说,定义了任意一对数据点之间的关系,以及目标函数在其定义域内变化的速度和平滑程度。

前面的每个选择都会导致生成的 GP 的行为发生 drastical 不同。例如,线性平均函数将导致 GP 预测中的线性行为,而二次平均函数将导致二次行为。通过在协方差函数中使用不同的参数,我们还可以控制 GP 的变异性。

3.3 使用平均函数定义函数行为

首先,我们介绍了 GP 的平均函数,它定义了 GP 的预期行为,或者说我们相信函数在所有可能的情况下平均情况下的样子。正如我们将要看到的,这有助于我们指定与函数的一般行为和形状相关的任何先验知识。本节中使用的代码包含在 CH03/01 - Mean functions.ipynb 中。为了使我们的讨论具体化,我们使用了一个房价数据集,其中表 3.1 中有五个数据点。

表 3.1 示例训练数据集。预测目标(价格)随特征(居住面积)的增加而增加。

居住面积(以 1000 平方英尺为单位)价格(以 10 万美元为单位)
0.50.0625
10.25
1.50.375
32.25
44

在这个数据集中,我们建模的函数值是房价,它们是严格正的,并且随着居住面积的增加而增加。这些性质是直观的,即使不知道未观察到的房屋的价格,我们也确信这些未见价格也具有这些性质。

我们的目标是将这些性质纳入我们的平均函数中,因为它们描述了我们对函数行为的预期。在我们开始建模之前,我们首先编写一个辅助函数,该函数接受一个 GP 模型(以及其似然函数),并在范围从 0 到 10(即,一个 10,000 平方英尺的居住面积)内可视化其预测。实现如下:

def visualize_gp_belief(model, likelihood):
    with torch.no_grad():                                ❶
        predictive_distribution = likelihood(model(xs))  ❶
        predictive_mean = predictive_distribution.mean   ❶
        predictive_upper, predictive_lower =             ❶
        ➥ predictive_distribution .confidence_region()  ❶

    plt.figure(figsize=(8, 6))

    plt.plot(xs, ys, label="objective", c="r")
    plt.scatter(train_x, train_y, marker="x", c="k", label="observations")

    plt.plot(xs, predictive_mean, label="mean")          ❷
    plt.fill_between(
        xs.flatten(), predictive_upper, predictive_lower, alpha=0.3,
        label="95% CI"
    )                                                    ❸

    plt.legend(fontsize=15);

❶ 计算预测

❷ 绘制平均线

❸ 绘制 95% CI 区域

我们在第 2.4.3 节中看到了这段代码是如何工作的,现在我们将其放入一个方便的函数中。有了这个,我们就准备好实现我们的 GP 模型,并看看我们的选择如何影响所产生的预测。

3.3.1 使用零均值函数作为基本策略

均值的最简单形式是一个在零处的常数函数。在没有数据的情况下,此函数将产生零作为其默认预测。当没有关于我们可能将其作为先验知识合并到 GP 中的目标函数的额外信息时,将使用零均值函数。

使用零均值函数实现的 GP 如下所示:

class ConstantMeanGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
       super().__init__(train_x, train_y, likelihood)
       self.mean_module = gpytorch.means.ConstantMean()   ❶
       self.covar_module = gpytorch.kernels.RBFKernel()

    def forward(self, x):
       mean_x = self.mean_module(x)
       covar_x = self.covar_module(x)
       return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

❶ 默认值为零的常数均值函数

请记住,根据第 2.4.2 节的内容,要构建一个使用 GPyTorch 的 GP 模型,我们实现 __init__()forward() 方法。在第一个方法中,我们初始化我们的均值和协方差函数;在第二个方法中,我们将输入 x 通过这些函数,并返回相应的多元高斯分布。

注意:在我们从第 2.4.2 节中的实现中,我们使用 gpytorch.means.ZeroMean 类,而在这里,我们使用 gpytorch.means.ConstantMean 类来初始化我们的均值函数。然而,这个常数均值函数的默认值是零,因此实际上,我们仍然在实现相同的 GP 模型。尽管这两种选择目前导致相同的模型,但在本章中,我们将展示如何使用 gpytorch.means.ConstantMean 来调整常数均值值,以获得更好的模型。

现在让我们初始化该类的一个对象,在我们的训练数据上对其进行训练,并可视化其预测。我们用以下代码来实现这个:

lengthscale = 1
noise = 1e-4

likelihood = gpytorch.likelihoods.GaussianLikelihood()     ❶
model = ConstantMeanGPModel(train_x, train_y, likelihood)  ❶

model.covar_module.lengthscale = lengthscale               ❷
model.likelihood.noise = noise                             ❷

model.eval()
likelihood.eval()

visualize_gp_belief(model, likelihood)

❶ 声明 GP

❷ 修复超参数

在这里,我们初始化 GP 模型并设置其超参数——长度尺度和噪声方差分别为 1 和 0.0001。我们将在本章后面看到如何适当设置这些超参数的值;现在,让我们只使用这些值。最后,我们在我们的 GP 模型上调用我们刚刚编写的辅助函数 visualize_gp_belief(),它将生成图 3.5。

所有我们在第 2.4.4 节中指出的 GP 的良好性质仍然存在:

  • 后验均值函数平滑地插值出我们的训练数据点。

  • 95% CI 在这些数据点周围消失,表示了一个良好校准的不确定性量化。

图 3.5 使用零均值函数的 GP 的预测。后验均值函数插值了观察到的数据点,并在远离这些观测点的区域回归到零。

从这个图表中我们还注意到,一旦我们足够远离我们的训练数据点(图表的右侧),我们的后验均值函数就会恢复到先验均值,即零。这实际上是高斯过程的一个重要特征:在没有数据的情况下(在没有观测到的区域),先验均值函数是推断过程的主要驱动力。这在直觉上是有道理的,因为没有实际观察,预测模型能做的最好的事情就是简单地依赖于其均值函数中编码的先验知识。

注意:在这一点上,我们看到为什么将先验知识明确定义地编码到先验高斯过程中是如此重要:在没有数据的情况下,预测的唯一驱动因素就是先验高斯过程。

然后自然地提出了一个问题:是否可以使用非零均值函数来诱导我们的高斯过程在这些未探索的区域中产生不同的行为,如果可以,我们有哪些选择?本节的剩余部分旨在回答这个问题。我们首先使用一个不为零的常数均值函数。

使用梯度下降法使用常数函数 3.3.2

如果一个常数均值函数不为零,那么当我们期望我们正在建模的目标函数具有我们先验知道的一些值范围时,这是合适的。由于我们正在建模房价,使用一个常数均值函数,其常数大于零,是有意义的,因为我们确实期望价格是正的。

当然,在许多情况下,我们不可能知道目标函数的平均值是多少。那么,我们应该如何为我们的均值函数找到一个合适的值呢?我们使用的策略是依赖于一个特定的数量:给定我们均值函数的值时,训练数据集有多大可能性。粗略地说,这个数量衡量了我们的模型解释其训练数据的能力。我们展示了如何在这一小节中使用这个数量来选择我们的高斯过程的最佳均值函数。

如果给定一些值c[1]时训练数据的似然性高于给定另一个值c[2]时的情况,那么我们更喜欢使用c[1]而不是使用c[2]。这量化了我们先前关于使用非零均值函数来建模正函数的直觉:一个常数均值函数,其值为正,比值为零(或负值)的函数更好地解释了来自完全正函数的观测。

我们如何计算这个似然性呢?GPyTorch 提供了一个方便的类,gpytorch.mlls.ExactMarginalLogLikelihood,它接受一个高斯过程模型并计算其训练数据的边际对数似然性,给定模型的超参数。

要看到这个似然数量在量化数据拟合方面的有效性,考虑图 3.6。这个图可视化了两个不同 GP 模型的预测:一个我们在前面子段中看到的零均值 GP,在左边,以及均值函数值为 2 的 GP,在右边。注意在第二个面板中,均值函数在图的右侧恢复到 2 而不是 0。在这里,第二个 GP 的(对数)似然性比第一个高,这意味着值 2 比值 0 更好地解释了我们的训练数据。

图 3.6 GP 预测,给出两个不同常数均值函数。值 2 比值 0 给出了更高的似然值,表明前者的均值函数比后者更合适。

有了这个对数似然计算,我们的最后一步就是简单地找到我们的均值函数的值,使得对数似然被最大化。换句话说,我们的目标是寻找最能解释我们训练数据的均值。由于我们可以访问对数似然计算,我们可以使用基于梯度的优化算法,例如梯度下降,来迭代地优化我们的均值。当收敛时,我们将得到一个很好的均值,它给出了高数据似然度。如果您需要恢复一下梯度下降的工作原理,我推荐 Luis Serrano 的《Grokking Machine Learning》(Manning,2021)的附录 B,它很好地解释了这个概念。

现在,让我们看看如何用代码实现这个过程。由于我们用gpytorch.means.ConstantMean类为我们的均值函数实现了 GP 模型,所以我们这里不需要做任何更改。所以,现在,让我们再次初始化我们的 GP 模型:

# declare the GP
lengthscale = 1
noise = 1e-4

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ConstantMeanGPModel(train_x, train_y, likelihood)

# fix the hyperparameters
model.covar_module.lengthscale = lengthscale
model.likelihood.noise = noise

这个过程的核心步骤是定义对数似然函数以及梯度下降算法。如前所述,前者是gpytorch.mlls.ExactMarginalLogLikelihood类的一个实例,实现如下:

mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

对于梯度下降算法,我们使用 Adam,这是一种在许多 ML 任务中取得了很大成功的最先进算法,尤其是 DL。我们用 PyTorch 声明如下:

optimizer = torch.optim.Adam([model.mean_module.constant], lr=0.01)

请注意,我们正在传递给torch.optim.Adam类的是model.mean_module.constant,这是我们希望优化的均值。当我们运行梯度下降过程时,Adam 算法会迭代地更新model.mean_module.constant的值,以改善似然函数。

现在我们需要做的最后一件事是运行梯度下降,其实现如下:

model.train()                         ❶
likelihood.train()                    ❶

losses = []
constants = []
for i in tqdm(range(500)):
    optimizer.zero_grad()

    output = model(train_x)
    loss = -mll(output, train_y)      ❷

    loss.backward()                   ❸

    losses.append(loss.item())
    constants.append(model.mean_module.constant.item())

    optimizer.step()                  ❸

model.eval()                          ❹
likelihood.eval()                     ❹

❶ 启用训练模式

❷ 损失作为负边缘对数似然

❸ 损失上的梯度下降

❹ 启用预测模式

开始时的train()调用和最后的eval()调用是我们总是需要进行的记录步骤,分别启用我们 GP 模型的训练模式和预测模式。每一步都使用optimizer.zero_grad()重置梯度是另一个记录任务,以确保我们不会错误地计算梯度。

在中间,我们有一个 500 步的梯度下降过程,通过迭代计算损失(即我们的负对数似然)并根据梯度下降来降低这个损失。在这个 for循环过程中,我们跟踪所获得的负对数似然值以及在每一步中调整的平均值。这是为了在训练之后,我们可以可视化检查这些值以确定它们是否收敛。

图 3.7 展示了这种可视化,显示了我们希望最小化的负对数似然的运行值以及 GP 的均值函数值。我们的损失一直降低(这是件好事!)随着均值常数的增加,显示出正常数比零更有可能性。两条曲线在 500 次迭代后都趋于平稳,表明我们已经收敛到均值常数的最优值。

图 3.7 展示了梯度下降过程中的负对数似然(较低的值为更好)和均值值。在这两个面板中,这些值已经收敛,表明我们已经到达了最优解。

注意,在使用梯度下降时,我建议您始终绘制出渐进损失的图像,就像我们刚才在这里所做的那样,以查看您是否已经收敛到最优值。在收敛之前停止可能会导致模型性能不佳。虽然在本章节中我们不会再展示这些渐进损失图,但附带的代码中包含了它们。

到目前为止,我们已经学会了在我们的 GP 模型中使用零均值函数作为默认值,并优化相对于数据似然的平均常数值。然而,在许多使用情况下,您可能对目标函数的行为有先验知识,因此更喜欢将更多结构纳入您的均值函数中。

例如,我们如何实现房屋价格随着更大的居住面积而增加的想法?向前迈进,我们学会了使用线性或二次均值函数来实现这一点。

3.3.3 使用带梯度下降的线性函数

我们继续使用线性均值函数,其形式为μ = w^T**x + b。这里,μ是测试点x上的预测均值,而w是连接x中每个特征的系数的权重向量,b是一个常数偏置项。

通过使用线性平均函数,我们对我们的目标函数的预期行为进行编码,即它等于数据点x的特征的线性组合。对于我们的住房价格示例,我们只有一个特征,即生活区域,我们期望它有一个正权重,因此通过增加生活区域,我们的模型将预测出更高的价格。

另一种思考线性平均模型的方式是,我们有一个线性回归模型(也假设目标标签为特征的线性组合),然后我们在预测上加上一个概率信念,一个 GP 模型。这给予我们线性回归模型的能力,同时保持使用 GP 建模的所有好处,即不确定性量化。

注意 在一个常数均值函数下,权重向量w固定为零向量,偏差b是我们学会在上一小节优化的平均值。换句话说,线性函数比常数均值函数是一个更一般的模型。

关于实施,使用具有线性平均函数的 GP 模型相当简单。我们只需用gpytorch.means.LinearMean实例取代我们的常数均值,这样(我们的forward()方法保持不变):

class LinearMeanGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.LinearMean(1)    ❶
        self.covar_module = gpytorch.kernels.RBFKernel()

❶ 线性平均

在这里,我们使用1来初始化我们的均值模块,表示我们正在处理一维目标函数。当处理高维函数时,你可以简单地在这里传递该函数的维度。除此之外,我们模型的其它部分与之前的相似。在我们的三个点数据集上拟合和训练这个新模型,我们得到图 3.8 中的预测。

图 3.8 GP 带有线性平均函数的预测。GP 有一个上升趋势,这直接是线性平均函数的正斜率的结果。

与我们到目前为止看到的恒定平均值不同,我们在这里使用的线性平均函数驱使整个 GP 模型具有上升趋势。这是因为我们在训练数据中的五个数据点的最佳拟合线是具有正斜率的线,这正是我们希望作为生活区域和价格之间关系的建模。

3.3.4 通过实施自定义均值函数使用二次函数

我们这里的线性平均函数成功捕捉到了价格的上涨趋势,但它假定价格增长率是恒定的。也就是说,对生活区域增加一个额外的平方英尺,预期上会导致价格的恒定增加。

然而,在许多情况下,我们可能有先验知识,即我们的目标函数以非恒定的速率增长,而线性均值无法建模。事实上,我们使用的数据点是这样生成的,即价格是生活区域的二次函数。这就是为什么我们看到较大的房子比较小的房子更快地变得更贵。在本小节中,我们将将我们的 GP 均值实现为一个二次函数。

在撰写本文时,GPyTorch 仅提供了常数和线性均值函数的实现。但是,正如我们将在本书中一次又一次地看到的那样,这个软件包的美妙之处在于其模块化:GP 模型的所有组件,如均值函数、协方差函数、预测策略,甚至边缘对数似然函数,都被实现为模块,因此,它们可以以面向对象的方式进行修改、重新实现和扩展。当我们实现自己的二次均值函数时,我们首先亲自体会到这一点。

我们要做的第一件事是定义一个均值函数类:

class QuadraticMean(gpytorch.means.Mean):
    def __init__(self, batch_shape=torch.Size(), bias=True):
        ...

    def forward(self, x):
        ...

这个类扩展了 gpytorch.means.Mean 类,它是所有 GPyTorch 均值函数实现的基类。为了实现我们的自定义逻辑,我们需要重新编写两个方法:__init__()forward(),这与我们实现 GP 模型时完全相同!

__init__() 中,我们需要声明我们的均值函数包含哪些参数。这个过程称为 参数注册

虽然线性函数有两个参数,斜率和截距,而二次函数有三个参数:一个二阶项 x[2] 的系数;一个一阶项 x 的系数;和一个零阶项的系数,通常称为偏差。这在图 3.9 中有所说明。

图 3.9 线性函数和二次函数的函数形式。线性函数有两个参数,而二次函数有三个。当这些函数用作 GP 的均值函数时,相应的参数是 GP 的超参数。

在这个基础上,我们这样实现 __init__() 方法:

class QuadraticMean(gpytorch.means.Mean):
    def __init__(self, batch_shape=torch.Size(), bias=True):
        super().__init__()
        self.register_parameter(
            name="second",
            parameter=torch.nn.Parameter(torch.randn(*batch_shape, 1, 1))
        )                                                                 ❶

        self.register_parameter(
            name="first",
            parameter=torch.nn.Parameter(torch.randn(*batch_shape, 1, 1))
        )                                                                 ❷

        if bias:
            self.register_parameter(
                name="bias",
                parameter=torch.nn.Parameter(torch.randn(*batch_shape, 1))
            )                                                             ❸
        else:
            self.bias = None

❶ 二阶系数

❷ 一阶系数

❸ 偏差

我们按顺序调用 register_parameter() 来注册二阶系数、一阶系数和偏差。由于我们不清楚这些系数应该取什么值,我们只需使用 torch.randn() 随机初始化它们。

注意 我们需要将这些参数注册为 torch.nn .Parameter 类的实例,这允许它们的值在梯度下降期间进行调整(训练)。

对于 forward() 方法,我们需要定义我们的均值函数如何处理输入。正如我们所说的,二次函数的形式为 ax² + bx + c,其中 abc 分别是二阶系数、一阶系数和偏差。所以我们只需要实现这个逻辑,如下所示:

class QuadraticMean(gpytorch.means.Mean):
    def __init__(self, train_x, train_y, likelihood):
        ...                                               ❶

    def forward(self, x):
        res = x.pow(2).matmul(self.second).squeeze(-1) \
            + x.matmul(self.first).squeeze(-1)            ❷
        if self.bias is not None:
            res = res + self.bias
        return res

❶ 已省略

❷ 二次函数的公式

有了这个二次均值函数,我们现在可以编写一个初始化其均值模块的高斯过程模型,使用的是我们刚刚实现的自定义 QuadraticMean 类:

class QuadraticMeanGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = QuadraticMean()
        self.covar_module = gpytorch.kernels.RBFKernel()

    def forward(self, x):
        ...                 ❶

❶ 省略

通过梯度下降重新运行整个训练过程,我们得到了图 3.10 中的预测。

图 3.10 由具有二次均值函数的高斯过程进行的预测。该高斯过程预测随着较大的居住面积,价格增长速度更快。

在这里,我们成功地模拟了房屋价格相对于居住面积的非恒定增长率。我们的预测在图的右侧增长得比左侧快得多。

我们可以对我们希望假设的任何功能形式进行这样的操作,例如更高阶的多项式或次线性函数。我们需要做的是实现一个具有适当参数的均值函数类,并使用梯度下降为这些参数分配值,以便为我们的训练数据提供良好的拟合。

到目前为止,我们的讨论展示了高斯过程模型的数学灵活性,即它们可以利用任何结构的均值函数,并且仍然能够产生概率预测。这种灵活性激发并推动了 GPyTorch 的设计,其对模块化的强调帮助我们轻松地扩展和实现自己的自定义均值函数。在接下来我们将讨论的 GPyTorch 的协方差函数中,我们看到了同样的灵活性和模块化。

3.4 用协方差函数定义变异性和平滑度

虽然均值函数定义了我们对目标函数整体行为的期望,但高斯过程的协方差函数或核函数扮演着更复杂的角色:定义了域内数据点之间的关系,并控制了高斯过程的结构和平滑度。在本节中,我们比较了高斯过程对模型不同组成部分进行更改时的预测。通过这样做,我们能够实际洞察如何为高斯过程模型选择适当的协方差函数。我们使用的代码在 CH03/02 - Covariance functions.ipynb 中。

在这些示例中,我们使用了在第 2.4.1 节中看到的 Forrester 函数作为我们的目标函数。我们再次在 -3 到 3 之间随机采样三个数据点,并将它们用作我们的训练数据集。本节中可视化的所有预测都来自于在这三个点上训练的高斯过程。

3.4.1 设置协方差函数的尺度

通过其协方差函数来控制高斯过程行为的第一种方法是设置长度尺度和输出尺度。这些尺度,就像均值函数中的常数或系数一样,是协方差函数的超参数:

  • 长度尺度控制着高斯过程输入的尺度,因此,高斯过程沿轴变化的速度有多快—也就是说,我们相信目标函数在输入维度上的变化程度有多大。

  • 输出尺度定义了 GP 的输出范围或者说是其预测范围。

通过设置这些尺度的不同值,我们可以增加或减少 GP 预测的不确定性,以及缩放我们的预测范围。我们使用以下代码实现:

class ScaleGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ZeroMean()
        self.covar_module =
        ➥ gpytorch.kernels.ScaleKernel(
          gpytorch.kernels.RBFKernel())     ❶

    def forward(self, x):
        ...                                 ❷

❶ gpytorch.kernels.ScaleKernel 实现了输出尺度。

❷ 忽略

注意,这里的 covar_module 属性代码与之前不同:我们将一个 gpytorch.kernels.ScaleKernel 对象放在了通常的 RBF 内核之外。这实际上实现了输出尺度,它通过某个定值因子对 RBF 内核的输出进行了缩放。而长度尺度则已经包含在了 gpytorch.kernels.RBFKernel 内核中。

在我们迄今为止使用的代码中,我们有一行代码用于设置我们内核的长度尺度,即 model.covar_module.base_kernel.lengthscale = lengthscale。这就是长度尺度的值所保存的位置。使用类似的 API,我们可以使用 model.covar_module.outputscale = outputscale 来设置内核的输出尺度。现在,为了验证长度尺度确实可以控制函数变化的速度,我们将比较两个 GP 的预测,一个长度尺度为 1,另一个为 0.3,如图 3.11 所示。

图 3.11 显示了 GP 以长度尺度为 1(左边)和 0.3(右边)进行预测的结果。长度尺度越小,GP 的预测就越不确定,变异性越高。

这两个面板的巨大差异清晰地表明了长度尺度的效果:

  • 一个较短的长度尺度相当于给定输入常量变化情况下客观函数更多的变异性。

  • 相反,一个较长的长度尺度强制函数更加光滑,也就是说给定输入更少的变化会使其变异性减少。

例如,沿着 x 轴移动一个单位,在图 3.11 左侧面板中的样本的变化幅度小于右侧面板的样本。

那么输出尺度呢?我们之前说这个参数将协方差函数的输出值缩放到一个不同的范围。这是通过将协方差输出结果乘以这个参数来实现的。因此,较大的输出尺度会使 GP 的预测范围更宽,而较小的输出尺度则会使预测范围缩小。为了验证这一点,让我们再次运行代码并重新生成预测,这次将输出尺度设置为 3。生成的输出结果如图 3.12 所示。

图 3.12 显示了 GP 的输出尺度为 3 的预测结果。较大的输出尺度使 GP 模拟的函数范围更宽,也容许更多的预测不确定性。

尽管图 3.11 和图 3.12 的左侧面板看起来相同,因为 GP 及其样本在两个图中具有相同的形状,但我们注意到图 3.12 的 y-轴更大,因为它的预测值和样本值都取得了较大的值(负值和正值都有)。这是使用具有较大输出尺度的 RBF 核对协方差值进行缩放的直接结果。

通过我们的协方差函数的两个超参数,我们已经看到我们可以解释由我们的 GP 模型建模的各种功能行为,我们在表 3.2 中总结了这些行为。我邀请您使用不同的长度和输出尺度重新运行此代码,以查看其效果并验证表格!

表 3.2 GP 的长度和输出尺度的角色总结

参数大值小值
长度尺度更平滑的预测,较少的不确定性更多的变异性,更多的不确定性
输出尺度较大的输出值,更多的不确定性较窄的输出范围,较少的不确定性

注意 这种建模的灵活性引发了一个自然的问题:你应该如何适当地设置这些超参数的值?幸运的是,我们已经知道了一种设置 GP 模型超参数的好方法。我们可以通过选择最好地解释我们的数据的值,或者换句话说,通过梯度下降来最大化边缘对数似然,从而实现这一目标。

就像我们希望优化均值函数的超参数一样,我们现在只需将我们希望优化的变量 - 协方差函数的参数 - 传递给 Adam 即可实现:

optimizer = torch.optim.Adam(model.covar_module.parameters(), lr=0.01)

通过运行梯度下降,我们可以获得这些参数的良好值。具体来说,我获得了大约 1.3 的长度尺度和大约 2.1 的输出尺度。也就是说,为了很好地拟合这三点训练数据集,我们希望 GP 稍微更平滑一些(具有大于 1 的长度尺度),并且还希望我们的预测范围更大一些(具有较大的输出尺度)。这无疑是一个让人放心的结果,因为我们的目标函数具有广泛的数值范围 - 在输入 3 处,它的值将达到 -2,这远超出了具有输出尺度 1 的 CI。

3.4.2 使用不同的协方差函数控制平滑度

到目前为止,我们仅仅使用了 RBF 核作为我们的协方差函数。然而,如果 RBF 不合适,完全可以使用不同的核函数用于我们的 GP。在本小节中,我们将学习使用另一种核函数家族,即马特恩核,并看看这种核函数对我们的 GP 会产生什么影响。

注意 通过使用马特恩核,我们正在指定 GP 模型函数的平滑度。这里的 "平滑度" 是一个技术术语,指的是函数的可微性;函数可微性越多次,它就越平滑。我们可以大致将其视为函数值以曲折方式 "跳动" 的程度。

RBF 核模拟的函数具有无限可微性,这是现实世界中很少有的函数特性。与此同时,Matérn 核生成的函数是有限可微的,这些函数可以被微分的次数(即这些函数的平滑度)由可设置的参数控制,我们马上就会讨论到。

要看到 Matérn 核的实际效果,我们首先重新实现我们的 GP 模型类:

class MaternGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood, nu):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ZeroMean()
        self.covar_module = gpytorch.kernels.MaternKernel(nu)

    def forward(self, x):
        ...                  ❶

❶ 省略

在这里,我们的 covar_module 属性被初始化为 gpytorch .kernels.MaternKernel 类的实例。这个初始化接受一个参数 nu,定义了我们的 GP 将具有的平滑程度,这也是我们 __init__() 方法的一个参数。

重要提示:在撰写本文时,GPyTorch 支持三个 nu 值,1/2、3/2 和 5/2,对应函数分别是不可微分、一次可微分和两次可微分的。换句话说,这个 nu 参数越大,我们的 GP 越平滑。

让我们首先尝试 nu = 0.5,通过在初始化 GP 时设置该值来实现:

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = MaternGPModel(train_x, train_y, likelihood, 0.5)

...     ❶

visualize_gp_belief(model, likelihood)

❶ 修正超参数并启用预测模式

这段代码生成图 3.13。

图 3.13 中 Matérn 1/2 核的 GP 预测,这表明目标函数不可微分,对应非常粗糙的样本

与我们之前在 RBF 中看到的情况不同,这个 Matérn 核的样本都非常参差不齐。事实上,它们都不可微分。当建模时间序列数据时,比如股票价格,nu = 0.5 是 Matérn 核的一个好值。

然而,在 BayesOpt 中通常不使用这个值,因为像图 3.13 中那样的参差不齐的函数非常不稳定(它们可以以不可预测的方式上下跳动),通常不是自动优化技术的目标。我们需要目标函数具有一定的平滑度,以便进行优化;否则,有效的优化是不切实际的目标。

Matérn 5/2 核通常是首选。它的预测结果与 Matérn 3/2 生成的结果在图 3.14 中可视化。

图 3.14 中 Matérn 5/2(左)和 Matérn 3/2(右)核的 GP 预测。这里的样本足够平滑,以便 GP 有效地从数据中学习,但也足够参差不齐,以真实地模拟现实生活中的过程。

我们看到这个 5/2 核的样本要平滑得多,这导致 GP 更有效地学习。然而,这些样本也足够粗糙,以至于它们类似于我们在现实世界中可能看到的函数。因此,BayesOpt 中的大多数工作,无论是研究还是应用,都使用这个 Matérn 5/2 核。在未来的章节中,当我们讨论 BayesOpt 的决策时,我们将相应地默认使用这个核。

注意 虽然我们在这里没有包含相应的细节,但 Matérn 核函数有自己的长度尺度和输出尺度,可以像前面的小节一样指定,以进一步定制生成的 GP 的行为。

通过将均值函数与核函数配对,我们可以在 GP 的预测中诱发复杂的行为。就像我们的先验影响了我们朋友圈中每个人在看到有人正确猜出一个秘密数字 100 次后的结论一样,我们对均值函数和核函数的选择决定了 GP 的预测。图 3.15 展示了三个例子,其中每种均值函数和核函数的组合导致了截然不同的行为。

图 3.15 展示了在相同数据集上训练时,均值函数和核函数的三种不同选择以及它们各自后验 GP 的预测,每种选择都导致不同的预测行为。

3.4.3 用多个长度尺度建模不同级别的变异性

因为我们一直只考虑一维目标函数(其输入具有一个特征),所以我们只需要考虑一个长度尺度。然而,我们可以想象到一种情况,在这种情况下,高维目标函数(其输入具有多个特征)在某些维度上具有更多的变异性,在其他维度上较为平滑。也就是说,某些维度具有小的长度尺度,而其他维度具有大的长度尺度。还记得本章开头的启发性例子吗:房屋价格的预测增加一个楼层的幅度比增加一个平方英尺的生活面积更大。在本小节中,我们探讨了如何在 GP 中维护多个长度尺度以对这些函数进行建模。

如果我们只使用一个长度尺度来处理所有维度,那么我们将无法忠实地对目标函数进行建模。这种情况需要 GP 模型为每个维度维护一个单独的长度尺度,以完全捕获它们各自的变异性。在本章的最后一节中,我们学习如何在 GPyTorch 中实现这一点。

为了帮助我们的讨论,我们使用一个具体的二维目标函数,称为 Ackley,它可以修改为在不同维度中具有各种级别的变异性。我们将该函数实现如下:

def ackley(x):
    # a modification of https:/ /www.sfu.ca/~ssurjano/ackley.xhtml
    return -20 * torch.exp(
        -0.2 * torch.sqrt((x[:, 0] ** 2 + x[:, 1] ** 2) / 2)
    )
    ➥ - torch.exp(torch.cos(2 * pi * x[:, 0] / 3)
    ➥ + torch.cos(2 * pi * x[:, 1]))

我们特别将该函数的定义域限制为两个维度上的方形区域,即 -3 到 3,通常表示为 [-3, 3]²。为了可视化这个目标函数,我们在图 3.16 中使用热图。

图 3.16 作为我们的目标使用的二维 Ackley 函数。在这里,x 轴的变异性比 y 轴小(变化较少),需要不同的长度尺度。

热图中的每个暗斑点都可以被看作是目标函数表面上具有低值的山谷。在这里,沿y轴有更多的山谷,而沿x轴则没有那么多,这表明第二个维度的变异性更大——也就是说,目标函数沿y轴上下波动的次数比沿x轴多得多。

再次强调,这意味着仅使用一个长度尺度来描述两个维度并不是一个好选择。相反,我们应该为每个维度(在本例中为两个)有一个长度尺度。然后,可以使用梯度下降独立地优化每个长度尺度。

使用每个维度的长度尺度的内核称为自动相关确定(ARD)。这个术语表示,在使用梯度下降优化这些长度尺度之后,我们可以推断出目标函数的每个维度与函数值相关的程度。具有较大长度尺度的维度具有较低的变异性,因此,在建模目标函数值时比具有较小长度尺度的维度不太相关。

使用 GPyTorch 实现 ARD 非常容易:我们只需在初始化协方差函数时将 ard_num_dims 参数指定为我们的目标函数具有的维度数即可。像这样使用 RBF 内核:

class ARDGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ZeroMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel(ard_num_dims=2)
        )

    def forward(self, x):
        ...                 ❶

❶ 省略

让我们看看,当在我们的 Ackley 函数上训练时,这个模型是否为两个维度给出了不同的长度尺度。为此,我们首先构造一个包含 100 个点的随机抽样训练数据集:

torch.manual_seed(0)
train_x = torch.rand(size=(100, 2)) * 6 - 3
train_y = ackley(train_x)

在使用梯度下降训练模型后,我们可以通过打印出优化后的长度尺度的值来检查它们。

>>> model.covar_module.base_kernel.lengthscale
tensor([[0.7175, 0.4117]])

这确实是我们预期的结果:第一个维度的长度尺度较大,函数值的变异性较低,而第二个维度的长度尺度较小,函数值的变异性较高。

更多关于内核的阅读

内核本身已经受到了 ML 社区的极大关注。除了我们迄今所涵盖的内容之外,还有一点需要注意的是,内核还可以编码复杂的结构,如周期性、线性和噪声。对于更全面和技术性的内核讨论,感兴趣的读者可以参考 David Duvenaud 的 内核菜谱 (www.cs.toronto.edu/~duvenaud/cookbook/)。

这个讨论标志着第三章的结束。在本章中,我们广泛地研究了我们的 GP 模型是如何受到均值和协方差函数的影响,特别是它们的各个参数。我们将这视为将我们对目标函数的了解(即先验信息)融入到我们的 GP 模型中的一种方式。我们还学会了使用梯度下降来估计这些参数的值,以获得最佳解释我们数据的 GP 模型。

这也标志着本书的第一部分的结束,我们将重点放在了 GP 上。从下一章开始,我们开始学习 BayesOpt 框架的第二个组成部分:决策制定。我们从两种最常用的 BayesOpt 策略开始,这两种策略旨在改善已见到的最佳点:概率改进和期望改进。

3.5 练习

这个练习是为了练习使用 ARD 实现 GP 模型。为此,我们创建一个目标函数,沿一个轴变化的程度比沿另一个轴变化的程度更大。然后,我们使用来自该函数的数据点训练一个有或没有 ARD 的 GP 模型,并比较学习到的长度尺度值。解决方案包含在 CH03/03 - Exercise.ipynb 中。

这个过程有多个步骤:

  1. 使用 PyTorch 在 Python 中实现以下二维函数:

    这个函数模拟了一个支持向量机(SVM)模型在超参数调整任务中的准确性曲面。x轴表示惩罚参数c的值,y轴表示 RBF 核参数γ的值。(我们在未来的章节中也将使用该函数作为我们的目标函数。)

  2. 在区间[0, 2]²上可视化该函数。热图应该看起来像图 3.17。

    图 3.17 SVM 模型在测试数据集上的准确率作为惩罚参数c和 RBF 核参数γ的函数。函数对γ的变化比对c的变化更快。

  3. 在区间[0, 2]²中随机选择 100 个数据点,作为我们的训练数据。

  4. 使用常数均值函数和 Matérn 5/2 卷积核来实现一个 GP 模型,其中输出规模作为gpytorch.kernels.ScaleKernel对象实现。

  5. 在初始化核对象时不要指定ard_num_dims参数,或将该参数设置为None。这将创建一个没有 ARD 的 GP 模型。

  6. 使用梯度下降训练 GP 模型的超参数,并在训练后检查长度尺度。

  7. 重新定义 GP 模型类,这次设置ard_num_dims = 2。使用梯度下降重新训练 GP 模型,并验证两个长度尺度具有显著不同的值。

总结

  • 先验知识在贝叶斯模型中起着重要作用,可以极大地影响模型的后验预测结果。

  • 使用 GP 模型可以通过均值和协方差函数来指定先验知识。

  • 均值函数描述了高斯过程模型的预期行为。在没有数据的情况下,高斯过程的后验均值预测回归到先验均值。

  • 高斯过程的均值函数可以采用任何函数形式,包括常数、线性函数和二次函数,这可以通过 GPyTorch 实现。

  • 高斯过程的协方差函数控制了高斯过程模型的平滑度。

  • 长度尺度指定了输出与函数输入之间的变异性水平。较大的长度尺度导致更加平滑,因此预测的不确定性较小。

  • 高斯过程中的每个维度都可以有自己的长度尺度。这被称为自动相关性确定(ARD),用于模拟在不同维度上具有不同变异程度的目标函数。

  • 输出尺度指定了函数输出的范围。较大的输出尺度导致更大的输出范围,因此在预测中有更多的不确定性。

  • Matérn 核类是 RBF 核类的泛化。通过指定其参数 nu,我们可以模拟高斯过程预测中的各种平滑程度。

  • 高斯过程的超参数可以通过最大化使用梯度下降的数据的边际似然来进行优化。

第二部分:用贝叶斯优化做决策

GP 只是方程式的一部分。为了充分实现 BayesOpt 技术,我们需要方程式的第二部分:决策策略,这些策略规定了如何进行函数评估以尽快优化目标函数。本部分列举了最流行的 BayesOpt 策略,包括它们的动机、数学直觉和实现。虽然不同的策略受到不同目标的驱使,但它们都旨在平衡勘探和开发之间的权衡——这是 BayesOpt 特别是以及不确定性问题下的决策制定的核心挑战。

第四章通过引入获取分数的概念来开始介绍事物,这是一种量化进行函数评估价值的方法。该章还描述了一种试图从迄今为止看到的最佳点改进的启发式方法,这导致了两种流行的 BayesOpt 策略:改进概率和预期改进。

第五章将 BayesOpt 与一个密切相关的问题联系起来:多臂赌博机。我们探索了流行的上置信界限策略,该策略使用乐观主义下的不确定性启发式,以及 Thompson 抽样策略,该策略利用了 GP 的概率性质来辅助决策。

第六章向我们介绍了信息理论,这是数学的一个子领域,在决策问题中有许多应用。利用信息理论的一个核心概念熵,我们设计了一个名为 BayesOpt 策略,该策略旨在获取关于我们搜索目标的最多信息。

在这一部分中,我们了解了解决勘探-开发权衡的不同方法,并建立了多种优化方法的多样化工具集。虽然我们已经学会了在前一部分中使用 GPyTorch 实现 GPs,但 BoTorch,首要的 BayesOpt 库,是我们在这一部分的重点。我们学习如何声明 BayesOpt 策略,使用它们来促进优化循环,并比较它们在各种任务中的性能。到本部分结束时,您将获得关于实施和运行 BayesOpt 策略的实践知识。

第五章:通过基于改进的政策优化最佳结果

本章包括

  • BayesOpt 循环

  • 在 BayesOpt 政策中开发开发和探索之间的权衡

  • 以改进为标准来寻找新数据点。

  • 使用改进的 BayesOpt 政策

在本章中,我们首先提醒自己 BayesOpt 的迭代性质:我们在已收集的数据上交替训练高斯过程(GP)并使用 BayesOpt 政策找到下一个要标记的数据点。这形成了一个良性循环,在这个循环中,我们的过去数据指导未来的决策。然后我们谈论我们在 BayesOpt 政策中寻找什么:一个决策算法,它决定标记哪个数据点。一个好的 BayesOpt 政策需要平衡足够探索搜索空间并集中在高性能区域。

最后,我们了解了两种政策,这两种政策旨在改进到目前为止 BayesOpt 循环中见过的最佳数据点:改进概率和最常用的 BayesOpt 政策之一,即期望改进。例如,如果我们有一个超参数调整应用程序,我们想要识别在数据集上提供最高验证准确性的神经网络,并且到目前为止我们见过的最高准确性为 90%,那么我们很可能想要改进这个 90%的阈值。本章中学到的策略试图创建这种改进。在我们在表 1.2 中看到的材料发现任务中,我们想要搜索混合温度低(对应高稳定性)的金属合金,而我们找到的最低温度为 187.24,上述两种政策将寻求找到低于这个 187.24 基准的值。

令人惊讶的是,由于我们对目标函数的信念的高斯性质,我们可以期望从最佳观察点的改进程度可以通过封闭形式计算。也就是说,虽然我们不知道未见位置的目标函数的样子,但是在 GP 下,仍然可以轻松地计算基于改进的数量。到本章结束时,我们深入了解了 BayesOpt 政策需要做什么以及如何使用这两个基于改进的政策来完成此操作。我们还学习了如何集成 BoTorch,即我们从本章到本书结束时使用的 Python 中的 BayesOpt 库(botorch.org/docs/introduction),以实现 BayesOpt 政策。

4.1 在 BayesOpt 中导航搜索空间

我们如何确保我们正确利用过去的数据来指导未来的决策?我们在 BayesOpt 政策中寻找什么作为自动决策程序?本节回答了这些问题,并为我们清楚地说明了在使用 GP 时 BayesOpt 的工作原理。

具体来说,在下一小节中,我们将重新检查在第 1.2.2 节中简要介绍的贝叶斯优化循环,以了解如何通过贝叶斯优化策略与高斯过程的训练同时进行决策。然后,我们将讨论贝叶斯优化策略需要解决的主要挑战:在具有高不确定性的区域和在搜索空间中利用良好区域之间的平衡。

以图 4.1 为例,该图显示了在 1 和 2 处训练的两个数据点的高斯过程。在这里,贝叶斯优化策略需要决定我们应该在 -5 和 5 之间的哪个点评估下一个目标函数。探索和利用的权衡是明显的:我们需要决定是否在搜索空间的左右极端(大约在 -5 和 5 附近)检查目标函数,在这些位置,我们的预测存在相当大的不确定性,或者留在平均预测值最高的区域周围。探索和利用的权衡将为我们讨论的不同贝叶斯优化策略奠定基础,这些策略分布在本章和后续章节中。

图 4.1 贝叶斯优化中的探索和利用的权衡。每个策略都需要决定是否查询具有高不确定性的区域(探索),还是查询具有高预测平均值的区域(利用)。

4.1.1 贝叶斯优化循环和策略

首先,让我们回顾一下贝叶斯优化循环的外观以及贝叶斯优化策略在此过程中的作用。在本章中,我们还实现了这个循环的框架,我们将在以后的章节中用来检查贝叶斯优化策略。回顾图 4.2,它是图 1.6 的重复,显示了贝叶斯优化在高层次上的工作方式。

图 4.2 贝叶斯优化循环,结合了模拟高斯过程和决策制定策略。此完整工作流现在可以用于优化黑盒函数。

具体来说,贝叶斯优化是通过一个循环完成的,该循环在以下几个方面交替进行:

  • 对当前训练集进行高斯过程(GP)训练。我们已经在以前的章节中彻底讨论了如何做到这一点。

  • 使用训练有素的高斯过程(GP)对搜索空间中的数据点进行评分,评估它们在帮助我们确定目标函数最优值时的价值(图 4.2 中的步骤 2)。选择最大化此分数的点进行标记,并添加到训练数据中(图 4.2 中的步骤 3)。如何进行此评分由使用的贝叶斯优化策略决定。我们在本章和第 5 和第六章中了解更多有关不同策略的信息。

我们重复此循环,直到达到终止条件,通常是一旦我们已经评估了目标函数达到目标循环迭代次数。这个过程是端到端的,因为我们不仅仅是使用高斯过程进行预测,而是利用这些预测来决定下一步收集哪些数据点,进而推动未来预测的生成。

定义 BayesOpt 循环是模型训练(GP)和数据收集(策略)的良性循环,彼此互相帮助和受益,目标是找到目标函数的最优解。BayesOpt 的良性循环是一个反馈循环,向着具有期望性质的平衡状态迭代;其组件协同作用以实现期望的结果,而不是在导致不良结果的恶性循环中相互恶化。

决定数据点如何根据其在帮助我们实现这一目标的价值方面得分的规则由 BayesOpt 策略决定,因此对于优化性能至关重要。一个好的策略会将高分分配给对优化真正有价值的数据点,从而更快、更高效地指引我们走向目标函数的最优点,而一个设计不良的策略可能会误导我们的实验,并浪费宝贵的资源。

定义 BayesOpt 策略 根据其价值对每个潜在的查询进行评分,从而决定我们应该在下一步查询目标函数的位置(评分最高的地方)。由策略计算的此评分被称为获取分,因为我们将其用作数据获取的方法。

连接到强化学习策略

如果你有强化学习(RL)的经验,你可能会注意到 BayesOpt 策略和 RL 策略之间的联系。在这两种技术中,策略告诉我们在决策问题中应该采取哪些行动。而在 RL 中,策略可能为每个动作分配一个分数,然后我们选择分数最高的动作,或者策略可能只是输出我们应该采取的动作。在 BayesOpt 中,它始终是前者,策略输出一个量化每个可能查询价值的分数,因此我们的工作是确定最大化此分数的查询。

不幸的是,设计一个好的 BayesOpt 策略的问题没有完美的解决方案。也就是说,并没有一种单一的 BayesOpt 策略能够在所有目标函数上始终表现出色。正如我们在本章和后续章节中所见,不同的策略使用不同的焦点和启发式方法。虽然某些启发式方法在某种类型的目标函数上效果良好,但其他启发式方法可能在不同类型的函数上有效。这意味着我们需要接触广泛的 BayesOpt 策略,并了解它们的目的,以便将它们应用于适当的情况——这正是我们将在第 4 至 6 章中所做的。

什么是策略?

每个 BayesOpt 策略都是一个决策规则,根据给定的标准或启发式评分数据点,以确定其在优化中的有用性。不同的标准和启发式导致不同的策略,而没有预先确定的一组 BayesOpt 策略。事实上,BayesOpt 研究人员仍然发布提出新策略的论文。在本书中,我们只讨论实践中最流行和常用的策略。

现在让我们花点时间来实现一个占位的贝叶斯优化循环,从现在开始我们将使用它来检查各种贝叶斯优化策略。这段代码在 CH04/01 - BayesOpt loop.ipynb 中实现。我们需要的第一个组件是一个我们想要使用贝叶斯优化来优化的目标函数。在这里,我们使用熟悉的一维 Forrester 函数作为要最大化的目标函数,它被定义在 -5 到 5 之间。我们还使用 xsys 作为基本事实,在其定义域 [–5, 5] 内计算 Forrester 函数的值:

def forrester_1d(x):                                              ❶
    y = -((x + 1) ** 2) * torch.sin(2 * x + 2) / 5 + 1return y.squeeze(-1)                                          ❶

bound = 5                                                         ❷

xs = torch.linspace(-bound, bound, bound * 100 + 1).unsqueeze(1)  ❷
ys = forrester_1d(xs)                                             ❷

❶ 目标函数的形式,假设未知

❷ 在 -5 到 5 之间的网格上计算的测试数据

我们需要做的另一件事是修改 GP 模型的实现方式,以便它们可以与 BoTorch 中的贝叶斯优化策略一起使用。实现 GP 构成了我们的贝叶斯优化循环的第一步。

由于 BoTorch 就建立在 GPyTorch 之上,因此只需要进行最小的修改。具体来说,我们使用以下 GP 实现,其中除了我们通常的 gpytorch.models.ExactGP 类之外,我们还继承了 botorch.models.gpytorch.GPyTorchModel 类。此外,我们声明了类特定属性 num_outputs 并将其设置为 1。这些是我们需要进行的最小修改,以便将我们的 GPyTorch 模型与 BoTorch 一起使用,后者实现了本章后面我们使用的贝叶斯优化策略:

class GPModel(gpytorch.models.ExactGP,
  botorch.models.gpytorch.GPyTorchModel):  ❶
    num_outputs = 1def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel()
        )
    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

❶ 用于 BoTorch 集成的修改

除此之外,我们 GP 实现中的其他一切都保持不变。现在我们编写一个帮助函数,用于在我们的训练数据上训练 GP:

def fit_gp_model(train_x, train_y, num_train_iters=500):
    noise = 1e-4                                     ❶

    likelihood = gpytorch.likelihoods
    ➥.GaussianLikelihood()                          ❷
    model = GPModel(train_x, train_y, likelihood)    ❷
    model.likelihood.noise = noise                   ❷

    optimizer = torch.optim.Adam(model.parameters(),
    ➥lr=0.01)                                       ❸
    mll = gpytorch.mlls.ExactMarginalLogLikelihood
    ➥(likelihood, model)                            ❸

    model.train()                                    ❸
    likelihood.train()                               ❸

    for i in tqdm(range(num_train_iters)):           ❸
        optimizer.zero_grad()                        ❸

        output = model(train_x)                      ❸
        loss = -mll(output, train_y)                 ❸

        loss.backward()                              ❸
        optimizer.step()                             ❸

    model.eval()                                     ❸
    likelihood.eval()                                ❸

    return model, likelihood

❶ 使用梯度下降训练 GP

❷ 声明 GP

❸ 使用梯度下降训练 GP

注意 我们在前几章中使用了所有先前的代码。如果您对某段代码有困难理解,请参考第 3.3.2 节以获取更多细节。

这涵盖了图 4.2 的第 1 步。现在,我们跳过第 2 步,即实现贝叶斯优化策略,并将其留待下一节和未来章节。要实现的下一个组件是可视化到目前为止收集的数据,当前 GP 信念以及一个贝叶斯优化策略如何评分其余数据点。该可视化的目标显示在图 4.3 中,我们在第一章中见过。具体来说,图的顶部面板显示了 GP 模型对真实目标函数的预测,而底部面板显示了由贝叶斯优化策略计算的收购分数。

图 4.3 贝叶斯优化进展的典型可视化。顶部面板显示了 GP 预测和真实的目标函数,而底部面板显示了一个名为期望改进(Expected Improvement)的贝叶斯优化策略所做的收购分数,我们将在第 4.3 节中学习到。

我们已经熟悉如何生成顶部面板,生成底部面板同样简单。这将使用类似于我们在第 3.3 节中使用的辅助函数来完成。该函数接受一个 GP 模型及其似然函数以及两个可选输入:

  1. policy 指的是 BayesOpt 策略对象,可以像调用任何 PyTorch 模块一样调用。在这里,我们将它调用在代表我们搜索空间的网格 xs 上,以获取整个空间的获取分数。我们将在下一节中讨论如何使用 BoTorch 实现这些策略对象,但我们现在不需要更多了解这些对象。

  2. next_x 是使获取分数最大化的数据点的位置,将其添加到正在运行的训练数据中:

def visualize_gp_belief_and_policy(
    model, likelihood, policy=None, next_x=None
):
    with torch.no_grad():
        predictive_distribution = likelihood(model(xs))  ❶
        predictive_mean = predictive_distribution.mean   ❶
        predictive_upper, predictive_lower =             ❶
          ➥predictive_distribution.confidence_region()  ❶

        if policy is not None:                           ❷
            acquisition_score = policy(xs.unsqueeze(1))  ❷

    ...                                                  ❸

❶ GP 预测

❷ 获取分数

❸ 省略

在这里,我们从 GP 和测试数据 xs 中生成预测。请注意,如果未传入 policy,我们不会计算获取分数,在这种情况下,我们也以我们已经熟悉的方式可视化 GP 预测-散点表示训练数据,平均预测的实线,95% CI 的阴影区域:

    if policy is None:
        plt.figure(figsize=(8, 3))

        plt.plot(xs, ys, label="objective", c="r")         ❶
        plt.scatter(train_x, train_y, marker="x", c="k",
            label="observations")                          ❷

        plt.plot(xs, predictive_mean, label="mean")        ❸
        plt.fill_between(                                  ❸
            xs.flatten(),                                  ❸
            predictive_upper,                              ❸
            predictive_lower,                              ❸
            alpha=0.3,                                     ❸
            label="95% CI",                                ❸
        )                                                  ❸

        plt.legend()
        plt.show()

❶ 真实值

❷ 训练数据

❸ 平均预测和 95% CI

请参考第 2.4.4 节以了解这个可视化的基础知识。

另一方面,如果传入了策略对象,我们将创建另一个子图以显示搜索空间中的获取分数:

    else:
        fig, ax = plt.subplots(
            2,
            1,
            figsize=(8, 6),
            sharex=True,
            gridspec_kw={"height_ratios": [2, 1]}
        )

        ...                                                    ❶

        if next_x is not None:                                 ❷
            ax[0].axvline(next_x, linestyle="dotted", c="k")   ❷

        ax[1].plot(xs, acquisition_score, c="g")               ❸
        ax[1].fill_between(                                    ❸
          xs.flatten(),                                        ❸
          acquisition_score,                                   ❸
          0,                                                   ❸
          color="g",                                           ❸
          alpha=0.5                                            ❸
        )                                                      ❸

        if next_x is not None:                                 ❷
            ax[1].axvline(next_x, linestyle="dotted", c="k")   ❷

        ax[1].set_ylabel("acquisition score")

        plt.show()

❶ GP 预测(与以前相同)

❷ 最大化获取分数的点,使用虚线垂直线进行可视化

❸ 获取分数

当传入 policynext_x 时,此函数将创建一个显示根据 BayesOpt 策略的获取分数的较低面板。最后,我们需要实现图 4.2 中 BayesOpt 循环的第 3 步,即(1)找到具有最高获取分数的点,并(2)将其添加到训练数据并更新 GP。对于识别给出最高获取分数的点的第一个任务,虽然在我们的 Forrester 示例中可能在一维搜索空间上进行扫描,但随着目标函数维数的增加,穷举搜索变得越来越昂贵。

请注意,我们可以使用 BoTorch 的辅助函数 botorch.optim.optimize .optimize_acqf(),该函数找到最大化任何 BayesOpt 策略得分的点。辅助函数使用 L-BFGS,一种准牛顿优化方法,通常比梯度下降方法更有效。

我们这样做:

  • policy 是 BayesOpt 策略对象,我们很快就会了解更多。

  • bounds 存储了我们搜索空间的边界,在本例中为-5 和 5。

  • q = 1 指定我们希望辅助函数返回的点数,这是一个。 (在第七章中,我们学习了允许同时对目标函数进行多次查询的设置。)

  • num_restartsraw_samples 分别表示在搜索给出最高获取分数的最佳候选项时 L-BFGS 使用的重复次数和初始数据点数。一般来说,我建议分别使用这些参数的维数的 20 倍和 50 倍。

  • 返回的值,next_xacq_val,分别是给出最高获取分数的点的位置和相应的最大化获取分数:

next_x, acq_val = botorch.optim.optimize_acqf(
    policy,
    bounds=torch.tensor([[-bound * 1.0], [bound * 1.0]]),
    q=1,
    num_restarts=20,
    raw_samples=50,
)

设置重新启动和原始样本的数量

num_restartsraw_samples 的值越高,当搜索最大化获取分数的最佳候选项时,L-BFGS 的穷举程度就越高。这也意味着 L-BFGS 算法运行的时间将更长。如果发现 L-BFGS 在最大化获取分数时失败,或者算法运行时间太长,可以增加这两个数字;反之,可以减少它们。

作为最后一步,我们在贝叶斯优化循环中汇总了我们到目前为止实现的内容。在该循环的每次迭代中,我们执行以下操作:

  1. 我们首先打印出迄今为止我们看到的最佳值(train_y.max()),这显示了优化的进展情况。

  2. 然后,我们重新在当前训练数据上训练 GP,并重新声明贝叶斯优化策略。

  3. 使用 BoTorch 中的辅助函数 botorch.optim.optimize_acqf(),我们确定在搜索空间中最大化获取分数的点。

  4. 我们调用辅助函数 visualize_gp_belief_and_policy(),该函数可视化我们当前的 GP 信念和优化进展。

  5. 最后,我们在确定的点(next_x)查询函数值并更新我们观察到的数据。

整个过程总结在图 4.4 中,该图显示了贝叶斯优化循环中的步骤及实现这些步骤的相应代码。每个步骤都由我们的辅助函数或 BoTorch 的模块化代码实现,使得整个过程易于跟踪。

图 4.4 贝叶斯优化循环中的步骤及相应的代码。每个步骤的代码是模块化的,这使得整个循环易于跟踪。

实际的代码实现如下:

num_queries = 10for i in range(num_queries):
    print("iteration", i)
    print("incumbent", train_x[train_y.argmax()], train_y.max())

    model, likelihood = fit_gp_model(train_x, train_y)   ❷

    policy = ...                                         ❸

    next_x, acq_val = botorch.optim.optimize_acqf(       ❹
        policy,                                          ❹
        bounds=torch.tensor([[-bound * 1.0],             ❹
        ➥[bound * 1.0]]),                               ❹
        q=1,                                             ❹
        num_restarts=20,                                 ❹
        raw_samples=50,                                  ❹
    )                                                    ❹

    visualize_gp_belief_and_policy(model, likelihood, policy,
        next_x=next_x)                                   ❺

    next_y = forrester_1d(next_x)                        ❻

    train_x = torch.cat([train_x, next_x])               ❻
    train_y = torch.cat([train_y, next_y])               ❻

❶ 可以进行的目标函数评估次数

❷ 更新当前数据上的模型

❸ 初始化贝叶斯优化策略,稍后讨论

❹ 找到给出最高获取分数的点

❺ 可视化当前 GP 模型和获取分数

❻ 在确定的点观察并更新训练数据

有了这一点,我们已经实现了一个 BayesOpt 循环的框架。现在唯一要做的就是用我们想要使用的实际 BayesOpt 策略来填充policy的初始化,然后笔记本就能在 Forrester 函数上运行 BayesOpt 了。请注意,虽然调用visualize_gp_belief_and_policy()不是必需的(也就是说,之前的 BayesOpt 循环仍然能够运行而不需要那一行代码),但是这个函数对我们观察 BayesOpt 策略的行为和特性以及诊断任何潜在问题是有用的,正如我们后面在本章中讨论的那样。

BayesOpt 策略最重要的特征之一是探索和利用之间的平衡,这是许多人工智能和机器学习问题中的一个经典权衡。在这里,发现我们目前不知道的高性能区域(探索)的可能性与集中在已知的良好区域(利用)的机会之间进行权衡。我们将在下一小节中更详细地讨论这种权衡。

4.1.2 平衡探索和利用

在这一小节中,我们讨论了决策过程中固有的一个问题,包括 BayesOpt 在内:在充分探索整个搜索空间和及时利用产生良好结果的区域之间的平衡。这一讨论将帮助我们形成对什么是一个好的 BayesOpt 策略的理解,并让我们意识到我们所学到的每个策略如何解决这种权衡。

为了说明探索和利用的权衡,想象一下你正在一家你之前只去过几次的餐厅用餐(见图 4.5)。你知道这家餐厅的汉堡很棒,但你不确定他们的鱼和牛排是否好吃。在这里,你面临着一个探索与利用的问题,你需要在尝试点可能是一道很好的菜肴(探索)和点你经常吃但可靠的餐点(利用)之间做出选择。

图 4.5 在餐厅点菜具有固有的探索(尝试新事物)与利用(点常吃的)的权衡。

过度的探索可能会导致你点到你不喜欢的东西,而持续的利用则可能会导致你错过你真正喜欢的一道菜。因此,两者之间的合理平衡至关重要。

这个无处不在的问题不仅存在于点餐中,也存在于人工智能的常见问题中,比如强化学习、产品推荐和科学发现。在 BayesOpt 中,我们面临着同样的权衡:我们需要充分探索搜索空间,以便不错过一个好的区域,但我们也应该专注于具有高客观价值的区域,以确保我们在优化方面取得进展。

注意:“具有高目标值的区域”指的是由输入 x 产生高输出 f(x) 值的区域,这是我们优化(特别是最大化)任务的目标。

让我们回到我们的代码示例,并假设当我们刚开始时,我们的训练数据集包含 Forrester 目标函数的两个观察值,分别为 x = 1 和 x = 2:

train_x = torch.tensor([
    [1.],
    [2.]
])
train_y = forrester_1d(train_x)

model, likelihood = fit_gp_model(train_x, train_y)

print(torch.hstack([train_x, train_y.unsqueeze(1)]))

这给出了输出

tensor([[1.0000, 1.6054],
        [2.0000, 1.5029]])

这表明点 1 处的评估大约为 1.6,而点 2 处的评估为 1.5。通过可视化训练好的 GP 所做出的预测,我们获得了图 4.6 中熟悉的样式。该图显示了我们面临的探索与利用之间的权衡:我们应该在不确定性较高的地方评估目标函数,还是应该留在平均预测较高的区域?

图 4.6:由 GP 训练的两个数据点的预测来自 Forrester 函数

每个 BayesOpt 策略都有一种不同的方式来处理这种权衡,因此,在如何最好地探索搜索空间方面提供了不同的建议。在图 4.6 中,一些策略可能会导致我们进一步探索未知区域,而另一些策略可能会建议我们聚焦于已知的高价值区域。同样,通常没有一种一刀切的方法(即,没有一种策略始终表现良好)。

4.2 寻找 BayesOpt 中的改进

我们几乎已经准备好在给定目标函数上运行 BayesOpt 了。现在我们需要一个具有关于我们搜索目标中每个潜在标记数据点的价值的评分规则的策略。同样,我们将看到的每个策略都提供了一个不同的评分规则,这些规则是受到不同的优化启发式的驱使的。

在这一部分,我们学习了一种在优化目标时具有直观意义的启发式方法。当我们寻求优化时,这种方法的形式是寻求从迄今为止所见过的最佳点开始改进。在即将到来的小节中,我们将了解到 GPs 有助于促进这一改进度量的计算。然后,我们将介绍不同定义改进的方式如何导致两种最常见的 BayesOpt 策略:提高概率和期望改进。

4.2.1 用 GP 测量改进

在本小节中,我们将讨论改进 BayesOpt 的定义、它如何构成一个良好的效用度量,并且使用正态分布处理与改进相关的量是直接的。因为我们在 BayesOpt 中的最终目标是确定目标函数的全局最优值——给出最高目标值的点——所以我们在评估目标函数时所观察到的值越高,我们的效用就应该越高。假设当我们在某个点x[1]处评估目标函数时,我们观察到值为 2。在另一种情况下,当我们评估另一个点x[2]时,我们观察到值为 10。直观地说,我们应该更重视第二个点,因为它给了我们一个更高的函数值。

但是,如果在观察到x[1]和x[2]之前,我们已经看到一个x[0]的点,它的值为 20 呢?在这种情况下,自然会认为即使x[2]比x[1]好,但两个点都不会带来任何额外的效用,因为我们在x[0]上已经有了一个更好的观察结果。另一方面,如果我们有一个x[3]的点,它的值为 21,那么我们会更高兴,因为我们已经找到了比x[0]更好的值。这在图 4.7 中有所说明。

图 4.7 寻求从观察到的最佳点改进。虽然点x[2]优于点x[1],但两者都“不好”意味着它们没有改进点x[0]。

这些比较指出了在 BayesOpt 中,我们关心的不仅仅是我们观察到的观察值的原始值,还有新发现的观察是否比我们的观察更好。在这种情况下,由于x[0]在函数值方面设定了一个非常高的标准,无论x[1]还是x[2]都不构成改进——至少不是我们在优化中关心的那种改进。换句话说,在优化中一个合理的目标是寻求从迄今为止我们见过的最佳点改进,因为只要我们从观察到的最佳点改进,我们就在取得进步。

定义 观察到的最佳点,或者产生了迄今为止我们找到的最高函数值的点,通常被称为现任。这个术语表示,这个点目前在我们搜索期间查询的所有点中持有最高值。

假设我们对目标函数有一个 GP 的信念,那么检查我们可以期望从观察到的最佳点改进多少可能会很容易。让我们从我们正在运行的 Forrester 函数的示例开始,我们的当前 GP 如图 4.6 所示。

在我们的训练集中的两个数据点中,(x = 1,y = 1.6) 是更好的一个,因为它有一个更高的函数值。这是我们当前的现任。同样,我们正在专注于从这个 1.6 的阈值改进;也就是说,我们希望找到的数据点会产生高于 1.6 的函数值。

从视觉上来看,我们可以将这个基于改进的方法想象为将高斯过程水平切断在现有解处,如图 4.8 所示。较暗颜色突出显示的部分对应于“改进现有解”(在 1.6 处)。任何低于这条线的点都不能构成改进,因此不会给我们带来额外的优化效益。虽然我们不知道某个点-比如x=0-是否会产生更高的函数值,但我们仍然可以尝试通过从高斯过程模型中得知的信息来推理它的可能性。

图 4.8 从高斯过程的角度看改进自现有解。高斯过程预测对应的改进自现有解的部分在较深的颜色中突出显示。

推理 x = 0 的概率将产生更高的函数值很容易做到,因为通过查询点 x = 0,我们观察到的改进自现有解正好对应于一个被部分截断的正态分布,如图 4.9 所示。

图 4.9 的左面板包含与图 4.8 相同的高斯过程,该过程在现有解处被切断,并额外显示了在x=0 处的正态分布预测的 CI。在此点 0 处垂直切割高斯过程,我们获得了图 4.9 的右面板,两个面板中的 CI 相同。我们看到,在右面板中正态分布的只有突出显示的部分代表我们可以从现有解中观察到的改进,这是我们关心的。这突出显示的部分是正态分布的一部分,正如我们在接下来的小节中所介绍的那样,这导致了许多数学的便利。

图 4.9 在 0 处对现有解的改进,以较深的颜色突出显示。左面板显示整个高斯过程,而右面板仅显示与 0 处的预测相对应的正态分布(误差栏在两个面板上相同)。在这里,对现有解的改进遵循一个被截断的正态分布。

这些便利不仅适用于 x = 0 。由于高斯过程在任何给定点处的预测是正态分布,因此在任何点处的改进也遵循被截断的正态分布,如图 4.10 所示。

图 4.10 在–2 和 3 处对现有解的改进,以较深的颜色突出显示。左面板显示整个高斯过程,中心面板显示-2 处的预测,右面板显示 3 处的预测。突出显示的部分显示可能的改进,这取决于给定点上的正态分布。

我们看到与图 4.9 相比,在 0 处突出显示了大部分正态分布为可能的改进,以下是正确的:

  • -2 处的预测(中心面板)更糟,因为它只有一半被突出显示为可能改进。这是因为 -2 处的平均预测值大致等于现有值,因此,在 -2 处的函数值改进为 1.6 的可能性是 50-50。

  • 作为另一个例子,右侧面板显示了在 3 处的预测,根据我们关于客观情况的高斯过程的信念,几乎不可能从现有情况改进,因为几乎整个正态分布在 3 处都低于现有门槛。

这表明不同的点将导致不同的可能改进,这取决于高斯过程的预测。

4.2.2 计算改进的概率

在 BayesOpt 中,我们的目标特别明确,即从当前的现有情况中提高。现在,我们终于准备开始讨论旨在实现这一目标的 BayesOpt 策略。在本小节中,我们了解改进概率(PoI),这是一种衡量候选点可能改进的策略。

衡量候选点从现有情况改进的可能性的想法对应于图 4.7 中一个点是“好”的概率的概念。这在第 4.2.1 节中与高斯过程中也有所提及,我们说:

  1. 0 处的点(图 4.9)的正态分布的大部分被突出显示为可能改进。换句话说,它有很高的改进可能性。

  2. -2 处的点(图 4.10 的中心面板)有 0.5 的改进概率,因为其正态分布的一半超过了现有状态。

  3. 3 处的点(图 4.10 的右侧面板)的正态分布大部分在门槛以下,因此它几乎不可能改进。

通过注意到从现有情况改进的概率等于图 4.9 和 4.10 中突出显示部分的面积,我们可以使这一计算更具体。

定义 任何正态曲线下的整个面积都为 1,因此图 4.9 和 4.10 中突出显示部分的面积恰好衡量了该正态随机变量(给定点的函数值)超过现有情况的可能性。

与突出显示的区域的面积对应的量与累积密度函数(CDF)有关,它被定义为一个随机变量取得小于或等于目标值的概率。换句话说,CDF 衡量了突出显示的区域的面积,即突出显示的区域的面积为 1 减去突出显示的区域的面积。

多亏了数学上的便利性,正态分布和高斯过程,我们可以轻松地使用累积分布函数计算此突出区域的面积,这需要相关正态分布的均值和标准差。在计算上,我们可以利用 PyTorch 的torch.distributions.Normal类,该类实现了正态分布并提供了有用的cdf()方法。具体来说,假设我们有兴趣计算 0 点的点能够改进现有情况的可能性。我们将按照图 4.11 中描述的过程进行操作:

  1. 首先,我们使用高斯过程来计算 0 点的均值和标准差预测。

  2. 接着,我们计算以先前的均值和标准差定义的正态曲线下的面积,以现有值为截止点。我们使用累积分布函数进行此计算。

  3. 最后,我们从 1 中减去 CDF 值,以获得候选点的 PoI。

图 4.11

图 4.11 PoI 分数计算的流程图。通过按照此过程,我们可以计算任何候选点能够从现有情况中改进的可能性。

注释 技术上,累积分布函数计算正态分布左侧的部分的面积,因此我们需要从 1 中减去 CDF 的输出,以获取阈值右侧的部分的面积,这对应于可能的改进。

我们首先生成该点处高斯过程的预测:

with torch.no_grad():
    predictive_distribution = likelihood(model(torch.tensor([[0.]])))
    predictive_mean = predictive_distribution.mean     ❶
    predictive_sd = predictive_distribution.stddev     ❷

❶ 0 点处的预测均值

❷ 0 点处的预测标准差

首先,我们使用相应的均值和标准差初始化一个一维正态分布:

normal = torch.distributions.Normal(predictive_mean, predictive_sd)

此正态分布是图 4.9 右侧面板中可视化的正态分布。最后,为了计算突出显示部分的面积,我们调用带有现有值的cdf()方法(即train_y.max(),我们的训练数据的最大值),并从 1 中减去结果:

>>> 1 - normal.cdf(train_y.max())

tensor([0.8305])

在这里,我们的代码显示,在 0 点处,我们有超过 80%的机会从现有情况中改进,这与图 4.9 中突出显示的正态分布的大部分一致。使用相同的计算,我们可以发现-2 点的 PoI 为 0.4948,3 点的 PoI 为 0.0036。通过观察图 4.10,我们可以看到这些数字是合理的。除了这三个点(0、-2 和 3)之外,我们还可以使用相同的公式计算给定点能够改进的可能性,即给定点的 PoI,跨越我们的搜索空间。

定义 图 4.11 中的过程给出了 PoI 策略的评分规则,其中搜索空间中每个点的得分等于该点能够改进现有情况的可能性。同样,这个分数也称为收集分数,因为我们将其用作数据收集的方法。

我们可以看到,此 PoI 策略使用了cdf()方法,但使用BoTorch更干净,这是一个实现 BayesOpt 策略的 Python 库。BoTorch 建立在 PyTorch 和 GPyTorch 之上,并与其无缝合作。正如我们之前所看到的,我们只需更改 GP 类中的两行代码,即可使模型与 BoTorch 兼容。此外,BoTorch 将其策略实现为模块,使我们能够以模块化的方式在 BayesOpt 循环中交换不同的策略。

BoTorch 策略的模块化

通过模块化,我们指的是我们可以在 BayesOpt 循环中用另一个策略替换当前正在使用的策略,只需更改策略的初始化。BayesOpt 循环的其余部分(训练 GP 模型、可视化优化进度和更新训练数据)不必更改。我们在 3.3 和 3.4.2 节中也观察到了与 GPyTorch 的均值函数和核函数类似的模块化。

要使用 BoTorch 实现 PoI 策略,我们执行以下操作:

policy = botorch.acquisition.analytic.ProbabilityOfImprovement( ❶
    model, best_f=train_y.max()                                 ❶
)                                                               ❶

with torch.no_grad():                                           ❷
    scores = policy(xs.unsqueeze(1))                            ❷

❶ 声明 PoI 策略

❷ 计算分数

BoTorch 类ProbabilityOfImprovement将 PoI 实现为 PyTorch 模块,将 GP 作为第一个参数,现有结果值作为第二个参数。变量scores现在存储xs中点的 PoI 分数,xs是-5 到 5 之间的密集网格。

再次,每个点的获取分数等于该点从现有结果改善的概率,根据我们的高斯过程信念。图 4.12 显示了我们的搜索空间中的这个分数,以及我们的高斯过程。

图 4.12 高斯过程预测(顶部)和 PoI(底部),虚线表示最大化 PoI 分数的点。这一点是我们在优化的下一次迭代中查询目标函数的地方。

我们观察到一些有趣的 PoI 分数行为:

  • 在现有结果的左侧区域(从 0 到 1)的 PoI 分数相对较高。这对应于此区域的高均值预测。

  • 绘图左侧区域的得分稍低。这是因为均值预测不那么高,但在这个区域存在足够的不确定性,仍然有相当大的改善概率。

  • 2 周围的区域的 PoI 接近 0。正如我们所见,此处的点的预测正态分布大多位于现有结果阈值下方。

现在,我们要做的就是确定在-5 到 5 之间具有最高 PoI 分数的点,即最大化从现有结果改善的概率的点。如前所述,我们利用 BoTorch 的辅助函数botorch.optim.optimize.optimize_acqf(),该函数找到最大化任何 BayesOpt 策略得分的点。我们使用以下代码执行此操作,该代码是实现 BayesOpt 循环的代码的一部分:

next_x, acq_val = botorch.optim.optimize.optimize_acqf(
    policy,
    bounds=torch.tensor([[-bound], [bound]], dtype=torch.float),
    q=1,
    num_restarts=10,
    raw_samples=20,
)

返回的值是 L-BFGS 找到的给出最高获取分数的点的位置以及相应的最大化获取分数。经过检查,我们有以下结果:

>>> next_x, acq_val

(tensor([[0.5985]]), tensor(0.9129))

此输出表明,最大化 PoI 分数的候选点在 0.91 PoI 处大约在 0.6 附近,这对应于图 4.12 中的点线。这个点就是我们将查询目标函数(即评估函数)以收集 BayesOpt 中的下一个数据点的地方。

具有最高预测均值的候选点

有趣的是,我们选择查询的点(大约在 0.6 附近)并不是具有最高预测均值的点(大约在 -0.5 左右)。后者的 PoI 稍低一些,因为此点的不确定性很高,因此事实上,与现有点相比,它不太可能得到改善,尽管其预测均值很高。

这就是在 BayesOpt 的单次迭代中决定使用 PoI 策略查询哪个点的全部内容。但是请记住图 4.2 中的 BayesOpt 循环,在该循环中我们在寻找下一个要查询的数据点(步骤 2)和使用新数据更新我们的 GP(步骤 1 和 3)之间交替。我们将在第 4.2.3 节中执行此操作。

运行 PoI 策略

在本小节中,我们最终运行 PoI 策略并分析其行为。我们再次重复整个过程——训练模型、声明 PoI 策略,并使用 optimize_acqf() 多次找到最佳点,直到达到终止条件为止。正如我们所见,这个循环在 CH04/01 - BayesOpt loop.ipynb 笔记本中实现。现在,我们需要在适当的 for 循环中初始化 PoI 策略。

此代码生成了一系列由辅助函数 visualize_gp_belief_and_policy() 生成的图,每个图显示了我们 BayesOpt 循环在我们进行的 10 次查询中的当前状态。这些图看起来类似于图 4.12,但增加了我们参考的目标函数:

num_queries = 10

for i in range(num_queries):
    print("iteration", i)
    print("incumbent", train_x[train_y.argmax()], train_y.max())

    model, likelihood = fit_gp_model(train_x, train_y)

    policy = botorch.acquisition.analytic.ProbabilityOfImprovement( ❶
        model, best_f=train_y.max()                                 ❶
    )                                                               ❶

    next_x, acq_val = botorch.optim.optimize_acqf(
    ...                                                             ❷

❶ 我们的 PoI 策略

❷ 已省略

BayesOpt 中的函数评估次数

我们在 BayesOpt 中使用的查询次数完全取决于我们能够承担的函数评估次数。第 1.1 节定义了昂贵的黑盒优化问题,这假设我们可以进行的查询次数相对较低,因为函数评估的成本很高。

还有其他标准可以用来确定何时终止 BayesOpt 循环。例如,当我们达到目标值或最近的 5 或 10 次查询中没有显著改进时,我们可以停止。在整本书中,我们坚持假设我们有一定数量的函数评估是可以进行的。

我们使用 10 次查询作为默认值来运行一维 Forrester 函数的 BayesOpt 策略,并检查策略的行为。本章的练习 2 处理二维函数并使用 20 次查询。

图 4.13 显示了第一、第五和最后一次迭代的绘图。我们发现 PoI 策略始终保持在 0 和 2 之间的区域内,在第 10 次和最后一次迭代中,我们已经收敛到局部最优点。这意味着我们未能充分探索搜索空间,因此错过了右侧区域中目标函数的全局最优点,即约为 4。

图 4.13 是 PoI 策略的进展情况。由于该策略旨在追求任何数量的改进,进展停滞在接近 2 的局部最优点,我们未能探索其他搜索空间区域。

当使用优化辅助函数 optimize_acqf()时,BoTorch 会发出警告。

当在上一页中运行 BayesOpt 与 PoI 的代码时,您可能会收到 BoTorch 发出的以下警告:

RuntimeWarning: Optimization failed in
`gen_candidates_scipy` with the following warning(s):
[OptimizationWarning('Optimization failed within `scipy.optimize.minimize`
with status 2 and message ABNORMAL_TERMINATION_IN_LNSRCH.')]
Trying again with a new set of initial conditions.
  warnings.warn(first_warn_msg, RuntimeWarning)

当帮助函数 optimize_acqf()(具体来说是线搜索子程序)未能成功优化收集得分(在这种情况下为 PoI 得分)时,会显示此警告。当收集得分函数高度不平滑时(例如,在图 4.13 的最后一个面板中,x=1.5 周围出现了一个尖峰),数值优化不稳定,这种故障经常发生。

不用详细了解优化例行程序,我们可以在使用 optimize_acqf()时增加重启次数(num_restarts 参数)和原始样本数(raw_samples 参数),这样可以增加发现拥有最高收集得分的数据点的机会。

为了便于讲解,从现在开始,在我们的代码中运行帮助函数 optimize_acqf()时,我们关闭此警告,使用警告模块中的上下文管理器:

with warnings.catch_warnings():
    warnings.filterwarnings('ignore', category=RuntimeWarning)
    next_x, acq_val = botorch.optim.optimize_acqf(...)

注意 尽管图 4.13 中的 PoI 表现可能令人失望(毕竟,我们已经花了很多时间构建这个看起来过度利用的 PoI 策略),但分析发生的情况将为我们提供改进性能的见解。

我们注意到,虽然 PoI 停留在局部最优点,但它正在实现其所需的功能。具体来说,由于 PoI 旨在改善当前的 incumbent,因此策略发现缓慢向右移动将以较高的概率实现这一目标。虽然 PoI 通过缓慢移动不断发现更多的改进,但我们将这种行为视为过度利用,因为该策略没有充分探索其他区域。

重要 换句话说,即使 PoI 正在做的事情符合我们最初想要实现的目标——即从在职者中进行改进,但结果的行为并不是我们想要的。这意味着,单纯追求从在职者处的任何改进都不是我们所关心的。

修复这种过度利用行为有两种方法。第一种是限制我们所说的 改进 的含义。我们对 PoI 的实验表明,在每次迭代中,策略只会通过缓慢朝着 GP 认为函数向上移动的方向移动,从现有状况中找到微小的改进。

如果我们重新定义所谓的 改进,即只有当改进比当前现有状况值至少大 ε 时才有效,并相应修改 PoI,那么策略将更有可能更有效地探索搜索空间。这是因为 GP 将知道停留在局部最优点不会导致与现有状况的显著改进。图 4.14 说明了这个想法。

图 4.14

图 4.14 通过要求改进至少为 ε = 0(左)和 ε = 2(右)来定义对 改进 更严格的定义。要求越大,PoI 策略就越具有探索性。更多细节请参见练习 1。

我们不会在这里深入讨论,但练习 1 探讨了这种方法。有趣的是,我们会观察到要求 PoI 观察到的改进越多,策略就越具有探索性。

4.3 优化改进的预期值

如前一节所示,天真地试图从现有状况中寻求改进会导致来自 PoI 的过度利用。这是因为简单地沿着适当方向微移离开现有状况就能获得高的 PoI。因此,优化这个 PoI 并不 是我们想要做的。在本节中,我们将学习进一步考虑我们可能观察到的改进的 幅度。换句话说,我们也关心我们能从现有状况中获得多少改进。这将引导我们进入最流行的贝叶斯优化策略之一:期望改进(EI)。

寻求考虑改进幅度的动机是明确的。考虑图 4.15 中的例子。

图 4.15

图 4.15 PoI(左)和 EI(右)之间的差异。前者只关心我们是否从现有状况改进,而后者考虑了改进的幅度。

左侧面板显示了 PoI 策略的计算,它只考虑每个候选数据点是否从现有状况改进。因此,稍微改进和显着从现有状况改进的点被平等地对待。

另一方面,右侧面板展示了如果我们还考虑可能改进的幅度会发生什么。在这里,虽然点 x[1] 和 x[2] 仍然被视为不理想(因为它们并未从现有状况 x[0] 改进),但 x[4] 被认为比 x[3] 更好,因为前者提供了更大的改进。同样,x[5] 被认为是五个候选点中最好的。

当然,这并不意味着我们现在可以简单地设计选择从现有指标获得最大改进的策略。我们仍然需要知道我们将观察到多少(如果有的话)改进,这只有当我们实际查询目标函数时才能发现。

注意 尽管我们不知道我们将观察到的改进的确切值,但我们可以以概率的方式推断每个候选点的改进量的大小。

回想一下,在图 4.9 和 4.10 中,我们有一个表示在给定点观察到的改进的截断正态分布。通过计算突出显示区域的面积,我们获得了一个点将从现有指标改进的概率,从而得到了 PoI 策略。然而,我们也可以执行其他计算。

定义 除了 PoI 之外,我们还可以计算与突出显示区域对应的随机变量的期望值。我们处理截断正态分布的事实使我们能够以封闭形式计算这个期望值。使用这个度量来评分数据点的 BayesOpt 策略被称为期望改进(EI)。

虽然 EI 得分的闭式公式不像 PoI 那样简单,但 EI 的评分公式仍然很容易计算。直观地说,使用改进的期望值可能比 PoI 更好地平衡探索和利用。毕竟,周围的点可能以很高的概率改进,但它们的改进可能很小(这是我们在实验中经验性地观察到的)。

一个远离我们不太了解的点可能会给出较低的改进概率,但因为有机会这个点会有较大的改进,EI 可能会给它分配较高的分数。换句话说,虽然 PoI 可能被认为是一种风险规避的 BayesOpt 策略,它关心的是从现有指标改进,无论改进有多小,但 EI 在风险和回报之间平衡,找到最好平衡权衡的点。图 4.16 对比了使用相同数据集和训练的 GP 的 PoI 和 EI。

图片

图 4.16 显示了 PoI(左)和 EI(右)之间的区别。EI 更好地平衡了探索和利用。

我们看到 PoI 选择的候选数据点(大约在 0.6 附近)与 EI 选择的数据点(大约在 –0.7 附近)不同。前者接近当前的指标,因此查询它有助于我们改进的可能性较大。然而,EI 看到远离现有指标的其他区域存在更大的不确定性,这可能导致更大的改进。由于这种推理,EI 更倾向于提供更好的探索和利用平衡的数据点。

EI 的另一个好的特性是它如何给具有相同预测均值或标准差的数据点分配获取分数。具体而言,它通过以下方式来实现:

  • 如果两个数据点具有相同的预测均值但不同的预测标准差,那么具有更高不确定性的数据点将获得更高的分数。因此,该策略奖励探索。

  • 如果两个数据点具有相同的预测标准差但不同的预测均值,那么具有更高均值的数据点将获得更高的分数。因此,该策略奖励开发。

这是 BayesOpt 策略的一个期望,因为它表达了我们在一切相等时对探索的偏好(即当预测均值相等时)但也表达了我们在一切相等时对开发的偏好(即不确定性相等时)。我们在第五章中再次可以看到这一特性,另一种 BayesOpt 策略称为置信上界

从计算上讲,我们可以使用几乎与 PoI 代码相同的代码将 EI 初始化为 BayesOpt 策略对象:

policy = botorch.acquisition.analytic.ExpectedImprovement(
    model, best_f=train_y.max()
)

现在,让我们使用 EI 策略重新运行整个 BayesOpt 循环,确保我们从相同的初始数据集开始。这产生了图 4.17,与 PoI 的图 4.13 进行比较。

图 4.17 显示了 EI 策略的进展。与 PoI 相比,该策略在探索和开发之间取得了更好的平衡,并在最后找到了全局最优解。

在这里,我们可以看到,尽管 EI 最初仍然集中在 2 附近的局部最优区域,但策略很快探索了搜索空间中的其他区域,以寻找对现有情况的更大改进。在第五次迭代中,我们可以看到我们现在正在检查左侧的区域。最后,在使用了所有的 10 个查询之后,EI 成功地确定了目标函数的全局最优解,优于上一节的 PoI。

注意:由于其简单性和在探索和开发之间的自然平衡,EI 是 BayesOpt 中最常用的策略之一,如果没有其他原因可以优先选择其他策略,该策略是一个很好的默认选择。

4.4 练习

本章中有两个练习:

  1. 第一个练习涵盖了使用 PoI 进行探索,通过改变我们对改进的定义。

  2. 第二个示例涵盖了使用 BayesOpt 进行超参数调整,目标函数模拟了 SVM 模型的精度曲面。

4.4.1 练习 1:通过 PoI 促进探索

解决 PoI 过度开发的一个方法是设定更高的标准来确定什么是改进。具体来说,我们发现单纯地找到最大化从现状改进的可能性的点会阻止我们摆脱局部最优解。

作为解决方案,我们可以修改策略以指定我们仅接受至少 ε 的改进。这将指导 PoI 在局部区域已经足够覆盖后,在搜索空间中寻找改进的其他区域。在 CH04/02 - Exercise 1.ipynb 中实现此练习,并演示其对 PoI 的积极影响。其步骤如下:

  1. 重新创建 CH04/01 - BayesOpt loop.ipynb 中的 BayesOpt 循环,其中将一维的 Forrester 函数作为优化目标。

  2. 在实现 BayesOpt 的 for 循环之前,声明一个名为 epsilon 的变量。此变量将作为鼓励探索的最小改进阈值。暂时将此变量设置为 0.1。

  3. for 循环内,像以前一样初始化 PoI 策略,但这次指定由 best_f 参数设置的现任阈值为现任值 加上 存储在 epsilon 中的值。

  4. 重新运行笔记本,并观察此修改是否比原始 PoI 策略更好地促进了更多的探索。

  5. PoI 变得更具探索性的程度在很大程度上取决于存储在 epsilon 中的最小改进阈值。将此变量设置为 0.001,观察当改进阈值不够大时是否能成功促进探索。当此值设置为 0.5 时会发生什么?

  6. 在上一步中,我们看到将改进阈值设置为适当的值对 PoI 非常关键。然而,如何在多个应用和目标函数中进行这样的设置并不明显。一个合理的启发式方法是将其动态设置为现任值的某个 α 百分比,指定我们希望看到现任值增加 1 + α。在代码中实现这一点,并设置 110% 的改进要求。

4.4.2 练习 2:用于超参数调整的 BayesOpt

此练习,在 CH04/03 - Exercise 2.ipynb 中实现,将 BayesOpt 应用于模拟 SVM 模型超参数调整任务的准确度表面的目标函数。x-轴表示惩罚参数 C 的值,而 y-轴表示 RBF 核参数 γ 的值。有关更多详细信息,请参阅第三章中的练习。步骤如下:

  1. 重新创建 CH04/01 - BayesOpt loop.ipynb 中的 BayesOpt 循环:

    1. 我们不再需要 Forrester 函数;相反,复制第三章中练习描述的二维函数的代码,并将其用作目标函数。

    2. 请注意,此函数的定义域为 [0, 2]²。

  2. 使用 xs 声明相应的测试数据,表示域的二维网格,ys 表示 xs 的函数值。

  3. 修改可视化优化进展的辅助函数。对于一维目标函数,很容易可视化 GP 预测以及收购积分。对于这个二维目标,该辅助函数应生成两个面板的绘图:一个显示基本事实,另一个显示获得分数。两个面板也应显示标记数据。绘图应类似于图 4.18。

    图 4.18 一个参考,显示可视化 BayesOpt 进展的辅助函数应该是什么样子。左面板显示真实的目标函数,而右面板显示收购分数。

  4. 从第三章的练习中复制 GP 类,该类使用 ARD 实现 Matérn 2.5 核。进一步修改此类,使其与 BoTorch 集成。

  5. 复用辅助函数 fit_gp_model() 和实现 BayesOpt 的 for 循环:

    1. 初始培训数据集应包含域中心的点:(1, 1)。

    2. 由于我们的搜索空间是二维的,请在 botorch.optim.optimize_acqf() 中设置 num_restarts=40raw_samples=100,以更详尽地搜索最大化收购分数的点。

    3. 将我们可以查询的数量设置为 20(我们可以评估目标函数的次数)。

  6. 对这个目标函数运行 PoI 策略。注意到该策略再次陷入了局部最优。

  7. 运行修改后的 PoI,其中最小改进阈值设置为 0.1:

    1. 参见练习 1,了解有关为 PoI 设置最小改进阈值的更多详细信息。

    2. 注意到此修改再次导致更好的优化性能。

    3. 在哪一轮迭代中达到了至少 90% 准确度?实现此准确度的模型参数是什么?

  8. 对这个目标函数运行 EI 策略:

    1. 注意到这种策略优于 PoI。

    2. 在达到至少 90% 的准确度的第一轮迭代是什么?实现此准确度的模型参数是什么?

  9. 检查基于单次 BayesOpt 运行的策略的性能可能会引导错误结论。更好地进行多次 BayesOpt 实验,使用不同的起始数据重复实验:

    1. 实现重复实验的想法,并可视化 10 个实验中的平均入职价值和误差线。

    2. 每个实验都应从搜索空间中均匀抽样的单个数据点开始。

    3. 运行我们列出的策略,比较它们的性能。

这标志着我们关于改进的基于 BayesOpt 策略的章节的结束。重要的是要记住在这里使用 Forrester 函数实现 BayesOpt 循环的代码,因为我们将在未来的章节中再次使用它来基准测试其他策略。具体来说,在第五章中,我们将了解受多臂老虎机问题启发的 BayesOpt 策略。

总结

  • 贝叶斯优化策略使用训练后的高斯过程来评分每个数据点在我们寻找目标函数最优解的过程中的价值。策略计算的得分被称为获取得分。

  • 在贝叶斯优化循环的每次迭代中,基于观察到的数据来训练高斯过程,策略建议一个新的数据点来查询,并且将此点的标签添加到训练集中。这个过程重复进行,直到我们不能再进行任何函数评估为止。

  • GPyTorch 模型只需要进行最小的修改就可以集成到实现贝叶斯优化策略的 BoTorch 中。

  • BoTorch 提供了一个名为optimize_acqf()的帮助函数,该函数来自优化模块optim .optimize,它接受一个策略对象,并返回最大化获取得分的数据点。

  • 一个好的贝叶斯优化策略需要平衡探索(学习高度不确定的区域)和利用(缩小高绩效区域)。

  • 不同的贝叶斯优化策略以不同的方式处理探索-利用的权衡。检查优化进度以分析和调整所使用策略的性能是很重要的。

  • 在贝叶斯优化中可能使用的一种启发式方法是找到比最佳值更好的点。

  • 寻找能够最有可能从最佳值改进的点得到了 PoI 策略。

  • 寻找从最佳值改进的期望提高最高的点得到了 EI 策略。

  • PoI 可以被认为是一种过度的利用和风险规避的策略,因为该策略仅仅旨在从最佳值改进,无论改进多小。没有任何进一步的修改,EI 往往比 PoI 更好地平衡探索和利用。

  • 由于对函数值的高斯信念,通过 PoI 和 EI 计算分数可以在闭合形式下完成。因此,我们可以轻松地计算和优化这些策略定义的分数。