3. More forward model examples - Uncertainty Quantification in Inverse Problems with CUQIpy

In 1. Forward models in CUQIpy and data generation, we introduced forward models in CUQIpy, here we discuss in some detail two more examples of forward models, a model for the 1D heat conduction problem described by a Poisson problem, and a model for 1D convolution. The latter we already introduced in 1. Forward models in CUQIpy and data generation but we elaborate on it here with more details and exercises. These two models are also used in other sections of the book to demonstrate various concepts in Bayesian inversion.

Table of contents¶

3.1. Learning objectives
3.2. Forward model: 1D Poisson
3.3. Forward model: 1D convolution

3.1. Learning objectives ¶

Create (or load pre-existing) linear and non-linear forward models in CUQIpy and use them
Create and plot input for the forward models and compute and plot the corresponding output
Mathematically define two forward models: 1D Poisson and 1D Convolution and use them in CUQIpy

3.2. Forward model: 1D Poisson ¶

Here we consider a heat conduction problem for a conductive rod of length $L = \pi$ with a varying conductivity (the conductivity of the rod changes from point to point). We assume that the rod length is large relative to its thickness and thus can be modeled as a 1D domain, we refer the reader to Kakaç et al. (2018) for further reading on heat conduction. We fix the temperature at the end-points of the rod and apply a heat source distributed along the length of the rod. We wait until the rod reaches an equilibrium temperature distribution. The equilibrium temperature of the rod is modelled using the Poisson equation as

