Abstract

Multi-output Gaussian processes (MOGPs) leverage the flexibility and interpretability of GPs while capturing structure across outputs, which is desirable, for example, in spatio-temporal modelling. The key problem with MOGPs is the cubic computational scaling in the number of both inputs (e.g., time points or locations), $n$, and outputs, $p$. Current methods reduce this to $O(n^3{m}^3)$, where is the desired degrees of freedom. This computational cost, however, is still prohibitive in many applications. To address this limitation, we present the Orthogonal Linear Mixing Model (OLMM), an MOGP in which exact inference scales linearly in $m$: $O(n^{3}m)$. This advance opens up a wide range of real-world tasks and can be combined with existing GP approximations in a plug-and-play way as demonstrated in the paper. Additionally, the paper organises the existing disparate literature on MOGP models into a simple taxonomy called the Mixing Model Hierarchy (MMH).

remresetThe remreset package *[parskip]latexCommand compat=newest calc, positioning, fit, tikzmark, arrows.meta, shapes, decorations.pathreplacing line/.style = thick, ->, > = Triangle[length=2.0mm, width=2.0mm] , arrow/.style = line , hidden node/.style = circle, minimum size = 0.9cm, draw = white, thick , latent node/.style = hidden node, draw = black, , factor node/.style = hidden node, rectangle, draw = black, , observed node/.style = latent node, fill = gray!15 , plate/.style = draw, label=[anchor=north west]south west:#1, rounded corners=2pt, shape=rectangle, inner sep=10pt, thick every picture/.append style = xscale = 2, yscale = 2 + assumptionAssumptionAssumptions assumptionAssumptionAssumptions corollaryCorollaryCorollaries corollaryCorollaryCorollaries definitionDefinitionDefinitions definitionDefinitionDefinitions exampleExampleExamples exampleExampleExamples factFactFacts factFactFacts lemmaLemmaLemmas lemmaLemmaLemmas modelModelModels modelModelModels propositionPropositionPropositions propositionPropositionPropositions questionQuestionQuestions questionQuestionQuestions remarkRemarkRemarks remarkRemarkRemarks theoremTheoremTheorems theoremTheoremTheorems asslistiAssumptionAssumptions asslistiAssumptionAssumptions corlistiCorollaryCorollaries corlistiCorollaryCorollaries deflistiDefinitionDefinitions deflistiDefinitionDefinitions exlistiExampleExamples exlistiExampleExamples factlistiFactFacts factlistiFactFacts lemlistiLemmaLemmas lemlistiLemmaLemmas modlistiModelModels modlistiModelModels proplistiPropositionPropositions proplistiPropositionPropositions qlistiQuestionQuestions qlistiQuestionQuestions remlistiRemarkRemarks remlistiRemarkRemarks thmlistiTheoremTheorems thmlistiTheoremTheorems sectionSec.Sects. propositionProp.Props. lemmaLem.Lems. modelMod.Mods. appendixApp.Apps.

 

Scalable Exact Inference in Multi-Output Gaussian Processes

 

Wessel P. Bruinsma${}^{\mathrm{1}\mathrm{,}\mathrm{2}}$ Eric Perim${}^{\mathrm{2}}$ Will Tebbutt${}^{\mathrm{1}}$ J. Scott Hosking${}^{\mathrm{3}}$ Arno Solin${}^4$ Richard E. Turner${}^{1,5}$

${}^1$University of Cambridge, ${}^2$Invenia Labs, ${}^3$British Antarctic Survey, ${}^4$Aalto University, ${}^5$Microsoft Research

1 INTRODUCTION

Application areas such as geostatistics, sensor networks, psychometrics, neuroscience, brain imaging, and finance all utilise multi-output statistical models. A rich, but disparate literature has deployed Gaussian processes (GPs, Rasmussen:2006:Gaussian_Processes) in such modelling problems, the aim being to leverage their interpretability, modularity, and tractability. Indeed, GPs have been successfully applied in a wide variety of single-output contexts: they can automatically discover structure in signals (Duvenaud:2014:Automatic_Construction), achieve state-of-the-art performance in regression tasks (Bui:2016:Deep_Gaussian_Processes_for_Regression), enable data-efficient models in reinforcement learning (Deisenroth:2011:PILCO_A_Model-Based_and_Data-Efficient), and support many applications in probabilistic numerics (Hennig:2015:Probabilistic_Numerics_and_Uncertainty_in), such as in optimisation (Brochu:2010:A_Tutorial_on_Bayesian_Optimization) and quadrature (Minka:2000:Quadrature_GP).

One of the first applications of GPs to multiple outputs was in geostatistics (Matheron:1969:Le_Krigeage_Universel). Today, multi-output GP (MOGP) models can be found in various areas, including geostatistics (Wackernagel:2003:Multivariate_Geostatistics), factor analysis (Teh:2005:Semiparametric_Latent_Factor; Yu:2009:Gaussian-Process_Factor_Analysis_for_Low-Dimensional), dependent or multi-task learning (Boyle:2005:Dependent_Gaussian_Processes; Bonilla:2007:Kernel_Multi-Task_Learning_Using_Task-Specific; Bonilla:2008:Multi-Task_Gaussian_Process; Osborne:2008:Towards_Real-Time_Information_Processing_of), latent force models (Alvarez:2009:Latent_Force_Models; Alvarez:2009:Sparse_Convolved_Gaussian_Processes_for; Alvarez:2010:Efficient_Multioutput_Gaussian_Processes_Through; Alvarez:2011:Computationally_Efficient_Convolved), state space modelling (Sarkka:2013:Spatiotemporal_Learning_via), regression networks (Wilson:2012:GP_Regression_Networks; Nguyen:2014:Collaborative_Multi-Output; Dezfouli:2017:Semi-Parametric_Network_Structure_Discovery_Models), and mixture models (Ulrich:2015:Cross_Spectrum; Bruinsma:2016:GGPCM; Parra:2017:Spectral_Mixture_Kernels_for_Multi-Output; Requeima:2019:The_Gaussian_Process_Autoregressive_Regression). Although many of these models are constructed in a similar spirit, they differ in their details and terminology, which makes them difficult to compare and organise.

A key practical problem with existing MOGPs is their computational complexity. For $n$ input points each having $p$ outputs, inference and learning in general MOGPs take $O(n^3{p}^3)$ time and $O(n^2{p}^2)$ memory, although these may be alleviated by a wide range of approximations (Quinonero:2005:Unifying_View; Titsias:2009:Variational_Learning; Lazaro-Gredilla:2010:Sparse_Spectrum_Gaussian_Process_Regression; Hensman:2013:Gaussian_Processes_for_Big_Data; Wilson:2015:Kernel_Interpolation_for_Scalable_Structured; Bui:2016:A_Unifying_Framework_for_Gaussian; Cheng:2017:Variational_Inference_for_Gaussian_Process; Hensman:2018:Variational_Fourier_Features_for_Gaussian). To mitigate these unfavourable scalings, a particular class of MOGPs, called the Linear Mixing Model (LMM, creftype 2.2), exploits low-rank structure of their covariance via the matrix inversion and determinant lemmas to reduce the complexity of inference and learning to $O(n^3{m}^3)$ time and $O(n^2{m}^2)$ memory, where $m$ is the desired degrees of freedom of the observations, typically much less than $p$. However, many multi-output regression problems cannot be solved adequately with only a small number of degrees of freedom, making the LMM prohibitively expensive in such cases. For example, consider temperatures measured across the Earth’s surface concatenated into a vector. Due to the complexity of the climate, we anticipate behaviour at many spatial length scales, so a model that supposes that this vector lives in a particularly low-dimensional “pancake”${}^{\text{1}}$ will not accurately reflect reality.

We present a unifying view of MOGPs to simplify the complicated web of models in this class. We call this the Mixing Model Hierarchy (MMH). Building on the Linear Mixing Model (LMM), which is the class of models at the bottom of the MMH, we introduce the Orthogonal Linear Mixing Model (OLMM). The OLMM is a particular LMM in which inference and learning take $O(n^{3}m)$ time and $O(n^{2}m)$ memory, without significantly sacrificing expressivity nor requiring any approximations. The OLMM is simple to implement and trivially compatible with existing single-output GP scaling techniques to scale to large numbers of data points $n$. In particular, for a large class of multi-output problems, e.g. many spatio-temporal problems, the time and memory complexity can be further reduced to $O(nm)$, which is linear in both the number of data points $n$ and degrees of freedom $m$. We demonstrate the OLMM in experiments on synthetic and real-world data sets.

