论文阅读:Is $L^2$ Physics-Informed Loss Always Suitable for Training Physics-Informed Neural Network
Is L 2 L^2 L2 Physics-Informed Loss Always Suitable for Training Physics-Informed Neural Network
问题分析
稳定性
在神经网络的训练过程中,当且仅当损失项 ℓ Ω , p ( u ) \ell_{\Omega,p}(u) ℓΩ,p(u) 和 ℓ ∂ Ω , p ( u ) \ell_{\partial\Omega,p}(u) ℓ∂Ω,p(u)均为零时, u ( x ) u(x) u(x) 才是 PDE 的精确解。然而,在实践中,由于优化过程的随机性或神经网络的表达能力,通常只能获得很小但非零的损失值。在这种情况下,自然会出现一个问题:具有较小损失的学习结果 u ( x ) u(x) u(x) 是否对应于精确解 u ∗ ( x ) u*(x) u∗(x) 的良好逼近?这种性质与 PDE 中的稳定性高度相关。
假设 Z 1 、 Z 2 、 Z 3 Z_1、Z_2、Z_3 Z1、Z2、Z3 是三个 Banach 空间。如果当 ∥ L u ˉ ( x ) − φ ˙ ( x ) ∥ Z 1 , ∥ B u ( x ) − ψ ˉ ( x ) ∥ Z 2 → 0. \lVert\bar{\mathcal{L}u}(x)-\dot{\varphi}(x)\rVert_{Z_1},\lVert\mathcal{B}u(x)-\bar{\psi}(x)\rVert_{Z_2}\to0. ∥Luˉ(x)−φ˙(x)∥Z1,∥Bu(x)−ψˉ(x)∥Z2→0. 时,对任意 u u u 存在 ∥ u ∗ ( x ) − u ( x ) ∥ Z 3 = O ( ∥ L u ( x ) − φ ( x ) ∥ Z 1 + ∥ B u ( x ) − ψ ( x ) ∥ Z 2 ) \|u^*(x)-u(x)\|_{Z_3}=O(\|\mathcal{L}u(x)-\varphi(x)\|_{Z_1}+\|\mathcal{B}u(x)-\psi(x)\|_{Z_2}) ∥u∗(x)−u(x)∥Z3=O(∥Lu(x)−φ(x)∥Z1+∥Bu(x)−ψ(x)∥Z2) ,则认为PDE 定义为是 ( Z 1 , Z 2 , Z 3 ) (Z_1, Z_2, Z_3) (Z1,Z2,Z3)-稳定的。
通过上述定义可以看出,如果一个PDE是 ( L 2 ( Ω ) , L 2 ( ∂ Ω ) , Z ) (L^2(\Omega),L^2(\partial\Omega),Z) (L2(Ω),L2(∂Ω),Z) -稳定的,那么就可以使用 L 2 L^2 L2 范数来最小化物理信息损失 ; ∥ L u ( x ) − φ ( x ) ∥ L 2 ( Ω ) 2 ;\|\mathcal{L}u(x)-\varphi(x)\|_{L^{2}(\Omega)}^{2} ;∥Lu(x)−φ(x)∥L2(Ω)2 和 ∥ B ˚ u ( x ) − ψ ( x ) ∥ L 2 ( ∂ Ω ) 2 \|\mathring{\mathcal{B}}u(x)-\psi(x)\|_{L^{2}(\partial\Omega)}^{2} ∥B˚u(x)−ψ(x)∥L2(∂Ω)2 并且,当损失趋近于零时,结果必然是精确解。
但是,这种稳定性并不广泛存在于PDE中,有很多不稳定的方程,比如:逆热方程。同时,即使方程是稳定的,也不一定是 ( L 2 ( Ω ) , L 2 ( ∂ Ω ) , Z ) (L^2(\Omega),L^2(\partial\Omega),Z) (L2(Ω),L2(∂Ω),Z) -稳定的。比如一些实际的高维HJB方程是稳定的,但不是 ( L 2 ( Ω ) , L 2 ( ∂ Ω ) , Z ) (L^2(\Omega),L^2(\partial\Omega),Z) (L2(Ω),L2(∂Ω),Z) -稳定的,使用 L 2 L^2 L2 范数将无法在实践中找到近似解。
以一类HJB方程为例,其成本率函数为:
r
(
x
,
m
)
=
a
1
∣
m
1
∣
α
1
+
⋯
+
a
n
∣
m
n
∣
α
n
−
φ
(
x
,
t
)
.
r(x,m)=a_1|m_1|^{\alpha_1}+\cdots+a_n|m_n|^{\alpha_n}-\varphi(x,t).
r(x,m)=a1∣m1∣α1+⋯+an∣mn∣αn−φ(x,t). 对应的 Hamilton-Jacobi-Bellman 方程可以改写为:
{
L
H
J
B
u
:
=
∂
t
u
(
x
,
t
)
+
1
2
σ
2
Δ
u
(
x
,
t
)
−
∑
i
=
1
n
A
i
∣
∂
x
i
u
∣
c
i
=
φ
(
x
,
t
)
(
x
,
t
)
∈
R
n
×
[
0
,
T
]
B
H
J
B
u
:
=
u
(
x
,
T
)
=
g
(
x
)
x
∈
R
n
,
\begin{cases}\mathcal{L}_{HJB}u:=\partial_t u(x,t)+\frac{1}{2}\sigma^2\Delta u(x,t)-\sum\limits_{i=1}^n A_i|\partial_{x_i}u|^{c_i}=\varphi(x,t)&(x,t)\in\mathbb{R}^n\times[0,T]\\ \mathcal{B}_{HJB}u:=u(x,T)=g(x)\quad&x\in\mathbb{R}^n\end{cases},
⎩
⎨
⎧LHJBu:=∂tu(x,t)+21σ2Δu(x,t)−i=1∑nAi∣∂xiu∣ci=φ(x,t)BHJBu:=u(x,T)=g(x)(x,t)∈Rn×[0,T]x∈Rn,
其中, A i = ( a i α i ) − 1 α i − 1 − a i ( a i α i ) − α i α i − 1 ∈ ( 0 , + ∞ ) A_i=(a_i\alpha_i)^{-\frac{1}{\alpha_i-1}}-a_i(a_i\alpha_i)^{-\frac{\alpha_i}{\alpha_i-1}}\in(0,+\infty) Ai=(aiαi)−αi−11−ai(aiαi)−αi−1αi∈(0,+∞) , c i = α i α i − 1 ∈ ( 1 , ∞ ) c_i=\frac{\alpha_i}{\alpha_i-1}\in(1,\infty) ci=αi−1αi∈(1,∞) . 对于函数 f : X → R f:X\to\mathbb{R} f:X→R ,其中 X X X 是一个可测空间,用 supp f \text{supp} f suppf 表示 f f f 的支集,即 { x ∈ X : f ( x ) ≠ 0 } \{x\in X:f(x)\neq0\} {x∈X:f(x)=0} 的闭包。
Sobolev 空间
分析 PDE 的一个重要概念是 Sobolev 空间,其定义如下:
对于
m
∈
N
,
p
∈
[
1
,
+
∞
)
m\in\mathbb{N},p\in[1,+\infty)
m∈N,p∈[1,+∞) 和一个开集
Ω
⊂
R
n
\Omega \subset \mathbb{R}^n
Ω⊂Rn , Sobolev 空间
W
m
,
p
(
Ω
)
W^{m,p}(\Omega)
Wm,p(Ω) 定义为
{
f
(
x
)
∈
L
p
(
Ω
)
:
D
α
f
∈
L
p
(
Ω
)
,
∀
α
∈
N
n
,
∣
α
∣
≤
m
}
\{f(x)\in L^p(\Omega):D^\alpha f\in L^p(\Omega),\forall\alpha\in\mathbb{N}^n,|\alpha|\le m\}
{f(x)∈Lp(Ω):Dαf∈Lp(Ω),∀α∈Nn,∣α∣≤m} . 函数空间
W
m
,
p
(
Ω
)
W^{m,p}(\Omega)
Wm,p(Ω) 具备 Sobolev 范数:
∥
f
∥
W
m
,
p
(
Ω
)
=
(
∑
∣
α
∣
≤
m
∥
D
α
f
∥
L
p
(
Ω
)
p
)
1
p
.
\begin{aligned} & \\ &\|f\|_{W^{m,p}(\Omega)}& =\left(\sum\limits_{|\alpha|\leq m}\|D^\alpha f\|_{L^p(\Omega)}^p\right)^\frac{1}{p}. \end{aligned}
∥f∥Wm,p(Ω)=
∣α∣≤m∑∥Dαf∥Lp(Ω)p
p1.
上述定义同样适用于时空域上的函数
Q
⊆
R
n
×
[
0
,
T
]
Q \subseteq \mathbb{R}^n \times [0,T]
Q⊆Rn×[0,T] ,稍微更改符号即可得到
:
W
m
,
p
(
Q
)
=
.
{
f
(
x
,
t
)
∈
L
p
(
Q
)
:
.
D
α
f
∈
L
p
(
Q
)
,
∀
α
∈
N
n
,
∣
α
∣
≤
m
}
:W^{m,p}(Q)\stackrel{.}{=}\{f(x,t)\in L^{p}(Q)\stackrel{.}{:}{{D}}^{\alpha}f\in L^p(Q),\forall\alpha\in\mathbb{N}^n,|\alpha|\leq m\}
:Wm,p(Q)=.{f(x,t)∈Lp(Q):.Dαf∈Lp(Q),∀α∈Nn,∣α∣≤m}
其中, D α D^\alpha Dα 仅对空间变量 x x x 进行运算。范数 ∥ ⋅ ∥ W m , p ( Q ) \|\cdot\|_{W^{m,p}(Q)} ∥⋅∥Wm,p(Q) 也可以类似进行定义。
HJB方程的稳定性
值函数 u(x, t) 是以下偏微分方程的唯一解,称为 Hamilton-Jacobi-Bellman 方程:
{
∂
t
u
(
x
,
t
)
+
1
2
σ
2
Δ
u
(
x
,
t
)
+
min
m
∈
M
[
r
(
x
,
m
(
t
,
x
)
,
t
)
+
∇
u
⋅
m
t
]
=
0
u
(
x
,
T
)
=
g
(
x
)
.
\begin{cases}\partial_t u(x,t)+\frac{1}{2}\sigma^2\Delta u(x,t)+\min\limits_{m\in\mathcal{M}}[r(x,m(t,x),t)+\nabla u\cdot m_t]=0\\ u(x,T)=g(x).\end{cases}
{∂tu(x,t)+21σ2Δu(x,t)+m∈Mmin[r(x,m(t,x),t)+∇u⋅mt]=0u(x,T)=g(x).
其中
M
\mathcal{M}
M 表示一组可能的控制函数。
对于
p
,
q
≥
1
p,q \ge 1
p,q≥1 令
r
0
=
(
n
+
2
)
q
n
+
q
r_0= \frac {(n+2)q} {n+q}
r0=n+q(n+2)q ,假设下述不等式对
p
,
q
p,q
p,q 以及
r
0
r_0
r0 成立:
p
≥
max
{
2
,
(
1
−
1
c
ˉ
)
n
}
;
q
>
(
c
ˉ
−
1
)
n
2
(
2
−
c
ˉ
)
n
+
2
;
1
r
0
≥
1
p
−
1
n
p\ge\max\left\{2,\left(1-\frac{1}{\bar{c}}\right)n\right\};q>\frac{(\bar{c}-1)n^2}{(2-\bar{c})n+2};\frac{1}{r_0}\ge\frac{1}{p}-\frac{1}{n}
p≥max{2,(1−cˉ1)n};q>(2−cˉ)n+2(cˉ−1)n2;r01≥p1−n1
其中 $\bar c= \max\limits_{1\leq i\leq n}c_i $ ,于是,对于任意的
r
∈
[
1
,
r
0
)
r \in [1,r_0)
r∈[1,r0) 以及任何有界开集
Q
⊆
R
n
×
[
0
,
T
]
Q \subseteq \mathbb{R}^n \times [0,T]
Q⊆Rn×[0,T] ,上述公式就是 $(Lp(\mathbb{R}n\times[0,T]),Lq(\mathbb{R}n),W^{1,r}(Q)) $ -稳定的,对于
c
ˉ
≤
2.
\bar{c}\leq2.
cˉ≤2.
上述证明指出,当 p , q = Ω ( n ) p,q = \Omega (n) p,q=Ω(n) 时,上文提到的HJB方程是 ( L p , L q , W 1 , r ) (L^p,L^q,W^{1,r}) (Lp,Lq,W1,r)-稳定的
因此,当状态函数 n 的维数很大时,如果 p 和 q 很小,则上文中提到的 HJB 方程是 ( L p , L q , W 1 , r ) (L^p,L^q,W^{1,r}) (Lp,Lq,W1,r)-不稳定的。此外,由于根据定义 L r = W 0 , r L^{r}{=}W^{0,r} Lr=W0,r ,这也意味着上文中提到的 HJB 方程甚至不是 ( L p , L q , L r ) (L^p,L^q,L^r) (Lp,Lq,Lr)-稳定的。因此,对于高维 HJB 问题,如果使用经典的 L 2 L^2 L2 物理信息损失来训练 PINN,即使损失非常小,学习到的解也可能与 u ∗ u^* u∗ 相距任意远。
解决方法
上述结果表明,我们应该在损失项
ℓ
Ω
,
p
(
u
)
\ell_{\Omega,p}(u)
ℓΩ,p(u) 和
ℓ
∂
Ω
,
p
(
u
)
\ell_{\partial\Omega,p}(u)
ℓ∂Ω,p(u) 中使用较大的
p
p
p 和
q
q
q 值,以保证对于高维 HJB 的近似解
u
u
u 足够接近
u
∗
u^∗
u∗ 。由于当
p
p
p 很大时,
L
p
L^p
Lp范数 和
L
∞
L^\infty
L∞范数 的行为相似。因此可以用
L
∞
L^\infty
L∞范数代替
L
p
L^p
Lp范数,直接优化
ℓ
Ω
,
∞
(
u
)
\ell_{\Omega,\infty}(u)
ℓΩ,∞(u) 和
ℓ
∂
Ω
,
∞
(
u
)
\ell_{\partial\Omega,\infty}(u)
ℓ∂Ω,∞(u) 。总的来说,训练目标可以表述为:
min
u
ℓ
∞
(
u
)
=
sup
x
∈
Ω
∣
L
u
(
x
)
−
φ
(
x
)
∣
+
λ
sup
x
∈
∂
Ω
∣
B
u
(
x
)
−
ψ
(
x
)
∣
\min\limits_{u}\ell_{\infty}(u)=\sup\limits_{x\in\Omega}|\mathcal{L}u(x)-\varphi(x)|+\lambda\sup\limits_{x\in\partial\Omega}|\mathcal{B}u(x)-\psi(x)|
uminℓ∞(u)=x∈Ωsup∣Lu(x)−φ(x)∣+λx∈∂Ωsup∣Bu(x)−ψ(x)∣
上式可以看作是一个最小-最大优化问题。内环是一个最大化问题,用于在
Ω
\Omega
Ω 和
∂
Ω
\partial\Omega
∂Ω 上找到
u
u
u 最违反 PDE 的数据点,外环是一个最小化问题,用于找到
u
u
u(即神经网络参数),使这些点上的损失最小化。
在深度学习中,这种最小-最大优化问题得到了深入研究,而对抗训练是许多应用中最有效的学习方法之一。本文利用对抗性训练实现。在每个训练步骤中,模型参数和数据点都会迭代更新。首先固定模型
u
u
u 并随机采样数据点
x
(
1
)
,
…
,
x
(
N
1
)
∈
Ω
x^{(1)},\dots,x^{(N_1)} \in \Omega
x(1),…,x(N1)∈Ω 和
x
(
1
)
,
…
,
x
(
N
2
)
∈
∂
Ω
x^{(1)},\dots,x^{(N_2)} \in \partial \Omega
x(1),…,x(N2)∈∂Ω,作为随机内循环优化的初始化。然后执行基于梯度的方法来获得具有大的逐点物理信息损失的数据点,内循环更新规则如下:
x
(
k
)
←
Project
Ω
(
x
(
k
)
+
η
sign
∇
x
(
L
u
θ
(
x
(
k
)
)
−
φ
(
x
(
k
)
)
)
2
)
x
~
(
k
)
←
Project
∂
Ω
(
x
~
(
k
)
+
η
sign
∇
x
(
B
u
θ
(
x
~
(
k
)
)
−
ψ
(
x
~
(
k
)
)
)
2
)
\begin{gathered} x^{(k)} \leftarrow\operatorname{Project}_{\Omega}\left(x^{(k)}+\eta\operatorname{sign}\nabla_{x}\left({\cal L}u_{\theta}(x^{(k)})-\varphi(x^{(k)})\right)^{2}\right) \\ \tilde{x}^{(k)} \leftarrow\operatorname{Project}_{\partial\Omega}\left(\tilde{x}^{(k)}+\eta\operatorname{sign}\nabla_x\left(\mathcal{B}u_\theta(\tilde{x}^{(k)})-\psi(\tilde{x}^{(k)})\right)^2\right) \end{gathered}
x(k)←ProjectΩ(x(k)+ηsign∇x(Luθ(x(k))−φ(x(k)))2)x~(k)←Project∂Ω(x~(k)+ηsign∇x(Buθ(x~(k))−ψ(x~(k)))2)
其中
Project
Ω
\operatorname{Project}_{\Omega}
ProjectΩ 和
Project
∂
Ω
\operatorname{Project}_{\partial\Omega}
Project∂Ω 用于将点投影到域内。当内循环结束后,就固定生成的数据点并计算梯度更新参数:
g
←
∇
θ
(
1
N
1
∑
i
=
1
N
1
(
L
u
θ
(
x
(
i
)
)
−
φ
(
x
(
i
)
)
)
2
+
λ
⋅
1
N
2
∑
i
=
1
N
2
(
B
u
θ
(
x
~
(
i
)
)
−
ψ
(
x
~
(
i
)
)
)
2
)
g\leftarrow\nabla_\theta\left(\frac{1}{N_1}\sum_{i=1}^{N_1}\left(\mathcal{L}u_\theta(x^{(i)})-\varphi (x^{(i)})\right)^2+\lambda\cdot\frac{1}{N_2}\sum_{i=1}^{N_2}\left(\mathcal{B}u_\theta(\tilde{x}^{(i)})-\psi(\tilde{x}^{(i)})\right)^2\right)
g←∇θ(N11i=1∑N1(Luθ(x(i))−φ(x(i)))2+λ⋅N21i=1∑N2(Buθ(x~(i))−ψ(x~(i)))2)
完整算法
实验结果
作者在高维线性二次高斯控制问题上进行了试验,并作了消融实验。
高维线性二次高斯控制问题
HJB方程:
{
∂
t
u
(
x
,
t
)
+
Δ
u
(
x
,
t
)
−
μ
∥
∇
x
u
(
x
,
t
)
∥
2
=
0
x
∈
R
n
,
t
∈
[
0
,
T
]
u
(
x
,
T
)
=
g
(
x
)
x
∈
R
n
,
\begin{cases}\partial_t u(x,t)+\Delta u(x,t)-\mu\|\nabla_x u(x,t)\|^2=0&x\in\mathbb{R}^n,t\in[0,T]\\ u(x,T)=g(x)&x\in\mathbb{R}^n,\end{cases}
{∂tu(x,t)+Δu(x,t)−μ∥∇xu(x,t)∥2=0u(x,T)=g(x)x∈Rn,t∈[0,T]x∈Rn,
精确解为:
u
(
x
,
t
)
=
−
1
μ
ln
(
∫
R
n
(
2
π
)
−
n
/
2
e
−
∥
y
∥
2
/
2
⋅
e
−
μ
g
(
x
−
2
(
T
−
t
)
y
)
d
y
)
u(x,t)=-\frac{1}{\mu}\ln\left(\int_{\mathbb{R}^n}(2\pi)^{-n/2}\mathrm{e}^{-\|y\|^2/2}\cdot\mathrm{e}^{-\mu g(x-\sqrt{2(T-t)}y)}\mathrm{d}y\right)
u(x,t)=−μ1ln(∫Rn(2π)−n/2e−∥y∥2/2⋅e−μg(x−2(T−t)y)dy)
其中,
μ
=
1
,
T
=
1
\mu =1, T=1
μ=1,T=1 , 成本函数为
g
(
x
)
=
ln
(
1
+
∥
x
∥
2
2
)
g(x)=\ln\left(\frac{1+\|x\|^2}{2}\right)
g(x)=ln(21+∥x∥2)
左侧为精确解 u ∗ u^∗ u∗;中间为 L 2 L^2 L2 损失的原始 PINN 方法和本文的对抗训练方法的学习解 u u u;右侧为逐点绝对误差 ∣ u − u ∗ ∣ |u − u^∗| ∣u−u∗∣。由于原方程是一个高维函数,因此这里显示的是其二维域上的可视化。具体来说,就是将函数 u ( x 1 , x 2 , 0 , ⋯ , 0 ; 0 ) u(x_1,x_2,0,\cdots,0;0) u(x1,x2,0,⋯,0;0) 中 x 1 , x 2 ∈ [ 0 , 1 ] x_1, x_2 \in [0, 1] x1,x2∈[0,1] 这两个变量可视化,其中水平轴和垂直轴分别对应于 x 1 x_1 x1 和 x 2 x_2 x2。
消融实验
总结
这篇文章从HJB方程的稳定性出发,证明了对于高维HJB方程, L p L^p Lp 范数只有当 p p p 足够大时才能确保神经网络近似解是逼近真实解的。随后文章提出使用 L ∞ L ^ \infty L∞ 范数来训练神经网络,并用对抗的思想,来近似 L ∞ L ^ \infty L∞ 范数。
这篇文章为PINN误差分析提供了新的视角,证明过程十分充足,附录有十多页的证明,可以仔细看看。
相关链接: