1. Solving PDE-based BIP using core CUQIpy

Here we build a Bayesian problem in which the forward model is a partial differential equation (PDE) model, the 1D heat problem in particular.

Try to at least run through part 1 to 4 before working on the optional exercises.

Table of contents¶

• 1.1. Learning objectives

• 1.2. Loading the PDE test problem

• 1.3. Building and solving the Bayesian inverse problem

• 1.4. Parametrizing the unknown parameters via step function expansion

• 1.5. ★ Observation on a part of the domain

• 1.6. ★ Parametrizing the unknown parameters via KL expansion

1.1. Learning objectives ¶

• Solve PDE-based Bayesian problem using CUQIpy.

• Use different parametrizations of the Bayesian parameters (e.g. KL expansion, non-linear maps).

★ Indicates optional section.

1.2. Loading the PDE test problem ¶

We first import the required python standard packages that we need:

From CUQIpy we import the classes that we use in this exercise:

We load the test problem Heat1D which provides a one dimensional (1D) time-dependent heat model with zero boundary conditions. The model is discretized using finite differences.

The PDE is given by:

where $u(\xi,\tau)$ is the temperature and $c^2$ is the thermal diffusivity (assumed to be 1 here). We assume the source term $f$ is zero. The unknown parameters (random variable) for this test problem is the initial temperature profile $g(\xi):=u(\xi,0)$.

The data $y$ is a random variable containing the temperature measurements everywhere in the domain $[0, L]$ at the final time $\tau^\mathrm{max}$ corrupted by Gaussian noise:

where $\mathcal{G}(g)$ is the forward model that maps the initial condition $g$ to the final time solution via solving the 1D time-dependent heat problem. $\eta$ is the measurement noise.

Given observed data $y^\text{obs}$ the task is to infer the initial temperature profile $g$.

Before we load the Heat1D problem, let us set the parameters: final time $\tau^\mathrm{max}$, number of finite difference nodes $N$, and the length of the domain $L$

We assume that the exact initial condition that we want to infer is a step function with three pieces. Here we use commonly used Python and NumPy functions to define the initial condition.

We plot the exact solution for each node $\xi_i$:

While it is possible to create a PDE model from scratch, for this notebook we use the Heat1D test problem that is already implemented in CUQIpy.

We may extract components of a Heat1D instance by calling the get_components method.

Let us take a look at what we obtain from the test problem. We view the model:

Note that the forward model parameter, named x here, is the unknown parameter which we want to infer. It represents the initial condition $g$ above, or a a parameterization of it.

We can look at the returned data:

And the problemInfo:

Now let us plot the exact solution (exact initial condition) of this inverse problem and the exact and noisy data (the final time solution before and after adding observation noise):

Note that the values of the initial temperature profile and the data at 0 and $L$ are not included in this plot.

In case the parameters were changed in the exercise, we reset them here

1.3. Building and solving the Bayesian inverse problem ¶

The joint distribution of the data $y$ and the parameter $x$ (where $x$ represents the unknown values of the initial condition $g$ at the grid nodes in this section) is given by

Where $p(x)$ is the prior pdf, $p(y|x)$ is the data distribution pdf. We start by defining the prior distribution $p(x)$:

Assuming an i.i.d. additive Gaussian noise and given the fact that we know the exact data in this case, we can compute the noise level as follows:

And then define the data distribution $p(y|x)$:

Now that we have all the components we need, we can create the joint distribution $p(x,y)$, from which the posterior distribution can be created by setting $y=y^\text{obs}$ (we named the observed data as data in our code):

First, we define the joint distribution $p(x,y)$:

The posterior distribution pdf is given by the Bayes rule:

By setting $y=\texttt{data}$ in the joint distribution we obtain the posterior distribution:

We can now sample the posterior. Let’s try the preconditioned Crank-Nicolson (pCN) sampler (~40 seconds):

Note here in the sampler warmup phase, we set the tuning frequency to 0.01, which means that the scale of the proposal distribution is adapted every 1% of the samples.

Let’s look at the $95\%$ credible interval:

We can see that the mean reconstruction of the initial temperature profile matches the general trend of the exact solution to some extent but it does not capture the piece-wise constant nature of the exact solution.

Also we note that since the heat problem has zero boundary conditions, the reconstruction tends to go to zero at the right boundary.

1.4. Parametrizing the unknown parameters via step function expansion ¶

One way to improve the solution of this Bayesian problem is to use better prior information. Here we assume the prior is a step function with three pieces (which is exactly how we created the exact solution). This also makes the Bayesian inverse problem simpler because now we only have three unknown parameters to infer.

To test this case we pass field_type='Step' to the constructor of Heat1D, which creates a StepExpansion domain geometry for the model during initializing the Heat1D test problem. It is also possible to change the domain geometry manually after initializing the test problem, but we will not do that here.

The parameter field_params is a dictionary that is used to pass keyword arguments that the underlying domain field geometry accepts. For example StepExpansion has a keyword argument n_steps and thus we can set field_params={'n_steps': 3}.

Let’s look at the model in this case:

We then continue to create the Bayesian inverse problem (prior, data distribution and then posterior) with a prior of dimension = n_steps.

And the posterior:

We then sample the posterior using Metropolis Hastings MH sampler (~40 seconds)