[ scale=1.5, labelbg/.style= fill=white, inner sep=0pt, opacity=0.7, text opacity=1 ] [thick, draw=black!30] (0, 0, -1) coordinate (backsw) – ++(0.8, 0, 0) coordinate (backse); [thick] (backse) – ++(0, 0.8, 0) coordinate (backne) – ++(-0.8, 0, 0) coordinate (backnw); [thick, draw=black!30] (backnw) – (backsw); [thick] (0, 0, 0) coordinate (frontsw) – ++(1, 0, 0) coordinate (frontse) – ++(0, 1, 0) coordinate (frontne) – ++(-1, 0, 0) coordinate (frontnw) – cycle; [thick, draw=black!30] (frontsw) – (backsw); [thick] (frontse) – (backse) (frontnw) – (backnw) (frontne) – (backne); [arrow] ($(frontsw)+(0,-0.07,0)$) – ($(frontse)+(0,-0.07,0)$) node [pos=0.475, fill=white, inner sep=.5pt, ] convolutional; [arrow] ($(frontsw)+(-0.07,0.07,0)$) – ($(frontnw)+(-0.07,-0.07,0)$) node [pos=0.475, rotate=90, fill=white, inner sep=.5pt, ] time-varying; [arrow] ($(frontnw)+(-0.1,0,-0.17)$) – ($(backnw)+(-0.1,0.025,0)$) node [pos=0.475, rotate=22, align=center, anchor=south ] prior
on $H$; [anchor=west] (labelfrontsw) at ($(frontsw)+(-110:0.25)$) $f(t)=Hx(t)$; [anchor=west] (labelfrontnw) at ($(frontnw)+(110:0.45)$) $f(t)=H(t)x(t)$; [anchor=west, labelbg] (labelfrontne) at ($(frontne)+(25:0.35)$) $f(t)={\displaystyle \int H(t,\tau )x(\tau )d\tau }$; [anchor=west] (labelfrontse) at ($(frontse)+(-25:0.1)$) $f(t)={\displaystyle \int H(t-\tau )x(\tau )d\tau }$; [dashed] (labelfrontsw.west) – (frontsw); [dashed] (labelfrontse.west) – (frontse); [dashed] (labelfrontnw.west) – (frontnw); [dashed] (labelfrontne.west) – (frontne); [anchor=south west, align=left, labelbg, xshift=3pt, yshift=3pt] at (frontsw) [1–3, 5, 6,
[8, 17–19]; [anchor=south east, yshift=3pt, xshift=-3pt, labelbg] at (frontse) [4, 7, 9, 10]; [anchor=south west] at (backnw) [11, 15, 16]; [anchor=south east, color=black] at (backse) [13]; [anchor=south west] at (backsw) [14, 15]; [anchor=west, align=center, text width=7cm] (Hs) at (1.0, 0.565, 0) Form of $H$ Form of $K$ Mixing [1, 5] $H$ $kI$ Inst. [2] $\left[H_1\mathrm{\cdots }{H}_Q\right]$ $\mathrm{diag}(k_{1}I,\mathrm{\dots },k_{Q}I)$ Inst. [3, 6, 8, 17–19] $H$ $\mathrm{diag}(k_1,\mathrm{\dots },k_Q)$ Inst. [4, 9, 10, 13] $H(t-t^{\prime })$ $\mathrm{diag}(\delta ,\mathrm{\dots },\delta )$ Conv. [7] Green’s function $\mathrm{diag}(k_1,\mathrm{\dots },k_Q)$ Conv. [11, 16] $H(t)$ $\mathrm{diag}(k_1,\mathrm{\dots },k_Q)$ Inst. [12, 14] $\left[HI\right]$ $\mathrm{diag}(k_1,\mathrm{\dots },k_{Q+P})$ Inst. [15] Lower triangular $\mathrm{diag}(k_1,\mathrm{\dots },k_Q)$ Inst. ; [anchor=north east, align=center, text width=.25] at ($(4.5,1.35,0)$) Models:
0.9 [1] Intrinstic Coregionalisation Model (Goovaerts:1997:Geostatistics_for_Natural_Resources_Evaluation) [2] Linear Model of Coregionalisation (Goovaerts:1997:Geostatistics_for_Natural_Resources_Evaluation) [3] Semiparametric Latent Factor Model (Teh:2005:Semiparametric_Latent_Factor) [4] Dependent Gaussian Processes (Boyle:2005:Dependent_Gaussian_Processes) [5] Multi-Task Gaussian Processes (Bonilla:2008:Multi-Task_Gaussian_Process) [6] Osborne:2008:Towards_Real-Time_Information_Processing_of [7] Latent Force Models (Alvarez:2009:Latent_Force_Models) [8] Gaussian Process Factor Analysis (Yu:2009:Gaussian-Process_Factor_Analysis_for_Low-Dimensional) [9] Efficient Multioutput Gaussian Processes Through Variational Inducing Kernels (Alvarez:2010:Efficient_Multioutput_Gaussian_Processes_Through) [10] Convolved Multiple Output Gaussian Processes (Alvarez:2011:Computationally_Efficient_Convolved)

; [anchor=north west, align=center, text width=.25] at ($(4.2,1.35,0)$) Models:
0.9 [11] Gaussian Process Regression Network (Wilson:2012:GP_Regression_Networks) [12] Collaborative Multi-Output Gaussian Processes (Nguyen:2014:Collaborative_Multi-Output) [13] Generalised Gaussian Process Convolution Model (Bruinsma:2016:GGPCM) [14] Semi-Parametric Network Structure Discovery Models (Dezfouli:2017:Semi-Parametric_Network_Structure_Discovery_Models) [15] The Gaussian Process Autoregressive Regression Model (Requeima:2019:The_Gaussian_Process_Autoregressive_Regression) [16] Grouped Gaussian Processes (Dahl:2018:Grouped_Gaussian_Processes_for_Solar) [17] Spatio-temporal Bayesian filtering and smoothing (Sarkka:2013:Spatiotemporal_Learning_via) [18] Linear Mixing Model (creftype 1) [19] Orthogonal Linear Mixing Model (creftype 2)

;

Figure 1: The Mixing Model Hierarchy, which organises multi-output GPs from the machine learning and geostatistics literature according to their modelling assumptions. In the table, $k_i=k_i(t,t^{\prime })$ denotes a kernel, $\delta =\delta (t-t^{\prime })$ is the Dirac delta function, “inst.” abbreviates “instantaneous”, and “conv.” abbreviates “convolutional”.

2 MULTI-OUTPUT GP MODELS

For tasks with $p$ outputs, multi-output Gaussian processes induce a prior distribution over vector-valued functions $f:𝒯\to {ℝ}^p$ by requiring that any finite number of function values $f_{p_1}(t_1),\mathrm{\dots },f_{p_n}(t_n)$ with ${(p_i)}_{i=1}^n\subseteq \{1,\mathrm{\dots },p\}$ are multivariate Gaussian distributed. Often, the input space $𝒯=ℝ$ is time, or $𝒯={ℝ}^d$ is some feature space. An MOGP $f\sim 𝒢𝒫(m,K)$ is described by a vector-valued mean function

and a matrix-valued covariance function

$K(t,t^{\prime })=𝔼[f(t)f^{\text{𝖳}}(t^{\prime })]-𝔼[f(t)]𝔼[f^{\text{𝖳}}(t^{\prime })]$

For $n$ observations $y(t_1),\mathrm{\dots },y(t_n)\in {ℝ}^p$, inference and learning take $O(n^3{p}^3)$ time and $O(n^2{p}^2)$ memory.

2.1 Mixing Model Hierarchy

Our unifying Mixing Model Hierarchy begins with Linear Mixing Models (LMMs), a simple class of MOGP models typically characterised by low-rank covariance structure. An LMM decomposes a signal $f(t)$ comprising $p$ outputs into a fixed basis $h_1,\mathrm{\dots },h_m\in {ℝ}^p$ whose coefficients $x_1(t),\mathrm{\dots },x_m(t)$ are time-varying and modelled independently with Gaussian processes:

$f(t)$ $=h_1{x}_1(t)+\mathrm{\dots }+h_m{x}_m(t)$

The noisy signal $y(t)$ is, then, generated by adding $𝒩(0,\mathrm{\Lambda })$-distributed noise to $f(t)$. Intuitively, this means that the $p$-dimensional data lives in a “pancake” around the $m$-dimensional column space of $H$, where typically $m\ll p$.

Model 1 (Linear Mixing Model).

Let $K$ be an $m\times m$ diagonal multi-output kernel with $K(t,t)=I_m$, $H$ a $p\times m$ matrix, and $\mathrm{\Lambda }=\mathrm{diag}({\sigma }_1^2,\mathrm{\dots },{\sigma }_p^2)$. Then the LMM is given by the following generative model:

$x$ $\sim 𝒢𝒫(0,K(t,t^{\prime })),$ (latent processes)

$f(t)|H,x$ $=Hx(t),$ (likelihood)

$y|f$ $\sim 𝒢𝒫(f(t),\delta [t-t^{\prime }]\mathrm{\Lambda }).$ (noise model)

We call $x$ the latent processes and $H$ the mixing matrix or basis. If we marginalise out $f$ and $x$, we find $y\sim 𝒢𝒫(0,HK(t,t^{\prime })H^{\text{𝖳}}+\mathrm{\Lambda })$, which reveals the low-rank covariance structure of the LMM. It also shows that the LMM is a time-varying generalisation of factor analysis (FA): choosing $K(t,t^{\prime })=\delta [t-t^{\prime }]I_m$ recovers FA exactly.

The graphical model of the LMM is illustrated in the top-left corner of creftype 6 in creftype A, which highlights two restrictions of the LMM compared to a general MOGP: (i) the instantaneous spatial covariance of $f$, $𝔼[f(t)f^{\text{𝖳}}(t)]=H{H}^{\text{𝖳}}$, does not vary with time, because $H$ does not vary with time; and (ii) the noise-free observation $f(t)$ is generated from $x(t^{\prime })$ for $t^{\prime }=t$ only, meaning that, for example, $f$ cannot be $x$ with a delay or a smoothed version of $x$. We hence call the LMM a time-invariant (due to (i)) and instantaneous (due to (ii)) MOGP.

The LMM can be generalised in three ways. First, the mixing matrix $H$ may vary with time. Then $H\in {ℝ}^{p\times m}$ becomes a matrix-valued function $H:𝒯\to {ℝ}^{p\times m}$, and the likelihood becomes $f(t)|H,x=H(t)x(t)$. We call such MOGP models time-varying (see creftype 6, top-right). Second, $f(t)$ may depend on $x(t^{\prime })$ for all $t^{\prime }\in 𝒯$. Then the mixing matrix $H\in {ℝ}^{p\times m}$ becomes a matrix-valued time-invariant filter $H:𝒯\to {ℝ}^{p\times m}$, and the likelihood becomes $f(t)|H,x=\int H(t-\tau )x(\tau )d\tau$. We call such MOGP models convolutional (see creftype 6, bottom-left). Finally, $f(t)$ may depend on $x(t^{\prime })$ for all $t^{\prime }\in 𝒯$ and this relationship may vary with time. Then the mixing matrix $H\in {ℝ}^{p\times m}$ becomes a matrix-valued time-varying filter $H:𝒯\times 𝒯\to {ℝ}^{p\times m}$, and the likelihood becomes $f(t)|H,x=\int H(t,\tau )x(\tau )d\tau$. We call such MOGP models time-varying and convolutional (see creftype 6, bottom-right). The graphical models corresponding to these generalisations of the LMM are depicted in creftype 6 in the creftype A.

The LMM can be extended in one other way, which is to include a prior distribution on $H$. This extension and the three previously proposed generalisations together form the Mixing Model Hierarchy (MMH), which is depicted in creftype 1. The MMH organises multi-output Gaussian process models according to their distinctive modelling assumptions. creftype 1 shows how sixteen MOGP models from the machine learning and geostatistics literature can be recovered as special cases of the various generalisations of the LMM. Not all multi-output Gaussian process models are covered by the MMH, however. For example, Deep GPs (Damianou:2015:Deep_Gaussian_Processes_and_Variational) and variations thereon (Kaiser:2017:Bayesian_Alignments_of_Warped_Multi-Output) are excluded because they transform the latent processes nonlinearly to generate the observations.

2.2 Inference and Learning in the LMM

In MOGPs, the complexities of inference and learning can often be alleviated by exploiting structure in the kernel. This is also the case for the Linear Mixing Model: the $p$-dimensional observations have only $m$ degrees of freedom, where typically $m\ll p$, and it turns out that the $O(n^3{p}^3)$ time and $O(n^2{p}^2)$ memory complexities can be brought down to $O(n^3{m}^3)$ and $O(n^2{m}^2)$. In particular, in the following we show that the covariance structure of the LMM allows us to summarise the $p$-dimensional observations in an $m$-dimensional space, without loss of information.

We identify a low-dimensional projection of observations $y$ that is a sufficient statistic for the model and which can therefore be used to accelerate inference. This projection of $y$ is given by the maximum likelihood estimate (MLE) of $x$ under the likelihood of the LMM. Denote $p(y|x)=𝒩(y;Hx,\mathrm{\Lambda })$, and let $T$ be the $m\times p$ matrix ${(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}$. creftype 3 in creftype B shows that $Ty={\arg \max }_{x}p(y|x)$ and that $Ty$ is an unbiased estimate of $x$. Hence, $Ty$ is a MLE for $x$, which means that it is a function of a sufficient statistic for $x$, if one exists. The following proposition shows that $Ty$ is actually minimally sufficient itself.

Prop. 1.

The MLE $Ty$ of $x$ is a minimal sufficient statistic for $x$. Proof. See creftype C.

For any prior $p(x)$ over $x$, sufficiency of $T$ gives that $p(x|y)=p(x|Ty)$; i.e., conditioning on $y$ is equivalent to conditioning on $Ty$, where $Ty$ can be interpreted as a “summary” or “projection” of $y$. This idea is formalised in the following proposition, which is an immediate consequence of creftypeplural 11 in creftype C.

Prop. 2.

Let $p(x)$ be a model for $x:𝒯\to {ℝ}^m$, not necessarily Gaussian, $H$ a $p\times m$ matrix, and $\mathrm{\Lambda }=\mathrm{diag}({\sigma }_1^2,\mathrm{\dots },{\sigma }_p^2)$. Then consider the following generative model:

$x$ $\sim p(x),$ (latent processes)

$f(t)|H,x$ $=Hx(t),$ (likelihood)

$y|f$ $\sim 𝒢𝒫(f(t),\delta [t-t^{\prime }]\mathrm{\Lambda }).$ (noise model)

Consider a $p\times n$ matrix $Y$ of observations of $y$. Then $p(f|Y)=p(f|TY)$. Furthermore, the distribution of the projected observed signal $Ty$ is given by

$$Ty|x\sim 𝒢𝒫(x(t),\delta [t-t^{\prime }]{\mathrm{\Lambda }}_T)\text{ with }{\mathrm{\Lambda }}_T={(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1},$$

and the likelihood of $Y={(y_i)}_{i=1}^n\subseteq {ℝ}^p$ is given by

$$p(Y)=\left[\prod _{i=1}^n\frac{𝒩(y_i;0,\mathrm{\Lambda })}{𝒩(y_i;0,{\mathrm{\Lambda }}_T)}\right]\int p(x)\prod _{i=1}^n𝒩(T{y}_i;x_i,{\mathrm{\Lambda }}_T)\mathrm{d}x.$$

Crucially, $Y$ are $p$-dimensional observations, $TY$ are $m$-dimensional summaries, and typically $m\ll p$, so conditioning on $TY$ is often much cheaper. In particular, if we apply creftype 2 to the LMM by letting $x\sim 𝒢𝒫(0,K(t,t^{\prime }))$, we immediately get the claimed reduction in complexities: whereas conditioning on $Y$ takes $O(n^3{p}^3)$ time and $O(n^2{p}^2)$ memory, we may equivalently condition on $TY$, which takes $O(n^3{m}^3)$ time and $O(n^2{m}^2)$ memory instead. This important observation is depicted in creftype 1(a).

In creftype 2, we call ${\mathrm{\Lambda }}_T=T\mathrm{\Lambda }{T}^{\text{𝖳}}={(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}$ the projected observation noise. The noise term ${\mathrm{\Lambda }}_T$ is important, because it couples the latent processes upon observing data. In particular, if the latent processes are independent under the prior and ${\mathrm{\Lambda }}_T$ is diagonal, then the latent processes remain independent when data is observed. This observation forms the basis of the computational gains achieved by the Orthogonal Linear Mixing Model, which we present next.

[anchor=center] (Y) $Y$; [below=1.5cm of Y] (TY) $TY$; [right=2.5cm of Y] (pyY) $p(f|Y)$; [below=1.5cm of pyY] (pTY) $p(f|TY)$; [arrow] (Y) – (TY); [arrow] (Y) – node [pos=0.5, anchor=south] inference node [pos=0.5, anchor=north] $O(n^3{p}^3)$ (pyY); [arrow] (TY) – node [pos=0.5, anchor=south] inference node [pos=0.5, anchor=north] $O(n^3{m}^3)$ (pTY); [arrow] (pTY) – (pyY); (-.5cm,.5cm) rectangle (2.25cm,-1.5cm); (a) Conditioning on $Y$ in the LMM

[anchor=center] (Y) $Y$; [below=1.5cm of Y, anchor=center] (TY) $TY$; [right=3cm of Y, anchor=center] (pyY) $p(f|Y)$; [right=3cm of TY, yshift=.75cm, anchor=center] (pu1TY) $p(x_1|{(TY)}_{1:})$; [right=3cm of TY, yshift=-.75cm, anchor=center] (pumTY) $p(x_m|{(TY)}_{m:})$; [arrow] (Y) – (TY); [arrow] (Y) – node [pos=0.5, anchor=south] inference node [pos=0.5, anchor=north] $O(n^3{p}^3)$ (pyY); [arrow] (TY) – node [pos=0.5, anchor=south, rotate=20] inference node [pos=0.5, anchor=north, rotate=20] $O(n^3)$ (pu1TY.west); [right=3cm of TY, anchor=center, yshift=.1cm] $\mathrm{⋮}$; [arrow] (TY) – node [pos=0.5, anchor=south, rotate=-20] inference node [pos=0.5, anchor=north, rotate=-20] $O(n^3)$ (pumTY.west); [arrow, ->] (pu1TY.east) to[bend right=40] ([yshift=-2pt]pyY.east); [arrow, ->] (pumTY.east) to[out=60, in=-15] ([yshift=2pt]pyY.east); (-.25cm,.5cm) rectangle (2.5cm,-1.5cm); (b) Conditioning on $Y$ in the OLMM

[ scale=1.5, arrowlabel/.style= fill=white, opacity=0.7, text opacity=1, inner sep=0 , rotate=-30 ] [black] (0, 0) circle (0.02); at (0, 0) [anchor=north east] $0$; [dashed, thick] plot [smooth, tension=0.5] coordinates (0, 0.75) (0.125, 0.875) (0.25, 1.1) (0.375, 1.25) (0.5, 1.38) (0.75, 1.58) (1, 1.625) (1.25, 1.55) (1.375, 1.45) (1.45, 1.25) (1.375, 1) (1.25, 0.9) (1, 0.95) (0.75, 0.75) (0.8, 0.5) (1, 0.25) (1.05, 0) ; [black] (1.0, 0.95) circle (0.035) node [arrowlabel, anchor=south west, yshift=0.3em, xshift=0.5em] $f(t)=(f_1(t),f_2(t))$; [mplorange, arrow] (0, 0) – node [midway, anchor=north west, arrowlabel, yshift=-0.1em, xshift=-0.1em] $h_1{x}_1(t)$ (1.175, 0.55); [mplorange, arrow] (1.175, 0.55) – node [pos=.5, anchor=west, arrowlabel, xshift=0.2em] $h_2{x}_2(t)$ (1.0, 0.95); [thick, mplorange!50] (1.084, 0.508) – (1.044, 0.599) – (1.135, 0.642); [mplblue, arrow] (0, 0) – node [pos=0.75, anchor=east, arrowlabel, xshift=-0.2em, yshift=0.2em] $h_1{x}_1(t)$ (1.2, 1.9); [mplblue, arrow] (1.2, 1.9) – node [pos=.5, anchor=west, arrowlabel, xshift=0.2em] $h_2{x}_2(t)$ (1.0, 0.95); (c) Difference between LMM and OLMM

Figure 2: (a–b) Commutative diagram depicting that conditioning on $Y$ in the LMM and in the OLMM is equivalent to independently conditioning every $x_i$ on ${(TY)}_{i:}$, but yields different computational complexities. (c) Illustration of the difference between the LMM and OLMM. The trajectory of a particle in two dimensions is modelled by the LMM ( blue) and OLMM ( orange). The noise-free position $f(t)$ is modelled as a linear combination of basis vectors $h_1$ and $h_2$ with coefficients $x_1(t)$ and $x_2(t)$ (two independent GPs). In the OLMM, the basis vectors $h_1$ and $h_2$ are constrained to be orthogonal, whereas in the LMM, $h_1$ and $h_2$ are unconstrained.

3 THE ORTHOGONAL LINEAR MIXING MODEL

The Orthogonal Linear Mixing Model (OLMM) is the special case of the LMM where the basis, the columns of $H$, is orthogonal. In particular, $H=U{S}^{\frac{1}{2}}$ where $U$ is a truncated orthogonal matrix and $S>0$ a diagonal scaling. We define this model as follows.

Model 2 (Orthogonal Linear Mixing Model).

Let $K$ be an $m\times m$ diagonal multi-output kernel with $K(t,t)=I_m$, $H$ a $p\times m$ matrix of the form $H=U{S}^{\frac{1}{2}}$ with $U$ a truncated orthogonal matrix and $S>0$ diagonal, and $\mathrm{\Lambda }={\sigma }^2{I}_p+HD{H}^{\text{𝖳}}$ a $p\times p$ matrix with $D>0$ diagonal. Then the OLMM is given by the following generative model:

$x$ $\sim 𝒢𝒫(0,K(t,t^{\prime })),$ (latent processes)

$f(t)|H,x$ $=Hx(t),$ (likelihood)

$y|f$ $\sim 𝒢𝒫(f(t),\delta [t-t^{\prime }]\mathrm{\Lambda }).$ (noise model)

The difference between the LMM and the OLMM is illustrated in creftype 1(c). In the OLMM, the number of latent processes must be smaller than or equal to the number of outputs, since the number of $p$-dimensional vectors that can be mutually orthogonal is at most $p$. Also, $D$ in $\mathrm{\Lambda }$ can be interpreted as heterogeneous noise deriving from the latent processes. Moreover, although $H$ and $\mathrm{\Lambda }$ do not depend on time, our analysis and results also apply to the case where $H_t$ and ${\mathrm{\Lambda }}_t$ do vary with time. Finally, for the OLMM, creftype 5 in creftype E shows that $T=S^{-\frac{1}{2}}{U}^{\text{𝖳}}$ and ${\mathrm{\Lambda }}_T={\sigma }^2{S}^{-1}+D.$

Whereas the LMM is a time-varying generalisation of FA, the OLMM is a time-varying generalisation of probabilistic principal component analysis (PPCA, Tipping:1999:Probabilistic_Principal_Component_Analysis): $D=0$ and $K(t,t^{\prime })=\delta [t-t^{\prime }]I_m$ recovers PPCA exactly. The OLMM is also similar to Gaussian Process Factor Analysis (GPFA, Yu:2009:Gaussian-Process_Factor_Analysis_for_Low-Dimensional), with the crucial difference being that in GPFA orthogonalisation of the columns of $H$ is done as a post-processing step, whereas in the OLMM orthogonality of the columns of $H$ is built into the model.

Generality of the OLMM.

In the separable case, where $K(t,t^{\prime })=k(t,t^{\prime })I_m$ for a scalar-valued kernel $k$, for every LMM with homogeneous observation noise ($\mathrm{\Lambda }={\sigma }^2{I}_p$), there exists an OLMM with $D=0$ that is equal in distribution of $y$. To see this, note that

$y^{\text{(LMM)}}$ $\sim 𝒢𝒫(0,k(t,t^{\prime })H{H}^{\text{𝖳}}+{\sigma }^2\delta [t-t^{\prime }]I_p),$

$y^{\text{(OLMM)}}$ $\sim 𝒢𝒫(0,k(t,t^{\prime })US{U}^{\text{𝖳}}+{\sigma }^2\delta [t-t^{\prime }]I_p).$

Hence, letting $US{U}^{\text{𝖳}}$ be the eigendecomposition of $H{H}^{\text{𝖳}}$ gives an OLMM equal in distribution of $y$. In the nonseparable case, where diagonal elements of $K$ are linearly independent, in general only the distribution of $y(t)$ at every $t$ can be recovered by an OLMM, but the correlation between $y(t)$ and $y(t^{\prime })$ for $t^{\prime }\ne t$ may be different. In terms of the joint distribution over $x$ and $y$, which is important for interpretability of the latent processes, creftype G shows that the Kullback–Leibler (KL) divergence between two LMMs with bases $H$ and $\widehat{H}$ is proportional to ${\parallel H-\widehat{H}\parallel }_F^2$, hence symmetric, where $\parallel \bullet {\parallel }_F$ denotes the Frobenius norm. As a consequence (creftype G), the KL between an LMM with basis $H$ and the OLMM closest in KL is upper bounded by ${\parallel I-V\parallel }_F^2$ where $V$ are the right singular vectors of $H$. This makes sense: $V=I$ implies that $H$ is of the form $U{S}^{\frac{1}{2}}$ with $U$ a truncated orthogonal matrix and $S>0$ diagonal. It also shows that an LMM is close to an OLMM if $V$ is close to $I$ in the sense of the Frobenius norm.

Choice of $H$.

The basis $H$ can be initialised and learned in several ways. (i) $H$ can be initialised randomly and learned through gradient-based optimisation of the marginal likelihood. (ii) Using that $𝔼[f(t)f^{\text{𝖳}}(t)]=H{H}^{\text{𝖳}}$, the basis $H$ can be initialised to $\widehat{U}{\widehat{S}}^{\frac{1}{2}}$ where $\widehat{\mathrm{\Sigma }}=\widehat{U}\widehat{S}{\widehat{U}}^{\text{𝖳}}$ is the eigendecomposition of an estimate $\widehat{\mathrm{\Sigma }}$ of the spatial covariance. (iii) In the case that there is a kernel over the outputs, e.g. in separable spatio-temporal GP models, $H$ can be set to $U{S}^{\frac{1}{2}}$ where $US{U}^{\text{𝖳}}$ is an eigendecomposition of the kernel matrix over the outputs. The hyperparameters of the kernel over the outputs can then be learned with gradient-based optimisation by differentiating through the eigendecomposition. We discuss this choice in more detail in creftype 4.

Diagonal projected noise.

As alluded to in creftype 2.2, the most important property of the OLMM is that the projected noise ${\mathrm{\Lambda }}_T$ from creftype 2 is diagonal: ${\mathrm{\Lambda }}_T={\sigma }^2{S}^{-1}+D.$ creftype 4 in the creftype D characterises exactly when this is the case. Specifically, creftype 4, says that the projected noise is diagonal if and only if $H$ is of the form $H={\mathrm{\Lambda }}^{\frac{1}{2}}U{S}^{\frac{1}{2}}$ with $U$ a truncated orthogonal matrix and $S>0$ diagonal. This condition is awkward, as it couples $H$ and $\mathrm{\Lambda }$. Fortunately, creftype 4 also shows that we may drop $H$’s dependency on $\mathrm{\Lambda }$ if and only if every column of $U$ is an eigenvector of $\mathrm{\Lambda }$. In the OLMM, every column of $U$ is indeed an eigenvector of $\mathrm{\Lambda }$. We therefore identify the OLMM as a general parametrisation of the LMM for which the projected noise is diagonal. Compared to the LMM, this brings computational advantages, which we discuss next.

Inference.

Since the projected noise is diagonal, the latent processes remain independent when data is observed. We may hence treat the latent processes independently, conditioning the $i$^th latent process $x_i$ on ${(TY)}_{i:}=Y^{\text{𝖳}}{u}_i/\sqrt{S_{ii}}$ under noise ${({\mathrm{\Lambda }}_T)}_{ii}={\sigma }^2/S_{ii}+D_{ii}$, where $u_i$ is the $i$^th column of $U$. Therefore, inference takes $O(n^{3}m)$ time and $O(n^{2}m)$ memory, which are linear in the number of degrees of freedom in the data $m$. This decoupled inference procedure is depicted in creftype 1(b). Remark that the decoupled problems can be treated in parallel to achieve sublinear wall time, and that in the separable case further speedups are possible.

Learning.

For computing the likelihood of observed data, the OLMM also offers computational benefits. As the projected noise is diagonal, by absorbing $D$ into $K$, the expression for $p(Y)$ from creftype 1 further simplifies to the following:

$p(Y)$ $=𝒩(\mathrm{vec}Y;0,{\sigma }^2{I}_{np})$

$\times {\displaystyle \prod _{i=1}^m}{\displaystyle \frac{𝒩({(TY)}_{i:};0,K_i+({\sigma }^2/S_{ii}+D_{ii})I_n)}{𝒩({(TY)}_{i:};0,({\sigma }^2/S_{ii})I_n)}},$

where $K_i$ is the $n\times n$ kernel matrix for the $i$^th latent process $x_i$. We conclude that learning also takes $O(n^{3}m)$ time and $O(n^{2}m)$ memory, again linear in the number of latent processes. Furthermore, we observe that $p(Y)$ is convex in $U$, which significantly aids learning; in comparison, for the LMM, $p(Y)$ is highly nonconvex in the mixing matrix $H$.

Interpretability.

Besides computational benefits, that the OLMM breaks down into independent problems for the latent processes also promotes interpretability. Namely, the independent problems can be separately inspected to interpret, diagnose, and improve the model. This is much easier than directly working with predictions for the data, which are high dimensional and often strongly correlated between outputs. In particular, if $H$ is sufficiently rich and the important independent problems are solved sufficiently well, then the OLMM achieves good mean squared error (MSE) for the data:

$$\underset{\text{MSE}}{\underset{﹇}{{\parallel y-Hx\parallel }^2}}=\underset{\begin{array}{c}\text{variance not}\\ \text{captured by basis}\end{array}}{\underset{﹇}{{\parallel P_{U^{⟂}}y\parallel }^2}}+\sum _{i=1}^m{S}_{ii}\underset{\begin{array}{c}\text{MSE of}\\ i\text{<sup class="ltx\_sup"><span class="ltx\_text" style="font-size:70\%;">th</span></sup>}\text{ latent process}\end{array}}{\underset{﹇}{{({(Ty)}_i-x_i)}^2}}$$

where $P_{U^{⟂}}$ is the orthogonal projection onto the orthogonal complement of $U$. See creftype F for a proof.

Scaling.

For both learning and inference, the problem decouples into $m$ independent single-output problems. Therefore, to scale to a large number of data points $n$, off-the-shelf single-output GP scaling techniques can be trivially applied to these independent problems. For example, if the variational inducing point method by Titsias:2009:Variational_Learning is used with $r\ll n$ inducing points, then inference and learning are further reduced to $O(n{r}^{2}m)$ time and $O(nrm)$ memory. Most importantly, if $k(t,t^{\prime })$ is Markovian (e.g. of the Matérn class) one can leverage state-space methods to exactly solve the $m$ independent problems in closed-form (Hartikainen:2010:Kalman_filtering_and_smoothing_solutions; Sarkka:2019:Applied_Stochastic_Differential_Equations). This brings down the scaling to $O(nm{d}^3)$ time and $O(nm{d}^2)$ memory, where $d$ is the state dimension (typically $d\ll m,n$). We further discuss this approach in creftype 4.

Missing data and heterogeneous observation noise.

Missing data is troublesome for the OLMM, because it is not possible to take away a subset of the rows of $H$ and retain orthogonality of the columns. Possible methods to circumvent this issue are discussed in creftype H, but are not explored further in this work. Furthermore, the specification of the OLMM does not allow for heterogeneous observation noise. This issue is discussed in creftype I, which shows that the OLMM can have automatic whitening of the data built in to allow for heterogeneous observation noise.

4 APPLICATION TO SEPARABLE SPATIO-TEMPORAL MODELS

Separable spatio-temporal GPs, which are of the form $f\sim 𝒢𝒫(0,k_t(t,t^{\prime })k_r(r,r^{\prime }))$, fall neatly within the OLMM framework if observed at a fixed number of locations in space. Specifically, the vector-valued process $f(t)={(f(t,r_i))}_{r=1}^p$ has distribution $𝒢𝒫(0,k_t(t,t^{\prime })K_r)$ where $K_r$ is the $p\times p$ matrix with ${(K_r)}_{ij}=k_r(r_i,r_j)$. Letting $US{U}^{\text{𝖳}}$ be the eigendecomposition of $K_r$, $f(t)$ is an OLMM with $H=U{S}^{\frac{1}{2}}$, which we note is square, and $K(t,t^{\prime })=k_t(t,t^{\prime })I_p$. Note that, in this case, projecting the data takes $O(n{p}^2)$ time. By relaxing $K$ to be a general diagonal multi-output kernel with $K(t,t)=I_p$, we obtain a nonseparable relaxation of the original separable model whilst retaining exact, computationally efficient inference. In particular, the orthogonal basis for the OLMM are the eigenvectors of a kernel matrix whose hyperparameters can be optimised.

Combining the OLMM framework with efficient state-space scaling techniques (Hartikainen:2010:Kalman_filtering_and_smoothing_solutions; Sarkka:2019:Applied_Stochastic_Differential_Equations; Solin:2018:Infinite-Horizon_Gaussian_Processes; nickisch2018state), which are either exact or arbitrarily good approximations, the complexities are reduced to $O(np+p^3)$ time and $O(np+p^2)$ memory for the entire problem, which are linear in $n$. This compares favourably with the filtering approach by Sarkka:2013:Spatiotemporal_Learning_via, which takes $O(n{p}^3)$ time and $O(n{p}^2)$ memory. The sparse approximation by hartikainen2011sparse provides a reduced-rank approximation with $O(n{q}^3)$ time and $O(n{q}^2)$ memory for some $q\ll p$, although it comes at the expense of accuracy and yields no asymptotic savings if $q$ must increase linearly in $p$ to maintain a reasonable accuracy. It also compares favourably with the Kronecker product decomposition (saatcci2012scalable, Ch. 5) approach, which requires $O(p^3+n^3)$ time and $O(p^2+n^2)$ memory complexity.

5 EXPERIMENTS

We test the OLMM in experiments on synthetic and real-world data sets. An implementation is available at https://github.com/wesselb/olmm.

Computational scaling.

We demonstrate that exact inference scales favourably in $m$ for the OLMM, whereas the LMM quickly becomes computationally infeasible as $m$ increases. We use a highly optimised implementation of exact inference for the LMM, kindly made available by Invenia Labs. creftype 3 shows the runtime and the memory usage of the LMM and OLMM. Observe that the LMM scales $O(m^3)$ in time and $O(m^2)$ in memory, whereas the OLMM scales $O(m)$ in both time and memory. For $m=25$, the LMM takes nearly 10 minutes to compute the evidence, whereas the OLMM only requires a couple seconds.

Temperature extrapolation.

We consider a simple spatio-temporal temperature prediction problem over Europe. Approximately 30 years worth of the ERA-Interim reanalysis temperature data (dee2011era) is smoothed in time with a Hamming window of width 31 and sub-sampled once every 31 days to produce a data set comprising $13\times 19=247$ outputs and approximately $350$ months worth of data. We train the OLMM and independent GPs (both models use Matérn–$5/2$ kernels with a periodic component) on the first 250 months of the data and test on the next 100 months. The basis for the OLMM is given by the eigenvectors of the kernel matrix over the points in space (Matérn–$5/2$ with a different length scale for latitude and longitude). creftype 1 shows the log-marginal likelihood (LML) of the training data, the posterior predictive log-probability (PPLP) and the root-mean-square error (RMSE) of the test data. Although the RMSE is equal to that of the independent GPs—the data is highly periodic and the predictions are highly accurate for both models—the LML and PPLP are much better, which indicates that that the OLMM is able to capture meaningful structure between temperatures at different points in space.

	LML	PPLP	RMSE
IGP	$-0.247$	$-0.219$	$1.99$
OLMM	$1.22$	$1.19$	$1.99$

Table 1: Log-marginal likelihood (LML) of the training data, posterior predictive log-probability (PPLP) and root-mean-square error (RMSE) of held-out test data of independent GPs (IGP) and the OLMM for the temperature extrapolation experiments.

Figure 3: Runtime (left) and memory usage (right) of the LMM and OLMM for computing the evidence of $n=1500$ observations for $p=200$ outputs

every picture/.append style = xscale = .5, yscale = .5 width=height=0.5xlabel=Number of bins $n$ along $x_1$, ylabel=${\displaystyle \frac{\text{Time Kron.}}{\text{Time OLMM}}}$, grid=major, xmin=5, xmax=15000, major grid style=densely dashed,black!25, legend style = at=(0.05, 0.95),anchor=north west, legend style=legend cell align=left, align=left, draw=none, fill=black!10, fill opacity=0.35, text opacity=1, draw opacity=1, scale only axis, axis on top, axis x line*=bottom, axis y line*=left, every axis plot/.append style=thick, axis line style = thick, xtick align=outside, xtick style=black, thick, ytick align=outside, ytick style=black, thick, clip=false [xmode = log] +[color=mycolor1, mark options=fill=mycolor1, opacity=1, mark size=1.5pt] coordinates (10000.0, 4.5624386332067655) (2000.0, 3.7446352996685253) (1000.0, 3.1158154627754544) (200.0, 2.0277622456016386) (100.0, 1.6077527986229443) (40.0, 1.0629959708548498) (20.0, 0.5306449960733431) (10.0, 0.18566913734215784) ; LML +[color=mycolor2, mark options=fill=mycolor2, opacity=1, mark size=1.5pt] coordinates (10000.0, 4.879464511037345) (2000.0, 3.349257213603934) (1000.0, 2.983019409435243) (200.0, 1.8322661062389873) (100.0, 1.3563551273191006) (40.0, 0.9694708480653853) (20.0, 0.5251749491862987) (10.0, 0.18835279833487512) ; RNG [color=black, dashed] coordinates (15000.0, 1) (5.0, 1) ;

Figure 4: Ratio of timings of the Kronecker approach (saatcci2012scalable, Ch. 5) and the OLMM to compute the log marginal likelihood (LML) and generate a single prior sample (RNG)

(a) Empirical and learned correlations

(b) Clustering

(c) Latent processes

Figure 5: Results of the large-scale climate simulator experiment, showing (a) the empirical and learned correlation between the simulators, (b) a dendrogram deriving from hierarhically clustering the simulators based on the learned correlations, and (c) predictions for the latent processes for the first two eigenvectors of the covariance matrix between simulators $i_s=1,2$ and the first two eigenvectors of the spatial covariance $i_r=1,2$, for the last 1000 days.

Rainforest tree point process modelling.

We consider a subset of the extensive rainforest data set credited to rainforest; Condit:2005; Hubbell:1999. The data features a 1000$$m $\times$ 500$$m rainforest dynamics plot in Barro Colorado Island, Panama. In the survey area, the locations of all Trichilia tuberculata (a tree species of the Mahogany family) have been measured (see creftype 7 in creftype J). We tackle this spatial point pattern with a log-Gaussian Cox process model, which is an inhomogeneous Poisson process model for count data. The unknown intensity function $\lambda (x)$ is modelled with a Gaussian process such that $f(x)=\log \lambda (x)$. We model locally constant intensity in subregions by discretising the region into $np$ bins (Moller+Syversveen+Waagepetersen:1998). This leads to having a Poisson observation model in each bin (see creftype J). This model reaches posterior consistency in the limit of bin width going to zero (Tokdar+Ghosh:2007). The accuracy thus improves with tighter binning. We use a separable Matérn-$5/2$ GP prior over $f(x_1,x_2)$, and discretise the area into a $n\times p=200\times 100$ grid with $np=20000$ grid bins, and treat the first dimension as time.

We use Elliptical Slice Sampling (Murray:2010:Elliptical_Slice_Sampling; murray2010slice) and Metropolis–Hastings (hastings1970monte) to perform inference in this problem, which require many evaluations of the log marginal likelihood and draws from the prior. The proposed methodology renders these operations much less computationally intensive than the general case. We find approximately an hour is required to perform ${10}^5$ MCMC iterations on a MacBook Pro. See creftype J for more details.

As previously discussed, the Kronecker product factorisation techniques (saatcci2012scalable, Ch. 5) are applicable in separable spatio-temporal scenarios. While our method in combination with state-space inference offers superior performance asymptotically, we demonstrate that we also obtain benefits on practical problems. creftype 4 shows for different number of bins $n$ (with $p=n/2$) the time for log marginal likelihood evaluation and sampling from the GP prior for our method and for the Kronecker product decomposition approach. Observe that the linear scaling of our methodology quickly renders it competitive. While this simple experiment does not provide a comprehensive assessment of the relative merits of these techniques—for example, only a single CPU core was used in this experiment—it clearly establishes the practical relevance of our technique in the separable setting.

Large-scale climate model calibration.

In this final experiment, we consider a large-scale climate modelling task over Europe. In particular, we use the OLMM to find relationships between 28 climate simulators4 (see taylor2012overview, for background) by letting $H=H_s\otimes H_r$, where $H_s$ are the first $m_s=5$ eigenvectors of a $28\times 28$ covariance matrix $K_s$ between the simulators, and $H_r$ are the first $m_r=10$ eigenvectors of the kernel matrix over the points in space (Matérn–$5/2$ with a different length scale for latitude and longitude). This means the $m_r{m}_s=50$ latent processes are indexed by two indices $i_s$ and $i_m$, one corresponding to the eigenvector of the simulator covariance and one correponding to the eigenvector of the spatial covariance. In the above, the covariance matrix $K_s$ between the simulators is a parameter that we learn. The kernels for the latent processes are Matérn–$5/2$. We consider $n=10000$ time points for all 28 simulators, each with 247 outputs, giving a total of roughly 70 million data points. For the independent problems, we use the variational inducing point method by Titsias:2009:Variational_Learning, where the positions of the inducing points are initialised to one every two months. All hyperparameters and the locations of the inducing points are optimised until convergence using scipy’s implementation of the L-BFGS-B algorithm (Nocedal:2006:Numerical_Optimisation), which takes about 4 hours on a MacBook Pro5. The learned length scales were ${23.3}^{\circ }$ for latitude and ${43.6}^{\circ }$ for longitude.

creftype 4(a) shows the empirical correlations and the correlations learned by the OLMM (derived from $K_s$). Observe that the empirical correlation is not informative, as it ignores all temporal structure. The OLMM, however, does not, and is consequently able to identify a rich correlation structure. In order to get insight into the learned correlations, we hierarchically cluster the models using farthest point linkage with $1-|\text{corr.}|$ as the distance. creftype 4(b) shows the resulting dendrogram, in which models are grouped by their similarity. For two models, the further to the right the branch connecting them is, the less similar the models are. In creftypeplural 4(b)4(a), HadGEM2 is clearly singled out: it is one of the simplest models, not including several processes that can be found in others, such as ocean & sea-ice, terrestrial carbon cycle, stratosphere, and ocean biogeochemistry (bellouin2011hadgem2). Furthermore, if we inspect the names of the simulators in the groups in creftype 4(b), we observe that often simulators of the same family are grouped together. We observe some interesting cases. (i) Although IPSL-CM5A-LR and IPSL-CM5A-MR are close, IPSL-CM5B-LR is grouped far apart. It turns out that IPSL-CM5A-LR and IPSL-CM5A-MR are different-resolution versions of the same model, while IPSL-CM5B-LR employs a different atmospheric model. (ii) ACCESS1.0 and ACCESS1.3 have a similar name, but differ greatly in their implementation: ACCESS1.0 is the basic model, while ACCESS1.3 is much more aspirational, including experimental atmospheric physics models and a particular land surface model (bi2013access). (iii) The distance between BCC_CSM1.1(m) and BCC_CSM1.1 can be explained by the more realistic surface air temperature predictions obtained by the former (wu2014overview), which is exactly the quantity we study. Finally, creftype 4(c) shows predictions for four latent processes ($i_s=1,2$ with $i_r=1,2$). The first spatial eigenvector ($i_r=1$) is constant in space; combined with the strongest eigenvector of $K_s$ ($i_s=1$), we obtain a strong signal constituting seasonal temperature changes.

6 DISCUSSION AND CONCLUSION

In this paper, we introduced the Orthogonal Linear Mixing Model as a multi-output GP in the Linear Mixing Model class where exact inference and learning scales linearly in the desired number $m$ of degrees of freedom of the observations. The OLMM is simple to implement and trivially compatible with off-the-shelf single-output GP scaling techniques to scale to large numbers of observations $n$. We introduced the Mixing Model Hierarchy, which organises MOGPs according to their generative models, and placed the OLMM within this scheme.

The OLMM is exciting for several reasons. It provides a way to do efficient inference in separable spatio-temporal models that is trivially compatible with single-output scaling techniques. Exact inference is linear in number of outputs, and, with state space methods, it is also linear in number of time points. The OLMM even extends separable spatio-temporal models to non-separable models whilst retaining exact, computationally efficient inference. Furthermore, the OLMM is interpretable: the latent processes are independent, which means that they can be inspected separately to diagnose the model. Finally, as demonstrated in the real-world experiments in the paper, the OLMM is applicable to a wide range of modelling tasks. An interesting future direction for the method would be targeting sub-linear time complexity by parallelisation (see, e.g., Sarkka:2019:Temporal_Parallelization_of_Bayesian_Filters).

 

Supplementary Material:

Scalable Exact Inference in Multi-Output Gaussian Processes

 

[ xscale=0.85, yscale=0.85, latent node/.append style= align=center, inner sep=0pt, text width=0.75cm ] [latent node] (xt) at (0, 0) $x_n$; [latent node] (ft) at ($(xt)+(0,1.5)$) $f_n$; [latent node] (xt2) at ($(xt)+(1,0)$) $x_{n+1}$; [latent node] (ft2) at ($(xt2)+(0,1.5)$) $f_{n+1}$; [latent node, dashed] (xtpre) at ($(xt)-(1,0)$) ; [latent node, dashed] (xtpost) at ($(xt2)+(1,0)$) ; [latent node, dashed] (ftpre) at ($(xtpre)+(0,1.5)$) ; [latent node, dashed] (ftpost) at ($(xtpost)+(0,1.5)$) ; [arrow, ->] (xt) – (ft) node [pos=.45, fill=white] $H$; [arrow, ->] (xt2) – (ft2) node [pos=.45, fill=white] $H$; [arrow, ->] (xt) – (xt2); [arrow, ->, dashed] (xtpre) – (xt); [arrow, ->, dashed] (xt2) – (xtpost); [anchor=south, rotate=90, draw=black, thick] at (-1.5, 0.9) instantaneous mixing; [anchor=south, draw=black, thick] at (0.5, 2.3) time invariant; [anchor=south] (desc) at (0.5, 1.85) $H\in {ℝ}^{p\times m}$ is a fixed matrix; [plate, fit=(xtpre) (desc) (ftpost), inner sep=5pt, inner ysep=3pt, yshift=-2pt, label=[anchor=south west]north west:LMM] () ; [ xscale=0.85, yscale=0.85, latent node/.append style= align=center, inner sep=0pt, text width=0.75cm ] [latent node] (xt) at (0, 0) $x_n$; [latent node] (ft) at ($(xt)+(0,1.5)$) $f_n$; [latent node] (xt2) at ($(xt)+(1,0)$) $x_{n+1}$; [latent node] (ft2) at ($(xt2)+(0,1.5)$) $f_{n+1}$; [latent node, dashed] (xtpre) at ($(xt)-(1,0)$) ; [latent node, dashed] (xtpost) at ($(xt2)+(1,0)$) ; [latent node, dashed] (ftpre) at ($(xtpre)+(0,1.5)$) ; [latent node, dashed] (ftpost) at ($(xtpost)+(0,1.5)$) ; [arrow, ->] (xt) – (ft) node [pos=.45, fill=white, fill opacity=.85] $H_n$; [arrow, ->] (xt2) – (ft2) node [pos=.45, fill=white, fill opacity=.85] $H_{n+1}$; [arrow, ->] (xt) – (xt2); [arrow, ->, dashed] (xtpre) – (xt); [arrow, ->, dashed] (xt2) – (xtpost); [anchor=south, draw=black, thick] at (0.5, 2.3) time-varying; [anchor=south] at (0.5, 1.85) $H:𝒯\to {ℝ}^{p\times m}$ is a time-varying matrix; [plate, fit=(xtpre) (desc) (ftpost), inner sep=5pt, inner ysep=3pt, yshift=-2pt, draw=none] () ;
[ xscale=0.85, yscale=0.85, latent node/.append style= align=center, inner sep=0pt, text width=0.75cm ] [latent node] (xt) at (0, 0) $x_n$; [latent node] (ft) at ($(xt)+(0,1.5)$) $f_n$; [latent node] (xt2) at ($(xt)+(1,0)$) $x_{n+1}$; [latent node] (ft2) at ($(xt2)+(0,1.5)$) $f_{n+1}$; [latent node, dashed] (xtpre) at ($(xt)-(1,0)$) ; [latent node, dashed] (xtpost) at ($(xt2)+(1,0)$) ; [latent node, dashed] (ftpre) at ($(xtpre)+(0,1.5)$) ; [latent node, dashed] (ftpost) at ($(xtpost)+(0,1.5)$) ; [arrow, ->, dashed] (xtpre) – (xt); [arrow, ->, dashed] (xt2) – (xtpost); [arrow, ->, dashed] (xt) – (ftpre); [arrow, ->, dashed] (xt) – (ftpost); [arrow, ->, dashed] (xt2) – (ftpre); [arrow, ->, dashed] (xt2) – (ftpost); [arrow, ->] (xt) – (ft) node [pos=0.45, fill=white, fill opacity=.85, inner sep=1pt] $H_0$; [arrow, ->] (xt) – (ft2) node [pos=0.2, fill=white, fill opacity=.85, inner sep=1pt] $H_{-1}$; [arrow, ->] (xt2) – (ft) node [pos=0.7, fill=white, fill opacity=.85, inner sep=1pt] $H_1$; [arrow, ->] (xt2) – (ft2) node [pos=0.45, fill=white, fill opacity=.85, inner sep=1pt] $H_0$; [arrow, ->] (xt) – (xt2); [anchor=south, rotate=90, draw=black, thick] at (-1.5, 0.9) convolutional mixing; [anchor=south] at (0.5, 1.85) $H:𝒯\to {ℝ}^{p\times m}$ is a time-invariant filter; [plate, fit=(xtpre) (desc) (ftpost), inner sep=5pt, inner ysep=3pt, yshift=-2pt, draw=none] () ; [ xscale=0.85, yscale=0.85, latent node/.append style= align=center, inner sep=0pt, text width=0.75cm ] [latent node] (xt) at (0, 0) $x_n$; [latent node] (ft) at ($(xt)+(0,1.5)$) $f_n$; [latent node] (xt2) at ($(xt)+(1,0)$) $x_{n+1}$; [latent node] (ft2) at ($(xt2)+(0,1.5)$) $f_{n+1}$; [latent node, dashed] (xtpre) at ($(xt)-(1,0)$) ; [latent node, dashed] (xtpost) at ($(xt2)+(1,0)$) ; [latent node, dashed] (ftpre) at ($(xtpre)+(0,1.5)$) ; [latent node, dashed] (ftpost) at ($(xtpost)+(0,1.5)$) ; [arrow, ->, dashed] (xtpre) – (xt); [arrow, ->, dashed] (xt2) – (xtpost); [arrow, ->, dashed] (xt) – (ftpre); [arrow, ->, dashed] (xt) – (ftpost); [arrow, ->, dashed] (xt2) – (ftpre); [arrow, ->, dashed] (xt2) – (ftpost); [arrow, ->] (xt) – (ft) node [pos=0.45, fill=white, fill opacity=.85, inner sep=1pt] $H_n^n$; [arrow, ->] (xt) – (ft2) node [pos=0.2, fill=white, fill opacity=.85, inner sep=1pt] $H_n^{n+1}$; [arrow, ->] (xt2) – (ft) node [pos=0.7, fill=white, fill opacity=.85, inner sep=1pt] $H_{n+1}^n$; [arrow, ->] (xt2) – (ft2) node [pos=0.45, fill=white, fill opacity=.85, inner sep=1pt] $H_{n+1}^{n+1}$; [arrow, ->] (xt) – (xt2); [anchor=south] at (0.5, 1.85) $H:𝒯\times 𝒯\to {ℝ}^{p\times m}$ is a time-varying filter; [plate, fit=(xtpre) (desc) (ftpost), inner sep=5pt, inner ysep=3pt, yshift=-2pt, draw=none] () ;

Figure 6: Graphical models illustrating the difference between time-invariant/time-varying and instantaneous/convolutional multi-output GP models, for data sampled at real-valued times $t_1,t_2,\mathrm{\dots }$ (sampling period $\mathrm{\Delta }t$). Abbreviations used: $x_n=x(t_n)$, $f_n=f(t_n)$, $H_n=H(n\mathrm{\Delta }t)$, and $H_n^m=H(t_m,t_n)$. For simplicity, the dynamics of $x$ are depicted as a Markov chain; since $x$ is modelled with a GP, $x_n$ actually depends on $x_{n^{\prime }}$ for all $n^{\prime }\le n$.

Appendix A Characterisation of MOGPs

creftypecap 6 illustrates the graphical models of time-invariant/time-varying and instantaneous/convolutional multi-output GP models.

Appendix B Maximum Likelihood Estimate

Throughout the appendix, we assume that the columns of $H$ are linearly independent and $\mathrm{\Lambda }>0$. As a consequence, $H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H$ is invertible and $H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H>0$.

Prop. 3.

Denote $p(y|x)=𝒩(y;Hx,\mathrm{\Lambda })$, and let $T$ be the $m\times p$ matrix ${(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}$. Then

$$Ty=\underset{x}{\arg \max }p(y|x)$$

and $Ty$ is an unbiased estimate of $x$: $𝔼[Ty|x]=x$.

Proof.

Note that

$$\log p(y|x)\simeq -\frac{1}{2}{(y-Hx)}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}(y-Hx)$$

Using invertibility of $H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H$, an elementary calculation then shows that the maximum with respect to $x$ is given by

$$x={(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}y=Ty.$$

To show that $Ty$ is an unbiased estimate of $x$, we use that $𝔼[y|x]=Hx$:

$$𝔼[Ty|x]=THx={(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)x=x.$$

$0.75\mathrm{\square }$

Appendix C Proof of Prop. 1

To prove creftype 1, we need the property of $T$ that it “preserves the signal-to-noise ratio”. This is characterised in the following lemma.

Lem. 1.

$$\frac{𝒩(y;Hx,\mathrm{\Lambda })}{𝒩(y;0,\mathrm{\Lambda })}=\frac{𝒩(Ty;x,{(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1})}{𝒩(Ty;0,{(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1})}.$$

Proof.

It is simple to check the equality by direct verification. We show, however, how the equality may be derived. To begin with, we have

$${(y-Hx)}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}(y-Hx)=y^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}y-2x^{\text{𝖳}}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}y+x^{\text{𝖳}}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}Hx.$$

Here

$$H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}y=(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H){(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}y=(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)Ty,$$

so

$${(y-Hx)}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}(y-Hx)=y^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}y-2x^{\text{𝖳}}(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)Ty+x^{\text{𝖳}}(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)x.$$

Adding and subtracting $y{T}^{\text{𝖳}}(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)Ty$, we find

$${(y-Hx)}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}(y-Hx)=y^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}y-y{T}^{\text{𝖳}}(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)Ty+{(x-Ty)}^{\text{𝖳}}(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)(x-Ty).$$

