Gauss markov random fields

1 Creating a stationary Markov SDE with desired covariance

The Gauss-Markov Random Field approach.

Warning: I’m taking crib notes for myself here, so I lazily switch between signal processing filter terminology and probabilist terminology. I assume Bochner’s and Yaglom’s Theorems as comprehensible methods for analysing covariance kernels.

Let’s start with stationary kernels. We consider an SDE $f:\mathbb{R}\to \mathbb{R}$ at stationarity. We will let its driving noise be some Wiener process. We are concerned with deriving the parameters of the SDE such that it has a given stationary covariance function $k$.

If there are no zeros in the spectral density, then there are no poles in the inverse transfer function, and we can model it with an all-pole SDE. This includes all the classic Matérn functions. This is covered in J. Hartikainen and Särkkä (2010), and Lindgren, Rue, and Lindström (2011). Worked examples starting from a discrete time formulation are given in a tutorial introduction Grigorievskiy and Karhunen (2016).

More generally, (quasi-)periodic covariances have zeros and we need to find a full rational function approximation. Särkkä, Solin, and Hartikainen (2013) introduces one such method. Bolin and Lindgren (2011) explores a slightly different class

Solin and Särkkä (2014) has a fancier method employing resonators a.k.a. filter banks, to address a concern of Steven Reece et al. (2014) that atomic spectral peaks in the Fourier transform are not well approximated by rational functions.

Bolin and Lindgren (2011) consider a general class of realisable systems, given by $${\mathcal{L}}_{1}X(\mathbf{s})={\mathcal{L}}_2\mathcal{W}(\mathbf{s})$$ for some linear operators ${\mathcal{L}}_1$ and ${\mathcal{L}}_2.$

In the case that ${\mathcal{L}}_1$ and ${\mathcal{L}}_2$ commute, this may be put in hierarchical form: $$\begin{array}{rl}{\mathcal{L}}_1{X}_0(\mathbf{s})& =\mathcal{W}(\mathbf{s})\\ X(\mathbf{s})& ={\mathcal{L}}_2{X}_0(\mathbf{s}).\end{array}$$

They explain

$X(\mathbf{s})$ is simply ${\mathcal{L}}_2$ applied to the solution one would get to if ${\mathcal{L}}_2$ was the identity operator.

They call this a nested PDE, although AFAICT you could also say ARMA. They are particularly interested in equations of this form: $${({\kappa }^2-\mathrm{\Delta })}^{\alpha /2}X(\mathbf{s})=(b+{\mathbf{B}}^{\mathrm{\top }}\mathrm{\nabla })\mathcal{W}(\mathbf{s})$$

The SPDE generating this class of models is $$\left(\prod _{i=1}^{n_1}{({\kappa }^2-\mathrm{\Delta })}^{{\alpha }_i/2}\right)X(\mathbf{s})=\left(\prod _{i=1}^{n_2}(b_i+{\mathbf{B}}_i^{\mathrm{\top }}\mathrm{\nabla })\right)\mathcal{W}(\mathbf{s})$$

They show that spectral density for such an $X(\mathbf{s})$ is given by $$S(\mathbf{k})=\frac{{\varphi }^2}{(2\pi {)}^d}\frac{\prod _{j=1}^{n_2}(b_j^2+{\mathbf{k}}^{\mathrm{\top }}{\mathbf{B}}_j{\mathbf{B}}_j^{\mathrm{\top }}\mathbf{k})}{\prod _{j=1}^{n_1}{({\kappa }_j^2+\parallel \mathbf{k}{\parallel }^2)}^{{\alpha }_j}}.$$

2 Convolution representations

See stochastic convolution or pragmatically, assume Gaussianity and see Gaussian convolution processes.

3 Covariance representation

Suppose there is a linear SDE on domain ${\mathbb{R}}^d$ whose measure has the desired covariance structure, and ignore all questions of existence and convergence for now. We define terms of the driving noise $\epsilon$ and a linear differential operator $\mathcal{L}$ such that $$\mathcal{L}f(\mathbf{x})=\epsilon (\mathbf{x}).$$

Assume there is a Green’s function for the PDE, i.e. that for any $\mathbf{s}\in {\mathbb{R}}^d$ we may find a function $G_{\mathbf{s}}(\mathbf{x})$ such that $$\mathcal{L}{G}_{\mathbf{s}}(\mathbf{x})={\delta }_{\mathbf{s}}(\mathbf{x}).$$

