Markov random fields in CUQIpy

Markov random fields in 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 [Siden20]. 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*} \]

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*} \]

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) [RAU+24].

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

Hide code cell content
from cuqi.distribution import GMRF
import numpy as np
# Define prior precision
d = 50
n = 200
# Define GMRF prior (zero boundary conditions are the default)
x_GMRF = GMRF(np.zeros(n), d)

Exercise:#

  1. Can you generate and plot a realization (sample) of x_GMRF? Does the realization show spatial correlation?

  2. 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?

  3. ★ 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:#

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 [RAU+24] illustrates and compares using the three Markov random fields in a 1D problem.