2. Probably the simplest BIP in the world (the long story) - Uncertainty Quantification in Inverse Problems with CUQIpy

Here we present the same BIP discussed in section 1. Probably the simplest BIP in the world (the short story), while providing additional details and exercises.

Table of contents¶

2.1. Learning objectives
2.2. The forward model
2.3. The prior
2.4. The noise distribution
2.5. The data distribution
2.6. The likelihood function
2.7. Maximum likelihood (ML) point estimate
2.8. The posterior distribution
2.9. Maximum a posteriori (MAP) estimate
2.10. Sampling from the posterior

2.1. Learning objectives ¶

Create a linear forward model object in CUQIpy and apply it to some parameters
Create a distribution object in CUQIpy that represents the prior and the noise, and sample, and visualize it
Create and design data distributions in CUQIpy for additive and multiplicative noise case and visualize it
Compute the MAP estimate in CUQIpy
Write the minimization problem that corresponds to finding the MAP estimate
Sample from the posterior distribution in CUQIpy and visualize the samples

2.2. The forward model ¶

Consider the following inverse problem: given observed data $b$ , determine $x_1$ , and $x_2$ :

b = x_1 + x_2 + e \;\;\mathrm{with}\;\; e \sim \mathrm{Gaussian}(0, 0.1)

(12)

We can write it as:

b = \mathbf{A}\mathbf{x} + e = \large(1,1\large)\binom{x_1}{x_2} + e

(13)

variable	description	dimension
$\mathbf{x}$	parameter to be inferred	2-dimensional
$\mathbf{A}$	forward model	1-by-2 matrix
$b$	data	1-dimensional
$e$	noise	1-dimensional

This problem is:

A Linear inverse problem since the forward model is linear.
Ill-posed (in the sense of Hadamard [1]) since the solution is not unique, i.e., for some given value of $b$ , e.g., $b=3$ , all points $(x_1, x_2)$ that satisfy $x_1 + x_2 = 3$ are solutions to the (noise-free) problem.

# Importing the required libraries
from cuqi.distribution import Gaussian
from cuqi.problem import BayesianProblem
from cuqi.model import LinearModel
import numpy as np
import matplotlib.pyplot as plt
from cuqi.utilities import plot_1D_density, plot_2D_density
np.random.seed(0)