Let’s take a look at the posterior:

We show the trace plot: a plot of the kernel density estimator (left) and chains (right) of the n_steps variables:

We show pair plot of 2D marginal posterior distributions:

We notice that there seems to be some burn-in, which is the initial part of the chain and appears as a trail leading to the high-density region. We show the pair plot after removing 1000 burn-in:

We can see that the burn-in, the trail part, is indeed removed. Another observation here is the correlation (or inverse correlation) between each pair of the variables.

We can also see the effect of removing the burn-in if we look at the $100\%$ credible interval before and after removing the burn-in. Let us look at the $100\%$ credible interval before:

And after removing the burn-in:

We compute the effective sample size (ESS) which approximately gives the number of independent samples in the chain, Roy (2020):

1.5. ★ Observation on a part of the domain ¶

Here we solve the same problem as in section 3 but with observing the data only on the right half of the domain.

We chose the number of steps to be 4:

Then we write the observation_nodes map which can be passed to the Heat1D. It is a lambda function that takes the forward model range grid (range_grid) as an input and generates a sub grid of the nodes where we have observations (data).

We load the Heat1D problem. Note in this case we do not pass an exactSolution. If no exactSolution is passed, the Heat1D test problem will create an exact solution.

Now let us plot the exact solution of this inverse problem and the exact and noisy data:

We then continue to create the Bayesian inverse problem (prior, data distribution and then posterior) with a prior of dimension = 4.

And the posterior:

We then sample the posterior using the Metropolis Hastings sampler (~50 seconds)

Let’s take a look at the posterior:

We see that the credible interval is wider on the side of the domain where data is not available (the left side) and narrower as we get to the right side of the domain.

1.6. ★ Parametrizing the unknown parameters via KL expansion ¶

Here we explore the Bayesian inversion for a more general initial temperature profile. We parametrize the unknown initial temperature profile using Karhunen–Loève (KL) expansion, which we introduced in section 3. More forward model examples. We denote by $\mathbf{x}$ now the expansion coefficients of the KL expansion of the initial temperature profile we are trying to infer.

Let’s load the Heat1D test case and pass field_type = 'KL', which behind the scenes will set the domain geometry of the model to be a KL expansion geometry (See KLExpansion documentation for details):

Now we inspect the model:

And the exact solution and the data:

Note that the exact solution here is a general signal that is not constructed from the basis functions. We define the prior $p(x)$ as an i.i.d. multivariate Gaussian with variance 3². We determined the choice of this variance value through trial and error approach to get a reasonable Bayesian reconstruction of the exact solution:

We define the data distribution:

And the posterior distribution:

We sample the posterior, here we use Component-wise Metropolis Hastings (CWMH) sampler, see section 1.2.2. Story 2 - Not all unknowns are created equal: Component-wise Metropolis Hastings CWMH for more details about CWMH. The sampling takes ~70 seconds:

Then we plot the $95\%$ credible interval (you can try plotting different credible intervals, e.g. $80\%$)

Using plot_ci with default options results in a credible interval with vanishing width at some locations where the upper and lower bounds of the credible interval seem to cross. To look into what actually happens here, let us first plot the upper and the lower bounds of the credible interval:

One might expect the upper bound to be always above the lower bound, but here we see that they intersect and switch order at some locations. This is because the credible interval is computed from the samples, which are KL coefficients samples, then converted to the original function values, the initial temperature profile over $[0, L]$. The resulting upper and lower bounds are not guaranteed to follow the expected order everywhere in $[0, L]$. Let us look at the upper and lower bounds of the credible interval in the parameter space (the KL coefficients) which are guaranteed to follow the expected order:

Note that, indeed, each coefficient in the lower bound is smaller than the corresponding coefficient in the upper bound. Let us use a different approach to plotting the credible interval. We can first convert the samples to function values first, i.e., the temperature profiles over $[0, L]$, then compute the credible interval from the resulting samples. To convert the samples to function values we run:

Then we plot the credible interval computed from the function values:

We can see that the width of the resulting credible interval is larger than zero everywhere over $[0, L]$. We can also plot the samples, converted to function values, to ensure that this credible interval is consistent with how the samples look like.

The samples indeed convey consistent information with the latter credible interval plot.

Let’s look at the effective sample size (ESS):

We note that the ESS varies considerably among the variables. We can view the trace plot for, let’s say, the first and the second variables:

We note that (as expected from the values of ESS) the chain quality of $x1$, which corresponds to the second coefficient in the KL expansion, is much better than that of the first variable $x0$. Low quality chain (where chain samples are highly correlated) indicates difficulty in exploring the corresponding parameter and possibly high sensitivity of the model to that particular parameter. Sampling methods that incorporate gradient information (which we do not explore in this notebook) are expected to work better in this situation.

A third way of looking at the credible intervals, is to look at the expansion coefficients $x_i$ credible intervals. We plot the credible intervals for these coefficients from both prior and posterior samples by passing the flag plot_par=True to plot_ci function:

The prior:

The posterior:

By comparing the two plots above, we see that the first few coefficients are inferred with higher certainty than the remaining coefficients. In parts, this is due to the nature of the Heat problem where high oscillatory initial temperature features (corresponding to the higher modes in the expansion) will be smoothed out (lost) faster and thus are harder to retrieve based on measurements from the final solution.