PINN 训练体会

3 PINN真的慢，真的不稳定，不用trick效果很烂，简单toy model很难超过传统数值

PINN天生就被大家寄予厚望，然而他在处理小规模问题取得了让人哄堂大笑的结果。真不知道当初Rassi是怎么说服审稿人这个方法有希望的。而关于PINN是不是没活硬整的算法已经讲了很多了，这里我仅仅介绍一些训练时的真实体会。

1为什么PINN这么慢？

因为神经网络非线性层+自动微分架构。在数据科学领域他是一把”精密手术刀“，但在计算数学的领域我们要的是“激光点打击”。一个两层的神经网络和一个1000阶的线性方程组相比：恕我直言，前者还是有待加强

2PINN到底是一个什么样的问题？

所有人看到Loss，这是一个优化问题，到底谁才是主项谁才是约束？这是一个很经典的“数学上等价但是数值上不等价的问题”

答案你可能意想不到：其实我们实际计算的过程中，是在PDE项增加调参系数，是将其视为一个PDE约束的逼近优化问题。

一般来讲，我们是用一个好项去balance一个坏项，例如说用L1正则项让问题稀疏化。但是确实从直觉的角度来看，我们要解的是方程，边界只是一个约束。但现在我们逼近的是边界，而方程只是约束条件。这个说法就听上去有点奇怪。

但如果我说，边界项是我们观察到的已知数据，方程是作为物理信息“嵌入”神经网络，也就是打包进Loss呢？这么说来看，PDE才是那个“来者不善”的来者。

所以，我这里采用分阶段训练居然意外的取得了很好的效果。从1e-1一直连降了三个数量级，虽然依然效果很烂就是了

而在此基础上，我们可以用控制的视角去看PINN，毕竟，将PDE作为约束就是那类问题所做的。

下期我或许会做这个？

3 方程艰险，数据如何求解？PINNs空耗算力，不如回归传统。数据规模太小，谁能做求解器？

PINN更多的是作为一种全新的范式，就从他发表的JCP之中，就可以看出，它实际上目标并非是收敛性的分析，而是一个大胆的尝试。用更编程的一些说法，PINN有点类似于Python，虽然它运行效率很慢，但是方便快捷，一列Loss的轮椅打法本身就是一大优势。甚至于我们嘲笑说自然语言编程，但现实就是在收敛速度分析之外，计算数学的初衷或许就是让更多人面对一个工程物理问题时能够算出来，将大家从痛苦的matlab画格子的痛苦解救出来，并且考虑到python的训练语料真的很多，至少远远多于Matlab，与其调包，不如加入AI coding的怀抱中

最后，所有人，立刻打开github,关注deepxde，里面有很多的代码

书接上回，我们尝试着训练一个网络，然而，可谓是拼尽全力无法战胜loss曲线。PINN的调参确实是玄学……吗？

纯新人的PINN训练指北7 赞同 · 2 评论文章

让我们掌声有请deepmind，深度学习高手为我们带来大名鼎鼎的PINN-NTK自适应调节方法。

在此之前，我们不得不问一件事情：

什么是NTK

NTK（Neural Tangent Kernel）神经正切核讨论的是一类极限情况，对于宽度趋近于无穷的神经网络：

初始化时：神经网络 ≈ 高斯过程（GP）

f(x,θ)∼N(0,Σ(x,x′))

训练时：神经网络 ≈ 核回归

fNTK(x)=Kx,training⋅Ktraining,training−1⋅y

而且核矩阵在训练过程中保持常数

为什么要用NTK视角分析PINN？

回到PINN问题本身，我们会发现PINN常常会出现这三大问题：

PINN训练为什么有时候会失败？
为什么会出现梯度不平衡？
如何科学地选择超参数？

而面对着三个挑战，NTK 理论为其提供了严谨的理论分析框架。

