欢迎来到《基于物理的深度学习书籍》（v0.2）！ 👋

简要说明：本文档包含在物理模拟背景下与深度学习相关的实用而全面的介绍。所有主题都尽可能以 Jupyter 笔记本的形式附带实践代码示例，以便快速入门。除了标准的数据监督学习外，我们还将研究物理损失约束、与可微分模拟更紧密耦合的学习算法、为物理问题量身定制的训练算法，以及强化学习和不确定性建模。我们生活在一个激动人心的时代：这些方法具有巨大的潜力，可以从根本上改变计算机模拟所能实现的目标。

说明：v0.2有哪些新内容？ 对于熟悉本书v0.1的读者来说，将DP集成到NN训练中的扩展章节和关于改进物理问题学习方法的全新章节（从标度-不变量和反演开始）是强烈推荐的起点。

即将推出

作为小预告，接下来的章节将有所展示：

• 如何训练网络来推断机翼等形状周围的流体流动，并估计预测的不确定性。这将提供一个替代模型，取代传统的数值模拟。
• 如何使用模型方程作为残差来训练代表解决方案的网络，以及如何通过使用可变模拟来改进这些残差约束。
• 如何针对逆问题与完整模拟器进行更紧密的交互。例如，我们将演示如何在训练循环中利用模拟器来规避标准强化学习技术的收敛问题。
• 我们还将讨论反演对于更新步骤的重要性，以及如何利用高阶信息来加快收敛速度，并获得更精确的神经网络。

在本文中，我们将介绍将物理模型引入深度学习的不同方法，即基于物理的深度学习（PBDL）方法。这些算法变体将按照整合紧密度递增的顺序进行介绍，并将讨论不同方法的利弊。重要的是，要知道每种不同的技术在哪些场景下特别有用。

可执行代码，就在此时此地

我们将重点介绍 Jupyter 笔记本，其主要优势在于所有代码示例都可以在浏览器中当场执行。您可以修改代码，并立即查看结果--请通过[在您的浏览器中运行这个预告示例]尝试一下。

此外，Jupyter 笔记本也很不错，因为它是一种有文化的编程。

意见和建议