/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/arviz/__init__.py:50: FutureWarning: 
ArviZ is undergoing a major refactor to improve flexibility and extensibility while maintaining a user-friendly interface.
Some upcoming changes may be backward incompatible.
For details and migration guidance, visit: https://python.arviz.org/en/latest/user_guide/migration_guide.html
  warn(

Let us define the forward model $\mathbf{A}$ , we first define the matrix:

A_matrix = np.array([[1.0, 1.0]])

Verify the dimension of the forward model matrix $\mathbf{A}$

A_matrix.shape

(1, 2)

Now we wrap A_matrix in a CUQIpy forward model object as follows:

A = LinearModel(A_matrix)

Let us test applying the forward model to some parameters:

some_x = np.array([1.0, 2.0])
print(A@some_x)

[3.]

# your code here

Geometries in CUQIpy¶

In CUQIpy, we use the concept of geometries to represent the interpretation of variables values (e.g. values of function on a 1D or 2D grid, discrete values, coefficients in some expansion, image pixels, etc). The Geometry object also defines the dimension of the variable and is equipped with methods of plotting the variable.

We notice that printing the forward model object A, for example, gives us

CUQI LinearModel: _DefaultGeometry1D[2] -> _DefaultGeometry1D[1].
    Forward parameters: ['x'].

The first geometry _DefaultGeometry1D[2] is the domain_geometry which represents the input space of the forward model, i.e., the space of the parameters $x_1$ and $x_2$ .
The second geometry _DefaultGeometry1D[1] is the range_geometry which represents the output space of the forward model, i.e., the space of the data $b$ .
The forward parameters ['x'] are the parameters that the forward model operates on, in this case, the parameters $x_1$ and $x_2$ .
You can access the domain and range geometries of the forward model object A by A.domain_geometry and A.range_geometry respectively.

print(A.domain_geometry)
print(A.range_geometry)

_DefaultGeometry1D[2]
_DefaultGeometry1D[1]

The _DefaultGeometry1D is a simple geometry that is used as a default geometry in CUQIpy if the user does not specify a geometry.
We will revisit this concept at a later stage depending on forward models needs.

2.3. The prior ¶

Bayesian approach: Use prior to express belief about solution¶

A common choice for simplicity is the zero mean Gaussian prior where the components are independent and identically distributed (i.i.d.):

\mathbf{x} \sim \mathrm{Gaussian}(\mathbf{0}, \delta^2 \mathbf{I})

(14)

The probability density function (PDF) of such a Gaussian prior is expressed as

\pi (\mathbf{x}) = \frac{1}{\sqrt{(2 \pi)^n \delta^{2n}}} \mathrm{exp}\left(-\frac{||\mathbf{x}||^2}{2\delta^2}\right)

(15)

where:

$n$ is the dimension of the parameter space (which is 2 in this specific case),
$\delta$ is the standard deviation of the prior distribution,
$\mathbf{I}$ is the identity matrix.

Let us define the prior distribution for the parameters $x_1$ and $x_2$ :

x = Gaussian(np.zeros(2), 2.5)

We can plot the prior PDF for the parameters $x_1$ and $x_2$ :

im = plot_2D_density(x, -5, 5, -5, 5)
plt.colorbar(im)

We can sample from the prior distribution:

x_samples = x.sample(1000)

We can visualize the samples from the prior distribution, one way to do this is to plot samples pair plot:

x_samples.plot_pair()

<Axes: xlabel='v0', ylabel='v1'>

Exercises

Can you generate one plot where you visualize the prior PDF and the pair plot together? (tip: after plotting the PDE, you can pass the pyplot current axis to plot_pair method by passing the keyword argument ax=plt.gca())
Does the samples seem to represent the prior distribution correctly?
x_1 and x_2 of the prior distribution are independent and identically distributed random variables. In this exercise: a. Define a new Gaussian distribution x_exercise where x_1 and x_2 are correlated with covariance of 0.7. Set x_1 variance to 1 and x_2 variance to 4 for. The covariance matrix is given by:
$\begin{bmatrix} 1 & 0.7 \\ 0.7 & 4 \end{bmatrix}$
(16)
Visualize the PDF of the new distribution x_exercise and the pair plot of the samples from the distribution in one plot. What do you observe?

# your code here

2.4. The noise distribution ¶

As mentioned earlier, we assume $e \sim \mathrm{Gaussian}(0, 0.1)$ . We can define the noise distribution as follows:

e = Gaussian(0, 0.1)

We print the noise distribution object:

print(e)

CUQI Gaussian.

We draw some samples from the noise distribution:

samples = e.sample(10000)

And visualize them. One way to do that is to use the trace plot in CUQIpy:

samples.plot_trace()

array([[<Axes: title={'center': 'v'}>, <Axes: title={'center': 'v'}>]],
      dtype=object)

On the left is the PDF of e estimated from e samples using the underlying Arviz kernel density estimation (KDE), and on the right is the chain plot of the samples.

Let us also plot the PDF of the noise distribution e, using the python function plot_pdf which uses the analytical expression of the PDF of the noise distribution directly:

plot_1D_density(e, -1.5, 1.5)

# your code here

2.5. The data distribution ¶

The noise in the measurement data follows $e \sim \mathrm{Gaussain}(0, 0.1)$ and due to the relation $b = \mathbf{A}\mathbf{x} + e$ , we can write

b | \mathbf{x} \sim \mathrm{Gaussian}(\mathbf{A}\mathbf{x}, \sigma^2\mathbf{I})

(17)

and in this case we specify $\sigma^2 = 0.1$ .

\pi (b | \mathbf{x}) = \frac{1}{\sqrt{(2 \pi)^m \sigma^{2m}}} \mathrm{exp}\left(-\frac{||\mathbf{A}\mathbf{x}- b||^2}{2\sigma^2}\right)

(18)

The data distribution is the conditional distribution of $b$ given $\mathbf{x}$ .
This PDF can only be evaluated for a given $\mathbf{x}$ .

We create the data distribution object as follows:

b = Gaussian(A@x, 0.1)

We print the data distribution object:

CUQI Gaussian. Conditioning variables ['x'].

Note that we can not sample from this distribution directly. If we try, we will get an error:

# Here we catch the error and print it
try:
    b.sample(10)
except Exception as e:
    print(e)

Cannot sample from conditional distribution. Missing conditioning variables: ['x']

Before sampling or evaluating the PDF of b, we need to specify the value of the parameter x, let us choose the following value:

particular_x = np.array([1.5, 1.5])

Then we condition the data distribution b on the given parameter particular_x:

b_given_particular_x = b(x=particular_x)

Now we have the distribution object b_given_particular_x that represents the data distribution given a particular value of the parameters x. We can now plot the PDF of this distribution:

plot_1D_density(b_given_particular_x, 1.5, 4.5)

We can use b_given_particular_x to simulate noisy data assuming that the true x parameters is particular_x:

b_obs = b_given_particular_x.sample()
print(b_obs)

3.2989960053953644

Exercises

Print the distribution object b_given_particular_x, how does the output differ from printing b.
What is the ratio of the noise in the data (the noise in b_obs) to the true data.
Explain how the data distribution is different from the noise distribution.
Assume we have a multiplicative noise case in which the level of the noise depends on the model output, i.e., $b = (\mathbf{A}\mathbf{x}) e$ where $e \sim \mathrm{Gaussain}(1, 0.1)$ .
- What is the corresponding data distribution?
- Define the data distribution in CUQIpy and call it b_mult. Plot its PDF (condition on particular_x). Tip: you can use lambda function to define the Gaussian covariance lambda x : 0.1*(A@x)**2.
- Plot the PDF of the data distribution b_mult given x=particular_x.
- In general, for the given x_particular, how does the noise to the exact data ratio change in the case of b_mult compared to b?
- Plot the PDF of the data distribution b_mult given a different value of x=np.array([0.4, 0.3]). What can we say about the noise ratio to the exact data in this case? and why?

# your code here

2.6. The likelihood function ¶

We obtain the likelihood function by fixing observed data $b^\mathrm{obs}$ in the data distribution and considering the function of $\mathbf{x}$ :

L (\mathbf{x} | b^\mathrm{obs}) \mathrel{\vcenter{:}}= \pi (b^\mathrm{obs} | \mathbf{x})

(19)

Her we have:

L (x_1, x_2 | b=b^\mathrm{obs}) = \frac{1}{\sqrt{2 \pi \cdot 0.1}} \mathrm{exp}\left(-\frac{(x_1+x_2- b^\mathrm{obs})^2}{2\cdot 0.1}\right)

(20)

In CUQIpy, we can define the likelihood function as follows:

likelihood = b(b=b_obs)
print(likelihood)

CUQI Gaussian Likelihood function. Parameters ['x'].

Note that the likelihood function is a density function and is not a distribution. If we try to compute pdf for example, we will get an error:

try:
    likelihood.pdf(x=particular_x)
except Exception as e:
    print(e)

'Likelihood' object has no attribute 'pdf'

while for example, we can evaluate the pdf for the distribution x:

x.pdf(particular_x)

array([0.02588303])

For the likelihood function, we can evaluate its log-density:

likelihood.logd(x=particular_x)

CUQIarray: NumPy array wrapped with geometry.
---------------------------------------------

Geometry:
 _DefaultGeometry1D[1]

Parameters:
 True

Array:
CUQIarray([-0.21463904])

We plot the likelihood function for the observed data b_obs:

x1_lim = np.array([-5, 5])
x2_lim = np.array([-5, 5])
plot_2D_density(
    likelihood,
    x1_lim[0], x1_lim[1],
    x2_lim[0], x2_lim[1])

2.7. Maximum likelihood (ML) point estimate ¶

The maximum likelihood (ML) estimate is equivalently the minimizer of the negative log of the likelihood. And in the case of Gaussian noise, it is the least-squares solution:

\mathbf{x}^* = \underset{\mathbf{x}}{\operatorname{argmin\;}} \frac{1}{2 \sigma^2} ||\mathbf{A}\mathbf{x}- b^\mathrm{obs}||_2^2

(21)

Again, we plot the likelihood function, but this time we add the line $x_2 = b^{obs}-x_1$ , shown as a red dashed line:

# Plot the likelihood
plot_2D_density(
    likelihood,
    x1_lim[0], x1_lim[1],
    x2_lim[0], x2_lim[1])

# Plot the line x2 = b_obs - x1
plt.plot(x1_lim, b_obs-x1_lim, '--r')
plt.ylim(x2_lim)

(-5.0, 5.0)

Note that all the points on the line $x_2 = b^{obs}-x_1$ have the same likelihood value and therefore there is no unique ML point. This is expected, since the problem we are solving is: $b^{obs} = x_1 + x_2$ . Combining the likelihood with the prior gives a unique maximum a posteriori (MAP) estimate as we will see next.

2.8. The posterior distribution ¶

Bayes’ rule¶

Bayes’ rule defines the so-called posterior distribution $\pi(\mathbf{x} | b)$ , which is the conditional distribution of the parameters $\mathbf{x}$ given the observed data $b$ :

\pi(\mathbf{x} | b) \propto \pi( b|\mathbf{x})\pi(\mathbf{x})

(22)

It shows that the posterior is proportional to the product of the likelihood $\pi( b|\mathbf{x})$ and the prior $\pi(\mathbf{x})$ . Note that $\pi( b|\mathbf{x})$ here denotes the likelihood and not the data distribution, despite often written that way.

In CUQIpy, we can use the class BayesianProblem to bundle the prior, the data distribution, and the data, then use it to explore the posterior distribution (e.g., find point estimates and sample the posterior):

BP = BayesianProblem(b, x)
print(BP)

BayesianProblem with target: 
 JointDistribution(
    Equation: 
	p(b,x) = p(b|x)p(x)
    Densities: 
	b ~ CUQI Gaussian. Conditioning variables ['x'].
	x ~ CUQI Gaussian.
 )

Now we pass the data:

BP.set_data(b=b_obs)
print(BP)

BayesianProblem with target: 
 Posterior(
    Equation:
	 p(x|b) ∝ L(x|b)p(x)
    Densities:
	b ~ CUQI Gaussian Likelihood function. Parameters ['x'].
 	x ~ CUQI Gaussian.
 )

Note the difference in the target of the BayesianProblem object before and after passing the data, where in the first case, the target is a JointDistribution of x and b, while in the second case, the target is a Posterior distribution of x given b.

2.9. Maximum a posteriori (MAP) estimate ¶

The MAP estimate is the maximizer of the posterior distribution:

\mathbf{x}^* = \underset{\mathbf{x}}{\operatorname{argmax\;}} \pi(\mathbf{x} | b)

(23)

The posterior maximizer is equivalent to the minimizer of the negative log of the posterior, which, in the case of a Gaussian noise and Gaussian prior, is the classic Tikhonov solution, see sections 3. Prior information and Bayesian inverse problems and 4. Gaussians priors for more details:

\mathbf{x}^* = \underset{\mathbf{x}}{\operatorname{argmin\;}} \frac{1}{2 \sigma^2} ||\mathbf{A}\mathbf{x}- b^\mathrm{obs}||_2^2 + \frac{1}{2\delta^2}||\mathbf{x} ||^2_2

(24)

# your code here

2.10. Sampling from the posterior ¶

The MAP estimate is a very useful point estimate, but it does not provide information about the uncertainty associated with the estimate. To quantify uncertainty in the solution, we can compute posterior variance or other statistics. In this example, we have a closed form expression of the posterior which we can readily compute these statistics from. However, in general, direct computation of such statistics might not be possible or feasible. A more general approach is to use sampling methods to explore the posterior distribution which we provide a guided exercise for, next.

Let us first plot the posterior distribution.

Note that in this example the posterior distribution is a multivariate distribution of two parameters only and it is easy to evaluate the PDF of the posterior distribution over a grid of points in the parameter space. However, typically, the posterior distribution is high-dimensional and evaluating the PDF over an n-dimensional grid is not feasible.

im = plot_2D_density(BP.posterior, x1_lim[0], x1_lim[1], x2_lim[0], x2_lim[1])
plt.colorbar(im)

Exercises

Use the result MAP estimate from the previous exercise and plot it over the posterior PDF. Are they consistent?
Use the BayesianProblem object to sample from the posterior distribution. Objects of the class BayesianProblem have a method sample_posterior that samples from the posterior distribution. You can use help(BP.sample_posterior) to see the documentation of the method.
The method sample_posterior returns a Samples object, which is the same type of object that is returned by the sample method of the Distribution class. Visualize the posterior samples using the plot_trace method of the Samples object.
Similar to what you did for the prior, plot the pair plot of the posterior samples over the prior PDF (in the same plot). Are the samples consistent with the posterior distribution PDF?

# your code here

2.11. References ¶

Latz, J. (2020). On the well-posedness of Bayesian inverse problems. SIAM/ASA Journal on Uncertainty Quantification, 8(1), 451-482.