\left\{ \begin{aligned} & \dfrac{\mathrm{d}}{\mathrm{d} \xi}\left(u(\xi) \dfrac{\mathrm{d} y(\xi)}{\mathrm{d} \xi}\right) = -f(\xi), \quad & \xi\in (0,L) \\ & y(0) = y(L) = 0. \end{aligned} \right.

(3)

Here, $\xi$ represents the spatial coordinate and $y$ represents the temperature distribution along the rod, $u(\xi)$ is the unknown conductivity of the rod and $f(\xi)$ is a deterministic heat source given by

\begin{aligned} f(\xi) = 10\exp( -\frac{ (\xi - L/2)^2} {0.02} ). \end{aligned}

(4)

To ensure that the conductivity of the rod is non-negative, we parameterize $u$ by the log conductivity $x$ as follows:

\begin{aligned} u( \cdot ) = \exp( x( \cdot ) ) \end{aligned}

(5)

where $x$ is not necessarily positive.

From this point on, we assume $x$ and $y$ denote the discretized versions of the log conductivity and temperature distribution, respectively, stemming from discretizing the Poisson equation using the finite difference method. Let us load the forward model that maps $x$ to the temperature distribution $y$ in CUQIpy. We will use the following parameters:

dim : number of equi-spaced discretization points (nodes) on the interval $[0,L]$
L : length of the rod
f : a function that represents the heat source

dim = 128
L = np.pi
f = lambda xi: 10*np.exp(-(xi-L/2)**2 / 0.02)

Then we can load the 1D Poisson forward model as follows:


A = Poisson1D(dim=dim, endpoint=L, source=f).model

A is a CUQIpy model that maps the log conductivity x to the measurements y via solving the Poisson equation above. We print the forward model to see its details.

CUQI PDEModel: Continuous1D[128] -> Continuous1D[127].
    Forward parameters: ['x'].
    PDE: SteadyStateLinearPDE.

We can look at the domain and range geometries of the forward model.

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

Continuous1D[128]
Continuous1D[127]

These geometries are of type Continuous1D which represents a 1D continuous signal/field defined on a grid. We can view the 1D grid which is stored as a numpy 1D array:

print(A.domain_geometry.grid)

[0.         0.02473695 0.0494739  0.07421085 0.0989478  0.12368475
 0.1484217  0.17315865 0.1978956  0.22263255 0.2473695  0.27210645
 0.2968434  0.32158035 0.3463173  0.37105425 0.3957912  0.42052815
 0.4452651  0.47000205 0.494739   0.51947595 0.5442129  0.56894985
 0.5936868  0.61842375 0.6431607  0.66789765 0.6926346  0.71737155
 0.7421085  0.76684545 0.7915824  0.81631935 0.8410563  0.86579325
 0.8905302  0.91526715 0.9400041  0.96474105 0.989478   1.01421495
 1.0389519  1.06368885 1.0884258  1.11316275 1.1378997  1.16263665
 1.1873736  1.21211055 1.2368475  1.26158445 1.2863214  1.31105835
 1.3357953  1.36053225 1.3852692  1.41000615 1.4347431  1.45948005
 1.484217   1.50895395 1.5336909  1.55842785 1.5831648  1.60790175
 1.6326387  1.65737565 1.6821126  1.70684955 1.7315865  1.75632345
 1.7810604  1.80579735 1.8305343  1.85527125 1.8800082  1.90474515
 1.9294821  1.95421905 1.978956   2.00369295 2.0284299  2.05316685
 2.0779038  2.10264075 2.1273777  2.15211465 2.1768516  2.20158855
 2.2263255  2.25106245 2.2757994  2.30053635 2.3252733  2.35001025
 2.3747472  2.39948415 2.4242211  2.44895805 2.473695   2.49843195
 2.5231689  2.54790585 2.5726428  2.59737975 2.6221167  2.64685365
 2.6715906  2.69632755 2.7210645  2.74580145 2.7705384  2.79527535
 2.8200123  2.84474925 2.8694862  2.89422315 2.9189601  2.94369705
 2.968434   2.99317095 3.0179079  3.04264485 3.0673818  3.09211875
 3.1168557  3.14159265]

Additionally, the properties domain_dim and range_dim of the forward model represent the dimension of the input and output of the forward model, respectively.

print(A.domain_dim)
print(A.range_dim)

128
127

Let us create an array representing a discretized constant conductivity provided at the grid nodes

some_x_array = 20*np.ones(A.domain_dim)

We can wrap the array in a CUQIarray object which is the main data structure in CUQIpy for variables (e.g. arrays and fields)

some_x = CUQIarray(some_x_array, geometry=A.domain_geometry)

Note that we pass geometry=A.domain_geometry to equip the CUQIarray object with the same geometry as the domain geometry of the forward model. This geometry will interpret the values of some_x_array as function values evaluated at grid points.

We can plot the conductivity using the plot method of the CUQIarray object

some_x.plot(marker='|')

We added the marker '|' for illustration purposes, to show the grid points of the interval [0, L] where the conductivity is evaluated. We can also evaluate the forward model at the conductivity and plot the solution, which is the temperature distribution along the rod as prescribed by the Poisson equation above.

some_y = A(some_x)
some_y.plot()

Exercises

Trying a different conductivity profile:
- Create another discretized constant conductivity another_x of value 10 as a CUQIarray object.
- Evaluate the forward model at the new conductivity profile and store the result in a variable another_y. Plot the solution.
- In one plot, compare the solutions some_y and another_y using the plot method of the CUQIarray object. What do you observe? and why?
Experimenting with setting up the map in Poisson1D:
- Execute help(Poisson1D)
- Note the map parameter of Poisson1D. We can use this parameter to transform the conductivity field before solving the PDE. Create the forward model again, name it A_map, with setting up map to be lambda x: np.exp(x) to ensure that the conductivity is always positive.
- Inspect the domain and range geometries of the new forward model A_map what is different this time, and why?
- Create a CUQIarray object x_with_map with a constant conductivity profile of value 0 and evaluate the forward model at x_with_map and store it in variable y_with_map. Make sure to use the right geometry for x_with_map.
- Plot x_with_map using x_with_map.plot(), What do you observe?
- Plot the solution y_with_map and some_y in the same plot. What do you observe? and why?
Experimenting with setting up the field_type in Poisson1D:
- Note the field_type parameter of Poisson1D. We can use this parameter to change the parameterization of the conductivity field. Create the forward model again, name it A_step, with setting both the map as in the previous exercise and the field_type to be "Step" to generate a step parameterization of the conductivity field.
- What is the domain and range geometries of the new forward model A_step?
- Create x_step to be x_step = CUQIarray(np.arange(A_step.domain_dim)*0.1, geometry=A_step.domain_geometry).
- Plot x_step using x_step.plot().
- What does the dimension A_step.domain_dim represent? and is it equal to the size of the domain geometry grid?
- Evaluate the forward model at x_step and store it in variable y_step. Plot the solution.

# your code here

3.3. Forward model: 1D convolution ¶

Convolution is a forward model that explains what happens to an input function $x$ when it goes through a system where it is convolved with another function $k$ that characterizes the system. Convolution arises, e.g., in models for certain measurement systems. Deconvolution is the “opposite”—it arises when we want to determine the input $x$ from the output of the system. A convolution of two arbitrary functions $x$ and $k$ always takes the form

y(\xi) = \int_D k(\xi - \xi') x(\xi') \, \mathrm{d} \xi'

(6)

where $y$ denotes the output (the convolved signal), $x$ is the input signal, and the function $k$ describes the system—note that $k$ is a function of a single variable. In our example, $k$ is a Gaussian function, but it can indeed be any function. In practice, a finite-dimensional representation of the convolution is employed. After discretizing the signal domain $D$ into $N$ points, the convolution model is expressed as a system of linear algebraic equations $\bm{y}=\mathbf{K}\bm{x}$ . See Bracewell (1999) for theory about convolution and its relation to the Fourier transform, and see Hansen (2002) for computational aspects of convolution and deconvolution.

Let us load the forward model that maps an input $\bm{x}$ to the convolved signal $\bm{y}$ in CUQIpy. We will use the following parameters to specify the forward model:

dim : is the number of discretization points for the signal, $N$
PSF : a function that represents the point spread function, i.e., the convolution kernel
PSF_size : the size of the PSF (less or equal to dim)
PSF_param : A parameter of the PSF, the larger the size the more the blur applied to the signal
BC : boundary conditions for the convolution

We set the parameters as follows:

dim = 201
PSF = 'gauss' # Gaussian PSF
PSF_size = np.round(dim/3)
PSF_param = np.round(dim/20)
BC='reflect' # Boundary condition

Then we can load the 1D convolution forward model as follows (note that we refer to the forward model as a convolution model which we obtain from the Deconvolution1D test problem, the latter represents the BIP of deconvolving the signal):

A_deconv = Deconvolution1D(dim=dim, PSF=PSF, PSF_size=PSF_size, PSF_param=PSF_param, BC=BC).model

Let us create a function representing a signal $\bm{x}$ that we want to convolve.

signal_function = lambda xi: (xi>np.round(dim*4/10))*(xi<np.round(dim*6/10))

We evaluate the function signal_function at the discretization grid points and wrap it in a CUQIarray object.

signal_array = signal_function(A_deconv.domain_geometry.grid)
signal = CUQIarray(signal_array, geometry=A_deconv.domain_geometry)

Let us plot the signal $\bm{x}$

signal.plot()

Now, we evaluate the forward model at the signal and plot the output, the convolved signal $\bm{y}$ :

A_deconv(signal).plot()

We see that the signal is blurred due to applying the convolution operation (the forward model).

Exercises

Type help(Deconvolution1D) to see the documentation of this test problem.
Experimenting with the convolution strength:
- Create another convolution model B with the same parameters as A_deconv but with a different PSF_param value, use PSF_param=PSF_param/2. Evaluate the forward model B at the signal and plot the solution along with A_deconv(signal) in the same plot. What do you observe?
Experimenting with setting up the PSF:
- Create a new convolution model C with the same parameters as A_deconv but with a different PSF, use PSF=np.ones(m)/m. Evaluate the forward model C at the signal and plot the solution along with A_deconv(signal) in the same plot (try for different values of m: 5, 10, 20). What do you observe?
Experimenting with setting up the BC:
- Create a new convolution model D with the same parameters as A_deconv but with a different BC, use BC="periodic". Evaluate the forward model D at the signal and plot the solution along with A_deconv(signal) in the same plot. What do you observe?
- Now create another signal another_signal as another_signal = signal + 0.1*np.exp((D.domain_geometry.grid - np.round(dim*7/8))/10). Plot the new signal. Evaluate the forward models A_deconv, and D at the new signal and plot the solutions in the same plot. What do you observe?

# your code here

References¶

Kakaç, S., Yener, Y., & Naveira-Cotta, C. P. (2018). Heat conduction. CRC press.
Bracewell, R. N. (1999). The Fourier Transform and Its Applications (3rd ed.). McGraw-Hill. https://books.google.com/books?id=ecH2KgAACAAJ
Hansen, P. C. (2002). Deconvolution and regularization with Toeplitz matrices. Numerical Algorithms, 29(4), 323–378.