作为深度学习方法与物理模拟更紧密、更通用结合的下一步,我们将把可微分数值模拟纳入到学习过程中。在下文中,我们将 "物理系统的可微分数值模拟 "简称为 "可微分物理"(DP)。

这些方法的核心目标是利用现有的数值求解器,并为其配备根据输入计算梯度的功能。一旦模拟的所有算子都实现了这一功能,我们就可以利用 DL 框架的反向传播自动微分功能,实现让梯度信息从模拟器流入 NN,反之亦然。这样做有很多好处,如改进学习反馈和泛化,我们将在下文概述。

与集成物理约束的损失函数方法相比, 它还能够处理更为复杂的解空间, 而不仅仅是一个反问题。例如, 在前一章中我们使用深度学习来解决单一的反问题, 而利用可微物理学, 我们可以训练神经网络以非常高效地解决更大类别的反问题。

7.1 可微算子

通过DP方法,我们在现有的数值求解器基础上进行建模。也就是说,这种方法在很大程度上依赖于计算方法领域为我们世界中各种物理效应开发的算法。首先,我们需要一个连续的数学模型来描述我们想要模拟的物理效应,如果没有这个模型,我们就会遇到麻烦。但幸运的是,我们可以利用现有的模型方程集合和离散化连续模型的已建立方法。

假设我们有一个关于感兴趣的物理量 的连续表达式,其中 是模型参数(例如扩散系数、粘度或电导率)。 的分量将用编号的下标表示,即

通常,我们对这样一个系统的时间演化感兴趣。离散化得到一个公式 ,我们重新排列以计算时间步长 后的未来状态。在 时刻,通过一系列操作 来计算,使得 ,其中 表示函数分解,即

注:为了将这个求解器整合到深度学习过程中,我们需要确保每个算子相对于其输入提供梯度,即在上述示例中为

请注意,我们通常不需要的所有参数的导数,例如,我们在下文中省略了,假设这是一个给定的模型参数,NN不应该与之交互。当然,它可以在我们感兴趣的解流形中变化,但 不会成为 NN 表示的输出。如果是这种情况,我们就可以在求解器中省略提供 。然而,下面的学习过程自然会将 作为一个自由度。

请注意,通常我们不需要计算的所有参数的导数,例如,在接下来的内容中我们省略了,假设这是一个给定的模型参数,神经网络不应该与之交互。当然,它可以在我们感兴趣的解空间内变化,但不会是神经网络表示的输出。如果是这种情况,我们可以在求解器中省略提供。然而,下面的学习过程自然地可以包含作为一个自由度。

7.2 雅可比矩阵

由于 通常是一个向量值函数,所以 表示雅各布矩阵 ,而不是单一值:

\begin{aligned} \frac{ \partial \mathcal P_i }{ \partial \mathbf{u} } = \begin{bmatrix} \partial \mathcal P_{i,1} / \partial u_{1} & \ \cdots \ & \partial \mathcal P_{i,1} / \partial u_{d} \ \vdots & \ & \ \ \partial \mathcal P_{i,d} / \partial u_{1} & \ \cdots \ & \partial \mathcal P_{i,d} / \partial u_{d} \end{bmatrix} \end{aligned}

如上所述, 表示向量 中的分量数。由于 的一个值映射到另一个值,因此在这里雅可比矩阵是方阵。当然,这并不一定是一般模型方程的情况,但对于可微分模拟来说,非方阵的雅可比矩阵并不会造成任何问题。

在实践中,我们依赖现代深度学习框架提供的_反向模式_微分,关注计算雅可比矩阵的转置与一个向量 的矩阵向量积,即表达式: 。如果我们需要在训练过程中构建和存储遇到的所有完整雅可比矩阵, 将导致巨大的内存开销并不必要地减慢训练速度。相反, 对于反向传播, 我们可以提供更快的运算来计算与雅可比转置的乘积, 因为链式规则的末端总是一个标量的损失函数。

考虑到上面的表述,我们需要通过链式规则来求解 在某个当前状态 下的函数组成链的导数。例如,对于其中的两个:

这只是“经典”链式法则 的向量版本,并直接扩展到更多级联函数的情况 ()。

这里,的导数仍然是雅可比矩阵,但需要注意的是,在链的“末端”,我们有标量损失(参见Overview),最右边的雅可比矩阵将始终是一个具有1列的矩阵,即一个向量。在反向模式下,我们从该向量开始,逐个计算与左侧雅可比矩阵的乘积。

前向和反向传播模式微分的详细内容, 请参阅诸如 Baydin等人的综述之类的外部材料。

7.3 基于 DP 算子的学习

因此,一旦我们模拟器的算子支持雅各布向量乘积的计算,我们就可以将它们集成到 DL 管道中,就像集成普通全连接层或 ReLU 激活一样。

