🔢 mathematics

Discretization-free Bayesian inverse problems in distribution spaces

This paper bridges infinite-dimensional Bayesian theory and computational practice for linear inverse problems in distribution spaces by demonstrating that discretizing the unknown is unnecessary, requiring only numerical quadratures independent of any discrete representation, as validated through an analysis of X-ray tomography.

Original authors: Daniela Calvetti, Erkki Somersalo

Published 2026-07-21

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Daniela Calvetti, Erkki Somersalo

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a detective trying to solve a mystery, but the clues you have are blurry, incomplete, and mixed with static noise. This is the daily reality of "inverse problems" in science and engineering. Whether it's figuring out what's inside a human body from X-rays, reconstructing a sound source from a few microphones, or guessing the shape of an object from its shadow, scientists are constantly trying to work backward from an effect to its cause. The tricky part is that many different causes could produce the same blurry effect, making the problem "ill-posed"—like trying to guess the exact ingredients of a cake just by tasting a single crumb. To solve this, scientists use a method called Bayesian inference. Think of it as a detective who doesn't just look at the clues but also brings a "prior belief" or a hunch about how the world usually works. By combining the noisy clues with these hunches, the detective can narrow down the list of suspects to the most likely culprit.

For decades, the standard way to do this detective work on a computer has been to chop the mystery into tiny, manageable pieces, like cutting a map into a grid of pixels. You turn the smooth, continuous world into a giant spreadsheet of numbers, solve the math for that spreadsheet, and hope the answer looks right. But this "pixelation" introduces errors and forces you to decide exactly how small those pieces should be before you even start. What if the answer you get depends more on how you cut the map than on the actual clues? This paper by Calvetti and Somersalo asks a bold question: Do we really need to chop the world into pixels at all? They propose a way to solve these mysteries using the language of "distributions"—a mathematical concept that treats things like smooth waves or sharp spikes as single, continuous entities—without ever forcing them into a rigid grid.

The authors show that for a specific class of problems involving linear relationships and Gaussian (bell-curve) uncertainties, you can skip the pixelation step entirely. Instead of forcing the unknown object into a grid of numbers, they treat the unknown as a continuous, wiggly shape and the measurements as "probes" that poke at this shape to ask questions. Imagine the unknown object is a giant, invisible cloud of fog. In the old way, you would try to describe the fog by counting how many drops are in every square inch of a grid. In this new way, you simply ask the fog, "How much fog is there in this specific shape?" by shining a light through a specific stencil. The math shows that you can calculate the most likely shape of the fog and how uncertain you are about it, using only these "stencil" questions and some clever integration (summing up areas), without ever defining a grid.

The paper demonstrates this idea using X-ray tomography, the technology used in CT scans. In a traditional scan, the computer divides the body into millions of tiny pixels and tries to guess the density of each one. The authors show that you can instead define the problem using smooth mathematical functions that represent the X-ray beams and the detectors. They prove that you can calculate the "posterior" (the updated, most likely picture of the object) directly in this smooth, continuous world. The best part? Once you've done the heavy math to solve the mystery in this continuous world, you can ask any question about the result without starting over. If you decide later that you only care about a specific tiny spot in the image, or if you want to look at the image with a different level of detail, you don't need to re-run the whole simulation. You just "probe" the continuous solution with a new set of questions.

The authors are careful to note that this isn't a magic wand that solves every problem instantly; it is a theoretical framework that works beautifully for linear problems with Gaussian noise. They don't claim to have solved every inverse problem in the universe, but they have shown that for these specific cases, the "discretize-then-analyze" approach (cutting things up first) is not a necessary evil. Their method, which they call "discretization-free," suggests that we can avoid the errors introduced by forcing smooth things into jagged grids. They provide a concrete example with X-ray tomography, showing that the complex numbers needed for the solution can be calculated using standard numerical integration (like adding up areas under a curve) rather than building a massive matrix of pixels. This means the solution is flexible: you can decide how much detail you want to see after you've solved the problem, rather than being locked into a specific resolution before you begin.

In essence, this paper offers a new way to think about solving mysteries in science. It argues that we don't need to force the continuous, flowing nature of reality into a rigid digital box to understand it. By using the tools of distribution spaces, we can keep the mystery smooth and continuous until the very last moment, only "pixelating" the final answer if we actually need to print it out or display it on a screen. It's a shift from "let's cut the world into pieces and solve the pieces" to "let's ask the world questions and listen to the answers," keeping the solution fluid and adaptable until the very end.

Technical Summary: Discretization-free Bayesian Inverse Problems in Distribution Spaces

Problem Statement
Inverse problems typically involve inferring an unknown, unobservable quantity from noisy measurements of a related observable. These problems are often ill-posed, requiring regularization through prior information. While the Bayesian approach provides a rigorous framework for this by combining prior distributions with likelihood models, a significant gap exists between theoretical formulations in infinite-dimensional function spaces and practical computational implementations.

