Markov random fields in CUQIpy - Uncertainty Quantification in Inverse Problems with CUQIpy

In some cases, we may want to generate samples that represent a field with some spatial correlation and smoothness properties.

One CUQIpy distribution object that can be achieved for this purpose is the Gaussian Markov random field (GMRF) distribution Sidén (2020). This distribution assumes a Gaussian distribution on the differences between neighboring elements of $\mathbf{x}$ , i.e. in 1D:

\begin{align*} x_i - x_{i-1} \sim \mathrm{Gaussian}(0, d^{-1}), \quad i=1, \ldots, n, \end{align*}

(1)

where we purposely leave out the details on the boundary conditions for this notebook.

To simplify the notation, we denote by $\mathrm{GMRF}$ the distribution that induces this property on a vector $\mathbf{x}$ defined by its mean and precision $d$ . That is, the above can be written as

\begin{align*} \mathbf{x} &\sim \mathrm{GMRF}(\mathbf{0}, d), \end{align*}

(2)

with some choice of the precision say $d=50$ . For more details on GMRF see the CUQIpy paper (the first part of the 2-part article) Riis et al. (2024).

The GMRF distribution is implemented in CUQIpy as GMRF class and can be used as follows:

# Define prior precision
d = 50
n = 200
# Define GMRF prior (zero boundary conditions are the default)
x_GMRF = GMRF(np.zeros(n), d)

Exercise:¶

Can you generate and plot a realization (sample) of x_GMRF? Does the realization show spatial correlation?
Create a Gaussian distribution x_Gaussian with mean np.zeros(n) and precision 50, and compare a sample from the GMRF distribution with the Gaussian distribution by plotting them on the same plot. What do you observe?
★ Generate 100000 samples of x_GMRF and store them in variable x_GMRF_samples. Verify the following about the distribution of the differences. Focus only on verifying the difference between elements 30 and 31 as a representative example.
- The mean of the difference between elements 30 and 31 is close to 0.
- The variance of the difference between elements 30 and 31 is close to $1/50$ .
- Hint: Use this line to create a Samples object of the differences diff_30_31_samples = Samples((x_GMRF_samples.samples[31] - x_GMRF_samples.samples[30]).reshape(1, -1)).

# your code here

Other Markov random fields in CUQIpy:¶

Cauchy Markov Random Field (CMRF): CMRF class
Laplace Markov Random Field (LMRF): LMRF class

CMRF and LMRF are similar to GMRF but with different distributions on the differences between neighboring elements in the signal, where CMRF assumes a Cauchy distribution and LMRF assumes a Laplace distribution. LMRF and CMRF are particularly useful in cases in which the signal to be inferred has sharp edges (jumps).

This 1D deconvolution example from Riis et al. (2024) illustrates and compares using the three Markov random fields in a 1D problem.

References¶

Sidén, P. (2020). Scalable Bayesian spatial analysis with Gaussian Markov random fields (Vol. 790). Linköping University Electronic Press.
Riis, N. A. B., Alghamdi, A. M. A., Uribe, F., Christensen, S. L., Afkham, B. M., Hansen, P. C., & Jørgensen, J. S. (2024). CUQIpy: I. Computational uncertainty quantification for inverse problems in Python. Inverse Problems, 40(4), 045009. 10.1088/1361-6420/ad22e7