Field dynamics inference for local and causal interactions

TL;DR

This paper presents a Bayesian method for inferring space-time fields and their correlation structures from noisy, incomplete data, based on physical assumptions like homogeneity, locality, and causality.

Contribution

It introduces a non-parametric prior model for correlation structures that leverages fundamental physical principles to improve field inference from limited data.

Findings

Successfully reconstructs fields and correlations from single noisy data realizations.

Demonstrates robustness with incomplete and noisy observational data.

Applicable across various scientific disciplines involving spatiotemporal data.

Abstract

Inference of fields defined in space and time from observational data is a core discipline in many scientific areas. This work approaches the problem in a Bayesian framework. The proposed method is based on statistically homogeneous random fields defined in space and time and demonstrates how to reconstruct the field together with its prior correlation structure from data. The prior model of the correlation structure is described in a non-parametric fashion and solely builds on fundamental physical assumptions such as space-time homogeneity, locality, and causality. These assumptions are sufficient to successfully infer the field and its prior correlation structure from noisy and incomplete data of a single realization of the process as demonstrated via multiple numerical examples.

Figures

Equations

$$S(x,y)\equiv {⟨s^x{\left(s^y\right)}^{\ast }⟩}_{P(s)},$$

$$(O_{S}s{)}^x\equiv {\int }_{I}S(x,y)s^y\mathrm{d}y,$$

$$(O_S{e}_k{)}^x={\int }_{I}S(x,y)e_k^y\mathrm{d}y={\lambda }_k{e}_k^x.$$

$$s^x=\sum _{k=-\infty }^{\infty }{e}_k^x{\tilde{s}}^k,$$

$${\tilde{s}}^k\equiv {\int }_I{\left(e_k^x\right)}^{\ast }{s}^x\mathrm{d}x\text{with}⟨{\tilde{s}}^k⟩=0,⟨{\tilde{s}}^k{\tilde{s}}^q⟩={\delta }_{kq}{\lambda }_k\forall k,q\in \text{mathds}Z.$$

$$S(x,y)=S(x-y)\forall x,y\in I.$$

$$s^{x+B}=s^x\text{and}S(x+B,y)=S(x,y),$$

$$s^x=\sum _{k=-\infty }^{\infty }{e}^{\frac{2\pi }{B}ikx}{g}^k{\xi }^k\text{with}{\left|g^k\right|}^2={\lambda }_k,{\xi }^k\sim \mathcal{G}({\xi }^k,1)\forall k\in \text{mathds}Z,$$

$$s^x=(\mathcal{F}\hat{G}\xi {)}^x\equiv \sum _{k=-\infty }^{\infty }\sum _{q=-\infty }^{\infty }{\mathcal{F}}_k^x{\hat{G}}_q^k{\xi }^q.$$

$${\tilde{s}}^k={\left({\mathcal{F}}^{†}s\right)}^k\equiv \int {\left({\mathcal{F}}_k^x\right)}^{\ast }{s}^x\mathrm{d}x.$$

$$d^i=(Rs{)}^i+n^i={\int }_I{R}_x^i{s}^x\mathrm{d}x\forall i\in \left\{1,...,M\right\}\text{with}n\sim \mathcal{G}(n,N),$$

$$d=R\mathcal{F}\hat{G}\xi +n.$$

$$\begin{array}{r}\begin{array}{rlrlrl}t]3& D_{\xi }& & ={\left({\hat{G}}^{†}{\mathcal{F}}^{†}{R}^{†}{N}^{-1}R\mathcal{F}\hat{G}+\text{mathds}1\right)}^{-1}& & =\text{mathds}1-{\hat{G}}^{†}{\mathcal{F}}^{†}{R}^{†}{\left(R\mathcal{F}\hat{G}{\hat{G}}^{†}{\mathcal{F}}^{†}{R}^{†}+N\right)}^{-1}R\mathcal{F}\hat{G}\\ & m_{\xi }& & =D_{\xi }{\hat{G}}^{†}{\mathcal{F}}^{†}{R}^{†}{N}^{-1}d& & ={\hat{G}}^{†}{\mathcal{F}}^{†}{R}^{†}{\left(R\mathcal{F}\hat{G}{\hat{G}}^{†}{\mathcal{F}}^{†}{R}^{†}+N\right)}^{-1}d\end{array}\end{array}$$

$$m=\mathcal{F}\hat{G}{m}_{\xi }\text{and}D=\mathcal{F}\hat{G}{D}_{\xi }{\hat{G}}^{†}{\mathcal{F}}^{†}.$$

$${\bar{s}}^x\equiv \sum _{k\in W}{g}^k{\xi }^k{e}^{\frac{2\pi }{B}ikx}\text{with}W\equiv \left\{-\text{nicefrac}K2,...,\text{nicefrac}K2\right\}.$$

$${\left({ϵ}^x\right)}^2\equiv ⟨{\left|s^x-{\bar{s}}^x\right|}^2⟩=\sum _{k\in \text{mathds}Z\setminus W}⟨{\left|g^k{\xi }^k{e}^{\frac{2\pi }{B}ikx}\right|}^2⟩=\sum _{k\in \text{mathds}Z\setminus W}{\left|g^k\right|}^2\forall x\in I,$$

$${\left({ϵ}_R^j\right)}^2\equiv ⟨{\left|{\left(Rs\right)}^j-{\left(R\bar{s}\right)}^j\right|}^2⟩=\sum _{k\in \text{mathds}Z\setminus W}{\left|g^k\right|}^2{\left|{\left(R{e}_k\right)}^i\right|}^2\text{with}{e}_k^x=e^{\frac{2\pi }{B}ikx}.$$

$$\Omega \equiv \left[0,T\right]\otimes \left(\underset{i=1}{\overset{d}{⨂}}\left[0,B_i\right]\right),$$

$$G(x,x^{\prime })=G(x-x^{\prime })\equiv \left(\mathcal{F}\hat{G}{\mathcal{F}}^{†}\right)(x,x^{\prime })=\sum _{\omega =-\infty }^{\infty }\sum _{\mathbf{x}=-\infty }^{\infty }{g}^{\omega \mathbf{k}}{e}^{2\pi i\left(\omega \frac{t-t^{\prime }}{2T}+\mathbf{k}\cdot \frac{\mathbf{x}-{\mathbf{x}}^{\prime }}{\mathbf{B}}\right)},$$