The solutions of the SDE, ignoring a whole bunch of existence stuff, are then given by the convolution of these Green’s functions with the driving noise, i.e. $f({\mathbf{x}}_p)\stackrel{\text{sorta}}{=}\int G_{\mathbf{s}}({\mathbf{x}}_p)\epsilon (\mathbf{s})d\mathbf{s}.$ We use this to find the covariance of the solutions in terms of inner products of these fundamental solutions. $$\begin{array}{rl}k({\mathbf{x}}_p,{\mathbf{x}}_q)& =\mathbb{E}[f({\mathbf{x}}_p)f({\mathbf{x}}_q)]\\ & =\mathbb{E}[\int G_{\mathbf{s}}({\mathbf{x}}_p)\epsilon (\mathbf{s})d\mathbf{s}\int G_{\mathbf{t}}({\mathbf{x}}_q)\epsilon (\mathbf{t})d\mathbf{t}]\\ & =\mathbb{E}[\iint G_{\mathbf{s}}({\mathbf{x}}_p)G_{\mathbf{t}}({\mathbf{x}}_q)\epsilon (\mathbf{s})\epsilon (\mathbf{t})d\mathbf{t}d\mathbf{s}]\\ & =\iint G_{\mathbf{s}}({\mathbf{x}}_p)G_{\mathbf{t}}({\mathbf{x}}_q)\mathbb{E}[\epsilon (\mathbf{s})\epsilon (\mathbf{t})]d\mathbf{t}d\mathbf{s}\\ & =\iint G_{\mathbf{s}}({\mathbf{x}}_p)G_{\mathbf{t}}({\mathbf{x}}_q){\sigma }_{\epsilon }^2{\delta }_{\mathbf{s}}(\mathbf{t})d\mathbf{t}d\mathbf{s}& \text{ whiteness}\\ & ={\sigma }_{\epsilon }^2\int G_{\mathbf{s}}({\mathbf{x}}_p)G_{\mathbf{s}}({\mathbf{x}}_q)d\mathbf{s}\\ & ={\sigma }_{\epsilon }^2⟨G_{\cdot }({\mathbf{x}}_p),G_{\cdot }({\mathbf{x}}_q)⟩\end{array}$$

After that, the question is, given a Green’s function can you produce a linear operator that realises it?

For example, the arc-cosine kernel of order $1$ corresponding to the ReLU is $$\begin{array}{rl}k({\mathbf{x}}_p,{\mathbf{x}}_q)& =\frac{{\sigma }_{\epsilon }^2\Vert {\mathbf{x}}_p\Vert \Vert {\mathbf{x}}_q\Vert }{2\pi }(\sin |\theta |+(\pi -|\theta |)\cos \theta )\end{array}$$ so for Green’s functions inducing this to exist we would want $$\begin{array}{rl}\int G_{\mathbf{s}}({\mathbf{x}}_p)G_{\mathbf{s}}({\mathbf{x}}_q)d\mathbf{s}& =\frac{\Vert {\mathbf{x}}_p\Vert \Vert {\mathbf{x}}_q\Vert }{2\pi }(\sin |\theta |+(\pi -|\theta |)\cos \theta )\end{array}$$ For this to work we would need $G_{\mathbf{s}}(\mathbf{x})\propto \Vert \mathbf{x}\Vert .$

4 Input measures

Warning: this is just a dump of some notes from a paper I was writing; It does not make much sense RN. The essential idea I want to get at is considering different enveloping strategies for the SDE; Enveloping the input noise, for example

Suppose ${\mathbf{x}}_p,{\mathbf{x}}_q\in {\mathbb{R}}^d$. The kernel satisfies $$\begin{array}{r}k({\mathbf{x}}_p,{\mathbf{x}}_q)=\sum _{j=1}^d\frac{\partial k}{\partial x_{pj}}{x}_{pj}.\end{array}$$ Let $f$ denote the Gaussian process with covariance function $k$ and let ${\mathcal{F}}_{\mu }[f]$ denote the Fourier transform of $f$ with respect to the finite measure $\mu$. Let the Fourier transform of $\mu$ be denoted $\mathcal{F}[\mu ](\omega )=\int e^{-i{\omega }^{\mathrm{\top }}\mathbf{x}}\mu (\mathrm{d}\mathbf{x})$, so that ${\mathcal{F}}_{\mu }[f]=\mathcal{F}[f(x){\partial }_x\mu (x)]=\mathcal{F}[f(x)]\ast \mathcal{F}[\mu ].$

