基于LM算法的多尺度神经网络求解椭圆界面问题

Solving Elliptic Interface Problems Using a Multi-Scale Neural Network Based on the LM Algorithm

摘要: 近几年以来，利用深度学习方法解决偏微分方程问题引起了广泛关注和研究。本文介绍了一种基于改进的Levenberg-Marquardt (LM)优化方法的多尺度神经网络，该方法在求解椭圆界面问题的精度方面显示出不俗的潜力。文章通过单个神经网络框架解决界面问题，并选用不同的激活函数进行对比，通过几个具有规则和不规则界面的数值实验，以验证神经网络效果及相关理论。

Abstract: In recent years, the application of deep learning methods to solve partial differential equations has garnered significant attention and research. This paper introduces a multi-scale neural network based on an improved Levenberg-Marquardt (LM) optimization method, which demonstrates remarkable potential in achieving high accuracy for solving elliptic interface problems. The paper addresses boundary problems within a single neural network framework, comparing different activation functions. Numerical experiments involving both regular and irregular boundaries validate the neural network’s performance and related theoretical foundations.

1. 引言

近年来，随着计算机技术的飞速进步和科学计算能力的显著提升，深度学习技术得以不断发展。基于神经网络(Neural Networks)的深度学习(Deep Learning)方法在多学科领域特别是在计算机视觉[1] [2] (Computer Vision)和自然语言处理[3] [4] (Natural Language Processing)等领域取得了令人瞩目的成就。这些成就如雨后春笋般层出不穷，吸引了众多学者深入探索。神经网络因其卓越的拟合能力而广受认可，被视为一种功能强大的通用函数逼近工具。最近，这一强大的工具也被引入到偏微分方程的求解中，成为了一个新兴的研究方向。应用深度学习求解偏微分方程的基本思想是通过训练神经网络而非求解代数方程，近似求解偏微分方程的解函数或解算子。目前，该领域已涌现出众多研究进展，并取得了令人瞩目的成果。然而，在利用神经网络求解带有界面条件的偏微分方程时，仍存在精度不足的问题。这类方程因界面上的不连续性，使得普通的神经网络方法难以准确捕捉到界面上的跳跃性特征，这成为求解过程中的一大难点。因此，探索使用更高精度的神经网络方法来求解带有界面条件的偏微分方程，成为了一个值得深入研究的方向。

2018年，Jiequn H等[5]将深度学习方法用于求解高维偏微分方程，19年Cai W等人[6]运用多尺度深度神经网络求解高维偏微分方程，文[7]说明深度时间神经网络(DTNN)是求解高维偏微分方程的有效深度学习方法。而[8]-[11]分别介绍了求解低维非线性偏微分方程的深度学习初始化迭代方法(Int-Deep)、深度最小二乘法这一基于无监督学习的求解椭圆偏微分方程的数值方法、用于完全非线性偏微分方程的深分支求解器及使用深度神经网络求解椭圆偏微分方程的非梯度方法。2024年，Paola F等人[12]介绍用于加速3D椭圆偏微分方程的有限元求解器中的代数多重网格方法，[13] [14]中提出了单位网络自适应分区(APUNet)这一用于解决偏微分方程的本地化深度学习方法及基于改进欧拉法的非线性偏微分方程神经网络求解器。彭杰，张玉武[15]对基于自适应神经网络的偏微分方程进行研究。

与全连接神经网络和ResNet的结果相比，多尺度融合网络[16] [17]作为一种新的深度学习结构来解决椭圆界面问题能够更好地捕捉到界面上的跳跃性特征，从而提高了准确性。此外，它的数值解可以保持解的连续性，同时保持磁通量跳跃通过不同的界面，从而保持微分方程的物理特性。基于此背景，我们利用多尺度神经网络来求解椭圆界面方程。在本文中，我们使用单个网络来解决椭圆界面问题，从而规避了使用两个或多个网络的复杂性。相较于[16] [17]中所使用的Adam优化算法而言，我们借鉴[18]中的LM优化算法，将多尺度融合网络与LM算法结合，并考虑该情况下不同激活函数对求解椭圆界面问题精度的影响效果，并通过数值实验验证我们所采用方法的有效性与价值性。

本文的内容安排如下。在第二节部分介绍我们的多尺度神经网络及LM优化算法。在第三节基于我们所使用的神经网络对规则和不规则界面的二维问题进行了数值实验，以验证我们的理论结果，随后进一步将其用于一个三维问题，展示了该方法的有效性。最后，第四节给出我们的相关结论及展望。

2. 深度学习方法

在这一节，我们以椭圆偏微分方程为例，介绍使用多尺度神经网络求解界面问题的具体流程。考虑以下具有不连续系数的椭圆偏微分方程：

$\{\begin{array}{cc}-\nabla \cdot \left(a\left(x\right)\nabla u\left(x\right)\right)=f\left(x\right),& x\in \Omega \\ u\left(s^{+}\right)=u\left(s^{-}\right)+p(s),& s\in \Gamma \\ \alpha \left(s^{+}\right)\frac{\partial u\left(s^{+}\right)}{\partial n\left(s\right)}=\alpha \left(s^{-}\right)\frac{\partial u\left(s^{-}\right)}{\partial n\left(s\right)}+q\left(s\right),& s\in \Gamma \\ u\left(s\right)=g\left(s\right),& s\in \partial \Omega \end{array}$ (1)

其中$\Omega ={\Omega }_p\cup {\Omega }_s\cup \Gamma$ ，界面$\Gamma$ 为${\Omega }_p$ 和${\Omega }_s$ 的交界处，是${ℝ}^d(d\ge 2$ 为给定的正整数)上的有界凸域，$g\left(s\right)$ 和$p\left(s\right)$ 及$q\left(s\right)$ 分别是边界$\partial \Omega$ 和界面$\Gamma$ 上的给定函数，设系数$\alpha \left(x\right)$ 为跨越界面$\Gamma$ 的分段常数函数，

$\alpha \left(x\right)=\{\begin{array}{ll}{\alpha }_1,\hfill & x\in {\Omega }_p\hfill \\ {\alpha }_2,\hfill & x\in {\Omega }_s\hfill \end{array}$ (2)

在给定的边界条件和界面条件下近似偏微分方程(1)的解。对于具有界面条件的椭圆问题，直接应用前面描述的深度学习方法无法得到准确的数值解。为了克服这个精度问题，学者们提出了几种技术，例如域分解法、扩展和投影法。在这里，我们没有使用域分解方法，而是采用了扩展和投影方法，使用基于这种扩展和投影方法的一个网络可以比使用多个网络获得更好的结果，同时在必要时它可以保持底层量的连续性。给定一个神经网络，为了保证准确性，根据扩展投影法，我们首先用指示函数$z\left(x\right)$ 将之扩展到更高的空间维度。也就是说，我们尝试在$\tilde{\Omega }\equiv \Omega \times ℝ\in {ℝ}^{d+1}$ 上定义$\tilde{u}\left(x,z\right)$ 满足$u\left(x\right)=\tilde{u}\left(x,z\left(x\right)\right)$ 。其中，指示函数$z\left(x\right)$ 用于区分不同的子域${\Omega }_p$ 和${\Omega }_s$ 。对于界面上的情况，我们可以将指标函数$z\left(x\right)$ 视为以下分段常数

$z\left(x\right)=\{\begin{array}{ll}-1,\hfill & x\in {\Omega }_p\hfill \\ 1,\hfill & x\in {\Omega }_s\hfill \end{array}$ (3)

通过使用以下最小二乘损失函数来描述椭圆界面问题的解$u\left(x\right)$ ：

$\underset{{}_{v\in H_g^2\left({\Omega }_p\cup {\Omega }_s\right)}}{\min }L\left(v\right)=\underset{{}_{v\in H_g^2\left({\Omega }_p\cup {\Omega }_s\right)}}{\min }{\displaystyle \sum _{i=1}^5{\beta }_i}{L}_i\left(v\right)$ (4)

其中${\beta }_i,i=1,2,3,4,5$ 是给定的正权重常数，$L_i\left(v\right),i=1,2,3,4,5$ 是对应的泛函，分别定义如下：

$L_1\left(v\right)={\displaystyle {\int }_{{\Omega }_p}{\left|-\nabla \cdot \left({\alpha }_1\nabla v_1\left(x\right)\right)-f\left(x\right)\right|}^2}$

$L_2\left(v\right)={\displaystyle {\int }_{{\Omega }_{\text{s}}}{\left|-\nabla \cdot \left({\alpha }_2\nabla v_2\left(x\right)\right)-f\left(x\right)\right|}^2}$

$L_3\left(v\right)={\displaystyle {\int }_{\Gamma }{\left|v_1\left(s\right)+p\left(s\right)-v_2\left(s\right)\right|}^2}$ (5)

$L_4\left(v\right)={\displaystyle {\int }_{\Gamma }{\left|{\alpha }_2\frac{\partial v_2\left(s\right)}{\partial n\left(s\right)}-{\alpha }_1\frac{\partial v_1\left(s\right)}{\partial n\left(s\right)}-q\left(s\right)\right|}^2}$

$L_5\left(v\right)={\displaystyle {\int }_{\partial \Omega }{\left|v_2\left(s\right)-g\left(s\right)\right|}^2}$

2.1. 多尺度神经网络

深度学习方法训练过程的本质是寻找最优参数使损失函数最小。在深度学习的过程中，神经网络需要对损失函数中的5个积分进行求值，通过蒙特卡洛算法对相应域中的随机点求和，可以有效且高效地逼近损失函数。在每个区域中引入均匀分布的随机点${\left\{x_i\right\}}_{i=1}^{M_p}\in {\Omega }_p,{\left\{x_i\right\}}_{i=1}^{M_s}\in {\Omega }_s,{\left\{x_i\right\}}_{i=1}^{M_{\Gamma }}\in \Gamma ,{\left\{x_i\right\}}_{i=1}^{M_b}\in \partial \Omega$ ，相应的优化任务变为：

$\underset{\Theta =\left({\Theta }_p,{\Theta }_S\right)}{\min }{\displaystyle \sum _{i=1}^5{\beta }_i}{l}_i$ (6)

函数$l_i,\left(i=1,2,3,4,5\right)$ 是$L_i$ 的离散形式，形式如下：

$l_1\left(\tilde{v}\right)=\frac{\left|{\Omega }_p\right|}{M_p}{\displaystyle \sum _{i=1}^{M_p}{\left|-\nabla \cdot \left({\alpha }_1\nabla {\tilde{v}}_1\left(x_i\right)\right)-f\left(x_i\right)\right|}^2}$

$l_2\left(\tilde{v}\right)=\frac{\left|{\Omega }_s\right|}{M_s}{\displaystyle \sum _{i=1}^{M_s}{\left|-\nabla \cdot \left({\alpha }_2\nabla {\tilde{v}}_2\left(x_i\right)\right)-f\left(x_i\right)\right|}^2}$

$l_3\left(\tilde{v}\right)=\frac{\left|\Gamma \right|}{M_{\Gamma }}{\displaystyle \sum _{i=1}^{M_{\Gamma }}{\left|{\tilde{v}}_1\left(s_i\right)+p\left(s_i\right)-{\tilde{v}}_2\left(s_i\right)\right|}^2}$ (7)

$l_4\left(\tilde{v}\right)=\frac{\left|\Gamma \right|}{M_{\Gamma }}{\displaystyle \sum _{i=1}^{M_{\Gamma }}{\left|{\alpha }_2\frac{\partial {\tilde{v}}_2\left(s_i\right)}{\partial n\left(s\right)}-{\alpha }_1\frac{\partial {\tilde{v}}_1\left(s_i\right)}{\partial n\left(s\right)}-q\left(s_i\right)\right|}^2}$

$l_5\left(\tilde{v}\right)=\frac{|\partial \Omega |}{M_b}{\displaystyle \sum _{i=1}^{M_b}{\left|{\tilde{v}}_2\left(s_i\right)-g\left(s_i\right)\right|}^2}$

其中${\tilde{v}}_i\left(x\right),\left(i=1,2\right)$ 是函数$v_i\left(x\right)$ 的近似值，$\tilde{v}\left(x\right)$ 是函数$v\left(x\right)$ 的近似值，$|\cdot |$ 代表的是测度。

由于逼近能力不同，神经网络中不同的拓扑或架构会影响所提出的深度学习方法的性能。在这里，我们没有使用常用的ResNet结构或前馈神经网络，而是使用了[16] [17]中的多尺度融合神经网络(MSFN)，它已被数值证明能够更好地捕捉“急转弯”，从而提高准确性。作为前馈神经网络的推广，MSFN结构基于小波理论中的多尺度基函数思想提出，并结合改进的信息融合模块，大大增强了其逼近能力。神经网络结构见图1：

Figure 1. Multi-scale neural network architecture diagram

图1. 多尺度神经网络结构图

MSFN在输入层接收到信息后，将其划分成多个并行子网络，并在每个子网络中还进一步加入了改进后的信息融合模块，以增强其逼近能力。最后，MSFN结构的输出是子网络的输出向量之和，加权平均作为可训练参数。基于多尺度神经网络求解椭圆界面问题，选择模型对参数进行更新时，使用LM优化算法。选择不同的激活函数对多尺度神经网络进行对比实验，并进行相关应用。

2.2. LM优化算法

在神经网络的参数训练过程中，必须借助优化方法解决目标函数的最小化问题。当前，Adam优化器因良好的稳定性与收敛效率，成为深度学习领域的主流选择。然而，这一选择在精度提升方面存在明显局限，若将LM [18]方法作为优化器引入深度学习框架，模型精度可实现显著突破，其相对误差能够达到10⁻⁶量级甚至更高，且在求解精度层面，表现可能优于传统的有限元方法。LM这类近似二阶方法在求解PDE问题时精度优势的根本在于：它通过利用损失函数的曲率信息，能够更准确地模拟损失曲面的局部几何形状，从而计算出比一阶梯度方向更优的更新方向和步长。这对于由平方和构成、且常呈病态条件的PDE损失函数尤为有效。

为进一步提升模型性能，本文选择LM方法作为核心优化器。考虑到内容完整性，下文将在具体实现部分简要阐述LM方法的基本原理与关键技术细节。针对深度学习场景中常见的函数最小化问题，我们考虑以下通用损失函数：

$\mathrm{Loss}\left(P\right)=\frac{1}{M}{\displaystyle \sum _{i=1}^M{\left(\varphi \left(x_i\right)-{\varphi }_{\text{aug }}\left(x_i;P\right)\right)}^2}$ (8)

其中$P\in {ℝ}^{N_p}$ 表示指定神经网络中的超参数向量，$x_i\left(i=1,\cdots ,M\right)$ 是采样点。我们利用LM算法得到其最小值，则有如下迭代方案：

$P^{\left(k+1\right)}=P^{\left(k\right)}+\delta$ (9)

其中$P^{\left(k\right)}$ 表示$P$ 的当前近似值，$\delta$ 满足以下修正的正态方程问题

$\left(J^{T}J+\mu I\right)\delta =J^T\left(\Phi -{\Phi }_{\text{aug }}\left(P^{\left(k\right)}\right)\right)$ (10)

通过归一化处理，我们将损失函数重写如下：

$\mathrm{Loss}\left(P\right)={\displaystyle \sum _{i=1}^M{\left(\frac{1}{\sqrt{M}}\varphi \left(x_i\right)-\frac{1}{\sqrt{M}}{\varphi }_{\text{aug }}\left(x_i;P\right)\right)}^2}$ (11)

接下来我们介绍LM算法的步骤，从而找到上述损失函数的最小点，具体如下：

LM算法：

步骤一：初始化容差$ϵ$ 和一个较大的$\mu$ ($\mu ={\mu }_{\max }={10}^8$ )，以便快速下降；设定${\mu }_{div}=3$ 和${\mu }_{\text{mul }}=2$ ；初始化超参数$P$ 、loss等；

步骤二：使用LM方法更新超参数$P$ 的值，计算新的损失，雅可比矩阵$J$ ，并设$\mu =\max \left\{\mu /{\mu }_{div},{10}^{-15}\right\}$ ；

步骤三：对于当前近似值$P^k$ ，通过LM算法更新超参数$P^{\left(k+1\right)}$ ，计算增量$\delta$ 和新的损失；

1. 如果新的损失小于之前的值，则设置$\mu =\max \left\{\mu /{\mu }_{div},{10}^{-15}\right\}$ ，然后返回到步骤三，得到新的近似值$P^{(k+1)}$ ；

2. 否则计算$\cos \theta =\frac{\langle \delta ,{\delta }_{\text{old}}\rangle }{\Vert {\delta }_{\text{old }}\Vert \cdot \Vert \delta \Vert }$ ，即前一个和当前搜索方向之间夹角的余弦值：

(a). 如果$\left(1-\cos \theta \right)$ 乘以当前损失大于先前的损失，则放弃当前搜索方向并设置$\mu =\max \left\{\mu \ast {\mu }_{\text{mul }},{\mu }_{\text{max }}\right\}$ 。使用新的$\mu$ 求解修正的方程获得新的搜索方向$\delta$ 并返回到步骤三；

(b). 否则，接受当前搜索方向，然后更新相关信息，设置$\mu =$ $\max \left\{\mu /{\mu }_{\text{div }},{10}^{-15}\right\}$ 并返回步骤三；

(c). 当损失小于给定的容差$ϵ$ 时停止迭代。

通过上述的算法步骤从而找到损失函数最小值。

3. 数值实验

在这部分中，我们进行数值测试来验证理论结果。为此，考虑了三种数值实验，一种是规则界面的情况，一种是不规则界面的情况，以及一种三维椭球界面的问题。这里我们使用Python来实现我们的数值实验，并对每种界面使用了三种不同的激活函数进行对比实验。我们采用相对$L^2$ 误差来测量神经网络的精确解$u^{\ast }$ 和近似解$\widehat{u}$ 之间的差异，其定义如下：

${\epsilon }_{R{L}^2}=\frac{{\Vert u-\widehat{u}\Vert }_{L^2}}{{\Vert u\Vert }_{L^2}}=\frac{\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left|u^{*}\left(x_i\right)-\widehat{u}\left(x_i\right)\right|}^2}}}{\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left|u^{*}\left(x_i\right)\right|}^2}}}$ (12)