$$s^x={\left(\mathcal{F}\hat{G}\xi \right)}^x={\int }_{\Omega }G(x-x^{\prime }){\xi }^{x^{\prime }}\mathrm{d}{x}^{\prime }\text{with}{\xi }^{x^{\prime }}\equiv {\left(\mathcal{F}\xi \right)}^{x^{\prime }}.$$

$$G_{\mathrm{c}}(x-x^{\prime })=\Theta (t-t^{\prime })G(x-x^{\prime }),$$

$$G_{\mathrm{c}}(x-x^{\prime })=\Theta (t-t^{\prime })\Theta (T-(t-t^{\prime }))G(x-x^{\prime }).$$

$$G_{\mathrm{lc}}(x-x^{\prime })=\Theta (l^2)G_{\mathrm{c}}(x-x^{\prime }),$$

$$l^2=(t-t^{\prime }{)}^2-(\mathbf{x}-{\mathbf{x}}^{\prime }{)}^{†}{\mathbf{C}}^{-1}(\mathbf{x}-{\mathbf{x}}^{\prime }).$$

$$\mathcal{L}s=0,$$

$$\left(\mathcal{L}G\right)(x-x^{\prime })=\delta (x-x^{\prime }).$$

$$\mathcal{L}(x-x^{\prime })=\left(\mathcal{F}\hat{f}{\mathcal{F}}^{†}\right),$$

$$G_{k^{\prime }}^k=\frac{1}{f^k}{\delta }_{k^{\prime }}^k,$$

$$g^k=\frac{1}{f^k}.$$

$$f^k=f_r(k)+i{f}_i(k)\text{with}k\in \text{mathds}Z,f_r,f_i\in L^2(\text{mathds}R).$$

Full text

Philipp Frank* Reimar Leike Torsten A. Enßlin

{affiliations}

Ph. Frank, R. Leike, PD. Dr. T. A. Enßlin

Max-Planck Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85748, Garching, Germany

Ludwig-Maximilians-Universität München, Geschwister-Scholl-Platz 1, 80539, München, Germany

philipp@mpa-garching.mpg.de

1 Introduction

Modeling as well as inferring quantities defined in space and time on the basis of observational data thereof has always been at the very core of many scientific areas. In recent years, astrophysical imaging began to become sensitive to the temporal dimension, in addition to the spatial ones. This is due to the fact that although large astrophysical objects such as galaxies appear to be static on observational timescales, small objects such as stars and binary black holes exhibit transient, periodic, and quasi-periodic modulations of the emission on observable timescales.

In addition, the spatio-temporal correlation structure of non-astronomical systems plays a central role in the calibration of modern telescopes. In particular the temporal variability of these systems is used in order to identify and distinguish them from the typically static astronomical object of interest. As a prominent example, modern radio telescopes such as LOFAR [1] or the upcoming SKA are limited in resolution by the deflection of incoming radio signals due to the ionosphere. The strength of these distortions ultimately depends on the electron density of the ionosphere. As this density is not known for all observed locations $x$ at time $t$, it has to be inferred along with the incoming flux. Typically, the electron density is probed via observing a calibration target with known flux at location $x^{\prime}$ and time $t^{\prime}$. Therefore, it is also necessary to make a statement about the correlation structure of the electron density in order to extrapolate the information gained at $(t^{\prime},x^{\prime})$ to the space-time location $(t,x)$ where the actual observation is made.

In order to tackle these as well as other inference problems of this kind in space and time, we have to perform inference of continuous quantities, or fields, from a finite set of measurement data. This problem is in general ill-posed, as we aim to constrain infinite degrees of freedoms (dofs) on finite, usually also noisy, measurements. Consequently we rely on Bayesian inference, more precisely on Information Field Theory (IFT) [2], and use this language to encode prior knowledge about the system under consideration. Typically, this prior knowledge is incomplete, and there exist a set of unknown hyper-parameters, such as the spatio-temporal prior correlation structure, which have to be inferred along with the field of interest. We outline in this paper how physically motivated concepts such as spatio-temporal homogeneity, locality as well as causality can be encoded into the prior correlation structure. Furthermore, we demonstrate that the resulting hierarchical prior model is restrictive enough to perform inference on the basis of noisy and incomplete data of a single realization of the process, while still being flexible enough to capture complicated correlation structures.

Traditionally there exist two different pictures on random fields in space-time, one results in space and time being treated separately [3, 4], while the other models the field as defined over a single space, namely space-time [5, 6] (see e.g. [7] for an extensive discussion). In this work we rely on the latter picture. Consequently the corresponding inference problems can be regarded as the task of inferring a field defined on a single space, given a finite amount of measurements in this space.

To this end, in section 2 we start with a brief introduction to IFT and the notation used in this work. In section 3 we discuss how to encode our prior knowledge about the field and its correlation structure into a joint prior distribution thereof. This prior is then used in section 4 in order to solve the corresponding inference problem and the performance of the resulting algorithm is demonstrated in section 5. Ultimately, in section 6, we conclude the paper with a brief summary of the proposed concepts.

2 Information Field Theory and Gaussian processes

Information Field Theory is a statistical field theory that aims to describe Bayesian inference of fields defined over some continuous space, or space-time. For simplicity, we first consider a one dimensional random process, and provide an extension to space-time at the end of this section. To this end, consider a zero mean square-integrable random process $s$ defined over a closed interval $I=[0,B]$, i.E. $s^{x}\in L^{2}(I)$ with probability $P(s)$. We define the covariance function $S$ as

[TABLE]

where ∗ denotes complex conjugation. If we associate with $S$ a linear operator $O_{S}:L^{2}(I)\rightarrow L^{2}(I)$ and define its application via

[TABLE]

we may define the eigenvalues $\lambda_{k}$ and eigenfunctions $e_{k}$ of the linear operator $O_{S}$ via

[TABLE]

Since the eigenfunctions of $O_{S}$ form an orthonormal basis and the random process $s$ lies within the span of $e_{k}$, the Karhunen-Loève theorem [8, 9] states that $s$ may be represented in this basis as

[TABLE]

with the modes $\tilde{s}^{k}$ defined via

[TABLE]

Specifically all $\tilde{s}^{k}$ become zero mean and uncorrelated random variables with variance $\lambda_{k}$. Consequently, an inference problem of $s^{x}$ can be reduced to an inference problem in $\tilde{s}^{k}$

2.1 Statistically homogeneous Gaussian processes

Statistically homogeneous Gaussian processes are a special, but very useful, process for prior modeling of physical processes as the statistical homogeneity implies that a priori no specific location in $I$ is singled out. We may again define a zero mean random process $s^{x}$ with the additional requirement that the covariance takes the form

[TABLE]

If we additionally require the space to obey periodic boundary conditions such that

[TABLE]

the Wiener-Khinchin theorem [10] implies that the eigenbasis of the linear operator associated with this covariance function is the Fourier basis and its spectrum is the Fourier power spectrum. This allows for a representation of $s$ as

[TABLE]

where $i$ denotes the imaginary unit and $\xi^{k}$ are independent and identically distributed Gaussian random variables with mean zero and variance one. For a compact notation we may define the Fourier transformation $\mathcal{F}^{k}_{x}=e^{\frac{2\pi}{B}ikx}$ and an infinite dimensional diagonal matrix $\widehat{G}^{k}_{q}=\delta_{kq}g^{k}$ to write

[TABLE]

The application of its adjoint, abbreviated as $\mathcal{F}^{\dagger}$, is defined as

[TABLE]

2.2 Linear Measurements and the Wiener Filter

In order to perform Bayesian inference of $s$ given some observational data $d$, we require a data generating model (or data-model) that describes how $d$ is obtained from $s$ and possibly additional nuisance parameters $n$ that describe the measurement noise. A simple but very powerful idealized model is a linear measurement of $s$ with additive Gaussian distributed noise $n$, independent of $s$, defined as

[TABLE]

where $R:L^{2}(I)\rightarrow\mathds{R}^{M}$ is a linear operator that maps onto a discrete $M$-dimensional space, called data space, and $\mathcal{G}(n,N)$ denotes a Gaussian distribution with zero mean and covariance $N$. A linear measurement operator may represent common scenarios such as measurements of values at single locations, integrated measurements over a specific area, partially masked areas, convolution with a point spread function, and linear combinations thereof.

We can represent the data-model (Eq. (11)) as a generating process in terms of $\xi$ by inserting the Fourier basis representation of $s$ (Eq. (9)) to get

[TABLE]

The inference problem may thus be regarded as the task of constructing the posterior distribution of $\xi$, given $d$ and the background information i.E. the specific form of $R$ and the prior spectrum $g$ which in turn defines $\widehat{G}$. This problem allows for a closed form solution by means of quadratic completion (see e.g. [2]) and the posterior remains a Gaussian distribution with mean $m_{\xi}$ and covariance $D_{\xi}$ given as

[TABLE]

where the first part of the equations is the common representation of the Wiener Filter, now for a infinite number random variables $\xi$ which are the coefficients of the random process $s$ in the eigenbasis of its prior. The right hand side can be obtained by straightforward manipulation of the expressions. It has the convenient property that the only matrix inversion involved appears in the finite dimensional data-space and thus entirely avoids inversion of infinite dimensional matrices.

The posterior of the coefficients $\xi$ can be used to construct the posterior of $s$ by insertion of the modes into the expansion of $s$ (Eq. (9)). Therefore the posterior mean $m$ and the covariance $D$ of $s$ are denoted as

[TABLE]

This concludes the description of the Wiener Filter theory applied to square integrable random processes in terms of the eigenbasis of the linear operator associated with the prior covariance. A mathematically more rigorous and coordinate free discussion of these concepts is beyond the scope of this work, but is described in great detail by e.g. [11].

2.3 Consistent discretization

For many physically relevant choices of $R$ and $\widehat{G}$ the Fourier integrals involved in Eqs. (13) and (14) may be difficult to solve or may not have a closed form representation analytically. Therefore, for practical applications, the inference problem is often discretized and the discrete problem is solved instead. However, as shown by e.g. [12], care must be taken when defining a discretization in order to ensure that the finite dimensional approximation is consistent with the infinite dimensional inference problem. In this work, we achieve a discrete representation by truncating the Fourier series at a maximal / minimal value $\pm\nicefrac{{K}}{{2}}$. Specifically

[TABLE]

A measure for the discretization error may be defined by means of the expected squared difference between $\bar{s}$ and $s$ as

[TABLE]

which quantifies the difference between the infinite dimensional process and the finite dimensional approximation of the quantity of interest $s$. In order to ensure that inference is consistent, the discretization error $\epsilon_{R}$ of the observed quantity $Rs$ has to be considered. It is given as

[TABLE]

A small $\epsilon_{R}$ is sufficient to ensure that the discrete approximation of the inference problem is close to the continuous one as it ensures that the contribution of modes not in $W$ to the observed quantity is small and therefore the information gained about these modes via the observation is also small compared to the information gain about the modes in $W$. For posterior analysis of $s$ it is also relevant to have a good discrete approximation and therefore in general also $\epsilon$ should be small.

A minimal requirement is that the Gaussian process is continuous, which implies that $\left|g^{k}\right|^{2}$ decays asymptotically at least with $\nicefrac{{1}}{{|k|^{2+\gamma}}}\ ,\gamma\geq 0$. This ensures that the series expansion of $s$ converges and that there exists a $K$ such that $\epsilon$ becomes small. Fortunately, the assumption of continuity is met by most physically relevant processes.

In general, the magnitude of $\epsilon_{R}$ can only be specified for a given measurement scenario as it depends on the specific form of the measurement operator $R$. However, as the properties of $R$ can fully be defined via its action on the Fourier basis $e_{k}$, we may qualitatively discuss three distinct cases. First, consider the case where $\left|Re_{k}\right|^{2}\sim\mathcal{O}(1)$. Typical examples are the measurements of individual locations or sub intervals of $I$. In this case $\epsilon_{R}$ is comparable to $\epsilon$. The second case are measurements that suppress small scales, as for example integration over an interval or convolution with a spatially extended kernel such as a point spread function. In case of integration, we get that $\left|Re_{k}\right|^{2}\propto\nicefrac{{1}}{{|k|^{2}}}$ and thus $\epsilon_{R}$ becomes smaller then $\epsilon$. In these two cases the discretization error of the observable is comparable or smaller then $\epsilon$ and thus a small discretization error for the field is sufficient for a consistent reconstruction. The third case are measurement operators which amplify small scale structures. One important special example are measurement operations involving spatial derivatives. The action of the derivative on $e_{k}$ leads to a multiplicative factor proportional to $k$ and therefore $\left|Re_{k}\right|^{2}\propto|k|^{2}$. Care must be taken in this case since $\epsilon_{R}$ may become large or even infinite even though $\epsilon$ is small. This also shows that not all combinations of $g$ and $R$ lead to an inference problem that allows for a consistent finite dimensional representation.

