论文阅读:Achieving High Accuracy with PINNs via Energy Natural Gradients
Achieving High Accuracy with PINNs via Energy Natural Gradients
问题
据观察,偏微分方程残差、边界项和初始条件对梯度贡献的大小常常具有不平衡的大小。为了解决这个问题,针对损失的各个组成部分制定了不同的加权策略(Wang et al., 2022)。尽管改进了 PINN 训练,但上述工作均未报告相对 L2 误差低于 1 0 − 4 10^{−4} 10−4 。
PINN 损失离散化中配置点的选择已在多种工作中进行了研究(Lu et al., 2021; Wang et al., 2022a;Wu et al., 2023)。所有这些研究的共同点是观察到搭配点应集中在高 PDE 残差区域,本文参考(Wu et al., 2023)对不同提出的抽样策略进行了广泛比较。此外,据报道,对于时间相关问题,课程学习可以减轻与解决长期进化问题相关的训练病态(Wang et al., 2022a;Krishnapriyan et al., 2021)。同样,虽然所有上述工作都显着改善了 PINN 训练,但没有一个贡献误差可以达到低于 1 0 − 4 10^{−4} 10−4 。
在 PINN 的背景下,人们提出了不同的优化策略,这些策略在概念上不同于基于直接梯度的目标优化。例如,贪婪算法用于逐个神经元增量构建浅层神经神经元,这对于各种偏微分方程具有高精度,相对误差高达 1 0 − 8 10^{−8} 10−8(Hao 等人,2021)。然而,所提出的贪婪算法仅在计算上易于处理浅层神经网络。另一个设想是将二次 PINN 损失重新表述为鞍点问题,涉及用于逼近解的网络和惩罚非零残差的判别器网络。由此产生的鞍点公式可以通过竞争梯度下降来解决(Zeng et al., 2022),作者报告了高度准确的——高达 1 0 − 8 10^{−8} 10−8 的相对 L2 误差——针对许多示例问题的 PINN 解决方案。然而,这种方法的代价是训练两个神经网络并将最小化问题替换为鞍点问题。最后,粒子群优化方法在 PINN 的背景下提出,它们提高了标准优化器的准确性,但尽管有计算负担,但未能实现优于 1 0 − 3 10^{−3} 10−3 的准确性。
自然梯度法是一种成熟的优化算法,这里仅讨论与偏微分方程数值解相关的工作。事实上,在没有明确参考自然梯度文献和术语的情况下,自然梯度在有限元背景下的 PDE 约束优化中广泛使用。例如,在某些情况下,质量或刚度矩阵可以解释为格拉米矩阵,表明该方法确实是一种自然梯度方法。对于明确的例子,可以参考(Schwedes et al., 2016; 2017)。在基于神经网络的方法的背景下,由 Sobolev、Fisher-Rao 和 Wasserstein 几何形状引起的各种自然梯度已被提出并针对 PINN 进行了测试(Nurbekyan et al., 2022)。这项工作的重点是这些方法的有效实现,没有考虑基于能量的自然梯度,本文发现这是实现高精度所必需的。
本文用
L
p
(
Ω
)
L^{p}(\Omega)
Lp(Ω) 表示
Ω
⊆
R
d
\Omega\subseteq\mathbb{R}^d
Ω⊆Rd 上可积
p
p
p 次方的函数空间,并赋予它规范范数。对于足够平滑的函数
u
u
u,本文用
∂
i
u
=
∂
u
/
∂
x
i
\partial_iu=\partial u/\partial x_i
∂iu=∂u/∂xi 表示其偏导数,并用
(
D
l
u
)
i
1
,
.
.
.
,
i
l
:
=
∂
i
1
…
∂
i
l
u
(D^lu)_{i_1,...,i_l}:=\partial_{i_1}\ldots\partial_{i_l}u
(Dlu)i1,...,il:=∂i1…∂ilu 表示与
l
l
l 阶导数相关的张量。 本文用 $\nabla u=(\partial_1u,\ldots,\partial_du)^\top $ 表示足够平滑函数
u
u
u 的梯度,拉普拉斯算子
Δ
\Delta
Δ 定义为
Δ
u
:
=
∑
i
=
1
d
∂
i
2
u
\Delta u:=\sum_{i=1}^d\partial_i^2u
Δu:=∑i=1d∂i2u。本文将
L
p
(
Ω
)
L^{p}(\Omega)
Lp(Ω) 中具有
k
k
k 阶弱导数的函数的 Sobolev 空间表示为
W
k
,
p
(
Ω
)
W^{k,p}(\Omega)
Wk,p(Ω),这是一个具有如下范数的 Banach 空间:
∥
u
∥
W
k
,
p
(
Ω
)
p
:
=
∑
l
=
0
k
∥
D
l
u
∥
L
p
(
Ω
)
p
.
\|u\|_{W^{k,p}(\Omega)}^p:=\sum_{l=0}^k\|D^lu\|_{L^p(\Omega)}^p.
∥u∥Wk,p(Ω)p:=l=0∑k∥Dlu∥Lp(Ω)p.
下面本文主要处理
p
=
2
p = 2
p=2 的情况,并写为
H
k
(
Ω
)
H^{k}(\Omega)
Hk(Ω) 而不是
W
k
,
2
(
Ω
)
W^{k,2}(\Omega)
Wk,2(Ω)。
考虑自然数
d
、
m
、
L
、
N
0
,
.
.
.
,
N
L
d、m、L、N_0,..., N_L
d、m、L、N0,...,NL 并设
θ
=
(
(
A
1
,
b
1
)
,
.
.
.
,
(
A
L
,
b
L
)
)
\theta = ((A_1, b_1), ..., (A_L, b_L))
θ=((A1,b1),...,(AL,bL)) 为矩阵向量对的元组,其中
A
l
∈
R
N
l
×
N
l
−
1
,
b
l
∈
R
N
l
A_{l}\in\mathbb{R}^{N_{l}\times N_{l-1}} , b_l \in \mathbb{R}^{N_{l}}
Al∈RNl×Nl−1,bl∈RNl 且
N
0
=
d
,
N
L
=
m
N_0 = d, N_L =m
N0=d,NL=m 。每个矩阵向量对
(
A
l
,
b
l
)
(A_l, b_l)
(Al,bl) 都会产生一个仿射线性映射
T
l
:
R
N
l
−
1
→
R
N
l
T_l:\mathbb{R}^{N_{l-1}}\to\mathbb{R}^{N_l}
Tl:RNl−1→RNl 。一个具有参数
θ
\theta
θ 且相对于某个激活函数
ρ
:
R
→
R
\rho\colon\mathbb{R}\to\mathbb{R}
ρ:R→R 的神经网络函数可表示如下:
u
θ
:
R
d
→
R
m
,
x
↦
T
L
(
ρ
(
T
L
−
1
(
ρ
(
⋯
ρ
(
T
1
(
x
)
)
)
)
)
)
)
u_{\boldsymbol{\theta}}:\mathbb{R}^{\boldsymbol{d}}\to\mathbb{R}^{\boldsymbol{m}},\quad x\mapsto T_{L}(\rho(T_{L-1}(\rho(\cdots\rho(T_{1}(x)))))))
uθ:Rd→Rm,x↦TL(ρ(TL−1(ρ(⋯ρ(T1(x)))))))
这种网络的参数数量和神经元数量可以由
∑
l
=
0
L
−
1
(
n
l
+
1
)
n
l
+
1
\sum_{l=0}^{L-1}(n_l+1){n_{l+1}}
∑l=0L−1(nl+1)nl+1 给出。如果网络深度为 2,本文称其为浅网络,否则称为深网络。在其余部分,本文将限制在
m
=
1
m = 1
m=1 的情况,因为本文只考虑实值函数。进一步,在本文的实验中,选择 tanh 作为激活函数,以便假设网络函数
u
θ
u_\theta
uθ 和参数化
θ
↦
u
θ
\theta\mapsto u_\theta
θ↦uθ 所需的平滑度概念。对于
A
∈
R
n
×
m
A\in\mathbb{R}^{n\times m}
A∈Rn×m,本文使用
A
+
A^+
A+ 表示 A 的任何伪逆。
自然梯度
现有方法总结
人们提出了各种基于神经网络的偏微分方程近似解方法(Beck et al., 2020;Weinan et al., 2021;Kovachki et al., 2021)。其中大多数将偏微分方程的解视为某些函数空间上典型凸能量的最小值,并使用该能量来优化网络参数。本文提出两种重要的方法,并介绍稍后用于处理这两种方法的统一设置。
-
PINN
考虑以下形式的一般偏微分方程:
L u = f i n Ω B u = g o n ∂ Ω , \begin{aligned} &\mathcal{L}u=f\quad\mathrm{in~}\Omega\\ &\mathcal{B}u=g\quad\mathrm{on~}\partial\Omega, \end{aligned} Lu=fin ΩBu=gon ∂Ω,
其中 Ω ⊆ R d \Omega\subseteq\mathbb{R}^d Ω⊆Rd 是一个开集, L \mathcal{L} L 是一个(可能是非线性的)偏微分算子, B \mathcal{B} B 是一个边界值算子。假设解 u u u 在希尔伯特空间 X X X 中求得,并且右侧 f f f 和边界值 g g g 分别是 Ω \Omega Ω 和 ∂ Ω \partial \Omega ∂Ω 上的平方可积函数。在这种情况下,可以将上式重新表述为具有如下目标函数的最小化问题
E ( u ) = ∫ Ω ( L u − f ) 2 d x + τ ∫ ∂ Ω ( B u − g ) 2 d s E(u)=\int_{\Omega}(\mathcal{L}u-f)^2\mathrm{d}x+\tau\int_{\partial\Omega}(\mathcal{B}u-g)^2\mathrm{d}s E(u)=∫Ω(Lu−f)2dx+τ∫∂Ω(Bu−g)2ds
对于惩罚参数 τ > 0 \tau > 0 τ>0。当且仅当 E ( u ) = 0 E(u) = 0 E(u)=0 时,函数 u ∈ X u \in X u∈X 满足偏微分方程。为了获得近似解,可以通过神经网络参数化函数 u θ u_\theta uθ 并根据损失函数最小化网络参数 θ ∈ R p \theta\in\mathbb{R}^p θ∈Rp
L ( θ ) : = ∫ Ω ( L u θ − f ) 2 d x + τ ∫ ∂ Ω ( B u θ − g ) 2 d s L(\theta):=\int_{\Omega}(\mathcal{L}u_{\theta}-f)^2\mathrm{d}x+\tau\int_{\partial\Omega}(\mathcal{B}u_{\theta}-g)^2\mathrm{d}s L(θ):=∫Ω(Luθ−f)2dx+τ∫∂Ω(Buθ−g)2ds
这种将方程表述为最小化问题的通用方法称为残差最小化,在偏微分方程神经网络的背景下可以追溯到(Dissanayake & Phan-Thien,1994;Lagaris et al., 1998)。最近,这种方法以深度 Galerkin 方法或物理信息神经网络的名称得到普及,其中损失也可以被增强,以包含源自解决方案的现实世界测量的回归项(Sirignano & Spiliopoulos,2018;Raissi 等人) .,2019)。在实践中,目标函数中的积分必须以合适的方式离散化。 -
深度里兹方法
当使用偏微分方程的弱公式时,标准考虑变分公式,即考虑能量泛函,使得欧拉-拉格朗日方程是偏微分方程的弱公式。这个想法已经被 (Ritz, 1909) 利用来计算偏微分方程解的多项式逼近系数,并在 (E & Yu, 2018) 的神经网络背景下得到普及,他为这种方法创造了深度 Ritz 方法的名称。抽象地说,这种方法类似于残差公式。给定变分能量 E : X → R E : X \to R E:X→R,在希尔伯特空间 X X X 上,通过神经网络 u θ u_\theta uθ 参数化并得到损失函数 L ( θ ) : = E ( u θ ) L(\theta) := E(u_\theta) L(θ):=E(uθ) 。请注意,这种方法与 PINN 不同,例如对于泊松方程 − Δ u = f -\Delta u=f −Δu=f ,PINN剩余能量由 u ↦ ∥ Δ u + f ∥ L 2 ( Ω ) 2 u\mapsto\|\Delta u+f\|_{L^2(\Omega)}^2 u↦∥Δu+f∥L2(Ω)2 给出,其相应的变分能量则由下式给出 u ↦ 1 2 ∥ ∇ u ∥ L 2 ( Ω ) 2 − ∫ Ω f u d x u\mapsto\frac12\|\nabla u\|_{L^2(\Omega)}^2-\int_{\Omega}fu\mathrm{d}x u↦21∥∇u∥L2(Ω)2−∫Ωfudx。特别是,能量需要不同的函数平滑度,因此在不同的索博列夫空间上定义。在深度里兹方法中纳入基本边界值与 PINN 方法不同。在 PINN 中,对于任何 τ > 0 \tau > 0 τ>0,能量的唯一最小化器是 PDE 的解,而在深度里兹方法中,惩罚能量的最小化器解决了 Robin 边值问题,这可以解释为扰动问题。为了很好地近似原始问题,惩罚参数需要很大,这会导致病态问题(Muller & Zeinhofer,2022a;Courte & Zeinhofer,2023)
一般设置物理信息神经网络以及深度 Ritz 方法都属于最小化能量
E
:
X
→
R
E : X \to R
E:X→R 或更准确地说是参数空间上的相关目标函数
L
(
θ
)
:
=
E
(
u
θ
)
L(\theta) := E(u_\theta)
L(θ):=E(uθ) 的一般框架的神经网络。这里,假设
X
X
X 是函数的希尔伯特空间,并且由参数
θ
\theta
θ 的神经网络计算的函数
u
θ
u_\theta
uθ 位于
X
X
X 中,并假设
E
E
E 允许唯一的最小化器
u
⋆
∈
X
u^⋆ \in X
u⋆∈X。此外,假设参数化
P
:
R
p
→
X
,
θ
↦
u
θ
P : \mathbb R^p \to X, \theta \mapsto u_\theta
P:Rp→X,θ↦uθ 是可导的,其范围由
F
Θ
=
{
u
θ
:
θ
∈
R
p
}
\mathcal F_\Theta = \{u_\theta : \theta \in \mathbb R^p\}
FΘ={uθ:θ∈Rp} 表示。可以将该参数模型上的广义切线空间表示为:
T
θ
F
Θ
:
=
s
p
a
n
{
∂
θ
i
u
θ
:
i
=
1
,
…
,
p
}
T_\theta\mathcal{F}_\Theta:=\mathrm{span}\left\{\partial_{\theta_i}u_{\theta}:i=1,\ldots,p\right\}
TθFΘ:=span{∂θiuθ:i=1,…,p}
除了 PINN 训练过程的显着改进之外,原始 PINN 公式的基于梯度的优化迄今为止还无法打破一定的优化障碍,即使对于简单的情况也是如此。通常,在 L2 范数中测量的误差约为
1
0
−
3
10^{−3}
10−3。这种现象归因于 PINN 公式的刚度,如 (Wang et al., 2021) 中的实验验证。此外,偏微分方程的残差平方公式是条件数的平方——这在经典离散化方法中是众所周知的。由于离散偏微分方程会导致病态线性系统,这会恶化迭代求解器(例如标准梯度下降)的收敛性。另一方面,自然梯度下降通过保证更新方向遵循函数空间梯度信息来避免离散化的这种病态,其中偏微分方程问题通常具有更简单的结构。
方法提出
自然梯度的概念由 Amari 在监督学习和盲源分离的参数估计背景下推广(Amari,1998)。主要的想法是修改基于梯度的优化方案中的更新方向,以在参数的合适表示空间中模拟梯度。它通常归因于在表示空间上使用 Fisher 度量,而且 Fisher 度量的产品、Wasserstein 和 Sobolev 几何也已被成功使用。
自然梯度的一个微妙之处是函数空间中几何的定义。这可以通过公理化或通过潜在函数的 Hessian 矩阵来完成。本文遵循这样的想法,使用由凸函数空间目标的 Hessian 矩阵引起的自然梯度,其中自然梯度可以解释为广义高斯-牛顿方法,该方法已被建议用于监督学习任务的神经网络训练。与现有的工作相反,本文在应用中遇到无限维且非强凸的目标。
在这里,考虑在希尔伯特空间
X
X
X 上定义的凸能量
E
:
X
→
R
E : X \to \mathbb R
E:X→R 的最小化设置,它涵盖了物理信息神经网络和深度里兹法。作为优化网络参数的目标函数,像以前一样使用
L
(
θ
)
=
E
(
u
θ
)
L(\theta) = E(u_\theta)
L(θ)=E(uθ)。可以定义希尔伯特和能量格拉姆矩阵
G
H
(
θ
)
i
j
:
=
⟨
∂
θ
i
u
θ
,
∂
θ
j
u
θ
⟩
X
G_H(\theta)_{ij}:=\langle\partial_{\theta_i}u_{\theta},\partial_{\theta_j}u_{\theta}\rangle_X
GH(θ)ij:=⟨∂θiuθ,∂θjuθ⟩X
以及
G
E
(
θ
)
i
j
:
=
D
2
E
(
u
θ
)
(
∂
θ
i
u
θ
,
∂
θ
j
u
θ
)
G_E(\theta)_{\boldsymbol{i}j}:=D^2E(u_{\boldsymbol{\theta}})(\partial_{\boldsymbol{\theta}_i}u_{\boldsymbol{\theta}},\partial_{\boldsymbol{\theta}_j}u_{\boldsymbol{\theta}})
GE(θ)ij:=D2E(uθ)(∂θiuθ,∂θjuθ)
其中,
⟨
⋅
,
⋅
⟩
X
\langle\cdot,\cdot\rangle_X
⟨⋅,⋅⟩X 作为神经网络训练的 Sobolev 内积,Nurbekyan 等人提出了使用如下更新方式
∇
H
L
(
θ
)
=
G
H
(
θ
)
+
∇
L
(
θ
)
\nabla^HL(\theta)=G_H(\theta)^+\nabla L(\theta)
∇HL(θ)=GH(θ)+∇L(θ)
本文将其称为希尔伯特自然梯度(H-NG),在特殊情况下
X
X
X 是索博列夫空间索博列夫自然梯度。在有关自然梯度的文献中众所周知,
D
P
θ
∇
H
L
(
θ
)
=
Π
T
θ
F
Θ
(
∇
E
(
u
θ
)
)
.
DP_{\theta}\nabla^{H}L(\theta)=\Pi_{T_{\theta}\mathcal{F}_{\Theta}}(\nabla E(u_{\theta})).
DPθ∇HL(θ)=ΠTθFΘ(∇E(uθ)).
换句话说,遵循自然梯度相当于沿着希尔伯特空间梯度在函数空间中模型切线空间上的投影移动。通过 Hessian 矩阵识别函数空间梯度导致牛顿更新的观察激发了本文现在引入的能量自然梯度的概念。
定义1(能量自然梯度)。考虑问题
min
θ
∈
R
p
L
(
θ
)
\min_{\theta\in\mathbb{R}^p}L(\theta)
minθ∈RpL(θ),其中
L
(
θ
)
=
E
(
u
θ
)
L(\theta) = E(u_\theta)
L(θ)=E(uθ) 并用
∇
L
(
θ
)
\nabla L(\theta)
∇L(θ) 表示欧几里得梯度。然后可得到能量自然梯度(E-NG)定义如下:
∇
E
L
(
θ
)
:
=
G
E
+
(
θ
)
∇
L
(
θ
)
\nabla^EL(\theta):=G_E^+(\theta)\nabla L(\theta)
∇EL(θ):=GE+(θ)∇L(θ)
对于线性 PDE 算子
L
\mathcal L
L,残差产生二次能量,能量 Gram 矩阵采用以下形式:
G
E
(
θ
)
i
j
=
∫
Ω
L
(
∂
θ
i
u
θ
)
L
(
∂
θ
j
u
θ
)
d
x
+
τ
∫
∂
Ω
B
(
∂
θ
i
u
θ
)
B
(
∂
θ
j
u
θ
)
d
s
\begin{aligned} G_{E}(\theta)_{\boldsymbol{i}j}& =\int_{\Omega}\mathcal{L}(\partial_{\theta_i}u_\theta)\mathcal{L}(\partial_{\theta_j}u_\theta)\mathrm{d}x \\ &+\tau\int_{\partial\Omega}\mathcal{B}(\partial_{\theta_i}u_\theta)\mathcal{B}(\partial_{\theta_j}u_\theta)\mathrm{d}s \end{aligned}
GE(θ)ij=∫ΩL(∂θiuθ)L(∂θjuθ)dx+τ∫∂ΩB(∂θiuθ)B(∂θjuθ)ds
另一方面,深度里兹法的二次能量
E
(
u
)
=
1
2
a
(
u
,
u
)
−
f
(
u
)
E(u) = \frac 1 2 a(u, u) − f (u)
E(u)=21a(u,u)−f(u),其中 a 是对称且强制双线性形式,且
f
∈
X
∗
f \in X^*
f∈X∗ 可得出:
G
E
(
θ
)
i
j
=
a
(
∂
θ
i
u
θ
,
∂
θ
j
u
θ
)
G_E(\theta)_{ij}=a(\partial_{\theta_i}u_{\theta},\partial_{\theta_j}u_{\theta})
GE(θ)ij=a(∂θiuθ,∂θjuθ)
对于能量自然梯度,可以得到以下将能量自然梯度与牛顿更新相关的结果。
定理2(函数空间中的能量自然梯度)。如果假设
D
2
E
D^2E
D2E 在任何地方都是强制的,那么可以得到如下恒等式:
D
P
θ
∇
E
L
(
θ
)
=
Π
T
θ
F
Θ
D
2
E
(
u
θ
)
(
D
2
E
(
u
θ
)
−
1
∇
E
(
u
θ
)
)
DP_\theta\nabla^EL(\theta)=\Pi_{T_\theta\mathcal{F}_\Theta}^{D^2E(u_\theta)}(D^2E(u_\theta)^{-1}\nabla E(u_\theta))
DPθ∇EL(θ)=ΠTθFΘD2E(uθ)(D2E(uθ)−1∇E(uθ))
现在假设
E
E
E 是一个二次函数,具有有界且正定的二阶导数
D
2
E
=
a
D^2E = a
D2E=a,它允许极小值
u
∗
∈
X
u^* \in X
u∗∈X。那么可以得到:
D
P
θ
∇
E
L
(
θ
)
=
Π
T
θ
F
Θ
a
(
u
θ
−
u
∗
)
DP_{\theta}\nabla^{E}L(\theta)=\Pi_{T_{\theta}\mathcal{F}_{\Theta}}^{a}(u_{\theta}-u^{*})
DPθ∇EL(θ)=ΠTθFΘa(uθ−u∗)
证明想法,完整证明见论文附录。在
D
2
E
D^2E
D2E 具有强制性的情况下,它会在希尔伯特空间
X
X
X 上引入黎曼度量。由于相对于该度量的梯度由
D
2
E
(
u
)
−
1
∇
E
(
u
)
D^2E(u)^{−1} \nabla E(u)
D2E(u)−1∇E(u) 给出,因此前式类似于有限维情况或希尔伯特空间自然梯度的情况。在能量
E
E
E 是二次且
D
2
E
=
a
D^2E = a
D2E=a 有界且非简并的情况下,相对于内积
a
a
a 的梯度没有经典定义。然而,可以检查
a
(
u
−
u
∗
,
v
)
=
D
E
(
u
)
v
a(u − u^*, v) = DE(u)v
a(u−u∗,v)=DE(u)v,即误差
u
−
u
∗
u − u^*
u−u∗ 可以解释为相对于内积
a
a
a 的梯度,从而产生后式。
特别是,可以从上式中看到,在参数空间中使用能量自然梯度与函数空间中的牛顿更新密切相关,其中对于二次能量,牛顿方向由误差 u θ − u ⋆ u_\theta − u^⋆ uθ−u⋆ 给出。
H-NG 和 E-NG 的复杂性
H-NG 和 E-NG 的计算(取决于 Gram 矩阵 G H G_H GH 和 G E G_E GE 的组装)同样昂贵。幸运的是,格拉姆矩阵的计算成本通常同样昂贵。对于二次问题,Hessian 矩阵通常并不比希尔伯特内积更难评估,甚至对于非二次情况,内积方面的封闭形式表达式也经常可用。请注意,HNG 模拟 GD,E-NG 模拟 X 中的牛顿法。实际上,自然梯度的计算成本很高,因为它需要求解线性方程组,其复杂度为 O ( p 3 ) O(p^3) O(p3),其中 p p p 为参数维度。将其与计算梯度的 O ( p ) O(p) O(p) 成本进行比较。在本文的实验中,发现 E-NGD 与 GD 和 Adam 相比,即使后者被允许更多的计算时间,也能实现显着更高的精度。
实验结果
本文在四个问题上测试了能量自然梯度方法:二维泊松方程的 PINN 形式、五维泊松方程的 PINN 形式、一维热方程的 PINN 形式以及一维非线性椭圆方程的深度里兹形式。
方法描述
对于本文所有的数值实验,通过线性搜索实现能量自然梯度步骤,如下算法所述。
本文选择线性搜索的区间 [ 0 , 1 ] [0, 1] [0,1] 来确定学习率,因为学习率为 1 将对应于函数空间中的近似牛顿步。但由于模型的参数化是非线性的,有利于进行线搜索,不能简单地选择牛顿步长。在本文的实验中,在 [ 0 , 1 ] [0, 1] [0,1] 上的对数间隔网格上使用网格搜索来确定学习率 η ∗ \eta ^* η∗。虽然很简陋,但它可以很容易地并行化,并且在本文的实验中快速高效地执行。
Gram 矩阵
G
E
G_E
GE 的组装可以高效地并行完成,从而避免了索引对
(
i
,
j
)
(i, j)
(i,j) 上可能代价高昂的循环。本文不是计算 Gram 矩阵
G
E
(
θ
)
G_E(\theta)
GE(θ) 的伪逆,而是解决如下最小二乘问题:
∇
E
L
(
θ
)
∈
arg
min
ψ
∈
R
p
∥
G
E
(
θ
)
ψ
−
∇
L
(
θ
)
∥
2
2
\nabla^EL(\theta)\in\arg\min_{\psi\in\mathbb{R}^p}\|G_E(\theta)\psi-\nabla L(\theta)\|_2^2
∇EL(θ)∈argψ∈Rpmin∥GE(θ)ψ−∇L(θ)∥22
在具体实现中,最小二乘求解的 JAX 实现依赖于奇异值分解。
为了对损失函数以及 Gram 矩阵的条目中出现的积分进行数值评估,本文对规则网格上的固定积分点和重复随机绘制的积分点进行了实验。本文根据标准差为 0.1 且均值消失的高斯函数初始化网络的权重和偏差。
本文报告优化过程期间和之后的相对 L2 和 H1 误差。为此,本文使用了比优化期间多 10 倍的积分点。本文将能量 NG 的效率与以下优化器进行比较。首先,考虑在对数网格上进行线性搜索的普通梯度下降(在实验中表示为 GD)。然后,使用指数递减的学习率计划来测试 Adam 的性能,为防止振荡,从 1 0 − 3 10^{−3} 10−3 的初始学习率开始,在 15000 步之后开始每 10000 步减少 1 0 − 1 10^{−1} 10−1 倍直到达到最小学习率 1 0 − 7 10^{−7} 10−7 或完成最大迭代量。本文还与拟牛顿法 BFGS (Nocedal & Wright, 1999) 进行了比较。最后,用线搜索测试希尔伯特自然梯度下降(记为H-NGD)
泊松问题
考虑如下泊松问题:
−
Δ
u
(
x
,
y
)
=
f
(
x
,
y
)
=
2
π
2
sin
(
π
x
)
sin
(
π
y
)
-\Delta u(x,y)=f(x,y)=2\pi^2\sin(\pi x)\sin(\pi y)
−Δu(x,y)=f(x,y)=2π2sin(πx)sin(πy)
在边界值为零的单位正方形
[
0
,
1
]
2
[0, 1]^2
[0,1]2 上。解由下式给出:
u
∗
(
x
,
y
)
=
sin
(
π
x
)
sin
(
π
y
)
u^*(x,y)=\sin(\pi x)\sin(\pi y)
u∗(x,y)=sin(πx)sin(πy)
问题的 PINN 损失如下:
L
(
θ
)
=
1
N
Ω
∑
i
=
1
N
Ω
(
Δ
u
θ
(
x
i
,
y
i
)
+
f
(
x
i
,
y
i
)
)
2
+
1
N
∂
Ω
∑
i
=
1
N
∂
Ω
u
θ
(
x
i
b
,
y
i
b
)
2
\begin{aligned} L(\theta)=\frac1{N_{\Omega}}\sum_{i=1}^{N_{\Omega}}(\Delta u_{\theta}(x_{i},y_{i})+f(x_{i},y_{i}))^2\\ +\frac1{N_{\partial\Omega}}\sum_{i=1}^{N_{\partial\Omega}}u_{\theta}(x_{i}^{b},y_{i}^{b})^2 \end{aligned}
L(θ)=NΩ1i=1∑NΩ(Δuθ(xi,yi)+f(xi,yi))2+N∂Ω1i=1∑N∂Ωuθ(xib,yib)2
其中
{
(
x
i
,
y
i
)
}
i
=
1
,
.
.
.
,
N
Ω
\{(x_i, y_i)\}_{i=1,...,N_\Omega}
{(xi,yi)}i=1,...,NΩ 表示内部配置点,
{
(
x
i
b
,
y
i
b
)
}
i
=
1
,
.
.
.
,
N
∂
Ω
\{(x^b_i, y^b_i )\}_{i=1,...,N_{\partial \Omega}}
{(xib,yib)}i=1,...,N∂Ω 表示
∂
Ω
\partial \Omega
∂Ω 上的配置点。在这种情况下,
H
2
(
Ω
)
H^2(\Omega)
H2(Ω) 的能量内积由下式给出:
a
(
u
,
v
)
=
∫
Ω
Δ
u
Δ
v
d
x
+
∫
∂
Ω
u
v
d
s
a(u,v)=\int_{\Omega}\Delta u\Delta v\mathrm{d}x+\int_{\partial\Omega}uv\mathrm{d}s
a(u,v)=∫ΩΔuΔvdx+∫∂Ωuvds
请注意,此内积对
H
2
(
Ω
)
H^2(\Omega)
H2(Ω) 没有强制作用,并且与
H
2
(
Ω
)
H^2(\Omega)
H2(Ω) 内积不同。上式中的积分是使用与PINN 损失函数
L
L
L 的定义相同的配置点来计算的。为了逼近解
u
∗
u^*
u∗,使用浅层神经网络,以双曲正切作为激活函数,宽度为 64,因此有 257 个可训练权重。本文在
Ω
\Omega
Ω 内部选择 900 个等距配置点,在边界上选择 120 个等距配置点。能量自然梯度下降和希尔伯特自然梯度下降分别应用 500 次迭代,而对 GD 和 Adam 进行 200000 次迭代训练.
如上表和上图所示,可以观察到 E-NG 更新需要相对较少的迭代来产生泊松方程的高精度近似解。请注意,H-NG 下降并未收敛,这正强调了使用函数空间目标的 Hessian 几何信息的重要性,就像 E-NG 下降中所做的那样。
考虑的一阶优化器,即 Adam 和 vanilla 梯度下降可靠地降低了相对误差,但即使允许更多的迭代次数,也无法达到高于 6.9 ⋅ 1 0 − 4 6.9 · 10^{−4} 6.9⋅10−4 的精度。拟牛顿法 BFGS(Nocedal & Wright,1999)比一阶方法获得了更高的精度,但与表 1 相比,能量自然梯度法的精度仍然大约高两个数量级。
根据当前的实现和考虑的网络规模,一次自然梯度更新的成本仅为 Adam 算法一次迭代的两倍到三倍(见下表)。
使用 Adam 等优化器训练 PINN 模型轻松需要 100 倍自然梯度训练所需的迭代次数(无法产生高度准确的解决方案)使所提出的方法更快、更准确。请注意,使用自然梯度方法的优化时间不到一分钟,而使用 Adam 优化器的优化时间需要两个小时以上,这是两个数量级的改进。
为了说明能量自然梯度 ∇ E L ( θ ) \nabla^EL(\theta) ∇EL(θ)、标准梯度 ∇ L ( θ ) \nabla L(\theta) ∇L(θ) 和误差 u θ − u ∗ u_\theta − u^* uθ−u∗ 之间的差异,下图绘制了在初始化时函数空间中的有效更新方向 D P ( θ ) ∇ L ( θ ) DP(\theta)\nabla L(\theta) DP(θ)∇L(θ) 和 D P ( θ ) ∇ E L ( θ ) DP(\theta)\nabla^EL(\theta) DP(θ)∇EL(θ) 。
显然,能量自然梯度更新方向与误差的匹配程度比普通参数梯度好得多。
高维度的例子
作为更高维度的例子,再次考虑如下五个空间维度的泊松方程:
−
Δ
u
=
f
i
n
[
0
,
1
]
5
,
u
(
x
)
=
∑
k
=
1
5
sin
(
π
x
k
)
o
n
∂
[
0
,
1
]
5
.
\begin{aligned} &-\Delta u=f&&\mathrm{in~}[0,1]^5,\\ &u(x)=\sum_{k=1}^5\sin(\pi x_k)&&\mathrm{on~}\partial[0,1]^5. \end{aligned}
−Δu=fu(x)=k=1∑5sin(πxk)in [0,1]5,on ∂[0,1]5.
假设其真实解为:
u
∗
:
R
5
→
R
,
x
↦
∑
k
=
1
5
sin
(
π
x
k
)
u^*:\mathbb{R}^5\to\mathbb{R},\quad x\mapsto\sum_{k=1}^5\sin(\pi x_k)
u∗:R5→R,x↦k=1∑5sin(πxk)
因此
f
=
π
2
u
∗
f = \pi ^2u^*
f=π2u∗。对于给定的一组内部和边界搭配点
(
x
1
i
,
.
.
.
,
x
5
i
)
i
=
1
,
.
.
.
,
N
Ω
⊆
[
0
,
1
]
5
(x^i_1, . . . , x^i_5)_{i=1,...,N_\Omega} ⊆ [0, 1]^5
(x1i,...,x5i)i=1,...,NΩ⊆[0,1]5 和
(
x
1
b
,
i
,
.
.
.
,
x
5
b
,
i
)
i
=
1
,
.
.
.
,
N
∂
Ω
⊆
∂
[
0
,
1
]
5
(x^{b,i} _1 , . . . , x^{b,i} _5 )_{i=1,...,N_{\partial \Omega}} ⊆ \partial[0, 1]^5
(x1b,i,...,x5b,i)i=1,...,N∂Ω⊆∂[0,1]5 定义损失函数和能量内积定义同上。在此示例中,通过在训练过程的每次迭代中重复采集
N
Ω
=
3000
N_\Omega = 3000
NΩ=3000 个随机内部和
N
∂
Ω
=
500
N_{\partial \Omega} = 500
N∂Ω=500 个边界搭配点。使用具有输入维度 5 和 64 个隐藏神经元以及双曲正切激活的浅层神经网络。
如上图所示,能量自然梯度法在相对较短的时间内产生高精度的解。 Adam 和普通梯度下降的收敛行为再次显示出快速饱和行为,并且都无法产生具有竞争力的准确解决方案。在此示例中,拟牛顿法 BFGS 无法产生比 Adam 优化器更准确的解。同样,E-NGD 不仅是最准确的优化器,精度提高了两个数量级,而且在实现这一精度的同时,时间也减少了一个数量级。
热方程
考虑如下一维热方程:
∂
t
u
(
t
,
x
)
=
1
4
∂
x
2
u
(
t
,
x
)
f
o
r
(
t
,
x
)
∈
[
0
,
1
]
2
u
(
0
,
x
)
=
sin
(
π
x
)
f
o
r
x
∈
[
0
,
1
]
u
(
t
,
x
)
=
0
f
o
r
(
t
,
x
)
∈
[
0
,
1
]
×
{
0
,
1
}
.
\begin{aligned} \partial_tu(t,x)&=\frac14\partial_x^2u(t,x)\quad\mathrm{for~}(t,x)\in[0,1]^2\\ u(0,x)&=\sin(\pi x)\quad\mathrm{for~}x\in[0,1]\\ u(t,x)&=0\quad\mathrm{for~}(t,x)\in[0,1]\times\{0,1\}. \end{aligned}
∂tu(t,x)u(0,x)u(t,x)=41∂x2u(t,x)for (t,x)∈[0,1]2=sin(πx)for x∈[0,1]=0for (t,x)∈[0,1]×{0,1}.
解由下式给出:
u
∗
(
t
,
x
)
=
exp
(
−
π
2
t
4
)
sin
(
π
x
)
u^*(t,x)=\exp\left(-\frac{\pi^2t}4\right)\sin(\pi x)
u∗(t,x)=exp(−4π2t)sin(πx)
PINN 损失为:
L
(
θ
)
=
1
N
Ω
T
∑
i
=
1
N
Ω
T
(
∂
t
u
θ
(
t
i
,
x
i
)
−
1
4
∂
x
2
u
θ
(
t
i
,
x
i
)
)
2
+
1
N
i
n
∑
i
=
1
N
Ω
(
u
θ
(
0
,
x
i
i
n
)
−
sin
(
π
x
i
i
n
)
)
2
+
1
N
∂
Ω
∑
i
=
1
N
∂
Ω
u
θ
(
t
i
b
,
x
i
b
)
2
\begin{aligned} L(\theta)& =\frac1{N_{\Omega_T}}\sum_{\boldsymbol{i}=1}^{\boldsymbol{N_{\Omega_T}}}\left(\partial_{\boldsymbol{t}}u_{\boldsymbol{\theta}}(t_i,x_i)-\frac14\partial_{\boldsymbol{x}}^2u_{\boldsymbol{\theta}}(t_i,x_i)\right)^2 \\ &+\frac1{N_{\mathrm{in}}}\sum_{i=1}^{N_{\Omega}}\left(u_{\theta}(0,x_{i}^{\mathrm{in}})-\sin(\pi x_{i}^{\mathrm{in}})\right)^2 \\ &+\frac1{N_{\partial\Omega}}\sum_{i=1}^{N_{\partial\Omega}}u_\theta(t_i^b,x_i^b)^2 \end{aligned}
L(θ)=NΩT1i=1∑NΩT(∂tuθ(ti,xi)−41∂x2uθ(ti,xi))2+Nin1i=1∑NΩ(uθ(0,xiin)−sin(πxiin))2+N∂Ω1i=1∑N∂Ωuθ(tib,xib)2
其中
{
(
t
i
,
x
i
)
}
i
=
1
,
.
.
.
,
N
Ω
T
\{(t_i, x_i)\}_{i=1,...,N_{\Omega T}}
{(ti,xi)}i=1,...,NΩT 表示时空圆柱内部的配置点,
{
(
t
i
b
,
x
i
b
)
}
i
=
1
,
.
.
.
,
N
∂
Ω
\{(t^b _i , x^b _i )\}_{i=1,...,N_{\partial \Omega}}
{(tib,xib)}i=1,...,N∂Ω 表示空间边界上的配置点,
{
(
x
i
i
n
)
}
i
=
1
,
.
.
.
,
N
i
n
\{(x^{in} _i )\}_{i=1,...,N_{in}}
{(xiin)}i=1,...,Nin 表示初始条件的配置点。能量内积定义在空间上:
a
:
(
H
1
(
I
,
L
2
(
Ω
)
)
∩
L
2
(
I
,
H
2
(
Ω
)
)
)
2
→
R
a\colon\left(H^1(I,L^2(\Omega))\cap L^2(I,H^2(\Omega))\right)^2\to\mathbb{R}
a:(H1(I,L2(Ω))∩L2(I,H2(Ω)))2→R
即:
a
(
u
,
v
)
=
∫
0
1
∫
Ω
(
∂
t
u
−
1
4
∂
x
2
u
)
(
∂
t
v
−
1
4
∂
x
2
v
)
d
x
d
t
+
∫
Ω
u
(
0
,
x
)
v
(
0
,
x
)
d
x
+
∫
I
×
∂
Ω
u
v
d
s
d
t
.
\begin{gathered} a(u,v) =\int_0^1\int_\Omega\left(\partial_tu-\frac14\partial_x^2u\right)\left(\partial_tv-\frac14\partial_x^2v\right)\mathrm{d}x\mathrm{d}t \\ +\int_{\Omega}u(0,x)v(0,x)\mathrm{d}x+\int_{I\times\partial\Omega}uv\mathrm{d}s\mathrm{d}t. \end{gathered}
a(u,v)=∫01∫Ω(∂tu−41∂x2u)(∂tv−41∂x2v)dxdt+∫Ωu(0,x)v(0,x)dx+∫I×∂Ωuvdsdt.
在本文的实现中,内积由与损失函数定义中相同的正交点离散化。网络架构和训练过程与前面的泊松问题示例相同,运行两种 NG 方法进行 2000 次迭代。同样在这个例子中,能量自然梯度方法显示了其高精度和高效率。下表显示了训练后的相对 L2 误差
上图为训练过程的可视化
下表为运行时间
上图为不同梯度的视觉比较
深度里兹法的非线性示例
本文利用深度里兹法测试了非线性问题的能量自然梯度法。考虑如下寻找能量最小值的一维变分问题:
E
(
u
)
:
=
1
2
∫
Ω
∣
u
′
∣
2
d
x
+
1
4
∫
Ω
u
4
d
x
−
∫
Ω
f
u
d
x
E(u):=\frac12\int_{\Omega}|u^{\prime}|^2\mathrm{d}x+\frac14\int_{\Omega}u^4\mathrm{d}x-\int_{\Omega}fu\mathrm{~d}x
E(u):=21∫Ω∣u′∣2dx+41∫Ωu4dx−∫Ωfu dx
其中
Ω
=
[
−
1
,
1
]
\Omega = [−1, 1]
Ω=[−1,1],
f
(
x
)
=
π
2
cos
(
π
x
)
+
cos
3
(
π
x
)
f (x) = \pi^2 \cos(\pi x) + \cos^3(\pi x)
f(x)=π2cos(πx)+cos3(πx)。相关的欧拉拉格朗日方程产生如下非线性偏微分方程:
−
u
′
′
+
u
3
=
f
i
n
Ω
∂
n
u
=
0
o
n
∂
Ω
\begin{aligned} -u^{\prime\prime}+u^3&=f&\mathrm{~in~}\Omega\\ \partial_nu&=0&\mathrm{~on~}\partial\Omega \end{aligned}
−u′′+u3∂nu=f=0 in Ω on ∂Ω
因此最小化器由
u
∗
(
x
)
=
cos
(
π
x
)
u^*(x) = \cos(\pi x)
u∗(x)=cos(πx) 给出。由于能量不是二次方,因此能量内积取决于
u
∈
H
1
(
Ω
)
u \in H^1(\Omega)
u∈H1(Ω),由下式给出:
D
2
E
(
u
)
(
v
,
w
)
=
∫
Ω
v
′
w
′
d
x
+
3
∫
Ω
u
2
v
w
d
x
D^2E(u)(v,w)=\int_\Omega v^{\prime}w^{\prime}\operatorname{d}x+3\int_\Omega u^2vw\operatorname{d}x
D2E(u)(v,w)=∫Ωv′w′dx+3∫Ωu2vwdx
为了离散能量和内积,本文使用具有 20000 个等距正交点的梯形积分。使用宽度为 32 个神经元的浅层神经网络和双曲正切作为激活函数。
可以再次观察到,E-NG 更新有效地导致了解决方案的非常准确的近似,
上图为训练过程可视化,下表为获得的相对 L2 误差。
下表报告了各个方法的计算时间。
在这个例子中,Hilbert NG 下降同样有效。注意能量内积和希尔伯特空间内积在这种情况下非常相似。同样,Adam 和标准梯度下降较早饱和,其误差比自然梯度方法高得多。可以观察到,使用深度里兹法获得高精度需要高精度的积分,这就是使用相对精细的网格和梯形积分的原因。
总结
本文从牛顿下降法和自然梯度法获得灵感,提出使用能量自然梯度方法来对神经网络参数进行优化,并在多个样例上取得了显著的精度改善。
相关链接: