  A number of problems in probability and statistics can be addressed using the
multivariate normal (Gaussian) distribution. In the one-dimensional case,
computing the probability for a given mean and variance simply requires the
evaluation of the corresponding Gaussian density. In the $n$-dimensional
setting, however, it requires the inversion of an $n \times n$ covariance
matrix, $C$, as well as the evaluation of its determinant, $\det(C)$. In many
cases, such as regression using Gaussian processes, the covariance matrix is of
the form $C = \sigma^2 I + K$, where $K$ is computed using a specified
covariance kernel which depends on the data and additional parameters
(hyperparameters). The matrix $C$ is typically dense, causing standard direct
methods for inversion and determinant evaluation to require $\mathcal O(n^3)$
work. This cost is prohibitive for large-scale modeling. Here, we show that for
the most commonly used covariance functions, the matrix $C$ can be
hierarchically factored into a product of block low-rank updates of the
identity matrix, yielding an $\mathcal O (n\log^2 n) $ algorithm for inversion.
More importantly, we show that this factorization enables the evaluation of the
determinant $\det(C)$, permitting the direct calculation of probabilities in
high dimensions under fairly broad assumptions on the kernel defining $K$. Our
fast algorithm brings many problems in marginalization and the adaptation of
hyperparameters within practical reach using a single CPU core. The combination
of nearly optimal scaling in terms of problem size with high-performance
computing resources will permit the modeling of previously intractable
problems. We illustrate the performance of the scheme on standard covariance
kernels.