Nevertheless, given a consistent combination of $g$ and $R$, there always exists a cutoff $K$, that can be chosen prior to the reconstruction, for which the error between the discrete representation and the continuous inference problem becomes small.

2.3.1 Unknown prior spectrum

In this work, however, we aim to infer the prior correlation structure, specifically the form of $g$, in addition to $\xi$, from the observed data. This poses a problem, since we cannot deduce a sensible choice for $K$ a priori. In some cases, the measurement setup allows to provide an estimate for the small scale / asymptotic behaviour of the true spectrum from the observed data, and therefore allows to set $K$ accordingly. In many cases, however, such an estimate is not feasible without performing the full reconstruction. One approach to resolve this issue is that after a reconstruction with a chosen $K$, a new reconstruction using a larger cutoff $K^{\prime}>K$ is performed, and both results are compared. If the two results obtained this way are similar, specifically if the difference between the large and small reconstruction for $s$ as well as the observable $Rs$ are small, we may conclude that the chosen discretization sufficiently resolves the true underlying process. If this is not the case, the reconstruction has to be re-run with an even larger cutoff and the procedure has to be repeated, until the deviations become small enough. This ensures that the reconstruction is consistent with the infinite dimensional problem, assuming that the true observed process has a further decaying spectrum for modes above the cutoff.

2.4 Higher dimensional representation in space-time

In order perform inference in a space-time setting we define a $d+1$ dimensional space $\Omega$, where $d$ is the number of spatial dimensions. Specifically we define

[TABLE]

and impose periodic boundary conditions along each axis of $\Omega$. This allows for a multi-dimensional Fourier series expansion of random processes $s\in L^{2}(\Omega)$, that are statistically homogeneous in space and time, in direct analogy to the one dimensional setting. Furthermore we may label coordinates on $\Omega$ via $x=\left(t,\mathbf{x}\right)$, and additionally $k=\left(\omega,\mathbf{k}\right)$ labels the associated Fourier modes.

The periodic boundary conditions introduce possibly unwanted correlations between the boundaries, in particular along the time axis. To avoid this, we extend the time domain to $2T$, to be twice the size of the observed domain, and ensure by construction of the prior of $g$ that the process becomes uncorrelated for moments in time with a distance greater than $T$. Similar to the procedure of finding an appropriate discretization, we may ensure that the posterior of $g$ has sufficiently decreased within the interval $[0,T]$ after the reconstruction. If this requirement is not met, the time domain has to be enlarged even further and the reconstruction is performed again. For the sake of simplicity, in all examples of this work, we keep the periodic boundaries in the spatial dimensions. However, these may be omitted as well by an analogous extension of the space.

3 Prior

In the introduction we proposed three different concepts we aim to encode into our prior model, particularly into the prior of the correlation structure, defined via $g$. These are

• •

Statistical space-time homogeneity

• •

(Spatio-) temporal causality

• •

Locality.

The first concept is already satisfied by construction as we assume that the correlation structure is diagonal in the Fourier representation and thus fully specified via the Fourier spectrum $g$. We may associate a Linear operator $G$ with $g$, by Fourier transforming the diagonal Matrix $\widehat{G}$ defined in Eq. (9) as

[TABLE]

where $\mathbf{B}=\left(B_{1},...,B_{d}\right)$ denotes the length of the spatial domain along each axis and $\cdot$ denotes the scalar product. This is the representation of $G$ in space-time coordinates. Indeed we find that

[TABLE]

This equivalence allows for a convenient physical interpretation of $G$ and $\xi$. Assume that $s$ models the deviations of a physical system from its steady state, induced via an external force, an excitation. Furthermore assume that the response of the system to such an excitation is stationary in space and time, i.E. is the same irrespective of where the excitation happened. This implies that $\xi$ plays the role of the external excitation, while $G$ models the stationary response (also called Green’s function) of the system to $\xi$. This interpretation is purely artificial at this point, however it motivates the prior concepts that we aim to include into $G$ as they are fundamental for the response of a physical system.

The second concept, causality, is introduced via an additional constraint on $G$. As $G$ should model the response of the system to the excitations $\xi$, it should not contain a response at times before an excitation happens. This can be formulated mathematically as

[TABLE]

where $\Theta$ denotes a step function in time. This trivially ensures that no response happens before an excitation. Technically, as we imposed periodic boundary conditions in time, the constraint has to be modified such that

[TABLE]

The additional step function ensures that an excitation at $t^{\prime}$ can only cause a response within the interval $[t^{\prime},T+t^{\prime}]$ but leaves later times unaffected. Together with the expansion of the space by a factor of $2$ as discussed in section 2.4 it ensures that excitations at $T$ do not wrap around the space to cause a response at $t=0$.

In space-time, the physical concept of causality also implies a maximal finite propagation speed of interactions. This leads to propagation within a light cone, as depicted in figure 1. Therefore, given a maximal propagation speed $c$, we can restrict the reconstruction to Greens functions that are non-zero only within the light cone. We can implement this constraint by extending eq. 21 as

[TABLE]

where

[TABLE]

This ensures that the propagator is only non-zero within the light cone. In general we might expect different maximal propagation speeds in different directions, and consequently $\mathbf{C}$ becomes a symmetric tensor. In case propagation is isotropic in space with speed $c$ we have $\mathbf{C}_{ij}=c^{2}\delta_{ij}$. Note that even in cases where we do not know $\mathbf{C}$, for example if the medium in which the interaction is realized is unknown, such a constraint can also be useful if we elevate $\mathbf{C}$ to be an additional unknown parameter that has to be inferred. Since $\mathbf{C}$ is assumed to be the same for all scales, an inference algorithm can use large scale information, where the signal to noise ration (SNR) is usually higher, to determine $\mathbf{C}$, which then effectively increases the SNR for smaller scales. A detailed description how $\mathbf{C}$ enters the reconstruction can be found in appendix A.

In order to encode the last concept, locality, we first have to revisit some properties of Green’s functions. Consider a system, undisturbed by external excitation, that can be described via

[TABLE]

where $\mathcal{L}$ is a linear differential operator. The response $G$ of such a system to an external force has to fulfill

[TABLE]

In order to give rise to a homogeneous $G$, $\mathcal{L}$ also has to be homogeneous and therefore is diagonal in Fourier space

[TABLE]