其以神经网络隐藏层宽度趋于无限为核心假设，将复杂的 PINN 训练转化为可解析的核回归问题，让收敛速率、谱特性等关键训练规律变得定量且明确。

在这个基础上，我们能精准剖析 PINN 训练失败的根源、梯度不平衡的本质，并为超参数选择提供科学的理论依据，而非经验性试错。

具体操作：计算PINNs的NTK

首先我们来计算一下，对于这样的深度学习网络之中的NTK正切核矩阵到底是怎么样的

考虑一个简单的边值PDE问题：其中 L 表示一个微分算子

Lu(x)=f(x),x∈Ω

u(x)=g(x),x∈∂Ω

可以定义一个PINN损失函数：

L(θ)=Lb(θ)+Lr(θ)

其中边界损失：

Lb(θ)=12Nb∑i=1Nb|u(xib,θ)−g(xib)|2

PDE残差损失：

Lr(θ)=12Nr∑i=1Nr|Lu(xir,θ)−f(xir)|2

NTK矩阵定义：

K(t)=(Kuu(t)Kur(t)Kru(t)Krr(t))

其中：

(Kuu)ij=∂u(xib,θ)∂θ⋅∂u(xjb,θ)∂θ

(Krr)ij=∂Lu(xir,θ)∂θ⋅∂Lu(xjr,θ)∂θ

谱偏差：为什么PINN学得慢？

在我们计算出核矩阵之后，我们就可以基于NTK的有关结论来说明问题的来源

NTK核矩阵，其特征值分布为

K=QTΛQ,Λ=diag(λ1,λ2,…,λN)

此时，训练误差衰减：

‖error(t)‖≈‖e−Kt‖

不同的特征值会带来谱偏差现象：

大特征值 → 快速收敛（低频成分）
小特征值 → 缓慢收敛（高频成分）
特征值快速衰减 → 高频成分难以学习

PPT图示（文字描述）：

特征值衰减曲线：
λ ████████████░░░░░░░░░░░
  ▏
  ▏
  ▏
  ▏低频 ←→ 高频

因此：NTK的特征值分解揭示了PINN训练困难的根本原因——谱偏差。较大的特征值对应的模式会快速收敛，而较小的特征值（通常对应高频成分）收敛极慢。这解释了为什么含有高频特征的解将会使得核矩阵的特征值奇异，使得它难以用PINN准确学习。

损失函数失衡：为什么收敛速度很慢

除了谱偏差外，PINN的另一个核心问题是不同损失项之间的不平衡。边界条件和PDE残差的收敛速度可能差异巨大，导致训练过程中需要不断调整权重。

PINN的损失函数包含多个项：

L=λbLb+λrLr

Lb 和 Lr 的收敛速度可能相差数十个数量级
导致某些项已收敛而另一些项仍很大
传统的固定权重无法解决

PPT图示：

训练过程：
损失值
  ↑
  │    _______ 边界损失（快）
  │   /
  │  /
  │ /
  │/__________________ → 迭代次数
  │    ~~~~~~~ PDE损失（慢）

算法：动态权重调整

每 k 步更新权重：

λb(k+1)=λb(k)⋅tr(Krr)tr(Kuu)

λr(k+1)=λr(k)⋅tr(Kuu)tr(Krr)

或者使用梯度归一化：

g~b=gb‖gb‖,g~r=gr‖gr‖

Lnew=Lb+Lr

文案：基于NTK分析，作者提出了一种自适应训练算法。核心思想是根据NTK特征值的trace来动态调整不同损失项的权重，使得不同组件的收敛速度达到平衡。

3.8 数值实验

PPT内容： 实验1：一维Poisson方程

uxx=f(x),x∈[−1,1]

u(−1)=u(1)=0

真解：u(x)=sin⁡(πx)

结果对比：

方法	相对L2误差
标准PINN	1.2e-2
+ 动态权重	3.1e-4

实验2：多尺度问题 高频解 u(x)=sin⁡(10πx)

