Reconstruction of Gaussian and log-normal fields with spectral smoothness

We develop a method to infer log-normal random fields from measurement data affected by Gaussian noise. The log-normal model is well suited to describe strictly positive signals with fluctuations whose amplitude varies over several orders of magnitude. We use the formalism of the minimum Gibbs free energy to derive an algorithm that uses the signal's correlation structure to regularize the reconstruction. The correlation structure, described by the signal's power spectrum, is thereby reconstructed from the same data set. We show that the minimization of the Gibbs free energy, corresponding to a Gaussian approximation to the posterior marginalized over the power spectrum, is equivalent to the empirical Bayes ansatz, in which the power spectrum is fixed to its maximum a posteriori value. We further introduce a prior for the power spectrum that enforces spectral smoothness. The appropriateness of this prior in different scenarios is discussed and its effects on the reconstruction's results are demonstrated. We validate the performance of our reconstruction algorithm in a series of one- and two-dimensional test cases with varying degrees of nonlinearity and different noise levels.

Figure 1

Power spectra of the one-dimensional test case. The solid (black) line shows the theoretical power spectrum; the dashed (blue) line, the power in the randomly drawn signal realization studied here; and the dotted (green) line, the power spectrum reconstructed without smoothing. The dotted horizontal line indicates the noise power, given by ${\sigma }_n^2/100$.

Figure 2

Signal reconstruction in the one-dimensional test case. The solid (blue) line shows the randomly drawn signal realization, the crosses are noisy data points used to do the reconstruction, the dashed (green) line shows the signal reconstruction obtained without spectral smoothing, and the dotted (red) line is the Wiener filter reconstruction. The lower panel shows the same with the true signal substracted.

Figure 3

Power spectra of the one-dimensional test case with a smoothness prior. The solid (black) line shows the theoretical power spectrum; the dashed (blue) line, the power spectrum reconstructed with a smoothness prior with ${\sigma }_p^2=1000$; and the dotted (green) line, the reconstruction with a stiffer smoothness prior with ${\sigma }_p^2=10$. Hatched regions around the dashed and dotted lines are the corresponding $1\sigma$ uncertainty regions estimated from the inverse Hessian of the Hamiltonian. The dotted horizontal line indicates the noise level.

Figure 4

Inverse Hessian of the Hamiltonian with a smoothness prior with ${\sigma }_p^2=1000$. Plotted is the full matrix given by the inverse of Eq. (27), evaluated at $p=p^{\left(\mathrm{MAP}\right)}$. The diagonal of this matrix can be translated into the uncertainty interval plotted in Fig. 3, however, that does not take into account the correlations across different $k$ modes that are visible in this plot.

Figure 5

Residual power between true signals and their reconstructions averaged over 100 different realizations. The solid (black) line shows the theoretical power spectrum of the signals, given by Eq. (14) with $P_0=0.2$, $k_0=5$, and $\gamma =4$. The horizontal dotted line indicates the noise level, corresponding to ${\sigma }_n^2=0.1$. The remaining lines show the residual power for different reconstruction schemes. From top to bottom, these are reconstructions without any spectral smoothing [dotted (green) line]; with an ad hoc convolution of the power spectrum [dash-dotted (magenta) line], with the smoothness prior given by Eq. (24) with ${\sigma }_p^2=100$ [short-dashed (red) line]; and using the correct power spectrum [long-dashed (blue) line].

Figure 6

Power spectra of signals with correlation functions given by Eq. (A4) (top) and Eq. (A6) (bottom). Solid (black) lines show the theoretical power spectra; dashed (blue) lines, the power in the random field realization drawn from the theoretical power spectra; and dotted (green) lines, the reconstructed power spectra using the spectral smoothness prior given in Eq. (24) with ${\sigma }_p^2=100$. Parameters substituted into Eq. (A4) are $\sigma =0.2$ and $C_0=1/\left(4\sqrt{2\pi }\right)$ and those used in Eq. (A6) are $L=0.2$ and $C_0=0.25$. The noise variance in both cases is ${\sigma }_n^2=0.1$. The hatched area indicates the uncertainty of the reconstructed power spectrum estimated from the inverse Hessian, and the dotted horizontal line, the noise level.

