Inverse Problems and Parameter Estimation for PDE Models

An inverse problem asks a simple but fundamentally difficult question: given measurements or simulation outputs, what are the hidden parameters that produced them? In forward modelling, you specify initial and boundary conditions, material properties, and source terms, then integrate the governing equations to predict the system’s behaviour. In an inverse problem, you flip the process. You observe the outcome — a temperature field, a pressure profile, a displacement — and try to recover the unknown inputs that drove it.

This is parameter estimation for PDE models: the mathematical machinery for turning measurements into parameter estimates, and the computational tools that make it practical. It appears everywhere from geophysics (estimating subsurface permeability from well measurements) to materials science (recovering diffusion coefficients from observed concentration fields) to physics-informed machine learning (discovering governing equations from data).

The challenge is that inverse problems are almost always ill-posed. Small errors in your data can produce wildly wrong estimates. Without careful mathematical formulation and numerical regularization, your best-fitting parameters may be completely misleading.

Key Takeaways

• Every inverse problem starts with a forward model. The parameter-to-solution map $y = \mathcal{G}(\theta)$ defines how parameters $\theta$ produce observables $y$. Inverse problems recover $\theta$ from noisy measurements.
• Discretization matters more than you’d think. The order in which you discretize the forward model, impose regularization, and solve the inverse problem (DTI vs CTI vs ITD) determines whether your results depend on your mesh or reflect physical reality.
• Regularization tames ill-posedness. Tikhonov regularization, L-curve analysis, and the discrepancy principle give you practical ways to choose regularization parameters without overfitting.
• The adjoint method computes gradients efficiently. Instead of computing expensive finite-difference approximations for every parameter, adjoint equations solve for gradient information in a single backward sweep — critical when you have thousands of parameters.
• Bayesian inversion gives you distributions, not just point estimates. The posterior distribution $\mu^y$ characterizes the full parameter uncertainty conditional on data, and function-space formulations ensure discretization-invariant results.

What is an Inverse Problem?

Every inverse problem pairs a forward model with a dataset and asks a reconstruction question. The forward model is a mathematical object that maps parameters to observables. In PDE contexts, this is usually a solution operator: given a set of material properties, boundary conditions, or source terms, solve the governing equations and extract the quantities you care about.

For a linear PDE, the forward problem looks like

$$Au = f, \quad y = \mathcal{C}(u),$$

where $A$ is a differential operator encoding physics (diffusion, advection, elasticity), $u$ is the state variable, $f$ contains source terms and boundary conditions, and $\mathcal{C}$ is an observation operator that extracts measurements at specific locations or times.

The inverse problem inverts this relationship. You’re given noisy data $y^\delta$ and you want to find the parameters $\theta$ that produced it:

$$y^\delta = \mathcal{G}(\theta^\dagger) + \eta, \qquad \eta \sim \mathcal{N}(0, \Gamma),$$

where $\mathcal{G}$ is the observation operator that applies the forward model then the observation map, $\theta^\dagger$ is the true parameter, and $\eta$ represents measurement noise.

The fundamental issue is well-posedness. Hadamar’d well-posedness requires three properties: existence of a solution, uniqueness, and continuous dependence on data. Inverse problems often violate the third property. If you perturb your measurements by a tiny amount $\delta$, the corresponding change in your estimated parameters might blow up. This is the hallmark of an ill-posed problem.

Pro tip: Think of ill-posedness as the mathematical expression of “information loss.” When you observe a temperature distribution, you’ve lost information about the initial state, boundary conditions, or material properties. Reconstructing those hidden quantities is mathematically equivalent to running time backwards, which amplifies any noise present.

Why PDE Inverse Problems Are Hard

When your forward model is a PDE, the situation is more constrained than in simple inverse problems like deconvolution. The forward operator $\mathcal{G}$ is typically a compact operator on infinite-dimensional function spaces, which means it has eigenvalues that decay to zero. This decay is exactly what makes the problem ill-posed — small singular values mean that noise in those directions gets amplified without bound.

This distinguishes PDE-driven inverse problems from standard statistical estimation. Your parameters live in function space (spatially varying fields), not in $\mathbb{R}^n$, and your prior information about them should reflect that structure. The field of Bayesian inversion in function space — pioneered by Franka (1970), Mandelbaum (1984), and formalized by Lassas and Siltanen (2004) — provides the rigorous framework for handling this.

The key insight from function-space Bayesian inversion is that you should formulate your prior and posterior measures on infinite-dimensional spaces, then discretize them consistently. If you discretize first and then do Bayes, your posterior may depend on the mesh — a sign that you haven’t respected the continuum nature of the problem.

See also Verification vs Validation in Scientific Simulations for context on how parameter estimation relates to the broader V&V pipeline, and Uncertainty Quantification and Sensitivity Analysis in Scientific Simulation for the complementary topic of propagating parameter uncertainty through forward models.

Three Discretization Strategies

When you implement an inverse problem numerically, the order in which you discretize matters fundamentally. There are three established strategies, each with different trade-offs between theoretical cleanliness and practical complexity.

DTI (Discretize-Then-Inform) is the most common in practice. You take your standard PDE solver, discretize the forward model on a mesh, and then solve the resulting discrete inverse problem with standard regularization methods. This is straightforward but has a hidden cost: your solution will depend on the mesh resolution. If you refine the mesh, your regularized solution may change significantly. This is theoretically unsatisfying and practically risky.

CTI (Continuous-Then-Inform) inverts the order. You formulate the entire inverse problem — including prior distributions, likelihood, and regularization — in the infinite-dimensional setting. Then you discretize consistently. For example, you might place a Gaussian process prior on the infinite-dimensional parameter space, specify a likelihood on the observation space, and derive the posterior in function space. Only then do you discretize for computation. The advantage is that your posterior doesn’t depend on an arbitrary mesh — it reflects the continuum problem.

The seminal work by Lassas, Saksman, and Siltanen (2009) on discretization-invariant Bayesian inversion established that CTI approaches yield posterior measures that converge to a mesh-independent limit. Their framework uses Besov space priors that encode prior knowledge about regularity, and it shows that as long as your discrete spaces approximate the continuous ones, your posterior estimates converge.

ITD (Inform-Then-Discretize) is the most sophisticated but also the most flexible. Here you reformulate the problem itself — perhaps by changing the forward operator, modifying boundary conditions, or adding auxiliary constraints — to improve its well-posedness before discretizing. For example, you might replace a second-order PDE with a first-order system (Redheffer-Wign er form), which is better suited for certain inversion algorithms. This strategy is less commonly taught in textbooks but appears frequently in modern scientific computing.

⚠️ Common mistake: Using DTI and assuming your solution converges as you refine the mesh. For ill-posed problems, refining the mesh without tightening regularization will eventually diverge — your solution will oscillate more and more around the data. You need to show that your regularized sequence converges as $\delta \to 0$ and regularization parameter $\alpha \to 0$ simultaneously.

Regularization: Taming Ill-Posed Problems

The most common approach to making an ill-posed problem solvable is Tikhonov regularization. Instead of solving $Au = f$ directly, you solve a modified problem that penalizes “complex” solutions:

$$\min_{u} |Au – f|_{\Gamma}^2 + \alpha |Lu|^2, $$

where $L$ is a regularization operator (often a discrete approximation of the gradient or Laplacian), $\alpha$ is the regularization parameter, and $\Gamma$ is the noise covariance. The first term is the data misfit (how well your model explains the observations). The second term is the regularization penalty (how “smooth” or “simple” your solution is).

The regularization parameter $\alpha$ is the critical knob. Too large, and you oversmooth — your solution is close to the prior but ignores the data. Too small, and you under-regularize — you fit noise and get unstable solutions. Finding the right value is itself an inverse problem.