方法	相对L2误差
标准PINN	8.5e-1
+ 动态权重	2.3e-2

文案：数值实验验证了NTK理论的分析和自适应算法的有效性。在高频和多尺度问题上，动态权重调整显著提升了PINN的性能。

偏微分方程：

Cauchy-Riemann方程是一类什么样的方程：

∂u∂t=∂v∂x∂v∂t=−∂u∂x 其次，构造了如下的能量方法： E(t)=∑k=0∞(ρkk!)2∫(|∂xku|2+|∂xkv|2)dx

对E关于t求导，考虑 ρ=ρ(t) ，得到两项

E(t)=∑k=0∞2kρ(ρkk!)2∫(|∂xku|2+|∂xkv|2)dx+∑k=0∞(ρkk!)2ddt∫(|∂xku|2+|∂xkv|2)dx

这是两个部分组成的，可以考虑一个函数：

Y(t)=∑k=0∞2kρ(ρkk!)2∫(|∂xku|2+|∂xkv|2)dx

因此函数的第一项实际上就是 ρ′(t)Y(t)

第二项则考虑进行放缩

∑k=0∞ddt∫(|∂xku|2+|∂xkv|2)dx≤∑k=0∞∫(∂xku∂xk(∂tu)+∂xkv∂xk(∂tv))dx≤∑k=0∞∫(∂xku∂xk+1v)−∂xkv∂xk+1udx≤∑k=0∞(∫∂xku2)12(∫∂xk+1v2)12+∑k=0∞(∫∂xkv2)12(∫∂xk+1u2)12=A

这样的到了相似的两项，继续进行放缩：

∑k=0∞(ρkk!)2(∫∂xkv2)12(∫∂xk+1u2)12≤∑k=0∞kρ(ρkk!)2(∫∂xkv2)+∑k=0∞ρk(ρkk!)2(∫∂xk+1v2)=∑k=0∞kρ(ρkk!)2(∫∂xku2)+∑k=0∞k+1k(ρk+1k+1!)2(∫∂xk+1v2)

而对于另一项，求和同样可以得到：

A=Y(t)+∑k+1k(ρk+1k+1!)2(∫∂xk+1v2)≤Y(t)+2Y(t)

因此 E(t)≤(3+ρ′(t))Y(t)

我们先来考虑一个简单的二阶线性偏微分方程：

u-\Delta u = f

需要给出一个弱解u\in H_0^1(\Omega)，即在\Omega这个区间内存在紧支撑集，且一阶弱导数与原函数平方可积。

而弱解，实际上就是要满足对于任意的g \in H_0^1(\Omega)

\int_{\Omega} (ug-\Delta u g)dx = \int fg

这里就使用小名鼎鼎的Reisz表示定理，对于一个Hilbert空间上的内积(u,g)

(u,g)=Tg存在唯一解，当且仅当T是一个有界算子

那么，问题就转化为了：

1(u,g)是不是一个内积

2Tg是不是一个有界算子

首先我们回答第一个问题，H_0^1空间内这是不是一个内积，或者说，对于这样的Soblev空间如何定义内积。

要想证明内积，实际上是可以通过范数诱导得到。根据H_0^1空间的定义:

||u||H02=||u||L22+∑j=1||∂xju||L22

这样定义实际上希望满足Cauchy性，一个解的Cauchy列的极限希望是这个解，而对于偏微分方程就需要原函数接近的同时导数也接近，因此这样定义标准内积。

在这个基础上，用范数诱导内积。得到：

H01=∫Ωug+∇u⋅∇g

这实际上可以通过积分得到，考虑到g在H_0^1上，所以在边界处取零，无需考虑边界项，因此

我们确实可以把原问题转化为 ∫Ωug+∇u⋅∇g=∫Ωfg

接下来，我们就需要回答两个问题，这么诱导能不能诱导出一个内积，以及\intfg是不是有界的。

内积是满足三条性质的双线性函数:1对称，2双线性，3正定。

对称性=,双线性=a+b，正定 \ge 0确实都是满足的。

这三条性质同时满足，这确实是一个内积。而且在H_0^1的框架下，是一个有界的量。

接下来，我们验证右边是一个有界算子，将f视作一个固定的函数，那么

Tu:u→∫Ωfu 是一个从Soblev空间到实数的一个泛函，而且也确实满足：

∫Ωfg≤∫Ω|f|2dx∫|g|2≤∫Ω|f|2dx∫Ω(|g|2+|∇g|2)dx=C||g||H01

在凑成Reisz表示定理的羁绊之后我们直接套用结论，很好，结果是有界的。

不过我们还有个疑问，如果没有凑出这个特有的内积形式，那怎么办？举个例子，我们现在要面对这样的问题：

\int_{\Omega} \nabla u \cdot \nabla g dx= \int_{\Omega} fg dx

我们需要证明左边与原内积依然是等价的，所谓等价，就是只差一个常数

这就需要使用Poincre不等式证明：

∫f2dx≤C∫|∇f|2dx

对于一个有界区域这是成立的，因为我们可以先考虑一维的情形

∫f2dx≤∫ab∫axff′dsdx≤∫abfdx∫abf′≤(b−a)(∫abf2dx)1/2(∫ab(f′)2dx)1/2

而对于高维的情况，已知这样对于每一个变量都成立，即：

∫f2dx≤C∫∂xif2dx 通过取最大 sup{|xj−yj||x,y∈Ω} 的方法给出一个估计： ∫f2dx≤C∫|Df|2dx

但这样的常数选取实际上依赖于有界区间,在此基础上，介绍另一个不等式，Hardy不等式，适用于处理误解区间的情形。

通过这些证明原函数可以被高阶导数控制，可以证明C||\nabla u|| \ge ||u||_{H_0^1} 因此可以得到两种内积的等价性，证明 ∫Ω∇u⋅∇g=∫Ωfg 也存在唯一解

当然，我们还可以证明更一般的情形，也就是Lax-Milgram定理所阐述的情形

数值解：有限元方法

1偏微分方程的弱形式：

考虑一个经典的热方程： ∇⋅(c∇u)=−f

考虑Dirchlet边界条件：，u=g，x∈∂Ω

这个方程差不多已经出现过很多次，我们引入变分操作，实际上就是把能够乘进去的v全都乘一遍，然后将其做积分

∫∇⋅(c∇u)vdx=−∫fvdx

我们要求对于任意的v都能够满足方程，这样就能求得一个解。

所以什么是“能够乘进去的v”呢？可以考虑一个好函数空间 ϕ(x) ：

好函数要求是无限可微的，而且好函数空间具有紧支撑集。

你可能不需要知道这是什么意思，在此处就是要让我们不需要考虑边界项，而且无限光滑，也就是让它可以随便积分随便求导。

有一个“好函数”作伙伴，”不那么好函数“u(x)，只需要满足可积性条件，也可以求导了

∫v(x)ϕ(x)dx=−∫u(x)ϕ(x)dx

如果对任意的好函数\phi(x)都能满足这样的题意，那么就称此时的v是一个弱导数，此后我们如果不加以说明，所有的导数都是以弱导数为基础，包括梯度项

原本的可积性条件就是平方可积 ∫f2dx<∞ ,此后不加以说明，

回到原本的变分问题：

如果我们对于v加一些限制，它需要在\Omega上是一个紧支集，意味着它在边界处取零

这样就可以化简为一个对称的: ∫c∇u(x)⋅∇v(x)dx=∫fvdx

基于此，我们就得到了这个方程解的存在空间，需要一阶导数L^2与原函数L^2平方可积。

平方可积也并非是随便搞的，原函数的平方可积是基础，导数平方可积这实际上就是利用Cauchy不等式：

∫c(x)∇u⋅∇v≤||c(x)||∞∫|∇u(x)|2dx∫|∇v(x)|2dx<∞