此时,会出现一个非常合理的问题:“大多数物理求解器可以分解为一系列向量和矩阵操作。所有最先进的深度学习框架都支持这些操作,那么我们为什么不直接使用这些操作来实现我们的物理求解器呢?”

确实在理论上这是可能的。但问题是,TensorFlow和PyTorch中的每个向量和矩阵操作都是单独计算的,并且内部需要存储当前的前向计算状态以进行反向传播(上述的“”)。然而对于一个典型的模拟来讲,我们并不太关注求解器产生的每个中间结果。我们通常更关注诸如从 等重要更新步骤。

因此,在实践中,将求解过程分解为一系列有意义但单一的操作符是一个非常好的主意。这不仅通过防止计算不必要的中间结果节省了大量的工作,还允许我们选择计算这些操作符的更新(和导数)的最佳可能的数字方法。例如,由于这个过程与伴随法优化非常相似,我们可以重用在这个领域开发的许多技术,或者利用已有的数字方法。例如,我们可以利用多网格求解器的运行时来进行矩阵求逆。

这种方法的不利之处是,它需要对手头的问题和数值方法有一定的了解。此外,给定的求解器可能无法直接提供梯度计算。因此,如果我们想将深度学习应用于我们没有充分理解的模型方程,通过可微物理(Differentiable Physics,DP)方法进行学习可能不是一个好主意。然而,如果我们真的不理解我们的模型,我们最好还是回去再深入研究一下...

在实践中,我们应该对导数运算符采取贪婪的策略,并且只提供与学习任务相关的那些。例如,如果我们的网络在上述例子中从未产生参数,并且在我们的损失公式中也没有出现导数,那么我们在反向传播步骤中就永远不会遇到导数。

下图总结了基于 DP 的学习方法,并说明了一次 PDE 求解中通常要处理的运算序列。由于许多运算在实践中是非线性的,这通常会使神经网络面临具有挑战性的学习任务:

7.4 一个实际例子

举个简单的例子,我们考虑将被动标量密度在速度场中的对流作为物理模型

我们不直接将其作为残差方程使用 (如4 物理损失中的 v2), 而是用我们喜欢的网格和离散化方案对其进行离散化, 以获得一个随时间更新系统状态的表达式。这是一个标准的_前向_求解过程。为简化起见, 我们这里假设 仅仅是空间中的函数, 即随时间保持不变。我们稍后会讨论 的时间演化。

我们将这个重新表达式表示为 。它将 的状态映射到演化后的新状态,即:

作为逆问题和学习任务的简单示例,让我们考虑找到速度场的问题。这个速度应该将给定的初始标量密度状态在时间转变为由演化到稍后“结束”时间的状态,具有某种形状或配置。通俗地说,我们希望找到一个流场,通过PDE模型将变形为目标状态。表达这一目标的最简单方法是通过两种状态之间的L2损失。因此,我们想要最小化损失函数

作为逆问题和学习任务的一个简单例子,我们来考虑寻找速度场 的问题。这个速度应该能将时间 时的给定初始标量密度状态 转变为一个状态,这个状态通过 演化到之后的 "结束 "时间 ,并具有一定的形状或配置 。非正式地讲,我们希望找到一种流,通过 PDE 模型将 变形为目标状态。表达这一目标的最简单方法是两个状态之间的 损失。因此,我们希望最小化损失函数

请注意,这里描述的这个反问题是一个纯粹的优化任务:不涉及 NN,我们的目标是获得 。我们并不希望将这一速度应用到其他未见过的_测试数据_上,这在实际学习任务中是很常见的。

我们的标记密度 的最终状态完全由通过 开始的演化决定,这就产生了下面的最小化问题:

现在,我们想通过梯度下降(GD)找到这个目标的最小值,梯度由本章前面介绍的可微分物理方法决定。一旦使用梯度下降法,我们就可以比较容易地切换到更好的优化器,或将 NN 引入其中,因此它始终是一个很好的起点。为了方便阅读,我们在下文中将省略雅克比矩阵的转置。不幸的是, 雅克比因子是这样定义的,但我们实际上并不需要未转置的雅克比因子。请记住,实际上我们处理的是转置的雅克比矩阵 ,它被 "缩写 "了。

由于离散化的速度场 包含了所有的自由度,我们所需要做的就是通过 来更新速度,它被分解为 分量部分通常很简单:我们会得到

如果将 表示为一个向量,例如网格的每个单元都有一个条目,那么 同样也将是一个大小相当的列向量。这是因为 始终是一个标量损失函数,因此雅克比矩阵在 维度上的维数为 1。直观地说,这个向量将简单地包含结束时间的 与目标密度 之间的差值。

