Fast Multitask Gaussian Process Regression

Gaussian process (GP) regression is a powerful probabilistic modeling technique with built-in uncertainty quantification. When one has access to multiple correlated simulations (tasks), it is common to fit a multitask GP (MTGP) surrogate which is capable of capturing both inter-task and intra-task correlations. However, with a total of $N$ evaluations across all tasks, fitting an MTGP is often infeasible due to the $\mathcal{O}(N^2)$ storage and $\mathcal{O}(N^3)$ computations required to store, solve a linear system in, and compute the determinant of the $N \times N$ Gram matrix of pairwise kernel evaluations. In the single-task setting, one may reduce the required storage to $\mathcal{O}(N)$ and computations to $\mathcal{O}(N \log N)$ by fitting "fast GPs" which pair low-discrepancy design points from quasi-Monte Carlo to special kernel forms which yields nicely structured Gram matrices, e.g., circulant matrices. This article generalizes fast GPs to fast MTGPs which pair low-discrepancy design points for each task to special product kernel forms which yields nicely structured block Gram matrices, e.g., circulant block matrices. An algorithm is presented to efficiently store, invert, and compute the determinant of such Gram matrices with optionally different sampling nodes and different sample sizes for each task. Derivations for fast MTGP Bayesian cubature are also provided. A GPU-compatible, open-source Python implementation is made available in the FastGPs package (https://alegresor.github.io/fastgps/). We validate the efficiency of our algorithm and implementation compared to standard techniques on a range of problems with low numbers of tasks and large sample sizes.