Hence, rearranging,

$${(y-Hx)}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}(y-Hx)-y^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}y={(x-Ty)}^{\text{𝖳}}(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)(x-Ty)-y{T}^{\text{𝖳}}(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)Ty,$$

which yields the result. $0.75\mathrm{\square }$

Proof (creftype 1).

By a general characterisation of minimal sufficient statistics (see, e.g., Th. 6.2.13 in Casella:2001:Statistical_Inference), $Ty$ is a minimal sufficient statistic for $x$ if and only if it is true that $p(y_1|x)/p(y_2|x)$ is constant as a function of $x$ if and only if $T{y}_1=T{y}_2$. Indeed, by creftype 1,

$$\log \frac{p(y_1|x)}{p(y_2|x)}={(T{y}_1-T{y}_2)}^{\text{𝖳}}{(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}x+\text{const.}$$

which, by invertibility of $H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H$, does not depend on $x$ if and only if $T{y}_1=T{y}_2$. $0.75\mathrm{\square }$

Appendix D Characterisation of Diagonal Projected Noise

Lem. 2.

Let $U$ be a truncated orthogonal matrix and $D>0$ diagonal. Then every solution $M>0$ of $U^{\text{𝖳}}{M}^{-1}U=D^{-1}$ is of the form $M=UD{U}^{\text{𝖳}}+VE{V}^{\text{𝖳}}$ with $E>0$ diagonal and $V$ such that $\left[UV\right]$ is orthogonal.

Proof.

We equivalently show that $M^{-1}$ for every solution $M>0$ of $U^{\text{𝖳}}MU=D^{-1}$ is of the desired form. Denote

$$𝒟=\{M\ge 0:U^{\text{𝖳}}MU=D^{-1}\},𝒩=\{M\ge 0:U^{\text{𝖳}}MU=0\}.$$

If $M_1,M_2\in 𝒟$, then $U^{\text{𝖳}}(M_1-M_2)U=0$ shows that $M_1-M_2\in 𝒩$. Therefore, if $M\in 𝒟$, then $𝒟=M+𝒩$. It is clear that $U{D}^{-1}{U}^{\text{𝖳}}\in 𝒟$. Thus, $𝒟=U{D}^{-1}{U}^{\text{𝖳}}+𝒩.$ Let $M\in 𝒩$, and let $M=VE{V}^{\text{𝖳}}$ be the spectral decomposition of $M$. Assume that $E>0$ by dropping the appropriate columns from $V$. Then $U^{\text{𝖳}}MU=0$ says that $U$ and $V$ have pairwise orthogonal columns, so the columns of $V$ are in $U^{⟂}$. We conclude that every $M\in 𝒟$ is of the form $M=U{D}^{-1}{U}^{\text{𝖳}}+VE{V}^{\text{𝖳}}$ with $E>0$ diagonal and $V$ a truncated orthogonal matrix with columns in $U^{⟂}$. If $M>0$, then in particular $M$ has full rank, so it is clear that $\left[UV\right]$ must form an orthogonal matrix. Therefore, $M^{-1}=UD{U}^{\text{𝖳}}+V{E}^{-1}{V}^{\text{𝖳}}$, which is of the desired form. $0.75\mathrm{\square }$

Prop. 4.

The projected noise is diagonal if and only if $H$ is of the form $H={\mathrm{\Lambda }}^{\frac{1}{2}}U{S}^{\frac{1}{2}}$ with $U$ a truncated orthogonal matrix and $S>0$ diagonal. Suppose that this is the case, and fix such a $U$. Then $H$ is of the form $H=U{D}^{\frac{1}{2}}$ with $D>0$ diagonal if and only if every column of $U$ is an eigenvector of $\mathrm{\Lambda }$.

Proof.

The projected noise is diagonal if and only if $H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H=S$ for some $S>0$ diagonal. This condition is equivalent to

$$S^{-\frac{1}{2}}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-\frac{1}{2}}{\mathrm{\Lambda }}^{-\frac{1}{2}}H{S}^{-\frac{1}{2}}=I,$$