The L-Curve

The L-curve is the most widely used heuristic for choosing $\alpha$. You plot the log of the data misfit $|Au – f|$ against the log of the regularized solution norm $|Lu|$ as $\alpha$ varies. The resulting curve typically has an “L” shape: a vertical segment where regularization dominates, a horizontal segment where data misfit dominates, and a corner region where the two balance.

Pro tip: The corner of the L-curve is usually where $\alpha$ should sit, but not always. For problems with correlated noise (e.g., spatially smooth measurement errors), you may want to shift slightly toward the horizontal segment — letting the solution be slightly rougher is often better than oversmoothing.

The L-curve is easy to compute for small-scale problems but becomes expensive for large-scale PDE inverse problems. Each point on the curve requires solving the forward problem with a different regularization parameter. With hundreds of parameters, that’s impractical.

The Discrepancy Principle

A more theoretically grounded approach is the discrepancy principle (Moroz, 1968). If you know or can estimate the noise level $\delta = |\eta|$, the discrepancy principle selects $\alpha$ such that:

$$|Au^\alpha – f| = \eta \delta, \quad \eta > 1.$$

The safety factor $\eta > 1$ prevents overfitting: you stop fitting as soon as the residual reaches the estimated noise level. Anything beyond that is fitting noise, not signal. The principle is easy to implement — you solve for $\alpha$ by bisection or Newton’s method — and it has provable convergence properties.

See Time Integration Methods for PDE Solvers for context on how numerical stability relates to regularization choices, and Stiff Problems and Stiff Solvers when your regularized system exhibits stiffness.

The Adjoint Method: Computing Gradients Efficiently

Suppose your forward model maps 1,000 material parameters $\theta_1, \ldots, \theta_{1000}$ to observables. You want to minimize a misfit function — perhaps the sum of squared differences between observed and predicted data — and you need gradients $\partial F / \partial \theta_i$ for every $i$.

A naive approach uses finite differences: perturb $\theta_i$ by $\epsilon$, solve the forward problem again, and estimate the gradient. That’s 1,000 forward solves per gradient evaluation. For a PDE solver that takes minutes per solve, that’s hours just for one gradient step.

The adjoint method solves this in two sweeps regardless of the number of parameters. It exploits the fact that the gradient of a scalar loss function with respect to parameters can be computed by solving one additional PDE — the adjoint equation.

The adjoint method dates back to the work of Pironneau (1974) and Kontoleadis et al. (2013), and was formalized for PDE-constrained optimization by Gunther et al. (2013) in the context of dolfin-adjoint. The basic idea:

1. Solve the forward problem $A(\theta)u = f$ to get the state $u$.
2. Solve the adjoint equation $A(\theta)^* p = \frac{\partial R}{\partial u}$ to get the adjoint state $p$, where $R$ is your cost function and $A^*$ is the adjoint operator.
3. Compute the gradient $\nabla_\theta F = \frac{\partial F}{\partial \theta} + \text{terms involving } p \text{ and } u$.

The adjoint equation is essentially the forward model run backwards in time (if the PDE is time-dependent) with the cost function’s gradient as a source term. For a parabolic PDE, you solve from the final time back to the initial time, picking up contributions along the way.

The computational cost of the adjoint method is essentially one extra forward solve — not one per parameter. This makes it practical for thousands of parameters, where finite differences would be impossible.

Givoli (2021) provides an excellent pedagogical tutorial on the adjoint method for finite element problems, showing how the discrete adjoint is derived from the finite element formulation and how it relates to the continuous adjoint. The key insight is that the discrete adjoint of a discretized forward model is not necessarily the discretization of the continuous adjoint — this is where the discretization strategies from Section 2 become relevant again.

Pro tip: If you’re using FEniCS for your forward model, dolfin-adjoint is the canonical choice for automatic adjoint derivation. It automatically constructs the discrete adjoint from your forward code by differentiating the assembled system. The setup is simple — wrap your forward solve with DOLFIN_Adjoint(), specify your cost function, and call compute_gradient() — and it handles the adjoint derivation automatically.

For a hands-on introduction, the Stanford AMBAD adjoint tutorial (Bradić et al.) walks through the derivation of adjoint equations for Navier-Stokes and advection-diffusion problems, with clear finite element formulations. It’s an excellent practical reference.

Bayesian Inversion: Uncertainty-Aware Parameter Estimation

Point estimates — whether from Tikhonov regularization or adjoint-based optimization — give you a single “best” set of parameters. But they don’t tell you how confident you should be. In many applications, you need the full uncertainty distribution: What’s the 95% credible interval on the permeability? How likely is it that the true parameter lies outside the estimated range?

Bayesian inversion provides this by treating parameters as random variables. You specify a prior distribution $\mu_0$ encoding your prior knowledge about the parameters, then compute the posterior distribution $\mu^y$ using Bayes’ rule:

$$\mu^y(\theta) \propto \exp\left(-\frac{1}{2} | y – \mathcal{G}(\theta) |_{\Gamma}^2\right) \mu_0(\theta). $$

The exponential term is the likelihood (how likely the data are given the parameters), and $\mu_0$ is the prior. The posterior combines both: it balances the data against prior knowledge.

The Challenge: MCMC

Computing the posterior is tractable only in simple linear-Gaussian cases. For nonlinear PDEs, you need Markov Chain Monte Carlo (MCMC) methods to sample from $\mu^y$. The problem is that each MCMC proposal requires solving the forward model. With 10,000 proposals and a PDE solve that takes 5 minutes, you’re looking at days or weeks of computation.

This is why surrogate models are essential in Bayesian inverse problems. Instead of solving the PDE at every MCMC step, you build a fast approximation of the forward map. Common surrogates include:

• Gaussian Process (GP) regression — provides both predictions and uncertainty estimates. Effective for low-dimensional parameter spaces but computationally expensive for high-dimensional parameters ($\mathcal{O}(N^3)$ for $N$ training points).
• Neural network surrogates — faster evaluation but harder to quantify uncertainty. The DeepGaLA framework (Jimenez-Beltran et al., 2024) addresses this by combining the Deep Galerkin Method with Laplace Approximation, providing calibrated uncertainty estimates at the network’s last layer.

Discretization-Invariant Inversion

A critical insight from function-space Bayesian inversion (Lassas et al., 2009) is that your posterior should not depend on the mesh. In the DTI approach, you might place a prior on a discretized parameter vector, then refine the mesh and discover that your posterior has shifted. This is because the prior itself depends on the discretization.

The CTI approach avoids this by placing priors on infinite-dimensional function spaces (typically Gaussian processes or Besov spaces), specifying the likelihood on the observation space, and discretizing the posterior measure consistently. As the mesh refines, the discrete posterior converges to a mesh-independent limit.

For a practical walkthrough of Bayesian inverse problems in function space, Dan MacKinlay’s notebook provides an accessible introduction that bridges the theory and implementation, showing how measurement discretization and computational discretization interact in the Bayesian framework.

Python Implementation: TRIPs-Py and dolfin-adjoint

Two Python packages make practical inverse problems accessible. TRIPs-Py (Pasha et al., 2024) focuses on linear discrete inverse problems, while dolfin-adjoint handles PDE-constrained optimization with automatic adjoint derivation.

TRIPs-Py: Regularization for Linear Inverse Problems

TRIPs-Py is designed for linear inverse problems of the form $Ax = b$. It provides:

• Direct methods: TSVD, Tikhonov, truncated GSVD
• Iterative methods: GMRES, LSQR, CGLS, hybrid variants
• Regularization parameter selection: L-curve, discrepancy principle, GCV
• Test problems: 1D/2D deblurring, tomography

from trips_py import Deblurring1D, TSVD, Tikhonov, GCV

# Setup
deblurr = Deblurring1D()
nx = 200
x_true = deblurr.gen_xtrue(nx, test='curve2')
A = deblurr.forward_Op_1D(parameter=30, nx=nx)
b = deblurr.gen_data(x_true)
b_noisy, delta = deblurr.add_noise(b, 'Gaussian', noise_level=0.01)

# Solve with Tikhonov regularization
solver = Tikhonov()
x_reg = solver.solve(A, b_noisy, regularization_parameter=0.01)

# Or use GCV to find the regularization parameter automatically
solver = Tikhonov(regularization='auto')
solver.set_reg_param_method('GCV')
x_reg = solver.solve(A, b_noisy)

The CommitCrime option in TRIPs-Py is worth understanding. By default, TRIPs-Py avoids the inverse crime by using mismatched forward operators for data generation and solution. This is a practical safeguard: if your forward operator for solving exactly matches the one used to generate the data, your regularization may appear more effective than it actually is. Set CommitCrime=True only if you specifically want to study the inverse crime phenomenon.

dolfin-adjoint: Automatic Adjoint for PDE Systems

dolfin-adjoint integrates with FEniCS to provide automatic adjoint derivation for PDE-constrained optimization:

from fenics import *
from dolfin_adjoint import *

# Define forward model
mesh = RectangleMesh(50, 50, 1.0, 1.0)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
alpha = Function(V)  # Unknown parameter field

# Forward PDE: -∇·(α∇u) = f
f = Constant(1.0)
solve(Inner(grad(u), grad(v))*dx == f*v*dx, alpha, 
      adjoint=adj_solve(Inner(grad(u), grad(v))*dx == f*v*dx))

# Cost function
J = 0.5 * sum(v*v for v in V.sub(0).split())

# Compute gradient via adjoint method
compute_gradient(J, alpha)

The adj_solve function automatically differentiates the forward solve, generating the adjoint equation internally. You don’t need to derive or implement the adjoint manually — dolfin-adjoint handles it. This is particularly valuable when your PDE has complex boundary conditions or nonlinear operators.

For problems where you need uncertainty quantification alongside parameter estimation, the combination of dolfin-adjoint (for efficient gradient computation) with Bayesian samplers (like emcee or pymc) provides a powerful framework. The gradient information from dolfin-adjoint accelerates gradient-based MCMC methods like the No-U-Turn Sampler (NUTS).

The “Inverse Crime”: Why Your Synthetic Data May Be Lying to You

The inverse crime is a notorious pitfall in inverse problem research. It occurs when the forward operator used to generate synthetic test data is identical to the one used to solve the inverse problem. The result? Your regularization appears dramatically more effective than it actually is.

Here’s why: regularization methods like Tikhonov implicitly penalize certain solution features. If your synthetic data was generated with an operator that happens to align with your regularization prior, you’ll recover the ground truth effortlessly — not because your method is good, but because you set yourself up to succeed.

The consequence is that published benchmarks based on inverse crime data can be misleading. A method that achieves 99% accuracy on inverse crime data might drop to 60% on real-world data where the forward operator is genuinely mismatched (due to discretization differences, model simplifications, or measurement geometry).

⚠️ Common mistake: Generating test data with the same forward model you use for inversion, then claiming your method is “robust” based on the excellent recovery. This is the inverse crime — it inflates performance metrics by exploiting structural alignment between data generation and inversion.

TRIPs-Py addresses this by default: its Deblurring1D class uses different boundary conditions for data generation versus solution, avoiding the inverse crime. To study the inverse crime specifically, you must explicitly set CommitCrime=True. This default protection makes TRIPs-Py uniquely useful for honest benchmarking.

The mathematical foundation of the inverse crime was clarified by Kaipio and Sereno (2007), who showed that the inverse crime can be quantified as the alignment between the null spaces of the forward operator and the regularization operator. When these null spaces align, regularization becomes artificially effective. When they don’t, the true difficulty of the inverse problem emerges.

Related: PINNs and Machine Learning for Inverse Problems

Machine learning has entered the inverse problem space in two major ways: Physics-Informed Neural Networks (PINNs) and neural operator surrogates.

PINNs (Raissi et al., 2019) embed the PDE directly into the neural network’s loss function. The network learns to satisfy both the data and the physics simultaneously. For inverse problems, PINNs can discover unknown parameters by minimizing the residual of the governing equations alongside the data misfit. The beauty is that you don’t need to know the parameters to set up the network — they’re treated as trainable constants.

However, PINNs face challenges:

• Spectral bias: Neural networks learn low-frequency components first, which can cause them to miss fine-scale parameter variations.
• Loss balancing: Different terms in the loss (PDE residual, boundary conditions, data misfit) have different scales. Getting the right balance requires tuning or adaptive weighting.
• Convergence guarantees: Unlike Tikhonov regularization, there are no proven convergence results for PINNs in the noise-to-signal limit.

Neural operator surrogates (like DeepGaLA) are another emerging approach. Instead of learning the solution map directly, they learn a fast approximation of the forward model that can be evaluated millions of times during MCMC sampling. DeepGaLA (Jimenez-Beltran et al., 2024) provides a concrete framework: train a neural network surrogate with the Deep Galerkin Method, then use Laplace Approximation to quantify uncertainty at the last layer. The resulting random surrogate enables scalable Bayesian inverse problems without sacrificing reliability.

For deeper context on PINNs, see Physics-Informed Neural Networks (PINNs) for Scientific Simulations.

Summary + Next Steps

Inverse problems for PDE models sit at the intersection of three disciplines: numerical analysis (how to solve them efficiently), applied mathematics (how to guarantee convergence), and statistics (how to quantify uncertainty). The practical challenge is choosing the right combination of methods for your specific problem.

Here’s a practical workflow:

1. Start with a forward model. Implement or validate your PDE solver. If you haven’t done this, see Verification vs Validation in Scientific Simulations.
2. Formulate the inverse problem. Specify what parameters you’re recovering, what data you have, and what noise level to expect.
3. Choose a discretization strategy. DTI is simplest; CTI gives mesh independence. If you need both, ITD offers the most flexibility.
4. Pick a regularization method. Tikhonov with L-curve or discrepancy principle is the standard starting point. TRIPs-Py provides implementations.
5. Compute gradients efficiently. If you have many parameters, use the adjoint method. dolfin-adjoint handles the derivation automatically.
6. Quantify uncertainty. For high-dimensional parameters, consider neural surrogates like DeepGaLA. For lower dimensions, MCMC with a GP surrogate is effective.
7. Avoid the inverse crime. Always check that your forward operator for data generation differs from your inversion operator, or explicitly study the crime.

The field is evolving rapidly. Deep learning approaches are complementing traditional regularization methods, and automatic differentiation tools like dolfin-adjoint are making adjoint-based inversion accessible to practitioners who don’t want to derive adjoint equations manually. The key is grounding your choice in the mathematical properties of your specific problem — knowing when a problem is ill-posed, how regularization tames it, and what uncertainty quantification tells you about your parameters.

Related Guides

• Verification vs Validation in Scientific Simulations — Understanding how parameter estimation fits into the broader validation framework
• Uncertainty Quantification and Sensitivity Analysis in Scientific Simulation — Propagating parameter uncertainty through forward models
• Time Integration Methods for PDE Solvers — How numerical time stepping affects forward solve accuracy
• Stiff Problems and Stiff Solvers — When regularization introduces stiffness into your forward solves
• Physics-Informed Neural Networks (PINNs) for Scientific Simulations — ML approaches that embed PDEs directly into learning