  Gaussian random fields (GRF) are a fundamental stochastic model for
spatiotemporal data analysis. An essential ingredient of GRF is the covariance
function that characterizes the joint Gaussian distribution of the field.
Commonly used covariance functions give rise to fully dense and unstructured
covariance matrices, for which required calculations are notoriously expensive
to carry out for large data. In this work, we propose a construction of
covariance functions that result in matrices with a hierarchical structure.
Empowered by matrix algorithms that scale linearly with the matrix dimension,
the hierarchical structure is proved to be efficient for a variety of random
field computations, including sampling, kriging, and likelihood evaluation.
Specifically, with $n$ scattered sites, sampling and likelihood evaluation has
an $O(n)$ cost and kriging has an $O(\log n)$ cost after preprocessing,
particularly favorable for the kriging of an extremely large number of sites
(e.g., predicting on more sites than observed). We demonstrate comprehensive
numerical experiments to show the use of the constructed covariance functions
and their appealing computation time. Numerical examples on a laptop include
simulated data of size up to one million, as well as a climate data product
with over two million observations.
