PaPi: Parallel Picard-Duhamel Integration for Fast Neuro-Spectral Architectures
Authors: Arthur Bizzi, Márcio Marques, Kleyton Da Costa, Vitor Balestro, Tiago Novello, João M. Pereira, Daniel Yukimura
Organizations: Imperial College London London, United Kingdom · Instituto de Matemática Pura e Aplicada Rio de Janeiro, RJ, Brazil · Computer Science University College London London, United Kingdom · Instituto de Matemática e Estatística Universidade Federal Fluminense Niterói, RJ, Brazil · Department of Mathematics University of Georgia Athens, GA, USA · Instituto de Matemática Pura e Aplicada Rio de Janeiro, Brazil
Neuro-spectral architectures (NeuSA) are a class of physics-informed neural networks (PINNs) that have proven effective for solving time-dependent PDEs in the frequency domain. However, training and inference often rely on costly numerical integration of the learned vector field, typically using sequential Runge-Kutta time stepping, which limits temporal parallelism. We introduce PaPi, an integration scheme inspired by weakly nonlinear flows that enables time-Parallel inference through exponential integration and a Duhamel-Picard scheme. We show that the Duhamel-Picard iteration converges exponentially for all flows, and provide explicit error bounds that depend on how close the flow is to linear. Experiments show that PaPi accelerates training and inference by up to two orders of magnitude, achieving training times of under one second for low-dimensional PDE learning problems, while maintaining competitive predictive accuracy.
Figures & tables
2D Layers
Marmousi
3D Layers
2D Burgers
Sine-Gordon
KdV
Model
rMSE
TT
rMSE
TT
rMSE
TT
rMSE
TT
rMSE
TT
rMSE
TT
PINN
0.549
420
0.618
427
0.062
4000
0.183
605
0.130
569
6.1e-3
908
QRes
0.150
765
0.401
707
0.026
5496
0.073
1178
0.013
781
9.5e-3
1240
FLS
0.492
408
0.712
408
0.071
4050
0.174
578
0.129
594
1.5e-2
935
NeuSA
0.069
508
0.266
604
0.008
577
0.053
119
9.0e-4
234
3.6e-3
1322
PaPi (ours)
0.068
12
0.154
66
0.018
251
0.028
30
5.6e-4
0.3
1.2e-3
3.2
Table 1: Quantitative comparison. Mean rMSE and training time over five seeds for each model. TT denotes training time in seconds. PaPi has the shortest training time on all six problems and the lowest rMSE on five. Appendix C reports standard deviations and rMAE.
Figure 1: Wave equation in the Marmousi medium at t=2 (seed 0): the displacement u on the region of interest [−2,2]2 (top) and the absolute error against the reference (bottom), with the relative L2 error at this time in each error tile. PaPi resolves the reflected and refracted fronts that NeuSA blurs and the PINN misses.
Noise
Model
rMSE
∣λ1−1∣
∣λ2−0.1∣/0.1
TT (s)
0%
PINN (5,000 steps)
9.9e-2 ± 4.1e-2
4.9e-1 ± 2.0e-1
8.2e-2 ± 3.5e-2
419
QRes (5,000 steps)
4.0e-2 ± 2.1e-2
2.0e-1 ± 1.0e-1
3.7e-2 ± 1.9e-2
561
FLS (5,000 steps)
9.9e-2 ± 4.1e-2
4.9e-1 ± 2.0e-1
8.1e-2 ± 3.5e-2
411
PaPi (500 steps)
1.0e-2 ± 8.4e-3
4.7e-2 ± 4.1e-2
9.4e-3 ± 8.2e-3
0.6
PaPi (5,000 steps, float64)
3.8e-5 ± 1.2e-5
1.1e-4 ± 5.3e-5
3.6e-5 ± 2.2e-5
66
1%
PINN (5,000 steps)
9.8e-2 ± 4.1e-2
4.9e-1 ± 2.0e-1
8.0e-2 ± 3.5e-2
432
Table 2: Inverse Burgers’ problem: mean ± std over the ten starting pairs of the relative L2 error on [0,1] , the relative errors of the identified coefficients and the training time per run. The observation is the final frame only; noise is relative to the standard deviation of the data.
Appendix figures & tables11 assets
Supplementary material from the paper’s appendix.
Appendix
Problem
Basis
Modes
nt ( T )
ϵ
Steps
lr
K
Quadrature
Train nodes
Deploy
Params
2D Layers
cosine
2012
201 (2)
1
2,000
10−2
2
exp. trapezoid
101
8
3.2e5
Marmousi
cosine
2012
201 (2)
2
2,500 (2,000)
cosine
6
trapezoid
–
8
3.2e5
3D Layers
cosine
1013
111 (1)
10
4,000 (2,000)
10−2
2
exp. trapezoid
56
8
6.1e4
2D Burgers
Fourier
2012
201 (1)
0.1
1,000 (200)
5⋅10−3
2
trapezoid
–
8
8.1e5
Sine-Gordon
sine
201
201 (3)
0.1
1,000
2⋅10−2
2
trapezoid
–
8
9.7e5
KdV
Fourier
81
1,001 (1)
0.1
10,000 (1,000)
10−2
2
trapezoid
–
32
1.6e5
Appendix
Table 3: Settings per problem. nt is the number of time nodes on [0,T] ; “train nodes” the coarser grid the training rollout uses with the exponential trapezoidal rule (–: full grid); K the warm-started sweeps per optimiser step; “deploy” the number of Picard sweeps at inference. Steps are PaPi’s, with NeuSA’s in parentheses when different.
Figure 2: Training time against rMSE on every benchmark of Table 1 (means over five seeds, one RTX 5090). The arrow joins NeuSA to PaPi: the same model, trained through the parallel Picard–Duhamel rollout instead of the sequential RK4 one.
Figure 3: Relative L2 error at each time t (mean over five seeds, band from best to worst seed) for PaPi, NeuSA and the best MLP-based baseline, QRes.
Problem
Model
rMSE
rMAE
TT (s)
2D Layers
PINN
5.5e-1 ± 8.9e-2
7.8e-1 ± 1.3e-1
420 ± 56
QRes
1.5e-1 ± 7.2e-2
1.9e-1 ± 7.5e-2
765 ± 212
FLS
4.9e-1 ± 5.2e-2
7.2e-1 ± 8.3e-2
408 ± 60
NeuSA
6.9e-2 ± 2.0e-2
9.4e-2 ± 3.1e-2
508 ± 19
PaPi
6.8e-2 ± 1.8e-2
9.1e-2 ± 2.7e-2
11.8 ± 0.0
Marmousi
PINN
6.2e-1 ± 3.5e-2
9.2e-1 ± 8.1e-2
427 ± 68
Appendix
Table 4: Main benchmarks: mean ± std over five seeds, all on the same RTX 5090.
Figure 4: rMSE of the trained PaPi models after k cold Picard–Duhamel sweeps from the linear flow, mean over five seeds. Every problem but KdV is converged by eight sweeps.
NeuSA (PyTorch)
RK4 (JAX)
PaPi (JAX)
speedup / step
Problem
rMSE
steps
TT (s)
rMSE
steps
TT (s)
rMSE
steps
TT (s)
torch/JAX
RK4/PaPi
2D Layers
6.9e-2
2,000
508
9.1e-2
2,000
117
6.8e-2
2,000
11.8
4.4
9.9
Marmousi
2.7e-1
2,000
604
1.9e-1
5,000
290
1.5e-1
2,500
66
5.2
2.2
3D Layers
7.6e-3
2,000
577
5.3e-3
4,000
676
1.8e-2
4,000
251
1.7
2.7
2D Burgers
5.3e-2
200
119
2.5e-2
1,000
123
2.8e-2
1,000
29.6
4.8
4.2
Sine-Gordon
9.0e-4
1,000
234
7.3e-4
1,000
16.7
5.6e-4
1,000
0.3
14
56
Appendix
Table 5: Same model, three integrators, same RTX 5090: mean over five seeds of rMSE, optimiser steps, training time and time per step. “NeuSA” is the authors’ PyTorch code (torchdyn RK4); “RK4 (JAX)” our port of the same rollout; both ratios are of times per step.
Time nodes
Picard, K=16
Duhamel, K=1
PD, K=16
PD time (ms)
101
7.55×1017
5.025
8.63×10−4
3.50
201
7.04×1017
5.021
2.16×10−4
5.10
401
6.88×1017
5.018
5.39×10−5
7.76
Appendix
Table 6: Only PD converges. Relative error E at fixed sweep budgets; times are medians of 20 warmed float64 GPU solves, excluding compilation.
Figure 5: Each PD sweep extends the accurate time window, while Picard fails early at every budget. (a,b) Error at each time t on the 401-node grid (gray: above one); the bottom PD row is the Duhamel solution. (c) Once converged ( K=32 ), PD error falls at exactly second order in Δt .
Figure 6: Step size and optimizer family determine accuracy; regularization only adds error. (a) Mean error per optimizer and initial learning rate. (b) L2 penalty for Adam at η0=0.01 . (c) Polyak SGD step-size cap. (d) All configurations against training time; the star marks Adam at η0=0.01 . Error bars are standard deviations over three seeds.
Figure 7: The PD gain comes from the update, and part of the training error comes from time integration. (a) Both updates from both initial paths (401 nodes). (b) Exact-force RK4 in float32 and float64. (c) Frozen Adam and AdamW models evaluated with finer RK4. Panels (b,c) score the original 201 nodes.
Figure 8: Given enough updates Muon matches Adam, while float64 does not rescue L-BFGS. (a) Extended step-size ranges at 1,000 updates. (b) Residuals continued to 5,000 updates: seed median with min–max band, smoothed by a 51-update running median (dotted line: original budget). (c) L-BFGS in float32 and float64 from identical initial parameters.
Neural PDE solvers use trainable models to approximate solutions or solution operators for partial differential equations, often with physical constraints.
Physics-Informed Neural Networks (PINNs) build neural representations of time-dependent PDE solutions, naturally incorporating physics knowledge and observational data, which makes them well suited to both forward and inverse PDE problems. PINNs, however, are known to suffer from spectral bias and lack of causality. Neuro-Spectral Architectures (NeuSA), a recently proposed alternative to PINNs, mitigate both issues, but their numerical integration becomes unstable for stiff differential equations arising in many relevant physical problems. This study proposes Neuro-Spectral Exponential Time Differencing Architectures (NEXT), which combines the spectral representation of the PDE solution in NeuSA with high-order exponential integrators. Within this approach, the linear stiff part of the vector field induced by the PDE is integrated exactly through matrix exponentials, while the possibly nonlinear remainder is modeled by a neural network. The effectiveness of NEXT is verified through benchmark experiments on a set of stiff PDEs, in which NEXT is stable and accurate while NeuSA diverges numerically. It is also shown that NEXT can be applied to inverse problems, where the model has to learn unknown parameters or boundary conditions from sparse data. All code used in this work is publicly available at: https://github.com/marcioh2m/next.git .
Márcio Marques, Leonardo Mendonça, Leonardo M. Moreira +6
Instituto de Matemática Pura e Aplicada (IMPA), Rio de Janeiro, Brasil · Universidade do Estado do Rio de Janeiro (UERJ), Rio de Janeiro, Brasil · Universidade Federal Fluminense (UFF), Niterói, Brasil
In this paper, we introduce INEUS, a meshfree iterative neural solver for partial integro-differential equations (PIDEs). The method replaces the explicit evaluation of nonlocal jump integrals with single-jump sampling and reformulates PIDE solving as a sequence of recursive regression problems. Like Physics-Informed Neural Networks (PINNs), INEUS learns global solutions over the entire space-time domain, yet it offers a more efficient treatment of nonlocal terms and avoids the computationally expensive differentiation of full PIDE residuals. These features make INEUS particularly well suited for high-dimensional PDEs and PIDEs. Supported by a contraction-based convergence proof for linear PIDEs, our numerical experiments show that INEUS delivers accurate and scalable solutions for various high-dimensional linear and nonlinear examples.
Jean-Loup Dupret, Davide Gallon, Patrick Cheridito
Department of Mathematics and RiskLab, ETH Zurich, Switzerland · Applied Mathematics: Institute for Analysis and Numerics, University of Münster, Germany
Physics-informed neural networks (PINNs) have emerged as a promising mesh-free paradigm for solving partial differential equations, yet adoption in science and engineering is limited by slow training and modest accuracy relative to modern numerical solvers. We introduce the Sequential Correction Algorithm for Learning Efficient PINN (Scale-PINN), a learning strategy that bridges modern physics-informed learning with numerical algorithms. Scale-PINN incorporates the iterative residual-correction principle, a cornerstone of numerical solvers, directly into the loss formulation, marking a paradigm shift in how PINN losses can be conceived and constructed. This integration enables Scale-PINN to achieve unprecedented convergence speed across PDE problems from different physics domain, including reducing training time on a challenging fluid-dynamics problem for state-of-the-art PINN from hours to sub-2 minutes while maintaining superior accuracy, and enabling application to representative problems in aerodynamics and urban science. By uniting the rigor of numerical methods with the flexibility of deep learning, Scale-PINN marks a significant leap toward the practical adoption of PINNs in science and engineering through scalable, physics-informed learning. Codes are available at https://github.com/chiuph/SCALE-PINN.
Pao-Hsiung Chiu, Jian Cheng Wong, Chin Chun Ooi +3
Institute of High Performance Computing, A*STAR, Singapore. · Centre for Frontier AI Research, A*STAR, Singapore. · School of Mechanical Engineering, Tianjin University, China. +1