Normative Modeling of Neuroimaging Data using Scalable Multi-Task Gaussian Processes

For $N$ samples and $T$ tasks, the time and space complexity of MTGPR are $\mathcal{O}(N^{3}T^{3})$ and $\mathcal{O}(N^{2}T^{2})$, respectively. These high computational demands (compared to STGPR with $\mathcal{O}(N^{3}T)$ and $\mathcal{O}(N^{2}T)$) are mainly due to the need for computing the inverse cross-covariance matrix in learning and inference phases. In neuroimaging problems that we consider, these can both be relatively high where $N$ in generally in the order of $10^{2}-10^{4}$ and $T$ is in the order of $10^{4}-10^{5}$ or even higher. Therefore, improving the computational efficiency of MTGPR is crucial for certain problems, and there have been several approaches proposed for this in the machine learning literature [13, 3]. Here we briefly review two main directions to address the computational tractability issue of MTGPR.

In the first set of approaches, approximation techniques are used to improve estimation efficiency. Bonilla et al. [5] made one of the earliest efforts in this direction, in which they proposed to use Nyström approximation on $M$ inducing inputs [13] out of $N$ samples in combination with the probabilistic principal component analysis, in order to approximate reduced $M$-rank and $P$-rank sample and task covariance matrices, respectively. Their approximation reduced the time complexity of hyperparameter learning to $\mathcal{O}(NTM^{2}P^{2})$. Elsewhere, Alvarez and Lawrence [2] proposed to approximate a sparse version of MTGPR, assuming conditional independence between each output variable with all others given the input process. This assumption besides using $M$ out of $N$ input samples as inducing inputs reduces the computational complexity of MTGPR to $\mathcal{O}(N^{3}T+NTM^{2})$ and $\mathcal{O}(N^{2}T+NTM)$ in time and storage, where for $N=M$ is the same as a set of $T$ independent STGPRs. Alvarez et al. in [4] extended their previous work by developing the concept of inducing function rather than inducing input. Their new approach so-called variational inducing kernels achieves time complexity of $\mathcal{O}(NTM^{2})$.

The second set of approaches utilize properties of Kronecker product [8] to reduce the time and space complexity in computing the exact (and not approximated) inverse covariance matrix. Stegle et al. [15] proposed to use these properties in combination with eigenvalue decomposition of input and task covariance matrices for efficient parameter estimation, and likelihood evaluation/optimization in MTGPR. In this method, the joint covariance matrix is defined as a Kronecker product between the input and task covariance matrices. This approach reduces the time and space complexity of MTGPR to $\mathcal{O}(N^{3}+T^{3})$ and $\mathcal{O}(N^{2}+T^{2})$, respectively. To account also for structured noise, Rakitsch et al. [14] extended this method by using two separate Kronecker products for the signal and noise. Importantly, this provides a significant reduction in computational complexity using all samples (i.e., not just inducing inputs), and is exact in the sense that it does not require any approximation or relaxing assumptions.

Boldface capital letters, $\mathbf{A}$, and capital letters, $A$, are used to denote matrices and scalar numbers. We denote the vertical vector which is resulted from collapsing columns of a matrix $\mathbf{A}\in\mathbb{R}^{N\times T}$ with $vec(\mathbf{A})\in\mathbb{R}^{NT}$. In the remaining text, we use $\otimes$ and $\odot$ to respectively denote Kronecker and the element-wise matrix products. We denote an identity matrix by $\mathbf{I}$; and the determinant, diagonal elements, and the trace of matrix $\mathbf{A}$ with $\left|\mathbf{A}\right|$, $diag(\mathbf{A})$, and $Tr[\mathbf{A}]$, respectively.

In our experiments, we use a public fMRI dataset collected for reconstructing visual stimuli (black and white letters and symbols) from fMRI data [12]. In this dataset, fMRI responses were measured while $10\times 10$ checkerboard patch images were presented to subjects according to a blocked design. Checkerboard patches constituted random (1320 trials) and geometrically meaningful patterns (720 trials). We use the preprocessed data available in Nilearn package [1] wherein the fMRI data are detrended and masked for the occipital lobe (5438 voxels).²² 2 See http://nilearn.github.io/auto_examples/02_decoding/plot_miyawaki_reconstruction.html. Whilst our approach is quite general, we demonstrate S-MTGPR by simulating normative modeling for novelty detection. Therefore, we aim to predict the masked fMRI 3D-volume from the presented visual stimuli in an encoding setting. To this end, we randomly selected 600 random pattern trials, for training the encoding model. The model then learns to represent this reference or normative class such that anomalous or abnormal samples can be detected and characterised. The rest of non-random patterns (720 trials) and random patterns (720 trials) are used for evaluating the encoding model and testing anomaly-detection performance, achieved by fitting a generalised extreme value distribution to the most deviating voxels. In our experiments, we use PCA to transform the fMRI data in the training set from the voxel space to $\mathbf{Z}$, and the resulting $P=10,25,50,100,250,500,1000$ PCA components are used as basis matrix $\mathbf{B}$ in the optimization and inference.