因此，我们就得到了Soblev空间，这在偏微分方程界的地位相当于线性方程组之中一个欧式空间R^n

Reisz定理说明，这样的方程必然是存在唯一解的，但是R老师并没告诉怎么求解这个问题。

到头来，往往需要找一个有限维的子空间，通过在这个子空间上得到解。至于有限维怎么求解，当然是构造线性方程组然后求解线性方程组。

这样的子空间也就是大名鼎鼎的”有限元空间“，构造有限元空间，生成有限元基函数并且求解有限元方程组的整个过程就被称之为”有限元方法“。

回到原问题之中：如果我们已知函数处在一个有限维的有限元空间 Uh ，从中任意选取试探函数 vh ，都存在 uh

∫c(x)DuhDvhdx=∫fvhdx ,若将左项写成 a(uh,vh)=∫Ωc(x)uhvh 右项写成内积的形式

也可以写成，a(uh,vh)=(f,vh)，∀vh∈Uh ,

已知有限元空间是一个有限维的网络，因此我们可以选取一组基 {ϕ1,ϕ2,⋯,ϕn} , uh=∑uiϕi

其中 a(uh,ϕi)=∑uja(ϕi,ϕj)=(f,ϕj)

这实际上最终将会构成一个矩阵（）A=（a(ϕi,ϕj)）u=(ui),v=(f,ϕj)

这样实际上就构成了一个线性方程组Au = v

关键在于如何求解这个刚度矩阵A之中的每个元素了，而且既然基函数是我们亲自选的，这个A当然得是一个好矩阵：它是稀疏的，大部分元素都是零。

那么问题就在于，怎样选取基函数？怎样组装刚度矩阵?

我们举个最简单的1d网络剖分为例：[a，b]区间内，划分n个区域共0-N，总共N+1个点

这样就很自然的得到一组基所需要满足的条件： ϕj(xi)=δij ，n个点，n个基，可以确定一个n次方程组，得到一组拉格朗日基底。

但还需要考虑稀疏性：因此需要构造一个分片多项式的函数，例如线性函数就是一次元，二次函数就是二次元，以此类推。

举个线性元的例子：假设考虑

那就是 ϕj(x)=(x/h−(j−1))I[j−1h,jh]+((j+1)−x/h)I[jh,j+1h],j=1,2,3,⋯,N−1

（ϕ0(x)=（1−x/h)I[0,h],ϕN(x)=(N−1)h+x/h)I[N−1h,Nh]

这样，就可以考虑根据分片性质得到稀疏矩阵，仅仅在主对角线和副对角线上生成一个三对角矩阵。此时考虑从1到N+1开始编号

首先考虑i=2,3,\cdots,N

对角线上 Aii=1h2∫i−2hihc(x)dx,Aii−1=−1h2∫i−2hi−1hc(x)dx,Aii+1=−1h2∫i−1hihc(x)dx

i=1,N+1时需要考虑边界条件A_{1,1}=1/h^2\int_0^hc(x)dx,A_{N+1,N+1}=1/h^2\int_{N-1h}^{Nh}c(x)

（）bj=（f,ϕj）=∫(j−2)hj−1hf(x)(x/h−j+2)dx+∫(j−1)hjhf(x)(−x/h+j)dx,j=2,3,4,⋯,N=∫0hf(x)(1−x/h)dx,j=1=∫N−1hNhf(x)(x/h−(N−1)h)dx,j=N+1

与此同时，还需要对于边界项进行额外的考虑，对于1d网络，实际上就是要让

u1=ua,uN+1=ub

调整刚度矩阵

A[1,:]=0,A[1,1]=1,A[N+1,:]=0,A[N+1,N+1]=1,b[0]=ua,b[N+1]=ub

有限元空间划分网格-构造基函数-计算参考点-组装刚度矩阵-考虑边界条件-求解线性方程组，这样我们就求解一个有限元解。

