第二十八周周报
- 一、文献阅读
- 题目信息
- 摘要
- Abstract
- 网络架构
- 实验——Data-driven discovery of partial differential equations（偏微分方程的数据驱动发现）
- 1. Continuous time models（连续时间模型）
- 例子：(Navier–Stokes equation)
- 2. Discrete time models（离散时间模型）
- 例子：Korteweg–de Vries equation（KdV方程）
- 结论
- 缺点以及后续展望
- 二、动手深度学习
- 自动微分
- 从零实现线性回归
一、文献阅读
题目信息
- 题目：《Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations》
- 期刊： Journal of Computational Physics
- 作者： M. Raissi、P. Perdikaris、G.E. Karniadakis
- 发表时间： 2019
- 文章链接：https://www.sciencedirect.com/science/article/pii/S0021999118307125
摘要
文献主要讲述了作者是如何利用物理信息神经网络（Physics-informed neural networks，PINN）正反解求偏微分方程（PDEs）。首先作者提出在数据量稀缺背景下，多数机器学习技术缺乏鲁棒性且无法保证收敛，而利用PINN求解相关的偏微分方程可以解决这类问题。。因此，作者在论文中详细讲述了深度学习与物理规定律相结合将物理定律融入神经网络的损失函数中，并提出偏微分方程的数据驱动求解和数据驱动发现两类主要问题展开研究。在研究中，作者根据数据的性质设计了连续时间和离散时间模型，其中，连续时间模型使用一般形式的偏微分方式，离散时间模型使用Runge–Kutta方法。然后，作者通过流体、量子力学、非线性浅水波传播等实验证明了PINN的有效性。最后，作者也提出了PINN的不足，例如算法收敛性、梯度消失、泛化能力不理想等问题，后续需要进一步研究以及解决相关问题。
Abstract
The manuscript delineates the utilization of Physics-informed Neural Networks (PINNs) for both forward and inverse solutions of partial differential equations (PDEs). The author commences by positing that in the context of scarce data, the majority of machine learning techniques exhibit a lack of robustness and convergence guarantees, which can be circumvented by employing PINNs to address PDEs. Consequently, the paper provides an exhaustive exposition on the integration of deep learning with physical laws, embedding these laws into the loss functions of neural networks, and explores two principal issues: data-driven solutions and data-driven discovery of PDEs. In the research, the author designs models for both continuous and discrete time, with the continuous time model employing a general form of partial differential equations and the discrete time model leveraging the Runge-Kutta method. Subsequently, the efficacy of PINNs is substantiated through experiments in fluid dynamics, quantum mechanics, and the propagation of nonlinear shallow water waves. Ultimately, the author acknowledges the limitations of PINNs, including issues with algorithmic convergence, gradient vanishing, and suboptimal generalization capabilities, which necessitate further research and resolution.
网络架构
下图展示论文中作者设计的PINNs的架构
主要结构如下：
输入层: 输入数据 x（空间空间） 和 t （时间变量）通过神经网络(NN)。
N
(
x
,
t
,
θ
)
\mathcal{N}(\mathbf{x}, t, \theta)
N(x,t,θ)是一个神经网络模型，
θ
\theta
θ 表示网络的参数。
输出层： 网络的输出，是微分方程的解u，这是PINN的预测解，用于损失函数的优化。其中求导同过自动微分技术(AD)完成。
损失函数：: 模型的损失函数如下：
L
(
θ
)
=
L
data
+
L
bic
+
L
p
\mathcal{L}(\theta) = \mathcal{L}_{\text{data}} +\mathcal{L}_{\text{bic}} +\mathcal{L}_{\text{p}}
L(θ)=Ldata+Lbic+Lp
其由三部分组成：
Labeled Data Loss(预测值与真实值之间的损失) :
L
data
\mathcal{L}_{\text{data}}
Ldata 是处理有标签的数据，即已知的解。预测值与真实值之间的误差。
Boundary/Initial Conditions Loss（初始条件、边界条件损失） :
L
bic
\mathcal{L}_{\text{bic}}
Lbic用于确保网络输出满足边界条件或初始条件。
边界或初始条件是从边界或初始时间点上计算输出与预期条件之间的误差来实现。
方程的边界条件和初始条件是PINN求解偏微分方程（PDEs）问题的关键组成部分。
其中，初始条件被用来约束网络在初始时间点的输出，确保它与给定的初始状态相匹配。边界条件被用来约束网络的输出，确保它在边界上满足方程的约束。
Residual on PDE Equations Loss（残差损失）:
L
p
\mathcal{L}_{\text{p}}
Lp 也叫做残差项，它确保网络输出满足物理定律，通过计算网络输出在PDE中的残差来实现。
实验——Data-driven discovery of partial differential equations（偏微分方程的数据驱动发现）
作者在论文中的实验均已开源
github网址如下：https://github.com/maziarraissi/PINNs
上一周的周报我们已经解读了论文在解决正问题（即在指定的初始和边界条件下求解PDEs）时的方法。
其中
作者对连续时间模型提出一般形式的微分方程：
u
t
+
N
[
u
]
=
0
,
x
∈
Ω
,
t
∈
[
0
,
T
]
u_{t}+\mathcal{N}[u]=0, x \in \Omega, t \in[0, T]
ut+N[u]=0,x∈Ω,t∈[0,T]
作者对离散时间模型提出Runge-Kutta方法：
u
n
+
c
i
=
u
n
−
Δ
t
∑
j
=
1
q
a
i
j
N
[
u
n
+
c
j
]
,
i
=
1
,
…
,
q
u
n
+
1
=
u
n
−
Δ
t
∑
j
=
1
q
b
j
N
[
u
n
+
c
j
]
\begin{array}{l}u^{n+c_i} = u^n - \Delta t \sum_{j=1}^{q} a_{ij} \mathcal{N}[u^{n+c_j}], \quad i = 1, \ldots, q \\ u^{n+1} = u^n - \Delta t \sum_{j=1}^{q} b_j \mathcal{N}[u^{n+c_j}] \end{array}
un+ci=un−Δt∑j=1qaijN[un+cj],i=1,…,qun+1=un−Δt∑j=1qbjN[un+cj]
此外，作者还分别举Schrödinger方程与Allen-Cahn方程来验证PINN求解正问题的能力，并且都得到了不错的效果。
所以这一周我们来详细研究一下作者是如何使用PINNs解决反问题（即确定未知的参数、边界条件或偏微分方程本身的问题）。反问题通常是指使用PINN根据已知的解u来确定模型的未知条件，这些条件能包括模型的参数或边界条件，从而构建一个表现较好的模型。
简单来说，正问题是从已知条件出发找到结果，而反问题是从已知结果出发，推测出可能的条件。
作者同样使用了连续时间模型和离散时间模型进行实验验证PINN在反问题中的表现性能。
1. Continuous time models（连续时间模型）
在连续时间模型中，作者依旧使用的是一般形式的微分方程来解决。
例子：(Navier–Stokes equation)
纳维-斯托克斯方程描述了牛顿流体的动量守恒和质量守恒，并且考虑了压强、温度、密度等相关的状态方程。
纳维-斯托克斯方程的推导源于牛顿第二定律，假设流体中的应力是扩散粘性项和压力项之和，于是描述了粘性流。
它可以用来模拟天气、洋流、机翼附近的气流等流体力学体系，在工程中有极为重要的应用。
二维的Navier-Stokes方程如下：
u
t
+
λ
1
(
u
u
x
+
v
u
y
)
=
−
p
x
+
λ
2
(
u
x
x
+
u
y
y
)
,
v
t
+
λ
1
(
u
v
x
+
v
v
y
)
=
−
p
y
+
λ
2
(
v
x
x
+
v
y
y
)
，
\begin{array}{l}u_{t}+\lambda_{1}\left(u u_{x}+v u_{y}\right)=-p_{x}+\lambda_{2}\left(u_{x x}+u_{y y}\right), \\v_{t}+\lambda_{1}\left(u v_{x}+v v_{y}\right)=-p_{y}+\lambda_{2}\left(v_{x x}+v_{y y}\right)，\end{array}
ut+λ1(uux+vuy)=−px+λ2(uxx+uyy),vt+λ1(uvx+vvy)=−py+λ2(vxx+vyy)，
其中
u
(
t
,
x
,
y
)
u(t, x, y)
u(t,x,y)表示速度场的x分量，
v
(
t
,
x
,
y
)
v(t, x, y)
v(t,x,y)表示y分量，
p
(
t
,
x
,
y
)
p(t, x, y)
p(t,x,y)表示压力。
λ
=
(
λ
1
，
λ
2
)
λ= (λ 1， λ2)
λ=(λ1，λ2)为未知参数。
作者假设Navier-Stokes方程的解，即
u
x
+
u
y
=
0
u_x + u_y = 0
ux+uy=0，这个方程述了流体的质量守恒。
反问题就表现在作者通过对已知的
{
t
i
,
x
i
,
y
i
,
u
i
,
v
i
}
i
=
1
N
\left\{t^{i}, x^{i}, y^{i}, u^{i}, v^{i}\right\}_{i=1}^{N}
{ti,xi,yi,ui,vi}i=1N的测量值，通过PINN确定参数
λ
λ
λ 以及压力
p
(
t
,
x
,
y
)
p(t, x, y)
p(t,x,y)的最佳值。
作者将定义残差项为： f : = u t + λ 1 ( u u x + v u y ) + p x − λ 2 ( u x x + u y y ) , g : = v t + λ 1 ( u v x + v v y ) + p y − λ 2 ( v x x + v y y ) , \begin{array}{l}f:=u_{t}+\lambda_{1}\left(u u_{x}+v u_{y}\right)+p_{x}-\lambda_{2}\left(u_{x x}+u_{y y}\right), \\g:=v_{t}+\lambda_{1}\left(u v_{x}+v v_{y}\right)+p_{y}-\lambda_{2}\left(v_{x x}+v_{y y}\right),\end{array} f:=ut+λ1(uux+vuy)+px−λ2(uxx+uyy),g:=vt+λ1(uvx+vvy)+py−λ2(vxx+vyy),（其实就是把二维的Navier-Stokes方程移项得到的残差项，目的就是为了使f与g接近0）
所以将损失函数定义为：
M
S
E
:
=
1
N
∑
i
=
1
N
(
∣
u
(
t
i
,
x
i
,
y
i
)
−
u
i
∣
2
+
∣
v
(
t
i
,
x
i
,
y
i
)
−
v
i
∣
2
)
+
1
N
∑
i
=
1
N
(
∣
f
(
t
i
,
x
i
,
y
i
)
∣
2
+
∣
g
(
t
i
,
x
i
,
y
i
)
∣
2
)
MSE:=\frac{1}{N} \sum_{i=1}^{N}\left(\left|u\left(t^{i}, x^{i}, y^{i}\right)-u^{i}\right|^{2}+\left|v\left(t^{i}, x^{i}, y^{i}\right)-v^{i}\right|^{2}\right)+\frac{1}{N} \sum_{i=1}^{N}\left(\left|f\left(t^{i}, x^{i}, y^{i}\right)\right|^{2}+\left|g\left(t^{i}, x^{i}, y^{i}\right)\right|^{2}\right)
MSE:=N1i=1∑N(
u(ti,xi,yi)−ui
2+
v(ti,xi,yi)−vi
2)+N1i=1∑N(
f(ti,xi,yi)
2+
g(ti,xi,yi)
2)
前半部分是计算测量值与真实值之间的损失，后半部分是残差损失，起到物理约束的作用。
在实验中作者以圆柱绕流问题为例，用谱/hp - 元求解器生成高分辨率数据集。
作者随机抽取5000个数据点（占总数据1%）用于训练9层每层20个神经元的神经网络。
实验结果如下：
无噪声训练数据时，估计参数λ₁和λ₂的误差分别为0.078%和4.67%；
即使训练数据被1%的不相关高斯噪声破坏，误差分别为0.17%和5.70%，且在无压力训练数据情况下能定性准确预测压力场。
2. Discrete time models（离散时间模型）
例子：Korteweg–de Vries equation（KdV方程）
KdV方程研究浅水中小振幅长波运动时共同发现的一种单向运动浅水波偏微分方程
它描述了长一维波在许多物理环境中的演化，包括具有弱非线性恢复力的浅水波、密度分层海洋中的长内波、等离子体中的离子声波和晶格上的声波。
作者通过举这个例子强调了作者提出的PINN处理涉及高阶导数的计算偏微分方程的能力。
KdV方程如下：
u
t
+
λ
1
u
u
x
+
λ
2
u
x
x
x
=
0
u_{t}+\lambda_{1} u u_{x}+\lambda_{2} u_{x x x}=0
ut+λ1uux+λ2uxxx=0
其中
(
λ
1
,
λ
2
)
(λ 1, λ2)
(λ1,λ2)为未知参数。
结合作者提出离散时间模型的Runge-Kutta方法，得到如下方程：
N
[
u
n
+
c
j
]
=
λ
1
u
n
+
c
j
u
x
n
+
c
j
−
λ
2
u
x
x
x
n
+
c
j
\mathcal{N}\left[u^{n+c_{j}}\right]=\lambda_{1} u^{n+c_{j}} u_{x}^{n+c_{j}}-\lambda_{2} u_{x x x}^{n+c_{j}}
N[un+cj]=λ1un+cjuxn+cj−λ2uxxxn+cj
神经网络的共享参数（如下图中式子(23)、式子(24)和式子(25)）以及KdV方程的参数
λ
=
(
λ
1
，
λ
2
)
λ= (λ 1， λ2)
λ=(λ1，λ2) 可以通过最小化损失函数（下图中式子(26)）来学习。
作者用传统谱方法模拟得到数据集，提取两个时刻的解快照并随机抽样生成训练数据集。
在网络架构设计方面，作者采用4层每层50个神经元的离散时间物理信息神经网络训练。
实验结果如下：
在无噪声训练数据时，估计参数λ₁和λ₂的误差分别为0.023%和0.006%。
含1%噪声时，误差分别为0.057%和0.017%。
实验结果表明，无论训练数据是否被噪声破坏，PINN的方法都能够正确识别未知参数。
结论
作者在论文中提出了用物理信息神经网络（PINN）来处理数据和求解偏微分方程。PINN的核心在于其把物理定律嵌入到了神经网络架构中。经过作者的实验验证了这种方法在很多计算科学的问题上都显示出了很好的效果，给深度学习和数学物理的结合提供了新的研究方向。但是，其还是不能够取代那些传统的数值方法。相反，这两种方法可以互相补充，取长补短来共同解决复杂的PDEs正反问题。
缺点以及后续展望
作者在文章的最后还提出了PINNs的不足：神经网络的构造方面，比如要设置多少层和需要多少数据来训练才能让模型的表现最佳这些问题还没有得到解决；模型训练中产生的问题，如梯度消失、过拟合等问题依旧需要研究解决；PINN的初始设置，比如开始时权重怎么给定和数据的处理方式也会影响结果，影响模型的泛化能力。
针对这些缺陷，作者提出了未来研究PINN求解微分方程的一些展望，例如，需要进一步探索神经网络结构与数据量的关系和研究神经网络的结构和需要的数据量之间的关系。改进网络开始时的设置和数据的处理方式；找到更好的方法来衡量模型预测的准确性等问题，这都是值得我们作为PINN后续的方向以及需要解答的问题。
二、动手深度学习
自动微分
因为后续要进行PINN代码相关的实验复现，所以需要了解自动微分（automatic differentiation，AD）在pytorch是怎么实现的。
深度学习框架通过⾃动计算导数，即⾃动微分（automatic differentiation）来加快求导。
假设我们想对函数 y = 2 x ⊤ x y = 2x^⊤x y=2x⊤x关于列向量x求导。
- 我们先创建变量x并为其分配⼀个初始值。
import torch
x = torch.arange(4.0)
x
- 在我们计算y关于x的梯度之前，我们需要⼀个地方来存储梯度。
因为不会在每次对⼀个参数求导时都分配新的内存，又因为训练中经常会成千上万次地更新相同的参数，而每次都分配新的内存可能很快就会将内存耗尽。所以我们将梯度存储在x.grad中。
x.requires_grad_(True) # 等价于x=torch.arange(4.0,requires_grad=True)
x.grad # 默认值是None
- 然后使用点积
torch.dot()计算y，x是⼀个⻓度为4的向量，计算x和x的点积，得到了我们赋值给y的标量输出。
y = 2 * torch.dot(x, x) # dot是点积
y
- 接下来，通过反向传播函数来自动计算y关于x每个分量的梯度，并展示梯度
[0., 1., 2., 3.] *4 = [ 0., 4., 8., 12.]
y.backward()
x.grad
- 函数
y
=
2
x
⊤
x
y = 2x^⊤x
y=2x⊤x关于x的梯度应为4x(因为
2
x
2
2x^2
2x2对x求导为
4
x
4x
4x)。
验证这个梯度是否计算正确。
x.grad == 4 * x
- 然后计算x的另⼀个函数，需要注意将梯度清零。此时 y = ∑ i = 0 x y = \sum_{i=0}{x} y=∑i=0x
# 在默认情况下，PyTorch会累积梯度，我们需要清除之前的值
x.grad.zero_() #下划线代表重写内容
y = x.sum()
y.backward()
x.grad
使用动微分的⼀个好处是，即使构建函数的计算图需要通过Python控制流（例如，条件、循环或任意函数调用），可以计算得到的变量的梯度。
while循环的次数和if语句的结果都取决于输⼊a的值。
def f(a):
b = a * 2
while b.norm() < 1000:
b = b * 2
if b.sum() > 0:
c = b
else:
c = 100 * b
return c
计算梯度
a = torch.randn(size=(), requires_grad=True)
d = f(a)
d.backward()
我们现在可以分析上⾯定义的f函数。对于任何a，存在某个常量标量k，使得
f
(
a
)
=
k
∗
a
f(a)=k*a
f(a)=k∗a，其中k的值取决于输⼊a。
因此，我们可以⽤d/a验证梯度是否正确。
a.grad == d / a
从零实现线性回归
- ⽣成⼀个包含1000个样本的数据集，每个样本包含从标准正态分布中采样的2个特征。我们的合成数据集是
⼀个矩阵 X ∈ R 1000 × 2 X ∈ R^{1000×2} X∈R1000×2。
使⽤线性模型参数 w = [ 2 , − 3.4 ] ⊤、 b = 4.2 w = [2, −3.4]⊤、b = 4.2 w=[2,−3.4]⊤、b=4.2 和噪声项 ϵ ϵ ϵ⽣成数据集及其标签：
y = X w + b + ϵ . y = Xw + b + ϵ. y=Xw+b+ϵ.
将 ϵ ϵ ϵ视为模型预测和标签时的潜在观测误差。将标准差设为0.01。
%matplotlib inline
import random
import torch
from d2l import torch as d2l
def synthetic_data(w, b, num_examples): #@save
X = torch.normal(0, 1, (num_examples, len(w)))
y = torch.matmul(X, w) + b
y += torch.normal(0, 0.01, y.shape)
return X, y.reshape((-1, 1))
true_w = torch.tensor([2, -3.4])
true_b = 4.2
# features中的每⼀⾏都包含⼀个⼆维数据样本，labels中的每⼀⾏都包含⼀维标签值。
features, labels = synthetic_data(true_w, true_b, 1000)
print('features:', features[0],'\nlabel:', labels[0])
# 通过⽣成第⼆个特征features[:, 1]和labels的散点图，可以直观观察到两者之间的线性关系。
d2l.set_figsize()
d2l.plt.scatter(features[:, (1)].detach().numpy(), labels.detach().numpy(), 1); #散点大小,不加默认为20
- 定义⼀个函数，该函数能打乱数据集中的样本并以⼩批量⽅式获取数据。
- 在下⾯的代码中，我们定义⼀个data_iter函数，该函数接收批量⼤⼩、特征矩阵和标签向量作为输⼊，⽣成⼤⼩为batch_size的⼩批量。每个⼩批量包含⼀组特征和标签。
def data_iter(batch_size, features, labels):
num_examples = len(features)
indices = list(range(num_examples))
# 这些样本是随机读取的，没有特定的顺序
# 随机排序列表
random.shuffle(indices)
for i in range(0, num_examples, batch_size):
batch_indices = torch.tensor(
# num_examples防止超出列表
indices[i: min(i + batch_size, num_examples)])
#带yield的函数是一个生成器,可循环返回数据。它记住上一次返回时在函数体中的位置。对生成器函数的第二次(或第 n 次)调用跳转至该函数中间,而上次调用的所有局部变量都保持不变
yield features[batch_indices], labels[batch_indices]
#读取第⼀个⼩批量数据样本并打印。每个批量的特征维度显⽰批量⼤⼩和输⼊特征数。同样的，批量的标签形状与batch_size相等。
batch_size = 10
for X, y in data_iter(batch_size, features, labels):
print(X, '\n', y)
break
- 开始⽤⼩批量随机梯度下降优化我们的模型参数之前，我们需要先有⼀些参数。 下⾯的代码中，我们通过从均值为0、标准差为0.01的正态分布中采样随机数来初始化权重，并将偏置初始化为0。
w = torch.normal(0, 0.01, size=(2,1), requires_grad=True)
b = torch.zeros(1, requires_grad=True)
# 定义模型
def linreg(X, w, b): #@save
return torch.matmul(X, w) + b
# 因为需要计算损失函数的梯度，所以我们应该先定义损失函数。
def squared_loss(y_hat, y): #@save
"""均⽅损失"""
return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2
#在每⼀步中，使⽤从数据集中随机抽取的⼀个⼩批量，然后根据参数计算损失的梯度。接下来，朝着减少损失的⽅向更新我们的参数。下⾯的函数实现⼩批量随机梯度下降更新。该函数接受模型参数集合、学习速率和批量⼤⼩作为输⼊。每⼀步更新的⼤⼩由学习速率lr决定。因为我们计算的损失是⼀个批量样本的总和，所以我们⽤批量⼤⼩（batch_size）来规范化步⻓，这样步⻓⼤⼩就不会取决于我们对批量⼤⼩的选择。
def sgd(params, lr, batch_size): #@save
"""⼩批量随机梯度下降"""
with torch.no_grad():
for param in params:
param -= lr * param.grad / batch_size
param.grad.zero_()
# 在每个迭代周期（epoch）中，我们使⽤data_iter函数遍历整个数据集，并将训练数据集中所有样本都使⽤⼀次（假设样本数能够被批量⼤⼩整除）。这⾥的迭代周期个数num_epochs和学习率lr都是超参数，分别设为3和0.03。
lr = 0.03
num_epochs = 3
net = linreg
loss = squared_loss
注意：
-
with 语句适用于对资源进行访问的场合,确保不管使用过程中是否发生异常都会执行必要的“清理”操作,释放资源,比如文件使用后自动关闭/线程中锁的自动获取和释放等。
(1)紧跟with后面的语句被求值后,返回对象的“–enter–()”方法被调用,这个方法的返回值将被赋值给as后面的变量;
(2)当with后面的代码块全部被执行完之后,将调用前面返回对象的“–exit–()”方法。 -
with torch.no_grad的作用在该模块下,所有计算得出的tensor的requires_grad都自动设置为False。
for epoch in range(num_epochs):
for X, y in data_iter(batch_size, features, labels):
l = loss(net(X, w, b), y) # X和y的⼩批量损失
# 因为l形状是(batch_size,1)，⽽不是⼀个标量。l中的所有元素被加到⼀起，
# 并以此计算关于[w,b]的梯度
l.sum().backward()
sgd([w, b], lr, batch_size) # 使⽤参数的梯度更新参数
with torch.no_grad():
train_l = loss(net(features, w, b), labels)
print(f'epoch {epoch + 1}, loss {float(train_l.mean()):f}')
输出结果：
print(f'w的估计误差: {true_w - w.reshape(true_w.shape)}')
print(f'b的估计误差: {true_b - b}')