which, in turn, holds if and only if ${\mathrm{\Lambda }}^{-\frac{1}{2}}H{S}^{-\frac{1}{2}}=U$ is a truncated orthogonal matrix. Thus, the projected noise is diagonal if and only if $H$ is of the form $H={\mathrm{\Lambda }}^{\frac{1}{2}}U{S}^{\frac{1}{2}}$ with $U$ a truncated orthogonal matrix and $S>0$ diagonal.

For the second statement, note that the condition that every column of $U$ is an eigenvector of $\mathrm{\Lambda }$ is equivalent to the condition that $\mathrm{\Lambda }$ is of the form $\mathrm{\Lambda }=UD{U}^{\text{𝖳}}+VE{V}^{\text{𝖳}}$ with $D,E>0$ diagonal and $V$ such that $\left[UV\right]$ forms an orthogonal matrix. Suppose that $\mathrm{\Lambda }$ is of this form. Then

$${\mathrm{\Lambda }}^{\frac{1}{2}}=U{D}^{\frac{1}{2}}{U}^{\text{𝖳}}+V{E}^{\frac{1}{2}}{V}^{\text{𝖳}}$$

implies that

$$H=U{D}^{\frac{1}{2}}{U}^{\text{𝖳}}U{S}^{\frac{1}{2}}+V{E}^{\frac{1}{2}}{V}^{\text{𝖳}}U{S}^{\frac{1}{2}}=U{D}^{\frac{1}{2}}{S}^{\frac{1}{2}},$$