where we also introduced the diagonal Fourier space representation of $\mathcal{L}$ denoted via the complex Fourier coefficients $f^{k}$. Consequently $G$ takes the form

[TABLE]

and therefore the eigen-spectrum of $G$ reads

[TABLE]

It turns out that locality can be encoded more intuitively in terms of a prior for $f$. As we can see in Figure 2, the locality of a response is related to the order in the derivatives of the differential operator $\mathcal{L}$. Low order derivatives result in an almost instantaneous response of the system while higher order derivatives lead to an apparent non-local response. Therefore we seek to formulate a prior for $f$ such that lower order derivatives are a priori favoured. Nevertheless it should be possible that $f$ can deviate from this assumption if there is enough evidence in the data to support this. We impose a certain degree of smoothness for $f$, i.E. we want that two modes $f^{k}$, $f^{k^{\prime}}$ are correlated, where the correlation decays as $|k-k^{\prime}|$ increases. To this end we define the complex modes $f$ as

[TABLE]

In words, the real and imaginary part of the modes $f^{k}$ are defined via two square integrable functions $f_{\nicefrac{{r}}{{i}}}$ which get evaluated at the integer spaced Fourier locations $k$. We place a Gaussian process prior on both functions $f_{r}$ and $f_{i}$ of the form

[TABLE]

where the exponential factor reads

[TABLE]

where $\Delta_{k}$ denotes the Laplace operator, $\sigma$ is an overall scaling parameter that steers the strength of this prior, and $\gamma$ is a low frequency cutoff to ensure that the prior is proper. For further details see [13, 14]. The covariance function associated with $T$ takes the form

[TABLE]

Therefore the real and imaginary part of the Fourier modes $f^{k}$, which are just the functions $f_{\nicefrac{{r}}{{i}}}$ evaluated at $k\in\mathds{Z}$, are defined to be two independent, infinite dimensional Gaussian random vectors with zero mean and covariance

[TABLE]

All in all, we can combine the above concepts to end up with a generative description of our prior which takes the form

[TABLE]

where $\widehat{g_{\mathrm{lc}}}$ denotes a diagonal matrix in Fourier space with $(g_{\mathrm{lc}})^{k}$ on its diagonal and $a=\left(f,\xi,\mathbf{C}\right)$ denote the quantities of interest in the generative process. Furthermore

[TABLE]

with $f$ being a priori distributed according to eq. (31).

3.1 Comparison to Matèrn type and other parametric kernels

There exists a vast literature about Gaussian process priors with a stationary covariance [15, 16] which discuss a great variety of different forms of covariance functions. Two important classes are squared exponential, and the Matèrn class of covariance Functions. As a stationary covariance allows for a diagonal representation in Fourier space, it makes sense to compare the spectra of the associated operators. Squared exponential kernels imply that the spectrum takes the form

[TABLE]

While such kernels are very popular as they are particularly easy to implement and use in practice, the quadratic-exponential suppression of small scale structures often appears to be non-physical. A physically better motivated type of spectrum is provided by Matèrn covariances which give rise to spectra of the form

[TABLE]

As many physical processes can be well approximated as a power-law, this parametrization provides a more sensible statistical structure. In addition the large-scale cutoff ensures that the process is well defined as the variance remains finite for all $k$. We notice that the denominator of this spectrum is very well represented by the prior process we impose on $f$ and thus the Matèrn class of covariances is represented in our prior assumptions. In contrast to a fixed Matèrn covariance function, however, the form of $f$ remains unknown prior to the reconstruction and thus is inferred to match the observed data. In addition, a non-parametric process for $f$ appears to be more flexible in modeling deviations from a (possibly idealized) power-law shape of the spectrum.

3.2 Prior distributions for excitations

So far we restricted the discussion to independent Gaussian distributed excitations $\xi$ which, for a given fixed $G$, give rise to a Gaussian process for $s$.

From the physical perspective of $s$ being the result of a dynamic response $G$ to external excitations $\xi$, it is not necessary that $\xi$ is Gaussian distributed. In fact, in many applications it might be more realistic to define a different prior for the excitations. In order to demonstrate the implications of a non-Gaussian prior on $\xi$, in section 5, we show results for the inference problem in cases where $\xi$ is distributed according to an inverse-gamma distribution at each location in space-time. This prior is typically used as a sparsity prior, and in our case results in a system that is subject to sparse external excitations. Of course, many other prior distributions are also reasonable and important, but for the sake of simplicity we stick to $\xi$ being either Gaussian or inverse-gamma distributed in the examples. Note that even though physically motivated, exchanging the prior of $\xi$ to be an inverse-gamma distribution is non-trivial in the continuum limit. The goal of the examples is to demonstrate the applicability of this method also for non-Gaussian excitations, and therefore a rigorous mathematical treatment is beyond the scope of this work. In the case of inverse Gamma excitations, we therefore revert to the discrete representation of the process and leave the continuous treatment to future research.

4 Inference

In section 2.2 we already discussed the inference problem in the case of a linear measurement and a Gaussian prior with known covariance. Now, with the appropriate prior for the covariance at hand, we can set up the task of inferring the covariance together with the field, given observational data about the field. To this end, consider again a linear measurement as defined in eq. (11) which in terms of $a$ takes the form

[TABLE]

where $s(a)$ is given via Eq. (35). If we assume $n$ being Gaussian distributed with zero mean and covariance $N$ we get that the joint distribution of $(d,a)$ reads

[TABLE]

The posterior $P(a|d)$ is proportional to the joint distribution up to a factor that only depends on $d$ since

[TABLE]

This posterior is intractable, due to the non-linear dependency of $s$ on $a$. Consequently the corresponding inference problem cannot be solved analytically and we have to rely on a numerical approximation.

There exist a variety of different approximation techniques for Bayesian inverse problems ranging from point estimations such as the maximum a posterior (MAP) estimate, over variational approximations, to posterior sampling techniques such as Markov-Chain Monte Carlo (MCMC) [17, 18] and Hybrid Monte Carlo (HMC) [19]. As shown in [20], MAP tend to perform poorly in the task of reconstructing the excitations together with the prior correlation structure, as uncertainty information is vital to correctly estimate the prior statistics, which are missing in point estimates. While MCMC and HMC algorithms are very attractive due to their theoretical guarantees to converge to the true posterior statistics, they tend to become expensive for many astrophysical field inference applications compared to simpler, less expressive approaches. Therefore, in this work, we use a variational approximation algorithm called Metric Gaussian Variational Inference (MGVI) [20] where we approximate the true posterior with a Gaussian distribution in order to get an estimate for the mean and the covariance of the posterior. As shown in [20], MGVI provides an accurate estimate of the first moment (i.E. the posterior mean) as well as a tight lower bound on the second moment (posterior covariance) when compared against HMC techniques, while being substantially faster.

