openCypher* 针对任何关系数据库
原文：
towardsdatascience.com/opencypher-against-any-relational-database-a3b2388579df
关系数据库作为图形数据库 = Mindful (openCypher-2-SQL)
·发表于 数据科学前沿 ·8 分钟阅读·2023 年 7 月 25 日
–
图片由作者提供。阴阳月。由 Syed Ahmad 修改的免版税照片 Unsplash
在任何关系数据库上运行的 openCypher 图查询的有限子集是 Mindful 计划。此阶段的查询是只读的，且不包含元图查询。Mindful 是 微软的 openCypher 到 SQL 转换器 的封闭源代码修改版，遵循 MIT 许可证，其中 Mindful 生成 SQL 以操作任何关系/SQL 数据库。
考虑到这一点……让我们开始了解范围……
在 Mindful 的背景下，“任何关系数据库”意味着 openCypher 查询被转换为针对任何实际关系数据库的 SQL，而不是那些需要为图类型查询特别修改表格或将数据作为 JSON 注入字段并在该 JSON 数据上执行图形查询的关系数据库。
openCypher 查询被转换为 SQL，以便在任何标准关系数据库上运行。
解释视频：
对您业务的适用性 — 数据科学
您可能已经拥有一个现有的数据仓库、语义层或本质上是关系型的数据库，并且使用 SQL 作为主要查询语言……而您希望使用图形查询来查询您的数据资产。
相反，你可能急需从现有的图数据库迁移到关系/SQL 基于数据库，并需要数据迁移测试和实施工具。Mindful，一个 openCypher 到 SQL 的转译器，旨在成为实现你目标的工具。
现有的在关系数据库上进行图查询的实现需要特殊的表来有效表示节点类型和边类型（例如，具有单列主键的表）。Mindful 实现允许你在具有多列主键的表上运行 filter openCypher 查询。
在这篇文章中，我们展示了如何在不影响现有关系数据库栈的情况下实现这一点，通过采用一种数据科学策略，在 DDL 的评论部分使用 JSON 存储有关关系数据的同态图结构的元信息。例如，ORACLE、SQL Server、PostgreSQL 和大多数主流关系数据库都支持对元数据/DDL 进行注释。
实际上，通过使用现有数据库/数据仓库或语义层中提供的标准工具对数据库进行文档化，你将你的关系数据栈转换为图兼容的数据栈。
让我们深入探讨使这一切成为可能的数据科学和数据库理论……
关注细节
现代数据库管理系统显然具有以下特点，我们将对其进行探讨：
任何关系数据库都可以视为图数据库或关系数据库，取决于你的观点（如下所示）；以及
一旦数据库投入使用，数据库架构随处可见。例如，甚至连以前提倡无架构的图数据库业务如Neo4j 也开始谈论“架构”。
考虑到这一定位，我们可以想象，未来所有有价值的数据库都将提供图查询语言以及普通结构化查询语言（SQL）查询。
任何关系数据库都可以视作图数据库
图-关系范式、关系知识图谱 和多模型数据库并不新鲜，但为了设定背景，请注意实体-关系图与其对应的属性图模式及其反向映射：
同态映射——属性图模式到实体-关系图。图像作者提供。
我们的架构用于在电影院预订座位以观看电影。一个人可以为在电影院观看电影的一个或多个座位进行预订；其中座位位于该电影院的某个区域和行中。
Shutterstock 图片，2187621947. 授权给 Victor Morgante / FactEngine。 免版税标准许可证 2187621947在 Shutterstock 上。
数据库理论…
使 Mindful 成功的两个方面是：
-
外键关系和多对多表在其他关系数据库中与图数据库的边类型是同义的；
-
图查询语言如 openCypher 支持可以用类似于 SQL 查询的方式编写的查询。
让我们单独查看这些内容…
外键和多对多表
外键 <-> 边类型
关系数据库中的外键关系与属性图模式中的对应边类型有 1:1 的映射关系。
在我们的模式中，每个名为 Row 的表中的 Row 记录都有一个指向‘Cinema’表中 Cinema 记录的外键引用，表示该行‘IS IN’的影院。
实体关系图中的外键引用。图片由作者提供。
…在我们的属性图模式中有一个对应的边类型。
属性图模式中的边类型。图片由作者提供。
…在我们的 DDL（数据库定义语言）中，这可以通过一个表示属性图模式的“IS_IN”标签的 JSON（JavaScript 对象表示法）轻松表示，嵌入为 CREATE TABLE DDL 语句中‘Row’表的对应外键参考定义的注释/文档。
CREATE TABLE "Row" (Cinema_Id INTEGER REFERENCES [Cinema] NOT NULL
,RowNr NUMBER NOT NULL
, CONSTRAINT [Row_PK] PRIMARY KEY ([Cinema_Id],[RowNr])
, FOREIGN KEY ([Cinema_Id]) REFERENCES [Cinema] ([Cinema_Id])
ON DELETE CASCADE ON UPDATE CASCADE /* { Label:"IS_IN"} */
)
…我们可以将关系视图和图视图之间的同态视为一种变形动画：
外键关系与边类型之间的同态。图片由作者提供。
多对多表 <-> 边类型
类似地，关系数据库中的多对多连接表在传统图数据库的属性图模式中与其对应的边类型具有同态。
在我们的模式中，一个预订可能有多个座位，而一个座位（在其生命周期内）可能有多个预订。
注意，在我们的模式中，‘Booking’和‘Seat’表的主键是多列的，而我们在‘BookingHasSeat’多对多连接表上的主键有七列。
在 ER 图中的多对多连接表。图片由作者提供。
…在属性图模式视图中，Booking 和 Seat 节点类型之间有一个标记为‘HAS’的对应边类型。
属性图模式中的边类型。图片由作者提供。
…我们可以很容易地将‘HAS’标签作为 JSON 捕获在 DDL 中‘BookingHasSeat’表主键的注释中。
CREATE TABLE "BookingHasSeat" (CinemaId INTEGER
,Letter TEXT(1)
,DateTime DATETIME
,RowNr NUMBER
,Person_Id INTEGER
,Cinema_Id INTEGER
,Film_Id INTEGER
,CONSTRAINT [BookingHasSeat_PK] PRIMARY KEY
([CinemaId],[Letter],[DateTime],[RowNr],[Person_Id],[Cinema_Id],[Film_Id])
/* { Label:"HAS"} */
,FOREIGN KEY ([Person_Id],[Film_Id],[DateTime],[CinemaId]) REFERENCES
[Booking] ([Person_Id],[Film_Id],[DateTime],[Cinema_Id])
ON DELETE CASCADE ON UPDATE CASCADE /* { Label:"HAS"} */
,FOREIGN KEY ([RowNr],[Cinema_Id],[Letter]) REFERENCES
[Seat] ([RowNr],[Cinema_Id],[Letter])
ON DELETE CASCADE ON UPDATE CASCADE /* { Label:"HAS"} */
)
…我们也可以将关系视图和图视图之间的同态视为一种变形动画：
架构
Mindful 的实现变得相当简单易懂。通过将适当的同态从关系模式映射到图模式，并且我们可以将有关同态的元信息存储在关系数据库的 DDL 中，我们有一种非常简单的机制来编写 openCypher 查询……因为实际上，关系数据库的元模型与图数据库的元模型之间接近同构。即，从适当的视角来看，它们在概念上几乎是相同的。
同态 — openCypher 到 SQL
…两种查询类型的故事
图查询语言有两种查询类型：
-
过滤查询；以及
-
元图/图遍历查询。
元图查询是一种在属性图模式的“图”（或模型）中查找路径、关系或结构的查询类型。
另一方面，过滤查询从数据库的所有数据的超集中过滤数据。例如，用自然语言表达：“Peter 预订了哪些影院？”
在这一阶段，我们最感兴趣的是过滤查询，因为这是 Mindful 代码所处的阶段，但一个元图查询可能会提出这样的问题，其他形式的自然语言表述为：“从一个人 Peter 到他预订的影院的最短路径是什么？”。
从图像上看，我们可以看到，对于我们的模式，最短路径是从人到预订，再到会话，再到影院，如下所示（而不是人-预订-座位-行-影院）：
最短路径 — 人到影院 — 元图/图遍历。图像作者提供。
对于我们的目的……我们需要过滤查询，因为 openCypher 过滤查询与 SQL 查询有同态，而 标准 SQL 不支持元图查询。
让我们使用 openCypher 查询
假设我们想知道登录名为“Peter”的人在哪个日期和时间在什么影院观看了哪个电影，以及座位字母和行号……我们可以为我们的模式编写以下 openCypher 查询：
**MATCH (p:Person)<-[:IS_FOR]-(b:Booking)-[:HAS]->(seat:Seat)-[:IS_IN]->(:Row)-[:IS_IN]->(c:Cinema), (b)-[:IS_FOR]->(s:Session)-[:IS_FOR]->(f:Film)
WHERE p.LoginName = ‘Peter’
RETURN p.LoginName, f.Name, s.DateTime, seat.Letter, seat.RowNr, c.CinemaName
这是一种过滤类型的查询，因为我们正在将结果筛选到那些登录名为“Peter”的人。
…现在让我们在关系数据库上运行这个查询，使用连接到 Mindful 的软件：
openCypher 查询与关系数据库 — 结果。图像作者提供。
在这个阶段，微软的 openCypher 到 SQL 转译器（以及 Mindful）生成的 SQL 相当繁琐，我在这里不会展示，但足以说明，相比于 openCypher 查询，相应的 SQL 查询相当长。以下是我们模式的等效 SQL 查询：
SELECT Person.[LoginName],Film.[Name],Session.[DateTime],
Seat.Letter, Seat.RowNr,Cinema.[CinemaName]
FROM Person,
Booking,
Session,
Film,
Seat,
Row,
Cinema
,BookingHasSeat
WHERE Booking.Person_Id = Person.Person_Id
AND Booking.Film_Id = Session.Film_Id
AND Booking.DateTime = Session.DateTime
AND Booking.Cinema_Id = Session.Cinema_Id
AND Session.Film_Id = Film.Film_Id
AND BookingHasSeat.Person_Id = Booking.Person_Id
AND BookingHasSeat.Film_Id = Booking.Film_Id
AND BookingHasSeat.DateTime = Booking.DateTime
AND BookingHasSeat.CinemaId = Booking.Cinema_Id
AND BookingHasSeat.RowNr = Seat.RowNr
AND BookingHasSeat.Cinema_Id = Seat.Cinema_Id
AND BookingHasSeat.Letter = Seat.Letter
AND Seat.Cinema_Id = Row.Cinema_Id
AND Seat.RowNr = Row.RowNr
AND Row.Cinema_Id = Cinema.Cinema_Id
AND Person.LoginName = 'Peter'
总结
所以，这就是如此简单。我们展示了关系模式和图模式之间的同态，使得可以编写软件，以便在任何适当配置的关系数据库上执行 openCypher 查询（过滤类型），并且 openCypher 被转换为 SQL。
感谢阅读。时间允许的话，我将写更多关于图数据库、关系数据库、图查询和同态的内容。
— — — — — — —
Mindful openCypher 查询目前为只读，并且在此阶段没有元图查询。
归属：本文中的模式源于最初在 ActiveFacts examples github 页面 下以 MIT 许可证 开源共享的内容。
— — — — — — — — — — — — — 结束 — — — — — — — — — — —
通过物理启发的 DeepONet 进行算子学习：从头开始实现
深入探讨 DeepONets、物理启发的神经网络以及物理启发的 DeepONets
·发表于 Towards Data Science ·阅读时间 23 分钟·2023 年 7 月 7 日
–
图 1. ODE/PDEs 广泛用于描述系统过程。在许多情况下，这些 ODE/PDEs 接受一个函数（例如，强迫函数 u(t)）作为输入，并输出另一个函数（例如，s(t)）。传统上，数值求解器用于连接输入和输出。最近，神经算子的开发大大提高了处理效率。（图像由作者提供）
常微分方程（ODEs）和偏微分方程（PDEs）是许多科学和工程学科的基础，从物理学和生物学到经济学和气候科学。它们是描述物理系统和过程的基本工具，捕捉了数量随时间和空间的连续变化。
然而，这些方程的一个独特特点是它们不仅接受单一数值作为输入，还接受函数。例如，考虑预测建筑物因地震而产生的振动。地面的震动随时间变化，可以表示为一个函数，该函数作为描述建筑物运动的微分方程的输入。同样，在音乐厅中声波传播的情况下，乐器产生的声波可以是一个随时间变化的音量和音调的输入函数。这些变化的输入函数从根本上影响了结果输出函数——建筑物的振动和音乐厅的声学场。
传统上，这些 ODEs/PDEs 通常使用有限差分或有限元方法等数值求解器来解决。然而，这些方法存在一个瓶颈：每当有新的输入函数时，求解器必须重新运行一次。这个过程可能计算密集且缓慢，尤其是在复杂系统或高维输入情况下。
为了应对这一挑战，Deep Operator Network（简称DeepONet）的创新框架由Lu et al.于 2019 年提出。DeepONets 旨在学习将输入函数映射到输出函数的算子，本质上是学习预测这些 ODEs/PDEs 在任意给定输入函数下的输出，而不需要每次都重新运行数值求解器。
尽管 DeepONets 很强大，但它们继承了数据驱动方法面临的共同问题：我们如何确保网络的预测与包含在控制方程中的已知物理定律一致？
进入物理信息化学习领域。
物理信息化学习是一个迅速发展的机器学习分支，它将物理原理与数据科学结合起来，以增强对复杂物理系统的建模和理解。它涉及利用特定领域的知识和物理定律来指导学习过程，提高机器学习模型的准确性、泛化能力和可解释性。
在这个框架下，2021 年，Wang et al.提出了 DeepONets 的新变种：物理信息化 DeepONet。这种创新方法在 DeepONets 的基础上，通过将我们对物理定律的理解融入学习过程中，进行改进。我们不再只是让模型从数据中学习，而是用源于几个世纪科学探究的原理来指导它。
这看起来是一个非常有前景的方法！但是我们应该如何在实践中实现它？这正是我们今天要探讨的内容🤗
在这篇博客中，我们将讨论物理信息化 DeepONet 背后的理论，并逐步讲解如何从零开始实现它。我们还将把我们开发的模型付诸实践，通过实际案例展示其强大能力。
如果你也有兴趣使用物理信息化 DeepONet 解决逆问题，可以查看我的新博客：利用物理信息化 DeepONet 解决逆问题：带代码实现的实用指南
让我们开始吧！
内容表
· 1. 案例研究
· 2. 物理信息化 DeepONet
∘ 2.1 DeepONet：概述
∘ 2.2 物理信息化神经网络（PINNs）
∘ 2.3 物理信息化 DeepONet
· 3. 物理信息化 DeepONet 的实现
∘ 3.1 定义架构
∘ 3.2 定义 ODE 损失
∘ 3.3 定义梯度下降步骤
· 4. 数据生成与组织
∘ 4.1 u(·) 轮廓生成
∘ 4.2 数据集生成
∘ 4.3 数据集组织
· 5. 训练物理信息深度运算网络
· 6. 结果讨论
· 7. 重点总结
· 参考文献
1. 案例研究
让我们在一个具体的例子中扎根讨论。在这篇博客中，我们将重现 Wang et al. 原论文中考虑的第一个案例研究，即由以下常微分方程（ODE）描述的初值问题：
具有初始条件 s(0) = 0。
在这个方程中，u(t) 是随时间变化的输入函数，而 s(t) 是我们感兴趣的在时间 t 的系统状态。在物理场景中，u(t) 可能代表施加在系统上的力，而 s(t) 可能代表系统的响应，比如位移或速度，具体取决于上下文。我们这里的最终目标是学习强迫项 u(t) 与 ODE 解 s(t) 之间的映射关系。
传统的数值方法，如欧拉方法或龙格-库塔方法，可以有效地求解此方程。然而，请注意，强迫项 u(t) 可以采取各种轮廓，如下图所示：
图 2. u(t) 的示例轮廓。 （作者提供的图片）
因此，每当 u(t) 变化时，我们需要重新运行整个模拟以获取相应的 s(t)（如图 3 所示），这可能会计算密集且效率低下。
图 3. s(t) 的相应轮廓。它们是通过使用 RK45 算法求解 ODE 计算得出的。 （作者提供的图片）
那么，我们如何更高效地解决这种问题呢？
2. 物理信息深度运算网络
如介绍中所述，物理信息 DeepONet 是解决我们目标问题的有前途的解决方案。在这一部分，我们将详细解析其基本概念，使其更易于理解。
我们将首先讨论原始 DeepONet 的基础原则。接着，我们将探索物理信息神经网络的概念及其如何为问题解决提供额外的维度。最后，我们将展示如何将这两个想法无缝集成以构建物理信息 DeepONet。
2.1 DeepONet：概述
DeepONet，简而言之就是深度运算网络，代表了深度学习的新前沿。与传统的机器学习方法将一组输入值映射到输出值不同，DeepONet 旨在将整个函数映射到其他函数。这使得 DeepONet 在处理自然涉及函数输入和输出的问题时特别强大。那么它是如何实现这一目标的呢？
为了符号化我们想要实现的目标：
图 4. 我们的目标是训练一个神经网络，以近似将强迫项 u(·)映射到目标输出 s(·)的算子，这两者都是时间的函数。（图片由作者提供）
左边是将输入函数 u(·)映射到输出函数 s(·)的算子G。右边，我们希望使用神经网络来近似算子 G。一旦实现了这一点，我们可以利用训练好的神经网络在给定任何 u(·)的情况下快速计算 s(·)。
对于当前的案例研究，输入函数 u(·)和输出函数 s(·)都将时间坐标t作为唯一参数。因此，我们希望构建的神经网络的“输入”和“输出”应如下所示：
图 5. 我们希望训练的神经网络模型的输入和输出。（图片由作者提供）
实质上，我们的神经网络应接受 u(t)的整个轮廓作为第一个输入，以及一个特定时间点t作为第二个输入。随后，它应输出在时间点t评估的目标输出函数 s(·)，即 s(t)。
为了更好地理解这一设置，我们认识到 s(t)的值首先依赖于 s(·)的轮廓，而 s(·)的轮廓又依赖于 u(·)，其次依赖于 s(·)被评估的时间点。这也是时间坐标t需要作为输入之一的原因。
目前我们需要弄清楚两个问题：首先，我们应该如何将 u(·)的连续轮廓输入网络？其次，我们应该如何拼接这两个输入，即t和 u(·)。
1️⃣ 我们应该如何输入 u(·)的连续轮廓？
实际上，我们并不这样做。一种直接的解决方案是离散表示函数 u(·)。更具体地说，我们只是评估 u(·)在足够但有限的多个位置的值，然后将这些离散的 u(·)值输入到神经网络中：
图 6. 在被输入到神经网络之前，u(·)轮廓被离散化。（图片由作者提供）
这些位置在原始 DeepONet 论文中被称为传感器。
2️⃣ 我们应该如何将输入t和 u(·)拼接在一起？
初看之下，我们可能会想直接在输入层将它们拼接起来。然而，事实证明，这种简单的方法不仅会限制我们可以使用的神经网络类型，而且在实践中会导致次优的预测准确度。不过，还有更好的方法。现在是介绍DeepONet的时候了。
总之，DeepONet 提出了一种用于算子学习的新网络架构：它由两个主要组件组成：分支网络和主干网络。分支网络将离散函数值作为输入，并将其转换为特征向量。同时，主干网络将坐标（在我们当前的案例研究中，坐标仅为t。对于 PDE，将包括时间和空间坐标）作为输入，并将其也转换为相同维度的特征向量。这两个特征向量然后通过点积进行合并，最终结果用作在输入坐标处评估 s(·)的预测值。
图 7. DeepONet 包括一个分支网络来处理输入函数 u(·)和一个主干网络来处理时间/空间坐标。两个网络的输出具有相同的维度，并通过点积进行合并。可选地，可以在点积后添加一个偏置项以进一步提高模型的表达能力。（图片由作者提供）
在原始 DeepONet 论文中，作者指出，这种在“分支”和“主干”网络中体现的“分而治之”策略受到算子通用逼近定理的启发，旨在为算子学习引入强的归纳偏置。这也是作者声称使 DeepONet 成为有效解决方案的关键点。
如果你想了解更多关于 DeepONet 理论基础的内容，请参考原始论文的附录 A。
DeepONet 的主要优势之一是其效率。一旦训练完成，DeepONet 可以实时推断新的输入函数的输出函数，无需进一步训练，只要新的输入函数在其训练过的输入函数范围内。这使 DeepONet 成为需要实时推断的应用中的强大工具。
DeepONet 的另一个显著优势在于其灵活性和多功能性。虽然主干网络和分支网络最常见的选择是全连接层，但 DeepONet 框架允许高度的架构自定义。根据输入函数 u(·)和坐标的特征，可以采用各种神经网络架构，如 CNN、RNN 等。这种适应性使 DeepONet 成为一个高度多功能的工具。
然而，尽管存在这些优势，DeepONet 的局限性也很明显：作为一种纯数据驱动的方法，DeepONet 不能保证其预测结果会遵循描述所考虑物理系统的先验知识或控制方程。因此，DeepONet 可能无法很好地泛化，尤其是当面对位于训练数据分布之外的输入函数，即分布外（OOD）输入时。对此的一个常见解决方案是准备大量训练数据，但在实际中这可能并不总是可行，特别是在数据收集可能昂贵或耗时的科学和工程领域。
那么我们应该如何解决这些局限性呢？是时候讨论物理信息学习，更具体地说，是物理信息神经网络（PINNs）了。
2.2 物理信息神经网络（PINNs）
在传统的机器学习模型中，我们主要依赖数据来学习潜在的模式。然而，在许多科学和工程领域，捕捉我们对动态系统的先验知识的控制方程（ODE/PDE）是可用的，它们提供了除了观察数据之外的另一种信息来源。如果正确地将这一额外的知识源纳入模型中，它可能会改善模型的性能和泛化能力，特别是在处理有限或噪声数据时。这就是物理信息学习的作用所在。
当我们将物理信息学习与神经网络的概念结合时，我们将得到物理信息神经网络（PINNs）。
PINNs 是一种神经网络，其中网络不仅仅是拟合数据，还要尊重由微分方程描述的已知物理定律。这是通过引入ODE/PDE 损失来实现的，它测量了控制微分方程的违反程度。通过这种方式，我们将物理定律注入网络训练过程，使其物理信息化。
图 8. 物理信息神经网络的损失函数包括 PDE 损失的贡献项，这有效地测量了预测解是否满足控制微分方程。注意，由于自动微分的存在，相对于输入的输出的导数可以很容易地计算出来。（图片来源：作者）
尽管 PINNs 在许多应用中已被证明有效，但它们也不是没有局限性。PINNs 通常是针对特定的输入参数（例如边界和初始条件、外部强迫等）进行训练的。因此，每当输入参数发生变化时，我们就需要重新训练 PINN。因此，它们在不同操作条件下的实时推断效果不是特别好。
还记得哪个方法是专门用于处理变化的输入参数的吗？没错，就是 DeepONet！现在是将物理信息学习的理念与 DeepONet 结合的时候了。
2.3 物理信息 DeepONet
物理信息 DeepONet的主要思想是结合 DeepONets 和 PINNs 的优点。就像 DeepONet 一样，物理信息 DeepONet 能够将一个函数作为输入，并产生一个函数作为输出。这使得它在实时推断新输入函数时非常高效，无需重新训练。
另一方面，像 PINN 一样，物理信息 DeepONet 在学习过程中融入了已知的物理定律。这些定律作为额外的约束引入到训练过程中的损失函数中。这种方法使得模型即使在处理有限或嘈杂数据时，也能做出物理一致的预测。
我们如何实现这种整合呢？类似于 PINNs，我们增加一个额外的损失项，以衡量模型的预测如何符合已知的微分方程。通过优化这个损失函数，模型学会进行数据一致（如果在训练过程中提供了测量数据）和物理一致的预测。
图 10. 物理信息 DeepONet 使用 DeepONet 作为骨干架构，同时利用物理信息学习的概念来训练模型。这样，训练后的物理信息 DeepONet 既数据一致又物理一致。（图像由作者提供）
总结来说，物理信息 DeepONet 是一个强大的工具，结合了两者的优势：DeepONet 的高效性和物理信息学习的准确性。它代表了一种有前景的方法，用于解决那些实时推断和物理一致性都至关重要的复杂问题。
在下一部分，我们将开始进行案例研究，并将理论转化为实际代码。
3. 物理信息 DeepONet 的实现
在这一部分，我们将详细讲解如何定义一个物理信息 DeepONet 模型，以解决我们的目标案例研究。我们将使用 TensorFlow 来实现它。让我们先导入必要的库：
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras
tf.random.set_seed(42)
3.1 定义架构
如前所述，物理信息 DeepONet 与原始 DeepONet 具有相同的架构。以下函数定义了 DeepONet 的架构：
def create_model(mean, var, verbose=False):
"""Definition of a DeepONet with fully connected branch and trunk layers.
Args:
----
mean: dictionary, mean values of the inputs
var: dictionary, variance values of the inputs
verbose: boolean, indicate whether to show the model summary
Outputs:
--------
model: the DeepONet model
"""
# Branch net
branch_input = tf.keras.Input(shape=(len(mean['forcing'])), name="forcing")
branch = tf.keras.layers.Normalization(mean=mean['forcing'], variance=var['forcing'])(branch_input)
for i in range(3):
branch = tf.keras.layers.Dense(50, activation="tanh")(branch)
# Trunk net
trunk_input = tf.keras.Input(shape=(len(mean['time'])), name="time")
trunk = tf.keras.layers.Normalization(mean=mean['time'], variance=var['time'])(trunk_input)
for i in range(3):
trunk = tf.keras.layers.Dense(50, activation="tanh")(trunk)
# Compute the dot product between branch and trunk net
dot_product = tf.reduce_sum(tf.multiply(branch, trunk), axis=1, keepdims=True)
# Add the bias
output = BiasLayer()(dot_product)
# Create the model
model = tf.keras.models.Model(inputs=[branch_input, trunk_input], outputs=output)
if verbose:
model.summary()
return model
在上面的代码中：
-
我们假设主干网络和分支网络都是完全连接的网络，每个网络有 3 个隐藏层，每层包含 50 个神经元，并且使用 tanh 激活函数。这个架构是基于初步测试选择的，并且应该作为这个问题的一个良好的起点。然而，可以很容易地用其他架构（例如 CNN、RNN 等）和其他层超参数进行替换。
-
主干网络和分支网络的输出通过点积合并。正如原始 DeepONet 论文中建议的，我们添加了一个偏置项以提高预测准确性。
BiasLayer()是一个自定义定义的类，用于实现这个目标：
class BiasLayer(tf.keras.layers.Layer):
def build(self, input_shape):
self.bias = self.add_weight(shape=(1,),
initializer=tf.keras.initializers.Zeros(),
trainable=True)
def call(self, inputs):
return inputs + self.bias
3.2 定义 ODE 损失
接下来，我们定义一个函数来计算 ODE 损失。回顾一下我们的目标 ODE 是：
因此，我们可以按如下方式定义该函数：
@tf.function
def ODE_residual_calculator(t, u, u_t, model):
"""ODE residual calculation.
Args:
----
t: temporal coordinate
u: input function evaluated at discrete temporal coordinates
u_t: input function evaluated at t
model: DeepONet model
Outputs:
--------
ODE_residual: residual of the governing ODE
"""
with tf.GradientTape() as tape:
tape.watch(t)
s = model({"forcing": u, "time": t})
# Calculate gradients
ds_dt = tape.gradient(s, t)
# ODE residual
ODE_residual = ds_dt - u_t
return ODE_residual
在上面的代码中：
-
我们使用
tf.GradientTape()来计算 s(·)相对于t的梯度。请注意，在 TensorFlow 中，tf.GradientTape()作为上下文管理器使用，任何在 tape 上下文中执行的操作都会被 tape 记录。在这里，我们显式地观察变量t。因此，TensorFlow 会自动跟踪涉及t的所有操作，在这种情况下，它是 DeepONet 模型的前向传播。之后，我们使用 tape 的gradient()方法来计算 s(·)相对于t的梯度。 -
我们包括了一个额外的输入参数
u_t，它表示在t时刻评估的输入函数 u(·)的值。这构成了我们目标 ODE 的右侧项，并且它是计算 ODE 残差损失所需的。 -
我们使用
@tf.function装饰器将我们刚刚定义的常规 Python 函数转换为 TensorFlow 图。这是有用的，因为梯度计算可能非常昂贵，并且在图模式下执行可以显著加速计算。
3.3 定义梯度下降步骤
接下来，我们定义一个函数来编译总损失函数并计算总损失相对于网络模型参数的梯度：
@tf.function
def train_step(X, X_init, IC_weight, ODE_weight, model):
"""Calculate gradients of the total loss with respect to network model parameters.
Args:
----
X: training dataset for evaluating ODE residuals
X_init: training dataset for evaluating initial conditions
IC_weight: weight for initial condition loss
ODE_weight: weight for ODE loss
model: DeepONet model
Outputs:
--------
ODE_loss: calculated ODE loss
IC_loss: calculated initial condition loss
total_loss: weighted sum of ODE loss and initial condition loss
gradients: gradients of the total loss with respect to network model parameters.
"""
with tf.GradientTape() as tape:
tape.watch(model.trainable_weights)
# Initial condition prediction
y_pred_IC = model({"forcing": X_init[:, 1:-1], "time": X_init[:, :1]})
# Equation residual
ODE_residual = ODE_residual_calculator(t=X[:, :1], u=X[:, 1:-1], u_t=X[:, -1:], model=model)
# Calculate loss
IC_loss = tf.reduce_mean(keras.losses.mean_squared_error(0, y_pred_IC))
ODE_loss = tf.reduce_mean(tf.square(ODE_residual))
# Total loss
total_loss = IC_loss*IC_weight + ODE_loss*ODE_weight
gradients = tape.gradient(total_loss, model.trainable_variables)
return ODE_loss, IC_loss, total_loss, gradients
在上面的代码中：
-
我们只考虑两个损失项：与初始条件相关的损失
IC_loss和 ODE 残差损失ODE_loss。IC_loss通过将模型预测的 s(t=0)与已知的初始值 0 进行比较来计算，ODE_loss通过调用我们之前定义的ODE_residual_calculator函数来计算。如果有可用的测量 s(t)值（在上面的代码中未实现），数据损失也可以计算并添加到总损失中。 -
通常，总损失是
IC_loss和ODE_loss的加权和，其中权重控制在训练过程中对这些单独损失项的重视程度或优先级。在我们的案例研究中，将IC_weight和ODE_weight都设置为 1 就足够了。 -
类似于我们计算
ODE_loss的方式，我们也采用了tf.GradientTape()作为上下文管理器来计算梯度。然而，这里我们计算的是总损失相对于网络模型参数的梯度，这对于执行梯度下降是必要的。
在继续之前，让我们快速总结一下我们到目前为止所开发的内容：
1️⃣ 我们可以使用create_model()函数初始化一个 DeepONet 模型。
2️⃣ 我们可以计算 ODE 残差，以评估模型预测与所治理 ODE 的契合程度。这是通过ODE_residual_calculator函数实现的。
3️⃣ 我们可以使用train_step计算总损失及其相对于网络模型参数的梯度。
现在准备工作完成了一半🚀 在下一节中，我们将讨论数据生成和数据组织的问题（上述代码中的奇怪X[:, :1]会在那时变得清晰）。之后，我们终于可以训练模型并查看其表现。
4. 数据生成和组织
在本节中，我们讨论合成数据的生成及其在训练 Physics-informed DeepONet 模型中的组织方式。
4.1 生成 u(·)特征
用于训练、验证和测试的数据将是合成生成的。这样做的理由有两个：不仅方便，而且可以完全控制数据的特征。
在我们的案例研究中，我们将使用零均值的高斯过程生成输入函数 u(·)，并使用径向基函数（RBF）核。
高斯过程是一种强大的数学框架，常用于机器学习中建模函数。RBF 核是捕捉输入点之间相似性的热门选择。通过在高斯过程中使用 RBF 核，我们确保生成的合成数据表现出平滑和连续的模式，这在各种应用中通常是有利的。如需了解更多关于高斯过程的内容，请随时查看我之前的博客。
在 scikit-learn 中，这可以通过几行代码实现：
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF
def create_samples(length_scale, sample_num):
"""Create synthetic data for u(·)
Args:
----
length_scale: float, length scale for RNF kernel
sample_num: number of u(·) profiles to generate
Outputs:
--------
u_sample: generated u(·) profiles
"""
# Define kernel with given length scale
kernel = RBF(length_scale)
# Create Gaussian process regressor
gp = GaussianProcessRegressor(kernel=kernel)
# Collocation point locations
X_sample = np.linspace(0, 1, 100).reshape(-1, 1)
# Create samples
u_sample = np.zeros((sample_num, 100))
for i in range(sample_num):
# sampling from the prior directly
n = np.random.randint(0, 10000)
u_sample[i, :] = gp.sample_y(X_sample, random_state=n).flatten()
return u_sample
在上面的代码中：
-
我们使用
length_scale来控制生成函数的形状。对于 RBF 核，图 11 展示了不同核长度尺度下的 u(·)特征。 -
请记住，我们需要在将 u(·)输入 DeepONet 之前对其进行离散化。这是通过指定
X_sample变量来完成的，该变量在我们感兴趣的时间域内分配 100 个均匀分布的点。 -
在 scikit-learn 中，
GaussianProcessRegressor对象提供了一个sample_y方法，用于从具有长度尺度指定核的高斯过程抽取随机样本。注意，我们在使用GaussianProcessRegressor对象之前并没有调用.fit()，这与我们通常对其他 scikit-learn 回归器的做法不同。这是故意的，因为我们希望GaussianProcessRegressor使用我们提供的精确length_scale。如果你调用.fit()，length_scale将被优化为另一个值以更好地拟合给定的数据。 -
输出
u_sample是一个维度为 sample_num * 100 的矩阵。u_sample的每一行表示一个 u(·)的特征，其中包含 100 个离散值。
图 11. 在不同核长度尺度下的合成 u(·) 轮廓。（图片来源：作者）
4.2 数据集生成
现在我们已经生成了 u(·) 轮廓，让我们关注如何组织数据集，以便它可以输入到 DeepONet 模型中。
请记住，我们在上一节中开发的 DeepONet 模型需要 3 个输入：
-
时间坐标 t，这是介于 0 和 1 之间的标量（暂时不考虑批量大小）；
-
u(·) 的轮廓，这是一个由在预定义的、固定的时间坐标（介于 0 和 1 之间）下评估的 u(·) 值组成的向量；
-
u(t) 的值，这也是一个标量。这个 u(t) 值用于在时间坐标 t 下计算 ODE 损失。
因此，我们可以这样构建一个单一的样本：
（图片来源：作者）
当然，对于每个 u(·) 轮廓（在上图中标记为绿色），我们应考虑多个 t（及其对应的 u(t)）来评估 ODE 损失，以更好地施加物理约束。理论上，t 可以取考虑的时间域内的任何值（即在我们案例研究中为 0 和 1 之间）。然而，为了简化，我们只考虑在 u(·) 轮廓离散化的相同时间位置的 t。因此，我们更新后的数据集将如下所示：
（图片来源：作者）
请注意，上述讨论仅考虑了单一的 u(·) 轮廓。如果我们考虑所有的 u(·) 轮廓，我们的最终数据集将如下所示：
（图片来源：作者）
其中 N 代表 u(·) 轮廓的数量。现在有了这个前提，让我们看看一些代码：
from tqdm import tqdm
from scipy.integrate import solve_ivp
def generate_dataset(N, length_scale, ODE_solve=False):
"""Generate dataset for Physics-informed DeepONet training.
Args:
----
N: int, number of u(·) profiles
length_scale: float, length scale for RNF kernel
ODE_solve: boolean, indicate whether to compute the corresponding s(·)
Outputs:
--------
X: the dataset for t, u(·) profiles, and u(t)
y: the dataset for the corresponding ODE solution s(·)
"""
# Create random fields
random_field = create_samples(length_scale, N)
# Compile dataset
X = np.zeros((N*100, 100+2))
y = np.zeros((N*100, 1))
for i in tqdm(range(N)):
u = np.tile(random_field[i, :], (100, 1))
t = np.linspace(0, 1, 100).reshape(-1, 1)
# u(·) evaluated at t
u_t = np.diag(u).reshape(-1, 1)
# Update overall matrix
X[i*100:(i+1)*100, :] = np.concatenate((t, u, u_t), axis=1)
# Solve ODE
if ODE_solve:
sol = solve_ivp(lambda var_t, var_s: np.interp(var_t, t.flatten(), random_field[i, :]),
t_span=[0, 1], y0=[0], t_eval=t.flatten(), method='RK45')
y[i*100:(i+1)*100, :] = sol.y[0].reshape(-1, 1)
return X, y
在上述代码中，我们添加了一个选项，用于计算给定 u(·) 轮廓的相应 s(·)。虽然我们在训练中不会使用 s(·) 值，但我们仍然需要它们来测试模型性能。 s(·) 的计算是通过使用 scipy.integrate.solve_ivp 实现的，这是一个来自 SciPy 的 ODE 求解器，专门设计用于解决初值问题。
现在我们可以生成训练、验证和测试数据集。请注意，对于本案例研究，我们将使用 0.4 的长度尺度生成 u(·) 轮廓，并训练物理信息 DeepONet。
# Create training dataset
N_train = 2000
length_scale_train = 0.4
X_train, y_train = generate_dataset(N_train, length_scale_train)
# Create validation dataset
N_val = 100
length_scale_test = 0.4
X_val, y_val = generate_dataset(N_val, length_scale_test)
# Create testing dataset
N_test = 100
length_scale_test = 0.4
X_test, y_test = generate_dataset(N_test, length_scale_test, ODE_solve=True)
4.3 数据集组织
最后，我们将 NumPy 数组转换为 TensorFlow 数据集对象，以便于数据输入。
# Determine batch size
ini_batch_size = int(2000/100)
col_batch_size = 2000
# Create dataset object (initial conditions)
X_train_ini = tf.convert_to_tensor(X_train[X_train[:, 0]==0], dtype=tf.float32)
ini_ds = tf.data.Dataset.from_tensor_slices((X_train_ini))
ini_ds = ini_ds.shuffle(5000).batch(ini_batch_size)
# Create dataset object (collocation points)
X_train = tf.convert_to_tensor(X_train, dtype=tf.float32)
train_ds = tf.data.Dataset.from_tensor_slices((X_train))
train_ds = train_ds.shuffle(100000).batch(col_batch_size)
# Scaling
mean = {
'forcing': np.mean(X_train[:, 1:-1], axis=0),
'time': np.mean(X_train[:, :1], axis=0)
}
var = {
'forcing': np.var(X_train[:, 1:-1], axis=0),
'time': np.var(X_train[:, :1], axis=0)
}
在上面的代码中，我们创建了两个不同的数据集：一个用于评估 ODE 损失（train_ds），另一个用于评估初始条件损失（ini_ds）。我们还预先计算了 t 和 u(·) 的均值和方差。这些值将用于标准化输入。
数据生成和组织的部分就到这里。接下来，我们将启动模型训练并查看其表现。
5. 训练物理信息 DeepONet
作为第一步，让我们创建一个自定义类来跟踪损失演变：
from collections import defaultdict
class LossTracking:
def __init__(self):
self.mean_total_loss = keras.metrics.Mean()
self.mean_IC_loss = keras.metrics.Mean()
self.mean_ODE_loss = keras.metrics.Mean()
self.loss_history = defaultdict(list)
def update(self, total_loss, IC_loss, ODE_loss):
self.mean_total_loss(total_loss)
self.mean_IC_loss(IC_loss)
self.mean_ODE_loss(ODE_loss)
def reset(self):
self.mean_total_loss.reset_states()
self.mean_IC_loss.reset_states()
self.mean_ODE_loss.reset_states()
def print(self):
print(f"IC={self.mean_IC_loss.result().numpy():.4e}, \
ODE={self.mean_ODE_loss.result().numpy():.4e}, \
total_loss={self.mean_total_loss.result().numpy():.4e}")
def history(self):
self.loss_history['total_loss'].append(self.mean_total_loss.result().numpy())
self.loss_history['IC_loss'].append(self.mean_IC_loss.result().numpy())
self.loss_history['ODE_loss'].append(self.mean_ODE_loss.result().numpy())
然后，我们定义了主要的训练/验证逻辑：
# Set up training configurations
n_epochs = 300
IC_weight= tf.constant(1.0, dtype=tf.float32)
ODE_weight= tf.constant(1.0, dtype=tf.float32)
loss_tracker = LossTracking()
val_loss_hist = []
# Set up optimizer
optimizer = keras.optimizers.Adam(learning_rate=1e-3)
# Instantiate the PINN model
PI_DeepONet= create_model(mean, var)
PI_DeepONet.compile(optimizer=optimizer)
# Configure callbacks
_callbacks = [keras.callbacks.ReduceLROnPlateau(factor=0.5, patience=30),
tf.keras.callbacks.ModelCheckpoint('NN_model.h5', monitor='val_loss', save_best_only=True)]
callbacks = tf.keras.callbacks.CallbackList(
_callbacks, add_history=False, model=PI_DeepONet)
# Start training process
for epoch in range(1, n_epochs + 1):
print(f"Epoch {epoch}:")
for X_init, X in zip(ini_ds, train_ds):
# Calculate gradients
ODE_loss, IC_loss, total_loss, gradients = train_step(X, X_init,
IC_weight, ODE_weight,
PI_DeepONet)
# Gradient descent
PI_DeepONet.optimizer.apply_gradients(zip(gradients, PI_DeepONet.trainable_variables))
# Loss tracking
loss_tracker.update(total_loss, IC_loss, ODE_loss)
# Loss summary
loss_tracker.history()
loss_tracker.print()
loss_tracker.reset()
####### Validation
val_res = ODE_residual_calculator(X_val[:, :1], X_val[:, 1:-1], X_val[:, -1:], PI_DeepONet)
val_ODE = tf.cast(tf.reduce_mean(tf.square(val_res)), tf.float32)
X_val_ini = X_val[X_val[:, 0]==0]
pred_ini_valid = PI_DeepONet.predict({"forcing": X_val_ini[:, 1:-1], "time": X_val_ini[:, :1]}, batch_size=12800)
val_IC = tf.reduce_mean(keras.losses.mean_squared_error(0, pred_ini_valid))
print(f"val_IC: {val_IC.numpy():.4e}, val_ODE: {val_ODE.numpy():.4e}, lr: {PI_DeepONet.optimizer.lr.numpy():.2e}")
# Callback at the end of epoch
callbacks.on_epoch_end(epoch, logs={'val_loss': val_IC+val_ODE})
val_loss_hist.append(val_IC+val_ODE)
# Re-shuffle dataset
ini_ds = tf.data.Dataset.from_tensor_slices((X_train_ini))
ini_ds = ini_ds.shuffle(5000).batch(ini_batch_size)
train_ds = tf.data.Dataset.from_tensor_slices((X_train))
train_ds = train_ds.shuffle(100000).batch(col_batch_size)
这是一段相当长的代码，但它应该是自解释的，因为我们已经覆盖了所有重要部分。
为了可视化训练性能，我们可以绘制损失收敛曲线：
# History
fig, ax = plt.subplots(1, 3, figsize=(12, 4))
ax[0].plot(range(n_epochs), loss_tracker.loss_history['IC_loss'])
ax[1].plot(range(n_epochs), loss_tracker.loss_history['ODE_loss'])
ax[2].plot(range(n_epochs), val_loss_hist)
ax[0].set_title('IC Loss')
ax[1].set_title('ODE Loss')
ax[2].set_title('Val Loss')
for axs in ax:
axs.set_yscale('log')
训练结果如下所示：
图 12. 损失收敛图。（图像由作者提供）
此外，我们还可以看到在训练过程中某一特定目标 s(·)的预测准确性如何变化：
在训练开始时，我们可以看到模型预测与真实值之间存在明显的差异。然而，到训练结束时，预测的 s(·)已经收敛到真实值。这表明我们的物理信息深度网络学习得很充分。
6. 结果讨论
一旦训练完成，我们可以重新加载保存的权重并评估性能。
在这里，我们随机挑选了三个 u(·)轮廓从测试数据集中，并将其对应的 s(·)与我们的物理信息深度网络预测的结果以及由数值 ODE 求解器计算的真实值进行了比较。我们可以看到，预测结果与真实值几乎无法区分。
图 13. 从测试数据集中随机选择了三个 u(·)轮廓，这些轮廓显示在上排。下排显示了对应的 s(·)轮廓。我们可以看到，物理信息深度网络预测的结果与由数值 ODE 求解器计算的真实值几乎无法区分。（图像由作者提供）
这些结果相当令人惊讶，考虑到我们甚至没有使用任何 s(·)的观测数据（除了初始条件）来训练 DeepONet。这表明控制 ODE 本身为模型提供了足够的“监督”信号，以做出准确的预测。
另一个有趣的评估点是所谓的“分布外”预测能力。由于我们在训练 DeepONet 时强制执行了控制方程，我们可以预期训练得到的物理信息深度网络在 u(·)轮廓超出训练 u(·)分布时仍能做出不错的预测。
为了测试这一点，我们可以使用不同的长度尺度生成 u(·)轮廓。以下结果显示了使用 0.6 长度尺度生成的三个 u(·)轮廓，以及预测的 s(·)。这些结果相当不错，考虑到物理信息深度网络是用 0.4 的长度尺度训练的。
图 14. 训练得到的物理信息深度网络显示了一定程度的分布外预测能力。（图像由作者提供）
然而，如果我们继续将长度尺度减小到 0.2，我们会注意到明显的差异开始出现。这表明训练得到的物理信息深度网络（DeepONet）在泛化能力上存在限制。
图 15. 物理信息深度 ONet 的泛化能力是有限的。（作者提供的图像）
较小的长度尺度通常会导致更复杂的 u(·)轮廓，这与用于训练的 u(·)轮廓可能有很大不同。这可以解释为何训练后的模型在较小长度尺度区域预测准确性遇到挑战。
图 16. 我们训练的模型在泛化到较小长度尺度区域时面临挑战，因为 u(·)轮廓更复杂，与训练数据有较大区别。（作者提供的图像）
总的来说，我们可以说开发的物理信息深度 ONet 能够在仅有 ODE 约束的情况下正确学习系统动态和从输入函数到输出函数的映射。此外，物理信息深度 ONet 在处理“超分布”预测方面显示出一定的能力，这表明训练模型与控制 ODE 对齐可以提高模型的泛化能力。
7. 收获
我们在探索物理信息深度 ONet 的过程中走了很长一段路。从理解深度 ONet 和物理信息学习的基本概念，到通过代码实现看到它们的实际应用，我们已经详细讲解了这一强大方法在求解微分方程中的应用。
这里有几点关键的收获：
1️⃣ 深度 ONet是一个强大的操作符学习框架，得益于其创新的分支和主干网络架构。
2️⃣ 物理信息学习明确地将动态系统的控制微分方程纳入学习过程，从而具有提高模型解释性和泛化能力的潜力。
3️⃣ 物理信息深度 ONet结合了深度 ONet 和物理信息学习的优势，呈现出作为学习功能映射的有前景工具，同时遵循相关的控制方程。
希望你喜欢这次对物理信息深度 ONet的深入探讨。接下来，我们将转向使用物理信息深度 ONet 解决逆问题。敬请关注！
如果你觉得我的内容有用，可以在这里请我喝杯咖啡🤗 非常感谢你的支持！
你可以在这里找到包含完整代码的辅助笔记本 💻
如果你也对使用物理信息深度 ONet 解决逆问题感兴趣，请随时查看我的新博客：使用物理信息深度 ONet 解决逆问题：带代码实现的实用指南
如果你想了解最新的物理知识学习最佳实践，请查看我目前正在进行的设计模式系列：揭示物理知识驱动神经网络的设计模式
参考文献
[1] Lu 等人，DeepONet：基于算子通用近似定理的非线性算子学习，用于识别微分方程。arXiv, 2019。
[2] Wang 等人，学习参数偏微分方程的解算符，基于物理知识的 DeepOnets。arXiv, 2021。
Optical Flow with RAFT: 第一部分
深入探讨光流的深度学习
·
关注 发布于 Towards Data Science ·14 分钟阅读·2023 年 10 月 3 日
–
照片由 Zdeněk Macháček 提供，发布于 Unsplash
在这篇文章中，我们将学习一种旗舰级深度学习方法，这种方法在 2020 年获得了ECCV最佳论文奖，并被引用超过 1000 次。它也是许多顶级模型在KITTI 基准测试中的基础。这个模型叫做RAFT: Recurrent All-Pairs Field Transforms for Optical Flow，可以在PyTorch或GitHub上轻松获取。实现使其非常易于获取，但模型复杂，理解起来可能会令人困惑。在这篇文章中，我们将把 RAFT 分解为其基本组成部分，并详细了解它们。然后，我们将学习如何在 Python 中使用它来估计光流。在第二部分中，我们将探索隐秘的细节并可视化不同的模块，以便深入理解它们的工作原理。
-
介绍
-
RAFT 的基础
-
视觉相似性
-
迭代更新
-
如何使用 RAFT
-
结论
介绍
光流
光流是图像序列中像素的表观运动。为了估计光流，场景中物体的运动必须有相应的亮度位移。这意味着图像中的一个移动的红球在下一张图像中应具有相同的亮度和颜色，这使我们能够确定它在像素上的移动量。图 1 展示了一个光流的例子，其中一个逆时针旋转的天花板风扇被一系列图像捕捉到。
图 1. 图像序列的光流估计。帧 1，帧 2，帧 1 和帧 2 之间计算的光流。来源：作者。
最右边的彩色图像包含了从帧 1 到帧 2 的每个像素的表观运动，它的颜色编码方式使得不同的颜色表示像素运动的不同水平和垂直方向。这是一个密集光流估计的例子。
对稠密光流的估计为每个像素分配一个二维流向量，描述其在时间间隔内的水平和垂直位移。在稀疏光流中，这个向量仅分配给对应于强特征（如角点和边缘）的像素。为了使流向量存在，该像素在时间 t 的亮度必须与时间 t+1 时相同，这被称为 亮度一致性假设。位置(x,y)* 在时间t 的图像强度或亮度由 I(x,y,t) 给出。下面的图 2 展示了已知像素位移的示例，其中 dx 和 dy 是水平和垂直图像位移，dt 是帧之间的时间差。
图 2。像素从时间 t 到 t+dt 的位移。亮度一致性假设意味着该像素在两个帧中具有相同的颜色和强度。来源：作者。
亮度一致性假设意味着在 (x,y,t) 处的像素在 (x+dx, y+dy, t+dy) 处将具有相同的强度。因此：I(x, y, t) = I(x+dx, y+dy, t+dt)。
从亮度一致性假设出发，我们可以通过在*(x, y, t)* 周围展开右侧的 1ˢᵗ阶泰勒近似来推导光流方程[1]。
光流方程的推导。来源：作者。
水平和垂直梯度 Iₓ 和 Iᵧ 可以通过 Sobel 算子 进行近似，时间梯度 Iₜ 是已知的，因为我们有时间 t 和 t+1 的图像。流方程有两个未知数 u 和 v，分别是时间 dt 内的水平和垂直位移。一个方程中的两个未知数使这个问题成为一个 欠定 问题，许多尝试都旨在解决 u 和 v。RAFT 是一种深度学习方法，用于估计 u 和 v，但实际上，它比仅仅基于两帧预测光流要复杂得多。它经过精心设计，以准确估计光流场，下一节我们将深入探讨它的复杂细节。
RAFT 的基础
RAFT 是一个深度神经网络，能够估计给定一对连续图像的密集光流I₁和I₂。它估计一个流位移场*(f¹,* f²*)，将每个像素(u, v)在I₁中映射到I₂中对应的像素(u’, v’)，其中(u’, v’) = (u +* f¹*(u), v +* f²*(v))*。它通过提取特征、寻找其相关性，然后以模拟优化算法的方式迭代更新流。初始流要么初始化为全 0，要么可以使用向前投影的先前流估计来初始化，这被称为温启动。整体架构如下所示。
图 3. RAFT 的架构。修改自源。
请注意，它包含三个主要模块：特征编码器模块、视觉相似性模块和迭代更新模块。RAFT 架构有两个版本，一个大版本有 480 万参数，一个小版本有 100 万参数，在这篇文章中我们将重点关注大版本，但理解小版本在理解大版本后意义不大。
特征提取
RAFT 通过一个包含六个残差块的卷积神经网络（CNN）对两个输入图像进行特征提取，并将每个图像下采样到 1/8 分辨率，具有 D 个特征图。
图 4. RAFT 的编码块。修改自源。
特征编码器网络g在两个图像上使用共享权重进行操作，而上下文编码器网络f仅在I₁上操作，并提取作为流估计的主要参考的特征。除了细微的差异外，两个网络的整体架构几乎相同。上下文网络使用批归一化，而特征网络使用实例归一化，上下文网络提取C = c + h特征图，其中c是上下文特征图的数量，h是将初始化迭代更新模块隐藏状态的隐藏特征图数量。
特征网络 f 和上下文网络 g 的函数映射。来源：作者。
注意：原始论文中经常使用特征图大小 H/8xW/8 的简写符号：HxW。这可能会令人困惑，因此我们将遵循 H’ = H/8 的约定，使特征图大小为 H’xW’。我们也将提及从 I₁ as g¹中提取的特征图张量，I₂亦然。
视觉相似性
相关体积
视觉相似性是一个 4D H’xW’xH’xW’的全对关联体积C，通过计算特征图的点积得到。
4D 相关体积的计算。修改自源。
在相关体积中，来自特征图g¹的每个像素与特征图g²中的每个像素都有一个计算得到的相关性，我们称这些相关性中的每一个为2D 响应映射 (见图 5)。想象在 4D 空间中可能有些挑战，所以可以将体积的前两个维度展平：(H’xW’)xH’xW’，现在我们得到一个 3D 体积，其中g¹的每个像素都有自己的 2D 响应映射，显示其与g²的每个像素位置的相关性。由于特征是从图像中提取的，响应映射实际上指示了I₁的给定像素与I₂的每个像素的相关程度。
视觉相似性是一种全对 Correlation Volume，通过计算每个像素位置处的每个特征图的相关性，将I₁的像素与I₂的每个单一像素联系起来
相关金字塔
相关体积有效地提供了关于小像素位移的信息，但可能难以捕捉较大的位移。为了捕捉大和小的像素位移，需要多个级别的相关性。为此，我们构建了一个包含多个相关体积级别的相关金字塔，其中不同级别的相关体积通过对相关体积的最后两个维度进行平均池化来生成。平均池化操作在体积的最后两个维度产生了I₂的粗略相关特征，这使得I₁的精细特征能够与I₂的逐渐粗略的特征相关联。每个金字塔级别包含越来越小的 2D 响应映射。
图 5。左：I₁中单个像素与I₂所有像素的关系。右：相关金字塔中各种相关体积的 2D 响应映射。来源。
图 5 显示了不同平均池化级别的不同 2D 响应映射。相应的相关体积的尺寸被堆叠到一个 5D 相关金字塔中，其中包含四个级别的核大小：1、2、4 和 8。该金字塔提供了关于大和小位移的强大信息，同时保持对I₁的高分辨率。
相关查找
相关查找运算符 L꜀ 通过在每个级别的相关金字塔中索引特征来生成新的特征图。给定当前的光流估计*(f¹,* f²*)，I₁的每个像素：x = (u, v)映射到其在I₂中估计的对应关系：x’ = (u + f¹(u) + v + f²(v))。我们定义了x’*周围的局部邻域：
以像素x’ = (u’, v’)为中心的半径r的邻域。来源：作者。
对应关系是基于其流估计的I₂像素的新位置
所有金字塔层级上的常量半径意味着更大的上下文将被纳入到较低层级中。 即，半径为 4 对应于原始分辨率下的 256 像素。
实际上，这个邻域是围绕每个细分分辨率像素中心的正方形网格，r = 4 时，我们在每个像素周围得到一个 9x9 的网格，其中每个维度的长度为 (2r + 1)。我们通过 双线性重采样 在网格定义的位置周围对每个像素的相关性特征进行重采样（边缘位置使用零填充）。由于流偏移和平均池化，邻域网格值可能是浮点数，双线性重采样通过对附近像素的 2x2 子邻域进行加权平均来处理这一点。换句话说，重采样将提供 亚像素 精度。我们在金字塔的每一层的所有像素位置进行重采样，这可以通过 PyTorch 的 F.grid_sample() 高效完成。这些重采样后的特征被称为 相关性特征，并输入到更新块中。
高效的相关性查找（可选）
相关性查找的复杂度为 O(N²)，其中 N 是像素数量，这可能会成为大图像的瓶颈，但存在一种等效操作，其复杂度为 O(NM)，其中 M 是金字塔层数。该操作将相关性金字塔与查找相结合，利用了内积和平均池化的线性特性。下图显示了在 2ᵐx2ᵐ 网格上的平均相关性响应 Cᵐ（金字塔层级 m）。
等效的相关性实现。来源。
对于给定的金字塔层级 m，我们不需要对特征图 g¹ 进行求和，这意味着可以通过将特征图 g¹ 与平均池化后的特征图 g² 进行内积来计算相关性，这具有 O(N) 的复杂度。由于这仅适用于单个金字塔层级 m，我们必须为每一层计算这个内积，使其复杂度为 O(M)，总复杂度为 O(NM)。我们不是预计算金字塔的相关性，而是只预计算池化特征图，并在查找发生时按需计算相关性值。
迭代更新
更新操作符估计一系列光流：{f₀, f₁ ,…, fₙ}** 从初始起点 f₀ 开始，该起点可以是全 0 或向前投影的先前光流估计（热启动）。在每次迭代 k 中，它产生一个光流更新方向 Δf，该方向被加到当前估计中：fₖ₊₁ = fₖ + Δfₖ。更新操作符模仿优化算法，并经过训练以提供更新，使得估计的光流序列收敛到一个固定点：fₖ → f*。
更新块
更新块的输入包括：相关特征、当前光流估计、上下文特征和隐藏特征。其结构及突出显示的子块如下所示。
图 6. 大型架构的 RAFT 更新块，不同子块突出显示。蓝色-特征提取块，红色 — 递归更新块，绿色 — 光流头。改编自 来源。
更新块中的子块包括：
-
特征提取块 — 从相关性、光流和 I₁（上下文网络）中提取运动特征。
-
递归更新块 — 递归计算光流更新
-
光流头 — 最终卷积层，将光流估计重新调整为 H/8 x W/8 x 2
如图 6 所示，递归更新块的输入是光流、相关性和上下文特征的连接。潜在的隐藏状态由上下文网络中的隐藏特征初始化。（上下文网络提取了一堆 2D 特征图，然后将其分离为上下文特征图和隐藏特征图）。递归更新块由 2 个可分离的 ConvGRU 组成，这些 ConvGRU 可以在不显著增加网络规模的情况下增加感受野。在每次更新时，递归更新块中的隐藏状态被传递到光流头，以获得尺寸为 H/8 x W/8 x 2 的光流估计。该估计随后使用凸上采样进行上采样。
凸上采样
RAFT 的作者实验了 双线性 和凸上采样，并发现凸上采样提供了显著的性能提升。
图 7. 双线性 VS 凸上采样的比较。 来源。
凸上采样 估计每个细像素为其相邻 3x3 粗像素的 凸组合
让我们分解一下凸上采样的工作原理，下面的图 8 提供了一个很好的视觉示例。
图 8. 单个全分辨率像素（紫色）的凸上采样示例。 来源。
首先，我们假设一个细分分辨率的像素是其最近的粗分辨率邻居的凸组合。这一假设意味着粗分辨率像素的加权和必须等于真实的细分分辨率像素，且权重之和为一且非负。由于我们是按八倍因子上采样的，每个粗分辨率像素必须分解成 64 个(8x8)的细分像素（图 8 中的视觉效果不按比例）。我们还注意到 3x3 网格中心的每个 64 个像素都需要自己的权重集，总共需要的权重数为：(H/8 x W/8 x (8x8x9))。
实际上，权重由神经网络参数化，凸上采样块使用两个卷积层来预测一个(H/8 x W/8 x (8x8x9))的掩码，然后对九个邻居的权重进行 softmax，得到形状为(H/8 x W/8 x (8x8))的掩码。然后我们使用这个掩码来获得邻域的加权组合，并重新调整以得到 HxWx2 的流场。
训练
RAFT 的目标函数能够捕捉所有迭代的流预测。形式上，它是流预测和真实值之间加权的l1距离的总和，权重以指数形式增长。
RAFT 的损失，γ = 0.8。 来源。
如何使用 RAFT
我们可以使用 RAFT 来估计我们自己图像上的密集光流。首先，我们需要克隆GitHub 仓库并下载模型。此教程的代码在GitHub上。
!git clone https://github.com/princeton-vl/RAFT.git
%cd RAFT
!./download_models.sh
%cd ..
预训练的 RAFT 模型有几种不同的版本，根据作者，它们是：
-
raft-chairs — 在 FlyingChairs 上训练
-
raft-things — 在 FlyingChairs + FlyingThings 上训练
-
raft-sintel — 在 FlyingChairs + FlyingThings + Sintel + KITTI 上训练（用于提交的模型）
-
raft-kitti — raft-sintel 在仅 KITTI 上微调
-
raft-small — 在 FlyingChairs + FlyingThings 上训练
接下来，我们将 RAFT 的核心添加到路径中
sys.path.append('RAFT/core')
现在，我们需要一些辅助函数来与 RAFT 类接口。注意：这些辅助函数仅为 CUDA 编写，但你可以通过Colab轻松访问 GPU。
import torch
from raft import RAFT
from utils import flow_viz
from utils.utils import InputPadder
def process_img(img, device='cuda'):
return torch.from_numpy(img).permute(2, 0, 1).float()[None].to(device)
def load_model(weights_path, args):
""" Loads model to CUDA only """
model = RAFT(args)
pretrained_weights = torch.load(weights_path, map_location=torch.device("cpu"))
model = torch.nn.DataParallel(model)
model.load_state_dict(pretrained_weights)
model.to("cuda")
return model
def inference(model, frame1, frame2, device='cuda', pad_mode='sintel',
iters=12, flow_init=None, upsample=True, test_mode=True):
model.eval()
with torch.no_grad():
# preprocess
frame1 = process_img(frame1, device)
frame2 = process_img(frame2, device)
padder = InputPadder(frame1.shape, mode=pad_mode)
frame1, frame2 = padder.pad(frame1, frame2)
# predict flow
if test_mode:
flow_low, flow_up = model(frame1,
frame2,
iters=iters,
flow_init=flow_init,
upsample=upsample,
test_mode=test_mode)
return flow_low, flow_up
else:
flow_iters = model(frame1,
frame2,
iters=iters,
flow_init=flow_init,
upsample=upsample,
test_mode=test_mode)
return flow_iters
def get_viz(flo):
flo = flo[0].permute(1,2,0).cpu().numpy()
return flow_viz.flow_to_image(flo)
注意到*inference()*中的输入填充，我们需要确保所有图像的尺寸都能被 8 整除。raft.py 代码可以从命令行轻松访问，但如果我们想要与之接口，我们需要重写部分代码，或者可以创建一个特殊的类来传递参数。
# class to interface with RAFT
class Args():
def __init__(self, model='', path='', small=False,
mixed_precision=True, alternate_corr=False):
self.model = model
self.path = path
self.small = small
self.mixed_precision = mixed_precision
self.alternate_corr = alternate_corr
""" Sketchy hack to pretend to iterate through the class objects """
def __iter__(self):
return self
def __next__(self):
raise StopIteration
Args 类的默认初始化将直接与任何大型 RAFT 模型进行接口。为了演示 RAFT，我们将使用一个慢速旋转的天花板风扇视频的帧。现在我们可以加载一个模型并估算光流。
model = load_model("RAFT/models/raft-sintel.pth", args=Args())
flow_low, flow_up = inference(model, frame1, frame2, device='cuda', test_mode=True)
测试模式将返回 1/8 分辨率光流以及凸性上采样光流。
图 9. 上：RAFT 的输入图像序列。下：1/8 分辨率和上采样的光流估计。图片来源于作者。来源：作者。
再次，我们可以看到凸性上采样的显著优势，现在让我们来看看额外迭代的优势。图 10 展示了一个风扇图像上的 20 次迭代的 GIF 动画。
flow_iters = inference(model, frame1, frame2, device='cuda', pad_mode=None, iters=20, test_mode=False)
图 10. 光流估计的渐进迭代。来源：作者。
我们可以看到前几次迭代的明显好处，在这种情况下，模型能够在约 10 次迭代中收敛。现在我们将尝试使用温初始化，将前一个 1/8 分辨率的流估计传递给推理函数。
# get previous estimate at 1/8 res
flow_lo, flow_up = inference(model, frame1, frame2, device='cuda', pad_mode=None, iters=20, test_mode=True)
# 0 initialization
flow_lo_cold, flow_up_cold = inference(model, frame2, frame3, device='cuda', pad_mode=None, flow_init=None, iters=20, test_mode=True)
# warm initialization
flow_lo_warm, flow_up_warm = inference(model, frame2, frame3, device='cuda', pad_mode=None, flow_init=flow_lo, iters=20, test_mode=True)
图 11. 在第 2 和第 3 帧之间进行 0 VS 温初始化的光流估计。来源：作者。
在这种情况下，我们并未看到任何改善，右侧的温初始化实际上看起来比初始化为 0 的流还要糟糕。对于这个视频序列来说，温启动的好处似乎微乎其微，但在不同的环境中可能会有用。
结论
在本文中，我们了解了 RAFT，一个能够估算准确流场的先进模型。RAFT 能够通过计算从提取的特征图中的所有像素的全对相关体积来捕捉所有像素之间的关系。建立相关金字塔以捕捉大和小的像素位移。查找运算符基于当前流估计从相关金字塔中提取新的相关特征。更新块使用相关特征和当前流估计提供迭代更新，收敛到最终的流估计，然后使用凸性上采样进行上采样。在第二部分，我们将详细探讨网络并了解一些关键块的工作方式。
参考资料
[1] Horn, B. K. P., & Schunck, B. G. (1981). 确定光流。人工智能, 17(1–3), 185–203. doi.org/10.1016/0004-3702(81)90024-2
[2] Teed, Z., & Deng, J. (2020). Raft: 循环全对场变换用于光流。计算机视觉 — ECCV 2020, 402–419. doi.org/10.1007/978-3-030-58536-5_24
RAFT 中的光流：第二部分
解密深度光流模型
·
关注 发表在Towards Data Science ·10 分钟阅读·2023 年 10 月 3 日
–
图片由Kevin Hansen提供，来源于Unsplash
在这篇文章中，我们将以另一种方式来看待RAFT。第一部分中的直接方法能够详细解析网络的细节，但在这里我们将可视化这些细节并建立有价值的直觉。在第一部分中，我们的目标是理解 RAFT，以便我们可以直接使用它；而在第二部分，我们将旨在以一种方式理解 RAFT，使我们能够利用其架构的不同部分来构建自己的模型。
这里是本文的概览：
-
动机
-
RAFT 架构
-
查找操作符
-
迭代更新
-
结论
动机
RAFT 的概念被许多后续工作利用，理解 RAFT 是理解这些新方法的关键。你如何知道 RAFT 的哪些部分可以或应该被利用？为什么许多后续工作使用相关体？这些问题的答案来自于掌握 RAFT 的内部工作原理，仅仅了解论文表面内容通常是不够的，有时我们需要深入探讨，RAFT 也不例外。
RAFT 架构
首先，快速回顾一下 RAFT，其架构可以分解为三个基本模块，如下所示。
图 1. RAFT 的架构。修改自 来源。
特征提取
特征编码器是一个卷积神经网络（CNN），利用共享权重从图像 I₁ 和 I₂ 中提取特征。上下文编码器提取上下文和隐藏特征，这些特征都被输入到迭代更新模块中。
特征网络 f 和上下文网络 g 的函数映射。来源：作者。
4D 相关体的计算。修改自 来源。
视觉相似性
两个特征图的点积形成一个 4D 全对相关体，其中 g¹ 的每个像素映射到 g² 的所有像素，这些映射被称为 2D 响应图。其中 g¹ 和 g² 分别是从 I₁ 和 I₂ 中提取的特征图张量。 在相关体的最后两个维度上执行平均池化，卷积核大小为：1、2、4、8。
图 2. 左：I₁ 中单个像素与 I₂ 中所有像素的关系。右：相关金字塔中各种相关体的 2D 响应图。来源。
我们将这些相关体积堆叠成一个 5D 相关金字塔，每一层将 g¹ 的精细像素与 g² 的逐渐粗糙的像素特征相关联。这使我们能够捕获关于大规模和小规模像素位移的信息。查找操作符 从相关体积中提取相关特征。它获取 g¹ 的一个精细特征像素及其对应的光流场，并计算其新的明显位置，这称为其 对应关系。然后，它在预定义半径 r 周围形成一个 2D 网格，并随后沿网格执行 亚像素 双线性重采样 以获得新的网格值。这个重采样的特征网格包含流（下降）方向信息，查找操作符对金字塔中每层的每个像素执行此操作。这些逐像素特征网格称为 相关特征，然后被重新塑形并输入到 迭代更新块。
迭代更新
迭代更新块接受四个输入：上下文特征、相关特征、当前流估计和隐藏特征。流和相关特征被一起编码为运动特征，因为它们都描述了特征像素的相对运动。上下文特征保持不变，作为更新块的稳定参考。该块本身由一个 ConvGRU 组成，该网络计算一定数量的递归更新，随后由一个包含卷积层的流头将隐藏状态转换为原始输入分辨率的 1/8 的流估计。
图 3. RAFT 更新块用于大型架构，突出显示了不同的子块。蓝色—特征提取块，红色—递归更新块，绿色—流头。修改自 Source。
更新操作符的功能类似于优化算法，这意味着它从初始流 f₀ 开始，并迭代计算新的流值 Δfₖ，这些值被添加到先前的流估计中：fₖ₊₁ = fₖ + Δfₖ，直到收敛到一个固定值：fₖ → f*。在执行迭代更新时，流估计和相关特征（不是相关金字塔）会持续被精炼。一旦迭代耗尽，流估计将从 1/8 上采样到原始分辨率。
凸上采样
凸上采样 通过将每个精细像素估计为其邻近的 3x3 粗像素网格的 凸组合。权重由一个神经网络参数化，该网络能够为每个精细像素学习最佳权重。下方显示了一个示例。
图 4. 上：双线性 VS 凸采样的比较。下：单个全分辨率像素（紫色）的凸采样示例。 来源。
学习光流偏移量
重要的是要记住，RAFT 不一定估计光流，它估计的是从起始点的 光流偏移量，输出的是这些光流偏移量的累积。第一次估计是对 I₁ 像素位置上 0 的先前光流的更新。I₁ 的信息来自初始隐藏状态和上下文特征，这为更新块提供了持续的反馈，指导学习，使 RAFT 能够估计光流偏移量，而不是重新估计光流。
RAFT 不直接估计光流，它的更新块从起始点估计光流偏移量，模型输出的是光流偏移量的累积。
在训练过程中，网络的递归更新模拟了优化算法的步骤，其中每个新的光流估计 fₖ₊₁ = fₖ + Δfₖ 被目标函数越来越严格地审视，迫使网络随着迭代次数的增加学习更保守的 Δfₖ 估计。目标函数捕捉所有光流更新，是光流预测与真实值之间的加权 l1 距离的总和，权重呈指数增长。
RAFT 的损失，γ = 0.8。 来源。
这种优化算法的模仿以及查找操作符的半径协同作用，限制了每次更新的搜索空间，从而减少了过拟合的风险，加快了训练速度，并且提高了泛化能力。
解读 RAFT
现在我们可以开始解读 RAFT，并深入了解它如何进行预测。这个教程的代码位于 GitHub，在这里 RAFT 代码已被修改，以在每个主要块处输出其中间特征图。测试图像是逆时针旋转的吊扇，这会产生几乎所有方向的光流。为了获得较大的光流位移，我们将跳过序列中的几个图像，这将使光流特征更加明显，便于研究。
图 5. 预测的光流。来源：作者。
现在让我们检查来自特征编码器块的图。
图 6. 来自特征编码器块的特征图（0–127）。来源：作者。
特征编码器生成了fmap 1和fmap 2，它们看起来像噪声，但实际上包含了对流动估计至关重要的信息。隐藏的特征图类似于输入图像I₁，它们直接初始化了 ConvGRU 更新块的隐藏状态，提供了关于I₁中像素的宝贵信息。上下文特征图与输入图像有些相似，突出显示了诸如角落和边缘等强特征。原始论文建议，上下文特征提高了网络准确识别空间运动边界的能力。
相关性
相关金字塔对于计算准确的光流至关重要，因为它能够捕捉从I₁到I₂的像素对应关系。正如我们将看到的，相关体积是 RAFT 的核心，我们将可视化其出色的像素位移捕捉能力。我们将通过检查几个测试像素来接近这一点，并查看它们的二维响应图如何捕捉相对位移。我们只能看到特征，但像素位移的信息仍然会显现。
相关金字塔捕捉了图像序列中像素在多个分辨率级别上的对应关系。
下图展示了一个测试图像的估计流动，并标注了一些测试像素。
图 7. 带注释测试像素的流动图像。来源：作者。
每个像素的水平和垂直流场分量的估计值是：
-
像素 0: (-49.4, -4.3)
-
像素 1: (-5.8, -26.4)
-
像素 2: (23.5, -9.3)
访问相关特征
要访问给定测试像素的相关图，我们将以下函数添加到 corr.py 脚本中，以获得任何金字塔级别下给定测试像素的整数索引。
def get_corr_idx(loc, lvl, w=71, h=40):
""" Obtains index of test pixel location in correlation volume.
loc - test pixel location
lvl - Pyramid level
w - 1/8 of padded horizontal image width
h - 1/8 of padded vertical image height
"""
u = np.clip(np.round(loc[2]/(lvl*8)), 0, (w-1))
v = np.clip(np.round(loc[1]/(lvl*8)), 0, (h-1))
return int(u + w*v)
一旦我们获得了测试像素索引，我们可以获取其二维响应图、重新采样的相关特征和每个金字塔级别的对应关系。
test_pixel_idx = get_corr_idx(test_pixel, lvl=(2**i))
# get the 2D response map
corr_response = corr[test_pixel_idx, 0, :, :].detach().cpu().numpy()
# get the resampled correlation feature grids
resampled_corr = bilinear_sampler(corr, coords_lvl)
resampled_corr_response = resampled_corr[test_pixel_idx, 0, :, :].detach().cpu().numpy()
# get the correspondence
pixel_loc = centroid_lvl[test_pixel_idx, :, :, :].cpu().numpy().squeeze()
可视化相关特征
对于每个像素，图 8-10 中展示了前 15 次更新的第一个金字塔级别的二维响应图 GIF。对应关系由红色方框标记。大的相关值（亮点）指示了I₂中的相对像素位置。注意，随着网络学习流动偏移，对应关系如何围绕高相关值收敛。尽管这些是大位移，但在第一个金字塔级别上的所有对相关性仍能捕捉到它们。
图 8. 金字塔级别 1，像素 0 的二维相关响应图。来源：作者。
图 9. 金字塔级别 1，像素 1 的二维相关响应图。来源：作者。
图 10. 金字塔级别 1，像素 2 的二维相关响应图。来源：作者。
相关金字塔能够捕捉所有层级的相关性，但这并不总是显而易见的。随着金字塔层级的提升，事物开始变得更加抽象，确定 RAFT 实际在做什么变得越来越困难，此外，我们看到的是提取特征的相关性，使得事情更加模糊。
检索下降方向
匹配点通过查找操作符来提出下降方向。相关查找操作符在新的匹配点位置周围放置一个半径为r的网格：x’ = (u + f¹(u) + v + f²(v))，其中*(f¹,* f²*)*是当前的流场估计。围绕 x’的网格用于从相关体中进行双线性重采样。这些重采样网格是输入到更新操作符中的相关特征，以预测下一个流估计。下图显示了第一个金字塔层级 2D 相关响应在像素 1 的情况以及其对应的双线性重采样网格；上排是第一次迭代，流初始化为零；下排是第二次迭代。
图 11. 像素 1 放大后的相关响应和双线性重采样网格。上排：迭代 0，下排：迭代 1。匹配点在(行，列)。来源：作者。
请注意上排中，左侧的相关响应与右侧的重采样网格相同，这是由于零流初始化造成的。我们还注意到关于右上角的双线性重采样网格的一个非常重要的点，最大值直接位于中心的左侧。如果我们向右移动三像素并向上移动一像素，即(3, -1)，那么我们就会落在这个大值上。这是从相关体中检索出的建议下降方向。在迭代更新块中，网络利用这些信息来制定实际的下降方向Δfₖ。
在下排中，我们可以看到匹配点大致从(39, 34)移动到(41.93, 33.3)，这是一个(2.93, -0.7)的位移，显示网络确实利用了建议的下降方向。在右下角的重采样网格中，我们看到最大值位于中心并与匹配点对齐，这表明网络已经有一个接近收敛的流预测。
运动特征
运动特征是相关特征和当前流估计的卷积编码。它们提供了需要由更新块细化的像素流信息。以下展示了每次迭代的一些运动特征图。
图 12. 不同迭代下的运动特征图。是的，我挑选了有趣的特征。来源：作者。
看起来，运动特征与具有大幅度移动的像素对应，不同的特征图似乎对应于不同的像素流，这在图 12 右侧的特征 126 和 127 中表现得尤为明显。它们都以类似的方式收敛到实际流动预测。
结论
在这篇文章中，我们学习了 RAFT 及其内部工作原理。我们看到提取的隐藏特征提供了关于 I₁ 的有用信息，而提取的上下文特征则提供了关于 I₁ 强特征的参考信息。我们已经可视化了相关体如何捕捉到小范围和大范围像素位移的信息。事实证明，RAFT 中的相关概念被许多后续工作所采用，从可视化中建立的直觉加强了这种模式。如果你已经读到这里，恭喜你！你现在对 RAFT 有了比表面层次更深入的理解。
参考文献
[1] Teed, Z., & Deng, J. (2020). Raft: Recurrent all-pairs field transforms for optical flow. 计算机视觉 — ECCV 2020, 402–419. doi.org/10.1007/978-3-030-58536-5_24
优化需求满足：行业方法
·
关注 发表在 Towards Data Science · 8 分钟阅读 · 2023 年 2 月 6 日
–
在电信领域，数据科学有许多应用，涵盖从站点规划到预算和管理会计决策、客户生命周期管理以及市场营销。在这篇文章中，我们将讨论一个这样的用例，即网络站点规划。我们将从两个场景开始，以便您可以理解它们。
场景 1： 最近，电信决策者发现由于未规划的站点规划导致了巨大的开销。网络站点当然有各种成本，包括固定建设成本、可变成本和其他可变开销。公司的管理者发现了导致低效的各种差异。现在，管理者想要重新规划一切，关闭不必要的站点，同时确保剩余的站点完全覆盖所有需求/区域。这就是一个用例。目标是减少基站的总数量，同时仍然覆盖所有区域。
场景 2： 另一个用例是当管理者想要限制用于覆盖最大区域数量的站点数量，这可能是由于财务不足或某个季度或年份的预算限制。因此，在这种情况下，我们希望在选择不超过 p 个基站（上限）的情况下最大化我们的覆盖区域。在这种情况下，可能无法覆盖所有区域。
概述
让我们看看如何使用优化技术来轻松解决这两个用例。需求并不总是指用户或区域的数量，但可以从整体上考虑多个因素，这超出了本文的范围。
场景 1
让我们从第一个用例开始，即在满足所有需求的情况下，尽可能少地选择基站。我们将使用一个简单的例子来演示这些概念；然而，在现实世界中，你可能会有不同类型的数据集，包含成千上万的基站，并且每个基站的邻域或覆盖范围可以通过一定的距离阈值来检查。
让我们从一个例子开始。我们有四个基站，区域被分为四个区域。大面积通常会在地图上分成不同的网格。我们不会深入讨论这一点。例如，基站 #1 可以覆盖区域 1 和 2，而基站 #2 可以覆盖区域 3。基站 #2 可以覆盖区域 1、2 和 3。基站 #3 可以覆盖区域 3 和 4。基站 #4 可以覆盖区域 2 和 4。下图展示了这一点。
图 1：区域与网络塔的连接性（作者图）
在这种情况下，我们想要找到满足所有要求的最少数量的基站。例如，我们可能不需要基站 #4，因为需求已经由基站 #1、#2 和 #3 满足。因此，我们可能要排除基站 #4。这正是我们想要完成的：找到最佳的基站组合，以满足所有需求，同时消除不必要的基站。
方程
让我们用 J 表示位置/基站，用 I 表示需求节点/区域。我们还定义 aij，如果需求节点 i 被位置节点 j 覆盖，则 aij 等于 1。aij 可以使用某个距离阈值来定义。
表 1：场景 1 的样本问题表述和变量符号（作者绘制）
每个基站位置由一个二进制变量 xj 表示，该变量指示该位置是否被选中。如果基站被选中，则 xj=1，否则 xj=0。
目标函数是通过最小化选定基站的总和来进行优化。
方程 1：目标函数——场景 1 中最小化选定基站的总和（作者绘制）
约束条件是确保每个区域至少被一个选定的基站覆盖。因此，我们希望每个需求节点 i（在我们案例中是区域）至少由一个位置 j（在我们案例中是基站）服务。
方程 2：约束条件——场景 1 中每个区域必须由至少一个基站/位置覆盖（作者绘制）
因此，对于每个区域 i，我们查看所有基站，看看这些基站是否覆盖区域 i。每个区域应该至少由一个基站覆盖。请注意，在这个特定模型中，区域可以由多个基站服务。
通过解决这个优化问题，电信公司可以在最小化所需基站数量的同时确定最佳的基站位置，以覆盖整个国家/城市/区域。
Code
from pyomo.environ import *
model = ConcreteModel()
# Define parameters regions and cell sites
model.I = Set(initialize=[1,2,3,4]) # regions/demands
model.J = Set(initialize=[1,2,3,4]) # cell sites
model.a = Param(model.I, model.J, initialize={(1,1):1, (1,2):1, (1,3):0, (1,4):0, (2,1):1,(2,2):1,(2,3):1, (2,4):0, (3,1):0,(3,2):0,(3,3):1,(3,4):1,(4,1):0,(4,2):1,(4,3):0, (4,4):1})
# Define variables
model.x = Var(model.J, within=Binary)
# Define objective function
def objective(model):
return sum(model.x[j] for j in model.J)
model.obj = Objective(rule=objective, sense=minimize)
# Define constraints (sumproduct), for each Region i, we look at all cell sites, and see whether this cell site covers region i or not. sum must be greator than equal to 1
def demand_rule(model, i):
return sum(model.a[i,j]*model.x[j] for j in model.J) >= 1
model.demand_constraint = Constraint(model.I, rule=demand_rule)
# Solve the model
solver = SolverFactory('glpk')
results = solver.solve(model)
# Print the solution
print("Optimal solution:")
for j in model.J:
if model.x[j].value == 1:
print(f"Location {j} is open")
print(f"Number of open locations: {model.obj()}")
图 1：场景 1 的结果（作者绘制）
正如我们所见，我们可以用两个基站覆盖所有需求点/区域，你可以使用上面的图表自行验证结果。
场景 2
现在我们进入第二个用例，讨论如何对可用的站点位置设置上限，以满足最大需求。考虑到可用的基站数量有限，我们希望最大化覆盖范围。例如，公司可能决定设置一个最大为两个基站的上限，并关闭其他基站以节省费用。现在的问题是我们能用两个基站覆盖多少个区域。当然，不同的组合是可能的；例如，如果我们只使用基站 1 和 2，并关闭基站 3 和 4，那么我们实际上会在区域 3 失去覆盖。也许有某些两个基站的组合可以覆盖整个区域，这将通过我们的优化方程来解决。可能存在不止一个可行的解决方案。
方程
表 2：场景 2 的样本问题表述和变量符号（作者绘制）
在这种情况下，我们希望基本上最大化我们的覆盖范围，我们用 y 变量定义覆盖的概念。在这种情况下， yi 表示我们的区域是否会被选中。如果被选中，则 yi 将为 1，否则为 0。我们基本上希望最大化这个总和，也就是说，我们希望每个区域都被覆盖。这可以通过以下方程表示。
方程 3：目标函数——最大化场景 1 的覆盖/区域（图由作者提供）
然而，我们有一些限制条件，因此可能无法覆盖所有区域。让我们从第一个限制条件开始，即覆盖区域的基站数量必须大于等于 yi。例如，如果某个区域被选中，则 yi=1，那么我们基本上希望覆盖该区域的基站数量大于等于 yi。
方程 4：限制条件——每个选定的区域（yi）需要至少由 1 个基站/位置覆盖，用于场景 2（图由作者提供）
方程与第一部分相同，只是将 1 替换为 yi。因此，对于每个区域 I，我们查看它是否被基站覆盖，并将其加起来以确保它大于等于 yi。
除此之外，最重要的限制条件是公司决定将我们限制在最多两个基站。我们基本上希望基站总数小于等于 p。（在我们的情况下是 2）。
方程 5：限制条件——场景 2 的 p 个基站上限（图由作者提供）
现在，让我们将所有方程组合在一起，再次写出完整的内容。
方程 6：结合上述所有内容的完整方程，用于场景 2（图由作者提供）
代码
from pyomo.environ import *
model = ConcreteModel()
# Define parameters regions and cell sites
model.I = Set(initialize=[1,2,3,4]) # regions/demands
model.J = Set(initialize=[1,2,3,4]) # cell sites
model.a = Param(model.I, model.J, initialize={(1,1):1, (1,2):1, (1,3):0, (1,4):0, (2,1):1,(2,2):1,(2,3):1, (2,4):0, (3,1):0,(3,2):0,(3,3):1,(3,4):1,(4,1):0,(4,2):1,(4,3):0, (4,4):1})
model.p = Param(initialize=1)
model.x = Var(model.J, within=Binary)
model.y = Var(model.I, within=Binary)
def objective(model):
return sum(model.y[i] for i in model.I)
model.obj = Objective(rule=objective, sense=maximize)
##constraint
def coverage_rule(model, i):
return sum(model.a[i,j]*model.x[j] for j in model.J) >= model.y[i]
model.coverage_constraint = Constraint(model.I, rule=coverage_rule)
def cellsitelocation_rule(model):
return sum(model.x[j] for j in model.J) <= model.p
model.cellsite_constraint = Constraint(rule=cellsitelocation_rule)
opt = SolverFactory("glpk")
results = opt.solve(model)
print("Optimal solution:")
for j in model.J:
if model.x[j].value == 1:
print(f"Cellsite {j} is on")
print("Optimal solution:")
for i in model.I:
if model.y[i].value == 1:
print(f"Region {i} is covered")
print(f"Number of regions covered: {model.obj()}")
图 2：场景 2 的结果（图由作者提供）
我们仍然可以用有限数量的基站覆盖所有需求/区域。你现在可以尝试 p=1 来查看可能会遗漏哪些区域。重要的是要注意，使用上限时，你可能会遗漏一些区域，但仍能最大化覆盖范围。
结论
这些是我们在本文中探讨的非常简单的问题。然而，当涉及到这样的问题时，这可能是一个非常强大的工具。我们可以有一个更高级的方程，还带有加权部分，现在由你进一步思考。对于这些类型的问题还有许多其他变体，例如最小化从位置到需求点的运输成本。这些类型的问题有许多有趣的变体。
参考文献
[1] optimization.cbe.cornell.edu/index.php?title=Set_covering_problem
[2] www.im.ntu.edu.tw/~lckung/courses/OR15/slides/OR-Sp15_09_IPapplication.pdf
[3] www.pyomo.org/documentation
深度学习中的神经网络优化
原文：
towardsdatascience.com/optimisation-algorithms-neural-networks-101-256e16a88412
如何超越“普通”梯度下降算法改进训练
·发表于 Towards Data Science ·8 分钟阅读·2023 年 11 月 24 日
–
www.flaticon.com/free-icons/neural-network.神经网络图标。神经网络图标由 andinur 创作 — Flaticon.
背景
在我上一篇文章中，我们讨论了如何通过超参数调整来提高神经网络的性能：
如何通过调整超参数来改善神经网络的“学习”和“训练”
towardsdatascience.com
这是一个过程，通过对学习率和隐藏层数量等最佳超参数进行“调整”，以找到对我们的网络性能最优的参数。
不幸的是，对于大型深度神经网络（深度学习），这种调整过程极其缓慢。改进的一种方法是使用比传统“普通”梯度下降方法更快速的优化器。在这篇文章中，我们将深入探讨最流行的优化器和梯度下降的变体，这些可以提升训练速度以及收敛性，并在 PyTorch 中进行比较！
回顾：梯度下降
在深入之前，让我们快速复习一下梯度下降及其背后的理论。
梯度下降的目标是通过减去参数关于损失函数的梯度（偏导数）来更新模型的参数。学习率 α 用于调节此过程，以确保参数的更新在合理范围内进行，避免过度或不足地达到最优值。
梯度下降。方程由作者提供。
-
θ 是模型的参数。
-
J(θ) 是损失函数。
-
∇J(θ) 是损失函数的梯度。∇ 是梯度算子，也称为 nabla。
-
α 是学习率。
我曾写过一篇关于梯度下降及其工作原理的文章，如果你想对它有更多了解，可以参考一下：
解释了为什么梯度下降在数据科学中经常使用，并提供了 C 语言的实现
towardsdatascience.com
动量
梯度下降通常被可视化为一个球体在山坡上滚动。当它到达山谷底部时，它已经收敛到最小值，即最优值。一个持续向下滚动的球体会获得一定的 动量，然而，普通梯度下降是在每次迭代基础上进行的，并不了解之前的更新。
通过在梯度下降更新中包含动量，它为算法提供了关于先前计算的梯度的信息。
从数学上讲，我们得到的是：
动量梯度下降。方程由作者提供。
其中：
-
v_t 是当前的速度。
-
β 是动量系数，值在 0 和 1 之间。这有时也被解释为“摩擦力”。你需要找到最佳的 β 值，但通常 0.9 是一个不错的基准***。***
-
t 是当前的时间步长或迭代次数。
-
v_{t−1}* 是前一步的速度（上一次计算的值）。
其余术语与之前对普通梯度下降的定义相同！
注意，我们利用了先前梯度的信息来‘加速’当前梯度的方向。这提高了收敛速度，并减少了普通梯度下降中可能出现的任何振荡。
动量在 PyTorch 中也很容易实现。
optimizer = torch.optim.SGD([theta], lr=learning_rate, momentum=momentum)
Nesterov 加速梯度
Nesterov 加速梯度 (NAG)，或称为 Nesterov 动量，是对动量算法的轻微修改，通常能导致更好的收敛效果。
NAG 测量相对于损失函数的梯度时略微超前于 θ. 这改善了收敛性，因为 动量 值通常会朝着最优点前进。因此，每次允许算法略微前进一步可以使其更快收敛。
Nesterov 加速梯度下降。方程由作者编写。
其中 ∇J(θ+βv_{t−1}) 是在当前 θ 前稍微一点的损失函数的梯度。
上述方程中的所有术语与之前的优化器相同，因此我不会再次列出所有术语！
Nesterov 加速梯度也可以在 PyTorch 中轻松实现***：***
optimizer = torch.optim.SGD([theta], lr=learning_rate, momentum=momentum, nesterov=True)
AdaGrad
自适应梯度算法（Adagrad） 是一种梯度下降算法，它使用自适应学习率，如果特征/参数更新得更频繁，学习率会变得更小。换句话说，它对更陡峭的梯度比对较浅的梯度衰减更多。
Adagrad。方程由作者编写。
这里：
-
G 是一个对角矩阵，积累了每个参数在时间步长内所有梯度的平方。
-
ϵ 是一个小的平滑项，用于避免当 G 非常小时的除零问题。
-
⊙ 表示逐元素乘法。这是 Hadamard 乘积
上述方程中的其余术语与之前的优化器相同，因此我不会再次列出所有术语！
元素级矩阵乘法的一个例子，假设 A 和 B 都是 2x2：
Hadamard 乘积的一个例子。方程由作者用 LaTeX 编写。
正如我们所看到的，G 的值越大，对参数的更新就越小。它基本上是平方梯度的移动平均。这确保了学习过程变慢，不会超过最优点。
Adagrad 的一个问题是，它有时会使学习率衰减得太多，导致神经网络过早停止学习并陷入停滞。因此，通常不推荐在训练神经网络时使用 Adagrad。
optimizer = torch.optim.Adagrad([theta], lr=learning_rate)
RMSProp
RMSProp（均方根传播） 通过只考虑最近的梯度来解决 Adagrad 过早结束训练的问题。它通过引入另一个超参数 β 来做到这一点，从而缩小对对角矩阵 G 内部值的影响：
RMSProp。方程由作者用 LaTeX 编写。
上述方程中的所有项都与之前优化器的相同，所以我不会再次列出它们！
像其他优化器一样，RMSProp 在 PyTorch 中实现起来很简单：
optimizer = torch.optim.RMSprop([theta], lr=learning_rate, alpha=alpha, eps=eps)
Adam
我们将要看的最终优化器是自适应矩估计，更为人知的是Adam。该算法结合了动量和 RMSProp，因此可以说是两者的最佳结合。不过，它有几个额外的步骤：
Adam 优化器。公式由作者以 LaTeX 呈现。
前两个步骤和最后一步几乎是我们之前展示的动量和 RMSProp 算法。第三步和第四步是修正v_t和G_t的偏差，因为它们在开始时被初始化为 0。
Adam 是自适应学习率算法，类似于 RMSProp，因此使用此优化器时不一定需要调节学习率。
上述方程中的其他项与之前优化器相同，所以我不会再次列出它们！
以下是如何在 PyTorch 中应用 Adam：
optimizer = torch.optim.Adam([theta], lr=learning_rate)
其他优化器
这里有许多其他梯度下降优化器，我们考虑的只是一阶导数，这些被称为雅可比矩阵。还有一种二阶导数，称为赫西矩阵, 但其计算复杂度为O²，而一阶导数的计算复杂度仅为O。
实际上，深度神经网络有数万到数百万行数据，因此赫西梯度下降方法很少使用。大多数情况下，金标准确实是 Adam 或 Nestorov。
还有批量、小批量和随机梯度下降，这些会影响网络的计算速度。我在之前的文章中写过这些算法。
其他一些常用优化器包括：
完整的综合列表可以在这里找到。
性能比较
下面的代码是对我们之前讨论的不同优化器的比较，针对J(θ) = θ²损失函数。最小值在θ = 0:
import torch
import plotly.graph_objects as go
# Function to perform optimisation and log theta
def optimize(optimizer_class, theta_init, lr, iterations, **kwargs):
theta_values = []
theta = torch.tensor([theta_init], requires_grad=True)
optimizer = optimizer_class([theta], lr=lr, **kwargs)
for _ in range(iterations):
optimizer.zero_grad()
loss = theta.pow(2)
loss.backward()
optimizer.step()
theta_values.append(theta.item())
return theta_values
# Initial values
theta_init = 3.0
learning_rate = 0.01
iterations = 1000
# Optimiser configurations
optim_configs = {
"Momentum": {"optimizer_class": torch.optim.SGD, "lr": learning_rate, "momentum": 0.9},
"Nesterov": {"optimizer_class": torch.optim.SGD, "lr": learning_rate, "momentum": 0.9, "nesterov": True},
"Adagrad": {"optimizer_class": torch.optim.Adagrad, "lr": learning_rate},
"Adam": {"optimizer_class": torch.optim.Adam, "lr": learning_rate},
"RMSprop": {"optimizer_class": torch.optim.RMSprop, "lr": learning_rate}
}
# Run optimization for each optimizer and collect theta values
results = {}
for name, config in optim_configs.items():
results[name] = optimize(**config, theta_init=theta_init, iterations=iterations)
# Plot the results
fig = go.Figure()
for optimiser, theta_values in results.items():
fig.add_trace(go.Scatter(x=list(range(iterations)), y=theta_values, mode='lines', name=optimiser))
fig.update_layout(title="Optimiser Performance Comparison",
xaxis_title="Iteration Number",
yaxis_title="Value of Theta",
legend_title="Optimisers",
template="plotly_white",
width=900,
height=600,
font=dict(size=18),
xaxis=dict(tickfont=dict(size=16)),
yaxis=dict(tickfont=dict(size=16)),
title_font_size=24)
fig.show()
优化器比较。由作者使用 Python 绘制的图。
这个图很有趣，有几个关键点要指出：
-
动量和 Nestorov 都超出了 θ. 的最优值。
-
Adagrad 非常慢。这与我们之前讨论的情况一致，即学习率迅速衰减，导致训练过早停止和学习停滞。
-
Adam 和 RMSProp 似乎是最好的，其中 RMSProp 更快地达到最优值。
当然，这只是一个简单的示例，在实际问题中，最佳的优化器可能会有所不同。因此，尝试各种不同的优化器并选择表现最佳的往往是非常值得的。
这段代码可以在我的 GitHub 上找到：
[## Medium-Articles/Neural Networks/optimisers.py at main · egorhowell/Medium-Articles
在我的中等博客/文章中使用的代码。通过创建帐户来贡献开发…
github.com](https://github.com/egorhowell/Medium-Articles/blob/main/Neural%20Networks/optimisers.py?source=post_page-----256e16a88412--------------------------------)
摘要与进一步的思考
在这篇文章中，我们看到了几种加速和提高普通梯度下降性能的方法。这两种方法类型是基于动量的，使用来自先前梯度的信息，以及基于自适应的，依据计算出的梯度调整学习率。在文献中，Adam 优化器通常是最推荐和最常用于研究的优化器。然而，尝试不同的优化器总是值得的，以确定哪种最适合你的模型。
另一个话题！
我有一个免费的新闻通讯，分析数据，我每周分享成为更好的数据科学家的技巧。没有“虚 fluff”或“点击诱饵”，只有来自实践数据科学家的纯粹可操作的见解。
[## 分析数据 | Egor Howell | Substack
如何成为更好的数据科学家。点击阅读《分析数据》，由 Egor Howell 编写，Substack 出版…
与我联系！
参考文献与进一步阅读
-
关于优化神经网络的优秀博客
这是我关于神经网络的其他一些可能感兴趣的博客：
解释神经网络为何能学习（几乎）任何事物和一切
[towardsdatascience.com ## 前向传播与反向传播：神经网络基础
通过手动和代码（使用 PyTorch）解释神经网络如何“训练”和“学习”数据中的模式
[towardsdatascience.com
优化：Python 中的容量限制设施选址问题
原文：
towardsdatascience.com/optimization-capacitated-facility-location-problem-in-python-57c08f259fe0
查找最佳的仓库数量和位置以降低成本并满足需求
·发表于 Towards Data Science ·阅读时间 12 分钟·2023 年 2 月 28 日
–
图片由作者提供。
目录
-
简介
-
问题陈述
-
实现
3.1. 数据集
3.2. 客户、仓库和需求
3.3. 供应和固定成本
3.4. 运输成本
3.5. 优化
-
探索结果
-
结论
1. 引言
设施选址问题（FLPs） 是经典的优化任务。它们的目标是确定仓库或工厂的最佳潜在位置。
仓库可能有或没有容量限制。这将 有容量（CFLP） 与 无容量（UFLP） 问题变体区分开来。
业务目标是找到一组能够最小化成本的仓库位置。原始问题定义由 Balinski (1965) 提出了两个（年度）成本因素之和的最小化：
-
运输成本
-
仓库固定成本
运输成本指的是从仓库位置到客户的费用。仓库固定成本是特定于位置的。它可能包括如租金、税费、电费和维护等费用。
设施选址 是一个众所周知的主题，具有相当丰富的文献。因此，存在许多问题变体以及方法。这篇文章介绍了经典的 CFLP 公式，并分享了一个使用 PuLP 的实际 Python 示例。
2. 问题陈述
CFLP 的目标是确定能够满足客户需求的仓库数量和位置，同时降低固定和运输成本。因此，我们可以将问题表述为以下目标函数的最小化：
在前一个表达式中：
-
N是客户地点的集合。 -
M是候选仓库地点的集合。 -
fⱼ表示仓库j的年固定成本。 -
tᵢⱼ表示从仓库j到客户i的运输成本。 -
xᵢⱼ是从仓库j到客户i的单位数。 -
yⱼ是一个二进制变量yⱼ ∈ {0,1}，表示是否在位置j建立仓库（yⱼ = 1）或不建立（yⱼ = 0）。
现在让我们考虑将约束添加到目标函数中。
由于我们正在建模一个有容量限制的问题，每个设施j可以供应的年最大容量为Cⱼ。因此，交付给客户的单位数xᵢⱼ不能超过此值：
从仓库j到客户i的年交付单位数必须在零到dᵢ之间，其中dᵢ是客户i的年需求：
最后但同样重要的是，我们必须满足客户的需求。在这个示例中，我们规定每个服务客户地点的仓库必须完全满足其需求：
总之，我们可以如下定义问题：
图片由作者提供。
3. 实施
让我们导入所需的库：
-
NumPy、Pandas用于数据处理。 -
math用于特定的数学函数。 -
GeoPandas用于地理空间表示。 -
Matplotlib用于数据可视化。 -
PuLP用于优化。
import numpy as np
import pandas as pd
import geopandas
from math import sin, cos, asin, acos, radians
from pulp import LpProblem, LpMinimize, LpVariable, LpBinary, lpSum, LpStatus, value
import matplotlib.pyplot as plt
plt.style.use('ggplot')
3.1. 数据集
我们在意大利设定优化问题。
起始数据集可以在 simplemaps.com 上获得。我们可以从这里下载输入 csv 文件，并在MIT 许可证下个人和商业用途均可自由使用。
# Load dataframe
df = pd.read_csv(
'./it.csv',
usecols = ['city', 'lat', 'lng', 'population', 'capital', 'admin_name'])
我们感兴趣的是以下列：
-
city: 城镇名称； -
lat: 纬度； -
lng: 经度； -
population: 居民数量； -
capital: 表示城市是否是主要城市或行政中心； -
admin_name: 最高级别行政区域的名称。
3.2. 客户、仓库和需求
在创建客户、设施和需求时，我们假设：
-
客户是输入城市的一部分（30%）。
-
设施只能在行政中心建立。作为起始条件，我们假设我们可以在意大利 80%的主要城市建立仓库。
-
需求是恒定且全年已知的。它等于客户城镇人口的一个部分（2%）加上一个误差项。
RANDOM_STATE = 2 # For reproducibility
FRACTION_CUSTOMERS = 0.3 # Fraction of cities we want to keep as customers
FRACTION_WAREHOUSES = 0.8 # Fraction of cities we want to keep as warehouse locations
FRACTION_DEMAND = 0.02 # Fraction of citizens of a city that may order a product
# List of the 20 regions of Italy
REGION_LIST = [
'Lombardy', 'Veneto', 'Emilia-Romagna', 'Sicilia', 'Campania', 'Piedmont', 'Puglia',
'Lazio', 'Calabria', 'Tuscany', 'Sardegna', 'Marche', 'Friuli-Venezia Giulia', 'Abruzzo',
'Umbria', 'Trentino-Alto Adige', 'Liguria', 'Basilicata', 'Molise', 'Valle d’Aosta']
# Demand is composed of:
# 1\. A fraction of the population
# 2\. An error term of uniform distribution
# Note: demand is approximated to the closest int
# as its physical meaning denies decimals
df['demand'] = np.floor(
FRACTION_DEMAND * df.population + np.random.uniform(-10, 10, size=(df.shape[0],)))
# Create the warehouses dataframe:
# 1\. Filter the 20 regions of Italy
# 2\. Filter capitals as candidate warehouse locations
# 3\. Sample a fraction of the original cities
facility_df = df.\
loc[df.admin_name.isin(REGION_LIST)].\
loc[df.capital.isin(['admin', 'minor'])].\
sample(frac=FRACTION_WAREHOUSES, random_state=RANDOM_STATE, ignore_index=True)
# Create the customers dataframe:
# 1\. Filter the 20 regions of Italy
# 2\. Sample a fraction of the original cities
customer_df = df.\
loc[df.admin_name.isin(REGION_LIST)].\
sample(frac=FRACTION_CUSTOMERS, random_state=RANDOM_STATE, ignore_index=True)
# Customers IDs list
customer_df['customer_id'] = range(1, 1 + customer_df.shape[0])
注意：在在线数据集中，区域名称Valle d'Aosta包含的是排版（弯曲）撇号（U+2019），而不是打字机（直线）撇号（U+0027）。如果复制此代码，请考虑到这一点。
尽管这对于优化任务不是必需的，但我们可能希望在地图上观察我们的地点。geopandas简化了这一任务。可以使用points_from_xy方法轻松创建一个充满地理空间信息的GeoDataFrame：
def add_geocoordinates(df, lat='lat', lng='lng'):
'''
Add column "geometry" with <shapely.geometry.point.Point> objects
built from latitude and longitude values in the input dataframe
Args:
- df: input dataframe
- lat: name of the column containing the latitude (default: lat)
- lng: name of the column containing the longitude (default: lng)
Out:
- df: same dataframe enriched with a geo-coordinate column
'''
assert pd.Series([lat, lng]).isin(df.columns).all(),\
f'Cannot find columns "{lat}" and/or "{lng}" in the input dataframe.'
return geopandas.GeoDataFrame(
df, geometry=geopandas.points_from_xy(df.lng, df.lat))
customer_df = add_geocoordinates(customer_df)
facility_df = add_geocoordinates(facility_df)
我们可以通过geopandas访问意大利的地图，并绘制客户和潜在的仓库位置：
# Prepare country plot
world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
# Extract and plot the shape of Italy
italy = world[world.name == 'Italy']
ax = italy.plot(color='white', edgecolor='black', figsize=(10, 10))
# Plot customers as points
customer_df.\
plot(ax=ax, marker='X', color='red', markersize=30, alpha=0.5, label='Customer')
# Plot potential facility locations as points
facility_df.\
plot(ax=ax, marker='D', color='blue', markersize=30, alpha=0.5, label='Potential warehouse')
# Add legend
plt.legend(facecolor='white', title='Location')
# Add title
plt.title('Customer and potential warehouses')
# Remove ticks from axis
plt.xticks([])
plt.yticks([])
# Show plot
plt.show()
图片由作者提供。
同样，我们可以观察到 20 个意大利地区的平均需求：
# Prepare region dataframe:
# 1\. Filter the 20 regions of Italy
# 2\. Group by region
# 3\. Calculate:
# - Mean regional latitude
# - Mean regional longitude
# - Sum of regional demand
region_df = df.\
loc[df.admin_name.isin(REGION_LIST)].\
groupby(['admin_name']).\
agg({'lat': 'mean', 'lng': 'mean', 'demand': 'sum'}).\
reset_index()
# Add geo-coordinates
region_df = add_geocoordinates(region_df)
# Plot the shape of Italy
ax = italy.plot(color='white', edgecolor='black', figsize=(10, 10))
# Plot region area colored based on demand
region_df.\
plot(ax=ax, column='demand', marker='o', c='demand', cmap='plasma', markersize=2500, alpha=0.6)
# Add region 'center' as red dots
region_df.\
plot(ax=ax, marker='o', c='red', markersize=25, alpha=0.8, label='Customer location')
# Add region name above the center
for i, row in region_df.iterrows():
plt.annotate(
row.admin_name, xy=(row.lng, row.lat+0.2), horizontalalignment='center')
# Add color bar with demand scale
plt.colorbar(ax.get_children()[1], ax=ax, label='Annual Demand', fraction=0.04, pad=0.04)
# Add title
plt.title('Annual demand by region')
# Remove ticks from axis
plt.xticks([])
plt.yticks([])
# Show plot
plt.show()
图片由作者提供。
为了方便以后使用PuLP，让我们将需求数据存储在customer-demand对的字典中：
# Dictionary of cutomer id (id) and demand (value)
demand_dict = { customer : customer_df['demand'][i] for i, customer in enumerate(customer_df['customer_id']) }
3.3. 供应和固定成本
为了建模供应和固定成本，我们假设：
-
每个仓库可以满足的最大年度供应量等于平均区域需求的 3 倍。
-
每个仓库都有一个固定的年成本为 100.000,00 €，与其位置无关。
与需求相同，我们将供应和固定成本存储在字典中：
# Assumptions:
# 1\. Each warehouse has an annual cost of 100.000,00 euros: rent, electricity, ...
# 2\. Each warehouse can meet 3 times the regional average annual demand
COST_PER_WAREHOUSE = 100_000
SUPPLY_FACTOR_PER_WAREHOUSE = 3
SUPPLY_PER_WAREHOUSE = region_df.demand.mean() * SUPPLY_FACTOR_PER_WAREHOUSE
# Warehouses list
facility_df['warehouse_id'] = ['Warehouse ' + str(i) for i in range(1, 1 + facility_df.shape[0])]
# Dictionary of warehouse id (id) and max supply (value)
annual_supply_dict = { warehouse : SUPPLY_PER_WAREHOUSE for warehouse in facility_df['warehouse_id'] }
# Dictionary of warehouse id (id) and fixed costs (value)
annual_cost_dict = { warehouse : COST_PER_WAREHOUSE for warehouse in facility_df['warehouse_id'] }
3.4. 运输成本
运输成本的估算需要：
-
不同位置之间的距离，以及
-
每单位距离的成本函数。
我们可以使用哈弗森公式来近似两个位置之间的距离：
def haversine_distance(lat1, lon1, lat2, lon2):
'''
Calculate distance between two locations given latitude and longitude.
Args:
- lat1: latitude of the first location
- lon1: longitude of the first location
- lat2: latitude of the second location
- lon2: longitude of the second location
Out:
- Distance in Km
Ref:
- https://en.wikipedia.org/wiki/Haversine_formula
'''
return 6371.01 *\
acos(sin(radians(lat1))*sin(radians(lat2)) +\
cos(radians(lat1))*cos(radians(lat2))*cos(radians(lon1)-radians(lon2)))
让我们在两个城市上尝试一下：
-
米兰（纬度：45.4654219，经度：9.18854）
-
贝尔加莫（纬度：45.695000，经度：9.670000）
haversine_distance(45.4654219, 9.1859243, 45.695000, 9.670000)
45.508144765533906
我们得到的距离为 45.5 公里。不幸的是，这个度量与我们在汽车导航系统上看到的距离不一致，因为我们没有考虑路线：
米兰和贝尔加莫之间的哈弗森距离和路线距离。图片由作者提供。
尽管如此，我们可以使用我们的估算作为任务的合理近似。
最后，我们需要将距离转换为成本度量。目前在写作时，意大利的平均汽油价格为 1.87 €/L（来源）。一辆 EURO VI 卡车的平均耗油量约为 0.38 L/Km（来源）。通过一个简单但合理的近似，我们可以估算在意大利土地上每公里的平均成本为 0.71 €：
def traveling_cost(distance_in_km):
'''
Return traveling cost in euros given a distance in Km.
Args:
- distance_in_km: travel distance in Km
Out:
- cost of the trip in euros
'''
return 0.71 * distance_in_km
现在我们可以计算每个仓库-客户对的旅行成本，并将其存储在字典中：
# Dict to store the distances between all warehouses and customers
transport_costs_dict = {}
# For each warehouse location
for i in range(0, facility_df.shape[0]):
# Dict to store the distances between the i-th warehouse and all customers
warehouse_transport_costs_dict = {}
# For each customer location
for j in range(0, customer_df.shape[0]):
# Distance in Km between warehouse i and customer j
d = 0 if facility_df.city[i]==customer_df.city[j] else haversine_distance(
facility_df.lat[i], facility_df.lng[i], customer_df.lat[j], customer_df.lng[j])
# Update costs for warehouse i
warehouse_transport_costs_dict.update({customer_df.customer_id[j]: traveling_cost(d)})
# Final dictionary with all costs for all warehouses
transport_costs_dict.update({facility_df.warehouse_id[i]: warehouse_transport_costs_dict})
3.5. 优化
让我们回顾一下优化问题：
图片由作者提供。
我们可以定义两个决策变量xᵢⱼ和yⱼ，目标函数和约束条件如下：
# Define linear problem
lp_problem = LpProblem('CFLP', LpMinimize)
# Variable: y_j (constraint: it is binary)
created_facility = LpVariable.dicts(
'Create_facility', facility_df['warehouse_id'], 0, 1, LpBinary)
# Variable: x_ij
served_customer = LpVariable.dicts(
'Link', [(i,j) for i in customer_df['customer_id'] for j in facility_df['warehouse_id']], 0)
# Objective function
objective = lpSum(annual_cost_dict[j]*created_facility[j] for j in facility_df['warehouse_id']) +\
lpSum(transport_costs_dict[j][i]*served_customer[(i,j)] \
for j in facility_df['warehouse_id'] for i in customer_df['customer_id'])
lp_problem += objective
# Costraint: the demand must be met
for i in customer_df['customer_id']:
lp_problem += lpSum(served_customer[(i,j)] for j in facility_df['warehouse_id']) == demand_dict[i]
# Constraint: a warehouse cannot deliver more than its capacity limit
for j in facility_df['warehouse_id']:
lp_problem += lpSum(served_customer[(i,j)] for i in customer_df['customer_id']) <= annual_supply_dict[j] * created_facility[j]
# Constraint: a warehouse cannot give a customer more than its demand
for i in customer_df['customer_id']:
for j in facility_df['warehouse_id']:
lp_problem += served_customer[(i,j)] <= demand_dict[i] * created_facility[j]
我们可以解决优化问题：
lp_problem.solve()
我们可以如下检查结果：
print('Solution: ', LpStatus[lp_problem.status])
Solution: Optimal
我们现在对探索决策变量感兴趣：我们需要多少个仓库？它们的位置在哪里？
4. 探索结果
首先，让我们考虑商业目标：最小化成本。我们可以检查目标函数的值：
value(lp_problem.objective)
8964323.323646087
这是在给定约束条件下我们可以实现的最低成本。任何其他仓库数量或位置的选择都会导致目标函数值更高。
我们可以通过varValue属性访问决策变量。例如，我们可以查看yⱼ在j = 仓库 1时的值：
created_facility['Warehouse 1'].varValue
1.0
由于yⱼ = 1，我们应在该位置建立一个仓库。我们可以轻松地操作变量并计算所需设施的数量：
# List of the values assumed by the binary variable created_facility
facility_values = [i.varValue for i in created_facility.values()]
# Count of each distinct value of the list
[[i, facility_values.count(i)] for i in set(facility_values)]
[[0.0, 59], [1.0, 32]]
只需建造最初预算的 91 个地点中的 32 个即可。35.1%（32 / 91）的潜在仓库足以满足在给定约束条件下的需求。
我们可以将决策变量保存到初始数据框中，并观察所选择的位置：
# Create dataframe column to store whether to build the warehouse or not
facility_df['build_warehouse'] = ''
# Assign Yes/No to the dataframe column based on the optimization binary variable
for i in facility_df['warehouse_id']:
if created_facility[i].varValue == 1:
print('Build site at: ', i)
facility_df.loc[facility_df['warehouse_id'] == i, 'build_warehouse'] = 'Yes'
else:
facility_df.loc[facility_df['warehouse_id'] == i, 'build_warehouse'] = 'No'
Build site at: Warehouse 1
Build site at: Warehouse 2
Build site at: Warehouse 3
Build site at: Warehouse 4
Build site at: Warehouse 7
Build site at: Warehouse 8
Build site at: Warehouse 16
Build site at: Warehouse 18
Build site at: Warehouse 20
Build site at: Warehouse 21
Build site at: Warehouse 22
Build site at: Warehouse 23
Build site at: Warehouse 25
Build site at: Warehouse 26
Build site at: Warehouse 27
Build site at: Warehouse 29
Build site at: Warehouse 33
Build site at: Warehouse 35
Build site at: Warehouse 38
Build site at: Warehouse 48
Build site at: Warehouse 49
Build site at: Warehouse 55
Build site at: Warehouse 56
Build site at: Warehouse 57
Build site at: Warehouse 58
Build site at: Warehouse 63
Build site at: Warehouse 66
Build site at: Warehouse 70
Build site at: Warehouse 74
Build site at: Warehouse 82
Build site at: Warehouse 83
Build site at: Warehouse 84
colors = ['#990000', '#0059b3']
facility_df.build_warehouse.value_counts().plot.barh(
title='Warehouse sites to be established', xlabel='Number of sites', color=colors, ylabel='Establish', figsize=(7,6))
for i, v in enumerate(facility_df.build_warehouse.value_counts()):
plt.text(v, i, ' '+str(round(v,3)), color=colors[i], va='center', fontweight='bold')
图片来源于作者。
# Plot the shape of Italy
ax = italy.plot(color='white', edgecolor='black', figsize=(10, 10))
# Plot sites to establish
facility_df.\
loc[facility_df.build_warehouse =='Yes'].\
plot(ax=ax, marker='o', c='#0059b3', markersize=50, label='Build')
# Plot sites to discard
facility_df.\
loc[facility_df.build_warehouse =='No'].\
plot(ax=ax, marker='X', c='#990000', markersize=40, label='Discard')
# Add title
plt.title('Optimized Warehouse Sites')
# Add legend
plt.legend(title='Warehouse Site', facecolor='white')
# Remove ticks from axis
plt.xticks([])
plt.yticks([])
# Show plot
plt.show()
图片来源于作者。
同样，我们可以遍历决策变量xᵢⱼ，找到优化解中每个仓库服务的客户：
def get_linked_customers(input_warehouse):
'''
Find customer ids that are served by the input warehouse.
Args:
- input_warehouse: string (example: <Warehouse 21>)
Out:
- List of customers ids connected to the warehouse
'''
# Initialize empty list
linked_customers = []
# Iterate through the xij decision variable
for (k, v) in served_customer.items():
# Filter the input warehouse and positive variable values
if k[1]==input_warehouse and v.varValue>0:
# Customer is served by the input warehouse
linked_customers.append(k[0])
return linked_customers
# Warehouses to establish
establish = facility_df.loc[facility_df.build_warehouse =='Yes']
# Plot the shape of Italy
ax = italy.plot(color='white', edgecolor='black', figsize=(30, 30))
# Plot sites to establish
establish.\
plot(ax=ax, marker='o', c='#0059b3', markersize=100, label='Warehouse')
# Plot customers
customer_df.\
plot(ax=ax, marker='X', color='#990000', markersize=80, alpha=0.8, label='Customer')
# For each warehouse to build
for w in establish.warehouse_id:
# Extract list of customers served by the warehouse
linked_customers = get_linked_customers(w)
# For each served customer
for c in linked_customers:
# Plot connection between warehouse and the served customer
ax.plot(
[establish.loc[establish.warehouse_id==w].lng, customer_df.loc[customer_df.customer_id==c].lng],
[establish.loc[establish.warehouse_id==w].lat, customer_df.loc[customer_df.customer_id==c].lat],
linewidth=0.8, linestyle='--', color='#0059b3')
# Add title
plt.title('Optimized Customers Supply', fontsize = 35)
# Add legend
plt.legend(facecolor='white', fontsize=30)
# Remove ticks from axis
plt.xticks([])
plt.yticks([])
# Show plot
plt.show()
图片来源于作者。
5. 结论
在这篇文章中，我们介绍了一个经典的优化挑战：容量受限设施选址问题（CFLP）。我们描述了它的推导过程，并分享了一个实用的 Python 示例。特别地，由于我们从一个原始的地理位置数据集开始，我们覆盖了所有框架问题和求解问题所需的必要步骤和假设。
扫一扫
720
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