which is of the desired form. Conversely, suppose that $H=U{D}^{\frac{1}{2}}$ with $D>0$ diagonal. Then

$$S^{-\frac{1}{2}}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H{S}^{-\frac{1}{2}}=I=S^{-\frac{1}{2}}{D}^{\frac{1}{2}}{U}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}U{D}^{\frac{1}{2}}{S}^{-\frac{1}{2}}⟹U^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}U=S{D}^{-1}.$$

Therefore, by creftype 2,

$$\mathrm{\Lambda }=U{D}^{-1}S{U}^{\text{𝖳}}+VE{V}^{\text{𝖳}}$$

with $E>0$ diagonal and $V$ such that $\left[UV\right]$ is orthogonal, which again is of the desired form. $0.75\mathrm{\square }$

Appendix E Form of Projection and Projected Noise

Prop. 5.

Let $H$ and $\mathrm{\Lambda }$ be of the form for the OLMM (creftype 2). Then

$$T=S^{-\frac{1}{2}}{U}^{\text{𝖳}}\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{\mathrm{\Lambda }}_T={\sigma }^2{S}^{-1}+D.$$

Proof.

To begin with, using the matrix inversion lemma once and once in reverse,

$H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H$ $=S^{\frac{1}{2}}{U}^{\text{𝖳}}{({\sigma }^{2}I+U{S}^{\frac{1}{2}}D{S}^{\frac{1}{2}}{U}^{\text{𝖳}})}^{-1}U{S}^{\frac{1}{2}}$