本身的演化是由我们的离散物理模型 给出的,我们可以互换使用 。因此,更有趣的部分是雅可比矩阵 来计算完整的 。幸运的是,我们不需要 作为一个完整的矩阵,而只需要它与的乘积。

那么,的实际雅可比矩阵是什么呢?为了计算它,我们首先需要完善偏微分方程模型 ,这样可以得到一个可以求导的表达式。在下一章节中,我们将选择一个特定的对流格式和离散化方法,以便我们可以更具体地讨论。

7.4.1 引入一个具体的对流格式

下面我们将在 1 D 笛卡尔网格上使用一个简单的一阶迎风格式,其中单元 具有标量密度 和速度 。为简洁起见, 我们省略时间 处量的标注 ,即下面 简写为 。从上面可以看到, 我们使用_物理模型_更新标量密度 ,其中:

因此, 对于负的 ,我们使用 沿速度相反的方向查看, 即运动学意义上的_后向_。在这种情况下 将为零。对于正的 则反之, 我们将得到一个零值的 ,并通过 得到一个后向差分模板。为了选择前一种情况, 对于负的 我们得到:

因此 给出了 。直观地看,速度的变化取决于密度的空间导数。由于采用了一阶迎风格式,我们只考虑了两个相邻点(高阶方法会依赖的更多项)。

实际上,这一步等同于计算转置矩阵乘法。如果我们把上面的计算改写成 ,那么 。然而,在许多实际情况下,这种乘法的无矩阵实现可能比实际构造 更可取。

对于对流格式,我们可以考虑的另一个导数是与之前密度状态相关的导数,即 ,简写为 。对于单元 ,从上面的得出。然而,为了得到完整的梯度,我们需要加上来自 单元的潜在贡献,这取决于它们速度的符号。这种导数将在下一节中发挥作用。

7.4.2 时间演化

到目前为止,我们只处理了从时间 的一个 更新步,但我们当然可以有任意数量的这样的步骤。毕竟,我们在上文提出的目标是将初始标记状态 提前到时间 时的目标状态,这可能涵盖很长的时间间隔。

在上面的表达式中,每个都依次取决于时的速度和密度状态,即。因此,我们必须追溯损失 的影响,一直追溯到 如何影响初始标记状态。这可能需要通过 对对流格式进行大量评估。

这初听起来很有挑战性:例如,我们可以尝试将时间为的方程(1)插入时间为的方程,然后递归地重复这一过程,直到我们得到一个将与目标相关联的单一表达式。然而,由于雅可比矩阵的线性性质,我们将每个对流步骤,即我们的 PDE 的每次调用视为单独的模块化操作。每一次调用都遵循上一节描述的程序。

给定上述机制, 反向追踪非常简单: 对于 中的每个对流步骤, 我们计算雅可比矩阵与来自损失 或前一个对流步骤的导数的_输入_向量的乘积。我们重复这个过程直到追踪链条从具有 的损失一直回溯到 。从理论上讲, 速度 可以像 一样是时间的函数, 在这种情况下, 我们在每个时间步 都会得到一个梯度 。然而, 为了简化下面的内容, 我们假设我们有一个随时间不变的场, 即我们对每个通过 的对流都重用相同的速度 。现在, 每个时间步骤都会为我们提供一个对 的贡献,我们对所有步骤进行累积。

这里最后一项包含了从时间 到目标时间的完整标量密度回溯。乍看之下, 和式中的各项很令人困惑, 但仔细看, 左边每一行都只是添加了一个额外的时间步的雅可比矩阵。这遵循上面两个算子情况下的链式法则。因此各项包含许多相似的雅可比矩阵, 在实践中可以通过反向遍历 PDE 前向求解生成的计算步骤序列来有效计算。(如上所述, 这里我们省略了雅可比矩阵的转置。)

这种结构也清楚地表明,这一过程与 NN 的常规训练过程非常相似:通过嵌套函数调用对这些雅可比向量积进行评估,这正是深度学习框架训练 NN 的过程(我们只是用权重 代替了速度场)。因此,我们在实践中需要做的就是为 提供一个雅可比向量积的自定义函数。

7.5 隐式梯度计算

作为一个稍微复杂一点的例子,让我们考虑一下泊松方程 ,其中 是感兴趣的量,而 是给定的值。这是一个非常基本的椭圆多项式方程,对从静电到引力场等各种物理问题都很重要。它也出现在流体中,其中 是流体中的标量压力场,而右边的 则由流体速度 的散度给出。

对于流体,我们通常有 ,其中 。这里, 表示新的、无散度的速度场。这一步通常对执行硬约束 至关重要,也被称为 Chorin Projection,或 Helmholtz decomposition。它是向量微积分基本定理的直接结果。