Figure 7

Signal fields corresponding to the power spectra shown in Fig. 6. Top: A signal with a Gaussian two-point correlation function, Eq. (A4). Bottom: A signalwith a triangular correlation function, Eq. (A6). Shown in each plot are the signal realization $s$ [solid (blue) line, the data $d$ (crosses), the signal reconstruction $m$ [dashed (green) line], and an estimate of the local $1\sigma$ uncertainty of the reconstruction, given by $m\pm \mathrm{diag}{\left(D\right)}^{1/2}$ (hatched area). At the bottom of each panel, the same is plotted with the true signal subtracted.

Figure 8

One-dimensional log-normal reconstruction in a mildly nonlinear low-noise case. The power spectrum from which the signal realization was drawn is given by Eq. (14) with $P_0=0.1$, $k_0=5$, and $\gamma =4$. The noise variance is ${\sigma }_n^2=0.1$. Top: The exponentiated signal field $e^s$ [solid (blue) line], the data $d$ (crosses), the reconstruction $e^{m+\frac{1}{2}\mathrm{diag}\left(D\right)}$ [dashed (green) line], and the uncertainty interval of the reconstruction, given by $e^{m+\frac{1}{2}\mathrm{diag}\left(D\right)\pm {\left(\mathrm{diag}\left(D\right)\right)}^{1/2}}$ (hatched region). Bottom: The signal field $s$ [solid (blue) line], the logarithm of the data $\log \left(d\right)$ (crosses), the reconstruction $m$ [dashed (green) line], and the uncertainty interval for $m$, given by $m\pm {\left(\mathrm{diag}\left(D\right)\right)}^{1/2}$. In the lower panel, only data points for which $\log d$ is greater than $-1.5$ are shown.

Figure 9

One-dimensional log-normal reconstruction in a highly nonlinear low-noise case. Plotted quantities are the same as in Fig. 8. Parameters describing the power spectrum are $P_0=0.3$, $k_0=5$, and $\gamma =4$ and the noise variance is ${\sigma }_n^2=1$. In the lower panel, only data points for which $\log d$ is greater than $-3$ are shown.

Figure 10

One-dimensional log-normal reconstruction in a highly nonlinear high-noise case. Plotted quantities are the same as in Fig. 8. Parameters describing the power spectrum are $P_0=0.3$, $k_0=5$, and $\gamma =4$ and the noise variance is ${\sigma }_n^2=25$. In the lower panel, only data points for which $\log d$ is greater than $-3$ are shown.

Figure 11

Log-normal reconstruction example in a toroidal setting. The signal's power spectrum is given by Eq. (14) with $P_0=0.3$, $k_0=2$, and $\gamma =5$. The noise variance is ${\sigma }_n^2=10$. Top row: The data set $d$ (left) and its logarithm (right). In the logarithmic version, pixels with negative data values are plotted in dark blue. Middle row: The exponentiated signal field $e^s$ (left) and its reconstruction, given by $e^{m+\frac{1}{2}\mathrm{diag}\left(D\right)}$ (right). Bottom row: The corresponding nonexponentiated quantities $s$ (left) and $m$ (right).

Figure 12

Log-normal reconstruction example on a two-sphere. Quantities plotted in the top three rows are the same as in Fig. 11. Bottom row: The uncertainty of the signal reconstruction given by $\mathrm{diag}{\left(D\right)}^{1/2}$ (left) and the fractional uncertainty of the reconstruction of the exponentiated field given by $e^{\mathrm{diag}{\left(D\right)}^{1/2}}-1$ (right). The power spectrum parameters are $P_0=0.3$, $k_0=5$, and $\gamma =4$ and the noise variance is ${\sigma }_n^2=10$. In the top row, pixels without measurements are plotted in white.