本书由TUM的基于物理的仿真小组维护，其中 "书 "代表数字文本和代码示例集。如果您有任何意见，请随时通过老式电子邮件与我们联系。如果您发现错误，也请告知我们！我们知道这份文件远非完美，我们渴望改进它。在此先表示感谢！另外，我们还维护了一个链接集，其中包含最新的研究论文。

](http://physicsbaseddeeplearning.org/_images/divider-mult.jpg)

图 1 数值模拟时间序列的一些可视化示例。在本书中，我们将解释如何实现使用神经网络和数值求解器的算法。

谢谢!

如果没有许多人的帮助，这个项目就不可能完成。感谢大家 以下是按字母顺序排列的名单：

• Philipp Holl
• Maximilian Mueller
• Patrick Schnell
• Felix Trost
• Nils Thuerey
• Kiwon Um

此外，还要感谢 Georg Kohl 提供漂亮的分隔线图像（参见 [KUT20] ）、Li-Wei Chen 提供机翼数据图像，以及 Chloe Paillard 对本文部分内容的校对。

引用

如果您觉得本书有用，请通过以下方式引用：

@book{thuerey2021pbdl,
  title={Physics-based Deep Learning},
  author={Nils Thuerey and Philipp Holl and Maximilian Mueller and Patrick Schnell and Felix Trost and Kiwon Um},
  url={https://physicsbaseddeeplearning.org},
  year={2021},
  publisher={WWW}
}

让我们从一个非常简化的例子开始，突出物理学习方法的一些关键能力。假设我们的物理模型是一个非常简单的方程：沿着正x轴的抛物线。

虽然非常简单，但每个x点都有两个解，即我们有两种模式，一个在x轴上方，另一个在下方，如下图左侧所示。如果我们不注意，传统的学习方法会给出一个完全错误的近似解，如下面中间图中红色线所示。通过改进学习设置，理想情况下通过使用离散化的数值求解器，我们至少可以准确表示解的一种模式（如右侧绿色线所示）。

1.1 可微物理学

下面章节的一个关键概念我们将其称之为可微物理学(Differentiable Physics，DP)。这意味着我们可以使用模型方程的领域知识，然后将这些模型的离散化版本整合到训练过程中。正如其名称所暗示的那样，具有可微的表达式对于支持神经网络的训练至关重要。

让我们通过以下示例来说明利用可微分物理进行深度学习的性质：我们希望找到一个未知函数$f^{\ast }$，其接受来自$X$的输入，并在空间$Y$中生成解，即$f^{\ast }:X\to Y$。在接下来的内容中，我们经常用$\ast$上标表示理想化的、未知的函数，其与没有这个上标的离散化、可实现的函数相对应。

另外假设我们有一个通用的微分方程${\mathcal{P}}^{\ast }:Y\to Z$（我们的“模型”方程），它编码了解的某种属性，例如我们希望匹配的某种真实世界的行为。接下来，$P^{\ast }$通常代表时间演化，但它也可以是质量守恒的约束条件（那么${\mathcal{P}}^{\ast }$将衡量散度）。但为了尽可能保持简单，在接下来的内容中，我们将关注一个将解映射回输入空间$X$的模型，即${\mathcal{P}}^{\ast }:Y\to X$。

通过使用神经网络$f$来学习未知的理想函数$f^{\ast }$，我们可以采用经典的“监督式”训练方法来通过收集数据获取$f$。这种经典的设置需要通过从$X$中采样$x$，并添加相应的解$y$来获得数据集。我们可以通过经典的数值技术来获取这些数据。然后，我们按照通常的方式利用这个数据集来训练神经网络。

与这种监督方法相比，采用可微分物理方法利用了一个事实：即我们通常可以使用物理模型 $\mathcal{P}$ 的离散化版本，并将其用于指导 $f$ 的训练。也就是说，我们希望 $f$ 知道我们的模拟器 $\mathcal{P}$，并与其进行交互。这可以极大地改善学习效果，正如我们将在下面通过一个非常简单的例子进行说明的那样（更复杂的例子将在后面介绍）。

需要注意的是，为了使可微分物理方法起作用，正如其名称所暗示的那样，$\mathcal{P}$必须是可微分的。这些微分以梯度的形式出现，是推动学习过程的关键。

1.2 寻找抛物线的反函数

为了说明监督和可微分物理方法的区别，我们考虑以下简化的场景：给定函数 $\mathcal{P}:y\to y^2$，其中 $y$ 在区间 $[0,1]$ 内，需要找到未知函数 $f$，使得对于所有 $[0,1]$ 内的 $x$，都满足 $\mathcal{P}(f(x))=x$。注意：为了使事情更有趣，我们在这里使用 $y^2$ 来表示 $\mathcal{P}$ ，而不是更常见的抛物线 $x^2$，并且对于这个简单的情况，“离散化”只需要通过在计算机中用浮点数表示 $x$ 和 $y$ 来实现。

我们知道，$f$的解可以是正或负的平方根函数（它们的分段组合也是可能的）。知道这并不是太困难，一个可行的解决方案是训练一个神经网络来近似这个逆映射$f$。以经典有监督的方式（即纯粹基于数据）来实现这一点，显然是一个不错的起点，毕竟该方法已被证明是其他各种应用（如计算机视觉）的有力工具。

import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt

对于监督训练，我们可以使用求解器$\mathcal{P}$来预先计算训练所需的解：随机选择正或负平方根中的其中一个。这类似于一般情况，我们会收集所有可用的数据（例如使用优化技术来计算求解）。这种数据收集通常不会偏向于多模式解中的某个特定模式。

# X-Data
N = 200
X = np.random.random(N)
# Generation Y-Data
sign = (- np.ones((N,)))**np.random.randint(2,size=N)
Y = np.sqrt(X) * sign

利用代码显示图形：

plt.scatter(X,Y,s=9,color='red')

此时生成的点如图所示。

现在我们将使用一个简单的keras架构定义一个网络、损失函数和训练配置，该网络使用ReLU激活函数。

# Neural network
act = tf.keras.layers.ReLU()
nn_sv = tf.keras.models.Sequential([
tf.keras.layers.Dense(10, activation=act, input_shape=(1,)),
tf.keras.layers.Dense(10, activation=act),
tf.keras.layers.Dense(1,activation='linear')])

我们可以通过一个简单的均方误差损失函数开始训练，使用Keras的fit函数进行训练：

# Loss function
loss_sv = tf.keras.losses.MeanSquaredError()
optimizer_sv = tf.keras.optimizers.Adam(learning_rate=0.001)
nn_sv.compile(optimizer=optimizer_sv, loss=loss_sv)
# Training
results_sv = nn_sv.fit(X, Y, epochs=5, batch_size= 5, verbose=1)

训练过程如下所示。

Epoch 1/5
40/40 [==============================] - 0s 1ms/step - loss: 0.5084
Epoch 2/5
40/40 [==============================] - 0s 1ms/step - loss: 0.5022
Epoch 3/5
40/40 [==============================] - 0s 1ms/step - loss: 0.5011
Epoch 4/5
40/40 [==============================] - 0s 1ms/step - loss: 0.5002
Epoch 5/5
40/40 [==============================] - 0s 1ms/step - loss: 0.5007

由于神经网络和数据集规模都非常小，训练计算非常迅速。然而如果我们检查网络的预测结果，我们会发现它与我们希望找到的解相去甚远：它在x轴两侧的数据点之间取平均值，因此无法找到满意的问题解。

下面的图表清楚地展示了这一点：它显示了原始数据（绿色）和监督解（红色）。

# Results
def draw(X,Y,Y_pre):
    plt.figure(figsize=(8,5))
    plt.plot(X,Y,'.',label = 'Data Point',color='green',zorder= 1)
    plt.plot(X,Y_pre,'.',label='Supervised',color = 'red',zorder =2)
    plt.xlabel('y',fontsize = 16)
    plt.ylabel('x',fontsize = 16)
    plt.title('standard apporach',fontsize = 16)
    plt.yticks(fontproperties='Arial',size=12)
    plt.xticks(fontproperties = 'Arial',size = 12)
    plt.legend(fontsize = 12)      

y_pre = nn_sv.predict(X)
draw(X,Y,y_pre)

结果如下图所示。

这显然是完全错误的！红线与绿色解相去甚远。

请注意，红线通常不是完全接近零，而在连续的情况下，解的两个模式应该取平均值。这是由于在这个例子中只有200个采样点，采样相对粗糙造成的。

1.3 可微物理方法

现在让我们应用可微分物理方法来找到$f$：直接将离散模型$\mathcal{P}$纳入训练中。这里没有真实数据生成步骤；我们只需要从$[0,1]$区间进行采样。

保持与先前情况相同的$x$位置，并使用与之前架构相同的新的神经网络实例nn_dp：

# X-Data
# X = X , we can directly re-use the X from above, nothing has changed...
# Y is evaluated on the fly
# Model
nn_dp = tf.keras.models.Sequential([
tf.keras.layers.Dense(10, activation=act, input_shape=(1,)),
tf.keras.layers.Dense(10, activation=act),
tf.keras.layers.Dense(1, activation='linear')])

损失函数是训练的关键点：我们直接将函数$f$融入到损失函数中。 在这种简单的情况下，损失函数loss_dp只是简单地计算预测值y_pred的平方。

后来，这里可能会发生更多的事情：我们可以对预测解进行评估有限差分模板，或者计算求解器的整个隐式时间积分步长。 这里我们有一个简单的均方误差项，其形式为$|y_{pred}^2-y_{true}{|}^2$，我们在训练期间将其最小化。 没有必要让它变得如此简单：我们可以结合的知识和数值方法越多，我们就能更好地指导训练过程。

#Loss
mse = tf.keras.losses.MeanSquaredError()
def loss_dp(y_true, y_pred):
    return mse(y_true,y_pred**2)

optimizer_dp = tf.keras.optimizers.Adam(learning_rate=0.001)
nn_dp.compile(optimizer= optimizer_dp,loss=loss_dp)

进行训练并显示结果：

#Training
results_dp = nn_dp.fit(X,X,epochs=5,batch_size = 5,verbose = 1)
y_pre = nn_dp.predict(X)
draw(X,Y,y_pre)

训练过程如下：

Epoch 1/5
40/40 [==============================] - 0s 656us/step - loss: 0.2814
Epoch 2/5
40/40 [==============================] - 0s 1ms/step - loss: 0.1259
Epoch 3/5
40/40 [==============================] - 0s 962us/step - loss: 0.0038
Epoch 4/5
40/40 [==============================] - 0s 949us/step - loss: 0.0014
Epoch 5/5
40/40 [==============================] - 0s 645us/step - loss: 0.0012

现在，网络实际上已经学会了抛物线函数的一个很好的反函数！下图用红色点显示计算解。

注意，红色的线可能是在上面，也可能是在下面，多次训练可能会得到不同的结果。比如出现下图的结果。

现在得到的结果要好多了。

这里发生了什么？

• 我们根据当前的网络预测值来评估离散模型，而不是使用预先计算好的解，从而避免了在求解中出现不希望出现的多种模式平均化现象。这样，我们就能在网络预测附近找到最佳模式，并防止求解流形中存在的模式平均化。

• 我们仍然只能得到曲线的一边！这是意料之中的，因为我们使用的是一个确定性函数来表示解。因此我们只能表示单一模式。有趣的是，是顶部模式还是底部模式由 $f$ 中权重的随机初始化决定的，多运行几次示例就能看到这种效果。要捕捉多种模式，我们需要扩展 NN 以捕捉输出的完整分布，并用额外的维度对其进行参数化。

• 在本例中，$x$ 接近零的区域通常仍然是偏离的。在这里，网络基本上学习的是抛物线一半的线性近似值。造成这种情况的部分原因是神经网络的弱点：它非常小且浅。此外沿 x 轴均匀分布的样本点会使神经网络偏向于较大的 y 值。这些点对损失的贡献更大，因此网络会投入大部分资源来减少这一区域的误差。

1.4 讨论

这是一个非常简单的例子，但却非常清晰地展示了监督学习的失败案例。虽然乍一看似乎很不真实，但许多实际的 PDE 都会表现出各种各样的这些模式，而且我们通常并不清楚我们感兴趣的解空间中存在哪些模式以及它们的位置。在这种情况下使用监督学习是非常危险的，我们可能会在不知不觉中得到这些不同模式的平均值。

流体流动中的分岔就是一个明显的好例子。蜡烛上方升起的烟雾一开始是直的，然后由于其运动中的微小扰动，开始向随机方向摆动。下面的图片通过数值扰动说明了这种情况：完全对称的设置会开始向左或向右转动，这取决于近似误差如何积累。如果将这两种模式平均化，就会产生与上述抛物线示例类似的非物理直线流。

同样，我们在许多数值解中都有不同的模式，通常重要的是恢复它们，而不是将它们平均化。因此在接下来的章节中，我们将展示如何通过可微分物理来处理更实际、更复杂的情况。

图2显示浮力驱动的流体流动中的分叉：绿色显示的“烟雾”开始以完全直线的方式上升，但随着时间的推移，微小的数值误差增长，导致涡旋向一侧（右上角）或相反方向（右下角）交替的不稳定性。

1.5 下一步

下面的每本笔记本都有一个 "下一步 "部分，比如下面这个部分，其中包含了关于从哪里开始修改代码的建议。毕竟，这些笔记本的全部意义就在于提供可随时执行的程序，作为自己实验的基础。为了减少笔记本的运行时间，示例的数据集和 NN 大小通常都很小，但它们仍然是潜在复杂大型项目的良好起点。

对于上述简单的 DP 例子：

• 本笔记本有意使用非常简单的设置。更改上述训练设置和 NN，可以获得更高质量的解，如顶部第一张图片中显示的绿色解。

• 或者尝试将设置扩展到二维情况，即抛物面。给定函数 $\mathcal{P}:(y_1,y_2)\to y_1^2+y_2^2$，找出一个反函数 $f$，使得 $[0,1]$ 中的所有 $x$ 都能得到 $\mathcal{P}(f(x))=x$。

• 如果您想在不安装任何软件的情况下进行实验，也可以[在 colab 中运行此笔记本]。

本书名为 Physics-Based Deep Learning (基于物理的深度学习)，意指将物理建模和数值模拟与基于人工神经网络的方法相结合。基于物理的深度学习的大方向代表了一个非常活跃、快速发展和令人兴奋的研究领域。接下来的章节将更详细地介绍这个主题，并为后续章节奠定基础。

图 4 了解我们的环境并预测其演化是人类面临的关键挑战之一。实现这些目标的一个重要工具是模拟，而下一代模拟将极大受益于整合深度学习组件，以便对我们的世界进行更准确的预测。

2.1 动机

从天气预报[Sto14]（见上图），到量子物理学[OMalleyBK+16]，再到等离子体聚变控制[MLA+19]，利用数值分析来获得物理模型的解已成为科学研究不可分割的一部分。近年来，机器学习技术尤其是深度神经网络在多个领域取得了令人瞩目的成就：从图像分类[KSH12]到自然语言处理[RWC+19]，以及最近的蛋白质折叠[Qur19]。该领域充满活力，发展迅速，前景广阔。

2.1.1 取代传统的模拟？

这些深度学习（DL）方法的成功案例引起了人们的关注，认为这种技术有可能取代传统的、模拟驱动的科学方法。例如，最近的研究表明，基于神经网络的代理模型达到了实际工业应用（如机翼流）所需的精度[CT21]，同时在运行时间方面优于传统求解器几个数量级。

我们是否可以不依赖于仔细设计的基于基础原理的模型，而是通过处理足够大的数据集来提供正确的答案呢？正如我们将在接下来的章节中展示的那样，这种担忧是没有根据的。相反，对于下一代模拟系统来说，将传统的数值技术与深度学习方法相结合是至关重要的。

这种组合的重要性有一个核心原因，那就是深度学习方法强大的同时，也能从物理模型形式的领域知识中获益匪浅。深度学习技术和神经网络都比较新颖，有时应用起来比较困难，诚然，将我们对物理过程的理解妥善地融入到深度学习算法中往往并非易事。

在过去的几十年中，人们已经开发出高度专业和精确的离散方法来求解如 Navier-Stokes、Maxwell 或 Schroedinger 方程等基本模型方程。对离散方法进行的看似微不足道的改动却能决定在计算结果中是否能够看到关键的物理现象。本书非但不会抛弃数值数学领域已开发出的强大方法，反而会说明在应用深度学习时尽可能使用这些方法大有裨益。

2.1.2 黑匣子和魔法？

不熟悉 DL 方法的人通常会将神经网络与 黑盒子 联系在一起，并将训练过程视为人类无法理解的东西。然而这些观点通常源于道听途说和对这一主题了解不够。

相反，这种情况在科学界很常见：我们正在面对一类新的方法，而 所有的细节 还没有完全解决。这种情况在各种科学进步中都很常见。数值方法本身就是一个很好的例子。1950 年前后，数值逼近和求解器的地位十分艰难。例如，引用 H. Goldstine 的说法，数值不稳定性被认为是 未来焦虑的持续来源[Gol90]。现在，我们已经很好地掌握了这些不稳定性，数值方法无处不在，而且已经非常成熟。

因此，我们必须意识到，从某种程度上说，深度学习方法并不神奇或超凡脱俗。它们只是另一套数值工具。尽管如此，它们显然还是相当新的工具，而且现在绝对是我们所拥有的用于解决非线性问题的最强大的工具集。不能因为所有的细节还没有完全解决，还没有写得很好，就阻止我们将这些强大的方法纳入我们的数值工具箱。

2.1.3 协调DL和模拟

退一步说，本书的目的是利用我们所掌握的所有强大的数值模拟技术，并将它们与深度学习结合起来使用。因此，本书的核心目标是将以数据为中心的观点与物理模拟重新结合起来。

我们将讨论的关键问题包括：

• 解释如何使用深度学习技术来求解PDE问题
• 如何将它们与现有的物理学知识结合起来
• 不抛弃我们已有的关于数值方法的知识。

与此同时，值得注意的是，我们不会涵盖的内容包括：

• 深度学习和数值模拟的介绍，我们不打算对这一领域的研究文章进行广泛的研究

由此产生的方法在改进数值方法方面具有巨大潜力：例如，在求解器重复求解某个定义明确的问题领域中的案例时，投入大量资源训练一个支持重复求解的神经网络就非常有意义。基于该网络的特定领域专用性，这种混合网络的性能会大大超过传统的通用求解器。尽管还有很多问题有待解决，但第一批出版物已经证明，这一目标并不遥远[KSA+21, UBH+20]。

另一种看法是，我们自然界的所有数学模型都是理想化的近似值，都包含误差。为了获得非常好的模型方程，我们已经付出了很多努力，但为了向前迈进一大步，DL 方法提供了一个非常强大的工具，可以缩小与现实的差距[AAC+19]。

2.2 分类

在基于物理的深度学习领域，我们可以区分出从目标约束、组合方法、优化到应用的各种不同方法。更具体地说，所有方法要么针对正向模拟（预测状态或时间演化），要么针对反向问题（例如，从观测结果中获取物理系统的参数）。

无论我们考虑的是正向问题还是反向问题，以下主题最关键的区别在于深度学习技术与领域知识（通常以偏微分方程（PDE）的形式）之间的整合方式的性质。基于物理学的深度学习（PBDL）技术大致可分为以下三类：

• 监督学习（Supervised:）：数据由物理系统（真实或模拟）产生，但不存在进一步的交互。这是典型的机器学习方法。
• 损失项（Loss-terms）：物理动力学（或其部分）通常以可微分算子的形式被编码在损失函数中。学习过程可以重复评估损失，并且通常从基于 PDE 的表达式中获得梯度。这些软约束有时也被称为“具有物理信息”的训练方法。
• 交织方法（Interleaved）：完整的物理模拟与深度神经网络的输出交织并结合在一起。这需要一个完全可微分的模拟器，代表了物理系统与学习过程之间最紧密的耦合。交织式可微分物理方法在时间演化中尤为重要，其可以为动力学的未来行为做出估计。

因此，可以根据正向求解与反向求解以及物理模型在训练深度神经网络的优化循环中的紧密程度进行分类。在这里，特别是利用“可微分物理”进行交织的方法允许深度学习和数值模拟方法之间的非常紧密的集成。

2.2.1 命名

值得指出的是，在接下来的内容中，我们所称谓的可微分物理学 (Differentiable Physics，DP）在其他资源和研究论文中有多种不同的名称。可微分物理学的名称受到深度学习中可微分编程范式的启发。在这里，例如我们还有可微分渲染方法，用于模拟光线如何形成我们作为人类所看到的图像。相比之下，我们将从现在开始专注于物理模拟，因此才有了这个名字。

然而，来自其他背景的人更常用其他名称。例如可微分物理方法等效于使用伴随方法（adjoint method），并将其与深度学习过程相结合。实际上，它也等效于将反向传播/反向模式微分应用于数值模拟。然而正如上面所提到的，出于对深度学习观点的考虑，从现在开始我们将统称这些方法为可微分物理方法。

2.3 展望未来

物理模拟是一个庞大的领域，我们无法涵盖所有可能的物理模型和模拟类型。

注意本书的重点在于：

• 基于场的模拟（没有拉格朗日方法）
• 与深度学习的结合（存在其他许多有趣的机器学习技术，但这里不讨论）
• 实验留作展望（即用真实世界的观测数据取代合成数据）

值得注意的是，我们正从一些非常基本的构建模块开始构建这些方法。以下是跳过后面章节的一些注意事项。

提示：如果满足以下条件，可以跳过...

• 非常熟悉数值方法和偏微分方程求解器，并且希望立即开始DL主题。那么监督训练章节是一个很好的起点。
• 另一方面，如果已经深入NNs&Co，并且想跳过研究相关主题，我们建议从物理损失项章节开始，这为下一章奠定了基础。

不过，在这两种情况下，简单了解一下我们在符号和缩写 章节中的 符号 并不会有什么坏处！

2.4 实现

这段文字还介绍了一系列深度学习和模拟API的概述。我们将使用流行的深度学习API，[如pytorch](https://pytorch.org]和 tensorflow，并且额外介绍了可微分模拟框架 Φ_Flow (phiflow) 。其中一些示例也使用了JAX。因此，通过学习这些示例，您应该对当前API中提供的功能有一个很好的概览，从而可以为新任务选择最适合的API。

由于我们（在大多数Jupyter笔记本示例中）处理的是随机优化问题，因此下面的许多代码示例每次运行时都会产生稍微不同的结果。这在神经网络训练中非常常见，但在执行代码时要牢记这一点。这也意味着文本中讨论的数字可能与您重新运行示例后看到的数字并不完全一致。

2.5 模型与方程

下面我们将简要介绍深度学习（真的是非常非常简要！），主要是为了介绍符号表示法。

此外，我们将在下面讨论一些模型方程。请注意，与其他一些文本和API不同，我们将避免使用Model来表示经过训练的神经网络。它们将被称为NNs（神经网络）或networks（网络）。这里的Model（模型）通常表示一组描述物理效应的模型方程，它们通常是偏微分方程（PDEs）。

2.5.1 深度学习与神经网络

在本书中，我们专注于与物理模型的联系，而深度学习有很多很棒的介绍。因此，我们将保持简短：深度学习的目标是近似一个未知函数： $$f^{\ast }(x)=y^{\ast }\text{\hspace{1em}(1)}$$ 其中，$y^{\ast }$ 表示参考解或真实解。$f^{\ast }(x)$应该用NN表示$f(x;\theta )$来近似。我们通常借助损失函数 $L(y,y^{\ast })$的某种变体来确定 $f$，其中 $y=f(x;\theta )$ 是 NN 的输出。这就给出了一个最小化问题，即找到 $f(x;\theta )$使得 $e$ 最小化。在最简单的情况下，我们可以使用 $L^2$ 误差，得出 $$\arg \underset{\theta }{\min }{\left|f(x;\theta )-y^{\ast }\right|}_2^2\text{\hspace{1em}(2)}$$ 通常，我们使用随机梯度下降（SGD）优化器（例如Adam [KB14]）对其进行优化，即训练。我们将依赖于自动微分来计算损失函数$L$相对于权重$\theta$的梯度$\partial L/\partial \theta$。对于梯度计算来说，这个函数是标量且非常重要，损失函数通常也被称为误差函数、成本函数或目标函数。

对于训练，我们要区分：从某种分布中提取的训练数据集、验证数据集（来自相同的分布，但数据不同），以及与训练数据集有某种不同分布的测试数据集。后一种区别非常重要。对于测试集，我们需要 超出分布范围（OOD）的数据，以检查训练模型的泛化效果。需要注意的是，这就为测试数据集提供了一个巨大的可能性范围：从肯定会起作用的微小变化，到基本上保证会失败的完全不同的输入。虽然没有黄金标准，但在生成测试数据时仍需谨慎。

如果以上内容对你来说还不够明显，我们强烈建议阅读《深度学习》这本书的第6章至第9章，尤其是关于多层感知器（MLPs）和卷积神经网络（CNNs）的部分。你可以在深度学习的链接中找到这些章节，分别是MLPs和Conv-Nets，即卷积神经网络（CNNs）。

分类问题和回归问题之间的经典 ML 区别在这里并不重要：下面的内容中，我们只讨论回归问题。

2.5.2 作为物理模型的偏微分方程

接下来的部分将简要介绍稍后在深度学习示例中所使用到的一些模型方程。我们的目标通常是连续的偏微分方程，用 ${\mathcal{P}}^{\ast }$ 表示，其解处在 $d\in 1,2,3$ 维的空间域 $\Omega \subset {\mathbb{R}}^d$ 中。此外我们还经常考虑方程在有限时间间隔 $t\in {\mathbb{R}}^{+}$ 内的演化过程。相应的物理场可以是 d 维向量场，例如 $\mathbf{u}：{\mathbb{R}}^d\times {\mathbb{R}}^{+}\to {\mathbb{R}}^d$，也可以是标量场 $\mathbf{p}：{\mathbb{R}}^d\times {\mathbb{R}}^{+}\to \mathbb{R}$。向量的分量通常用$x，y，z$下标表示，即$\mathbf{v}=(v_x，v_y，v_z{)}^T$（当$d=3$时），而位置则用$x\in \Omega$表示。

为了得到 ${\mathcal{P}}^{\ast }$ 的唯一解，我们需要指定合适的初始条件，通常是在 $t=0$ 时的所有相关物理量，以及 $\Omega$ 边界的边界条件，在下文中用 $\Gamma$ 表示。

${\mathcal{P}}^{\ast }$ 表示一个连续的公式，我们对其连续性做了轻微的假设，通常我们会假设其存在一阶导数和二阶导数。

然后，我们可以使用数值方法，通过离散化获得平滑函数的近似值，如 ${\mathcal{P}}^{\ast }$。这必然会引入离散化误差，我们希望误差越小越好。这些误差可以用与精确分析解的偏差来衡量，对于离散模拟 PDEs，它们通常表示为截断误差 $O(\Delta x^k)$ 的函数，其中 $\Delta x$ 表示离散化的空间步长。同样，我们通常通过时间步长 $\Delta t$ 进行时间离散化。

我们通过执行大小为 $\Delta t$ 的步长来求解一个离散的偏微分方程$\mathcal{P}$。解可以表示为 $\mathbf{u}$ 及其导数的函数：$\mathbf{u}(\mathbf{x}，t+\Delta t)=\mathcal{P}\left({\mathbf{u}}_x，{\mathbf{u}}_{xx}，\dots {\mathbf{u}}_{xx\dots x}\right)$，其中 ${\mathbf{u}}_x$ 表示空间导数 $\partial \mathbf{u}(\mathbf{x},t)/\partial \mathbf{x}$。

接下来，我们将对模型方程进行概述，然后在后续章节及性能实际模拟和实现示例。

2.5.3 一些PDE示例

以下PDEs是很好的例子，我们将在不同的设置中使用它们来展示如何将它们纳入到深度学习中。

Burgers方程

我们经常以一维或二维的Burgers方程作为起点。这是一个经过广泛研究的偏微分方程，与Navier-Stokes方程不同，该方程不包含任何额外的约束条件（如质量守恒）。因此，它会导致有趣的激波形成。Burgers方程包含一个对流项（运动/输运）和一个扩散项（由热力学第二定律引起的耗散）。

二维Burgers方程表示为： $$\begin{array}{l}\frac{\partial u_x}{\partial t}+\mathbf{u}\cdot \nabla u_x=\nu \nabla \cdot \nabla u_x+g_x\\ \frac{\partial u_y}{\partial t}+\mathbf{u}\cdot \nabla u_y=\nu \nabla \cdot \nabla u_y+g_y\end{array}\text{\hspace{1em}(3)}$$ 式中，$\nu$及$\mathbf{g}$分别表示扩散常数和外力。

在没有外力的情况下，，将单个一维速度分量表示为$u=u_x$，可以得到burgers方程在一维空间中简单的变体： $$\frac{\partial u}{\partial t}+u\nabla u=\nu \nabla \cdot \nabla u\text{\hspace{1em}(4)}$$ Navier-Stokes方程

就复杂性而言，Navier-Stokes方程是一个很好的下一步，它是一个成熟的流体模型。除了动量守恒方程（类似于Burgers方程）外，该方程还包括质量守恒方程。这可以防止激波的形成，但同时为数值方法引入了新的挑战，即添加了无散度运动的硬约束条件。

在二维情况下，没有任何外部力的Navier-Stokes方程可以写成： $$\begin{array}{rl}\frac{\partial u_x}{\partial t}+\mathbf{u}\cdot \nabla u_x& =-\frac{\Delta t}{\rho }\nabla p+\nu \nabla \cdot \nabla u_x\\ \frac{\partial u_y}{\partial t}+\mathbf{u}\cdot \nabla u_y& =-\frac{\Delta t}{\rho }\nabla p+\nu \nabla \cdot \nabla u_y\\ \nabla \cdot \mathbf{u}& =0\end{array}\text{\hspace{1em}(5)}$$ 其中，像前面一样，𝜈 为粘度。

NS方程一个有趣的变体是通过包含Boussinesq近似来考虑温度对流体密度的影响。使用一个标记场$v$来指示高温区域，可以得到以下方程组： $$\begin{array}{rl}\frac{\partial u_x}{\partial t}+\mathbf{u}\cdot \nabla u_x& =-\frac{\Delta t}{\rho }\nabla p\\ \frac{\partial u_y}{\partial t}+\mathbf{u}\cdot \nabla u_y& =-\frac{\Delta t}{\rho }\nabla p+\xi v\\ \mathrm{subject}\mathrm{to}\nabla \cdot \mathbf{u}& =0\\ \frac{\partial v}{\partial t}+\mathbf{u}\cdot \nabla v& =0\end{array}\text{\hspace{1em}(6)}$$ 其中𝜉表示浮力的强度。

最后，3D中的Navier-Stokes模型给出了以下方程组： $$\begin{array}{rl}\frac{\partial u_x}{\partial t}+\mathbf{u}\cdot \nabla u_x& =-\frac{\Delta t}{\rho }\nabla p+\nu \nabla \cdot \nabla u_x\\ \frac{\partial u_y}{\partial t}+\mathbf{u}\cdot \nabla u_y& =-\frac{\Delta t}{\rho }\nabla p+\nu \nabla \cdot \nabla u_y\\ \frac{\partial u_z}{\partial t}+\mathbf{u}\cdot \nabla u_z& =-\frac{\Delta t}{\rho }\nabla p+\nu \nabla \cdot \nabla u_z\\ \mathrm{subject}\mathrm{to}\nabla \cdot \mathbf{u}& =0\end{array}\text{\hspace{1em}(7)}$$

2.5.4 正向模拟

在我们真正开始学习方法之前，有必要介绍一下使用上述模型方程的最基本的变体：常规的 正向 模拟，即从一组初始条件开始，随着时间的推移，用模型方程的离散版本演化系统的状态。我们将展示如何对一维burgers方程和二维的Navier-Stokes模拟进行这种正向模拟。

2.6 用phiflow对Burgers方程进行正演模拟

本章将介绍如何运行“正向”即常规模拟，即从给定的初始状态开始，通过数值逼近来近似后续状态，并介绍Φ_Flow框架。Φ_Flow提供了一组可微分的构建模块，可以直接与深度学习框架进行接口，因此非常适合本书的主题。在深入和更复杂的集成之前，本章（以及下一章）将展示如何使用Φ_Flow进行常规模拟。随后，我们将展示如何轻松将这些模拟与神经网络进行耦合。

ΦFlow的主要代码库（以下简称为"phiflow"）位于https://github.com/tum-pbs/PhiFlow，您可以在https://tum-pbs.github.io/PhiFlow/找到其他API文档和示例。

对于这个jupyter笔记本（以及后续的笔记本），你可以在第一个段落的末尾找到一个“[在Colab中运行]”的链接（或者你可以使用页面顶部的启动按钮）。这将在Colab中加载PBDL GitHub仓库中的最新版本的笔记本，你可以立即执行它：在Colab中运行。

2.6.1 Model

我们使用Burgers方程作为物理模型。这个方程是一个非常简单、非线性且非平凡的模型方程，可以导致有趣的冲击波形成。因此它是我们进行实验的一个很好的起点，该方程的一维版本（方程(4)）如下所示： $$\frac{\partial u}{\partial t}+u\nabla u=\nu \nabla \cdot \nabla u$$

2.6.2 导入和加载Phiflow

让我们先做一些准备工作：首先，我们将导入phiflow库，更具体地说是用于流体模拟的numpy运算符：phi.flow（用于DL框架X的可微分版本则通过phi.X.flow加载）。

注意： 下面的第一条带有"!"前缀的命令将通过pip在您的Python环境中安装来自GitHub的phiflow Python包。我们假设您尚未安装phiflow，但如果您已经安装了，请注释掉第一行（对于后续的所有笔记本也是如此）。

!pip install --upgrade --quiet phiflow==2.2
from phi.flow import *

注意：这里如果提示dtype错误，则需要修改C:\ProgramData\anaconda3\Lib\site-packages\phi\math\backend_dtype.py文件，修改第135行为 DType(object): np.dtype('O'),，修改第139行为_FROM_NUMPY[np.bool_] = DType(bool)。主要原因是numpy升级后抛弃了np.object及np.bool。如果是低版本的numpy则不会报错。

接下来，可以定义并初始化必要的常量（用大写字母表示）：模拟域具有N=128个单元格作为1D速度$u$的离散点，这些离散点位于在周期性域$\Omega$上的区间$[-1,1]$内。我们将使用32个时间步长STEPS，时间间隔为1，因此我们得到DT=1/32。此外，我们将使用粘度NU为$\nu =0.01/\pi$。

我们还将定义一个初始状态，由numpy数组INITIAL_NUMPY中的$-\text{sin}(\pi x)$给出，我们将在下一个单元格中使用它来初始化模拟中的速度$u$。此初始化将在计算域的中心产生一个漂亮的冲击波。

Phiflow是面向对象的，以场数据形式（由张量对象在内部表示）为中心。也就是说，您通过构建一些网格来组装模拟，并在时间步长内更新它们。

Phiflow在内部使用带有命名维度的张量。这在后面处理具有附加批处理和通道维度的2D模拟时将非常方便，但现在我们只需将1D数组简单地转换为具有单个空间维度'x'的phiflow张量即可。

N = 128
DX = 2./N
STEPS = 32
DT = 1./STEPS
NU = 0.01/(N*np.pi)
# initialization of velocities， cell centers of a CenteredGrid have DX/2 offsets for linspace()
INITIAL_NUMPY = np.asarray( [-np.sin(np.pi * x) for x in np.linspace(-1+DX/2，1-DX/2，N)] ) # 1D numpy array
INITIAL = math.tensor(INITIAL_NUMPY， spatial('x') ) # convert to phiflow tensor

接下来，我们从转换为张量的INITIAL numpy数组中初始化一个1D velocity 网格。 我们通过bounds参数$[-1,1]$来指定域$\Omega$的范围，并且网格使用周期性边界条件（extrapolation.PERIODIC）。这两个属性是张量和网格之间的主要区别：后者具有边界条件和物理范围。

为了便于说明，我们还将输出有关速度对象的一些信息：它是一个大小为128的phi.math张量。请注意，实际的网格内容包含在网格的values中。下面我们使用numpy()函数将phiflow张量的内容转换为numpy数组，并输出其中的五个条目。对于具有更多维度的张量，我们需要在这里指定额外的维度，例如，对于一个2D速度场，我们会使用'y,x,vector'（对于只有单个维度的张量，我们可以省略它）。

velocity = CenteredGrid(INITIAL， extrapolation.PERIODIC， x=N， bounds=Box('x'，-1:1))
vt = advect.semi_lagrangian(velocity， velocity， DT)
print(`Velocity tensor shape: ` + format( velocity.shape )) # == velocity.values.shape
print(`Velocity tensor type: ` + format( type(velocity.values) ))
print(`Velocity tensor entries 10 to 14: ` + format( velocity.values.numpy('x')[10:15] ))

输出结果为：

Velocity tensor shape: (x=128)
Velocity tensor type: <class 'phi.math._tensors.CollapsedTensor'> 
Velocity tensor entries 10 to 14: [0.49289819 0.53499762 0.57580819 0.61523159 0.65317284]

2.6.3 运行模拟

现在我们可以准备运行模拟了。为了计算模型方程中的扩散和对流分量，我们可以简单地调用phiflow中现有的diffusion和semi_lagrangian运算符：diffuse.explicit(u,...)通过中心差分计算模型中项$\nu \nabla \cdot \nabla u$的显式扩散。接下来，advect.semi_lagrangian(f,u)用于通过速度u对任意场f进行稳定的一阶逼近输运。在模型中有$\partial u/\partial t+u\nabla f$，因此可以使用semi_lagrangian函数在实现中将速度与自身进行输运：

velocities = [velocity]
age = 0.
for i in range(STEPS):
    v1 = diffuse.explicit(velocities[-1]， NU， DT)
    v2 = advect.semi_lagrangian(v1， v1， DT)
    age += DT
    velocities.append(v2)
print(`New velocity content at t={}: {}`.format( age， velocities[-1].values.numpy('x，vector')[0:5] ))

输出结果为：

New velocity content at t=1.0: [[0.0057228 ]
[0.01716715]
[0.02861034]
[0.040052 ]
[0.05149214]]

在这里，我们实际上将所有时间步收集到列表velocities中。一般来说，这并不是必需的（并且对于长时间运行的模拟可能会消耗大量内存），但在这里可以用于后续速度状态演变的绘制。

print语句输出了一些速度值，已经表明我们的模拟正在发生一些事情，但很难从这些数字中直观地了解 PDE 的行为。因此，让我们将随时间变化的状态可视化，以显示正在发生的事情。

2.6.4 可视化

我们可以很容易地通过图表直观地看到这种一维情况：下面的代码用蓝色显示初始状态，然后用绿色、青色和紫色显示 $10/32,20/32,1$ 的时间。

# get `velocity.values` from each phiflow state with a channel dimensions， i.e. `vector`
vels = [v.values.numpy('x，vector') for v in velocities] # gives a list of 2D arrays

import pylab

fig = pylab.figure().gca()
fig.plot(np.linspace(-1，1，len(vels[ 0].flatten()))， vels[ 0].flatten()， lw=2， color='blue'， label=`t=0`)
fig.plot(np.linspace(-1，1，len(vels[10].flatten()))， vels[10].flatten()， lw=2， color='green'， label=`t=0.3125`)
fig.plot(np.linspace(-1，1，len(vels[20].flatten()))， vels[20].flatten()， lw=2， color='cyan'， label=`t=0.625`)
fig.plot(np.linspace(-1，1，len(vels[32].flatten()))， vels[32].flatten()， lw=2， color='purple'，label=`t=1`)
pylab.xlabel('x'); pylab.ylabel('u'); pylab.legend()

输出结果如图所示。

这很好地显示了在我们的区域中心形成的激波，它是由两个初始速度 凸起（左侧的正凸起（向右移动）和中心右侧的负凸起（向左移动））碰撞形成的。

由于这些线条可能会有很多重叠，我们在接下来的章节中还将使用另一种可视化方法，即在二维图像中显示所有时间步长的演变过程。我们的一维域将沿 Y 轴显示，沿 X 轴的每个点代表一个时间步长。

下面的代码将我们的速度状态集合转换为二维数组，将单个时间步长重复 8 次，以使图像更宽。当然，这完全是可有可无的，但却能让我们更容易地看到汉堡模拟中发生的事情。

import os
def show_state(a， title):

    # we only have 33 time steps， blow up by a factor of 2^4 to make it easier to see
    # (could also be done with more evaluations of network)
    a = np.expand_dims(a， axis=2)
    for i in range(4):
        a = np.concatenate([a， a]， axis=2)
    a = np.reshape(a， [a.shape[0]， a.shape[1]*a.shape[2]])
    #print(`Resulting image size` +format(a.shape))
    fig， axes = pylab.subplots(1， 1， figsize=(16， 5))
    im = axes.imshow(a， origin='upper'， cmap='inferno')
    pylab.colorbar(im)
    pylab.xlabel('time')
    pylab.ylabel('x')
    pylab.title(title)

vels_img = np.asarray(np.concatenate(vels， axis=-1)， dtype=np.float32)
# save for comparison with reconstructions later on
os.makedirs(`./temp`， exist_ok=True)
np.savez_compressed(`./temp/burgers-groundtruth-solution.npz`，np.reshape(vels_img， [N， STEPS+1]))

show_state(vels_img， `Velocity`)

运行结果为：

至此，phiflow 的首次模拟结束。它并不太复杂，但正因为如此，它为下一章评估和比较不同的基于物理的深度学习方法提供了一个良好的起点。不过在此之前，我们将在下一节针对更复杂的模拟类型进行研究。

2.6.5 后续步骤

在此模拟设置的基础上进行一些尝试：

• 请随意尝试 - 上述设置非常简单，您可以更改模拟参数或初始化。例如，可以通过 Noise() 使用噪声场来获得更混乱的结果（参见上文 速度 单元中的注释）。
• 更复杂一点：将模拟扩展到 2D（或更高）。这需要对整个程序进行修改，但上述所有运算符都支持更高维度。在尝试之前，您可能需要查看下一个示例，它涵盖了一个 2D Navier-Stokes 案例。

2.7 Navier-Stokes 正向模拟

现在让我们举一个更复杂的例子：基于纳维-斯托克斯方程的流体模拟。使用 Φ_Flow (phiflow)仍然非常简单，因为这里存在所有步骤的可微分算子。纳维-斯托克斯方程（不可压缩形式）引入了额外的压力场 $p$，以及质量守恒约束。我们还将随着气流移动一个标记场，这里用 $d$ 表示。它表示温度较高的区域，并通过浮力因子$\xi$产生作用力： $$\begin{array}{rl}\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}& =-\frac{1}{\rho }\nabla p+\nu \nabla \cdot \nabla \mathbf{u}+(0,1{)}^T\xi d\text{s.t.}\nabla \cdot \mathbf{u}=0,\\ \frac{\partial d}{\partial t}+\mathbf{u}\cdot \nabla d& =0\end{array}$$

在这里，我们的目标是无压缩流 (即 $\rho =\text{const}$)， 并通过项 $(0，1{)}^T\xi d$ 使用一个简单的浮力模型 (Boussinesq 近似)。这在不显式计算 $\rho$ 的情况下为无压缩求解器近似了密度的变化。我们假设一个沿 y 方向作用的重力力， 通过向量 $(0，1{)}^T$。

我们将在一个封闭域上求解这个 PDE，该域的速度边界条件为 $\mathbf{u}=0$，压力边界条件为 $\frac{\partial p}{\partial x}=0$，该域的物理尺寸为 $100\times 80$ 单位。[run in colab]

2.7.1 实现

与上一节一样，第一条带!前缀的命令是在 python 环境中通过pip安装来自 GitHub 的 phiflow python 软件包。(如有必要，可跳过或修改此命令）。

from phi.flow import *  # The Dash GUI is not supported on Google colab, ignore the warning
import pylab

2.7.2 设置模拟

下面的代码会设置一些常数，这些常数用大写字母表示。我们将使用 $40\times 32$ 单元来离散域，通过 $\nu$ 引入轻微的粘性，并定义时间步长为 $\Delta t=1.5$。

我们首先创建一个 CenteredGrid，它由一个 Sphere 几何对象初始化。这将表示热烟雾生成的入流区域 INFLOW。

DT = 1.5
NU = 0.01

INFLOW = CenteredGrid(Sphere(center=tensor([30,15], channel(vector='x,y')), radius=10), extrapolation.BOUNDARY, x=32, y=40, bounds=Box(x=(0,80),y=(0,100))) * 0.2

流入的烟雾将被注入第二个居中的网格 smoke 中，该网格代表上面的标记字段 $d$。请注意，我们在上面定义了一个大小为 $100\times 80$ 的 Box。这是在我们的模拟中以空间单位表示的物理尺度，也就是说，$1$ 的速度将使烟雾密度每 1 个时间单位移动 1 个单位，这可能大于或小于离散网格中的一个单元，具体取决于 x,y 的设置。您可以将模拟网格参数设置为直接类似于现实世界的单位，或者考虑适当的换算因子。

上面的入流球体已经使用“世界”坐标: 它位于第一个轴上的 $x=30$， $y=15$ (在 $100\times 80$ 的域框内)。

接下来，我们为要模拟的量创建网格。在本例中，我们需要速度场和烟雾密度场。

smoke = CenteredGrid(0, extrapolation.BOUNDARY, x=32, y=40, bounds=Box(x=(0,80),y=(0,100)))  # sampled at cell centers
velocity = StaggeredGrid(0, extrapolation.ZERO, x=32, y=40, bounds=Box(x=(0,80),y=(0,100)))  # sampled in staggered form at face centers 

我们以交错形式对单元中心的烟场和速度进行采样。交错网格内部包含 2 个不同维度的居中网格，可通过 unstack 函数转换为居中网格（或简单的 numpy 数组），如上文链接中所述。

接下来，我们定义模拟的更新步骤，调用必要的函数将流体系统的状态向前推进 dt。下一个单元将计算一个这样的步骤，并绘制一个模拟帧后的标记密度图。

def step(velocity, smoke, pressure, dt=1.0, buoyancy_factor=1.0):
    smoke = advect.semi_lagrangian(smoke, velocity, dt) + INFLOW
    buoyancy_force = (smoke * (0, buoyancy_factor)).at(velocity)  # resamples smoke to velocity sample points
    velocity = advect.semi_lagrangian(velocity, velocity, dt) + dt * buoyancy_force
    velocity = diffuse.explicit(velocity, NU, dt)
    velocity, pressure = fluid.make_incompressible(velocity)
    return velocity, smoke, pressure

velocity, smoke, pressure = step(velocity, smoke, None, dt=DT)

print(`Max. velocity and mean marker density: ` + format( [ math.max(velocity.values) , math.mean(smoke.values) ] ))

pylab.imshow(np.asarray(smoke.values.numpy('y,x')), origin='lower', cmap='magma')

在这个 step() 调用中发生了很多事情：我们平移了烟场，通过Boussinesq 模型添加了一个向上的力，平移了速度场，最后通过压力求解使其无发散。

布尔辛斯克模型使用(0, buoyancy_factor) 元组乘法将烟场转化为交错的 2 分力场，并通过 at() 函数在速度分量的位置进行采样。该函数可确保交错速度分量的各个力分量被正确插值。请注意，这也会直接确保保留原始网格的边界条件。它还会在内部对生成的力网格执行 StaggeredGrid(..., extrapolation.ZERO,...) 内插。

在 make_incompressible 中的压力投影步骤通常是上述序列中计算成本最高的步骤。它根据域的边界条件求解泊松方程，并根据计算出的压力梯度更新速度场。

为了测试，我们还打印了更新后的速度平均值和最大密度。从生成的图像中可以看到，第一轮烟雾区域有轻微的向上运动（此处尚未显示）。

2.7.3 数据类型和尺寸

我们在这里为模拟字段创建的变量是Grid 类的实例。与张量一样，网格也有 shape 属性，该属性列出了所有批次、空间和通道维度。phiflow中的形状不仅存储了维度的大小，还存储了它们的名称和类型。

print(f"Smoke: {smoke.shape}")
print(f"Velocity: {velocity.shape}")
print(f"Inflow: {INFLOW.shape}, spatial only: {INFLOW.shape.spatial}")

请注意，这里的 phiflow 输出显示了维度的类型，例如，${}^S$ 表示空间维度，${}^V$ 表示向量维度。在后面的学习中，我们还将引入批量维度。

我们可以通过 .sizes获取形状对象的实际内容，或者通过 .get_size('dim')查询特定维度 dim 的大小。下面是两个示例：

print(f"Shape content: {velocity.shape.sizes}")
print(f"Vector dimension: {velocity.shape.get_size('vector')}")

可以使用 values 属性访问网格值。这是与 phiflow 张量对象的一个重要区别，后者没有values，如下代码示例所示。

print("Statistics of the different simulation grids:")
print(smoke.values)
print(velocity.values)

# in contrast to a simple tensor:
test_tensor = math.tensor(numpy.zeros([3, 5, 2]), spatial('x,y'), channel(vector="x,y"))
print("Reordered test tensor shape: " + format(test_tensor.numpy('vector,y,x').shape) ) # reorder to vector,y,x 
#print(test_tensor.values.numpy('y,x')) # error! tensors don't return their content via ".values"

网格还有更多的特性，这里 有详细记录。还要注意的是，交错网格具有非均匀形状，因为网格面的数量与单元格的数量不相等（在本例中，x 部分的单元格数量是 40 的 31 次方，而 y 部分的单元格数量是 39 的 32 次方）。INFLOW "网格的尺寸自然与 "smoke "网格相同。

2.7.4 时间演化

有了这种设置，我们就可以通过反复调用 step 函数，轻松地将模拟时间向前推进一些。

for time_step in range(10):
    velocity, smoke, pressure = step(velocity, smoke, pressure, dt=DT)
    print('Computed frame {}, max velocity {}'.format(time_step , np.asarray(math.max(velocity.values)) ))

现在，热羽流开始上升：

pylab.imshow(smoke.values.numpy('y,x'), origin='lower', cmap='magma')

让我们再计算并展示几个模拟步骤。由于流入的水流偏离中心向左（X 位置为 30），因此当羽流撞击到域的顶壁时会向右弯曲。

steps = [[ smoke.values, velocity.values.vector[0], velocity.values.vector[1] ]]
for time_step in range(20):
  if time_step<3 or time_step%10==0: 
    print('Computing time step %d' % time_step)
  velocity, smoke, pressure = step(velocity, smoke, pressure, dt=DT)
  if time_step%5==0:
    steps.append( [smoke.values, velocity.values.vector[0], velocity.values.vector[1]] )

fig, axes = pylab.subplots(1, len(steps), figsize=(16, 5))
for i in range(len(steps)):
    axes[i].imshow(steps[i][0].numpy('y,x'), origin='lower', cmap='magma')
    axes[i].set_title(f"d at t={i*5}")

我们还可以看看速度。上面的 steps 列表已经将速度的 vector[0] 和 vector[1] 分量存储为 numpy 数组，接下来我们可以展示一下。

fig, axes = pylab.subplots(1, len(steps), figsize=(16, 5))
for i in range(len(steps)):
    axes[i].imshow(steps[i][1].numpy('y,x'), origin='lower', cmap='magma')
    axes[i].set_title(f"u_x at t={i*5}")
    
fig, axes = pylab.subplots(1, len(steps), figsize=(16, 5))
for i in range(len(steps)):
    axes[i].imshow(steps[i][2].numpy('y,x'), origin='lower', cmap='magma')
    axes[i].set_title(f"u_y at t={i*5}")

这里看起来很简单，但这个模拟设置却是一个功能强大的工具。仿真可以很容易地扩展到更复杂的情况或三维空间，而且它已经与深度学习框架的反向传播管道完全兼容。

在接下来的章节中，我们将展示如何使用这些模拟来训练 NN，以及如何通过训练好的 NN 来引导和修改它们。这将说明通过在循环中加入求解器，我们可以在多大程度上改进训练过程，尤其是当求解器是可变的时候。在进入这些更复杂的训练过程之前，我们将在下一章介绍一种更简单的监督方法。这一点非常重要：即使是针对基于物理的高级学习设置，有效的监督训练始终是第一步。

2.7.5 下一步

您可以在此基础上创建各种流体模拟。例如，尝试改变空间分辨率、浮力系数和模拟运行的总长度。

2.8 优化和收敛性

本章将概述不同优化算法的推导过程。与其他文章不同的是，我们将从牛顿法这一经典优化算法入手，推导出几种广泛使用的变体，然后再回到深度学习（DL）优化器。我们的主要目标是将 DL 纳入这些经典方法的范畴。虽然我们将专注于 DL，但我们也会在本书稍后部分重温经典算法，以改进学习算法。物理模拟夸大了神经网络造成的困难，这就是为什么下面的主题与基于物理的学习任务特别相关。

2.8.1 符号

本章使用一种自定义符号，这种符号经过精心挑选，可以清晰简要地表示所考虑的所有方法。我们有一个标量损失函数 $L(x)：{\mathbb{R}}^N\to \mathbb{R}$，最佳值（$L$ 的最小值）位于位置 $x^{\ast }$，$\Delta$ 表示在 $x$ 中的一个步。$x$中不同的中间更新步用下标表示，如$x_n$或$x_k$。

本章使用经过仔细选择的自定义符号， 以清晰简洁地表示正在考虑的所有方法。我们有一个标量损失函数 $L(x):{\mathbb{R}}^N\to \mathbb{R}$，最优点 (L 的最小值)在 $x^{\ast }$ 处， $\Delta$ 表示 $x$ 中的一步。 $x$ 中的不同中间更新步骤用下标表示， 例如 $x_n$ 或 $x_k$。

在下文中，我们经常需要取反，即除以某个量。对于矩阵 $A$ 和 $B$，我们定义 $\frac{A}{B}\equiv B^{-1}A$。当 $a$ 和 $b$ 是向量时，结果就是用下面两种公式之一得到的矩阵。我们将具体说明使用哪一种： $$\frac{a}{b}\equiv \frac{a{a}^T}{a^{T}b}\text{ 或 }\frac{a}{b}\equiv \frac{a{b}^T}{b^{T}b}$$

将 $/partial//partialx$ 应用于 $L$ 一次，可以得到雅各布值 $J(x)$。由于 $L$ 是标量，$J$ 是行向量，梯度（列向量）$\nabla L$ 由 $J^T$ 给出。再次应用 $\partial /\partial x$ 可以得到 Hessian 矩阵 $H(x)$，再次应用 $\partial /\partial x$ 可以得到第三导数张量，用 $K(x)$ 表示。幸运的是，我们从来不需要把 $K$ 作为一个完整的张量来计算，但在下面的一些推导中需要用到它。为了缩短下面的符号，当函数或导数在 $x$ 处求值时，我们通常会去掉 $(x)$，例如，$J$ 将表示为 $J(x)$。

下图概述了一些常用量的矩阵形状。之后我们就不需要它了，但在本图中，$N$表示$x$的维度，即$x\in {\mathbb{R}}^N$。

2.8.2 准备工作

我们需要一些工具来进行下面的推导，现将其总结如下，以供参考。

毫不奇怪，我们需要一些泰勒级数展开式。用上面的符号表示为 $$L(x+\Delta )=L+J\Delta +\frac{1}{2}H{\Delta }^2+\cdots$$

然后，我们还需要拉格朗日形式，它可以从区间 $[x,x+\Delta ]$ 中得到 $\xi$ 的精确解：

$$L(x+\Delta )=L+J\Delta +\frac{1}{2}H(\xi ){\Delta }^2$$

在一些情况下，我们会用到微积分基本定理，为了完整起见，在此重复一遍：

$$f(x+\Delta )=f(x)+{\int }_0^1\text{d}s{f}^{\prime }(x+s\Delta )\Delta .$$

此外，我们还将利用具有常数 $\mathcal{L}$ 的 Lipschitz连续性： $|f(x+\Delta )+f(x)|\le \mathcal{L}\Delta$，以及著名的 Cauchy-Schwartz 不等式：$u^{T}v<|u|\cdot |v|$。

2.8.3 牛顿法

现在，我们可以从最经典的优化算法开始：牛顿法它是通过将我们感兴趣的函数近似为抛物线而得出的。这是因为几乎所有的最小值近看都像抛物线。

因此，我们可以用抛物线的形式 $\frac{1}{2}H(x-x^{\ast }{)}^2+c$ 来表示最佳值 $x^{\ast }$ 附近的 $L$ 。其中$c$表示常数偏移。在位置 $x$ 处，我们观察到 $H$ 和 $J^T=H\cdot (x_k-x^{\ast })$ 。将其重新排列后，可以直接得出计算最小值的方程：$x^{\ast }=x_k-\frac{J^T}{H}$。牛顿法默认以单步计算 $x^{\ast }$ ，因此牛顿法中的 $x$ 更新量给定为： $$\Delta =-\frac{J^T}{H}$$

让我们来看看牛顿法的收敛阶次。对于 $L$ 的最优值 $x^{\ast }$，让 ${\Delta }_n^{\ast }=x^{\ast }-x_n$ 表示从当前的 $x_n$ 到最优值的步长，如下图所示。

假设 $J$ 可微分，我们可以在 $x^{\ast }$ 处对 $J^T$ 进行拉格朗日展开：

$$\begin{array}{rl}0& =J^T(x^{\ast })=J^T(x_n)+H(x_n){\Delta }_n^{\ast }+\frac{1}{2}K({\xi }_n){{\Delta }_n^{\ast }}^2\\ \frac{J^T}{H}& =-\frac{K}{2H}{{\Delta }_n^{\ast }}^2-{\Delta }_n^{\ast }\end{array}$$

在第二行中，我们已经除以 $H$，并以缩短符号去掉了 $(x_n)$ 和 $({\xi }_n)$ 。当我们把这一行插入 ${\Delta }_n^{\ast }$ 时，我们会得到 $$\begin{array}{rl}{\Delta }_{n+1}^{\ast }& =x^{\ast }-x_{n+1}\\ & =x^{\ast }-(x_n-\frac{J^T}{H})\\ & ={\Delta }_n^{\ast }-\frac{K}{2H}{{\Delta }_n^{\ast }}^2-{\Delta }_n^{\ast }\\ & =-\frac{K}{2H}{{\Delta }_n^{\ast }}^2\end{array}$$

因此，与最优值的距离会因 ${{\Delta }_n^{\ast }}^2$ 的变化而改变，这意味着一旦我们足够接近最优值，就会出现二次收敛。这当然很好，但它仍然取决于前系数 $\frac{K}{2H}$，如果前系数 $>1$ 就会发散。请注意，这是一个精确的表达式，由于拉格朗日展开的缘故，没有截断。到目前为止，我们已经实现了二次收敛，但并不能保证收敛到最优。为此，我们必须允许步长可变。

2.8.4 自适应步长

因此，作为牛顿法的下一步，我们引入了步长可变的$\lambda$，从而得到迭代$x_{n+1}=x_n+\lambda \Delta =x_n-\lambda \frac{J^T}{H}$ 。如下图所示，如果 $L$ 不完全是抛物线，而且小的 $H$ 可能会以不理想的方式超调，那么这一点就特别有用。本例中最左侧的例子：

为了说明收敛性，我们需要一些基本假设：损失函数的凸性和平滑性。然后，我们将专注于证明损失会减少，并且我们会沿着一连串损失值较低的较小集合 $\forall xL(x)<L(x_n)$ 移动。

首先， 我们对 $L$ 应用基本定理：

$$L(x+\lambda \Delta )=L+{\int }_0^1\text{ds}J(x+s\lambda \Delta )\lambda \Delta ,$$

类似地， 使用它围绕此位置表达 $J$：

$$\begin{array}{r}J(x+s\lambda \Delta )=J+s\lambda {\Delta }^T{\int }_0^1\text{dt}H(x+st\lambda \Delta )\end{array}$$

将这个 $J^T$ 插入 $L$ 中会产生：

$$\begin{array}{rl}L(x+\lambda \Delta )-L(x)& ={\int }_0^1\text{d}s[J+s\lambda {\Delta }^T{\int }_0^1\text{d}tH(x+st\lambda \Delta ){]}^T\lambda \Delta \\ & ={\int }_0^1\text{d}s[-H\Delta +{\int }_0^1\text{d}tH(x+st\lambda \Delta )s\lambda \Delta {]}^T\lambda \Delta \\ & ={\int }_0^1\text{d}s{\int }_0^1\text{d}t[-H\Delta +H(x+st\lambda \Delta )s\lambda \Delta {]}^T\lambda \Delta \\ & ={\int }_0^1\text{d}s{\int }_0^1\text{d}t[-H\Delta (1+\lambda s)+[H(x+st\lambda \Delta )-H]s\lambda \Delta {]}^T\lambda \Delta \\ & =-(H\Delta {)}^T(1+\frac{\lambda }{2})\lambda \Delta +{\lambda }^2{\int }_0^1\text{d}ss{\int }_0^1\text{d}t[[H(x+st\lambda \Delta )-H]\Delta {]}^T\Delta \end{array}$$

在这里，我们首先使用了 $\Delta =-J^T/H$ ，然后把它移到第三行的积分中，与 $s\lambda \Delta$ 项一起。第四行将 $H\Delta$ 因数化。这样，我们就可以对上式（10）最后一行中的第一项 $-H\Delta (1+\lambda s)$ 进行积分。

在接下来的一系列步骤中，我们将首先使用 $(H\Delta {)}^T\Delta =\Vert H^{\frac{1}{2}}\Delta {\Vert }^2$。这一术语将简称为 $epsilon\equiv \Vert H^{\frac{1}{2}}\Delta {\Vert }^2$。由于 $H$ 的特性，这个 $epsilon$ "只是 "代表一个小的正因子，它会一直存在到最后。之后，在下面的第 3 行和第 4 行中，我们可以开始寻找损失变化的上限。我们首先使用考希-施瓦茨不等式，然后利用仿射共轭矩阵的特殊 Lipschitz 条件。对于 $H$，它的形式是 $\Vert H(x{)}^{-\frac{1}{2}}(H(y)-H(x))\Vert \le \mathcal{L}\Vert H(x{)}^{\frac{1}{2}}(y-x){\Vert }^2$。这就要求 $H$ 是对称的正定式，这在实际中并不是太不合理。承上启下，我们可以得到： $$\begin{array}{rl}\cdots & =-(\lambda +\frac{{\lambda }^2}{2})\cdot \Vert H^{\frac{1}{2}}\Delta {\Vert }^2+{\lambda }^2{\int }_0^1\text{d}ss{\int }_0^1\text{d}t[H(x+st\lambda \Delta )-H]\Delta {]}^T\Delta \\ & =-\lambda ϵ+\frac{{\lambda }^2ϵ}{2}+{\lambda }^2{\int }_0^1\text{d}ss{\int }_0^1\text{d}t[H^{-\frac{1}{2}}[H(x+st\lambda \Delta )-H]\Delta {]}^T{H}^{\frac{1}{2}}\Delta \\ & \le -\lambda ϵ+\frac{{\lambda }^2ϵ}{2}+{\lambda }^2{\int }_0^1\text{d}ss{\int }_0^1\text{d}t\Vert H^{-\frac{1}{2}}[H(x+st\lambda \Delta )-H]\Delta \Vert \cdot \Vert H^{\frac{1}{2}}\Delta \Vert \\ & \le -\lambda ϵ+\frac{{\lambda }^2ϵ}{2}+{\lambda }^3{\int }_0^1\text{d}s{s}^2{\int }_0^1\text{d}tt\mathcal{L}\Vert H^{\frac{1}{2}}\Delta {\Vert }^2\cdot \Vert H^{\frac{1}{2}}\Delta \Vert \\ & =-\lambda ϵ+\frac{{\lambda }^2ϵ}{2}+\mathcal{L}{\lambda }^3{\int }_0^1\text{d}s{s}^2{\int }_0^1\text{d}tt\Vert H^{\frac{1}{2}}\Delta {\Vert }^3\\ & =-\lambda ϵ+\frac{{\lambda }^2ϵ}{2}+\mathcal{L}{\lambda }^3{\int }_0^1\text{d}s{s}^2{\int }_0^1\text{d}tt{ϵ}^{\frac{3}{2}}\\ & =-\lambda ϵ+\frac{{\lambda }^2ϵ}{2}+\frac{{\lambda }^3\mathcal{L}{ϵ}^{\frac{3}{2}}}{6}\end{array}$$

由于 $H^T=H$，我们把 $H^{-\frac{1}{2}}$ 移到了第 2 行的积分内。在第 4 行中，除了应用特殊的 Lipschitz 条件外，我们尽可能地将 $s,t$ 和 $\lambda$ 从积分中移出。最后三行只是简化术语，用 $ϵ$ 表示 $H$ 的出现，并对积分进行评估。这样，我们就得到了以步长 $\lambda$ 为单位的三次方形式。最重要的是，第一个线性项是负的，因此对于较小的 $\lambda$ 将占主导地位。这样，我们就证明了对于足够小的 $\lambda$ 来说，步长将是负值： $$L(x+\lambda \Delta )=L(x)-\lambda ϵ+\frac{{\lambda }^2ϵ}{2}+\frac{{\lambda }^3\mathcal{L}{ϵ}^{\frac{3}{2}}}{6}$$

然而，这本质上要求我们自由选择 $\lambda$，因此该证明不适用于上述固定步长。它的一个很好的特性是，我们 "只 "要求 $H$ 的 Lipschitz 连续性，而不是 $J$ 或甚至 $L$。

总之，我们已经证明，牛顿方法的自适应步长可以收敛，这非常好。然而，它需要 Hessian $H$ 作为核心要素。不幸的是，$H$ 在实践中很难获得。这是一个真正的拦路虎，也是以下方法的动因。这些方法保留了牛顿法的基本步骤，但对 $H$ 进行了近似。

2.8.5 近似 Hessian

接下来， 我们将重新审视牛顿法导出的三种流行算法。

2.8.6 Broyden 方法

近似 $H$ 的第一种方法是做一个非常粗略的猜测，即从特征矩阵 $H_0=1$ 开始，然后通过有限差分近似迭代更新 $H_n$。对于布洛伊登方法，我们使用向量除法 $\frac{a}{b}\frac{a{b}^T}{b^{T}b}$。

为了简化有限差分，我们另外假设在当前位置 $x_n$ 时，已经达到 $J(x_n)=0$。当然，这并不一定正确，但在优化过程中，我们可以用下面这个简洁的表达式来修改 $H$：

$$H_{n+1}=H_n+\frac{J(x_{n+1}{)}^T}{\Delta }$$

和之前一样，我们对 $x$ 使用 $\Delta =-\frac{J^T}{H}$ 的步长，分母来自有限差分 $\frac{J(x_{n+1}{)}^T-J(x_n{)}^T}{\Delta }$，假设当前的雅各布值为零。请记住，这里的 $\Delta$ 是一个向量（见上文的向量分割），所以有限差分给出了一个大小为 $N\times N$ 的矩阵，可以添加到 $H_n$。

Broyden 方法的优势在于我们无需计算全Hessian，而且可以高效地评估$H$的更新。然而，上述两个假设使其成为一个非常粗糙的近似值，因此通过 $\frac{J^T}{H}$ 中的逆Hessian对更新步骤进行归一化也相应地不可靠。

2.8.7 BFGS

这就产生了 BFGS 算法（以_Broyden-Fletcher-Goldfarb-Shanno_命名），它引入了一些重要的改进：它不假定 $J$ 立即为零，并对更新的冗余部分进行补偿。这是必要的，因为有限差分 $(J(x_{n+1})-J(x_n))/{\Delta }_n$ 给出了 $H$ 的完整近似值。我们可以尝试执行某种平均程序，但这会严重破坏 $H_n$ 中的现有内容。因此，我们只沿着当前步长 $\Delta$ 减去 $H_n$ 中的现有条目。这样做是有道理的，因为有限差分近似正好可以得到沿 $\Delta$ 的估计值。结合使用向量分割 $\frac{a}{b}\frac{a{a}^T}{a^{T}b}$，这些变化给出了 $H$ 的更新步长： $$H_{n+1}=H_n+\frac{J(x_{n+1}{)}^T-J(x_n{)}^T}{{\Delta }_n}-\frac{H_n{\Delta }_n}{{\Delta }_n}$$

在实践中，BFGS 还利用直线搜索来确定步长 $\lambda$。由于 $H$ 的大小较大，通常使用的 BFGS 变体也会使用 $H$ 的缩小表示来节省内存。尽管如此，通过沿搜索方向的有限差分近似更新Hessian矩阵的核心步骤仍是 BFGS 的核心，并使其至少能部分补偿损失景观的缩放效应。目前，BFGS 变体是解决经典非线性优化问题最广泛使用的算法。

2.8.8 高斯牛顿法

通过限制 $L$ 为经典的 $L^2$ 损失，可以推导出牛顿法的另一个有吸引力的变体。这就是_高斯-牛顿（GN）算法。因此，我们仍然使用 $\Delta =-\frac{J^T}{H}$ ，但对于任意的 $f(x)$，则依赖于 $L=|f{|}^2$ 形式的平方损失。$f$的导数用$J_f,H_f$表示，而不是像之前那样用$J,H$表示$L$。由于链式规则，我们有 $J=2f^T{J}_f$。

二阶导数得出以下表达式。为了简化 GN，我们省略了下面第二行中的二阶项： $$\begin{array}{rl}H& =2J_f^T{J}_f+2f^T{H}_f\\ & \approx 2J_f^T{J}_f\\ & =2\frac{J^{T}J}{4|f{|}^2}\\ & =\frac{J^{T}J}{2L}\end{array}\text{\hspace{1em}(13)}$$

在这里，一阶近似的剩余 $J_f^T{J}_f$ 项可以简化，这要归功于我们对 $L^2$ 损失的关注：$J_f=J/(2f^T)$ 和 $|f{|}^2=L$。

等式 (13) 的最后一行意味着我们基本上是用 $J$ 的平方来逼近 Hessian。这在很多情况下都是合理的，将其插入到我们上面的更新步骤中，就得到了高斯-牛顿更新 ${\Delta }_{\text{GN}}\approx -\frac{J^T}{J^{T}J}$。

看这个更新， 它本质上采用了形式 ${\Delta }_{\text{GN}}\approx -\frac{1}{J}$ 的步骤， 即更新基于 $L$ 的雅可比的逼近逆。这为所有参数提供了大致相等的步长， 作为这样的提供了我们将在本书后面重新审视的有趣构建块。以上形式意味着我们仍然必须求一个大矩阵的逆， 这代价高昂， 且矩阵本身可能甚至不可逆。

2.8.9 回到深度学习

上文我们已经说明，牛顿方法是许多流行的非线性优化算法的核心。现在，高斯-牛顿法终于为我们提供了通向深度学习算法的阶梯，特别是通向亚当优化器的阶梯。

2.8.10 Adam

像往常一样，我们从牛顿步骤 $\Delta =-\lambda \frac{J^T}{H}$ 开始，但即使是高斯-牛顿对 $H\approx J^{T}J$ 的最简单近似，也需要求一个潜在的巨大矩阵的逆。这对于神经网络的权重来说是不可行的，因此一个合理的问题是，我们如何进一步简化这一步骤？对于Adam来说，答案就是：用对角线近似。具体来说，Adam使用： $$H\approx \sqrt{\text{diag}(J^{T}J)}$$

这只是真正的Hessian的一个非常粗略的近似值。我们在这里使用的只是平方一阶导数，当然，一般来说，${\left(\frac{\partial f}{\partial x}\right)}^2\ne \frac{{\partial }^{2}f}{\partial x^2}.$。这只适用于高斯-牛顿的一阶近似，即方程 (13) 的第一项。现在，亚当更进一步，只保留 $J^{T}J$ 的对角线。这个量在深度学习中以权重梯度的形式随时可用，并使 $H$ 的反演变得微不足道。因此，它至少提供了对单个权重曲率的一些估计，但忽略了它们之间的相关性。

有趣的是，亚当并没有通过 $\text{diag}(J^{T}J)$ 进行完全反转，而是使用了分量平方根。这就有效地得到了 $\sqrt{\text{diag}(J^{T}J)}\approx \sqrt{\text{diag}(J^2)}\approx \text{diag}(J)$。因此，亚当沿着 $\Delta \approx \text{sign}(-J)$ 移动，近似于在所有维度上执行一个固定大小的步长。

在实践中，Adam还引入了一些技巧，即计算梯度 $J^T$ 以及梯度平方动量，在优化迭代过程中平均这两个量。这使得估算结果更加稳健，这一点至关重要：如果梯度输入错误地过小，那么归一化可能会导致爆炸。Adam 还在除法时增加了一个小常数，而平方根同样有助于减轻超调。

总而言之：Adam 利用一阶更新和对角高斯-牛顿近似赫塞斯进行归一化。此外，它还利用动量进行稳定。

2.8.11 梯度下降

为了实现梯度下降（GD）优化，我们现在采取最后一步，在 $-\lambda \frac{J^T}{H}$ 中假设 $H=1$。这样，我们就有了一个由缩放梯度 $\Delta =-\lambda J^T$ 组成的更新。

有趣的是，在没有任何形式的反转的情况下，$J^T$ 本身具有添加到 $x$ 的正确形状，但它 "生活 "在错误的空间 $1/x$ 中，而不是 $x$ 本身。上述所有 "经典 "方案以某种形式使用的 $J^{-1}$ 与 $J^T$ 本质上是不同的，这一事实就说明了这一点。缺乏反演步骤的问题在涉及物理学的问题中尤为突出，因此也是本章的重要启示之一。我们将在以后的 Scale-Invariance and Inversion中更详细地讨论这个问题。现在，我们建议大家牢记不同算法之间的关系，即使 "仅仅 "使用 Adam 或 GD。

最后，让我们来看看 GD 的收敛性。我们再次假设 $L$ 具有凸性和可微性。用泰勒级数展开 $\Delta$ 的一步，用 Lipschitz 条件$J(x+\Delta )-J(x)\le \mathcal{L}\Delta$限定二阶项，我们得到

$$\begin{array}{rl}L(x+\Delta )-L(x)& \le J(x)\Delta +\frac{\mathcal{L}}{2}|\Delta {|}^2\\ & =J(-\lambda J^T)+\frac{\mathcal{L}}{2}|-\lambda J^T{|}^2\\ & =-\lambda J{J}^T+{\lambda }^2\frac{\mathcal{L}}{2}|J{|}^2\\ & =-\lambda |J{|}^2+{\lambda }^2\frac{\mathcal{L}}{2}|J{|}^2\end{array}$$

就像上面方程(11)中的牛顿方法一样，我们有一个负线性项，它在足够小的$\lambda$时支配着损失。结合起来，由于第一行中的 Lipschitz 条件，我们可以得到以下上界 $L(x+\Delta )\le L(x)-\lambda |J{|}^2+\frac{{\lambda }^2\mathcal{L}}{2}.|J{|}^2$.通过选择 $\lambda \le \frac{1}{\mathcal{L}}$，我们可以进一步简化这些项，并得到一个取决于 $J$ 平方的上界：$L(x+\Delta )\le L(x)-\frac{\lambda }{2}.|J{|}^2$，从而确保收敛。

遗憾的是，这一结果在实践中对我们帮助不大，因为在深度学习中，梯度的所有常见用法 $\mathcal{L}$ 都是未知的。不过，我们还是应该知道，梯度的 Lipschitz 常数在理论上可以为 GD 提供收敛保证。

至此，我们结束了经典优化器及其与深度学习方法的关系之旅。值得注意的是，为了清楚起见，我们在这里重点讨论了非随机算法，因为随机算法的证明会变得更加复杂。

有了这些背景知识，现在正是开始研究一些实际例子的好时机，这些例子从尽可能简单的开始，采用完全有监督的方法进行训练。

这里的监督（Supervised）主要是指："用老办法做事"。当然，在深度学习（DL）的语境中，"老套"还是相当新的。此外，"老套"并不总是意味着不好--我们稍后将讨论训练网络的方法，这些方法明显优于使用监督训练的方法。

尽管如此，"监督训练 "是人们在 DL 背景下遇到的所有项目的起点，因此值得研究。此外，虽然 "监督训练 "的结果通常不如基于物理学的方法，但在某些应用场景中，如果没有好的模型方程，"监督训练 "可能是唯一的选择。

3.1 问题设置

在监督训练中，我们面对一个未知函数 $f^{\ast }(x)=y^{\ast }$，收集大量数据对 $[x_0,y_0^{\ast }]，...[x_n,y_n^{\ast }]$（训练数据集），并直接训练一个 NN 来表示 $f^{\ast }$ 的近似值 $f$。

我们用这种方法得到的 $f$ 通常不是精确的，而是通过最小化问题得到的：通过调整表示 $f$ 的 NN 的权重 $\theta$，使得 $${\text{arg min}}_{\theta }\sum _i(f(x_i;\theta )-y_i^{\ast }{)}^2.$$

这样我们就能得到 $\theta$，从而在我们选择 $f$ 和超参数进行训练的情况下，尽可能精确地得到 $f(x;\theta )=y\approx y^{\ast }$ 的结果。请注意，上面我们假设的是最简单的 $L^2$ 损失。更一般的版本会在损失 $L$ 中使用一个误差度量 $e(x,y)$，通过 ${\text{arg min}}_{\theta }{\sum }_{i}e(f(x_i;\theta ),y_i^{\ast }))$ 来最小化。关于合适度量的选择，我们稍后再讨论。

无论我们选择哪种度量方法，这一表述都给出了监督式方法的实际 "学习 "过程。

训练数据通常需要相当大的规模，因此使用数值模拟求解物理模型 $\mathcal{P}$ 来产生大量可靠的输入输出对进行训练是很有吸引力的。这意味着训练过程使用一组模型方程，并对其进行数值近似，以训练 NN 表示 $f$。这样做有很多好处，例如，我们没有真实世界设备的测量噪声，也不需要人工标注大量样本来获取训练数据。

另一方面，这种方法也继承了以模拟代替实验的常见挑战：首先，我们需要确保所选模型有足够的能力来预测我们感兴趣的真实世界现象的行为。此外，数值近似会产生数值误差，需要将误差控制得足够小，以满足所选应用的需要。由于这些问题在经典模拟中都有深入研究，因此同样可以利用现有知识来设置 DL 训练任务。

3.2 代理模型

上述监督方法的核心优势之一是我们可以得到一个 代理模型(Surrogate Model)，即一个模仿原始 $\mathcal{P}$ 行为的新函数。针对现实世界现象的 PDE 模型的数值近似通常计算成本很高。另一方面，训练有素的 NN 每次评估的成本是固定的，而且在 GPU 或 NN 单元等专用硬件上进行评估通常是微不足道的。

尽管如此，还是要小心谨慎：神经网络可以迅速产生大量介于两者之间的结果。考虑一个具有 $128$ 特征的 CNN 层。如果我们将其应用于 $128^2$ 的输入，即大约 16k 个单元，我们会得到 $128^3$ 的中间值。这就超过了 200 万。所有这些值至少需要暂时存储在内存中，并由下一层进行处理。

尽管如此， 我们仍需谨慎: NN 可以快速生成大量中间结果。考虑一个具有 $128$ 个特征的 CNN 层。如果我们将其应用于 $128^2$，即大约 16 k 个单元的输入， 我们会得到 $128^3$ 个中间值。这已经超过了 200 万。这些值至少需要暂时存储在内存中，并由下一层处理。

然而，用快速、可学习的近似方法取代复杂、昂贵的求解器是一个非常有吸引力和有趣的方向。

3.3 给我看看代码!

最后，让我们来看一个训练神经网络的代码示例：我们将用[TWPH20]中的代用模型替换机翼周围湍流的完整求解器。

3.4 翼型周围 RANS 流动的监督训练

3.4.1 概述

在这个监督训练的例子中，我们有一个翼面周围的湍流气流，我们想知道在不同雷诺数和攻角下，该翼面周围的平均运动和压力分布。因此，在给定机翼形状、雷诺数和攻角的情况下，我们希望获得机翼周围的速度场和压力场。

经典的近似方法是_雷诺兹平均纳维-斯托克斯（RANS）模型，这种设置仍然是纳维-斯托克斯求解器在工业中最广泛的应用之一。不过，我们现在的目标不是依靠传统的数值方法来求解 RANS 方程，而是通过神经网络训练一个完全绕过数值求解器的代理模型，并以速度和压力的形式求解。

run in colab

3.4.2 准备工作

根据监督模型的描述， 我们的学习任务非常简单直接， 可以写成

$$\begin{array}{r}{\text{arg min}}_{\theta }\sum _i(f(x_i;\theta )-y_i^{\ast }{)}^2，\end{array}$$ 其中，$x$ 和 $y^{\ast }$ 分别由一组物理场组成，而索引 $i$ 则表示数据集所有离散点的差值。

我们的目标是在围绕翼面的计算域 $\Omega$ 中推断速度 $\mathbf{u}=u_x,u_y$ 和压力场 $p$。作为输入，我们有雷诺数$\text{Re}$。在$\mathbb{R}$ 中的雷诺数 $\text{Re}$、攻角 $\alpha \in \mathbb{R}$ 中的攻角 $\alpha$，以及机翼形状 $\mathbf{s}$，这些数据以 $128^2$ 的栅格编码。$\text{Re}$ 和 $\alpha$ 是以自由流流速 $\mathbf{f}$ 的形式提供的，其 x 和 y 分量表示为相同大小的常数域，在机翼区域包含零。因此，输入和输出具有相同的尺寸：$x,y^{\ast }\in {\mathbb{R}}^{3\times 128\times 128}$。输入包含 $[f_x,f_y,\text{mask}]$，而输出则存储通道 $[p,u_x,u_y]$。这正是我们下面要为 NN 指定的输入和输出维度。

这里需要注意的一点是，我们在 $y^{\ast }$ 中关注的量包含三个不同的物理场。虽然两个速度分量在本质上非常相似，但压力场通常具有与速度近似平方缩放的不同行为（参见 伯努利）。这意味着我们需要小心处理简单求和（如上文的最小化问题），并注意对数据进行归一化处理。

3.4.3 代码即将出现...

让我们开始实施吧。请注意，我们将跳过数据生成过程。下面的代码改变this codebase。在这里，我们只需下载一小部分在OpenFOAM中通过 Spalart-Almaras RANS 仿真生成的训练数据。

import numpy as np
import os.path, random
import torch
from torch.utils.data import Dataset
print("Torch version {}".format(torch.__version__))

# get training data
dir = "./"
if True:
    # download
    if not os.path.isfile('data-airfoils.npz'):
        import requests
        print("Downloading training data (300MB), this can take a few minutes the first time...")
        with open("data-airfoils.npz", 'wb') as datafile:
            resp = requests.get('https://dataserv.ub.tum.de/s/m1615239/download?path=%2F&files=dfp-data-400.npz', verify=False)
            datafile.write(resp.content)
else: 
    # alternative: load from google drive (upload there beforehand):
    from google.colab import drive
    drive.mount('/content/gdrive')
    dir = "./gdrive/My Drive/"

npfile=np.load(dir+'data-airfoils.npz')    
print("Loaded data, {} training, {} validation samples".format(len(npfile["inputs"]),len(npfile["vinputs"])))
print("Size of the inputs array: "+format(npfile["inputs"].shape))

Torch version 2.0.1+cu118
Downloading training data (300MB), this can take a few minutes the first time...
Loaded data, 320 training, 80 validation samples
Size of the inputs array: (320, 3, 128, 128)

如果您在 colab 中运行本笔记本，上面的 else语句（默认情况下未激活）可能会对您有帮助：与其每次都重新下载训练数据，您可以手动下载一次并将其存储在谷歌硬盘中。我们假设它存储在根目录下，名为 data-airfoils.npz。之后，您可以使用上面的代码从谷歌硬盘加载文件，这样通常会快很多。如果您想通过 colab 进行更广泛的实验，强烈建议使用此方法。

3.4.4 RANS 训练数据

现在我们有了一些训练数据。一般来说，尽可能多地了解我们正在使用的数据是非常重要的（对于任何 ML 任务来说，garbage-in-gargabe-out 原则绝对成立）。我们至少应该从维度和粗略统计的角度来理解数据，最好还能从内容的角度来理解数据。否则，我们将很难解释训练运行的结果。而且，尽管有所有的 DL 魔法：如果你无法在数据中找出任何模式，NNs 肯定也找不到任何有用的模式。

因此，让我们来看看其中一个训练样本...下面只是一些辅助代码，用于并排显示图像。

import pylab

# helper to show three target channels: normalized, with colormap, side by side
def showSbs(a1,a2, stats=False, bottom="NN Output", top="Reference", title=None): 
  c=[]
  for i in range(3):
    b = np.flipud( np.concatenate((a2[i],a1[i]),axis=1).transpose())
    min, mean, max = np.min(b), np.mean(b), np.max(b); 
    if stats: print("Stats %d: "%i + format([min,mean,max]))
    b -= min; b /= (max-min)
    c.append(b)
  fig, axes = pylab.subplots(1, 1, figsize=(16, 5))
  axes.set_xticks([]); axes.set_yticks([]); 
  im = axes.imshow(np.concatenate(c,axis=1), origin='upper', cmap='magma')

  pylab.colorbar(im); pylab.xlabel('p, ux, uy'); pylab.ylabel('%s           %s'%(bottom,top))
  if title is not None: pylab.title(title)

NUM=72
showSbs(npfile["inputs"][NUM],npfile["targets"][NUM], stats=False, bottom="Target Output", top="Inputs", title="3 inputs are shown at the top (free-ux, free-uy, mask), with the 3 output channels (p,ux,uy) at the bottom")

接下来，让我们定义一个小的辅助类 DfpDataset 来组织输入和目标。我们将把相应的数据传输到 pytorch 的 DataLoader类。

我们还设置了一些全局参数来控制训练参数，其中最重要的可能是：学习率 LR，即从上一次训练中获得的 $\eta$。当训练运行没有收敛时，这是第一个需要实验的参数。

在这里，我们将始终保持相对较小的学习率。(使用学习率衰减会更好，即有可能改善收敛性，但为了清晰起见，此处省略）。

# some global training constants

# number of training epochs
EPOCHS = 100
# batch size
BATCH_SIZE = 10
# learning rate
LR = 0.00002

class DfpDataset():
    def __init__(self, inputs,targets): 
        self.inputs  = inputs
        self.targets = targets

    def __len__(self):
        return len(self.inputs)

    def __getitem__(self, idx):
        return self.inputs[idx], self.targets[idx]

tdata = DfpDataset(npfile["inputs"],npfile["targets"])
vdata = DfpDataset(npfile["vinputs"],npfile["vtargets"])

trainLoader = torch.utils.data.DataLoader(tdata, batch_size=BATCH_SIZE, shuffle=True , drop_last=True) 
valiLoader  = torch.utils.data.DataLoader(vdata, batch_size=BATCH_SIZE, shuffle=False, drop_last=True) 

print("Training & validation batches: {} , {}".format(len(trainLoader),len(valiLoader) ))

输出结果为：

Training & validation batches: 32 , 8

3.4.5 网络设置

现在我们可以设置神经网络的架构，我们将使用完全卷积 U-net。这是一种广泛使用的架构，使用不同空间分辨率的卷积堆叠。与普通卷积网络的主要区别在于，从编码器到解码器部分引入了skip connection。这可以确保在特征提取过程中不会丢失任何信息。(请注意，这只有在网络作为一个整体使用时才有效。例如，如果我们想将解码器作为一个独立的组件来使用，这就行不通了）。

下面是体系结构概述:

首先，我们将定义一个辅助程序，用于在网络中设置一个卷积块，即 blockUNet 。注意，我们不使用任何池化！相反，我们按照 最佳实践，使用跨距和转置卷积（解码器部分需要对称，即内核大小不均）。完整的 pytroch 神经网络通过 DfpNet 类进行管理。

import os, sys, random
import numpy as np

import torch
import torch.nn as nn
import torch.optim as optim
import torch.autograd 
import torch.utils.data 

def blockUNet(in_c, out_c, name, size=4, pad=1, transposed=False, bn=True, activation=True, relu=True, dropout=0. ):
    block = nn.Sequential()

    if not transposed:
        block.add_module('%s_conv' % name, nn.Conv2d(in_c, out_c, kernel_size=size, stride=2, padding=pad, bias=True))
    else:
        block.add_module('%s_upsam' % name, nn.Upsample(scale_factor=2, mode='bilinear'))
        # reduce kernel size by one for the upsampling (ie decoder part)
        block.add_module('%s_tconv' % name, nn.Conv2d(in_c, out_c, kernel_size=(size-1), stride=1, padding=pad, bias=True))

    if bn:
        block.add_module('%s_bn' % name, nn.BatchNorm2d(out_c))
    if dropout>0.:
        block.add_module('%s_dropout' % name, nn.Dropout2d( dropout, inplace=True))

    if activation:
        if relu:
            block.add_module('%s_relu' % name, nn.ReLU(inplace=True))
        else:
            block.add_module('%s_leakyrelu' % name, nn.LeakyReLU(0.2, inplace=True))

    return block
    
class DfpNet(nn.Module):
    def __init__(self, channelExponent=6, dropout=0.):
        super(DfpNet, self).__init__()
        channels = int(2 ** channelExponent + 0.5)

        self.layer1 = blockUNet(3         , channels*1, 'enc_layer1', transposed=False, bn=True, relu=False, dropout=dropout )
        self.layer2 = blockUNet(channels  , channels*2, 'enc_layer2', transposed=False, bn=True, relu=False, dropout=dropout )
        self.layer3 = blockUNet(channels*2, channels*2, 'enc_layer3', transposed=False, bn=True, relu=False, dropout=dropout )
        self.layer4 = blockUNet(channels*2, channels*4, 'enc_layer4', transposed=False, bn=True, relu=False, dropout=dropout )
        self.layer5 = blockUNet(channels*4, channels*8, 'enc_layer5', transposed=False, bn=True, relu=False, dropout=dropout ) 
        self.layer6 = blockUNet(channels*8, channels*8, 'enc_layer6', transposed=False, bn=True, relu=False, dropout=dropout , size=2,pad=0)
        self.layer7 = blockUNet(channels*8, channels*8, 'enc_layer7', transposed=False, bn=True, relu=False, dropout=dropout , size=2,pad=0)
     
        # note, kernel size is internally reduced by one for the decoder part
        self.dlayer7 = blockUNet(channels*8, channels*8, 'dec_layer7', transposed=True, bn=True, relu=True, dropout=dropout , size=2,pad=0)
        self.dlayer6 = blockUNet(channels*16,channels*8, 'dec_layer6', transposed=True, bn=True, relu=True, dropout=dropout , size=2,pad=0)
        self.dlayer5 = blockUNet(channels*16,channels*4, 'dec_layer5', transposed=True, bn=True, relu=True, dropout=dropout ) 
        self.dlayer4 = blockUNet(channels*8, channels*2, 'dec_layer4', transposed=True, bn=True, relu=True, dropout=dropout )
        self.dlayer3 = blockUNet(channels*4, channels*2, 'dec_layer3', transposed=True, bn=True, relu=True, dropout=dropout )
        self.dlayer2 = blockUNet(channels*4, channels  , 'dec_layer2', transposed=True, bn=True, relu=True, dropout=dropout )
        self.dlayer1 = blockUNet(channels*2, 3         , 'dec_layer1', transposed=True, bn=False, activation=False, dropout=dropout )

    def forward(self, x):
        # note, this Unet stack could be allocated with a loop, of course... 
        out1 = self.layer1(x)
        out2 = self.layer2(out1)
        out3 = self.layer3(out2)
        out4 = self.layer4(out3)
        out5 = self.layer5(out4)
        out6 = self.layer6(out5)
        out7 = self.layer7(out6)
        # ... bottleneck ...
        dout6 = self.dlayer7(out7)
        dout6_out6 = torch.cat([dout6, out6], 1)
        dout6 = self.dlayer6(dout6_out6)
        dout6_out5 = torch.cat([dout6, out5], 1)
        dout5 = self.dlayer5(dout6_out5)
        dout5_out4 = torch.cat([dout5, out4], 1)
        dout4 = self.dlayer4(dout5_out4)
        dout4_out3 = torch.cat([dout4, out3], 1)
        dout3 = self.dlayer3(dout4_out3)
        dout3_out2 = torch.cat([dout3, out2], 1)
        dout2 = self.dlayer2(dout3_out2)
        dout2_out1 = torch.cat([dout2, out1], 1)
        dout1 = self.dlayer1(dout2_out1)
        return dout1

def weights_init(m):
    classname = m.__class__.__name__
    if classname.find('Conv') != -1:
        m.weight.data.normal_(0.0, 0.02)
    elif classname.find('BatchNorm') != -1:
        m.weight.data.normal_(1.0, 0.02)
        m.bias.data.fill_(0)

接下来，我们可以初始化一个 DfpNet实例。

下面，EXPO参数控制着我们的 Unet 特征图的指数：这直接决定了网络的大小（3 可以使网络拥有约 15 万个参数）。这对于一个有$3\times 128^2=\text{ca.}49k$ 输出的生成式 NN 而言，这个参数相对较小，但训练时间很快，而且可以防止过度拟合，因为我们在这里使用的数据集相对较小。因此，这是一个很好的起点。

# channel exponent to control network size
EXPO = 3

# setup network
net = DfpNet(channelExponent=EXPO)
#print(net) # to double check the details...

nn_parameters = filter(lambda p: p.requires_grad, net.parameters())
params = sum([np.prod(p.size()) for p in nn_parameters])

# crucial parameter to keep in view: how many parameters do we have?
print("Trainable params: {}   -> crucial! always keep in view... ".format(params)) 

net.apply(weights_init)

criterionL1 = nn.L1Loss()
optimizerG = optim.Adam(net.parameters(), lr=LR, betas=(0.5, 0.999), weight_decay=0.0)

targets = torch.autograd.Variable(torch.FloatTensor(BATCH_SIZE, 3, 128, 128))
inputs  = torch.autograd.Variable(torch.FloatTensor(BATCH_SIZE, 3, 128, 128))

Trainable params: 147363   -> crucial! always keep in view... 

指数为 3 时，该网络有 147555 个可训练参数。正如打印语句中的微妙提示所示，在训练 NN 时，这是一个需要时刻关注的关键数字。改变设置很容易得到一个拥有数百万个参数的网络，结果可能会出现各种收敛和过拟合问题。参数的数量肯定要与训练数据的数量相匹配，还应该与网络的深度成比例。不过，这三者之间的确切关系取决于问题的具体情况。

3.4.6 训练

最后，我们可以对 NN 进行训练。这一步可能要花费一些时间，因为训练会对所有 320 个样本运行 100 次，并不断评估验证样本，以跟踪 NN 的当前状态。

history_L1 = []
history_L1val = []

if os.path.isfile("network"):
  print("Found existing network, loading & skipping training")
  net.load_state_dict(torch.load("network")) # optionally, load existing network

else:
  print("Training from scratch")
  for epoch in range(EPOCHS):
      net.train()
      L1_accum = 0.0
      for i, traindata in enumerate(trainLoader, 0):
          inputs_curr, targets_curr = traindata
          inputs.data.copy_(inputs_curr.float())
          targets.data.copy_(targets_curr.float())

          net.zero_grad()
          gen_out = net(inputs)

          lossL1 = criterionL1(gen_out, targets)
          lossL1.backward()
          optimizerG.step()
          L1_accum += lossL1.item()

      # validation
      net.eval()
      L1val_accum = 0.0
      for i, validata in enumerate(valiLoader, 0):
          inputs_curr, targets_curr = validata
          inputs.data.copy_(inputs_curr.float())
          targets.data.copy_(targets_curr.float())

          outputs = net(inputs)
          outputs_curr = outputs.data.cpu().numpy()

          lossL1val = criterionL1(outputs, targets)
          L1val_accum += lossL1val.item()

      # data for graph plotting
      history_L1.append( L1_accum / len(trainLoader) )
      history_L1val.append( L1val_accum / len(valiLoader) )

      if epoch<3 or epoch%20==0:
          print( "Epoch: {}, L1 train: {:7.5f}, L1 vali: {:7.5f}".format(epoch, history_L1[-1], history_L1val[-1]) )

  torch.save(net.state_dict(), "network" )
  print("Training done, saved network")

程序输出结果为：

Training from scratch
Epoch: 0, L1 train: 0.29219, L1 vali: 0.23295
Epoch: 1, L1 train: 0.25406, L1 vali: 0.22507
Epoch: 2, L1 train: 0.22487, L1 vali: 0.21019
Epoch: 20, L1 train: 0.05228, L1 vali: 0.04134
Epoch: 40, L1 train: 0.03730, L1 vali: 0.03020
Epoch: 60, L1 train: 0.03236, L1 vali: 0.02523
Epoch: 80, L1 train: 0.03364, L1 vali: 0.02302
Training done, saved network

NN 终于训练完成！损失的绝对值应该已经很好地下降了：在标准设置下，验证损失的初始值约为 0.2，100 个epoch后约为 0.02。

让我们看看图表，直观地了解训练是如何随时间推移的。这对于识别训练的长期趋势非常重要。在实践中，很难发现命令行日志中 100 个左右的嘈杂数字的整体趋势是略微上升还是下降，而这在可视化中更容易发现。

import matplotlib.pyplot as plt

l1train = np.asarray(history_L1)
l1vali  = np.asarray(history_L1val)

plt.plot(np.arange(l1train.shape[0]),l1train,'b',label='Training loss')
plt.plot(np.arange(l1vali.shape[0] ),l1vali ,'g',label='Validation loss')
plt.legend()
plt.show()

程序输出如下：

您应该看到一条曲线，它在大约 40 个epoch中一直在下降，然后开始变得平缓。在最后一部分，它仍在缓慢下降，最重要的是，验证损失并没有增加。这肯定是过度拟合的迹象，我们应该避免。(试着人为减少训练数据量，就能有意识地造成过度拟合）。

3.4.7 训练进度和验证

如果仔细观察这张图，你会发现一些奇怪的现象：为什么验证损失低于训练损失？

当然，这些数据与训练数据类似，但从某种程度上说，它们略微 "艰难 "一些，因为网络在训练过程中肯定从未收到过任何验证样本。验证损失略微偏离训练损失是很自然的，但对于这些输入，L1 损失怎么会更低？

这是由于上述训练循环在 pytorch 中的实现方式造成的：虽然训练损失是在训练模式下通过 net.train() 进行评估的，但评估是在调用 net.eval() 之后进行的。这将关闭批量归一化，并禁用 dropout（如果激活）等功能。这会稍微改变评估结果。代码还会运行一个训练步骤，根据网络在一个epoch中的演化状态来测量图中每个点的损失。更新网络后，再运行验证样本。因此，所有验证样本所使用的状态都与历时的初始状态略有不同（希望更好一些）。由于这两个原因，验证损失可能会有偏差，在本例中通常会略低一些。

这里需要提醒的是：永远不要用训练数据来评估你的网络！因为过度拟合是一个非常常见的问题。至少要使用网络从未见过的数据，即验证数据，如果看起来不错，再尝试一些不同的（至少稍微超出分布范围的）输入，即_测试数据_。下一个单元格将在验证数据上运行训练有素的网络，并使用 showSbs 函数显示其中一个数据。

net.eval()
for i, validata in enumerate(valiLoader, 0):
    inputs_curr, targets_curr = validata
    inputs.data.copy_(inputs_curr.float())
    targets.data.copy_(targets_curr.float())
    
    outputs = net(inputs)
    outputs_curr = outputs.data.cpu().numpy()
    if i<1: showSbs(targets_curr[0] , outputs_curr[0], title="Validation sample %d"%(i*BATCH_SIZE))

输出结果为：

从外观上看，输入和网络输出之间至少应该有大致的相似之处。不过，我们还是留待测试数据时再进行更详细的评估吧。

3.4.8 测试评估

现在我们来看看实际的测试样本：在这种情况下，我们将使用新的机翼形状作为分布外（OOD）数据。这些形状是网络在任何训练样本中都从未见过的，因此它可以告诉我们网络对未见输入的泛化程度（验证数据不足以得出泛化结论）。

我们将使用与之前相同的可视化方式，正如伯努利方程所示，尤其是第一列中的_压力_对网络来说是一个具有挑战性的量。由于它与输入自由流速度和局部峰值呈立方比例关系，因此是网络最难推断的量。

下面的单元首先下载包含这些测试数据样本的较小档案，然后通过网络运行。评估循环还会计算累计 L1 误差，这样我们就可以量化网络在测试样本上的表现。

if not os.path.isfile('data-airfoils-test.npz'):
  import urllib.request
  url="https://physicsbaseddeeplearning.org/data/data_test.npz"
  print("Downloading test data, this should be fast...")
  urllib.request.urlretrieve(url, 'data-airfoils-test.npz')

nptfile=np.load('data-airfoils-test.npz')
print("Loaded {}/{} test samples\n".format(len(nptfile["test_inputs"]),len(nptfile["test_targets"])))

testdata = DfpDataset(nptfile["test_inputs"],nptfile["test_targets"])
testLoader  = torch.utils.data.DataLoader(testdata, batch_size=1, shuffle=False, drop_last=True) 

net.eval()
L1t_accum = 0.
for i, validata in enumerate(testLoader, 0):
    inputs_curr, targets_curr = validata
    inputs.data.copy_(inputs_curr.float())
    targets.data.copy_(targets_curr.float())

    outputs = net(inputs)
    outputs_curr = outputs.data.cpu().numpy()

    lossL1t = criterionL1(outputs, targets)
    L1t_accum += lossL1t.item()
    if i<3: showSbs(targets_curr[0] , outputs_curr[0],  title="Test sample %d"%(i))

print("\nAverage test error: {}".format( L1t_accum/len(testLoader) ))

Downloading test data, this should be fast...
Loaded 10/10 test samples
Average test error: 0.028802116494625808

默认设置下的平均测试误差应约为 0.03。由于输入已归一化，这意味着所有三个字段的平均误差为每个量的最大值的 3%。这对于新形状来说不算太差，但显然还有改进的余地。

通过观察可视化效果，您会发现输出结果中缺少特别高的压力峰值和较大的 yvelocities 小区。这主要是由于网络太小，没有足够的资源来重建细节。

尽管如此，我们还是成功地用一个非常小而快速的神经网络架构取代了相当复杂的 RANS 求解器。它 "开箱即用"（通过 pytorch），支持 GPU，可微分，引入的误差仅为百分之几。通过额外的改动和更多的数据，这个设置可以变得非常精确[CT21]。

3.4.9 下一步

有很多显而易见的方法可以尝试（见下面的建议），例如延长训练时间、扩大数据集、扩大网络规模等。

• 试验学习率、辍学率和网络规模，以减少测试集上的误差。在给定训练数据的情况下，你能把误差减到多小？
• 述设置使用的是归一化数据。相反，你可以恢复撤销归一化后的原始场 来检查网络在原始数量上的表现。
• 正如你所看到的，在这里你能从这个数据集得到的东西有点有限，你可以去这个项目的主 github 仓库下载更大的数据集，或者生成自己的数据。

3.5 监督方法讨论

前面的例子说明，我们可以很容易地使用监督训练来解决复杂的任务。主要的工作量是收集足够大的示例数据集。一旦有了这些数据集，我们就可以训练一个网络来近似这些解法采样的解流形，训练好的网络可以很快给出预测结果。在使用监督训练时，有几个要点需要注意。

3.5.1 一些需要注意的事项...

3.5.1.1 自然出发点

监督训练 是任何 DL项目的自然起点。在这里，我们的意思是总是，从使用尽可能少的数据进行完全监督测试开始是有意义的。这将是一个纯粹的过拟合测试，但如果你的网络不能快速收敛并在单个示例上提供非常好的性能，那么你的代码或数据就存在根本性的问题。因此，我们没有理由继续进行更复杂的设置，因为这会增加发现这些基本问题的难度。

总结前几节的零散评论，这里有一套建立 DL 项目的 "黄金法则"：

• 始终从 1 个样本的过拟合测试开始。- 检查网络有多少可训练参数。
• 缓慢增加训练数据量（以及可能的网络参数和深度）。
• 调整超参数（尤其是学习率）。
• 然后引入其他组件，如可变求解器或对抗训练。

3.5.1.2 稳定性

有监督训练的一个优点是非常稳定。当我们加入更复杂的物理模型或研究更复杂的 NN 架构时，情况也不会变得更好。

因此，在开始进行简单的过拟合测试时，请再次确保您能看到训练损失呈指数下降。这是一个很好的设置，可以找出作为最核心超参数的学习率的上限和合理范围。以后你可能还需要降低学习率，但至少应该能大致估算出 $\eta$ 的合适值。

3.5.1.3 了解你的数据

所有数据驱动型方法都遵循 "垃圾进垃圾出 "原则。因此，了解你所处理的数据非常重要。虽然没有放之四海而皆准的最佳方法，但我们强烈建议对数据集进行广泛的统计跟踪。一个好的起点是每个数量的平均值、标准偏差、最小值和最大值。如果其中某些值不正常，则表明数据集中存在不良样本。

这些值也可以很容易地通过直方图进行可视化，以追踪不需要的异常值。少量的异常值很容易使数据集出现偏差。

最后，检查不同数量之间的关系通常是个好主意，这样可以获得数据集所包含内容的一些直觉。下图给出了这一步骤的示例。

3.5.1.4 魔力在哪里?

在讨论 DL 方法时，尤其是在使用相对简单的训练方法时，你经常会听到这样的评论："这不就是对数据进行插值吗？"

嗯，是的！这正是NN应该做的。从某种程度上说，没有其他事情可做。这就是所有 DL 方法的目的。它们为我们提供训练时所看到的数据的平滑表示。即使我们以后在训练时使用了花哨的物理模型，NN 也只是调整它们的权重来代表它们接收到的信号，并将其再现出来。

由于炒作和大量成功案例，不熟悉 DL 的人往往会产生这样的印象：DL 就像人类的大脑一样工作，能够从数据集中提取基本和一般原理（"上帝的信息" 有人说过吗？）目前的技术水平并非如此。尽管如此，它仍是我们用来逼近复杂非线性函数的最强大工具。它是一个伟大的工具，但重要的是要记住，一旦我们正确设置了训练，我们从中得到的只是 NN 所训练的函数的近似值--不涉及任何魔法。

这意味着，你不应该指望网络能在它从未见过的数据上工作。从某种程度上说，神经网络之所以如此出色，正是因为它们能够准确地适应训练时接收到的信号，但与其他学习到的表征相比，神经网络实际上并不擅长外推。因此，我们不能指望神经网络能神奇地处理新的输入。相反，我们需要确保能够正确塑造输入空间，例如，通过归一化和关注不变量。

举个更具体的例子：如果您总是在$[0\dots 1]$的范围内训练您的网络输入，那么就不要指望它能在$[27\dots 39]$的输入下工作。在某些情况下，可以通过减去平均值来对输入和输出进行归一化，并通过标准偏差或合适的量化值进行归一化（确保这不会破坏数据中的重要相关性）。

经验法则：确保在实际训练 NN 时，输入尽可能与推理时要使用的输入相似。

这一点在接下来的章节中必须牢记：例如，如果我们希望 NN 与特定的模拟环境结合使用，那么在训练过程中就必须实际加入模拟器。否则，网络可能会专门处理预先计算的数据，而这些数据与将 NN 与求解器结合时产生的数据不同，即会出现_分布偏移_。

3.1.5.5 网格和栅格

前面的翼型示例使用了带有标准卷积的笛卡尔网格。就性能和稳定性而言，这些通常是最划算的。尽管如此，这里的整个讨论当然也适用于其他类型的卷积，例如，与图卷积结合的不太规则的网格，或具有连续卷积的基于粒子的数据（参见[非结构化网格和无网格方法](http://physicsbaseddeeplearning.org/others lagrangian.html)）。当切换到这些方法时，你通常会看到学习性能降低，以换取采样灵活性的提高。 最后，关于全连接层或一般的MLP：我们建议尽可能避免这些。对于任何结构化数据，如空间函数或一般的字段数据，卷积是优选的，并且不太可能过拟合。例如，你会注意到细胞神经网络通常不需要丢弃，因为它们通过构建很好地正则化了。对于MLP，您通常需要相当多的方法来避免过拟合。

3.5.2 监督训练简述

总之，监督训练具有以下特点。

✅ 优点：

• 训练速度非常快。
• 稳定简单。
• 是一个很好的起点。

❌ 缺点：需要大量数据。

• 性能、准确性和泛化能力不够优秀。
• 与外部“进程”（如嵌入到求解器中）的交互困难。

接下来的章节将解释如何缓解监督训练的这些缺点。首先，我们将通过软约束将模型方程引入图像中，然后我们将重新审视将数值模拟和学习方法结合起来所面临的挑战。

前几节中的有监督设置可以通过相当简单的训练过程快速获得近似解。然而，令人痛心的是，我们只是将物理模型和数值方法作为一种 "外部 "工具，来产生一大堆数据。我们人类拥有大量关于如何用数学方法描述物理过程的知识。正如以下几章所展示的，我们可以通过人类的物理知识来指导训练过程，从而改进训练过程。

4.1 使用物理模型

给定一个具有时间演化的 $\mathbf{u}(\mathbf{x},t)$ 的 PDE，我们通常可以用 $\mathcal{F}$ 的导数函数来表示，方法是 $${\mathbf{u}}_t=\mathcal{F}({\mathbf{u}}_x,{\mathbf{u}}_{xx},...{\mathbf{u}}_{xx...x}),$$

其中 $\mathbf{x}$ 的下标表示相对于一个空间维度的高阶空间导数（当然也可以包括相对于不同坐标轴的混合导数）。${\mathbf{u}}_t$ 表示随时间的变化。

在这种情况下，我们可以用神经网络来近似未知的 $\mathbf{u}$ 本身。如果近似值（我们称之为 $\tilde{\mathbf{u}}$ )是准确的，那么 PDE 应该自然满足。换句话说，残差 R 应该等于零： $$R={\mathbf{u}}_t-\mathcal{F}({\mathbf{u}}_x,{\mathbf{u}}_{xx},...{\mathbf{u}}_{xx...x})=0.$$

这与训练神经网络的目标很好地结合在一起：我们可以结合直接损失项来训练最小化残差。与之前类似，除了残差项之外，我们还可以使用预先计算出的解 $[x_0,y_0],...[x_n,y_n]$，将 $\mathbf{u}$ 与 $\mathbf{u}(\mathbf{x})=y$ 作为约束条件。这一点通常很重要，因为除非指定初始条件和边界条件，否则大多数实际的 PDE 都没有唯一解。因此，如果我们只考虑 $R$，可能会得到带有随机偏移或其他不理想成分的解。因此，监督采样点有助于在某些地方固定解。

现在我们的训练目标变成了：

$${\text{arg min}}_{\theta }\sum _i{\alpha }_0(f(x_i;\theta )-y_i^{\ast }{)}^2+{\alpha }_{1}R(x_i)，$$ 其中，${\alpha }_{0,1}$ 表示超参数，分别缩放监督项和残差项的贡献。当然，我们还可以在此添加适当比例系数的附加残差项。

注意方程(1)中两个不同项的含义很有启发：第一项是传统的、有监督的 L2 损失。如果我们只对这一损失进行优化，我们的网络就能很好地学习近似训练样本，但可能会在解中平均出多种模式，在样本点之间的区域表现不佳。相反，如果我们只优化第二项（物理残差），我们的神经网络可能会在局部满足 PDE，但产生的解仍可能与训练数据相去甚远。出现这种情况的原因可能是解中存在 "空位"，即不同的解均满足残差。因此，我们要同时优化这两个目标，以便在最佳情况下，网络在学习近似训练数据的特定解的同时，还能捕捉到底层 PDE 的知识。

需要注意的是，与用于监督训练的数据样本类似，我们无法保证残差项 $R$ 在训练过程中实际为零。训练过程的非线性优化将尽可能减少监督项和残差项，但这并不能保证。仍有可能存在大量非零残差项。我们将在接下来的代码示例中更详细地讨论这个问题，现在重要的是要记住，这种方式的物理约束只代表_软约束_，并不能保证最大限度地减少这些约束。

前面的概述并没有明确说明神经网络是如何产生 $\mathbf{u}$ 的。在这里，我们可以区分两种不同的方法：通过选择目标函数的显式表示（下文中的 v1），或通过使用全连接神经网络来表示解（v2）。例如，对于 v1，我们可以设置一个空间网格（或图形，或一组样本点），而在第二种情况下，不存在显式表示，而是由神经网络接收空间坐标，在查询位置生成解。下面我们将详细介绍这两种变体。

4.2 变体 1: 用于显式表示的残差导数

对于变式 1，我们选择离散化，并建立一个覆盖目标域的计算网格。在不失一般性的前提下，我们假设这是一个笛卡尔网格，以 $\mathbf{p}$ 的位置对空间进行采样。现在，训练一个 NN 来生成网格上的解：$\mathbf{u}(\mathbf{p})=f(\mathbf{x};\theta )$。对于规则网格，CNN 将是 $f$ 的最佳选择，而对于三角形网格，我们可以使用图网络，或使用点演算的粒子网络。

现在，我们可以在计算网格上将 $R$ 的方程离散化，并用我们选择的方法计算导数。唯一需要注意的是：为了将残差纳入训练，我们必须制定评估方法，使深度学习框架可以通过计算进行反向传播。由于我们的网络 $f()$ 产生了解 $\mathbf{u}$，而残差取决于它（$R(\mathbf{u})$），因此我们至少需要 $\partial R/\partial \mathbf{u}$，这样才能对权重 $\theta$ 进行梯度反向传播。幸运的是，如果我们用 DL 框架的操作来表示 $R$，那么该框架的反向传播功能就能解决这个问题。

这种变体在DL中有着相当悠久的 "传统"，例如，Tompson et al. [TSSP17]很早就提出了这种变体来学习无散度运动。举个具体例子：如果我们的目标是学习无散度的速度 $\mathbf{u}(t)$，我们就可以使用这种训练方法来训练 NN，而无需预先计算无发散速度场作为训练数据。为了简洁起见，我们在这里将去掉空间指数（$\mathbf{p}$），并将重点放在$t$上，我们同样可以将其简化：无散度必须在任何时候都保持不变，因此我们可以考虑从$t=0$开始的单一步骤，即从有散度的$\mathbf{u}(0)$到无散度的$\mathbf{u}(1)$的归一化步骤。对于正常求解器，我们必须计算一个压力 $p={\nabla }^{-2}\mathbf{u}(0)$，使得 $\mathbf{u}(1)=\mathbf{u}(0)-\nabla p$。这就是著名的矢量微积分基本定理，或Helmholtz decomposition，把矢量场分成_solenoidal_（无散度）和非旋转部分（压力梯度）。

为了学习这种分解，我们可以在计算网格上用 CNN 近似计算 $p$：$p=f(\mathbf{u}(0);\theta )$。学习目标变成最小化 $\mathbf{u}(0)$ 的发散，也就是最小化 $\nabla \cdot (\mathbf{u}(0)-\nabla f(\mathbf{u}(0);\theta ))$ 。为了实现这个残差，我们只需在计算网格上提供 $\mathbf{u}$ 的发散算子 $(\nabla \cdot )$ 即可。这通常可以通过 DL 框架中的卷积层轻松实现，该卷积层包含发散的有限差分权重。非常好的是，在这种情况下，我们甚至不需要额外的监督样本，通常只需使用这种残差公式进行训练即可。另外，与下面的变式 2 不同，我们可以直接处理相当大的解空间（我们并不局限于学习单个解）。我们可以在 代码库 中找到一个实现示例。

总的来说，变式 1 与可微分物理训练有很多共同之处（它基本上是一个子集）。由于我们将在 Introduction to Differentiable Physics和之后的文章中更详细地讨论可微分物理，从现在起，我们将专注于直接 NN 表示（变体 2）。

4.3 变体 2: 来自神经网络表示的导数

采用物理残差作为软约束的第二种变体是使用全连接的 NN 来表示 $\mathbf{u}$。这种_physics-infformed_方法由 Raissi 等人推广[RPK19]，它有一些有趣的优点和缺点，我们将在下文中概述。在下面的代码示例和讨论中，我们将以物理信息版本（变体 2）为目标。

这里的中心思想是，我们在学习问题中所追求的上述通用函数 $f$ 也可以用来获得物理场的表示，例如，满足 $R=0$ 的场 $\mathbf{u}$。这意味着 $\mathbf{u}(\mathbf{x})$ 将转化为 $\mathbf{u}(\mathbf{x},\theta )$，其中我们选择的 NN 参数 $\theta$ 可以尽可能精确地表示出所需的 $\mathbf{u}$。

这种观点的一个很好的副作用是，NN 表征本质上支持导数的计算。导数 $\partial f/\partial \theta$ 是通过梯度下降进行学习的关键基石，这在 Overview 中已有解释。现在，我们可以用同样的工具来计算空间导数，比如 $\partial \mathbf{u}/\partial x$，注意，上面对于 $R$，我们把这个导数写成了简写符号 ${\mathbf{u}}_x$。对于随时间变化的函数，这当然也适用于 $\partial \mathbf{u}/\partial t$，即上面符号中的 ${\mathbf{u}}_t$ 。

因此，对于由${\mathbf{u}}_t$和${\mathbf{u}}_x$项组成的通用$R$，一旦我们有了代表$\mathbf{u}$的NN，我们就可以依靠DL框架的反向传播算法来计算这些导数。从本质上讲，这给我们提供了一个函数（NN），该函数接收空间和时间坐标来生成 $\mathbf{u}$ 的解。因此，输入通常是相当低维的，例如三维情况下随时间变化的 3+1 个值，通常会产生一个标量值或空间向量。由于缺乏明确的空间采样点，这里选择了 MLP（即全连接 NN）架构。

举个简单的例子，一维Burgers方程， $\frac{\partial u}{\partial t}+u\nabla u=\nu \nabla \cdot \nabla u$ ，我们可以直接提出一个损失项 $R=\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}-\nu \frac{{\partial }^{2}u}{\partial x^2}$，在训练时应尽可能将其最小化。对于每个项，例如 $\frac{\partial u}{\partial x}$，我们可以简单地查询实现 $u$ 的 DL 框架，以获得相应的导数。对于高阶导数，例如 $\frac{{\partial }^{2}u}{\partial x^2}$，我们可以简单地多次查询框架的导函数。在下一节中，我们将举例说明如何在 tensorflow 中实现这一功能。

4.4 到目前为止的总结

上述方法为我们提供了一种将物理方程作为软约束纳入 DL 学习的方法：残差损失。通常情况下，这种设置适用于逆问题，即我们想要找到一个 PDE 解的某些测量或观测结果。由于重构的成本很高（将在下文中演示），解流形不应该过于复杂。例如，仅使用物理残余损失通常无法捕捉到广泛的解，如之前的监督机翼示例。

为了说明物理信息损失是如何在变体 2 中发挥作用的，让我们以重建任务作为反问题的例子。我们将使用Burgers方程 $\frac{\partial u}{\partial t}+u\nabla u=\nu \nabla \cdot \nabla u$ 作为一个简单但非线性的一维方程，我们在时间 $t=0.5$ 时有一系列_观测值_。解应该符合Burgers方程的残差公式，并与观测结果相匹配。此外，让我们在计算域两侧施加狄利克雷边界条件 $u=0$，并将解定义在时间区间 $t\in [0,1]$ 上。

请注意，与之前的正向模拟示例类似，我们仍将使用 128 个点（$n=128$）对解进行采样，但现在我们通过 NN 进行了离散化。因此，我们也可以对中间的点进行采样，而无需明确选择用于插值的基函数。现在通过 NN 进行离散化，可以在内部决定如何使用其自由度来安排激活函数作为基函数。因此，我们无法直接控制重构。

运行于 colab

5.1 原理

根据 2.5 模型与方程和前一节中的 $x,y^{\ast }$ 符号，这个重构问题意味着我们要求解的是： $${\text{arg min}}_{\theta }\sum _i(f(x_i;\theta )-y_i^{\ast }{)}^2+R(x_i),$$

其中 $x_i$ 表示时空点 $x_i=[p_i,t_i]$，参考解为 $y_i^{\ast }=y^{\ast }(x_i)$，索引 $i$ 表示数据集的不同采样点。$f$ 和 $y^{\ast }$ 都表示 $u$ 在空间和时间的不同位置的解，由于我们处理的是一维速度，所以 $f,y^{\ast }：{\mathbb{R}}^2\to \mathbb{R}$。在这个例子中，$y^{\ast }$ 表示布尔格斯方程 $\mathcal{P}$ 的参考值 $u$，在所有选定的时空点 $x_i=[p_i,t_i]$ 上，$f$ 应尽可能接近这个参考值。

上述第一项是 "监督 "数据项，第二项表示残差函数 $R$。它收集了$f(;\theta )$及其导数的额外评估，以计算$\mathcal{P}$的残差。这种方法--使用神经网络的导数来计算PDE残差--通常被称为_physics-informed_方法，产生一个physics-informed神经网络（PINN）[RPK19] 来表示逆重建问题的解。

上述残差函数$\mathcal{R}$收集$f(;\theta )$及其导数的其他评估值以构建$\mathcal{P}$的残差。这种方法——使用神经网络的导数来计算PDE残差——通常称为物理信息(physicals-informed)方法，产生一个物理信息神经网络(PINN)来表示反演重构问题的解。

因此，在上述公式中，$R$ 应简单地收敛为零。为简单起见，我们省略了目标函数中的缩放因子。请注意，实际上我们在这里处理的只是 $\mathcal{P}$ 的单个解 $u$ 的单点样本。

5.2 准备工作

现在我们先用 tensorflow 后端加载 phiflow，然后初始化随机采样。(注：本例使用的是旧版本 1.5.1 的 phiflow）。

!pip install --quiet phiflow==1.5.1

from phi.tf.flow import *
import numpy as np

#rnd = TF_BACKEND  # for phiflow: sample different points in the domain each iteration
rnd = math.choose_backend(1)  # use same random points for all iterations

注：在python3.9之后，优于collections.Iterable被废弃，第7行语句会报错。这里需要打开文件C:\ProgramData\anaconda3\Lib\site-packages\phi\backend\scipy_backend.py，修改第37行为 if isinstance(values, collections.abc.Iterable):，注意重启内核。

我们在这里导入了 phiflow，但不会像Simple Forward Simulation of Burgers Equation with phiflow中那样用它来计算 PDE 的解。相反，我们将使用 NN 的导数（如上一节所述）来为训练设置损失公式。

接下来，我们建立了一个简单的NN，其中包含8个完全连接的层和每个20个单元的tanh激活。

我们还将定义 boundary_tx 函数和 open_boundary 函数，前者给出解的约束条件数组（本例中所有约束条件均为 $=0.5$），后者存储 $x=\pm 1$ 为 0 的约束条件。

def network(x, t):
    """ Dense neural network with 8 hidden layers and 3021 parameters in total.
        Parameters will only be allocated once (auto reuse).
    """
    y = math.stack([x, t], axis=-1)
    for i in range(8):
        y = tf.layers.dense(y, 20, activation=tf.math.tanh, name='layer%d' % i, reuse=tf.AUTO_REUSE)
    return tf.layers.dense(y, 1, activation=None, name='layer_out', reuse=tf.AUTO_REUSE)

def boundary_tx(N):
    x = np.linspace(-1,1,128)
    # precomputed solution from forward simulation:
    u = np.asarray( [0.008612174447657694, 0.02584669669548606, 0.043136357266407785, 0.060491074685516746, 0.07793926183951633, 0.0954779141740818, 0.11311894389663882, 0.1308497114054023, 0.14867023658641343, 0.1665634396808965, 0.18452263429574314, 0.20253084411376132, 0.22057828799835133, 0.23865132431365316, 0.25673879161339097, 0.27483167307082423, 0.2929182325574904, 0.3109944766354339, 0.3290477753208284, 0.34707880794585116, 0.36507311960102307, 0.38303584302507954, 0.40094962955534186, 0.4188235294008765, 0.4366357052408043, 0.45439856841363885, 0.4720845505219581, 0.4897081943759776, 0.5072391070000235, 0.5247011051514834, 0.542067187709797, 0.5593576751669057, 0.5765465453632126, 0.5936507311857876, 0.6106452944663003, 0.6275435911624945, 0.6443221318186165, 0.6609900633731869, 0.67752574922899, 0.6939334022562877, 0.7101938106059631, 0.7263049537163667, 0.7422506131457406, 0.7580207366534812, 0.7736033721649875, 0.7889776974379873, 0.8041371279965555, 0.8190465276590387, 0.8337064887158392, 0.8480617965162781, 0.8621229412131242, 0.8758057344502199, 0.8891341984763013, 0.9019806505391214, 0.9143881632159129, 0.9261597966464793, 0.9373647624856912, 0.9476871303793314, 0.9572273019669029, 0.9654367940878237, 0.9724097482283165, 0.9767381835635638, 0.9669484658390122, 0.659083299684951, -0.659083180712816, -0.9669485121167052, -0.9767382069792288, -0.9724097635533602, -0.9654367970450167, -0.9572273263645859, -0.9476871280825523, -0.9373647681120841, -0.9261598056102645, -0.9143881718456056, -0.9019807055316369, -0.8891341634240081, -0.8758057205293912, -0.8621229450911845, -0.8480618138204272, -0.833706571569058, -0.8190466131476127, -0.8041372124868691, -0.7889777195422356, -0.7736033858767385, -0.758020740007683, -0.7422507481169578, -0.7263049162371344, -0.7101938950789042, -0.6939334061553678, -0.677525822052029, -0.6609901538934517, -0.6443222327338847, -0.6275436932970322, -0.6106454472814152, -0.5936507836778451, -0.5765466491708988, -0.5593578078967361, -0.5420672759411125, -0.5247011730988912, -0.5072391580614087, -0.4897082914472909, -0.47208460952428394, -0.4543985995006753, -0.4366355580500639, -0.41882350871539187, -0.40094955631843376, -0.38303594105786365, -0.36507302109186685, -0.3470786936847069, -0.3290476440540586, -0.31099441589505206, -0.2929180880304103, -0.27483158663081614, -0.2567388003912687, -0.2386513127155433, -0.22057831776499126, -0.20253089403524566, -0.18452269630486776, -0.1665634500729787, -0.14867027528284874, -0.13084990929476334, -0.1131191325854089, -0.09547794429803691, -0.07793928430794522, -0.06049114408297565, -0.0431364527809777, -0.025846763281087953, -0.00861212501518312] );
    t = np.asarray(rnd.ones_like(x)) * 0.5
    perm = np.random.permutation(128) 
    return (x[perm])[0:N], (t[perm])[0:N], (u[perm])[0:N]

def _ALT_t0(N): # alternative, impose original initial state at t=0
    x = rnd.random_uniform([N], -1, 1)
    t = rnd.zeros_like(x)
    u = - math.sin(np.pi * x)
    return x, t, u

def open_boundary(N):
    t = rnd.random_uniform([N], 0, 1)
    x = math.concat([math.zeros([N//2]) + 1, math.zeros([N//2]) - 1], axis=0)
    u = math.zeros([N])
    return x, t, u

最重要的是，我们现在还可以构建我们希望最小化的残差损失函数f，以便引导 NN 为我们的模型方程获取一个解。从最上面的方程中可以看出，我们需要与 $t$、$x$ 相关的导数以及 $x$ 的二次导数。下面 f 的前三行就是这样做的。

之后，我们只需将导数结合起来，就能形成Burgers方程。这里我们使用 phiflow 的 gradient函数：

def f(u, x, t):
    """ Physics-based loss function with Burgers equation """
    u_t = gradients(u, t)
    u_x = gradients(u, x)
    u_xx = gradients(u_x, x)
    return u_t + u*u_x - (0.01 / np.pi) * u_xx

接下来，让我们在内域中设置采样点，以便将解与之前在 phiflow 中进行的正向模拟进行比较。

下面的单元格分配两个张量:grid_x将使用128个单元覆盖我们的域，即-1到1的范围，而grid_t将使用33个时间点对时间区间$[0，1]$进行采样。

最后的math.expand_dims()调用只是添加了一个batch维度，这样得到的张量与下面的示例兼容。

N=128
grids_xt = np.meshgrid(np.linspace(-1, 1, N), np.linspace(0, 1, 33), indexing='ij')
grid_x, grid_t = [tf.convert_to_tensor(t, tf.float32) for t in grids_xt]

# create 4D tensor with batch and channel dimensions in addition to space and time
# in this case gives shape=(1, N, 33, 1)
grid_u = math.expand_dims(network(grid_x, grid_t))

现在，grid_u包含一个完整的图，用于在$128\times 33$个位置上评估我们的NN，一旦我们通过session.run运行它，就会以$[1，128，33，1]$数组的形式返回结果。让我们试一试：我们可以初始化一个 TF 会话，评估 grid_u 并将其显示在图像中，就像我们之前计算的 phiflow 解一样。

(注意，我们将使用Simple Forward Simulation of Burgers Equation with phiflow中的 show_state。因此，X 轴并不显示实际的模拟时间，而是显示被放大16倍的32个时间步，以使图像中随时间变化更容易看到。)

import pylab as plt
print("Size of grid_u: "+format(grid_u.shape))

session = Session(None)
session.initialize_variables()

def show_state(a, title):
    for i in range(4): a = np.concatenate( [a,a] , axis=3)
    a = np.reshape( a, [a.shape[1],a.shape[2]*a.shape[3]] )
    fig, axes = plt.subplots(1, 1, figsize=(16, 5))
    im = axes.imshow(a, origin='upper', cmap='inferno')
    plt.colorbar(im) ; plt.xlabel('time'); plt.ylabel('x'); plt.title(title)
    
print("Randomly initialized network state:")
show_state(session.run(grid_u),"Uninitialized NN")

运行输出为：

Size of grid_u: (1, 128, 33, 1)
Randomly initialized network state:

这种可视化效果已经显示出空间和时间上的平滑过渡。到目前为止，这纯粹是我们正在采样的 NN 的随机初始化。因此，到目前为止，它与我们基于 PDE 模型的解决方案毫无关系。接下来的步骤将根据数据（来自boundary函数）和来自f的模型约束条件对约束条件进行实际评估，以获取 PDE 的实际解。

5.3 损失函数和训练

作为学习过程的目标，我们现在可以将direct 约束条件，即 $t=0.5$ 时的解和狄利克雷 $u=0$ 边界条件与来自 PDE 残差的损失结合起来。对于这两个边界约束条件，我们将在下方使用 100 个点，然后在内部区域使用额外的 1000 个点对解进行采样。

直接约束通过 network(x, t)[:, 0] - u 进行评估，其中 x 和 t 是我们希望对解进行采样的时空位置，而 u 则提供相应的真值。

对于物理损失点，我们没有真值解，但我们只需通过 NN 导数评估 PDE 残差，看看解是否满足 PDE 模型。如果不满足，就会直接产生误差，需要在优化过程中通过更新步骤来减少误差。相应的表达式如下f(network(x, t)[:, 0], x, t)。请注意，对于数据和物理项，network()[:, 0]表达式不会从 $L^2$ 评估中删除任何数据，它们只是丢弃了网络返回的 $(n,1)$ 张量的最后一个 size-1 维。

# Boundary loss
N_SAMPLE_POINTS_BND = 100
x_bc, t_bc, u_bc = [math.concat([v_t0, v_x], axis=0) for v_t0, v_x in zip(boundary_tx(N_SAMPLE_POINTS_BND), open_boundary(N_SAMPLE_POINTS_BND))]
x_bc, t_bc, u_bc = np.asarray(x_bc,dtype=np.float32), np.asarray(t_bc,dtype=np.float32) ,np.asarray(u_bc,dtype=np.float32)
#with app.model_scope():
loss_u = math.l2_loss(network(x_bc, t_bc)[:, 0] - u_bc)  # normalizes by first dimension, N_bc

# Physics loss inside of domain
N_SAMPLE_POINTS_INNER = 1000
x_ph, t_ph = tf.convert_to_tensor(rnd.random_uniform([N_SAMPLE_POINTS_INNER], -1, 1)), tf.convert_to_tensor(rnd.random_uniform([N_SAMPLE_POINTS_INNER], 0, 1))
loss_ph = math.l2_loss(f(network(x_ph, t_ph)[:, 0], x_ph, t_ph))  # normalizes by first dimension, N_ph

# Combine
ph_factor = 1.
loss = loss_u + ph_factor * loss_ph # allows us to control the relative influence of loss_ph 

optim = tf.train.GradientDescentOptimizer(learning_rate=0.02).minimize(loss)
#optim = tf.train.AdamOptimizer(learning_rate=0.001).minimize(loss) # alternative, but not much benefit here

上面的代码只是初始化了损失的评估，我们仍然没有进行任何优化步骤，但我们终于可以开始着手处理这个问题了。

尽管方程很简单，但收敛速度通常很慢。迭代本身的计算速度很快，但这种设置需要大量的迭代。为了将运行时间控制在合理范围内，我们默认只进行 10k 次迭代 (ITERS)。你可以增加这个值，以获得更好的结果。

session.initialize_variables()

import time
start = time.time()

ITERS = 10000
for optim_step in range(ITERS+1):
  _, loss_value = session.run([optim, loss])
  if optim_step<3 or optim_step%1000==0: 
        print('Step %d, loss: %f' % (optim_step,loss_value))
        #show_state(grid_u)
        
end = time.time()
print("Runtime {:.2f}s".format(end-start))

Step 0, loss: 0.276599
Step 1, loss: 0.155847
Step 2, loss: 0.125085
Step 1000, loss: 0.053451
Step 2000, loss: 0.050352
Step 3000, loss: 0.047439
Step 4000, loss: 0.045276
Step 5000, loss: 0.042576
Step 6000, loss: 0.040332
Step 7000, loss: 0.036957
Step 8000, loss: 0.031184
Step 9000, loss: 0.028826
Step 10000, loss: 0.027606
Runtime 33.19

训练可能需要相当长的时间，在典型的notebook上大约需要2分钟，但至少误差显著下降(从约0.2降低到约0.03)，网络似乎成功地收敛到一个解。

让我们通过在规则网格的中心评估网络来显示网络的重构，这样我们就可以将解显示为图像。请注意，这实际上相当昂贵，我们必须对所有的$128\times 32$个采样点运行整个网络及其几千个权重。

不过，乍一看它看起来相当不错。与上面显示的随机初始化相比，发生了非常明显的变化:

5.4 评估

让我们更详细地比较一下解法。下面是用于约束解法的实际样本点（时间步长为 16，$t=1/2$）（灰色）和重构解法（蓝色）：

u = session.run(grid_u)

# solution is imposed at t=1/2 , which is 16 in the array
BC_TX = 16 
uT = u[0,:,BC_TX,0]

fig = plt.figure().gca()
fig.plot(np.linspace(-1,1,len(uT)), uT, lw=2, color='blue', label="NN")
fig.scatter(x_bc[0:100], u_bc[0:100], color='gray', label="Reference")
plt.title("Comparison at t=1/2")
plt.xlabel('x'); plt.ylabel('u'); plt.legend()

反馈结果为：

<matplotlib.legend.Legend at 0x7f8eca2a7a00>

域两侧的情况还不错（满足了 Dirichlet 边界条件 $u=0$），但中心（位于 $x=0$）的激波没有得到很好的表现。

让我们来看看 $t=0$ 时的初始状态重建得如何。这是最有趣、也是最棘手的部分（其余部分基本上都是根据模型方程和边界条件，给出第一种状态）。

事实证明，初始状态的精确度其实并不高：PINN 的蓝色曲线与参考数据（灰色显示）的约束条件相去甚远...随着迭代次数的增加，求解结果会越来越好，但对于这种相当简单的情况来说，迭代次数之多令人吃惊。

# ground truth solution at t0
t0gt = np.asarray( [ [-math.sin(np.pi * x) * 1.] for x in np.linspace(-1,1,N)] )
velP0 = u[0,:,0,0]

fig = plt.figure().gca()
fig.plot(np.linspace(-1,1,len(velP0)), velP0, lw=2, color='blue', label="NN")
fig.plot(np.linspace(-1,1,len(t0gt)), t0gt, lw=2, color='gray', label="Reference") 
plt.title("Comparison at t=0")
plt.xlabel('x'); plt.ylabel('u'); plt.legend()

运行结果为：

尤其是$x=\pm 1/2$处的最大值/最小值相差甚远，而$x=\pm 1$处的边界不符合要求：解不为零。

我们有用于该模拟的前向模拟器，因此我们可以使用网络的 $t=0$ 解来评估时间评估的重建效果。这可以衡量通过 PINN 损失的软约束捕捉模型方程时间演化的程度。

下图显示了蓝色的初始状态，以及在 $t=8/32$ 和 $t=15/32$ 时的两个演化状态。请注意，这些都是模拟版本，接下来我们将展示学习版本。

(注：下面的代码段还有一些可选代码，用于显示[STEPS//4]处的状态。默认情况下这些代码都已注释，如果你愿意，可以取消注释或添加额外的代码，以显示更多的时间演化过程）。

# re-simulate with phiflow from solution at t=0
DT = 1./32.
STEPS = 32-BC_TX # depends on where BCs were imposed
INITIAL = u[...,BC_TX:(BC_TX+1),0] # np.reshape(u0, [1,len(u0),1]) 
print(INITIAL.shape)

DOMAIN = Domain([N], boundaries=PERIODIC, box=box[-1:1])
state = [BurgersVelocity(DOMAIN, velocity=INITIAL, viscosity=0.01/np.pi)]
physics = Burgers()

for i in range(STEPS):
    state.append( physics.step(state[-1],dt=DT) )

# we only need "velocity.data" from each phiflow state
vel_resim = [x.velocity.data for x in state]

fig = plt.figure().gca()
pltx = np.linspace(-1,1,len(vel_resim[0].flatten()))
fig.plot(pltx, vel_resim[ 0].flatten(),       lw=2, color='blue',  label="t=0")
#fig.plot(pltx, vel_resim[STEPS//4].flatten(), lw=2, color='green', label="t=0.125")
fig.plot(pltx, vel_resim[STEPS//2].flatten(), lw=2, color='cyan',  label="t=0.25")
fig.plot(pltx, vel_resim[STEPS-1].flatten(),  lw=2, color='purple',label="t=0.5")
#fig.plot(pltx, t0gt, lw=2, color='gray', label="t=0 Reference") # optionally show GT, compare to blue
plt.title("Resimulated u from solution at t=0")
plt.xlabel('x'); plt.ylabel('u'); plt.legend()

注：如果报错np.object was a deprecated...则需要修改：

1. 文件C:\ProgramData\anaconda3\Lib\site-packages\phi\backend\scipy_backend.py第64行
2. 文件C:\ProgramData\anaconda3\Lib\site-packages\phi\struct\struct.py第306行
3. 文件C:\ProgramData\anaconda3\Lib\site-packages\phi\struct\struct.py第315行

修改np.object为np.dtype('O')

下面是 u 在相同时间步长下的 PINN 输出：

velP = [u[0,:,x,0] for x in range(33)]
print(velP[0].shape)

fig = plt.figure().gca()
fig.plot(pltx, velP[BC_TX+ 0].flatten(),       lw=2, color='blue',  label="t=0")
#fig.plot(pltx, velP[BC_TX+STEPS//4].flatten(), lw=2, color='green', label="t=0.125")
fig.plot(pltx, velP[BC_TX+STEPS//2].flatten(), lw=2, color='cyan',  label="t=0.25")
fig.plot(pltx, velP[BC_TX+STEPS-1].flatten(),  lw=2, color='purple',label="t=0.5")
plt.title("NN Output")
plt.xlabel('x'); plt.ylabel('u'); plt.legend()

根据目测标准判断，这两个版本的 $u$ 看起来非常相似，但随着时间的推移，误差会逐渐增大，存在显著差异，这并不奇怪。尤其是在$x=0$处激波附近解的陡峭化没有被很好地 "捕捉 "到。不过从这两幅图中很难看出来，让我们量化误差，显示实际差异：

error = np.sum( np.abs( np.asarray(vel_resim[0:16]).flatten() - np.asarray(velP[BC_TX:BC_TX+STEPS]).flatten() )) / (STEPS*N)
print("Mean absolute error for re-simulation across {} steps: {:7.5f}".format(STEPS,error))

fig = plt.figure().gca()
fig.plot(pltx, (vel_resim[0       ].flatten()-velP[BC_TX         ].flatten()), lw=2, color='blue',  label="t=5")
fig.plot(pltx, (vel_resim[STEPS//4].flatten()-velP[BC_TX+STEPS//4].flatten()), lw=2, color='green', label="t=0.625")
fig.plot(pltx, (vel_resim[STEPS//2].flatten()-velP[BC_TX+STEPS//2].flatten()), lw=2, color='cyan',  label="t=0.75")
fig.plot(pltx, (vel_resim[STEPS-1 ].flatten()-velP[BC_TX+STEPS-1 ].flatten()), lw=2, color='purple',label="t=1")
plt.title("u Error")
plt.xlabel('x'); plt.ylabel('MAE'); plt.legend()

上述代码将计算出地面实况再模拟与 PINN 演化之间的平均绝对误差约为 $1.5\cdot 10^{-2}$，这对于模拟的数值范围而言意义重大。为了与正向模拟和后续案例进行比较，这里还提供了所有随时间变化的步骤，并配有彩色地图。

# show re-simulated solution again as full image over time
sn = np.concatenate(vel_resim, axis=-1)
sn = np.reshape(sn, list(sn.shape)+[1] ) # print(sn.shape)
show_state(sn,"Re-simulated u")

接下来，我们将存储 $t=0\dots 1$ 时间间隔内的完整解，以便稍后将其与常规正向求解的完整解进行比较，并与微分物理解进行比较。因此，请继续关注完整的评估和比较。这将在Burgers Optimization with a Differentiable Physics Gradient中，在我们讨论完如何运行微分物理优化的细节之后进行。

vels = session.run(grid_u) # special for showing NN results, run through TF 
vels = np.reshape( vels, [vels.shape[1],vels.shape[2]] )

# save for comparison with other methods
import os; os.makedirs("./temp",exist_ok=True)
np.savez_compressed("./temp/burgers-pinn-solution.npz",vels) ; print("Vels array shape: "+format(vels.shape))

结果为：

Vels array shape: (128, 33)

5.5 后续步骤

当然，这种设置只是 PINN 和物理软约束的一个起点。设置参数的选择是为了相对快速地运行。正如我们在接下来的章节中所展示的，通过将求解器和学习更紧密地结合在一起，这种逆求解的性能可以得到大幅提升。

不过，上述 PINN 设置的解也可以直接加以改进。例如，可以尝试：

• 调整训练参数，以进一步减少误差，同时不使解发散。
• 调整 NN 架构以进一步改进（但要跟踪权重计数）。
• 激活不同的优化器，并观察其行为变化（这通常需要调整学习率）。请注意，在这个相对简单的例子中，更复杂的优化器并不一定做得更好。
• 或者修改设置，使测试案例更有趣：例如，将边界条件移到模拟时间的较后一点，使重建的时间间隔更大。

目前的好消息是我们有了一种 DL 方法，可以通过最小化残差，以软约束的形式包含物理定律。不过，正如前面那个非常简单的例子所示，这只是一个概念上的起点。

积极的一面是，我们可以利用反向传播的 DL 框架来计算模型的导数。与此同时，这也使我们在这些导数的可靠性方面受到所学表示法的支配。而且，每个导数都需要通过整个网络进行反向传播。这可能会非常昂贵，尤其是对于高阶导数。

虽然设置相对简单，但通常很难控制。NN 可以灵活地自行完善求解，但同时，当它不能专注于求解的正确区域时，就需要一些技巧。

6.1 这算是“机器学习”吗?

说到这里，大家可能还会想到一个问题：我们真的可以称它为机器学习吗？当然，这样的命名问题是肤浅的--如果一种算法是有用的，那么它叫什么名字并不重要。不过，在这里，这个问题有助于强调机器学习或优化等领域的算法通常具有的一些重要特性。

不把上一个笔记本的优化称为机器学习（ML）的一个主要原因是，我们测试和限制解的位置就是我们感兴趣的最终位置。因此，训练集、验证集和测试集之间没有真正的区别。为已知和给定的样本集计算解更类似于经典优化，而之前的Burgers例子等逆问题就源于经典优化。

对于机器学习而言，我们通常会假设模型的最终性能将在不同的、可能未知的输入集上进行评估。测试数据通常应能捕捉到这种超出分布范围（OOD）的行为，这样我们就能对模型在 "真实世界 "中的泛化程度做出估计，而这些 "真实世界 "是我们在应用中部署模型时会遇到的。

与此相反，在本文所述的 PINN 训练中，我们在已知和给定的时空区域内重建一个单一的解。因此，来自该区域的任何样本都遵循相同的分布，因此并不能真正代表测试或 OOD 样本。由于 NN 直接对解法进行编码，因此它几乎不可能产生不同的解法，也不可能在训练范围之外表现出色。如果我们对不同的解感兴趣，就必须从头开始训练 NN。

6.2 总结

因此，通过物理软约束，我们可以利用 NN 工具对 PDE 的解进行编码。这种变式 2 的一个固有缺点是只能得到单一的解，而且不能很好地与传统的数值技术相结合。例如，学习到的表征不适合用共轭梯度法等经典迭代求解器来完善。

这意味着过去几十年中开发的许多强大技术无法在此背景下使用。让这些数值方法重新发挥作用，将是接下来几节的核心目标之一。

✅ 优点:

• 使用物理模型。
• 可以通过反向传播方便地计算导数。

❌ 缺点:

• 相当缓慢......
• 物理约束仅作为软约束施加。
• 在很大程度上与经典数值方法不兼容。
• 导数的准确性依赖于学习到的表示。

为了解决这些问题，我们接下来将研究如何利用现有的数值方法，通过使用可微分求解器来改进 DL 流程。

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

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

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

7.1 可微算子

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

假设我们有一个关于感兴趣的物理量 $\mathbf{u}(\mathbf{x},t):{\mathbb{R}}^d\times {\mathbb{R}}^{+}\to {\mathbb{R}}^d$ 的连续表达式${\mathcal{P}}^{\ast }(\mathbf{x},\nu )$，其中 $\nu$ 是模型参数（例如扩散系数、粘度或电导率）。$\mathbf{u}$ 的分量将用编号的下标表示，即 $\mathbf{u}=(u_1,u_2,\dots ,u_d{)}^T$。

通常，我们对这样一个系统的时间演化感兴趣。离散化得到一个公式 $\mathcal{P}(\mathbf{x},\nu )$，我们重新排列以计算时间步长 $\Delta t$ 后的未来状态。在 $t+\Delta t$ 时刻，通过一系列操作 ${\mathcal{P}}_1,{\mathcal{P}}_2\dots {\mathcal{P}}_m$ 来计算，使得 $\mathbf{u}(t+\Delta t)={\mathcal{P}}_m\circ \dots {\mathcal{P}}_2\circ {\mathcal{P}}_1(\mathbf{u}(t),\nu )$，其中 $\circ$ 表示函数分解，即 $f(g(x))=f\circ g(x)$。

注：为了将这个求解器整合到深度学习过程中，我们需要确保每个算子${\mathcal{P}}_i$相对于其输入提供梯度，即在上述示例中为$\partial {\mathcal{P}}_i/\partial \mathbf{u}$。

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

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

7.2 雅可比矩阵

由于 $\mathbf{u}$ 通常是一个向量值函数，所以 $\partial {\mathcal{P}}_i/\partial \mathbf{u}$ 表示雅各布矩阵 $J$，而不是单一值：

\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}

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

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

考虑到上面的表述，我们需要通过链式规则来求解 ${\mathcal{P}}_i$ 在某个当前状态 ${\mathbf{u}}^n$ 下的函数组成链的导数。例如，对于其中的两个： $$\frac{\partial ({\mathcal{P}}_1\circ {\mathcal{P}}_2)}{\partial \mathbf{u}}{|}_{{\mathbf{u}}^n}=\frac{\partial {\mathcal{P}}_1}{\partial \mathbf{u}}{|}_{{\mathcal{P}}_2({\mathbf{u}}^n)}\frac{\partial {\mathcal{P}}_2}{\partial \mathbf{u}}{|}_{{\mathbf{u}}^n}，$$

这只是“经典”链式法则 $\frac{d}{dx}f(g(x))=f^{\prime }(g(x))g^{\prime }(x)$ 的向量版本，并直接扩展到更多级联函数的情况 ($i>2$)。

这里，${\mathcal{P}}_1$和${\mathcal{P}}_2$的导数仍然是雅可比矩阵，但需要注意的是，在链的“末端”，我们有标量损失（参见Overview），最右边的雅可比矩阵将始终是一个具有1列的矩阵，即一个向量。在反向模式下，我们从该向量开始，逐个计算与左侧雅可比矩阵$\frac{\partial {\mathcal{P}}_1}{\partial \mathbf{u}}$的乘积。

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

7.3 基于 DP 算子的学习

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

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

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

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

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

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

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

7.4 一个实际例子

举个简单的例子，我们考虑将被动标量密度$d(\mathbf{x},t)$在速度场$\mathbf{u}$中的对流作为物理模型${\mathcal{P}}^{\ast }$： $$\frac{\partial d}{\partial t}+\mathbf{u}\cdot \nabla d=0$$

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

我们将这个重新表达式表示为 $\mathcal{P}$。它将 $d(t)$ 的状态映射到演化后的新状态，即:

$$d(t+\Delta t)=\mathcal{P}(d(t)，\mathbf{u}，t+\Delta t)$$

作为逆问题和学习任务的简单示例，让我们考虑找到速度场$\mathbf{u}$的问题。这个速度应该将给定的初始标量密度状态$d^0$在时间$t^0$转变为由$\mathcal{P}$演化到稍后“结束”时间$t^e$的状态，具有某种形状或配置$d^{\text{target}}$。通俗地说，我们希望找到一个流场，通过PDE模型将$d^0$变形为目标状态。表达这一目标的最简单方法是通过两种状态之间的L2损失。因此，我们想要最小化损失函数$L=|d(t^e)-d^{\text{target}}{|}^2$。

作为逆问题和学习任务的一个简单例子，我们来考虑寻找速度场 $\mathbf{u}$ 的问题。这个速度应该能将时间 $t^0$ 时的给定初始标量密度状态 $d^0$ 转变为一个状态，这个状态通过 $\mathcal{P}$ 演化到之后的 "结束 "时间 $t^e$，并具有一定的形状或配置 $d^{\text{target}}$。非正式地讲，我们希望找到一种流，通过 PDE 模型将 $d^0$ 变形为目标状态。表达这一目标的最简单方法是两个状态之间的 $L^2$ 损失。因此，我们希望最小化损失函数 $L=|d(t^e)-d^{\text{target}}{|}^2$。

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

我们的标记密度 $d(t^e)$ 的最终状态完全由通过 $\mathbf{u}$ 从 $\mathcal{P}$ 开始的演化决定，这就产生了下面的最小化问题： $${\text{argmin}}_{\mathbf{u}}|\mathcal{P}(d^0，\mathbf{u}，t^e)-d^{\text{target}}{|}^2$$

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

由于离散化的速度场 $\mathbf{u}$ 包含了所有的自由度，我们所需要做的就是通过 $\Delta \mathbf{u}=\partial L/\partial \mathbf{u}$来更新速度，它被分解为 $\Delta \mathbf{u}=\frac{\partial d}{\partial \mathbf{u}}\frac{\partial L}{\partial d}$ 。$\frac{\partial L}{\partial d}$ 分量部分通常很简单：我们会得到 $$\frac{\partial L}{\partial d}=\frac{\partial |\mathcal{P}(d^0，\mathbf{u}，t^e)-d^{\text{target}}{|}^2}{\partial d}=2(d(t^e)-d^{\text{target}})$$

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

$d$ 本身的演化是由我们的离散物理模型 $\mathcal{P}$ 给出的，我们可以互换使用 $\mathcal{P}$ 和 $d$。因此，更有趣的部分是雅可比矩阵 $\partial d/\partial \mathbf{u}=\partial \mathcal{P}/\partial \mathbf{u}$ 来计算完整的 $\Delta \mathbf{u}=\frac{\partial d}{\partial \mathbf{u}}\frac{\partial L}{\partial d}$。幸运的是，我们不需要 $\partial d/\partial \mathbf{u}$ 作为一个完整的矩阵，而只需要它与$\frac{\partial L}{\partial d}$的乘积。

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

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

下面我们将在 1 D 笛卡尔网格上使用一个简单的一阶迎风格式，其中单元 $i$ 具有标量密度 $d_i$ 和速度 $u_i$。为简洁起见， 我们省略时间 $t$ 处量的标注 $(t)$，即下面 $d_i(t)$ 简写为 $d_i$。从上面可以看到， 我们使用_物理模型_更新标量密度 $d_i(t+\Delta t)=\mathcal{P}(d_i(t)，\mathbf{u}(t)，t+\Delta t)$，其中: $$\begin{array}{rl}& d_i(t+\Delta t)=d_i-\Delta t[u_i^{+}(d_{i+1}-d_i)+u_i^{-}(d_i-d_{i-1})]\\ & \\ & u_i^{+}=\text{min}(u_i/\Delta x，0)\\ & u_i^{-}=\text{max}(u_i/\Delta x，0)\\ & \end{array}$$

因此， 对于负的 $u_i$，我们使用 $u_i^{+}$ 沿速度相反的方向查看， 即运动学意义上的_后向_。在这种情况下 $u_i^{-}$ 将为零。对于正的 $u_i$ 则反之， 我们将得到一个零值的 $u_i^{+}$，并通过 $u_i^{-}$ 得到一个后向差分模板。为了选择前一种情况， 对于负的 $u_i$ 我们得到: $$\mathcal{P}(d_i(t)，\mathbf{u}(t)，t+\Delta t)=(1+\frac{u_i\Delta t}{\Delta x})d_i-\frac{u_i\Delta t}{\Delta x}{d}_{i+1}$$

因此$\partial \mathcal{P}/\partial u_i$ 给出了 $\frac{\Delta t}{\Delta x}{d}_i-\frac{\Delta t}{\Delta x}{d}_{i+1}$。直观地看，速度$u_i$的变化取决于密度的空间导数。由于采用了一阶迎风格式，我们只考虑了两个相邻点(高阶方法会依赖$d$的更多项)。

实际上，这一步等同于计算转置矩阵乘法。如果我们把上面的计算改写成 $\mathcal{P}(d_i(t),\mathbf{u}(t),t+\Delta t)=A\mathbf{u}$，那么 $(\partial \mathcal{P}/\partial \mathbf{u}{)}^T=A^T$。然而，在许多实际情况下，这种乘法的无矩阵实现可能比实际构造 $A$ 更可取。

对于对流格式，我们可以考虑的另一个导数是与之前密度状态相关的导数，即 $d_i(t)$，简写为 $d_i$。对于单元 $i$，从上面的$\partial \mathcal{P}/\partial d_i$得出$1+\frac{u_i\Delta t}{\Delta x}$。然而，为了得到完整的梯度，我们需要加上来自 $i+1$ 和 $i-1$ 单元的潜在贡献，这取决于它们速度的符号。这种导数将在下一节中发挥作用。

7.4.2 时间演化

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

在上面$d_i(t+\Delta t)$的表达式中，每个$d_i(t)$都依次取决于$t-\Delta t$时的速度和密度状态，即$d_i(t-\Delta t)$。因此，我们必须追溯损失 $L$ 的影响，一直追溯到 $\mathbf{u}$ 如何影响初始标记状态。这可能需要通过 $\mathcal{P}$ 对对流格式进行大量评估。

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

给定上述机制， 反向追踪非常简单: 对于 $\mathcal{P}$ 中的每个对流步骤， 我们计算雅可比矩阵与来自损失 $L$ 或前一个对流步骤的导数的_输入_向量的乘积。我们重复这个过程直到追踪链条从具有 $d^{\text{target}}$ 的损失一直回溯到 $d^0$。从理论上讲， 速度 $\mathbf{u}$ 可以像 $d$ 一样是时间的函数， 在这种情况下， 我们在每个时间步 $t$ 都会得到一个梯度 $\Delta \mathbf{u}(t)$。然而， 为了简化下面的内容， 我们假设我们有一个随时间不变的场， 即我们对每个通过 $\mathcal{P}$ 的对流都重用相同的速度 $\mathbf{u}$。现在， 每个时间步骤都会为我们提供一个对 $\Delta \mathbf{u}$ 的贡献，我们对所有步骤进行累积。 $$\begin{array}{rl}\Delta \mathbf{u}=& \frac{\partial d(t^e)}{\partial \mathbf{u}}\frac{\partial L}{\partial d(t^e)}\\ & +\frac{\partial d(t^e-\Delta t)}{\partial \mathbf{u}}\frac{\partial d(t^e)}{\partial d(t^e-\Delta t)}\frac{\partial L}{\partial d(t^e)}\\ & +\cdots \\ & +(\frac{\partial d(t^0)}{\partial \mathbf{u}}\cdots \frac{\partial d(t^e-\Delta t)}{\partial d(t^e-2\Delta t)}\frac{\partial d(t^e)}{\partial d(t^e-\Delta t)}\frac{\partial L}{\partial d(t^e)})\end{array}$$