$=S^{\frac{1}{2}}{U}^{\text{𝖳}}({\sigma }^{-2}I-{\sigma }^{-4}U{S}^{\frac{1}{2}}{(D^{-1}+S^{\frac{1}{2}}{\sigma }^{-2}{U}^{\text{𝖳}}U{S}^{\frac{1}{2}})}^{-1}{S}^{\frac{1}{2}}{U}^{\text{𝖳}})U{S}^{\frac{1}{2}}$

$={\sigma }^{-2}S-{\sigma }^{-4}S{(D^{-1}+{\sigma }^{-2}S)}^{-1}S$

$={({\sigma }^2{S}^{-1}+D)}^{-1},$

so

$${\mathrm{\Lambda }}_T=T\mathrm{\Lambda }{T}^{\text{𝖳}}={(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}={\sigma }^2{S}^{-1}+D.$$

Furthermore, by a similar calculation,

$T$ $={(H^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}H)}^{-1}{H}^{\text{𝖳}}{\mathrm{\Lambda }}^{-1}$

$=({\sigma }^2{S}^{-1}+D)S^{\frac{1}{2}}{U}^{\text{𝖳}}({\sigma }^{-2}I-{\sigma }^{-4}U{S}^{\frac{1}{2}}{(D^{-1}+{\sigma }^{-2}S)}^{-1}{S}^{\frac{1}{2}}{U}^{\text{𝖳}})$

