参数偏微分方程的傅里叶神经算子
论文链接:https://arxiv.org/abs/2010.08895
项目链接:https://github.com/neuraloperator/neuraloperator
作者博客:https://zongyi-li.github.io/blog/2020/fourier-pde/
- 参数偏微分方程的傅里叶神经算子
- ABSTRACT
- 1 INTRODUCTION
- 2 LEARNING OPERATORS
- 3 NEURAL OPERATOR
- 4 FOURIER NEURAL OPERATOR
- 5 NUMERICAL EXPERIMENTS
- 5.1 伯格斯方程BURGERS’ EQUATION
- 5.2 达西渗流DARCY FLOW
- 5.3 纳维-斯托克斯方程Navier-Stokes equations
- 5.4 Zero-shot超分辨率
- 5.5 贝叶斯逆问题
- 6 DISCUSSION AND CONCLUSION
- A APPENDIX
- A.1 符号表
- A.2 光谱分析
- A.3 数据生成
- A.4 伯格斯方程和达西渗流结果
- A.5 贝叶斯逆问题
ABSTRACT
神经网络的经典发展主要集中在有限维欧几里得空间之间的学习映射。最近,这被推广到学习函数空间之间映射的神经算子。对于偏微分方程,神经算子直接学习从任意函数参数依赖到解的映射。因此,他们学习了整个偏微分方程家族,而不像经典方法只解一个方程实例。在这项工作中,我们通过直接在傅里叶空间中参数化积分核来制定一个新的神经算子,允许一个表达和高效的架构。我们对Burgers方程、Darcy流和Navier-Stokes方程进行了实验。傅里叶神经算子是第一个成功模拟zero-shot超分辨率湍流的基于ML的方法。与传统的PDE求解器相比,它的速度快了三个数量级。此外,在固定分辨率下,与以往基于学习的求解器相比,它具有更高的精度。
1 INTRODUCTION
科学和工程中的许多问题涉及到对某些参数的不同值反复求解复杂的偏微分方程系统。例子出现在分子动力学、微观力学和湍流中。通常这样的系统需要精细的离散化,以便捕捉被建模的现象。因此,传统的数值求解器速度很慢,有时效率低下。例如,当设计诸如翼型之类的材料时,需要解决相关的逆问题,其中需要对正演模型进行数千次评估。一种快速的方法可以使这类问题变得可行。
传统解决方案vs.数据驱动的方法。传统的求解方法如有限元法(FEM)和有限差分法(FDM)是通过对空间进行离散来求解方程的。因此,它们对分辨率进行了权衡:粗网格速度快,但精度较低;精细网格是精确的,但速度很慢。如上所述,复杂的PDE系统通常需要非常精细的离散化,因此对于传统的求解器来说非常具有挑战性和耗时。另一方面,数据驱动的方法可以直接从数据中学习方程族的轨迹。因此,基于学习的方法可以比传统的求解器快几个数量级。
机器学习方法可能是科学学科革命的关键,因为它提供了近似或增强传统解决方案的快速解决方案(Raissi等人,2019;Jiang等,2020;Greenfeld等人,2019;Kochkov et al, 2021)。然而,经典的神经网络在有限维空间之间映射,因此只能学习与特定离散化相关的解决方案。这通常是实际应用的限制,因此需要发展网格不变神经网络。我们首先概述了两种主流的基于神经网络的偏微分方程求解方法——有限维算子和神经有限元法。
有限维空间算子。这些方法将解算子参数化为有限维欧几里得空间之间的深度卷积神经网络(Guo et al ., 2016);朱& Zabaras (2018);Adler & Oktem (2017);Bhatnagar等人(2019);Khoo等人(2017)。根据定义,这些方法依赖于网格,并且需要针对不同的分辨率和离散化进行修改和调整,以实现一致的误差(如果可能的话)。此外,这些方法仅限于训练数据的离散化大小和几何形状,因此不可能在域中的新点查询解。相反,对于我们的方法,我们展示了误差对网格分辨率的不变性,以及在网格之间传递解决方案的能力。
Neural-FEM。第二种方法直接将解函数参数化为神经网络(E & Yu, 2018;Raissi等人,2019;Bar & Sochen, 2019;Smith等人,2020;Pan & Duraisamy, 2020)。此方法旨在为PDE的一个特定实例建模,而不是为解决方案算子建模。它是网格无关的和准确的,但对于任何给定的新的函数参数/系数实例,它需要训练一个新的神经网络。该方法与有限元等经典方法非常相似,用神经网络空间代替局部基函数的有限集的线性跨度。神经有限元法与经典方法存在着相同的计算问题:需要对每一个新实例求解优化问题。此外,该方法仅限于底层PDE已知的设置。
神经算子。最近,一项新的研究提出了用神经网络学习无网格、无限维算子(Lu et al, 2019;Bhattacharya et al ., 2020;尼尔森&斯图尔特,2020;李等,2020b;a;Patel et al, 2021)。神经算子通过产生一组可用于不同离散化的网络参数,弥补了上述有限维算子方法的网格依赖性质。它具有在网格之间传递解的能力。此外,神经算子只需要训练一次。获得一个新的参数实例的解只需要网络的前向传递,从而减轻了神经有限元方法中产生的主要计算问题。最后,神经算子不需要底层PDE的知识,只需要数据。到目前为止,由于计算积分算子的成本,神经算子还没有产生有效的数值算法,可以在有限维环境中与卷积或循环神经网络的成功相媲美。通过快速傅里叶变换,我们的工作减轻了这个问题。
傅里叶变换。傅里叶变换经常用于求解微分方程的频谱方法,因为微分相当于傅里叶域中的乘法。傅里叶变换在深度学习的发展中也发挥了重要作用。理论上,它们出现在普遍近似定理的证明中(Hornik等人,1989),经验上,它们已被用于加速卷积神经网络(Mathieu等人,2013)。涉及傅里叶变换或使用正弦激活函数的神经网络架构也被提出和研究(Bengio等人,2007;Mingo et al, 2004;Sitzmann et al, 2020)。最近,一些谱方法已经扩展到神经网络(Fan et al .,20119a;b;Kashinath et al, 2020)。在这些工作的基础上,我们提出了一个直接在傅里叶空间中定义的神经算子体系结构,具有准线性时间复杂度和最先进的近似能力。
我们的贡献。我们引入傅里叶神经算子,这是一种新的深度学习架构,能够学习无穷维函数空间之间的映射;积分算子被限制为卷积,并通过傅里叶域中的线性变换实例化。
- 傅里叶神经算子是第一个学习湍流状态下Navier-Stokes方程家族的分辨率不变解算子的工作,其中以前基于图的神经算子不收敛。
- 通过构造,该方法共享相同的学习网络参数,而不管在输入和输出空间上使用的离散化。它可以实现零拍摄超分辨率:在较低分辨率上训练,直接在较高分辨率上评估,如图1所示。
- 所提出的方法始终优于所有现有的深度学习方法,即使将分辨率固定为64×64。它在Burgers方程上的错误率降低了30%,在Darcy Flow上降低了60%,在Navier Stokes(黏度ν = 1e−4的湍流状态)上降低了30%。在学习整个时间序列的映射时,当黏度ν = 1e−3时,该方法的误差< 1%,当黏度ν = 1e−4时,该方法的误差为8%。
- 在256 × 256网格上,傅里叶神经算子的推理时间仅为0.05 s,而用于求解Navier-Stokes的伪谱方法的推理时间为2.20 s。尽管该方法具有巨大的速度优势,但在下游应用(如解决Bayesian逆问题)中使用时,其精度不会下降,如图6所示。
我们观察到,所提出的框架可以近似于高度非线性、高频模式和缓慢能量衰减的偏微分方程中的复杂算子上升。神经算子的力量来自于线性、全局积分算子(通过傅里叶变换)和非线性、局部激活函数的结合。类似于标准神经网络通过结合线性乘法和非线性激活来近似高度非线性函数的方式,所提出的神经算子可以近似高度非线性算子。
2 LEARNING OPERATORS
我们的方法从观察到的输入输出对的有限集合中学习两个无限维空间之间的映射。设 D ⊂ R d D\subset\mathbb{R}^d D⊂Rd 为有界开集,且 A = A ( D ; R d a ) \mathcal{A}=\mathcal{A}(D;\mathbb{R}^{d_a}) A=A(D;Rda), U = U ( D ; R d u ) U=\mathcal{U}(D;\mathbb{R}^{d_u}) U=U(D;Rdu)是分别取 R d a \mathbb{R}^{d_a} Rda和 R d u \mathbb{R}^{d_u} Rdu值的函数的可分离巴拿赫空间。此外,让 G † : A → U G^\dagger:\mathcal{A}\to\mathcal{U} G†:A→U是一个(典型的)非线性映射。我们研究作为参数偏微分方程的解算子出现的映射 G † G^\dagger G†-参见第5节的例子。假设我们有观测值 { a j , u j } j = 1 N \left\{a_j,u_j\right\}_{j=1}^N {aj,uj}j=1N,其中 a j ∼ μ a_j\sim\mu aj∼μ是A上支持的概率测度 μ \mu μ的 A \mathcal{A} A序列, u j = G † ( a j ) u_j=G^\dagger(a_j) uj=G†(aj)可能被噪声破坏。我们的目标是通过构造一个参数映射来建立 G † G^\dagger G†的近似值:
G : A × Θ → U or equivalently, G θ : A → U , θ ∈ Θ (1) G:\mathcal{A}\times\Theta\to\mathcal{U}\quad\text{or equivalently,}\quad G_\theta:\mathcal{A}\to\mathcal{U},\quad\theta\in\Theta \tag{1} G:A×Θ→Uor equivalently,Gθ:A→U,θ∈Θ(1)
对于某个有限维参数空间 Θ \Theta Θ,通过选择 θ † ∈ Θ \theta^\dagger\in\Theta θ†∈Θ,使得 G ( ⋅ , θ † ) = G θ † ≈ G † G(\cdot,\theta^\dagger)=G_{\theta^\dagger}\approx G^\dagger G(⋅,θ†)=Gθ†≈G†。这是在无限维中学习的自然框架,就像我们可以定义一个代价函数 C : U × U → R C:\mathcal{U}\times\mathcal{U}\to\mathbb{R} C:U×U→R然后求问题的最小值。
min θ ∈ Θ E a ∼ μ [ C ( G ( a , θ ) , G † ( a ) ) ] (2) \min_{\theta\in\Theta}\mathbb{E}_{a\sim\mu}[C(G(a,\theta),G^\dagger(a))] \tag{2} θ∈ΘminEa∼μ[C(G(a,θ),G†(a))](2)
这与经典的有限维环境直接相似(Vapnik, 1998)。在无限维的情况下,如何证明最小值的存在,仍然是一个具有挑战性的开放性问题。我们将在测试列车设置中通过使用数据驱动的经验近似来解决这个问题,该近似用于确定θ并测试近似的准确性。因为我们在无限维的情况下概念化了我们的方法,所以所有有限维的近似都有一个共同的参数集,这些参数集在无限维中是一致的。表示法表见附录A.3。
可学习算子。近似运算符 G † G^{\dagger} G†是一个不同的任务,通常比为参数 a ∈ A a\in\mathcal{A} a∈A的单个实例找到PDE的解 u ∈ U u\in\mathcal{U} u∈U更具挑战性。大多数现有的方法,从经典的有限元素,有限差分,以及有限体积到现代机器学习方法,如物理信息神经网络(PINN)(Raissi et al.,2019)针对后者,因此计算成本可能很高。这对于需要为许多不同的参数实例提供PDE解决方案的应用程序来说是不切实际的。另一方面,我们的方法直接接近运算符,因此更便宜、更快,与传统的求解器相比,可以节省大量的计算量。有关贝叶斯逆问题的示例应用程序,请参见第5.5节。
离散化。由于我们的数据 a j a_j aj和 u j u_j uj通常都是函数,因此我们假设只能访问逐点求值。令 D j = { x 1 , … , x n } ⊂ D {D}_{j}= \{ x_{1}, \ldots , x_{n}\} \subset D Dj={x1,…,xn}⊂D是域 D D D的 n n n点离散化,并假设我们有观测值 a j ∣ D j ∈ R n × d a , u j ∣ D j ∈ R n × d v a_{j}| _{D_{j}}\in \mathbb{R} ^{n\times d_{a}}, u_{j}| _{D_{j}}\in \mathbb{R} ^{n\times d_{v}} aj∣Dj∈Rn×da,uj∣Dj∈Rn×dv,对于以 j j j为索引的输入输出对的有限集合。为了离散化不变量,神经算子可以对任意 x ∈ D x\in D x∈D产生答案 u ( x ) u(x) u(x),可能是 x ∉ D j x\notin D_j x∈/Dj。这样的性质是非常可取的,因为它允许在不同的网格几何和离散之间传递解。
3 NEURAL OPERATOR
在(Li et al ., 2020b)中提出的神经算子被表述为迭代架构 v 0 ↦ v_0\mapsto v0↦ v 1 ↦ … ↦ v T v_1\mapsto\ldots\mapsto v_T v1↦…↦vT 其中 v j v_j vj , j = 0 , 1 , … , T − 1 j=0,1,\ldots,T-1 j=0,1,…,T−1是一个取 R v d \mathbb{R}^d_v Rvd值的函数序列。如图2 (a)所示,输入 a ∈ A a\in\mathcal{A} a∈A首先通过局部变换 P P P提升到更高的维度表示 v 0 ( x ) = P ( a ( x ) ) v_0(x)=P(a(x)) v0(x)=P(a(x)),该变换通常由浅全连接神经网络参数化。然后我们应用几个迭代的更新 v t ↦ v t + 1 v_{t}\mapsto v_{t+ 1} vt↦vt+1(定义见下)。输出 u ( x ) = Q ( v T ( x ) ) u( x) = Q( v_{T}( x) ) u(x)=Q(vT(x))是 v T v_{T} vT通过局部变换 Q : R d v → R d u Q: \mathbb{R} ^{d_{v}}\to \mathbb{R} ^{d_{u}} Q:Rdv→Rdu的投影。在每次迭代中,更新 v t ↦ v t + 1 v_t\mapsto v_{t+1} vt↦vt+1被定义为非局部积分算子 K \mathcal{K} K和局部非线性激活函数 σ \sigma σ的组合。
定义1(迭代更新) 定义对表示 v t ↦ v t + 1 v_t\mapsto v_{t+1} vt↦vt+1的更新,
v t + 1 ( x ) : = σ ( W v t ( x ) + ( K ( a ; ϕ ) v t ) ( x ) ) , ∀ x ∈ D (2) v_{t+1}(x):=\sigma\Big(Wv_t(x)+\Big(K(a;\phi)v_t\Big)(x)\Big),\quad\forall x\in D \tag{2} vt+1(x):=σ(Wvt(x)+(K(a;ϕ)vt)(x)),∀x∈D(2)
其中 K : A × Θ K → L ( U ( D ; R d v ) , U ( D ; R d v ) ) K: \mathcal{A} \times \Theta _{\mathcal{K} }\to \mathcal{L} ( U( D; \mathbb{R} ^{d_v}) , \mathcal{U} ( D; \mathbb{R} ^{d_v}) ) K:A×ΘK→L(U(D;Rdv),U(D;Rdv)) 映射到 U ( D ; R d v ) \mathcal{U}(D;\mathbb{R}^{d_v}) U(D;Rdv),参数化为 ϕ ∈ Θ K , W : R d v → R d v \phi\in\Theta_{\mathcal{K}},W:\mathbb{R}^{d_v}\to\mathbb{R}^{d_v} ϕ∈ΘK,W:Rdv→Rdv是一个线性变换, σ : R → R \sigma:\mathbb{R}\to\mathbb{R} σ:R→R是一个非线性激活函数,其作用是按组件定义的。
我们选择 K ( a ; ϕ ) \mathcal{K}(a;\phi) K(a;ϕ)是由神经网络参数化的核积分变换。
定义2(核积分算子 K \mathcal{K} K) 定义映射到公式(2)中的核积分算子
( K ( a ; ϕ ) v t ) ( x ) : = ∫ D κ ( x , y , a ( x ) , a ( y ) ; ϕ ) v t ( y ) d y , ∀ x ∈ D (3) \big(\mathcal{K}(a;\phi)v_t\big)(x):=\int_D\kappa\big(x,y,a(x),a(y);\phi\big)v_t(y)\mathrm{d}y,\quad\forall x\in D \tag{3} (K(a;ϕ)vt)(x):=∫Dκ(x,y,a(x),a(y);ϕ)vt(y)dy,∀x∈D(3)
其中 κ ϕ : R 2 ( d + d a ) → R d v × d v \kappa_{\phi}:\mathbb{R}^{2(d+d_{a})}\to\mathbb{R}^{d_{v}\times d_{v}} κϕ:R2(d+da)→Rdv×dv是 ϕ ∈ Θ K \phi\in\Theta_{\mathcal{K}} ϕ∈ΘK参数化的神经网络。
这里 κ ϕ \kappa_{\phi} κϕ扮演核函数的角色,我们从数据中学习。定义1和定义2共同构成了Li等人(2020b)首次提出的将神经网络推广到无限维空间的方法。注意,即使积分算子是线性的,神经算子也可以通过将线性积分算子与非线性激活函数组合来学习高度非线性的算子,类似于标准神经网络。
如果我们去除对函数 a a a的依赖,并施加 κ ϕ ( x , y ) = κ ϕ ( x − y ) \kappa_\phi(x,y)=\kappa_\phi(x-y) κϕ(x,y)=κϕ(x−y),我们得到公式(3)是一个卷积算子,从基本解的角度来看,这是一个自然的选择。在下一节中,我们通过直接在傅里叶空间中参数化 κ ϕ \kappa_\phi κϕ并使用快速傅里叶变换(FFT)来有效地计算公式(3)来利用这一事实。这导致了一个快速架构,可以获得PDE问题的最新结果。
4 FOURIER NEURAL OPERATOR
我们建议用傅里叶空间中定义的卷积算子代替公式(3)中的核积分算子。设 F \mathcal{F} F表示函数 F F F的傅里叶变换 f : D → R d v {f}:D\to\mathbb{R}^{d_{v}} f:D→Rdv和 F − 1 {\mathcal{F}^{-1}} F−1是它的逆。
( F f ) j ( k ) = ∫ D f j ( x ) e − 2 i π ⟨ x , k ⟩ d x , ( F − 1 f ) j ( x ) = ∫ D f j ( k ) e 2 i π ⟨ x , k ⟩ d k (\mathcal{F}f)_j(k)=\int_Df_j(x)e^{-2i\pi\langle x,k\rangle}\mathrm{d}x,\quad(\mathcal{F}^{-1}f)_j(x)=\int_Df_j(k)e^{2i\pi\langle x,k\rangle}\mathrm{d}k (Ff)j(k)=∫Dfj(x)e−2iπ⟨x,k⟩dx,(F−1f)j(x)=∫Dfj(k)e2iπ⟨x,k⟩dk
对于 j = 1 , … , d v j=1,\ldots,d_v j=1,…,dv ,其中 i = − 1 i=\sqrt-1 i=−1是虚单位。令在公式(3)中的 κ ϕ ( x , y , a ( x ) , a ( y ) ) = κ ϕ ( x − y ) \kappa_\phi(x,y,a(x),a(y))=\kappa_\phi(x-y) κϕ(x,y,a(x),a(y))=κϕ(x−y),并应用卷积定理得到
( K ( a ; ϕ ) v t ) ( x ) = F − 1 ( F ( κ ϕ ) ⋅ F ( v t ) ) ( x ) , ∀ x ∈ D . \big(\mathcal{K}(a;\phi)v_t\big)(x)=\mathcal{F}^{-1}\big(\mathcal{F}(\kappa_\phi)\cdot\mathcal{F}(v_t)\big)(x),\quad\forall x\in D. (K(a;ϕ)vt)(x)=F−1(F(κϕ)⋅F(vt))(x),∀x∈D.
因此,我们建议在傅里叶空间中直接参数化 κ ϕ \kappa_\phi κϕ。
**定义3(傅里叶积分算子 ( K (\mathcal{K} (K) ** 定义傅里叶积分算子:
( K ( ϕ ) v t ) ( x ) = F − 1 ( R ϕ ⋅ ( F v t ) ) ( x ) ∀ x ∈ D (4) \big(\mathcal{K}(\phi)v_t\big)(x)=\mathcal{F}^{-1}\big(R_\phi\cdot(\mathcal{F}v_t)\big)(x)\quad\forall x\in D \tag{4} (K(ϕ)vt)(x)=F−1(Rϕ⋅(Fvt))(x)∀x∈D(4)
其中, R φ R_φ Rφ是周期函数 κ : D ˉ → R d v × d v \kappa:\bar{D}\to\mathbb{R}^{d_{v}\times d_{v}} κ:Dˉ→Rdv×dv由 ϕ ∈ Θ K \phi\in\Theta_{\mathcal{K}} ϕ∈ΘK参数化。图2 (b)给出了一个说明。
对于频率模式 k ∈ D k\in D k∈D,我们有 ( F v t ) ( k ) ∈ C v d (\mathcal{F}v_t)(k)\in\mathbb{C}^d_v (Fvt)(k)∈Cvd和 R ϕ ( k ) ∈ C d v × d v R_\phi(k)\in\mathbb{C}^{d_v\times d_v} Rϕ(k)∈Cdv×dv。注意,因为我们假设 κ \kappa κ是周期性的,它允许傅里叶级数展开,所以我们可以处理离散模式 k ∈ Z d k\in\mathbb{Z}^d k∈Zd。我们选择一个有限维的参数化通过截断在最大模态数的傅里叶级数 k max = ∣ Z k max ∣ = ∣ { k ∈ Z d : ∣ k j ∣ ≤ k max , j k_{\max }= | Z_{k_{\max }}| = | \{ k\in \mathbb{Z} ^d: | k_j| \leq k_{\max , j} kmax=∣Zkmax∣=∣{k∈Zd:∣kj∣≤kmax,j, for $j= 1, \ldots , d} | $ 。因此,我们直接参数化 R ϕ R_{\phi } Rϕ为复值 ( k m a x × d v × d v ) ( k_{\mathrm{max}}\times d_{v}\times d_{v}) (kmax×dv×dv)张量,包括截断的傅立叶模的集合,因此从我们的符号中去掉 ϕ \phi ϕ。由于 κ \kappa κ是实值,我们施加共轭对称性。我们注意到集合Zkmax并不是 v t v_{t} vt的低频模态的标准选择。事实上,低频模态通常是通过在 k ∈ Z d k\in\mathbb{Z}^{d} k∈Zd的1范数上放置一个上界来定义的。我们如上所述选择 Z k max Z_{k_{\max}} Zkmax,因为它允许有效的实现。
离散情况和FFT。假设域 D D D被 n ∈ N n\in\mathbb{N} n∈N个点离散,我们得到 v t ∈ R n × d v v_t\in\mathbb{R}^{n\times d_v} vt∈Rn×dv和 F ( v t ) ∈ C n × d v \mathcal{F}(v_t)\in\mathbb{C}^{n\times d_v} F(vt)∈Cn×dv。由于我们将 v t v_t vt与一个只有 k m a x k_{\mathrm{max}} kmax傅里叶模式的函数进行卷积,我们可以简单地截断较高的模式以获得 F ( v t ) ∈ C k m a x × d v \mathcal{F} ( v_{t}) \in \mathbb{C} ^{k_{\mathrm{max}}\times d_{v}} F(vt)∈Ckmax×dv。乘以权张量 R ∈ C k m a x × d v × d v R\in \mathbb{C} ^{k_{\mathrm{max}}\times d_{v}\times d_{v}} R∈Ckmax×dv×dv就是
( R ⋅ ( F v t ) ) k , l = ∑ j = 1 d v R k , l , j ( F v t ) k , j , k = 1 , … , k max , j = 1 , … , d v . (5) \left(R\cdot(\mathcal{F}v_t)\right)_{k,l}=\sum_{j=1}^{d_v}R_{k,l,j}(\mathcal{F}v_t)_{k,j},\quad k=1,\ldots,k_{\max},\quad j=1,\ldots,d_v. \tag{5} (R⋅(Fvt))k,l=j=1∑dvRk,l,j(Fvt)k,j,k=1,…,kmax,j=1,…,dv.(5)
当离散化均匀且分辨率为 s 1 × ⋯ × s d = n s_1\times\cdots\times s_d=n s1×⋯×sd=n时,可以用快速傅里叶变换代替 F \mathcal{F} F。对于 f ∈ R n × d v , k = ( k 1 , … , k d ) ∈ Z s 1 × ⋯ × Z s d f\in\mathbb{R}^{n\times d_v},k=(k_1,\ldots,k_d)\in\mathbb{Z}_{s_1}\times\cdots\times\mathbb{Z}_{s_d} f∈Rn×dv,k=(k1,…,kd)∈Zs1×⋯×Zsd, x = ( x 1 , … , x d ) ∈ D x=(x_1,\ldots,x_d)\in D x=(x1,…,xd)∈D, FFT F ^ \hat{\mathcal{F}} F^和它的逆 F ^ − 1 \hat{\mathcal{F}}^{-1} F^−1定义为
( F ^ f ) l ( k ) = ∑ x 1 = 0 s 1 − 1 ⋯ ∑ x d = 0 s d − 1 f l ( x 1 , … , x d ) e − 2 i π ∑ j = 1 d x j k j a j , ( F ^ − 1 f ) l ( x ) = ∑ k 1 = 0 s 1 − 1 ⋯ ∑ k d = 0 s d − 1 f l ( k 1 , … , k d ) e 2 i π ∑ j = 1 d x j k j s j \begin{align} (\hat{\mathcal{F}}f)_{l}(k)&=\sum_{x_{1}=0}^{s_{1}-1}\cdots\sum_{x_{d}=0}^{s_{d}-1}f_{l}(x_{1},\ldots,x_{d})e^{-2i\pi\sum_{j=1}^{d}\frac{x_{j}k_{j}}{a_{j}}},\\ (\hat{\mathcal{F}}^{-1}f)_l(x)&=\sum_{k_1=0}^{s_1-1}\cdots\sum_{k_d=0}^{s_d-1}f_l(k_1,\ldots,k_d)e^{2i\pi\sum_{j=1}^d\frac{x_jk_j}{s_j}} \end{align} (F^f)l(k)(F^−1f)l(x)=x1=0∑s1−1⋯xd=0∑sd−1fl(x1,…,xd)e−2iπ∑j=1dajxjkj,=k1=0∑s1−1⋯kd=0∑sd−1fl(k1,…,kd)e2iπ∑j=1dsjxjkj
对于 l = 1 , . . . , d v l = 1,...,d_v l=1,...,dv$。在这种情况下,截断模式集变为
Z k m a x = { ( k 1 , … , k d ) ∈ Z s 1 × ⋯ × Z s d ∣ k j ≤ k m a x , j o r s j − k j ≤ k m a x , j , f o r j = 1 , … , d } . Z_{k_{\mathrm{max}}}=\{(k_1,\ldots,k_d)\in\mathbb{Z}_{s_1}\times\cdots\times\mathbb{Z}_{s_d}\mid k_j\leq k_{\mathrm{max},j}\mathrm{~or~}s_j-k_j\leq k_{\mathrm{max},j},\mathrm{~for~}j=1,\ldots,d\}. Zkmax={(k1,…,kd)∈Zs1×⋯×Zsd∣kj≤kmax,j or sj−kj≤kmax,j, for j=1,…,d}.
当实现时, R R R被视为 ( s 1 × ⋯ × s d × d v × d v ) (s_1\times\cdots\times s_d\times d_v\times d_v) (s1×⋯×sd×dv×dv)张量, Z k m a x Z_{k_\mathrm{max}} Zkmax的上述定义对应于 R R R的“corners”,这允许通过矩阵向量乘法直接并行实现公式(5)。在实践中,我们发现选择 k m a x , j = 12 k_{\mathrm{max},j}=12 kmax,j=12,每个通道产生 k m a x = 1 2 d k_\mathrm{max}=12^d kmax=12d个参数,足以满足我们考虑的所有任务。
R R R的参数化。通常, R R R可以定义为依赖于 ( F a ) (\mathcal{F}a) (Fa)来并行公式(3)。实际上,我们可以定义 R ϕ : Z d × R d v → R d v × d v R_\phi:\mathbb{Z}^d\times\mathbb{R}^{d_v}\to\mathbb{R}^{d_v\times d_v} Rϕ:Zd×Rdv→Rdv×dv 作为一个参数函数,它映射 ( k , ( F a ) ( k ) ) (k,(\mathcal{F}a)(k)) (k,(Fa)(k))到适当的傅里叶模态的值。我们对 R ϕ R_\phi Rϕ的线性和神经网络参数化进行了实验。我们发现线性参数化与之前描述的直接参数化具有相似的性能,而神经网络的性能较差。这可能是由于空间 Z d \mathbb{Z}^d Zd的离散结构。我们在这项工作中的实验集中在上述直接参数化上。
离散化的不变性。傅里叶层是离散不变的因为它们可以从任意离散的函数中学习并求值。由于参数是直接在傅里叶空间中学习的,因此在物理空间中解析函数就等于在 R d {R}^d Rd上处处定义良好的 e 2 π i ⟨ x , k ⟩ e^{2\pi i\langle x,k\rangle} e2πi⟨x,k⟩的基础上进行投影。这使我们能够实现如5.4节所示的Zero-shot超分辨率。此外,我们的体系结构在输入和输出的任何分辨率上都有一致的错误。另一方面,请注意,在图3中,我们比较的标准CNN方法的误差随着分辨率的增加而增加。
准线性复杂度。权重张量 R R R包含 k m a x < n k_\mathrm{max}<n kmax<n个模式,因此内部乘法的复杂度为 O ( k max ) O(k_{\max}) O(kmax)。因此,大部分的计算成本在于计算傅里叶变换 F ( v t ) \mathcal{F}(v_t) F(vt)及其逆。一般傅里叶变换的复杂度为 O ( n 2 ) O(n^2) O(n2),然而,由于我们截断了序列,复杂度实际上是 O ( n k m a x ) O(nk_\mathrm{max}) O(nkmax),而FFT的复杂度为 O ( n log n ) O(n\log n) O(nlogn)。通常,我们发现使用FFT非常有效。然而,需要统一的离散化。
5 NUMERICAL EXPERIMENTS
在本节中,我们将提出的傅里叶神经算子与多个有限维架构以及基于算子的近似方法在一维Burgers方程、二维Darcy流问题和二维Navier-Stokes方程上进行比较。数据生成过程分别在附录A.3.1、A.3.2和A.3.3中讨论。我们不与传统的求解器(FEM/FDM)或神经-FEM类型的方法进行比较,因为我们的目标是产生可用于下游应用的有效算子近似。我们将在第5.5节中演示贝叶斯逆问题的一个这样的应用。
我们通过将(2)和(4)中指定的四个傅里叶积分算子层与ReLU激活以及批处理归一化叠加来构建傅里叶神经算子。除非另有说明,否则我们使用N = 1000个训练实例和200个测试实例。我们使用Adam优化器训练500个epoch,初始学习率为0.001,每100个epoch减半。我们设置一维问题 k max , j = 16 , d v = 64 k_{\max,j}=16,d_{v}=64 kmax,j=16,dv=64;对于二维问题, k max , j = 12 , d v = 32 k_{\max,j}=12,d_{v}=32 kmax,j=12,dv=32。低分辨率数据从高分辨率下采样。所有的计算都在一个具有16GB内存的Nvidia V100 GPU上进行。
分辨率评估。传统的PDE求解方法如FEM和FDM近似于单一函数,因此它们对连续体的误差随着分辨率的增加而减小。另一方面,只要解出所有相关信息,算子近似与数据离散的方式无关。分辨率不变运算符在不同分辨率下具有一致的错误率,如图3所示。此外,分辨率不变运算符可以实现Zero-shot超分辨率,如第5.4节所示。
时间无关问题的基准(Burgers和Darcy):
NN:一个简单的逐点前馈神经网络。
RBM:经典的约简基法(使用POD基) (DeVore, 2014)。
FCN:基于全卷积网络的最先进的神经网络架构(Zhu & Zabaras, 2018)。
PCANN:一种算子方法,使用PCA作为输入和输出数据的自编码器,并用神经网络对潜在空间进行插值(Bhattacharya等人,2020)。
GNO:原始图神经算子(Li et al ., 2020b)。
MGNO:多极图神经算子(Li et al ., 2020a)。
LNO:一种基于核低秩分解的神经算子方法 κ ( x , y ) : = ∑ i = 1 r ϕ j ( x ) ψ j ( y ) \kappa(x,y):=\sum_{i=1}^r\phi_j(x)\psi_j(y) κ(x,y):=∑i=1rϕj(x)ψj(y),类似于(Lu et al, 2019)中提出的未堆叠DeepONet。
FNO:新提出的傅里叶神经算子。
时间相关问题的基准(Navier-Stokes):
ResNet:带有残差连接的18层二维卷积(He et al, 2016)。
U-Net:图像到图像回归任务的流行选择,由四个块组成,具有二维卷积和反卷积(Ronneberger等人,2015)。
TF-Net:基于空间和时间卷积的组合,设计用于学习湍流的网络(Wang et al ., 2020)。
FNO-2d:具有时域RNN结构的二维傅里叶神经算子。
FNO-3d:在时空中直接卷积的三维傅立叶神经算子。
5.1 伯格斯方程BURGERS’ EQUATION
一维Burgers方程是一种非线性偏微分方程,具有多种应用,包括模拟粘性流体的一维流动。它的形式是:
∂ t u ( x , t ) + ∂ x ( u 2 ( x , t ) / 2 ) = ν ∂ x x u ( x , t ) , x ∈ ( 0 , 1 ) , t ∈ ( 0 , 1 ] u ( x , 0 ) = u 0 ( x ) , x ∈ ( 0 , 1 ) (6) \partial_{t}u(x,t)+\partial_{x}(u^{2}(x,t)/2)=\nu\partial_{xx}u(x,t),\quad x\in(0,1),t\in(0,1] \\ u(x,0)=u_{0}(x),\quad x\in(0,1)\tag{6} ∂tu(x,t)+∂x(u2(x,t)/2)=ν∂xxu(x,t),x∈(0,1),t∈(0,1]u(x,0)=u0(x),x∈(0,1)(6)
具有周期边界条件,其中 u 0 ∈ L p e r 2 ( ( 0 , 1 ) ; R ) u_0\in L_{\mathrm{per}}^2( ( 0, 1) ; \mathbb{R} ) u0∈Lper2((0,1);R)为初始条件, ν ∈ R + \nu\in\mathbb{R}_+ ν∈R+ 为粘度系数。我们的目标是学习将初始条件映射到时刻1的解的算子, G † : L p e r 2 ( ( 0 , 1 ) ; R ) → H p e r r ( ( 0 , 1 ) ; R ) G^\dagger : L_{\mathrm{per}}^2( ( 0, 1) ; \mathbb{R} ) \to H_{\mathrm{per}}^r( ( 0, 1) ; \mathbb{R} ) G†:Lper2((0,1);R)→Hperr((0,1);R)由 u 0 ↦ u ( ⋅ , 1 ) u_0\mapsto u(\cdot,1) u0↦u(⋅,1)定义,对于任意 r > 0 r>0 r>0。
我们的实验结果见图3 (a)和表3(附录A.3.1)。与任何基准测试相比,我们提出的方法获得了最低的相对误差。此外,误差随分辨率不变,而基于卷积神经网络的方法(FCN)的误差随分辨率增大。与在物理空间中使用Nystrom采样的其他神经算子方法(如GNO和MGNO)相比,傅里叶神经算子更精确,计算效率更高。
5.2 达西渗流DARCY FLOW
考虑单元上二维达西流方程的稳态,即二阶线性椭圆偏微分方程:
− ∇ ⋅ ( a ( x ) ∇ u ( x ) ) = f ( x ) x ∈ ( 0 , 1 ) 2 u ( x ) = 0 x ∈ ∂ ( 0 , 1 ) 2 (7) -\nabla\cdot(a(x)\nabla u(x))=f(x)\quad x\in(0,1)^{2}\\u(x)=0\quad x\in\partial(0,1)^{2}\tag{7} −∇⋅(a(x)∇u(x))=f(x)x∈(0,1)2u(x)=0x∈∂(0,1)2(7)
具有Dirichlet边界,其中 a ∈ L ∞ ( ( 0 , 1 ) 2 ; R + ) a\in L^\infty((0,1)^2;\mathbb{R}_+) a∈L∞((0,1)2;R+)为扩散系数, f ∈ f\in f∈ L 2 ( ( 0 , 1 ) 2 ; R ) L^2((0,1)^2;\mathbb{R}) L2((0,1)2;R)为强制函数。这种PDE有许多应用,包括模拟地下流动的压力,线性弹性材料的变形,以及导电材料的电势。我们感兴趣的是学习将扩散系数映射到解 G † : L ∞ ( ( 0 , 1 ) 2 ; R + ) → H 0 1 ( ( 0 , 1 ) 2 ; R + ) G^\dagger:L^\infty((0,1)^2;\mathbb{R}_+)\to H_0^1((0,1)^2;\mathbb{R}_+) G†:L∞((0,1)2;R+)→H01((0,1)2;R+)由 a ↦ u a\mapsto u a↦u定义。注意,虽然PDE是线性的,但算子 G † G^\dagger G†不是。
我们的实验结果见图3 (b)和表4(附录A.3.2)。所提出的傅里叶神经算子的相对误差比任何基准都低近一个数量级。我们再次观察到误差相对于分辨率的不变性。
5.3 纳维-斯托克斯方程Navier-Stokes equations
我们考虑在单位环面上以涡量形式存在的粘性不可压缩流体的二维Navier-Stokes方程:
∂ t w ( x , t ) + u ( x , t ) ⋅ ∇ w ( x , t ) = ν Δ w ( x , t ) + f ( x ) , x ∈ ( 0 , 1 ) 2 , t ∈ ( 0 , T ] ∇ ⋅ u ( x , t ) = 0 , x ∈ ( 0 , 1 ) 2 , t ∈ [ 0 , T ] w ( x , 0 ) = w 0 ( x ) , x ∈ ( 0 , 1 ) 2 (8) \begin{aligned} \partial_{t}w(x,t)+u(x,t)\cdot\nabla w(x,t)& =\nu\Delta w(x,t)+f(x), && x\in(0,1)^{2},t\in(0,T] \\ \nabla\cdot u(x,t)& =0, && x\in(0,1)^{2},t\in[0,T] \\ w(x,0)& =w_0(x), && x\in(0,1)^{2} \end{aligned}\tag{8} ∂tw(x,t)+u(x,t)⋅∇w(x,t)∇⋅u(x,t)w(x,0)=νΔw(x,t)+f(x),=0,=w0(x),x∈(0,1)2,t∈(0,T]x∈(0,1)2,t∈[0,T]x∈(0,1)2(8)
其中 u ∈ C ( [ 0 , T ] ; H p e r r ( ( 0 , 1 ) 2 ; R 2 ) ) u\in C([0,T];H_{\mathrm{per}}^{r}((0,1)^{2};\mathbb{R}^{2})) u∈C([0,T];Hperr((0,1)2;R2))对于任意 r > 0 r>0 r>0 为速度场, w = ∇ × u w=\nabla\times u w=∇×u为涡度, w 0 ∈ L p e r 2 ( ( 0 , 1 ) 2 ; R ) w_0\in L_{\mathrm{per}}^2( ( 0, 1) ^2; \mathbb{R} ) w0∈Lper2((0,1)2;R)为初始涡度, ν ∈ R + \nu\in\mathbb{R}_+ ν∈R+为黏度系数, f ∈ f\in f∈ L p e r 2 ( ( 0 , 1 ) 2 ; R ) L_{\mathrm{per}}^2( ( 0, 1) ^2; \mathbb{R} ) Lper2((0,1)2;R)为强制函数。我们感兴趣的是学习将时间为10的涡度映射到之后某个时间的涡度的算子 T > 10 , G † : C ( [ 0 , 10 ] ; H p e r r ( ( 0 , 1 ) 2 ; R ) ) → T> 10, G^\dagger : C( [ 0, 10] ; H_{\mathrm{per}}^r( ( 0, 1) ^2; \mathbb{R} ) ) \to T>10,G†:C([0,10];Hperr((0,1)2;R))→ C ( ( 10 , T ] ; H p e r r ( ( 0 , 1 ) 2 ; R ) ) C( ( 10, T] ; H_{\mathrm{per}}^r( ( 0, 1) ^2; \mathbb{R} ) ) C((10,T];Hperr((0,1)2;R)) 由 w ∣ ( 0 , 1 ) 2 × [ 0 , 10 ] ↦ w ∣ ( 0 , 1 ) 2 × ( 10 , T ] w|_(0,1)^2\times[0,10]\mapsto w|_{(0,1)^2\times(10,T]} w∣(0,1)2×[0,10]↦w∣(0,1)2×(10,T]定义。已知涡度,很容易求出速度。虽然与速度相比,涡度更难建模,但它提供了更多的信息。通过对涡度问题的表述,神经网络模型模拟了伪谱方法。我们实验粘度 ν = 1 e − 3 , 1 e − 4 , 1 e − 5 ν = 1e−3,1e−4,1e−5 ν=1e−3,1e−4,1e−5,随着动态变混沌,最终时间T减小。由于基线方法不是分辨率不变的,我们将训练和测试的分辨率都固定为64 × 64。
由表1可知,当数据充足时,FNO-3D的性能最佳( ν = 1 e − 3 , N = 1000 ν = 1e−3,N = 1000 ν=1e−3,N=1000和$ ν = 1e−4,N = 10000 ) 。对于数据量不足的配置 ( )。对于数据量不足的配置( )。对于数据量不足的配置(ν = 1e−4,N = 1000 和 和 和 ν = 1e−5,N = 1000$),所有方法的误差都>15%,其中FNO-2D达到最低。请注意,我们只给出空间分辨率为64 × 64的结果,因为我们比较的所有基准测试都是针对该分辨率设计的。增加它会降低它们的性能,而FNO会产生相同的错误。
2D和3D卷积。FNO-2D、U-Net、TF-Net和ResNet都在空间域中进行2D卷积,在时间域中循环传播(2D+RNN)。算子将前10个时间步的解映射到下一个时间步(2D函数到2D函数)。另一方面,FNO-3D在时空中进行卷积。它将初始时间步长直接映射到完整的轨迹(3D函数到3D函数)。2D+RNN结构可以将解以固定区间长度∆T的增量传播到任意时间T,而Conv3D结构则固定为区间[0,T],但可以将解转换为任意时间离散化。我们发现,与RNN结构的方法相比,3-d方法更具表现力,更容易训练。
5.4 Zero-shot超分辨率
神经算子是网格不变性的,因此它可以在较低分辨率下进行训练,并在较高分辨率下进行评估,而无需看到任何更高分辨率的数据(Zero-shot超分辨率)。图1显示了一个例子,我们在上面设置的64 × 64 × 20分辨率数据上训练FNO-3D模型,(ν = 1e−4;N = 10000),转换为256 × 256 × 80分辨率,展示时空超分辨率。傅里叶神经算子是基准(FNO-2D, U-Net, TF-Net和ResNet)中唯一可以实现Zero-shot超分辨率的模型。令人惊讶的是,它不仅可以在空间域而且可以在时间域实现超分辨率。
5.5 贝叶斯逆问题
在本实验中,我们使用函数空间马尔可夫链蒙特卡罗(MCMC)方法(Cotter et al, 2013)从给定时间T = 50的稀疏噪声观测值的Navier-Stokes初始涡度的后检分布中提取样本。我们将傅里叶神经算子作为替代模型与用于生成训练测试数据的传统求解器进行比较(两者都在GPU上运行)。我们从后验生成25,000个样本(具有5,000个样本老化期),需要对前向算子进行30,000次评估。
如图6(附录A.5)所示,FNO和传统求解器恢复几乎相同的后置均值,当向前推进时,可以很好地恢复Navier Stokes的后期动态。与之形成鲜明对比的是,FNO计算单个实例的时间为0.05 s,而传统的求解器在优化使用最大可能的内部时间步长(不会导致爆炸)后,需要2.2s。对于使用FNO的MCMC来说,这相当于2.5分钟,而对于传统求解器来说,这超过了18个小时。即使我们考虑到数据生成和训练时间(离线步骤)需要12个小时,使用FNO仍然更快!经过训练后,FNO可以针对不同的初始条件和观测值快速执行多个MCMC运行,而传统的求解器每个实例需要18个小时。此外,由于FNO是可微的,它可以很容易地应用于PDE约束优化问题,而不需要伴随方法。
光谱分析。由于我们参数化 R ϕ R_\phi Rϕ的方式,由公式(4)输出的函数每个通道最多有 k ^ max , j \hat{k}_{\max,j} k^max,j个傅里叶模式。然而,这并不意味着傅里叶神经算子只能近似到 k max , j k_{\max,j} kmax,j模式的函数。实际上,发生在积分算子和最终解码器网络 Q Q Q之间的激活函数恢复了高频模式。作为一个例子,考虑粘度为 ν = 1 \nu=1 ν=1e-3的Navier-Stokes方程的解。将该函数截断为20个傅立叶模式会产生约 2 % 2\% 2%的误差,而我们的傅立叶神经算子学习参数依赖性并产生误差 ≤ 1 % \leq1\% ≤1%的近似值,只有 k m a x , j = 12 k_{\mathrm{max},j}=12 kmax,j=12 个参数化模式。
非周期边界条件。传统的傅里叶方法只适用于周期边界条件。然而,傅里叶神经算子没有这种限制。这是由于线性变换 W W W(偏置项)保持了非周期边界的轨迹。以Darcy Flow和Navier-Stokes的时域具有非周期边界条件为例,傅里叶神经算子仍能以极好的精度学习解算子。
6 DISCUSSION AND CONCLUSION
数据要求。数据驱动的方法依赖于数据的质量和数量。为了学习黏度为ν = 1e−4的Navier-Stokes方程,我们需要生成N = 10000对训练对faj;用数值解算器。然而,对于更具挑战性的PDE,生成一些训练样本可能已经非常昂贵了。未来的发展方向是将神经算子与数值求解相结合,使对数据的要求悬浮起来。
复发性结构。神经算子具有迭代结构,可以自然地表述为循环网络,其中所有层共享相同的参数而不牺牲性能。(我们在实验中没有施加这个限制。)
计算机视觉。算子学习并不局限于偏微分方程。图像自然可以被视为二维域上的实值函数,视频只是添加了一个时间结构。因此,我们的方法是计算机视觉问题的自然选择,其中离散化的不变性至关重要(Chi等人,2020)。
A APPENDIX
A.1 符号表
A.2 光谱分析
Navier Stokes方程数据的谱衰减如图4所示。光谱衰减斜率为 k − 5 / 3 k^{-{5}/{3}} k−5/3,与湍流区的能谱相匹配。我们注意到能谱不随时间衰减。
我们还以截断kmax的形式给出了光谱衰减,如图5所示。我们注意到所有方程(Burgers, Darcy和Navier-Stokes ν≤1e−4)都表现出高频模式。即使我们在傅里叶层截断kmax = 12,傅里叶神经算子也可以恢复高频模式。
A.3 数据生成
A.4 伯格斯方程和达西渗流结果
Burgers方程和Darcy Flow的详细错误率如表3和表4所示。
A.5 贝叶斯逆问题
Navier-Stokes方程的Bayesian逆问题结果如图6所示。可以看出,采用傅里叶神经算子作为替代算子的求解结果与传统求解器的求解结果一样好。