## Learning GP covariance matrices by example …

Another astro ph submission that I’m reviewing a few weeks late here is this one by Garnett et al., entitled “Detecting Damped Lyman-$\alpha$ Absorbers with Gaussian Processes”.  The application here is pretty much what it says in the title; that is, probabilistic classification of QSO spectra into those containing DLAs and those without (plus, of course, an estimation of the redshift and column density of any candidate absorption systems).  Underneath the hood this is a straightforwards Bayesian model selection problem in which the underlying quasar spectrum is modelled as the realisation of a Gaussian process and any DLAs as Voigt profiles acting multiplicatively on top.

The strength of the analysis here is that a series of sensible decisions are made to yield a fast and effective implementation.  In particular, (i) a large training dataset of labelled SDSS quasar spectra is used to learn a reasonable mean and covariance function for quasars without DLAs, with the covariance represented via a low-rank approximation; (ii) the Woodbury identity and Sylvester determinant theorem are used to make efficient matrix operations on the Gaussian; and (iii) the two parameters of each absorber are integrated out numerically via a quasi-Monte Carlo ‘grid’ (FYI: Gerber & Chopin give some useful notes on the relative performance of MC and QMC methods at given dimension in the introduction to their QSMC read paper).  While it may not be common in purely observational data analysis problems to have a pre-labelled catalogue against which to learn a covariance function, it is often the case in astronomy that one has a reference set of simulated mock datasets against which a covariance can be learnt so this approach could well be useful in a number of other settings.  For example, to learn the (presumably non-stationary) spatial covariance function underlying a log-Gaussian Cox process of the galaxy distribution in a magnitude-limited survey.

Looking forwards, I wonder as to the validity of the assumed Gaussian likelihood (as proxy for Poisson photon counting statistics) as the method is rolled out to noisier spectra; of course, removing Gaussianity from the top layer of the likelihood function would mean analytic solutions are no longer available for key marginalisation steps, with presumably a massive increase in computational cost to the algorithm.

Minor notes:
– I was surprised to see an unweighted mean in Eqn 22 since I would have thought each quasar might have a different observational noise on its j-th wavelength bin for which usually one would use inverse variance weighting
– The term ‘unbiased’ is used improperly (from a statistical point of view) in describing the QMC integral of Eqn 57.