Now that I have some simulations under way I have had time to take a look through the recent(ish) astro-ph posting, Matrix-free Large Scale Bayesian Inference in Cosmology, by Jasche & Lavaux, in which the authors propose a “messenger field” approach to the “Wiener sampler” allowing for simulation from a high-dimensional posterior without the need for large, computationally-expensive matrix inversions. Translating this into statistical terminology we have an “auxiliary variable” method for posterior simulation from the “Normal linear model” with correlated observational errors … which will be computationally efficient* if (and only if) the covariance function and domain of the latent Gaussian are amenable to FFT-based orthogonalization*. Indeed, FFT-based methods have a long history of providing powerful solutions for ‘matrix-free’ (i.e. matrix inversion free) sampling of high-dimensional Gaussian posteriors, particularly under the Normal linear model with uncorrelated observational errors and for latent Normal models with, e.g., log Poisson link functions: see, for instance, Wilke 2002, Paciorek 2007, Lang & Potthoff 2013, and on a related note the least squares trick of Simpson 2008. The key contribution of Jasche & Lavaux then is to identify an auxiliary variable decomposition to allow application of the FFT approach to the (less common) case of the Normal linear model with correlated observational error term.

Whether this approach will find much application beyond the specific problem(s) considered by the authors (and their collaborators) is unclear, although my gut feeling is that it’s rather unlikely. My understanding from previous discussions with cosmologists working in this field of large-scale Bayesian inference of the power spectrum parameters is that FFT techniques are already well-known—owing to the circulant structure of the key covariance matrices—but that their application is hampered for much the reasons as in (e.g.) contemporary geostatistical problems (e.g. Gething et al. 2010): the data space simply cannot be readily transformed to the regular grid structure required for FFT (or at least not without exceeding limitations on computational memory). Personally, I’ll still be recommending INLA to any astronomers asking me how to fit large-scale GP models …

### Like this:

Like Loading...

*Related*

Pingback: Batch job … | Another Astrostatistics Blog