We have $$\begin{array}{rl}\mathbb{E}|{\mathcal{F}}_{\mu }[f](\omega ){|}^2& =\iint \mathbb{E}[f({\mathbf{x}}_p)f({\mathbf{x}}_q)]e^{-i{\omega }^{\mathrm{\top }}({\mathbf{x}}_p-{\mathbf{x}}_q)}\mu (\mathrm{d}{\mathbf{x}}_p)\mu (\mathrm{d}{\mathbf{x}}_q)\\ & =\iint k({\mathbf{x}}_p,{\mathbf{x}}_q)e^{-i{\omega }^{\mathrm{\top }}({\mathbf{x}}_p-{\mathbf{x}}_q)}\mu (\mathrm{d}{\mathbf{x}}_p)\mu (\mathrm{d}{\mathbf{x}}_q)\\ & =\iint \sum _{j=1}^d\frac{\partial k}{\partial x_{pj}}{x}_{pj}{e}^{-i{\omega }^{\mathrm{\top }}{\mathbf{x}}_p}\mu (\mathrm{d}{\mathbf{x}}_p)e^{i{\omega }^{\mathrm{\top }}{\mathbf{x}}_q}\mu (\mathrm{d}{\mathbf{x}}_q)\\ & =\int \sum _{j=1}^{d}i\frac{\partial }{\partial {\omega }_j}(\int \frac{\partial k}{\partial x_{pj}}{e}^{-i{\omega }^{\mathrm{\top }}{\mathbf{x}}_p}\mu (\mathrm{d}{\mathbf{x}}_p))e^{i{\omega }^{\mathrm{\top }}{\mathbf{x}}_q}\mu (\mathrm{d}{\mathbf{x}}_q)\\ & ={\mathcal{F}}_{\mu }^{{\mathbf{x}}_q}[\sum _{j=1}^{d}i\frac{\partial }{\partial {\omega }_j}(\int \frac{\partial k}{\partial x_{pj}}{e}^{-i{\omega }^{\mathrm{\top }}{\mathbf{x}}_p}\mu (\mathrm{d}{\mathbf{x}}_p))]\\ & ={\mathcal{F}}_{\mu }^{{\mathbf{x}}_q}\left[\sum _{j=1}^{d}i\frac{\partial }{\partial {\omega }_j}({\mathcal{F}}_{\mu }^{{\mathbf{x}}_p}\left[\frac{\partial k}{\partial x_{pj}}\right])\right].\end{array}$$ Then $$\begin{array}{rl}\mathbb{E}|{\mathcal{F}}_{\mu }[f](\omega ){|}^2& =\int \sum _{j=1}^{d}i\frac{\partial }{\partial {\omega }_j}(\int \frac{\partial k}{\partial x_{pj}}{e}^{-i{\omega }^{\mathrm{\top }}{\mathbf{x}}_p}\mu (\mathrm{d}{\mathbf{x}}_p))e^{i{\omega }^{\mathrm{\top }}{\mathbf{x}}_q}\mu (\mathrm{d}{\mathbf{x}}_q)\\ & =-\int \sum _{j=1}^d\frac{\partial }{\partial {\omega }_j}({\omega }_j\int k({\mathbf{x}}_p,{\mathbf{x}}_q)e^{-i{\omega }^{\mathrm{\top }}{\mathbf{x}}_p}\mu (\mathrm{d}{\mathbf{x}}_p))e^{i{\omega }^{\mathrm{\top }}{\mathbf{x}}_q}\mu (\mathrm{d}{\mathbf{x}}_q)\\ & =-\sum _{j=1}^d\int (\int k({\mathbf{x}}_p,{\mathbf{x}}_q)e^{-i{\omega }^{\mathrm{\top }}{\mathbf{x}}_p}\mu (\mathrm{d}{\mathbf{x}}_p))e^{i{\omega }^{\mathrm{\top }}{\mathbf{x}}_q}\mu (\mathrm{d}{\mathbf{x}}_q)\\ & \phantom{=}-\int ({\omega }_j\frac{\partial }{\partial {\omega }_j}\int k({\mathbf{x}}_p,{\mathbf{x}}_q)e^{-i{\omega }^{\mathrm{\top }}{\mathbf{x}}_p}\mu (\mathrm{d}{\mathbf{x}}_p))e^{i{\omega }^{\mathrm{\top }}{\mathbf{x}}_q}\mu (\mathrm{d}{\mathbf{x}}_q)\\ (d+1)\mathbb{E}|{\mathcal{F}}_{\mu }[f](\omega ){|}^2& =-\int ({\omega }_j\frac{\partial }{\partial {\omega }_j}\int k({\mathbf{x}}_p,{\mathbf{x}}_q)e^{-i{\omega }^{\mathrm{\top }}{\mathbf{x}}_p}\mu (\mathrm{d}{\mathbf{x}}_p))e^{i{\omega }^{\mathrm{\top }}{\mathbf{x}}_q}\mu (\mathrm{d}{\mathbf{x}}_q)\end{array}$$