4.1 Variational Inference

In general, variational inference can be described as the task of approximating one probability distribution $P(x)$ for some quantity $x$ with another distribution $Q_{\sigma}(x)$ which, in addition, is defined up to a set of parameters $\sigma$. Approximation is then achieved via minimizing the Forward Kullbach-Leibler divergence (KL) with respect to $\sigma$. The KL is defined as

[TABLE]

where we also introduced the so called Information Hamiltonian $H$ defined as

[TABLE]

In our case, the approximate distribution $Q$ is chosen to be a Gaussian distribution in $a$ with mean $m$ and covariance $A$. For many relevant inference problems, and also for the one studied in this work, this approximation cannot be performed analytically as typically $\left<H_{P}\right>_{Q_{\sigma}}$ cannot be calculated analytically. Therefore, as discussed in section 2.3, we choose an appropriate discretization for the space-time domain $\Omega$, and perform a variational approximation of the corresponding discrete problem. It turns out that in order to achieve a reasonable resolution in space-time, the inference problems can become very high dimensional. As an example, in section 5, we show an application for a discretized 1+1 dimensional space-time with a resolution of $256\times 200$ pixels. This yields that the number of dofs in $a$ is $\approx 2\times 256\times 200\propto 10^{5}$ and consequently the number of entries in $A$ are $\propto 10^{10}$ which renders an explicit representation of $A$ on a computer to be inefficient generally. Therefore, following [20], we avoid an explicit representation by setting it to be equal to the inverse Fisher Metric $\mathcal{M}^{-1}$ [21] of the posterior, evaluated at $m$. For the posterior distribution as defined in eqs. (40) and (41), together with a Gaussian prior for $\xi$ with zero mean and unit covariance, we get that $\mathcal{M}$ takes the form

[TABLE]

Using the definition of $s$ (eq. (35)) we get that

[TABLE]

Therefore, together with the definition of $T$ (eq. (32)) we see that $\mathcal{M}$ has an implicit representation, i.e. it can be applied solely using Fourier transformations and diagonal operations, avoiding the explicit storage of this matrix at any point. Consequently the application of $A=\mathcal{M}^{-1}$ is achieved via linear solvers such as the conjugate gradient method. In addition, the structure of $\mathcal{M}$ allows for an efficient sampling of the approximate distribution $Q$. Specifically we may draw a random realization of $Q$ as

[TABLE]

where $n^{*}$ and $b^{*}$ are independent samples drawn from the noise statistics and the joint prior distribution with covariance $S$, respectively.

This is another important property since the KL involves expectation values w.r.t $Q$, which can be approximated via samples from $Q$. Ultimately, the Fisher metric is also a measure for the local curvature of the KL and therefore enables us to use second order optimization schemes to solve the corresponding optimization problem in $m$.

As discussed in the previous section, in some cases space-time causality can only be imposed if we also infer the propagation speed $\mathbf{C}$. To do so, we notice that given a prior for $\mathbf{C}$, the above still applies with the extension that $s$ is now also a function of $\mathbf{C}$. Therefore the Fisher metric gets an additional entry with the same structure as in eq. (44), for the corresponding gradient $\nicefrac{{\partial s}}{{\partial\mathbf{C}}}$ (See appendix A for further details).

All in all, the solution strategy as defined in MGVI, starting from a random initialized $m$, can be summarized as:

• •

Using $\mathcal{M}$, evaluated at the current approximate mean $m$, draw a set of samples $\{a^{*}\}$ from the approximate distribution $Q$.

• •

With these samples, calculate an estimation for the current value of the KL and its gradient.

• •

Together with the metric $\mathcal{M}$ perform a second order Newton minimization in order to get a new estimate for $m$.

• •

Repeat this procedure by re-evaluating everything using the updated $m$, until convergence.

5 Application

To demonstrate the applicability of our method we apply it to several synthetic data examples. First, we demonstrate the performance of the algorithm in a one dimensional setting, where we only aim to infer the Green’s function $G$ of the system, given the excitations and data. In the second example we perform a reconstruction of excitations that are distributed according to a unit Gaussian, as well as the Green’s function from data alone, in a 1+1 dimensional setting. In the last example we add another level of complexity via non-Gaussian statistics in the excitations. Specifically we assume the excitations to follow an inverse-gamma distribution at each location in discretized space-time. Hereby we perform source detection, the task of inferring sparse excitations at various locations in space-time, in a case where also the Green’s function is unknown. Throughout all examples, the prior model for $G$ follows the one described in section 3. Specifically, we describe $G$ in terms of $f$ and the propagation speed $\mathbf{C}$. We place a flat prior on the logarithm of $\mathbf{C}$ to ensure positivity of the propagation speed.

5.1 Implementation details

All examples were implemented using the NIFTy software package in its version 5 [22]. We included an implementation of the here introduced prior structure for the Green’s function into this publicly available package.

Throughout all examples, to solve the optimization problem associated with the MGVI algorithm, we realize the algorithm described in 4.1 with 30 total steps of the entire loop, where we draw 10 samples from the approximation $Q$ to estimate the KL during optimization. For posterior analysis, we use 50 approximate posterior samples.

For the first, one dimensional example, the optimization can be performed within less then a minute, while for the other two examples the total runtime is around 10 minutes on a standard laptop. The runtime and scaling of MGVI with model complexity as well as dimensionality is described in great detail in [20].

5.2 Temporal evolution

In our first application, shown in figure 3, we aim to demonstrate the inference of the dynamics encoding field $f$ alone in a one dimensional setting where there is only a temporal evolution to be reconstructed. The excitation field $\xi$ is known during inference. The hyper-parameters of the prior for $f$ (see Eq. (31)) are set to $(\sigma,\gamma)=(2.9,1.7)$. We generate synthetic data according to eq. 39, with $R$ being the identity. The excitations were drawn from an inverse gamma distribution to model known, sparse excitations of the system (e.g. in a laboratory setting where the unknown system is driven via sparse excitations). The synthetic signal $s$ as well as corresponding data $d$ is shown in the top-left panel of figure 3. The dynamic operator used to generate signal and data is of the form:

[TABLE]

with $(m,\gamma,\tilde{m},\tilde{\gamma})=(0.6,0.23,0.2,0.03)$. The results of the reconstruction are shown in figure 3.

We see that the reconstruction of the Green’s function is in agreement with the ground truth in the temporal domain, within uncertainties. Due to the fact that the reconstructed dynamics is uncertain, the recovered signal also has uncertainty although the excitations are known. The reconstructed Green’s function indicates that it is indeed possible to reconstruct an apparent non-local response of the system (due to higher order derivatives in this setup) since the true propagator as well as its reconstruction show oscillations that grow in the beginning of the propagator before decaying exponentially. We also notice that there is relatively high uncertainty in the first timesteps of the response $G_{t}$. This is caused by the low initial response to excitations of the true propagator. The initial part of the reconstructed propagator is purely dominated by noise and thus only constrained up to the noise level (the standard deviation of the noise $\sigma_{n}$ is set to be $\sigma_{n}=15$ and the temporal domain is discretized via 512 equidistant pixels).

We also notice that the posterior solution for the spectrum levels out for high-frequency modes (large $\omega$) below the signal to noise ratio, while the true spectrum continues to decay (ultimately also the true spectrum levels out in the numerical example due to the finite size of the considered space). The uncertainty increases in this region, but not enough to capture the true solution. Here we notice the limitations of the variational inference, which provides a local approximation of the posterior with a Gaussian. Consequently the true uncertainty might be underestimated, as in this case. However this deviation occurs orders of magnitudes below the peak of the spectrum and therefore has only a barely visible effect on the reconstruction of the signal. One way to allow for a better extrapolation to higher frequencies would be to provide a more restrictive prior for the dynamics encoding field e.g. by defining it on a polynomial basis which is a more suitable basis for this particular setup as the true dynamics is also described in terms of a polynomial. However we aim to provide a general and less restrictive approach here capable of also reconstructing non-polynomial dynamics encoding functions and consequently being less able to extrapolate to regions where we have no information provided by data.

In addition, in figure 4, we depict the reconstructed real and imaginary part of the inverse of the propagator spectrum and compare it to the spectrum associated to the differential operator $\mathcal{L}$ (Eq. (48)) used to generate the mock data. We see that the true spectrum is in agreement with the reconstruction, within uncertainty. Furthermore we notice that the posterior uncertainty is small close to the two resonant frequencies $\omega_{r}$ (see figure 4) corresponding to the two peaks in the propagator spectrum depicted in figure 3. This is due to the fact that at these frequencies both, the real and imaginary part of the spectrum, are close to zero and thus the magnitude of its inverse is large. Therefore small deviations around these values result in large changes in the corresponding realization of the random process $s$ and therefore the posterior uncertainty has to be small in order to stay consistent with the data.

To further quantify the reconstruction error, we also investigate the residual as well as the corresponding posterior uncertainty for the signal $s_{t}$ and the time representation of the Greens function $G_{t}$ (see figure 5). The residual is defined as the difference between the true solution and the posterior mean of the reconstruction. In addition to the case of $\sigma_{n}=15$, we also show the residuals and uncertainties for various other noise levels. For a better comparison, no other changes where made during reconstruction. In particular also the random number generator used to generate the synthetic data as well as the approximate posterior samples during reconstruction was seeded with the same random seed for all runs. We notice that the posterior uncertainty appears to be on a reasonable scale as the residual is within the one or two sigma confidence interval for almost all cases. In addition, we notice that the posterior uncertainty of the signal $s$ is particularly high in regions right after a strong excitation happened. This is due to the fact that the reconstruction of the Green’s function $G$, which is the response to these excitations, is most uncertain for the first timesteps. This uncertainty propagates into the uncertainty of $s$.

5.3 Space-time evolution

In our second example, shown in figure 6, we aim to reconstruct the dynamics as well as the excitations in a spatio-temporal (1+1 dimensional) setting from incomplete and noisy observations. We aim to infer the dynamical field $f$, the propagation speed $c$, and the excitations $\xi$ from noisy and incomplete measurements $d$ of the field $s$ alone. The hyper-parameters of the prior for $f$ (see Eq. (31)) are set to $(\sigma,\gamma)=(1.8,0.5)$. The dynamical system used for the generation of synthetic data is a product of a damped harmonic oscillator and an advection-diffusion generating term. This reads

[TABLE]

Furthermore, the excitations are Gaussian distributed with zero mean and unit covariance from which a single realization was drawn and convolved with the synthetic Green’s function corresponding to $\mathcal{L}$. The signal $s$, the data $d$ and the reconstruction of $s$ are shown in figure 6 for a case with $(c,m,\alpha,\beta,\tilde{m},\gamma,\tilde{\beta})=(0.4,0.16,-0.19,-0.05,0.1,0.5,0.2)$. We generate the data according to eq. 39 with a linear measurement response, which partially ($\approx 25\%$) masks the observed region, and Gaussian distributed noise with $\sigma_{n}=10$. The space-time is discretized via a regular grid with $256\times 200$ pixels, respectively.

The reconstruction algorithm is capable of reconstructing the signal in regions where we have observations thereof, while being relatively blind in unobserved regions. Consequently the posterior uncertainty is higher there. In addition, we notice that unlike the posterior mean, posterior samples consistently fill unobserved regions. Although in these regions the samples deviate strongly from the true signal, the information on the statistical properties, inferred from the observed regions, propagates into the unobserved regions due to the assumed statistical homogeneity. Therefore the posterior samples are statistically consistent throughout the entire space-time interval, which is important for posterior analysis. Furthermore we notice that there also exists variance in the statistical properties of the posterior, as can be seen for example in the difference between small scale structures of the posterior mean and the sample displayed in figure 6. This is due to the fact that also the reconstruction of the statistical properties (described via $f$) is imperfect due to the noisy data and thus subject to uncertainty. This posterior uncertainty about the small scale properties of $f$ results in a variation between different posterior samples of $f$, which ultimately propagates into the statistical properties of the corresponding sample of $s$.