Traditionally, practical solutions rely on a "discretize-then-analyze" approach, where the unknown is approximated by a finite-dimensional vector (e.g., pixel values) before solving. This introduces truncation errors and basis-dependent artifacts. Conversely, infinite-dimensional theories (in Hilbert or Banach spaces) offer a consistent framework but often lack direct computational utility, as posterior formulas involving infinite-dimensional operators (like covariance inverses) are difficult to project onto finite dimensions without re-discretizing the problem. The central challenge addressed by this paper is how to bridge infinite-dimensional theory and practical computation without resorting to a discrete approximation of the unknown of interest.

Methodology
The authors propose a framework for Bayesian linear inverse problems formulated entirely within the space of distributions (generalized functions), avoiding the discretization of the unknown variable $X$ .

Distribution Space Framework:
- The unknown $X$ is modeled as a generalized random variable taking values in a space of distributions $H = U^*$ , where $U$ is a space of test functions (e.g., $C^\infty$ functions).
- The topology and $\sigma$ -algebra of $H$ are defined by cylinder sets generated by finite evaluations against test functions.
- Measurements are defined as dual evaluations. A finite-dimensional observation $Z$ is obtained by pairing the random variable $B$ (the noisy forward model output) with a set of test functions $\Psi = \{\psi_1, \dots, \psi_m\}$ .
Gaussian Linear Model:
- The forward model is $B = AX + E$, where $A$ is a continuous linear operator and $E$ is Gaussian noise.
- $X$ and $E$ are assumed to be independent Gaussian generalized random variables with mean and covariance operators defined via their action on test functions.
- The posterior distribution is derived using the properties of conditional Gaussian distributions. Crucially, the posterior mean and covariance operators are derived in the infinite-dimensional space before any specific measurement set for the unknown is chosen.
Posterior Computation via Dual Evaluations:
- The solution to the inverse problem is the posterior distribution of a specific finite-dimensional projection $\langle \Phi, X \rangle$ , where $\Phi = \{\phi_1, \dots, \phi_n\}$ is an arbitrary set of test functions chosen a posteriori.
- The posterior mean vector $z$ and covariance matrix $D$ are computed using finite-dimensional correlation matrices ( $C_{11}, C_{12}, C_{22}$ ) constructed from the covariance operators of $X$ and $E$ and the forward operator $A_0$ (the adjoint of $A$ ).
- These matrix entries are integrals involving the test functions and the covariance kernels. The paper emphasizes that these integrals can be approximated via numerical quadrature rules that are independent of any basis representation of $X$ .
Case Study: X-ray Tomography:
- The framework is applied to X-ray tomography (both parallel and fan beam).
- The forward operator $A_0$ maps test functions on the detector screen to "ridge functions" (constant along rays) in the imaging domain.
- The authors demonstrate that the necessary matrix entries (e.g., $C_{22}$ and $C_{12}$ ) can be computed as integrals of these ridge functions against the prior covariance kernel.
- In specific cases (e.g., white noise priors and characteristic function approximations), these integrals reduce to geometric areas of intersections, recovering standard tomography matrix entries without explicitly defining a pixel grid for the unknown.

Key Contributions and Results

Discretization-Free Formulation: The paper establishes that for Gaussian linear inverse problems, the discretization of the unknown is not a necessary step for numerical treatment. The only approximations required are numerical quadratures for integrals, which do not rely on a preselected basis for the unknown.
Posterior Flexibility: A primary result (Theorem 2.5) shows that the posterior mean and covariance operators are independent of the specific test functions used to probe the unknown. Once the inverse problem is solved in the distribution space, one can choose any finite set of test functions $\Phi$ to extract details (e.g., specific regions of interest) from the posterior density without recalculating the posterior mean or covariance or re-running the analysis.
Consistency with Discretized Models: The paper analyzes the relationship between the proposed approach and traditional discretized models. It shows that standard pixel-based approximations can be viewed as specific numerical quadrature rules applied to the discretization-free integrals. The error in such approximations is shown to converge to zero as the number of pixels increases.
Resolution of Matrix Incompatibility: In traditional discretized tomography, the forward and adjoint matrices are sometimes incompatible (one is not the transpose of the other). In the proposed framework, the relevant matrix ( $C_{22}$ ) is automatically symmetric and positive definite, and the "tomography matrix" appears only in the cross-correlation term, rendering the incompatibility issue moot.

Significance and Claims
The authors claim that this work bridges the gap between the theoretical elegance of infinite-dimensional Bayesian analysis and the practical necessity of finite-dimensional computation. By framing the problem in distribution spaces, the approach:

Avoids "basis selection" errors and truncation artifacts inherent in standard discretization.
Allows for "marginalization-free" querying of the unknown, enabling the extraction of specific spatial details without re-computing the entire solution or integrating out irrelevant variables.
Provides a consistent definition of tomography matrices that arises naturally from the geometry of the problem rather than arbitrary pixelization.

The paper maintains a conceptual focus, noting that while the framework yields computationally applicable formulas (via quadrature), detailed numerical examples and comparisons of computational efficiency are left for future investigation. The authors also note that the current formulation assumes smooth operators mapping smooth functions to smooth functions, with generalizations to non-smooth operators deferred to the appendix and future work.

More like this