5 Without stationarity via Green’s functions

Suppose our SDE may be specified in terms of a Gaussian white driving noise with variance ${\sigma }_w^2$ and an impulse response function/Green’s function, $g$ such that $$\begin{array}{r}f(x):=\int g(\mathbf{u},\mathbf{x})\mathrm{d}w(\mathbf{u}).\end{array}$$ We know that the kernel is an inner product kernel and therefore invariant to rotation about $\mathbf{0},$ i.e. for orthogonal $Q$, $k(Q{\mathbf{x}}_p,Q{\mathbf{x}}_q)=k({\mathbf{x}}_p,{\mathbf{x}}_q).$ It follows that $g(Q\mathbf{u},Q\mathbf{x})=g(\mathbf{u},\mathbf{x}).$ In fact, we may write each in dot-product form, i.e. $k({\mathbf{x}}_p,{\mathbf{x}}_q)=k({\mathbf{x}}_p\cdot {\mathbf{x}}_q)$ and $g(\mathbf{u},\mathbf{x})=g(\mathbf{u}\cdot \mathbf{x}).$ The kernel satisfies $$\begin{array}{rl}k({\mathbf{x}}_p,{\mathbf{x}}_q)& =\mathbb{E}[\int g(\mathbf{u},{\mathbf{x}}_p)\mathrm{d}w(\mathbf{u})\int g(\mathbf{v},{\mathbf{x}}_q)\mathrm{d}w(\mathbf{v})]\\ & =\mathbb{E}[\iint g(\mathbf{u},{\mathbf{x}}_p)g(\mathbf{v},{\mathbf{x}}_q)\mathrm{d}w(\mathbf{u})\mathrm{d}w(\mathbf{v})]\\ & =\iint g(\mathbf{u},{\mathbf{x}}_p)g(\mathbf{v},{\mathbf{x}}_q){\sigma }_w^2\delta (\mathbf{u},\mathbf{v})\mathrm{d}\mathbf{v}\mathrm{d}\mathbf{u}\\ & ={\sigma }_w^2\int g(\mathbf{u},{\mathbf{x}}_p)g(\mathbf{u},{\mathbf{x}}_q)\mathrm{d}\mathbf{u}\end{array}$$ Up to a scaling factor, the Green’s function is simply the covariance kernel under the assumption that the driving noise is white.

Recalling $k({\mathbf{x}}_p,{\mathbf{x}}_q)=\mathbb{E}[\psi ({\mathbf{W}}^{\mathrm{\top }}{\mathbf{x}}_q)\psi ({\mathbf{W}}^{\mathrm{\top }}{\mathbf{x}}_p)]=\mathbb{E}[\psi (Z_p)\psi (Z_q)]$ the Green’s function thus must satisfy $$\begin{array}{rl}{\sigma }_w^2\int g(\mathbf{u}\cdot {\mathbf{x}}_p)g(\mathbf{u}\cdot {\mathbf{x}}_q)\mathrm{d}\mathbf{u}& =\mathbb{E}[\psi ({\mathbf{W}}^{\mathrm{\top }}{\mathbf{x}}_q)\psi ({\mathbf{W}}^{\mathrm{\top }}{\mathbf{x}}_p)]\end{array}$$ Now we need to see how this works for individual kernels.

6 References