In figures 8 and 8 we study the posterior properties of the Green’s function in more detail. In particular we compare the reconstructed Green’s function as well as the corresponding spectrum with the underlying ground truth. The spectrum is comparable to the ground truth in regions with sufficient SNR while it levels out in regions where we have no information given via data. In addition, the reconstructed propagator also shows oscillations consistent with the true propagator. However we notice that modes that propagate “downwards” are reconstructed well while the weaker “upwards” propagating modes are not reconstructed due to the fact that they are below the noise level. In addition, we see that deviations in posterior samples of the propagator only occur within a “cone” and remain zero outside. This is due to the fact that we also reconstruct the maximal propagation speed of the process, which is reconstructed to be $c\approx 0.45\ (\pm 0.08)$ enclosing the correct value of $c=0.4$.

5.4 Source detection

In our last example we aim to perform source detection in the excitation field, in a case where also the dynamic response is unknown. To demonstrate this scenario we again generate a 1+1-dimensional synthetic example where in this case we assume that we are only able to measure the temporal evolution of the system at several locations. In particular we measure the temporal evolution at 50 randomly selected locations of the space under consideration. This results in $\approx 80\%$ of the discrete space-time being unobserved, as the resolution is the same as in the previous example. As before, we assume that the measurements are subject to additive Gaussian noise with $\sigma_{n}=0.3$ and also assume the system to be at rest at $t=0$. The resulting data is shown in the top-left panel of figure 9. The unknown excitations are inverse gamma distributed to model strong but sparse excitations. We infer those from measurements of the system at multiple locations together with the dynamics encoding function $f$. The hyper-parameters of the prior for $f$ (see Eq. (31)) are set to $(\sigma,\gamma)=(2.9,1.7)$. The system used to generate the data exhibits damped traveling waves described by

[TABLE]

with $(c,m,\alpha,\beta)=(1.2,0.04,-0.013,0.005)$. From an information theoretical point of view this setup is very similar to the previous one since we only changed the measurement response to describe measurements of the temporal evolution at several locations, as well as the prior for excitations to be an inverse gamma prior. Consequently also the inference can be performed in the same way as before.

The setup as well as the reconstruction of the field evolving in space-time are shown in figure 9. We see that the reconstruction recovers many sources and the corresponding propagation. The algorithm uses information from the response of strong excitations to reconstruct the Green’s function of the system. Due to the assumed homogeneity in space-time, this information helps to improve the overall reconstruction in other regions. The quality of the reconstruction of single excitations additionally depends on the surrounding measurement scenario.

In figures 11 and 11 we depict the dynamic propagator as well as the spectrum and reconstructions. We can validate that the Green’s function was indeed reconstructed correctly, within uncertainties. We conclude that the task of source detection is possible even in cases where the underlying dynamics is unknown, as long as the assumptions of spatio-temporal homogeneity, causality, and locality hold.

6 Conclusion

In this work we considered the problem of reconstructing a random field $s$, defined in space and time, together with its correlation structure $S$ from noisy and incomplete data $d$ about $s$. We have shown that this Bayesian hierarchical inference problem can be reformulated to a (theoretically) equivalent problem by means of a generative process, where we aim to infer an excitation field $\xi$ as well as the dynamic response $G$. Ultimately the eigen-spectrum of $G$ was encoded in the dynamics encoding field $f$ and the propagation-speed encoding parameter $\mathbf{C}$. Together with the excitations $\xi$ they denote the quantity of interest $a$ of the inference problem. We proposed a Gaussian process prior for $f$ which gives rise to a non-parametric description of the dynamic response $G$. This gives rise to a non-linear generative process for $s$. As the proposed method is also applicable for non-Gaussian prior distributions for $\xi$ it can also model a variety of other, physically motivated, prior distributions for $s$. This flexibility is discussed for the example of an inverse gamma distribution for the excitations.

To demonstrate its applicability, the proposed method is applied to several synthetic data examples. These include a one dimensional example where the excitations were known and only $G$ had to be inferred, a 1+1 dimensional example with unknown Gaussian distributed excitations as well as an example with inverse gamma distributed $\xi$ and a measurement response that is sparse in the spatial domain.

As we restricted the prior assumptions for $G$ to the physically motivated concepts of space-time homogeneity, locality, and causality, the method appears to be applicable in a wide range of problems. One particular strength is the non-parametric formulation of the Green’s function. This becomes important in scenarios where physical models cannot provide a simple parametric description of evolution so far, to describe the Green’s function. In addition, a probabilistic description of excitations is sufficient for inference. Consequently the method is still applicable in cases where external influences cannot be described completely.

All in all, we believe that the proposed method is capable of dealing with current as well as upcoming inference problems involving fields defined over space and time, arising from the context of astrophysical imaging.

Appendix A Light cone prior on a discretized space

As discussed in section 3 the concept of causality in space-time results in a restriction of propagation within a light cone. In a continuous description, this restriction is realized via a convolution with a step function of the form

[TABLE]

Due to the fact that our calculations are ultimately performed on a finite grid, this definition appears to be somewhat problematic, as it introduces boundary effects along the edges of the step function, when realized on a discretized space. In addition, if we aim to elevate $\mathbf{C}$ to be an unknown parameter of the problem that has to be inferred, gradient based methods are no longer applicable due to the fact that the gradient is zero almost everywhere in space-time (or not defined on the boundary). Therefore we seek to find a way to relax the sharp boundary introduced via the cone, without loosing its useful properties. To do so we borrow an idea from quantum field theory where it turns out that these sharp boundaries are “smeared” out when considering them on a quantum scale. In our case the “quantum” scale can be regarded as the resolution of the discrete representation of space-time, although this analogy is purely artificial.

To achieve this relaxation, consider the following quantity

[TABLE]

It has two useful properties: For causal (including time-like as well as light-like) points the real part of $\Delta$, $\Re(\Delta)$, is zero as the square root is taken from a negative number. For non-causal (space-like) points the real part becomes positive. Furthermore, for fixed $t$, $\Re(\Delta)$ is asymptotically linear in $x$. Therefore, if we consider a Gaussian in $\Re(\Delta)$,

[TABLE]

we notice that this quantity remains one within the light cone, while it asymptotically falls off like a Gaussian in $x$, for fixed $t$. Here $\sigma$ stands for an optional scaling parameter which controls the width of the Gaussian. In all applications of this paper we replace the function $\Theta$ (Eq. (51)) with $\theta$ and set $\sigma$ to the size of a few pixels.

Acknowledgements

We would like to thank Jakob Knollmüller, Philipp Arras and Margret Westerkamp for fruitful discussions, Martin Reinecke for his contributions to NIFTy, and two anonymous referees for numerous comments that significantly improved the mathematical and overall presentation of the subject.

Conflict of interest

The authors declare no conflicts of interest.