$=({\sigma }^2{S}^{-1}+D)({\sigma }^{-2}S-{\sigma }^{-4}S{(D^{-1}+{\sigma }^{-2}S)}^{-1}S)S^{-\frac{1}{2}}{U}^{\text{𝖳}}$

$=({\sigma }^2{S}^{-1}+D){({\sigma }^2{S}^{-1}+D)}^{-1}{S}^{-\frac{1}{2}}{U}^{\text{𝖳}}$

$=S^{-\frac{1}{2}}{U}^{\text{𝖳}}.$

$0.75\mathrm{\square }$

Appendix F Decomposition of the Mean Squared Error

Prop. 6.

Let $H=U{S}^{\frac{1}{2}}$ with $U$ a truncated orthogonal matrix and $S^{\frac{1}{2}}>0$ diagonal. Then

$$\underset{\text{MSE}}{\underset{﹇}{{\parallel y-Hx\parallel }^2}}=\underset{\begin{array}{c}\text{variance not}\\ \text{captured by basis}\end{array}}{\underset{﹇}{{\parallel P_{U^{⟂}}y\parallel }^2}}+\sum _{i=1}^m\stackrel{\begin{array}{c}\text{variance of}\\ i\text{<sup class="ltx\_sup"><span class="ltx\_text" style="font-size:70\%;">th</span></sup>}\text{ latent process}\end{array}}{\stackrel{﹈}{S_{ii}}}\underset{\begin{array}{c}\text{MSE of}\\ i\text{<sup class="ltx\_sup"><span class="ltx\_text" style="font-size:70\%;">th</span></sup>}\text{ latent process}\end{array}}{\underset{﹇}{{({(Ty)}_i-x_i)}^2}}$$

where $T=S^{-\frac{1}{2}}{U}^{\text{𝖳}}$ and $P_{U^{⟂}}$ is the orthogonal projection onto the orthogonal complement of $U$.

Proof.

By expanding and using orthogonality of $U$,

${\parallel y-Hx\parallel }^2$ $={\parallel y\parallel }^2-2⟨y,U{S}^{\frac{1}{2}}x⟩+{\parallel U{S}^{\frac{1}{2}}x\parallel }^2$

$={\parallel y\parallel }^2-{\parallel U^{\text{𝖳}}y\parallel }^2+{\parallel U^{\text{𝖳}}y\parallel }^2-2⟨U^{\text{𝖳}}y,S^{\frac{1}{2}}x⟩+{\parallel S^{\frac{1}{2}}x\parallel }^2$

$=⟨y,(I-U{U}^{\text{𝖳}})y⟩+{\displaystyle \sum _{i=1}^m}({⟨u_i,y⟩}^2-2⟨u_i,y⟩S_{ii}^{\frac{1}{2}}{x}_i+{(S_{ii}^{\frac{1}{2}}{x}_i)}^2)$

$=⟨y,(I-U{U}^{\text{𝖳}})y⟩+{\displaystyle \sum _{i=1}^m}{S}_{ii}(S_{ii}^{-1}{⟨u_i,y⟩}^2-2S_{ii}^{-\frac{1}{2}}⟨u_i,y⟩x_i+x_i^2)$

$=⟨y,(I-U{U}^{\text{𝖳}})y⟩+{\displaystyle \sum _{i=1}^m}{S}_{ii}{({(Ty)}_i-x_i)}^2,$

where $u_i$ is the $i$^th column of $U$. Note that $P_U=U{U}^{\text{𝖳}}$ is the orthogonal projection onto $U$, so $I-U{U}^{\text{𝖳}}=P_{U^{⟂}}$ is the orthogonal projection onto the orthogonal complement of $U$. Therefore,

$⟨y,(I-U{U}^{\text{𝖳}})y⟩=⟨y,P_{U^{⟂}}y⟩=⟨y,P_{U^{⟂}}^{2}y⟩=⟨P_{U^{⟂}}^{\text{𝖳}}y,P_{U^{⟂}}y⟩=⟨P_{U^{⟂}}y,P_{U^{⟂}}y⟩={\parallel P_{U^{⟂}}y\parallel }^2.$

$0.75\mathrm{\square }$

Appendix G Kullback–Leibler Divergence Between an LMM and OLMM

Prop. 7.

Consider two LMMs with equal $\mathrm{\Lambda }={\sigma }^2{I}_p$, equal $K(t,t^{\prime })$, but different bases $H$ and $\widehat{H}$. Let $t_1,\mathrm{\dots },t_n\in 𝒯$ and denote $x_i=x(t_i)$ and $y_i=y(t_i)$. It then holds that

$${\mathrm{D}}_{\text{KL}}(p(y_{1:n},x_{1:n})\parallel \widehat{p}(y_{1:n},x_{1:n}))={\mathrm{D}}_{\text{KL}}(\widehat{p}(y_{1:n},x_{1:n})\parallel p(y_{1:n},x_{1:n}))=n\frac{1}{2{\sigma }^2}{\parallel H-\widehat{H}\parallel }_F^2$$

and

$$\underset{\widehat{H}:\text{OLMM}}{inf}{\mathrm{D}}_{\text{KL}}(p(y_{1:n},x_{1:n})\parallel \widehat{p}(y_{1:n},x_{1:n}))\le n\frac{𝔼[{\parallel f(t)\parallel }^2]}{{\sigma }^2}\underset{i}{max}(1-V_{ii})\le n\frac{𝔼[{\parallel f(t)\parallel }^2]}{2{\sigma }^2}{\parallel I-V\parallel }_F^2$$

where $\widehat{H}$ ranges over matrices of the form $U{S}^{\frac{1}{2}}$ with $U$ a truncated orthogonal matrix and $S^{\frac{1}{2}}>0$ diagonal, $V$ is the orthogonal matrix collecting the right singular vectors of $H$, and $𝔼[{\parallel f(t)\parallel }^2]$ denotes the variance of the observations under the first LMM before adding noise.

Proof.

Start out by expanding the Kullback–Leibler divergence and noting that $p(x_{1:n})=\widehat{p}(x_{1:n})$:

${\mathrm{D}}_{\text{KL}}(p(y_{1:n},x_{1:n})\parallel \widehat{p}(y_{1:n},x_{1:n}))$ $=-{𝔼}_{p(y_{1:n},x_{1:n})}\log {\displaystyle \frac{\widehat{p}(y_{1:n}|x_{1:n})\cancel{\widehat{p}(x_{1:n})}}{p(y_{1:n}|x_{1:n})\cancel{p(x_{1:n})}}}$

$=-{\displaystyle \sum _{i=1}^n}{𝔼}_{p(y_i,x_i)}[\log \widehat{p}(y_i|x_i)-\log p(y_i|x_i)]$

$=-{\displaystyle \sum _{i=1}^n}{𝔼}_{p(y_i,x_i)}[\log 𝒩(y_i;\widehat{H}{x}_i,{\sigma }^2{I}_p)-\log 𝒩(y_i;H{x}_i,{\sigma }^2{I}_p)].$

Here

$${𝔼}_{p(y_i,x_i)}[\log 𝒩(y_i;\widehat{H}{x}_i,{\sigma }^2{I}_p)]=-\frac{p}{2}\log 2\pi {\sigma }^2-\frac{1}{2{\sigma }^2}{𝔼}_{p(y_i,x_i)}[{\parallel y_i-\widehat{H}{x}_i\parallel }^2]$$

where