3.1. 规则界面

首先我们考虑一个圆心为(1, 1)，半径为2的圆域，界面$\Gamma$ 为相同圆心且半径为1的圆，其坐标为$\left(x_1,x_2\right)=\left(1+\cos \theta ,1+\sin \theta \right),\theta \in \left[0,2\pi \right)$ 。其中系数$\alpha$ 为

$\alpha =\{\begin{array}{ll}1,\hfill & x\in {\Omega }_p\hfill \\ 2,\hfill & x\in {\Omega }_s\hfill \end{array}$ (13)

其精确解为

$u\left(x_1,x_2\right)=\{\begin{array}{r}\hfill 2\tanh \left(x_1+x_2\right),\left(x_1,x_2\right)\in {\Omega }_p\\ \hfill \tanh \left(x_1+x_2\right),\left(x_1,x_2\right)\in {\Omega }_s\end{array}$ (14)

在$\Omega$ 内选取400个采样点，界面$\Gamma$ 和边界$\partial \Omega$ 选取200个采样点进行训练，测试数据由$\Omega$ 内均匀分布的10,000个采样点生成。考虑Adam和LM优化算法下的多尺度神经网络及Tanh、Sigmoid、SiLU三种激活函数的误差见表1：

Table 1. Circular interface error

表1. 圆形界面误差

优化算法	激活函数	相对$L^2$ 误差
Adam	Tanh	2.1810e-04
LM	Tanh	1.3241e-11
LM	Sigmoid	4.3823e-12
LM	SiLU	3.9786e-12

可以看出来，相对于Adam算法而言，我们所采用的LM算法下的误差量级显著下降，不同的激活函数对误差也有着一定程度的影响。其中LM算法下，激活函数为SiLU时我们得到的精确解、近似解和误差分布见图2：

Figure 2. Circular interface diagram

图2. 圆形界面图

3.2. 不规则界面

我们进一步考虑不规则花型界面问题，其精确解如下：

$u\left(x\right)=\{\begin{array}{cc}{e}^{r^2},& x\in {\Omega }_p\\ 0.1r^4-0.01\log \left(2r\right),& x\in {\Omega }_s\end{array}$ (15)

$r$ 表示点与原点之间的距离，计算域$\Omega$ 为圆心为原点半径为1的圆，花形界面定义为：$r=0.5+\frac{\sin \left(5\theta \right)}{7}$ ，其中系数$\alpha$ 为

$\alpha =\{\begin{array}{ll}1,\hfill & x\in {\Omega }_p\hfill \\ 10,\hfill & x\in {\Omega }_s\hfill \end{array}$ (16)

在${\Omega }_p$ 选取160个采样点，${\Omega }_s$ 选取320个采样点，界面$\Gamma$ 和边界$\partial \Omega$ 选取128个采样点进行训练，在${[-1,1]}^2$ 的平方域上生成201 × 201的均匀网格用作测试。同样考虑两种优化算法及三种激活函数下的多尺度神经网络误差见表2：

Table 2. Floral interface error

表2. 花形界面误差

优化算法	激活函数	相对$L^2$ 误差
Adam	Tanh	1.4700e-03
LM	Tanh	4.1374e-05
LM	Sigmoid	7.4543e-06
LM	SiLU	9.9306e-08

可以发现，当界面不规则时，同为Tanh激活函数时我们所采用的LM优化算法依旧比Adam算法有不小的精度提升；不同激活函数对误差的影响较大，当激活函数为SiLU时的近似效果最好，此时我们得到的花形界面图形见图3：

Figure 3. Flower pattern interface diagram

图3. 花形界面图

3.3. 三维案例

最后，我们考虑一个三维的例子。在这个例子中，考虑一个系数为常数的3D椭球界面问题。计算域定义为${\left[-1,1\right]}^3$ 的立方体。界面采用以原点为中心的椭球体表面的形式，其特征是长度为0.7、0.5和0.3的半主轴。这个椭球体由水平集函数$\phi ={\left(x/0.7\right)}^2+{\left(y/0.5\right)}^2+{\left(z/0.3\right)}^2-1$ 表示。对于该问题，子域系数$\alpha$ 都取1，精确解为

$u\left(x,y,z\right)=\{\begin{array}{ll}\exp \left(x^2+y^2+z^2\right)\hfill & \left(x,y,z\right)\in {\Omega }_p\hfill \\ 0.1{\left(x^2+y^2+z^2\right)}^2-0.01\log \left(2\sqrt{x^2+y^2+z^2}\right),\hfill & \left(x,y,z\right)\in {\Omega }_s\hfill \end{array}$ (17)

在$\Omega$ 选取1500个采样点，界面$\Gamma$ 和边界$\partial \Omega$ 选取240个采样点进行训练，在${\left[-1,1\right]}^3$ 的平方域上生成128 × 128 × 128的均匀网格用作测试。和前面的数值实验采用一样的对比实验，得出三维案例下的相关误差见表3：

Table 3. Three-dimensional ellipsoid interface error

表3. 三维椭球界面误差

优化算法	激活函数	相对$L^2$ 误差
Adam	Tanh	3.1283e-03
LM	Tanh	3.1265e-05
LM	Sigmoid	4.1338e-07
LM	SiLU	2.3768e-07

其中，该三维案例在$z=0$ 切平面处的精确解、近似解和误差分布见图4：

Figure 4. Cross-section of a three-dimensional ellipsoid interface

图4. 三维椭球界面切面图

4. 结论与展望

本文基于LM算法的多尺度神经网络求解椭圆界面问题，并将其与常规的Adam优化算法展开对比分析。研究结果表明，在三个典型数值案例的求解过程中，所采用的LM算法展现出更优异的近似效果，相较于Adam算法，其求解精度的数量级实现了显著提升，这一结果充分验证了本研究的可行性与有效性。

然而，需要客观指出的是，LM算法在处理高精度优化问题时，也存在明显的性能瓶颈。当将该算法拓展至大规模椭圆界面问题求解场景时，其计算效率不足的问题逐渐凸显，对硬件计算设备的算力、存储容量等性能指标提出了更高要求——这不仅增加了实际工程应用中的部署成本，也在一定程度上限制了该方法的适用范围。因此，如何优化LM算法的计算流程以提升效率、降低设备依赖，成为当前方案亟待突破的核心问题，有待进一步深入探索与完善。