Solving Inverse Parametrized Problems via Finite Elements and Extreme Learning Networks

Erik Burman; Mats G. Larson; Karl Larsson; Jonatan Vallin

math.NAFeb 16, 2026

Solving Inverse Parametrized Problems via Finite Elements and Extreme Learning Networks

Authors: Erik Burman, Mats G. Larson, Karl Larsson, Jonatan Vallin

Abstract

We develop an interpolation-based modeling framework for parameter-dependent partial differential equations arising in control, inverse problems, and uncertainty quantification. The solution is discretized in the physical domain using finite element methods, while the dependence on a finite-dimensional parameter is approximated separately. We establish existence, uniqueness, and regularity of the parametric solution and derive rigorous error estimates that explicitly quantify the interplay between spatial discretization and parameter approximation. In low-dimensional parameter spaces, classical interpolation schemes yield algebraic convergence rates based on Sobolev regularity in the parameter variable. In higher-dimensional parameter spaces, we replace classical interpolation by extreme learning machine (ELM) surrogates and obtain error bounds under explicit approximation and stability assumptions. The proposed framework is applied to inverse problems in quantitative photoacoustic tomography, where we derive potential and parameter reconstruction error estimates and demonstrate substantial computational savings compared to standard approaches, without sacrificing accuracy.

Figures & tables

Panel of Figure 1 — Figure 1: Domains. (a) The spatial domain $\Omega_{x}$ is a bounded subset in $\mathbb{R}^{d}$ for $d\in\{2,3\}$ . (b) The parameter domain $\Omega_{t}$ is the hypercube $[0,1]^{N_{t}}$ where $N_{t}$ may be large.

Figure 2: Spatial domain discretization . The spatial domain

\Omega_{x}

is discretized using a finite element mesh

\mathcal{T}_{h_{x}}

of mesh size

h_{x}

Figure 3: Low-dimensional interpolation in $\Omega_{t}$ . For $N_{t}\leq 4$ , the parameter domain can be partitioned into simplices with interpolation points at the vertices. This partition constitutes a mesh $\mathcal{T}_{h_{t}}$ of the parameter domain with mesh size $h_{t}$ on which we may employ standard piecewise-linear interpolation.

Panel of Figure 4 — Figure 4: High-dimensional interpolation in $\Omega_{t}$ . For $N_{t}>4$ , an ELM surrogate is used for interpolation. (a) $J$ quasi-random (Sobol’) points in $\Omega_{t}$ are used as interpolation points. (b) Random ReLU feature hyperplanes for the ELM surrogate in $\Omega_{t}$ ; each hyperplane marks where one unit switches between active and inactive regions.

Panel of Figure 5 — Figure 5: Convergence in low-dimensional setting. (a) The heatmap shows approximation error for different values of the pair $(h_{x},h_{t})$ . (b) The figure shows the approximation error as a function of $h_{x}$ for different choices of $h_{t}$ including a reference line to indicate the theoretical convergence rate of $2.0$ with respect to $h_{x}$ . (c) The figure shows the approximation error as a function of $h_{t}$ for different choices of $h_{x}$ including a reference line to indicate the theoretical convergence rate of $2.0$ with respect to $h_{t}$ .

$\Omega_{x}$ -grid	$5\times 5$	$10\times 10$	$20\times 20$	$30\times 30$	$40\times 40$	$50\times 50$
$\Omega_{t}$ -grid	error	error	error	error	error	error
$5\times 5$	2.2e–1	1.02e–1	7.94e–2	7.59e–2	7.48e–2	7.43e–2
$10\times 10$	1.93e–1	5.36e–2	2.34e–2	1.88e–2	1.73e–2	1.67e–2
$20\times 20$	1.88e–1	4.52e–2	1.27e–2	7.26e–3	5.49e–3	4.74e–3
$30\times 30$	1.88e–1	4.39e–2	1.12e–2	5.52e–3	3.64e–3	2.81e–3
$100\times 100$	1.87e–1	4.31e–2	1.02e–2	4.46e–3	2.51e–3	1.63e–3

Table 1: Error in low-dimensional setting. The relative

L^{2}

-error for different grids on

\Omega_{x}

and

\Omega_{t}