We benchmark the proposed method against the STGPR (i.e., mass-univariate) and MT-Kronprod models in terms of their runtime, performance of the regression, and quality of resulting normative models. In all models, we use a summation of a linear, a squared exponential, and a diagonal isotropic covariance functions for sample and task covariance matrices in order to accommodate both linear and non-linear relationships. In all cases, we use an isotropic Gaussian likelihood function. This likelihood function has different functionality in the STGPR versus MTGPR settings. In STGPR, it is defined independently for each voxel, thus it handles heteroscedastic, i.e., spatially varying noise. While in MTGPR a single noise parameter is shared for all voxels, hence it merely considers homoscedastic, i.e., spatially stationary, noise. The truncated Newton algorithm is used for optimizing the parameters. Table 1 summarizes the time complexity and the number of parameters of three benchmarked methods in our experiments.

We use the coefficient of determination ($R^{2}$) to evaluate the explained variance by regression models. In normative modeling, the top 5% values in normative probability maps are used to fit the generalized extreme value distribution (see [9]). To evaluate resulting normative models, we employ area under the curve (AUC) to measure the performance of the model in distinguishing between normal (here random patterns) from abnormal samples (here non-random patterns). All the steps (random sampling, modeling, and evaluation) are repeated 10 times in order to estimate the mean and standard deviation of the runtime, $R^{2}$, and AUC. All experiments are performed on a system with Intel^®Xeon^®E5-1620 0 @3.60GHz CPU and 16GB of RAM.³³ 3 The experimental codes are available at https://github.com/smkia/MTNorm.

Fig. 1 compares the runtime, $R^{2}$, and AUC of STGPR and MT-Kronprod, with those of S-MTGPR for different number of bases. As illustrated in Fig. 1(a) S-MTGPR is faster than other approaches where the total runtime of MT-Kronprod (3 days) and STGPR (6 hours) can be reduced to 16 minutes for $P=25$. This difference in runtime is even more pronounced in case of the optimization time where S-MTGPR is at least (for $P=1000$) 33 and 89 times faster than STGPR and MT-Kronprod, respectively. The multi-task approaches are slower than STGPR in the prediction phase mainly due to the mini-batch implementation of the prediction variance computation (to avoid memory overflow). Fig. 1(b) shows this computational efficiency is achieved without penalty to the regression performance; where for certain number of bases the S-MTGPR shows equivalent and even better $R^{2}$ than STGPR and MT-Kronprod. Furthermore, Fig. 1(c) demonstrates that multi-task learning, by considering spatial structures, generally provides a more accurate normative model of fMRI data in that it more accurately detects samples that were derived from a different distribution to those used to train the model. This fact is well-reflected in higher AUC values for S-MTGPR at $P=25,100,250,500,1000$. It is worthwhile to emphasize that these improvements are achieved by reducing the degree-of-freedom of the normative model from 21752 for STGPR to 10 for S-MTGPR (see Table 1).

For $\mathbf{A}\in\mathbb{R}^{M\times N}$, $\mathbf{B}\in\mathbb{R}^{P\times Q}$, and $\mathbf{C}$, $\mathbf{D}$ (with appropriate size) we have:

1. 1.

   $\mathbf{A=U_{A}S_{A}U_{A}^{\top}}$ is the eigenvalue decomposition of $\mathbf{A}$,

2. 2.

   $\mathbf{(ACB)^{-1}=B^{-1}C^{-1}A^{-1}}$,

3. 3.

   $\mathbf{(A\otimes B)(C\otimes D)=AC\otimes BD}$,

4. 4.

   $\mathbf{(A\otimes B)^{-1}=A^{-1}\otimes B^{-1}}$,

5. 5.

   the eigenvalue decomposition of $\mathbf{A\otimes B+I}$ is:
   $\mathbf{(U_{A}\otimes U_{B})(S_{A}\otimes S_{B}+I)(U_{A}^{\top}\otimes U_{B}^{\top})}$,

6. 6.

   $(\mathbf{A\otimes B})vec(\mathbf{C})=vec(\mathbf{BCA}^{\top})$,

7. 7.

   $\ln\left|\mathbf{AC}\right|=\ln(\left|\mathbf{A}\right|\left|\mathbf{C}\right|)=\ln\left|\mathbf{A}\right|+\ln\left|\mathbf{C}\right|$,

8. 8.

   for $\mathbf{C}\in\mathbb{R}^{N\times N}$, $\frac{\mathrm{d}}{\mathrm{d}x}\ln\left|\mathbf{C}\right|=Tr[\mathbf{C}^{-1}\frac{\mathrm{d}\mathbf{C}}{\mathrm{d}x}]$,

9. 9.

   $Tr[\mathbf{ACBD}]=Tr[\mathbf{CBDA}]=Tr[\mathbf{BDAC}]=Tr[\mathbf{DACB}]$.

Eq. 6(a) is derived from Eq. 5(a) as follows:

Eq. 6(b) is derived from Eq. 5(b) as follows:

Eq. 7 is derived as follows:

In the optimization process, the derivatives of $\mathcal{L}$ with respect to $\theta_{\mathbf{C}}\in\Theta_{\mathbf{C}}$, $\theta_{\mathbf{R}}\in\Theta_{\mathbf{R}}$, and $\theta_{\sigma^{2}}\in\Theta_{\sigma^{2}}$ can be efficiently computed as follows: