  Maximum likelihood estimation for parameter-fitting given observations from a
Gaussian process in space is a computationally-demanding task that restricts
the use of such methods to moderately-sized datasets. We present a framework
for unstructured observations in two spatial dimensions that allows for
evaluation of the log-likelihood and its gradient (i.e., the score equations)
in $\tilde O(n^{3/2})$ time under certain assumptions, where $n$ is the number
of observations. Our method relies on the skeletonization procedure described
by Martinsson & Rokhlin in the form of the recursive skeletonization
factorization of Ho & Ying. Combining this with an adaptation of the matrix
peeling algorithm of Lin et al. for constructing $\mathcal{H}$-matrix
representations of black-box operators, we obtain a framework that can be used
in the context of any first-order optimization routine to quickly and
accurately compute maximum-likelihood estimates.