Panel of Figure 6 — Figure 6: Convergence in high-dimensional setting. (a) The heatmap shows approximation error for the ELMs for different values of the pair $(h_{x},M)$ . (b) The figure shows the approximation error as a function of $h_{x}$ for different choices of $M$ . A reference line is included to indicate the theoretical convergence rate of $2.0$ with respect to $h_{x}$ . The errors for the fully trained networks with $M=25600$ ReLU units are also included. (c) The figure shows the approximation error as a function of $M$ for different choices of $h_{x}$ . A reference line is included to indicate the theoretical convergence rate of $-0.5$ with respect to $M$ . The errors for the fully trained networks with $h_{x}=0.05$ are also included.

$\Omega_{x}$ -grid	$10\times 10$		$15\times 15$		$20\times 20$		$30\times 30$		$40\times 40$
$M$	error	std	error	std	error	std	error	std	error	std
400	0.148	6.2e–3	0.127	7.1e–3	0.124	8.4e–3	0.124	6.7e–3	0.124	7.7e–3
800	0.126	3.3e–3	0.100	4.0e–3	0.096	4.1e–3	0.095	4.4e–3	0.095	4.4e–3
1600	0.113	2.0e–3	0.080	2.8e–3	0.074	3.3e–3	0.072	3.7e–3	0.072	3.9e–3
3200	0.102	1.2e–3	0.063	2.0e–3	0.055	2.2e–3	0.052	2.3e–3	0.051	2.5e–3
6400	0.096	5.6e–4	0.052	1.0e–3	0.042	1.4e–3	0.037	1.6e–3	0.036	1.6e–3

Table 2: Error in high-dimensional setting. The average of the relative

L^{2}

-error and the estimated standard deviation for the ELMs for different

\Omega_{x}

-grids and numbers of ReLU units.

Panel of Figure 7 — Figure 7: Measurement operator. (a) The absorbed energy is given by $\mu u$ depicted in the figure. (b) By averaging over pixels we get the measurement $Q[\mu u]$ which is a piecewise-constant function on $\Omega_{x}$ subordinate to the pixel partition.

$\Omega_{x}$ -grid	$h_{x}$	$N_{t}$	$J$	$M$	$N_{m}$
$30\times 30$	0.069	50	10000	40000	$25\times 25$

Table 3: Experimental details. Parameter choices in the initial potential reconstruction experiment using a measurement operator with a fixed pixel resolution (

N_{m}=25\times 25

) and full coverage of the domain

\Omega_{x}

Panel of Figure 8 — Figure 8: Potential reconstruction . Illustration of how the potential $\mu$ is successively reconstructed at different stages during the optimization.

Panel of Figure 9 — Figure 9: Loss and error during optimization. (a) The loss over the iterations, (b) the parameter reconstruction error (measured in the Euclidean norm) over the iterations and (c) the potential reconstruction error over the iterations.

Panel of Figure 10 — Figure 10: Measurement operators with different resolutions. Four examples of the information over $\Omega_{x}$ seen by the measurement operator $Q$ using different pixel resolutions.

$\Omega_{x}$ -grid	$h_{x}$	$N_{t}$	$J$	$M$	$N_{m}$
$30\times 30$	0.069	20	10000	40000	varied

Table 4: Experimental details (varied resolution). Parameter choices for the potential reconstruction experiment when varying the pixel resolution of the measurement operator.

Panel of Figure 11 — Figure 11: Effect of pixel resolution. Errors for the potential reconstruction when varying the pixel resolution of the measurement operator. (a) The reconstruction error over the iterations for different pixel resolutions. (b) The final reconstruction error (measured in the Euclidean norm) as a function of the number of pixels.

Panel of Figure 12 — Figure 12: Measurement operators with different domain coverage. Four examples of partial pixel coverage for the measurement operator $Q$ , where only the given percentage of the pixels in $\Omega_{x}$ is seen by the operator.

$\Omega_{x}$ -grid	$h_{x}$	$N_{t}$	$J$	$M$	$N_{m}$
$30\times 30$	0.069	20	10000	40000	$25\times 25$

Table 5: Experimental details (varied domain coverage). Parameter choices for the potential reconstruction experiment when varying the domain coverage of the measurement operator

Q

Panel of Figure 13 — Figure 13: Effect of domain coverage. Errors for the potential reconstruction when varying the percentage of the domain seen by the measurement operator $Q$ . (a) The reconstruction error (measured in the Euclidean norm) over the iterations for different levels of coverage. (b) The final reconstruction error as a function of the coverage.

$\Omega_{x}$ -grid	$h_{x}$	$N_{t}$	$J$	$M$	$N_{m}$
$20\times 20$	0.11	20	2000	8000	$15\times 15$

Table 6: Experimental details (varied noise level). Parameter choices for the potential reconstruction experiment when varying the noise amplitude in the measurements.

Figure 14: Effect of noise level. The mean reconstruction error, measured in the Euclidean norm, after 200 iterations as a function of noise amplitude for the original loss functional

L(t)

and the weighted loss functional

\tilde{L}(t)

PDE Surrogate Modeling

PDE surrogate modeling builds fast approximations of partial differential equation solutions for simulation, prediction, and design.

See all 322 papers

Hierarchical Inference and Closure Learning via Adaptive Surrogates for ODEs and PDEs
Mar 4, 2026Pengyu Zhang, Arnaud Vadeboncoeur, Alex Glyn-Davies +1PDE Inverse Problems PDE Surrogate Modeling
Department of Engineering, University of Cambrdige, 7a JJ Thomson Ave, Cambridge, CB3 0FA, UK · The Alan Turing Institute, 96 Euston Rd., London, NW1 2DB, UK
PDEInvBench: A Comprehensive Dataset and Design Space Exploration of Neural Networks for PDE Inverse Problems
May 25, 2026Divyam Goel, Nithin Chalapathi, Sanjeev Raja +1Neural Network Optimization Synthetic Benchmark Generation
Department of Computer Science, UC Berkeley · Departments of Computer Science and Chemical EngineeringMay UC Berkeley; LBNL
Solving Inverse PDE Problems using Minimization Methods and AI
Mar 2, 2026Noura Al Helwani, Sophie Moufawad, Georges SakrPDE Solving PDE Inverse Problems
Department of Physics, American University of Beirut, Lebanon · Department of Mathematics, American University of Beirut, Lebanon · The authors would like to acknowledge

Mar 4, 2026J/K move · Enter open · S save

Hierarchical Inference and Closure Learning via Adaptive Surrogates for ODEs and PDEs

Pengyu Zhang, Arnaud Vadeboncoeur, Alex Glyn-Davies, Mark Girolami

Department of Engineering, University of Cambrdige, 7a JJ Thomson Ave, Cambridge, CB3 0FA, UK · The Alan Turing Institute, 96 Euston Rd., London, NW1 2DB, UK

PDE Inverse Problems PDE Surrogate Modeling

Open paper AskPDF arXiv

Inverse problems are the task of calibrating models to match data. They play a pivotal role in diverse engineering applications by allowing practitioners to align models with reality. In many applications, engineers and scientists do not have a complete picture of i) the detailed properties of a system (such as material properties, geometry, initial conditions, etc.); ii) the complete laws describing all dynamics at play (such as friction laws, complicated damping phenomena, and general nonlinear interactions). In this paper, we develop a principled methodology for leveraging data from collections of distinct yet related physical systems to jointly estimate the individual model parameters of each system, and learn the shared unknown dynamics in the form of an ML-based closure model. To robustly infer the unknown parameters for each system, we employ a hierarchical Bayesian framework, which allows for the joint inference of multiple systems and their population-level statistics. To learn the closures, we use a maximum marginal likelihood estimate of a neural network embeded within the ODE/PDE formulation of the problem. To realize this framework we utilize the ensemble Metropolis-Adjusted Langevin Algorithm (MALA) for stable and efficient sampling. To mitigate the computational bottleneck of repetitive forward evaluations in solving inverse problems, we introduce a bilevel optimization strategy to simultaneously train a surrogate forward model alongside the inference. Within this framework, we evaluate and compare distinct surrogate architectures, specifically Fourier Neural Operators (FNO) and parametric Physics-Informed Neural Network (PINNs).

Abstract

PDE Surrogate Modeling

Explore similar work

Hierarchical Inference and Closure Learning via Adaptive Surrogates for ODEs and PDEs

PDEInvBench: A Comprehensive Dataset and Design Space Exploration of Neural Networks for PDE Inverse Problems

Solving Inverse PDE Problems using Minimization Methods and AI