而对于2d和更高次的网络，往往需要考虑网络本身的性质，考虑选取一个划分

一方面可以使用三角形网格，另一方面可以采用矩形网络。一般情况下，不规则网络采用三角形划分，规则网络采用矩形划分。

我们首先考虑选取N个有限元单位,N_m个节点

E_n表示第n个有限元单位，Z_k为第k个节点。

P是一个N_m*N_m维信息矩阵，表示一个邻点矩阵。

T是一个N*N_m维信息矩阵，表示每个有限元单位上所包含节点

首先，依然是需要首先构造一个有限元函数空间，这包括构造有限元函数的基。这实际上依然是类似于一维时的操作，构造一组拉格朗日基底。

考虑X_k一系列的节点，共有N_b个节点，依然是可以构造对应的拉格朗日基底函数。

三角形网格上的线性元
1 参考单元上的基函数对于参考三角形 K^（顶点 (0,0),(1,0),(0,1)）：

节点基函数 ϕi(ξ,η) (0,0) ϕ1=1−ξ−η (1,0) ϕ2=ξ (0,1) ϕ3=η

2.2 实际单元的映射

实际单元 K 的顶点为 (x1,y1),(x2,y2),(x3,y3)，映射关系如下：

[xy]=[x1y1]+[x2−x1x3−x1y2−y1y3−y1][ξη]

雅可比矩阵 J 定义为：

J=[∂x∂ξ∂x∂η∂y∂ξ∂y∂η]

2.3 单元刚度矩阵计算

对于单元 En，其刚度矩阵元素定义为：

An,αβ=∫Enc(x)∇ϕnα⋅∇ϕnβdx

参考单元上的积分转换

通过坐标映射将实际单元积分转化为参考单元 K^ 上的计算，公式变为：

An,αβ=∫K^c(x(ξ,η))(∇ϕα)T(JJT)−1(∇ϕβ)|detJ|dξdη

线性元的梯度特性

对于三角形线性单元，形函数的梯度为常数，具体为：

∇ϕ1=[−1−1],∇ϕ2=[10],∇ϕ3=[01]

常数系数简化（c(x)=cn）

若扩散系数 c(x) 在单元 En 内近似为常数 cn，则单元刚度矩阵可简化为：

An=cn|detJ|⋅BT(JJT)−1B

其中，矩阵 B 由形函数梯度构成：

B=[∇ϕ1,∇ϕ2,∇ϕ3]

2.4 载荷向量计算

单元 En 的载荷向量元素定义为：

Fn,α=∫Enf(x)ϕnα(x)dx

参考单元上的积分转换

通过坐标映射将实际单元积分转化为参考单元 K^ 上的计算，公式变为：

Fn,α=∫K^f(x(ξ,η))ϕα(ξ,η)|detJ|dξdη

首先考虑三角形的基底，构造一个(x_i,y_i)的三个点，这将会对应三个基函数。

我们可以通过规范三角形，三个点分别为(0,0),(0,1),(1,0)，分别构造三个basis函数：

对应的三个基函数分别是1-x-y,x,y，取值限制在这个小三角形单元上。

对于第n个小三角形单元，可以得到三个基函数\phi_{p_s}(x),p表示节点标号。

在得到基函数之后，就可以得到刚度矩阵A,A刚度矩阵实际上有3*N行，其中上述的积分实际上是T()

本质上，狄利克雷边界条件*u*=*g*给出了所有边界有限元节点处的解。

由于有限元解*uh*=∑_{j=1}^{Nb}**u_j**φ_j*中的系数*uj，实际上就是有限元节点*Xj*（j=1,⋯,Nb）处的数值解，

因此我们其实已知那些对应于边界有限元节点的*uj*。

需注意，boundarynodes(2,:)存储了所有边界有限元节点的全局节点编号。

若*m*∈boundarynodes(2,:)，则第*m*个方程被称为边界节点方程。

设*nbn*为边界节点的数量；