${𝔼}_{p(y_i,x_i)}[{\parallel y_i-\widehat{H}{x}_i\parallel }^2]$ $={𝔼}_{p(y_i,x_i)}\mathrm{tr}[y_i{y}_i^{\text{𝖳}}-2y_i{x}_i^{\text{𝖳}}{\widehat{H}}^{\text{𝖳}}+x_i{x}_i^{\text{𝖳}}\widehat{H}{\widehat{H}}^{\text{𝖳}}]$

$=\mathrm{tr}[H{H}^{\text{𝖳}}+{\sigma }^{2}I-2H{\widehat{H}}^{\text{𝖳}}+\widehat{H}{\widehat{H}}^{\text{𝖳}}]$

$=p{\sigma }^2+\mathrm{tr}[(H-\widehat{H}){(H-\widehat{H})}^{\text{𝖳}}]$

$=p{\sigma }^2+{\parallel H-\widehat{H}\parallel }_F^2.$

Therefore,

$${\mathrm{D}}_{\text{KL}}(p(y_{1:n},x_{1:n})\parallel \widehat{p}(y_{1:n},x_{1:n}))=n\frac{1}{2{\sigma }^2}{\parallel H-\widehat{H}\parallel }_F^2.$$

Let $H=US{V}^{\text{𝖳}}$ be the SVD of $H$ where $U$ is a truncated orthogonal matrix with the same shape as $H$, $S>0$ is a square diagonal matrix, and $V$ is a square orthogonal matrix. Note that $U^{\text{𝖳}}U=I$, but $U{U}^{\text{𝖳}}\ne I$. Then, choosing $\widehat{H}=US$,

$$\underset{\widehat{H}:\text{OLMM}}{inf}{\mathrm{D}}_{\text{KL}}(p(y_{1:n},x_{1:n})\parallel \widehat{p}(y_{1:n},x_{1:n}))\le n\frac{1}{2{\sigma }^2}{\parallel U(S{V}^{\text{𝖳}}-S)\parallel }_F^2=n\frac{1}{2{\sigma }^2}{\parallel S{V}^{\text{𝖳}}-S\parallel }_F^2$$

since ${\parallel UA\parallel }_F^2=\mathrm{tr}[A^{\text{𝖳}}{U}^{\text{𝖳}}UA]=\mathrm{tr}[A^{\text{𝖳}}A]={\parallel A\parallel }_F^2$. We now further simplify:

$${\parallel S{V}^{\text{𝖳}}-S\parallel }_F^2=\mathrm{tr}[(S{V}^{\text{𝖳}}-S){(S{V}^{\text{𝖳}}-S)}^{\text{𝖳}}]=\mathrm{tr}[S{V}^{\text{𝖳}}VS-S{V}^{\text{𝖳}}S-SVS+SS]=2\mathrm{tr}[SS-SVS].$$

Hence, by definition of the trace and the fact that $S$ is diagonal,

$${\parallel S{V}^{\text{𝖳}}-S\parallel }_F^2=2\sum _{i=1}^m{S}_{ii}^2(1-V_{ii})\le 2\left(\sum _{i=1}^m{S}_{ii}^2\right)\underset{i}{max}(1-V_{ii})=2𝔼[{\parallel f\parallel }^2]\underset{i}{max}(1-V_{ii}),$$

since

$$𝔼[{\parallel f(t)\parallel }^2]=𝔼\mathrm{tr}[f(t)f^{\text{𝖳}}(t)]=\mathrm{tr}[H{H}^{\text{𝖳}}]=\mathrm{tr}[S^2].$$

Therefore,

$${\parallel S{V}^{\text{𝖳}}-S\parallel }_F^2\le 2𝔼[{\parallel f\parallel }^2]\underset{i}{max}(1-V_{ii})\le 2𝔼[{\parallel f\parallel }^2]\sum _{i=1}^m(1-V_{ii})=𝔼[{\parallel f\parallel }^2]{\parallel I-V\parallel }_F^2,$$

where the equality follows from a similar calculation:

$${\parallel I-V\parallel }_F^2=\mathrm{tr}[I-V^{\text{𝖳}}-V+V^{\text{𝖳}}V]=2\mathrm{tr}[I-V].$$

$0.75\mathrm{\square }$

Appendix H Missing Data

Missing data is troublesome for the OLMM, because it is not possible to take away a subset of the rows of $H$ and retain orthogonality of the columns. Suppose that data are missing for a fixed subset of outputs. Collect that subset of outputs into $y_1$, and let $y_2$ be the remainder of the outputs. The issue can then be resolved by considering a multi-headed version of the OLMM:

$y_1|H_1,x$ $\sim 𝒢𝒫(H_{1}x(t),\delta [t-t^{\prime }]{\mathrm{\Lambda }}_1),$

$y_2|H_2,x$ $\sim 𝒢𝒫(H_{2}x(t),\delta [t-t^{\prime }]{\mathrm{\Lambda }}_2).$

For $k$ heads, if $y_I={(y_i)}_{i\in I}$ for a subset of heads $I\subseteq \{1,\mathrm{\dots },k\}$ are observed, then

$$T_I{y}_I={\mathrm{\Lambda }}_{T_I}{\sum }_{i\in I}{\mathrm{\Lambda }}_{T_i}^{-1}{T}_i{y}_i,{\mathrm{\Lambda }}_{T_I}^{-1}={\sum }_{i\in I}{\mathrm{\Lambda }}_{T_i}^{-1}.$$

While this new model resolves the missing data problem, the assumption that each collection of outputs (each head) is a linear transformation of the same collection of latent processes will be overly restrictive for many contexts: it is not typically the case that our prior beliefs about the relationships between outputs have anything to do with the availability of data pertaining to them.

In the general case, where data can be missing arbitrarily, let $Y_{\text{o}}$ be the observed data. Complement $Y_{\text{o}}$ with missing data $Y_{\text{m}}$ such that $Y=Y_{\text{o}}\cup Y_{\text{m}}$ is complete. Then a simple solution is to use variational inference. Assume a Gaussian approximate posterior distribution $q(Y_{\text{m}})$ over $Y_{\text{m}}$, and maximise the evidence lower bound (ELBO) $ℒ$ using gradient-based optimisation:

$$\log p(Y_{\text{o}})\ge {𝔼}_{q(Y_{\text{m}})}[\log p(Y)]+H[q(Y_{\text{m}})]=ℒ[q(Y_{\text{m}})],$$

where the expectation can be approximated using the reparametrisation trick (Kingma:2013:Auto-Encoding_VB), $\log p(Y)$ can be computed efficiently because $Y$ is complete, and $H[q(Y_{\text{m}})]$ denotes the entropy of $q(Y_{\text{m}})$. This approach provides a tractable solution when the number of missing data are not too numerous.

Appendix I Heterogeneous Observation Noise

Although the specification of the observation noise $\mathrm{\Lambda }={\sigma }^2{I}_p+HD{H}^{\text{𝖳}}$ in the OLMM does not allow for heterogeneous observation noise, it is possible to set $\mathrm{\Lambda }=\mathrm{diag}({\sigma }_1^2,\mathrm{\dots },{\sigma }_p^2)$ and use creftype 4 to include $\mathrm{\Lambda }$ in the parametrisation of $H$: $H={\mathrm{\Lambda }}^{\frac{1}{2}}U{S}^{\frac{1}{2}}$. This parametrisation can be interpreted in two ways:

The model has a whitening transform built in. In the projection $T$, the (noise in the) data will first by whitened by ${\mathrm{\Lambda }}^{-\frac{1}{2}}$. Hence, this parametrisation can be used as a more principled substitute for the usual data normalisation where the outputs are divided by their empirical standard deviation prior to feeding them to the model.

The basis is orthogonal with respect to a weighted Euclidean inner product: ${⟨h_i,h_j⟩}_{{\mathrm{\Lambda }}^{-1}}={\sum }_{k=1}^p{h}_{ik}{h}_{jk}/{\sigma }_k^2=0$ for $i\ne j$. Intuitively, this means that the basis is orthogonal in the usual sense after stretching the $i$^th dimension by ${\sigma }_i^{-1}$. Although this construction provides additional flexiblity, it does require that $D=0$ to avoid a circular dependency between $\mathrm{\Lambda }$ and $H$.

Appendix J Point Process Experimental Details and Additional Results

We use a total of $np=200\times 100$ bins in the model, where each bin is $5\mathrm{m}$ $\times$ $5\mathrm{m}$. The likelihood for this experiment is given by

$$p(Y\mid f)\approx \prod _{i=1}^n\prod _{j=1}^p\mathrm{Poisson}(Y_{ij}\mid \exp (f(r_{ij})))$$

where $r_{ij}$ is the coordinate of the $ij$^th bin, $Y_{ij}$ is the number of data points in the $ij$^th bin, and $Y$ is the $n\times p$ matrix of counts.

[axis on top, xmin=0, xmax=1000, xlabel=$x_1$ (meters), y dir=reverse, ymin=0, ymax=500, ylabel=$x_2$ (meters), axis background/.style=fill=white, legend style=legend cell align=left,align=left,draw=white!15!black, width=height=] [forget plot] graphics [xmin=-0.394321766561514,xmax=1000.39432176656,ymin=-0.394944707740916,ymax=500.394944707741] ./figures/trees-1.png;

[point meta min=-3, point meta max=3, axis on top, xmin=0, xmax=1000, xlabel=$x_1$ (meters), y dir=reverse, ymin=0, ymax=500, ylabel=$x_2$ (meters), axis background/.style=fill=white, legend style=legend cell align=left,align=left,draw=white!15!black, width=height=] [forget plot] graphics [xmin=-5.05050505050505,xmax=1005.05050505051,ymin=-1.25628140703518,ymax=501.256281407035] ./figures/intensity-1.png;

Figure 7: Observations of the rainforest tree locations (left), and posterior mean log-intensity for the log-Gaussian Cox process model (right) with a grid of $np=20000$ observation bins.

We perform ${10}^5$ iterations of Gibbs sampling, alternating between Elliptical Slice Sampling for the Gaussian process given its hyperparameters, and Metropolis Hastings with proposal distribution $𝒩(\theta ,{0.05}^2)$ for the log of the hyperparameters given the latent GP-distributed function. The kernel is a product of two Matérn-$5/2$ kernels, one over the vertical dimension and the other over the horizontal one, each with their own length scale. A single process variance is utilised, and a nugget term is added. The $\log$ of each of the four hyperparameters was given a $𝒩(0,1)$ prior. creftypecap 8 shows the log joint of the entire state after each iteration, while creftype 9 shows the progress of each hyperparameter per iteration.

Figure 8: Log–joint probability per iteration

Figure 9: Hyperparameters per iteration. Shows the horizontal length scale (top-left), vertical length scale (top-right), nugget variance (bottom-right), and process variance (bottom-left).

The times in creftype 4 were obtained via BenchmarkTools.jl (BenchmarkTools.jl-2016). The implementation of the standard Kronecker product decomposition trick makes use of Kronecker.jl, and Julia’s (bezanson2017julia) standard linear algebra libraries, which make use of OpenBLAS and LAPACK to efficiently perform matrix-matrix products and compute eigendecompositions. The implementation of the state-space GP additionally makes use of StaticArrays.jl for efficient stack-allocated matrices, and Stheno.jl for GP-related functionality. Timing experiments were conducted on a single CPU core to make the implementations comparable.

The state-space implementation of the GP makes use of the trick introduced by Solin:2018:Infinite-Horizon_Gaussian_Processes when computing the log marginal likelihood. However, the infinite-horizon trick is exploited only once the filtering covariance has converged, which is determined by the point at which the Frobenius norm of the difference between the filtering covariance at the $t$^th and $(t-1)$^th iterations drops below ${10}^{-12}$. In practice, we find that this produces log-marginal likelihoods that are essentially exact, as expected.