PaPi: Parallel Picard-Duhamel Integration for Fast Neuro-Spectral Architectures

Arthur Bizzi; Márcio Marques; Kleyton Da Costa; Vitor Balestro; Tiago Novello; João M. Pereira; Daniel Yukimura

cs.LGOct 20, 2026

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 ·

Abstract

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

Model	rMSE	TT	rMSE	TT	rMSE	TT	rMSE	TT	rMSE	TT	rMSE	TT
	2D Layers		Marmousi		3D Layers		2D Burgers		Sine-Gordon		KdV
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 $L^{2}$ 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	$\|\lambda_{1}-1\|$	$\|\lambda_{2}-0.1\|/0.1$	TT (s)
0%	PINN (5,000 steps)	9.9e-2 $\pm$ 4.1e-2	4.9e-1 $\pm$ 2.0e-1	8.2e-2 $\pm$ 3.5e-2	419
	QRes (5,000 steps)	4.0e-2 $\pm$ 2.1e-2	2.0e-1 $\pm$ 1.0e-1	3.7e-2 $\pm$ 1.9e-2	561
	FLS (5,000 steps)	9.9e-2 $\pm$ 4.1e-2	4.9e-1 $\pm$ 2.0e-1	8.1e-2 $\pm$ 3.5e-2	411
	PaPi (500 steps)	1.0e-2 $\pm$ 8.4e-3	4.7e-2 $\pm$ 4.1e-2	9.4e-3 $\pm$ 8.2e-3	0.6
	PaPi (5,000 steps, float64)	3.8e-5 $\pm$ 1.2e-5	1.1e-4 $\pm$ 5.3e-5	3.6e-5 $\pm$ 2.2e-5	66
1%	PINN (5,000 steps)	9.8e-2 $\pm$ 4.1e-2	4.9e-1 $\pm$ 2.0e-1	8.0e-2 $\pm$ 3.5e-2	432

Table 2: Inverse Burgers’ problem: mean

\pm

std over the ten starting pairs of the relative

L^{2}

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	$n_{t}$ ( $T$ )	$\epsilon$	Steps	lr	$K$	Quadrature	Train nodes	Deploy	Params
2D Layers	cosine	$201^{2}$	201 (2)	1	2,000	$10^{-2}$	2	exp. trapezoid	101	8	3.2e5
Marmousi	cosine	$201^{2}$	201 (2)	2	2,500 (2,000)	cosine	6	trapezoid	–	8	3.2e5
3D Layers	cosine	$101^{3}$	111 (1)	10	4,000 (2,000)	$10^{-2}$	2	exp. trapezoid	56	8	6.1e4
2D Burgers	Fourier	$201^{2}$	201 (1)	0.1	1,000 (200)	$5\cdot 10^{-3}$	2	trapezoid	–	8	8.1e5
Sine-Gordon	sine	201	201 (3)	0.1	1,000	$2\cdot 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.

n_{t}

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

L^{2}

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 $\pm$ 8.9e-2	7.8e-1 $\pm$ 1.3e-1	420 $\pm$ 56
	QRes	1.5e-1 $\pm$ 7.2e-2	1.9e-1 $\pm$ 7.5e-2	765 $\pm$ 212
	FLS	4.9e-1 $\pm$ 5.2e-2	7.2e-1 $\pm$ 8.3e-2	408 $\pm$ 60
	NeuSA	6.9e-2 $\pm$ 2.0e-2	9.4e-2 $\pm$ 3.1e-2	508 $\pm$ 19
	PaPi	6.8e-2 $\pm$ 1.8e-2	9.1e-2 $\pm$ 2.7e-2	11.8 $\pm$ 0.0
Marmousi	PINN	6.2e-1 $\pm$ 3.5e-2	9.2e-1 $\pm$ 8.1e-2	427 $\pm$ 68

Appendix

Table 4: Main benchmarks: mean

\pm

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.

Problem	rMSE	steps	TT (s)	rMSE	steps	TT (s)	rMSE	steps	TT (s)	torch/JAX	RK4/PaPi
	NeuSA (PyTorch)			RK4 (JAX)			PaPi (JAX)			speedup / step
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\times 10^{17}$	5.025	$8.63\times 10^{-4}$	3.50
201	$7.04\times 10^{17}$	5.021	$2.16\times 10^{-4}$	5.10
401	$6.88\times 10^{17}$	5.018	$5.39\times 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 $\Delta 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 $\eta_{0}=0.01$ . (c) Polyak SGD step-size cap. (d) All configurations against training time; the star marks Adam at $\eta_{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

Neural PDE solvers use trainable models to approximate solutions or solution operators for partial differential equations, often with physical constraints.

See all 269 papers

NEXT: Physics-Informed Neuro-Spectral Exponential Time Differencing Architectures
Sep 25, 2026Márcio Marques, Leonardo Mendonça, Leonardo M. Moreira +6Neural PDE Solvers PDE Solving
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
INEUS: Iterative Neural Solver for High-Dimensional PIDEs
May 7, 2026Jean-Loup Dupret, Davide Gallon, Patrick CheriditoNeural PDE Solvers PDE Solving
Department of Mathematics and RiskLab, ETH Zurich, Switzerland · Applied Mathematics: Institute for Analysis and Numerics, University of Münster, Germany
Scale-PINN: Learning Efficient Physics-Informed Neural Networks Through Sequential Correction
Feb 23, 2026Pao-Hsiung Chiu, Jian Cheng Wong, Chin Chun Ooi +3Neural PDE Solvers Physics-Informed ML
Institute of High Performance Computing, A*STAR, Singapore. · Centre for Frontier AI Research, A*STAR, Singapore. · School of Mechanical Engineering, Tianjin University, China. +1

Sep 25, 2026J/K move · Enter open · S save

NEXT: Physics-Informed Neuro-Spectral Exponential Time Differencing Architectures

Márcio Marques, Leonardo Mendonça, Leonardo M. Moreira, Christian Júnior de Oliveira+5

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

Neural PDE Solvers PDE Solving Spectral Methods PDE Inverse Problems

Open paper AskPDF arXiv

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 .

Abstract

Neural PDE Solvers

Explore similar work

NEXT: Physics-Informed Neuro-Spectral Exponential Time Differencing Architectures

INEUS: Iterative Neural Solver for High-Dimensional PIDEs

Scale-PINN: Learning Efficient Physics-Informed Neural Networks Through Sequential Correction