如果我们现在在求解器中引入一个可以修改 的 NN,我们不可避免地需要通过泊松求解进行反向传播。也就是说,我们需要的梯度,在这个符号中,梯度的形式是.

结合起来,我们的目标是计算 。外部梯度(来自 )和内部散度()都是线性算子,它们的梯度很容易计算。主要的困难在于如何从泊松方程中获得逆矩阵 (我们在此将其简化一些,但它通常与时间相关,而且是非线性的)。

在实践中,对于的矩阵向量乘积并不是通过矩阵运算显式计算得到,而是用一个(可能是无矩阵的)迭代求解器来近似计算。例如,共轭梯度(CG)方法在这里非常受欢迎。因此,我们可以理论上将这个迭代求解器视为一个函数,其中。值得注意的是,矩阵求逆是一个非线性的过程,尽管矩阵本身是线性的。由于像CG这样的求解器也是基于矩阵和向量操作的,我们可以将分解为所有求解器迭代过程中一系列简单操作的序列,如,并通过每个操作进行反向传播。这当然是可行的,但不是一个好主意:它可能会引入数值问题,并且会非常慢。如上所述,DL框架默认存储像这个例子中的每一个可微算子(如本例中的)内部状态,因此我们需要组织和保留大量的中间状态在内存中。尽管这些状态对于我们原始的PDE来说完全没有意义,它们只是CG求解器的中间状态。

如果我们退后一步来看,它的梯度就是。在这种情况下,是一个对称矩阵,因此。这是我们在上面的原始方程中遇到的相同逆矩阵,因此我们重复使用未修改的迭代求解器来计算梯度。我们不需要拆开它并通过存储中间状态来减慢速度。然而,迭代求解器计算的矩阵向量乘积。那么在反向传播过程中是什么?在优化设置中,我们总是在正向传播链的末尾有损失函数。反向传播步骤将为输出给出一个梯度,假设这里是,需要传播到正向传播过程的较早操作。因此,我们只需在反向传播过程中调用我们的迭代求解器来计算。假设我们已经选择了一个良好的求解器用于正向传播,那么在反向传播中我们将得到完全相同的性能和精度。

如果您对代码示例感兴趣,phiflow的differentiate-pressure example正是使用这个过程进行通过压力投影步骤进行优化的:一个在右侧受限的流场,通过左侧的内容进行优化,以便在压力投影步骤后与右侧的目标匹配。

这里的主要启示是:在前向计算中,重要的是_不要盲目地进行反向传播_,而是要考虑在前向传递的解析方程中,为哪几步计算梯度。在上述情况下,我们通常可以找到梯度的改进解析表达式,然后对其进行数值逼近。

隐式函数定理和时间 IFT:上述过程本质上产生隐式导数。我们没有显式推导所有前向步骤,而是依靠隐函数定理来计算导数。 Time:我们实际上可以将迭代求解器的步骤视为虚拟的“时间”,并通过这些步骤进行反向传播。与其他Dp方法一样,这使得神经网络能够与迭代求解器交互。一个例子是从 [UBH+20] 中学习 CG 求解器的初始猜测。详细信息和代码可以在这里找到。

7.6 到目前为止可微物理模拟的总结

总结一下, 使用可微物理模拟为我们提供了一种工具, 可以将选择的离散化的物理方程引入深度学习中。与前一章的残差约束相比, 这使得神经网络可以无缝地与物理求解器进行交互。

以前我们会在神经网络训练结束后完全放弃物理模型和求解器: 在5 PINN优化Burgers方程 的例子中,神经网络直接给出解,绕过任何求解器或模型方程。DP 方法与5 PINN优化Burgers方程 中的基于物理的神经网络 (v 2)有本质区别,它与受控离散化 (v 1)更加接近。后者实际上是 DP 训练的一个子集或部分应用。

然而, 与这两种残差方法相反, DP 使得神经网络可以与数值求解器一起训练,因此我们可以在之后的推理时利用物理模型 (由求解器表示)。这使我们能够解决单一的反问题, 并使神经网络对新的输入具有非常强的泛化能力。让我们在 DP 的背景下重新审视5 PINN优化Burgers方程中的示例问题。

  1. 我们不再仅仅预测残差,而是训练一个神经网络来预测速度场。
  2. 该速度场不再直接用于计算损失,而是作为输入馈送到物理求解器中。
  3. 求解器使用预测的速度场进一步演化标量密度场。
  4. 最终的标量密度场与目标进行比较,计算损失。
  5. 梯度通过求解器反向传播, 更新神经网络的参数。

这样, 训练好的网络就学会了预测一个速度场, 这场速度场能够通过物理方程的演化得到目标状态。在部署阶段, 我们可以继续使用物理求解器, 只是用训练好的网络来替代速度场的生成。这种方法增强了对物理约束的遵循性, 并提高了泛化能力。