CUQIpy: Bayesian UQ in Inverse Problems

• CUQIpy is an open-source Python package for computational uncertainty quantification that integrates Bayesian frameworks with inverse problem modeling under standard and PDE-based scenarios.
• It features a modular design with high-level syntax, automatic sampler selection, and extensibility through plugins like FEniCS for finite element analysis.
• The framework supports both explicit and implicit priors, enabling advanced posterior inference, diagnostics, and visualization across various applications such as imaging, EIT, and PAT.

Searching arXiv for CUQIpy papers and related framework papers to ground the article. CUQIpy is an open-source Python package for computational uncertainty quantification in inverse problems under a Bayesian framework. It was introduced in a two-part series in which Part I presents a general framework for inverse problems in Python, while Part II extends the package to PDE-based Bayesian inverse problems through a general framework that allows the integration of PDEs expressed natively or through third-party libraries such as FEniCS. A subsequent paper introduces a computational framework and implementation of implicit priors within the same package. Across these works, CUQIpy is characterized by high-level modeling syntax that closely matches mathematical expressions, a modular architecture, efficient sampling strategies, automatic sampler selection, visualization tools, and extensibility through plugins (Riis et al., 2023, Alghamdi et al., 2023, Everink et al., 15 Sep 2025).

1. Scope, objectives, and positioning

CUQIpy is designed for computational UQ in inverse problems, including small and large-scale imaging inverse problems and inverse problems governed by partial differential equations. Its stated objectives combine accessibility for non-experts with advanced functionality for expert users. The package is intended to let users specify forward models, prior information, likelihood assumptions, and posterior inference workflows with concise syntax that mirrors the underlying mathematical formulation (Riis et al., 2023).

Part I emphasizes a general inverse-problem framework with built-in probability distributions, test problems, computational methods, visualization, and automatic sampler selection. Part II specializes this framework to PDE-based Bayesian inverse problems and introduces a generic API for integrating PDE models as native implementations, black-box solvers, or third-party finite-element software. The later implicit-priors paper broadens the package beyond explicit density-based priors by providing a unified computational framework for priors encoded through denoisers, proximal operators, optimization constraints, and transformations (Alghamdi et al., 2023, Everink et al., 15 Sep 2025).

A central theme across the package is the separation of statistical modeling and UQ logic from the numerical realization of the forward problem. This design supports a uniform user interface across linear and nonlinear inverse problems, hierarchical models, PDE-constrained models, and workflows that rely on explicit or implicit prior information. This suggests that CUQIpy is intended not merely as a collection of samplers, but as a modeling environment in which posterior inference, diagnostics, and visualization are coupled to the structure of inverse problems.

2. Bayesian formulation and modeling abstractions

CUQIpy adopts a fully Bayesian formulation. In the generic setting, an inverse problem is written as

$A(x) = y,$

or, in the noisy observational form used for PDE-based inverse problems,

$y = \mathcal{A}(x) + \epsilon.$

Here, $x$ denotes the unknown parameter, $y$ the observed data, $\mathcal{A}$ the forward or parameter-to-observation map, and $\epsilon$ the observational noise. The statistical model is then specified through a prior and likelihood,

$x \sim p(x), \qquad y \sim p(y \mid x),$

with posterior

$$p(x \mid y_{\mathrm{obs}}) \propto p(y_{\mathrm{obs}} \mid x)\,p(x).$$

For hierarchical models, CUQIpy also supports hyperpriors, for example

$$d \sim p(d), \qquad s \sim p(s), \qquad x \sim p(x \mid d), \qquad y \sim p(y \mid x, s),$$

with the corresponding joint posterior over state and hyperparameters (Riis et al., 2023, Alghamdi et al., 2023).

The software organizes this formulation through object-oriented abstractions including Distribution, Model, Geometry, Sampler, Samples, and BayesianProblem. The BayesianProblem class is used to formulate the statistical model, attach observed data, and invoke posterior UQ. The design is geometry-aware: models and distributions are associated with domain and range contexts that also guide plotting and post-processing. Supported geometries include Continuous1D/2D, Discrete, Image2D, and expansion-based representations such as KL, step, and mapped geometries (Riis et al., 2023).

The package includes a broad collection of explicit distributions. Examples listed in the source materials include Gaussian, Gamma, Beta, Cauchy, Laplace, Uniform, InverseGamma, and Lognormal, as well as Markov random field priors such as GMRF, LMRF, and CMRF. The latter support smoothness, edge-preserving, and heavy-tailed regularization effects. Because these distributions enter the same modeling interface as forward models and likelihoods, CUQIpy can express both standard and hierarchical Bayesian inverse problems in a uniform syntax (Riis et al., 2023).

3. Architecture, extensibility, and PDE integration

A defining architectural principle of CUQIpy is the decoupling of Bayesian modeling, probabilistic distributions, and sampling methods from the PDE modeling layer. Part II introduces a general framework for PDE integration based on an abstract PDE class and a PDEModel wrapper. The abstract PDE interface is organized around methods to assemble, solve, and observe, while PDEModel couples a PDE object with specific domain and range geometries to produce the forward operator $\mathcal{A}$ (Alghamdi et al., 2023).

This framework is explicitly intended to support three styles of PDE integration: native CUQIpy PDE descriptions, user-provided black-box solvers, and third-party libraries such as FEniCS. The plugin architecture is central to this design. The source materials mention plugins for PyTorch, CIL, and FEniCS, and specifically identify a CUQIpy-FEniCS plugin that enables PDE definition through FEniCS’s Unified Form Language while retaining a uniform UQ interface (Riis et al., 2023, Alghamdi et al., 2023).

For finite-element workflows, the PDE extension introduces geometry classes tailored to FEM contexts, including FEniCSContinuous, which encapsulates FEniCS function spaces for parameters and data, and MaternKLExpansion, which implements KL-based smoothness priors using eigenmodes and eigenvalues derived from Matérn covariances on the FEM space. More generally, geometry classes provide parameter-space representations such as StepExpansion, KL expansions, FEM spaces, and mapped transformations, allowing priors and constraints to be expressed at the level of coefficients, fields, or latent variables (Alghamdi et al., 2023, Everink et al., 15 Sep 2025).

At the array level, CUQIpy is described as array-agnostic. Its core data structures permit fallback to NumPy and can also use PyTorch for automatic differentiation and GPU acceleration. A plausible implication is that this design allows the same high-level modeling code to be paired with different computational backends, depending on whether the primary requirement is standard numerical computing, automatic differentiation, or acceleration hardware (Riis et al., 2023).

4. Inference workflow, sampler selection, and diagnostics

CUQIpy supports a broad family of posterior sampling methods. The listed samplers include Metropolis-Hastings, preconditioned Crank-Nicolson, ULA, MALA, NUTS, LinearRTO, Gibbs and conjugate samplers, Component-Wise MH, and UGLA. The package also supports direct sampling for Gaussian models and hybrid or conjugate Gibbs strategies when closed-form conditional distributions are available (Riis et al., 2023).

A distinctive feature is automatic sampler selection. After the Bayesian problem is formulated, CUQIpy analyzes the problem structure and chooses a suitable inference method without user intervention. The source materials provide examples of this mapping: linear Gaussian problems can use direct sampling or LinearRTO; hierarchical models with conjugate structure can use Gibbs variants; problems with LMRF or CMRF priors can use UGLA or NUTS; and user-defined nonlinear models can use gradient-based samplers such as NUTS when automatic differentiation is available (Riis et al., 2023).

The modeling syntax is intentionally aligned with mathematical notation. For example, the Gaussian prior-likelihood model

$y = \mathcal{A}(x) + \epsilon.$0

is represented in CUQIpy with corresponding Gaussian objects, a forward model A, and a BayesianProblem object whose UQ() method runs the posterior inference workflow. The same design principle is used in PDE-based examples, including FEniCS-based Poisson models and black-box wave-equation solvers (Alghamdi et al., 2023).

Posterior analysis is supported through integrated diagnostics and visualization. The reported diagnostics include effective sample size, integrated autocorrelation time, and the R-hat statistic via ArviZ. Samples objects provide plotting methods such as .plot_ci(), .plot_mean(), .plot_std(), .plot_trace(), and .plot_violin(). In PDE contexts, posterior means, variances, and credible intervals can also be visualized through FEniCS plotting for FEM-native quantities (Riis et al., 2023, Alghamdi et al., 2023).

5. PDE-based inverse problems and demonstrated use cases

Part II demonstrates the PDE framework on parabolic, elliptic, and hyperbolic inverse problems. The examples include the heat equation, the Poisson equation, electrical impedance tomography, and photo-acoustic tomography. These examples are presented as illustrations of the package’s efficiency, consistency, and intuitive interface in PDE-based Bayesian inverse problems (Alghamdi et al., 2023).

For the one-dimensional heat equation, the package is used to infer an unknown initial temperature profile from later temperature measurements. The PDE is posed as

$y = \mathcal{A}(x) + \epsilon.$1

with homogeneous boundary conditions and initial condition $y = \mathcal{A}(x) + \epsilon.$2. The unknown initial profile can be parameterized by a step expansion,

