图核网络-史上最详细解析
Neural Operator: Graph Kernel Network for Partial Differential Equations
摘要
神经网络的经典发展主要针对有限维欧几里得空间之间的映射或两个有限维欧几里得空间之间的映射。我们工作的目的是将神经网络推广到可以学习无限维空间之间(算子)的映射。我们工作的关键创新是,利用精心设计的网络架构中的一组网络参数,可以描述无限维空间之间的映射以及这些空间的不同有限维近似之间的映射。我们通过组合非线性激活函数和一类积分算子来形式化无限维映射的近似。核积分通过图网络上的消息传递来计算。这种方法具有重要的实际意义,我们将在输入数据与偏微分方程(PDEs)及其解之间的映射的背景下进行说明。在这种背景下,这种学习的网络可以在不同的 PDE 近似方法(例如有限差分法或有限元法)之间进行泛化,并且可以在不同的分辨率和离散化级别的近似之间进行泛化。实验确认了所提出的图核网络具有期望的属性,并且显示出与最先进求解器相媲美的性能。
解释:关于有限维无限维的解释请查看link:http://t.csdnimg.cn/GFqp1
1.引言
在许多应用中,学习 Banach 空间之间的映射是非常有用的。特别是,输入空间或输出空间,或两者,可能是无限维的。学习此类映射的可能性为神经网络的设计开辟了新领域,具有广泛的潜在应用性。新想法需要基于传统神经网络的构建,这些网络是从有限维欧几里得空间映射到类别,或者映射到另一个有限维欧几里得空间。我们研究了在输入和输出空间由定义在 R d \mathbb{R}^{d} Rd子集上的实值函数组成的情况下神经网络的发展。
1.1 文献回顾与背景
我们提出了一类新的神经网络,这类神经网络被定义为在一个有界的 D ⊂ R d D\subset\mathbb{R}^d D⊂Rd 集合上映射函数空间。此类神经网络一旦训练完成,具有一个重要特性,即它们对离散化具有不变性,能够在不同的离散化之间共享相同的网络参数。相比之下,传统神经网络架构在很大程度上依赖于离散化,难以在不同的网格表示之间进行泛化。
我们的方法具有一个基于 Nystrom 近似公式的底层结构(Nystrom, 1930),它将不同的网格链接到一组网络参数上。我们在偏微分方程 (PDE) 的背景下展示了这类神经网络的新概念性应用,以及在输入数据(以函数形式)与输出数据(求解 PDE 的函数)之间的映射。我们同时考虑了监督和半监督的设置。
在 PDE 应用中,定义方程通常是局部的,而求解算子具有非局部效应,尽管这些效应会衰减。这些非局部效应可以通过 Nystrom 类型的积分算子描述,它为通过不同网格或数据结构连接提供了一致的方法,解决了计算方法中的问题。由于这一原因,图神经网络在 PDE 求解算子中展现了巨大的潜力,这也是我们工作的出发点。
解释
- 离散化不变性:传统神经网络在处理不同的网格表示时表现不佳,特别是当输入数据以不同的离散程度方式呈现时,网络难以泛化。
- Nystrom 近似:Nystrom 近似是一种用于求解积分方程的数值方法,通过采样和权重来逼近积分。作者利用 Nystrom 近似来连接不同的网格,统一表示不同的离散化方法。这种方法使得神经网络可以在不同的网格结构之间应用相同的参数,提高了网络的适应性。
- 非局部效应和图神经网络:在 PDE 应用中,方程的局部性与解算子的非局部性之间存在差异。非局部效应是指解在空间中的某个点可能会受到其他远处点的影响。
偏微分方程 (PDEs)。许多重要的工程和物理问题由 PDEs 控制。在过去的几十年里,PDE 的制定(Gurtin, 1982)和求解(Johnson, 2012)方面取得了显著进展,这些 PDE 控制了从微观尺度问题(如量子和分子动力学)到宏观尺度应用(如土木和海洋工程)的许多领域。尽管在实际应用中使用 PDEs 取得了成功,但在解决这些问题的过程中仍然面临两个主要挑战。
首先,识别或制定适合特定问题建模的基本 PDE 通常需要对相关领域有广泛的先验知识,然后结合普遍的守恒定律来设计预测模型。例如,建模固体结构的变形和断裂需要详细了解组成材料中的应力和应变之间的关系。对于复杂系统(如活细胞)来说,获得这种知识往往是困难的,制定此类系统的控制 PDE 仍然是一项巨大的挑战;从数据中学习这种知识的可能性可能会彻底改变这些领域。
第二,求解复杂的非线性 PDE 系统(如湍流和塑性问题)在计算上非常耗费资源;再次,使用来自此类计算的数据设计快速近似求解器的可能性具有巨大的潜力。在这两个挑战中,如果神经网络能够充分利用日益增多的数据量,那么就需要将它们设计为能够很好地适应函数空间之间的映射。
我们首先概述了两种主要的基于神经网络的 PDE 方法。我们考虑如下形式的 PDE:
(
L
a
u
)
(
x
)
=
f
(
x
)
,
x
∈
D
u
(
x
)
=
0
,
x
∈
∂
D
(\mathcal{L}_{a}u)(x)=f(x),\quad x\in D\\u(x)=0,\quad x\in\partial D
(Lau)(x)=f(x),x∈Du(x)=0,x∈∂D
其中,解
u
u
u是定义在
D
⊆
R
D\subseteq\mathbb{R}
D⊆R上的实数值函数,参数
a
a
a 是定义在
D
⊆
R
D\subseteq\mathbb{R}
D⊆R上的实数值函数。域
D
D
D 被离散化为
K
K
K 个点(见第 2 节),并且
N
N
N 对参数函数和(近似的)解函数
{
a
j
,
u
j
}
j
=
1
N
\{a_j,u_j\}_{j=1}^{N}
{aj,uj}j=1N被用来设计神经网络。
第一种方法是将解算子参数化为有限欧几里得空间之间的深度卷积神经网络 F : R K × Θ → R K F:\mathbb{R}^K\times\Theta\to\mathbb{R}^K F:RK×Θ→RK(Guo et al., 2016; Zhu and Zabaras, 2018; Adler and Oktem, 2017; Bhatnagar et al., 2019)。这种方法本质上并不与网格无关,需要对架构进行修改以适应不同的分辨率和离散化,从而尽可能保持一致的误差。我们在第 4 节中使用 Zhu and Zabaras (2018) 提出的架构来展示这个问题,该架构是为在 64×64 网格上求解方程设计的。此外,这种方法仅限于用于离散化空间的解,因为它不能在新点处查询解。
第二种方法 直接将解 u u u 参数化为神经网络 F : D × Θ → R \mathcal{F}:D\times\Theta\to\mathbb{R} F:D×Θ→R(E and Yu, 2018; Raissi et al., 2019; Bar and Sochen, 2019)。这种方法显然是网格无关的,因为解是定义在物理域上的。然而,参数依赖性体现在与网格相关的方式上。实际上,对于给定的新方程和新系数函数 a a a,需要重新训练一个新的神经网络 F a \mathcal{F}_a Fa 。这种方法与经典方法(如有限元法)非常相似,后者用有限集局部基函数的线性组合替代神经网络空间。这种方法存在与经典方法相同的计算问题:对于每个新参数,都需要解决一个优化问题。此外,该方法仅限于在已知的 PDE 情形中使用,而当函数之间的映射是完全数据驱动的学习时,这种方法不可行。我们提出的方法解决了这些问题。
我们的方法最接近经典的简化基方法(DeVore, 2014)或 Cohen 和 DeVore 的方法(Cohen and DeVore, 2015)。与当代的相关工作(Bhattacharya et al., 2020)一起,据我们所知,这是第一个实用的深度学习方法,专门用于学习无限维空间之间的映射。它解决了 Guo et al., 2016 和 Zhu and Zabaras, 2018 提出的网格相关问题,并通过产生与函数分辨率无关的近似质量,从而使解能够在网格之间转移。此外,它只需要对方程组 { a j , u j } j = 1 N \{a_j,u_j\}_{j=1}^{N} {aj,uj}j=1N行一次训练;随后对于新参数 a ∼ μ a\sim\mu a∼μ 的解。只需要进行一次网络的前向传播,从而缓解了 E and Yu, 2018 和 Raissi et al., 2019 所面临的主要计算问题。
最后,我们的方法不需要了解潜在的 PDE;实际上,映射 F 可以被视为一个黑盒,可能基于实验数据或昂贵计算机模拟的结果进行训练,而不必是 PDE 的形式。
图神经网络 (GNNs):图神经网络 (GNNs) 是一类作用于图结构数据的神经网络,近年来得到了广泛的发展并应用于各种领域。图神经网络包含了诸如图卷积、边卷积、注意力机制和图池化等技术(Kipf and Welling, 2016;Hamilton et al., 2017)。
GNNs 还被应用于物理现象的建模,如分子(Chen et al., 2019)和刚体系统(Battaglia et al., 2018),因为这些问题天然具备图结构的解释:粒子是节点,交互作用是边。Alet et al., 2019 的工作进行了初步研究,使用图网络来学习泊松方程的解,以及其他物理应用。他们提出了一种编码器-解码器设置,在潜在空间中构建图,并使用编码器和解码器之间的消息传递。然而,他们的模型使用了最近邻结构,当网格大小增加时,无法捕捉到非局部依赖性。
与之相反,我们直接构造一个图,节点位于输出函数的空间域上。通过消息传递,我们可以直接学习近似 PDE 解的网络核函数。当查询一个新位置时,我们只需将一个新节点添加到我们的空间图中,并将其连接到现有节点,通过利用 Nystrom 扩展避免插值误差。
连续神经网络 (Continuous Neural Networks):在无限维空间中定义神经网络的概念是一个长期研究的核心问题(Williams, 1996;Neal, 1996;Roux and Bengio, 2007;Globerson and Livni, 2016;Guss, 2016)。通常的思路是取无限宽极限,这将产生一种非参数方法,并与高斯过程回归有关联(Neal, 1996;Matthews et al., 2018;Garriga-Alonso et al., 2018),引入了深度高斯过程(Damianou and Lawrence, 2013;Dunlop et al., 2018)。然而,这些方法尚未产生可以与卷积神经网络或循环神经网络的成功相媲美的高效数值算法。
另一种思路是定义每一层为无限维空间之间的映射,具有有限维的参数依赖。这是我们在本工作中采用的方法,通过在每层之间共享参数,进一步推进这一思想。
1.2 贡献
我们引入了神经算子的概念,并通过图核网络(graph kernel networks)实现它,这是一种新型的深度神经网络方法,用于学习定义在 R D \mathbb{R}^D RD 上有界开集的无限维函数空间之间的映射。
-
与现有方法不同,我们的方法能够在不同的近似和网格之间学习函数空间之间的映射,这在图 1.1 中有所展示。
-
我们使用了 Nystrom 扩展,将函数空间上的神经网络与基于任意、可能非结构化网格的 GNN 系列连接起来。
-
我们展示了神经算子方法具有与经典方法和深度学习方法相当的近似精度。
-
数值结果还表明,神经算子只需要在少量样本上训练,就能泛化到整个问题类。
-
我们展示了进行半监督学习的能力,能够从少量点的数据中学习并泛化到整个域。
这些概念在一类典型的椭圆型偏微分方程(PDE)中得到了体现,该类问题广泛出现在科学和工程领域。
2. 问题设置
我们的目标是通过有限的输入-输出对集合,学习两个无限维空间之间的映射:监督学习。令 A A A 和 U U U 是可分的 Banach 空间,且 F † : A → U \mathcal{F}^{\dagger}:A\to U F†:A→U是一个(通常)非线性映射。假设我们有观察值 { a j , u j } j = 1 N \{a_j,u_j\}_{j=1}^{N} {aj,uj}j=1N,其中a_{j}\sim\mu 是从定义在 A A A 上的概率测度 μ μ μ 中独立同分布的序列,且 u j = F † ( a j ) u_j=\mathcal{F}^\dagger(a_j) uj=F†(aj)可能被噪声污染。我们旨在通过构建一个参数化映射 F : A × Θ → U ( 2 ) \mathcal{F}:A\times\Theta\to U(2) F:A×Θ→U(2)来逼近 F † \mathcal{F}^\dagger F†,其中 Θ \Theta Θ是一个有限维参数空间,通过选择 θ ∗ ∈ Θ \theta^{*}\in\Theta θ∗∈Θ使得 F ( ⋅ , θ ∗ ) ≈ F † \mathcal{F}(\cdot,\theta^*)\approx\mathcal{F}^\dagger F(⋅,θ∗)≈F†。
这是一个在无限维空间中学习的自然框架,因为我们可以定义一个代价函数
C
:
U
×
U
→
R
C:U×U→R
C:U×U→R,并求解问题的最小化:
min
θ
∈
Θ
E
a
∼
μ
[
C
(
F
(
a
,
θ
)
,
F
†
(
a
)
)
]
\min_{\theta\in\Theta}\mathbb{E}_{a\sim\mu}[C(\mathcal{F}(a,\theta),\mathcal{F}^\dagger(a))]
θ∈ΘminEa∼μ[C(F(a,θ),F†(a))]
这一问题与经典的有限维问题完全对应(Vapnik, 1998)。在无限维环境中,展示最小化解的存在性仍然是一个具有挑战性的开放问题。我们将通过训练-测试框架来处理这个问题,其中经验近似被用来计算代价。我们将我们的方法概念化为允许所有有限维近似共享一组网络参数,并在无限维环境中定义(无近似自由的)问题。
具体而言,我们将考虑 Banach 空间上的映射,它们由定义在
R
d
R^d
Rd
上有界开集中的实值函数组成。这些映射将输入函数转换为 PDE 的输入,并将它们映射到 PDE 的解,其中输入和解都是定义在
R
d
R^d
Rd上的实值函数。
一个常见的此类问题是对二阶椭圆 PDE 的逼近:
− ∇ ⋅ ( a ( x ) ∇ u ( x ) ) = f ( x ) , x ∈ D ( 3 ) -\nabla \cdot (a(x) \nabla u(x)) = f(x), \quad x \in D (3) −∇⋅(a(x)∇u(x))=f(x),x∈D(3) u ( x ) = 0 , x ∈ ∂ D u(x) = 0, \quad x \in \partial D u(x)=0,x∈∂D
对于某个有界开集 D ⊂ R d D \subset \mathbb{R}^d D⊂Rd,和一个固定函数 f ∈ L 2 ( D ; R ) f \in L^2(D; \mathbb{R}) f∈L2(D;R),该方程是众多应用中典型的 PDE,包括水文学和弹性力学的步伐(Bear and Corapcioglu, 2012; Antman, 2005)。对于给定的 a ∈ A = L ∞ ( D ; R + ) ∩ L 2 ( D ; R + ) a \in A = L^\infty(D; \mathbb{R}^+) \cap L^2(D; \mathbb{R}^+) a∈A=L∞(D;R+)∩L2(D;R+),方程 (3) 有唯一的弱解 u ∈ U = H 0 1 ( D ; R ) u \in U = H^1_0(D; \mathbb{R}) u∈U=H01(D;R),因此我们可以将解算子 F † \mathcal{F}^{\dagger} F† 定义为映射 a ↦ u a \mapsto u a↦u。尽管 PDE (3) 是线性的,解算子 F † \mathcal{F}^{\dagger} F† 并非如此。
由于我们的数据 a j a_j aj 和 u j u_j uj 通常是函数,为了数值计算,我们假设只能访问点评估数据。为了说明这一点,我们将继续使用上一段的例子,定义 P K = x 1 , … , x K ⊂ D P_K = {x_1, \dots, x_K} \subset D PK=x1,…,xK⊂D,这是域 D D D 的一个 K 点离散化,且假设我们有观测值 a j ∣ P K , u j ∣ P K ∈ R K a_j|{P_K}, u_j|{P_K} \in \mathbb{R}^K aj∣PK,uj∣PK∈RK,这是通过有限的输入-输出对集合来获得的。在下一节中,我们提出了一种基于核的图神经网络架构,该架构在离散化数据上训练,但能够在给定的新输入 a ∼ μ a \sim \mu a∼μ 时产生 u ( x ) u(x) u(x) 的答案。也就是说,我们的方法与离散化 P K P_K PK 无关,并且是一个真正的函数空间方法;我们通过数值方式验证了这一点,展示了当 K → ∞ K→∞ K→∞ 时误差不变的性质。这种性质非常理想,因为它允许在不同网格几何和离散化大小之间传递解。
我们注意到,虽然我们的方法的应用基于函数的点评估,但并不限于此。例如,可以将函数数值化为有限集截断基函数,表示函数。表示的不变性将相对于这一集的大小。我们的方法原则上可以通过适当选择的架构进行修改以适应这种情况,但我们在本工作中不追求这一方向。
解读:
-
− ∇ ⋅ ( a ( x ) ∇ u ( x ) ) = f ( x ) , x ∈ D -\nabla \cdot (a(x) \nabla u(x)) = f(x), \quad x \in D −∇⋅(a(x)∇u(x))=f(x),x∈D 这个倒三角符号
∇ ∇ ∇ 被称为nabla 运算符,也叫做梯度运算符,在不同的上下文中它可以代表不同的运算。具体到你提供的这个公式,它表示的是散度和梯度的组合运算。具体含义如下:- 梯度
∇
u
(
x
)
∇u(x)
∇u(x):
梯度运算符 ∇ ∇ ∇ 应用于标量函数 u ( x ) u(x) u(x) 时,表示函数在各个方向上的变化率,生成的是一个向量场。它表示标量函数的变化率,指向函数增长最快的方向,大小表示变化的速率。 - 散度
∇
⋅
(
向量场
)
∇⋅(向量场)
∇⋅(向量场):
散度运算符 ∇ ⋅ ∇⋅ ∇⋅ 应用于一个向量场时,表示该向量场在某一点的发散或汇聚程度。它给出了一个标量值,用于描述向量场的体积膨胀率或收缩率。
- 梯度
∇
u
(
x
)
∇u(x)
∇u(x):
-
L ∞ ( D ; R + ) ∩ L 2 ( D ; R + ) L^\infty(D; \mathbb{R}^+) \cap L^2(D; \mathbb{R}^+) L∞(D;R+)∩L2(D;R+) 的意思是函数空间的交集,具体解释如下:
- L ∞ ( D ; R + ) L^\infty(D; \mathbb{R}^+) L∞(D;R+):表示定义在区域 D D D上、值域为非负实数(即 R + R^+ R+ )的本质有界函数的空间。即,函数在 D D D 上的值不会无限大,它的上界存在并且是有限的。
- L 2 ( D ; R + ) L^2(D; \mathbb{R}^+) L2(D;R+)表示定义在区域 D D D 上、值域为非负实数的平方可积函数的空间。即,函数在 D D D 上的平方是可积的,也就是满足: ∫ D ∣ f ( x ) ∣ 2 d x < ∞ \int_D|f(x)|^2 dx<\infty ∫D∣f(x)∣2dx<∞
-
数值计算中的点评估:由于函数 a j a_j aj和 u j u_j uj 是连续定义的,但我们通常只能访问这些函数在某些离散点的值,即所谓的点评估。
3.图核网络(Graph Kernel Network)
我们提出了一种图核神经网络,用于解决第2节中描述的问题。作为我们架构的指导原则,我们采用了如下例子。设
L
a
\mathcal{L}_{a}
La是一个基于参数
a
∈
A
a∈A
a∈A 的微分算子,考虑如下的 PDE:
(
L
a
u
)
(
x
)
=
f
(
x
)
,
x
∈
D
(
4
)
(\mathcal{L}_a u)(x) = f(x), \quad x \in D(4)
(Lau)(x)=f(x),x∈D(4)
u
(
x
)
=
0
,
x
∈
∂
D
u(x) = 0, \quad x \in \partial D
u(x)=0,x∈∂D
对于有界开集
D
⊂
R
d
D⊂R^d
D⊂Rd 和一个固定函数
f
f
f,其定义在由
L
a
\mathcal{L}_{a}
La 结构确定的适当函数空间中。来自公式 (3) 的椭圆算子 是一个例子。根据
L
a
\mathcal{L}_{a}
La 的通用条件(Evans, 2010),我们可以定义 Green 函数
G
:
D
×
D
→
R
G:D×D→R
G:D×D→R 作为问题的唯一解:
L
a
G
(
x
,
⋅
)
=
δ
x
\mathcal{L}_a G(x, \cdot) = \delta_x
LaG(x,⋅)=δx
其中,
δ
x
\delta_x
δx 是以
x
x
x 为中心的
R
d
\mathbb{R}^d
Rd 上的 Dirac 函数。注意,
G
G
G 将取决于参数
a
a
a,因此我们将其记为
G
a
G_a
Ga。公式 (4) 的解可以表示为:
u
(
x
)
=
∫
D
G
a
(
x
,
y
)
f
(
y
)
d
y
(
5
)
u(x) = \int_D G_a(x, y) f(y) \, dy(5)
u(x)=∫DGa(x,y)f(y)dy(5)
这可以通过以下的形式计算轻松看出:
(
L
a
u
)
(
x
)
=
∫
D
(
L
a
G
(
x
,
⋅
)
)
(
y
)
f
(
y
)
,
d
y
=
∫
D
δ
x
(
y
)
f
(
y
)
,
d
y
=
f
(
x
)
(\mathcal{L}_a u)(x) = \int_D (\mathcal{L}_a G(x, \cdot))(y) f(y) , dy = \int_D \delta_x(y) f(y) , dy = f(x)
(Lau)(x)=∫D(LaG(x,⋅))(y)f(y),dy=∫Dδx(y)f(y),dy=f(x)
通常,Green 函数在
x
≠
y
x \neq y
x=y 的点是连续的。例如,当
L
a
\mathcal{L}_a
La 是一致椭圆时(Gilbarg and Trudinger, 2015),因此将其通过神经网络建模是自然的。受公式 (5) 的启发,我们提出了如下的迭代架构
⋅
t
=
0
,
…
,
T
−
1.
\cdot t=0,\ldots,T-1.
⋅t=0,…,T−1.。
v
t
+
1
(
x
)
=
σ
(
W
v
t
(
x
)
+
∫
D
κ
ϕ
(
x
,
y
,
a
(
x
)
,
a
(
y
)
)
v
t
(
y
)
ν
x
(
d
y
)
)
(
6
)
v_{t+1}(x)=\sigma\left(Wv_t(x)+\int_D\kappa_\phi(x,y,a(x),a(y))v_t(y) \nu_x(dy)\right)(6)
vt+1(x)=σ(Wvt(x)+∫Dκϕ(x,y,a(x),a(y))vt(y)νx(dy))(6)
其中,
σ
:
R
→
R
\sigma: \mathbb{R} \to \mathbb{R}
σ:R→R 是逐元素应用的固定函数,
ν
x
\nu_x
νx 是每个
x
∈
D
x \in D
x∈D 上的固定 Borel 测度。矩阵
W
∈
R
n
×
n
W \in \mathbb{R}^{n \times n}
W∈Rn×n,以及进入核函数
κ
ϕ
\kappa_\phi
κϕ 的参数
ϕ
\phi
ϕ,定义为
κ
ϕ
:
R
2
(
d
+
1
)
→
R
n
×
n
\kappa_\phi: \mathbb{R}^{2(d+1)} \to \mathbb{R}^{n \times n}
κϕ:R2(d+1)→Rn×n,需要从数据中学习。我们将
κ
ϕ
\kappa_\phi
κϕ 建模为一个将
R
2
(
d
+
1
)
\mathbb{R}^{2(d+1)}
R2(d+1) 映射到
R
n
×
n
\mathbb{R}^{n \times n}
Rn×n 的神经网络。
离散化与参数共享 离散化的连续图景可以视为通过基于 K 点的经验逼近来替代 Borel 测度 ν x \nu_x νx。在这种设置下,我们可以将 κ ϕ \kappa_\phi κϕ 视为一个 K × K K \times K K×K 的核块矩阵,其中每个条目 κ ϕ ( x , y ) \kappa_\phi(x, y) κϕ(x,y) 本质上是一个 n × n n \times n n×n 矩阵。每个块共享相同的网络参数。这是使该方法共享独立于所用离散化的公共参数的关键。
最后,我们注意到,尽管我们专注于将神经网络映射从 a a a 映射到 u u u,但也可以进行广泛的推广,例如将 f f f 映射到 u u u,或者处理在 ∂ D \partial D ∂D 上具有非零边界数据 g g g 并将 g g g 映射到 u u u。为了展示我们的想法,我们将首先考虑从 f f f 映射到 u u u 的情况(线性且已知解析解),然后再研究从 a a a 映射到 u u u 的(非线性)映射。
解读:
- Dirac 函数(又称 Dirac Delta 函数,或称 狄拉克δ函数)是一种在数学和物理学中常用的特殊函数,通常表示为
δ
(
x
)
δ(x)
δ(x)。Dirac 函数的定义如下:
- 零值性:在除了 x = 0 x=0 x=0 之外的所有点上, δ ( x ) = 0 δ(x)=0 δ(x)=0。
- 积分为 1:Dirac 函数在
x
=
0
x=0
x=0 处有一个无限高的尖峰,并满足以下性质:
∫ − ∞ ∞ δ ( x ) , d x = 1 \int_{-\infty}^{\infty} \delta(x) , dx = 1 ∫−∞∞δ(x),dx=1
简而言之,Dirac 函数可以被看作是一个“单位脉冲”,它在 x = 0 x=0 x=0 处集中了所有的值,在该点之外为零。它常用于描述一个瞬间的冲击或极端集中的分布。
- Green 函数的作用:
- Green 函数 G a ( x , y ) G_a(x,y) Ga(x,y) 是一种核函数,它表示在点 y y y处的一个冲击(即 Dirac 函数 δ y δ_y δy )如何影响系统在点 x x x 处的响应。通过引入 Green 函数,复杂的偏微分方程问题可以转换为积分方程的形式。
- 含义:给定偏微分方程 L a u ( x ) = f ( x ) \mathcal{L}_a u(x) = f(x) Lau(x)=f(x),可以将其理解为求解输入 f ( x ) f(x) f(x) 对应的系统响应 u ( x ) u(x) u(x)。而通过 Green 函数,系统的响应可以通过对所有点的叠加效应计算出来,这些效应由积分表示。 - 具体地, G a ( x , y ) G_a(x, y) Ga(x,y) 表示系统在位置 x x x 对位于 y y y 处的输入 f ( y ) f(y) f(y) 的响应。通过对 D D D 上所有 y y y 点进行积分,可以累积所有位置对 x x x 的贡献,最终得到 u ( x ) u(x) u(x) 的值。
- 计算过程:利用 Green 函数的定义
L
a
G
(
x
,
y
)
=
δ
x
(
y
)
\mathcal{L}_a G(x, y) = \delta_x(y)
LaG(x,y)=δx(y),将偏微分方程
L
a
u
(
x
)
=
f
(
x
)
\mathcal{L}_a u(x) = f(x)
Lau(x)=f(x) 转化为积分形式。Green 函数的定义使得可以使用 Dirac 函数的性质:
( L a u ) ( x ) = ∫ D ( L a G ( x , ⋅ ) ) ( y ) f ( y ) , d y = ∫ D δ x ( y ) f ( y ) , d y = f ( x ) (\mathcal{L}_a u)(x) = \int_D (\mathcal{L}_a G(x, \cdot))(y) f(y) , dy = \int_D \delta_x(y) f(y) , dy = f(x) (Lau)(x)=∫D(LaG(x,⋅))(y)f(y),dy=∫Dδx(y)f(y),dy=f(x)
因此,Green 函数的积分形式解 u ( x ) u(x) u(x) 正是系统的响应,它将输入 f ( y ) f(y) f(y) 的影响通过 Green 函数 G a ( x , y ) G_a(x, y) Ga(x,y) 传播到输出 u ( x ) u(x) u(x)。
示例 1:泊松方程(Poisson Equation) 我们考虑一个简化的泊松方程映射的情况,其中我们将
f
f
f 映射到
u
u
u。为此,我们设定
v
0
(
x
)
=
f
(
x
)
v_0(x) = f(x)
v0(x)=f(x),
T
=
1
T = 1
T=1,
n
=
1
n = 1
n=1,
σ
(
x
)
=
x
\sigma(x) = x
σ(x)=x,
W
=
w
=
0
W = w = 0
W=w=0,以及
ν
x
(
d
y
)
=
d
y
\nu_x(dy) = dy
νx(dy)=dy(勒贝格测度)。我们然后通过神经网络
κ
ϕ
\kappa_\phi
κϕ 获取 Green 函数
G
a
G_a
Ga,其依赖于
a
(
x
)
a(x)
a(x) 和
a
(
y
)
a(y)
a(y)。现在考虑情景
D
=
[
0
,
1
]
D = [0, 1]
D=[0,1] 且
a
(
x
)
≡
1
a(x) \equiv 1
a(x)≡1,此时公式 (3) 简化为可以显式计算的 Green 函数:
G
(
x
,
y
)
=
1
2
(
x
+
y
−
∣
y
−
x
∣
)
−
x
y
G(x, y) = \frac{1}{2}(x + y - |y - x|) - xy
G(x,y)=21(x+y−∣y−x∣)−xy
请注意,虽然映射
f
↦
u
f \mapsto u
f↦u 在函数空间中是线性的,但 Green 函数本身在两个变量中都不是线性的。图 2 显示了在
N
=
2048
N = 2048
N=2048 个样本
f
j
∼
μ
=
N
(
0
,
(
−
Δ
+
I
)
−
1
)
f_j \sim \mu = \mathcal{N}(0, (-\Delta + I)^{-1})
fj∼μ=N(0,(−Δ+I)−1) 训练之后的
κ
ϕ
\kappa_\phi
κϕ,其中在算子
−
Δ
+
I
-\Delta + I
−Δ+I 上具有周期边界条件。(可以通过 Karhunen-Loeve 的随机傅里叶级数简单地实现该测度的样本——参见 [Lord et al., 2014])。
解读:
-
这里提到在算子 − Δ + I -\Delta + I −Δ+I 上有周期边界条件,这指的是训练过程中使用的偏微分算子条件。周期边界条件意味着在边界处的函数值在一定意义上是相互“连接”的,这种条件常用于物理模型中的循环系统。
-
使用了 2048 个样本 f j f_j fj,每个样本都从一个高斯分布 N ( 0 , ( − Δ + I ) − 1 ) \mathcal{N}(0, (-\Delta + I)^{-1}) N(0,(−Δ+I)−1) 中采样。这个分布带有特定的协方差结构,包含了算子 − Δ + I -\Delta + I −Δ+I 的逆。
-
Karhunen-Loeve 展开是一种在信号处理和随机过程中的工具,用于分解随机过程为一组正交函数(类似于傅里叶展开)。这里提到可以通过 Karhunen-Loeve 展开来实现这些随机样本的生成。这个方法能有效地将随机过程表示为一系列不相关的成分,使得从该测度生成样本变得简单可行。
示例 2:二维泊松方程:我们通过将第 1 个示例中研究的泊松方程扩展到二维 (2D) 情况,进一步展示了图核网络的强大功能。在这里我们逼近映射 f ( x ) ↦ u ( x ) f(x) \mapsto u(x) f(x)↦u(x),其中 x ∈ D = [ 0 , 1 ] × [ 0 , 1 ] x \in D = [0, 1] \times [0, 1] x∈D=[0,1]×[0,1]。我们考虑两种方法:1) 使用图核网络逼近二维 Green 函数 G ( x , y ) G(x, y) G(x,y);2) 直接设计神经网络,以 f ( x ) f(x) f(x) 作为输入, u ( x ) u(x) u(x) 作为输出,从而逼近映射 f ( x ) ↦ u ( x ) f(x) \mapsto u(x) f(x)↦u(x)。
这两个神经网络使用相同的训练集进行训练,训练集的大小从 1 个样本到 100 个样本不等,并使用 1000 个测试样本进行测试。训练样本的相对 l 2 l_2 l2 测试误差如图 A.4 所示。我们观察到,图核网络使用最小数量(约 5 个)训练样本即可逼近该映射,同时与使用 100 个样本训练的密集神经网络相比,测试误差更小。因此,图核网络显著减少了逼近映射所需的训练样本数量。该性质在实际应用中尤其重要,因为获取大量训练数据对某些工程/物理问题来说通常是困难的。
图核网络具有如此强大逼近能力的原因是,它能够学习到泊松方程的真实 Green 函数,正如我们在一维案例中已经展示的那样。我们将在附录 A.4 中进一步展示图核网络如何能够捕捉到二维泊松方程的真实核函数。
算法框架:我们将初始化 v 0 ( x ) v_0(x) v0(x) 传入我们的网络 (6),可以看作是我们对解 u ( x ) u(x) u(x) 的初步猜测,以及我们想要显式表达的任何其他依赖性。一个自然的选择是从系数 a ( x ) a(x) a(x) 本身以及物理空间中的位置 x x x 开始。该 (d + 1) 维向量场被提升为 n 维向量场的形式,这可以视为整体神经网络的第一层。这随后被用作核神经网络的初始化,核神经网络会进行 T 次迭代。在最后一层中,我们通过另一个神经网络层投影回我们感兴趣的标量场。
由于逆椭圆算子 (3) 对输入数据
a
a
a(实际上是对
f
f
f)的平滑效果,我们用高斯平滑版本的系数
a
c
(
x
)
a_c(x)
ac(x) 及其梯度
∇
a
c
(
x
)
\nabla a_c(x)
∇ac(x) 来扩展初始化 (
x
,
a
(
x
)
x, a(x)
x,a(x))。因此,我们使用一个 2(d + 1) 维向量场进行初始化。本文中高斯平滑是通过方差为 5 的居中各向同性高斯分布来执行的。Borel 测度
ν
x
\nu_x
νx 被选择为支持在以
x
x
x 为中心、半径为
r
r
r 的球上的 Lebesgue 测度。因此我们有:
v
0
(
x
)
=
P
(
x
,
a
(
x
)
,
a
c
(
x
)
,
∇
a
c
(
x
)
)
+
p
(
7
)
v_0(x) = P(x, a(x), a_c(x), \nabla a_c(x)) + p(7)
v0(x)=P(x,a(x),ac(x),∇ac(x))+p(7)
v
t
+
1
(
x
)
=
σ
(
W
v
t
(
x
)
+
∫
B
(
x
,
r
)
κ
ϕ
(
x
,
y
,
a
(
x
)
,
a
(
y
)
)
v
t
(
y
)
ν
x
(
d
y
)
)
(
8
)
v_{t+1}(x) = \sigma \left(W v_t(x) + \int_{B(x, r)} \kappa_\phi(x, y, a(x), a(y)) v_t(y) \, \nu_x(dy) \right)(8)
vt+1(x)=σ(Wvt(x)+∫B(x,r)κϕ(x,y,a(x),a(y))vt(y)νx(dy))(8)
u
(
x
)
=
Q
v
T
(
x
)
+
q
(
9
)
u(x) = Q v_T(x) + q(9)
u(x)=QvT(x)+q(9)
其中, P ∈ R n × 2 ( d + 1 ) P \in \mathbb{R}^{n \times 2(d+1)} P∈Rn×2(d+1), p ∈ R n p \in \mathbb{R}^n p∈Rn, v t ( x ) ∈ R n v_t(x) \in \mathbb{R}^n vt(x)∈Rn, Q ∈ R 1 × n Q \in \mathbb{R}^{1 \times n} Q∈R1×n, q ∈ R q \in \mathbb{R} q∈R。公式 (8) 中的积分通过蒙特卡罗求和近似,利用一个消息传递图网络,其中边权重为 ( x , y , a ( x ) , a ( y ) ) (x, y, a(x), a(y)) (x,y,a(x),a(y))。选择测度 ν x ( d y ) = 1 B ( x , r ) d y \nu_x(dy) = 1_{B(x,r)}dy νx(dy)=1B(x,r)dy 有两个作用:1)它允许更高效的计算,2)它利用了 Green 函数的衰减特性。请注意,如果我们对真实核有更多信息,则可以将其加入该测度中。例如,如果我们知道真实核具有高斯结构,我们可以定义 ν x ( d y ) = 1 B ( x , r ) ρ x ( y ) d y \nu_x(dy) = 1_{B(x,r)} \rho_x(y) dy νx(dy)=1B(x,r)ρx(y)dy,其中 ρ x ( y ) \rho_x(y) ρx(y) 是高斯密度。此时, κ ϕ \kappa_\phi κϕ 将只需学习一个相对简单的函数。然而,在本研究中我们没有深入探讨该方向。
消息传递图网络:
消息传递图网络采用了一种标准架构,使用边特征 [Gilmer et al., 2017]。如果我们在 PDE 的物理域
D
D
D 上正确构造了图,则核积分可以看作是消息的聚合。给定节点特征
v
t
(
x
)
∈
R
n
v_t(x) \in \mathbb{R}^n
vt(x)∈Rn,边特征
e
(
x
,
y
)
∈
R
n
e
e(x, y) \in \mathbb{R}^{n_e}
e(x,y)∈Rne 和图
G
G
G,消息传递神经网络的聚合形式为:
v t + 1 ( x ) = σ ( W v t ( x ) + 1 ∣ N ( x ) ∣ ∑ y ∈ N ( x ) κ ϕ ( e ( x , y ) ) v t ( y ) ) v_{t+1}(x) = \sigma \left( W v_t(x) + \frac{1}{|N(x)|} \sum_{y \in N(x)} \kappa_\phi (e(x, y)) v_t(y) \right) vt+1(x)=σ Wvt(x)+∣N(x)∣1y∈N(x)∑κϕ(e(x,y))vt(y)
其中 W ∈ R n × n W \in \mathbb{R}^{n \times n} W∈Rn×n, N ( x ) N(x) N(x) 是根据图的邻域, κ ϕ ( e ( x , y ) ) \kappa_\phi (e(x, y)) κϕ(e(x,y)) 是一个神经网络,将输入边特征作为输入并输出一个矩阵 R n × n \mathbb{R}^{n \times n} Rn×n。关于公式 (8),边特征 e ( x , y ) = ( x , y , a ( x ) , a ( y ) ) ∈ R 2 ( d + 1 ) e(x, y) = (x, y, a(x), a(y)) \in \mathbb{R}^{2(d+1)} e(x,y)=(x,y,a(x),a(y))∈R2(d+1)。
图的构建: 为了使用消息传递框架(公式 (10)),我们需要设计一个将 PDE 的物理域 D D D 连接的图。节点被选择为离散化后的 K K K 个空间位置。这里我们使用的是标准的均匀网格,但还有其他可能性,例如有限元网格和随机点数据获取。边的连通性则根据积分测度(公式 (8))进行选择,即 Lebesgue 测度限制在一个球体内。每个节点 x ∈ R d x \in \mathbb{R}^d x∈Rd 连接到距离其在球体 B ( x , r ) B(x, r) B(x,r) 内的节点,定义邻域 N ( x ) N(x) N(x)。然后对于每个邻居 y ∈ N ( x ) y \in N(x) y∈N(x),我们分配边权重 e ( x , y ) = ( x , y , a ( x ) , a ( y ) ) e(x, y) = (x, y, a(x), a(y)) e(x,y)=(x,y,a(x),a(y))。 公式 (10) 可以视为公式 (3) 的蒙特卡罗近似。这样的局部结构允许了更加高效的计算,同时仍保持对网格细化的独立性。事实上,由于半径参数 r r r 是在物理空间中选择的,随着离散化大小 K K K 的增加,邻域集合 N ( x ) N(x) N(x) 的大小也会增加。这是我们方法中使其网格独立的关键特性。
Nyström 核近似:
尽管前述的图结构极大地减少了在整个域
D
D
D 上积分(相当于一个完全连接的图)的计算负担,边的数量仍然随
K
2
K^2
K2 规模增长。为了解决这一问题,我们采用了随机 Nyström 类型的核近似。特别是,我们均匀地从原始图中采样远小于
K
K
K 的节点数,构造新的随机子图。这一过程被重复
l
=
N
l = N
l=N 次,每次生成
l
l
l 个带有
m
m
m 个节点的随机子图。这可以看作是一种减少估计方差的方式。我们在训练时使用这些子图来评估 (10),从而实现了更有利的缩放 (
l
(
m
2
)
l(m^2)
l(m2))。在评估阶段,当我们希望在特定网格上求解时,我们将网格划分为带有
m
m
m 个节点的子图并分别进行计算。
我们现在将在 RKHS(Reproducing Kernel Hilbert Space,重现核希尔伯特空间)设置中突出显示这种核近似。一个真实的重现核希尔伯特空间 H \mathcal{H} H 是一个函数空间 f : D → R f : D \to \mathbb{R} f:D→R,其中在每个点上的函数值可以通过连续线性泛函表示,即 ∣ f ( x ) ∣ ≤ C ∣ f ∣ |f(x)| \leq C |f| ∣f(x)∣≤C∣f∣,对于某些常数 C ≥ 0 C \geq 0 C≥0,且与 x x x 无关。对于每个 RKHS,都存在唯一的、对称的、正定核函数 κ : D × D → R \kappa : D \times D \to \mathbb{R} κ:D×D→R,它给出了表示 f ( x ) = ⟨ f , κ ( ⋅ , x ) ⟩ f(x) = \langle f, \kappa(\cdot, x) \rangle f(x)=⟨f,κ(⋅,x)⟩。
令
T
:
H
→
H
T : \mathcal{H} \to \mathcal{H}
T:H→H 为通过核作用在
H
\mathcal{H}
H 上的线性算子:
T
f
=
∫
B
(
x
,
r
)
κ
(
x
,
y
)
f
(
y
)
ν
(
d
y
)
T f = \int_{B(x,r)} \kappa(x, y) f(y) \nu(dy)
Tf=∫B(x,r)κ(x,y)f(y)ν(dy)
令 T m : H → H T_m : \mathcal{H} \to \mathcal{H} Tm:H→H 为它的 m 点经验近似:
T m f = ∫ B ( x , r ) κ ( x , y ) f ( y ) ν m ( d y ) T_m f = \int_{B(x,r)} \kappa(x, y) f(y) \nu_m(dy) Tmf=∫B(x,r)κ(x,y)f(y)νm(dy)
其中
ν m ( d y ) = 1 m ∑ k = 1 m δ y k ( d y ) \nu_m(dy) = \frac{1}{m} \sum_{k=1}^m \delta_{y_k}(dy) νm(dy)=m1k=1∑mδyk(dy)
从而
T m f = 1 m ∑ k = 1 m κ ( x , y k ) f ( y k ) T_m f = \frac{1}{m} \sum_{k=1}^m \kappa(x, y_k) f(y_k) Tmf=m1k=1∑mκ(x,yk)f(yk)
该近似的误差实现了蒙特卡罗率 ¥ O ( m − 1 / 2 ) ¥ O(m^{-1/2})¥ O(m−1/2)¥。
命题 1: 假设
E
y
∼
ν
[
κ
(
y
,
y
)
4
]
<
∞
\mathbb{E}_y \sim \nu [\kappa(y, y)^4] < \infty
Ey∼ν[κ(y,y)4]<∞,则存在常数
C
≥
0
C \geq 0
C≥0 使得:
E
∣
T
−
T
m
∣
HS
≤
C
m
¥
\mathbb{E} |T - T_m|_{\text{HS}} \leq \frac{C}{\sqrt{m}}¥
E∣T−Tm∣HS≤mC¥
其中
∣
⋅
∣
HS
|\cdot|_{\text{HS}}
∣⋅∣HS 表示作用于
H
\mathcal{H}
H 上的算子的 Hilbert-Schmidt 范数。证明请见附录 A.2。在更严格的假设下,类似的结果也可以通过高概率证明 [Rosasco et al., 2010]。
引理 2:
(引自 Belkin)
T
T
T 和
T
K
T_K
TK 是 Hilbert-Schmidt 算子。此外,以概率大于
1
−
2
e
−
τ
1 - 2 e^{-\tau}
1−2e−τ 满足:
∣
T
−
T
K
∣
HS
≤
2
2
k
τ
K
¥
|T - T_K|_{\text{HS}} \leq \frac{2 \sqrt{2k} \sqrt{\tau}}{\sqrt{K}}¥
∣T−TK∣HS≤K22kτ¥
其中
k
=
sup
x
∈
D
κ
(
x
,
x
)
k = \sup_{x \in D} \kappa(x, x)
k=supx∈Dκ(x,x)。
我们注意到,在我们的算法中,核 κ : D × D → R n × n \kappa : D \times D \to \mathbb{R}^{n \times n} κ:D×D→Rn×n,而前述结果仅在 n = 1 n = 1 n=1 的情形下证明;尽管如此,这些结果为我们方法中的近似提供了有用的直觉。
5 结论
讨论与未来工作
我们引入了神经算子的概念,并通过图核网络实现,它旨在逼近函数空间之间的映射。这些网络被构造为无网格的,我们的数值实验表明,它们具有能够在不同网格上进行训练和泛化的所需特性。这是因为网络学习了无限维函数空间之间的映射,这些映射可以在不同的离散化层次上共享逼近。此外,另一个优势是数据可以通过使用 Nyström 近似加入非结构化网格。我们展示了我们的方法能够与数值分析社区开发的其他无网格方法相竞争,并优于当前状态下神经网络在大网格上的表现,后者依赖于网格。
数值分析社区开发的方法相对较不灵活,因为它们依赖于散度形式的椭圆 PDE 的变分结构。我们新的无网格方法有很多应用,它有潜力成为一种更快的求解器,能够从物理空间中稀疏的观测中学习。它是唯一一种能够在半监督场景下工作的算法,即当我们只在网格的某些部分有测量值时。它也是唯一一种能够在不同几何形状之间进行转换的算法。例如,当计算多个不同翼型的流体动力学时,我们可以构建不同的图并一起训练。在不规则网格上学习时,查询新的位置,我们的方法不需要任何插值,从而避免了插值误差。
缺点: 图核网络的运行时间和存储规模与边的数量成比例, E = O ( K 2 ) E = O(K^2) E=O(K2)。而其他依赖网格的方法如 PCA+NN 和 RBM 只需要 O ( K ) O(K) O(K)。这是不可避免的,因为为了学习连续函数或核,我们需要捕捉每两个节点之间的成对信息,其复杂度为 O ( K 2 ) O(K^2) O(K2)。当离散化固定时,只需要捕捉逐点信息,其复杂度为 O ( K ) O(K) O(K)。因此,当网格较大时,训练和评估整个网格的代价较高。而另一方面,使用子采样来减轻成本会导致数据中的信息丢失,并产生错误,这使得我们的方法不如 PCA+NN 和 RBM 具有竞争力。
未来工作: 为了应对上述问题,我们提出了将诸如多重网格和快速多极方法(Gholami et al., 2016)的想法与我们的方法结合,以减少复杂性。特别是,多重网格方法将构建与不同分辨率对应的多图,这样,在每个图中,节点只连接到它们最近的邻居。边的数量将缩放为 O ( K ) O(K) O(K),而不是 O ( K 2 ) O(K^2) O(K2)。来自 Nyström 近似和局部截断的误差可以避免。另一个方向是将该框架扩展到时间相关的 PDEs。由于图核网络本质上是一个迭代求解器,每个时间步 t t t 都对应于 PDEs 的时间步长,因此自然可以将其框架化为一个 RNN。