  A Gaussian process (GP) is a powerful and widely used regression technique.
The main building block of a GP regression is the covariance kernel, which
characterizes the relationship between pairs in the random field. The
optimization to find the optimal kernel, however, requires several large-scale
and often unstructured matrix inversions. We tackle this challenge by
introducing a hierarchical matrix approach, named HMAT, which effectively
decomposes the matrix structure, in a recursive manner, into significantly
smaller matrices where a direct approach could be used for inversion. Our
matrix partitioning uses a particular aggregation strategy for data points,
which promotes the low-rank structure of off-diagonal blocks in the
hierarchical kernel matrix. We employ a randomized linear algebra method for
matrix reduction on the low-rank off-diagonal blocks without factorizing a
large matrix. We provide analytical error and cost estimates for the inversion
of the matrix, investigate them empirically with numerical computations, and
demonstrate the application of our approach on three numerical examples
involving GP regression for engineering problems and a large-scale real
dataset. We provide the computer implementation of GP-HMAT, HMAT adapted for GP
likelihood and derivative computations, and the implementation of the last
numerical example on a real dataset. We demonstrate superior scalability of the
HMAT approach compared to built-in $\backslash$ operator in MATLAB for
large-scale linear solves $\bf{A}\bf{x} = \bf{y}$ via a repeatable and
verifiable empirical study. An extension to hierarchical semiseparable (HSS)
matrices is discussed as future research.