$y = \mathcal{A}(x) + \epsilon.$3

or by a Karhunen–Loève expansion. The examples use synthetic data with varying noise levels and observation locations, and quantify uncertainty both in the coefficient space and over the physical domain (Alghamdi et al., 2023).

For elliptic problems, the paper includes a two-dimensional Poisson equation example implemented through FEniCS. The forward model is based on

$y = \mathcal{A}(x) + \epsilon.$4

where $y = \mathcal{A}(x) + \epsilon.$5 is the log-conductivity, ensuring positivity of the conductivity field. The conductivity is parameterized through a truncated KL expansion, and posterior inference is reported using NUTS, with posterior mean, variance, and credible intervals shown for both fields and coefficients (Alghamdi et al., 2023).

Electrical impedance tomography is treated as an elliptic inverse problem with multi-frequency boundary current measurements. The data summary describes a multi-layered parameterization that combines a Gaussian random field via KL expansion with a non-linear Heaviside level-set mapping, thereby supporting discontinuous inclusions. The case study emphasizes credible intervals for reconstructed conductivities and their boundaries, and also notes that multiple datasets and model configurations are handled within the same framework (Alghamdi et al., 2023).

Photo-acoustic tomography provides the hyperbolic example. Here the task is to reconstruct an initial pressure from time-dependent ultrasound measurements using the wave equation. The example is notable for its use of black-box or user-supplied forward solvers and for modeling both complete and partial measurement scenarios. Posterior inference is reported with the pCN sampler, and the uncertainty in the reconstructions is described as increasing where data are limited (Alghamdi et al., 2023).

PDE class	Parameterization / implementation	Demonstrated use case
Heat (parabolic)	Finite Difference; Step, KL	Recovering temperature profile
Poisson and EIT (elliptic)	FEM (FEniCS); KL, Step, Level-set (KL + Heaviside)	Conductivity inference; inclusion detection
PAT (hyperbolic)	User-defined black-box; KL, Step	Reconstructing pressure

6. Implicit priors, advanced prior modeling, and limitations

The 2025 paper extends CUQIpy from explicit density-based priors to a broader computational treatment of implicit priors. In this formulation, explicit priors are specified through a probability density function $y = \mathcal{A}(x) + \epsilon.$6, while implicit priors encode prior information through operators or procedures such as denoisers, proximal maps, optimization constraints, or transformations, without requiring the prior density to be available in tractable form (Everink et al., 15 Sep 2025).

The paper identifies several implementation patterns. RestorationPrior encapsulates a restoration operator such as a denoiser or projector. MoreauYoshidaPrior implements the Moreau–Yoshida regularization envelope used in algorithms such as MYULA. RegularizedGaussian combines a Gaussian prior with added regularization or constraints for RLRTO-style methods. MappedGeometry supports latent-variable or transformation priors of the form

$y = \mathcal{A}(x) + \epsilon.$7

so that posterior sampling is carried out in the latent variable and mapped back to the parameter space. The summary also mentions TweediePrior and RegularizedLinearRTO as components of the implicit-prior implementation (Everink et al., 15 Sep 2025).

The supported classes of implicit priors include latent-variable or low-dimensional transform priors, Langevin-based and Plug-and-Play priors, and RLRTO priors. For Plug-and-Play ULA, the proximal operator is replaced by a denoiser $y = \mathcal{A}(x) + \epsilon.$8, leading to the update

$y = \mathcal{A}(x) + \epsilon.$9

For RLRTO, samples are defined operationally through repeated solution of regularized randomized optimization problems, rather than by evaluating a closed-form prior density (Everink et al., 15 Sep 2025).

The application examples for implicit priors include image processing tasks such as inpainting and deblurring, using PnP-ULA with denoisers such as DnCNN or wavelet thresholding, and parameter estimation in PDEs such as Poisson problems using RLRTO with TV or nonnegativity constraints, or MYULA with a TV prior implemented through a proximal operator. The source materials also note that monotonicity and stepwise structure can be incorporated through custom projectors or regularizations (Everink et al., 15 Sep 2025).

The same paper states several limitations and interpretive cautions. Many implicit priors do not admit a tractable density, which precludes traditional MCMC methods such as Metropolis-Hastings. For arbitrary denoisers there may be no well-defined target posterior and convergence guarantees may fail. Some RLRTO-style priors may be data-dependent and therefore improper in the Bayesian sense. Diagnostic workflows can also become more difficult because standard convergence and density-based validation tools rely on knowledge of the target distribution. These issues do not negate the computational role of implicit priors, but they do place methodological emphasis on algorithm-specific interpretation and validation (Everink et al., 15 Sep 2025).