# Scalable Computations for Nonstationary Gaussian Processes

Paul G. Beckman
Courant Institute of Mathematical Sciences, New York University
Mathematics and Computer Science Division, Argonne National Laboratory
Christopher J. Geoga
Mathematics and Computer Science Division, Argonne National Laboratory
Department of Statistics, Rutgers University
Michael L. Stein
Department of Statistics, Rutgers University
Department of Statistics, University of Chicago
Mihai Anitescu22footnotemark: 2
Mathematics and Computer Science Division, Argonne National Laboratory
Department of Statistics, University of Chicago
corresponding author<EMAIL_ADDRESS>This material was based upon
work supported by the US Department of Energy, Office of Science, Office of
Advanced Scientific Computing Research (ASCR) under Contract DE-AC02-06CH11347

###### Abstract

Nonstationary Gaussian process models can capture complex spatially varying
dependence structures in spatial datasets. However, the large number of
observations in modern datasets makes fitting such models computationally
intractable with conventional dense linear algebra. In addition, derivative-
free or even first-order optimization methods can be slow to converge when
estimating many spatially varying parameters. We present here a computational
framework that couples an algebraic block-diagonal plus low-rank covariance
matrix approximation with stochastic trace estimation to facilitate the
efficient use of second-order solvers for maximum likelihood estimation of
Gaussian process models with many parameters. We demonstrate the effectiveness
of these methods by simultaneously fitting $192$ parameters in the popular
nonstationary model of Paciorek and Schervish using $107{,}600$ sea surface
temperature anomaly measurements.

Keywords: Nonstationary, Spatial Analysis, Optimization, Statistical Computing

##  1 Introduction

Gaussian processes are a prevalent class of models in spatial and
spatiotemporal statistics. This is due in part to the fact that the model is
completely specified by the first two moments of the Gaussian distribution, so
that a practitioner need select only a mean function and covariance function.
Let $Z(\bm{x})$ be a Gaussian process with mean function $\mu(\bm{x})=0$ and
covariance function

$\operatorname{Cov}(Z(\bm{x}),Z(\bm{x}^{\prime}))=k(\bm{x},\bm{x}^{\prime}|\bm{\theta})$
(1)

observed at locations $\\{\bm{x}_{i}\\}_{i=1}^{n}$ corresponding to
measurements $y_{i}=Z(\bm{x}_{i})$ with
$\bm{x}_{i}\in\Omega\subset\mathbb{R}^{d}$. Here $k(\cdot,\cdot)$ is a
parametric covariance function with parameters $\bm{\theta}$. Then defining
the vector $\bm{y}=[y_{1},...,y_{n}]^{\top}\in\mathbb{R}^{n}$, we have

$\bm{y}\sim\mathcal{N}\big{(}\bm{0},\bm{\Sigma}(\bm{\theta})\big{)},$ (2)

where the covariance matrix $\bm{\Sigma}(\bm{\theta})\in\mathbb{R}^{n\times
n}$ is defined by

$\bm{\Sigma}(\bm{\theta})_{ij}=\operatorname{Cov}(y_{i},y_{j})=k(\bm{x}_{i},\bm{x}_{j}|\bm{\theta}).$
(3)

A primary objective after selecting a parametric covariance model is to
estimate $\bm{\theta}$ from data. One can then predict the value of the
process at unobserved locations by treating the estimated parameters as the
true parameters. A standard parameter estimation method is to compute the
maximum likelihood estimator (MLE), denoted $\widehat{\bm{\theta}}$, which
minimizes the mean zero Gaussian negative log-likelihood function

$-\ell(\bm{\theta})=\frac{1}{2}\operatorname{logdet}\bm{\Sigma}(\bm{\theta})+\frac{1}{2}\bm{y}^{\top}\bm{\Sigma}(\bm{\theta})^{-1}\bm{y}+\frac{n}{2}\log(2\pi).$
(4)

For the remainder of this paper we will suppress the dependence of
$\bm{\Sigma}$ on $\bm{\theta}$ for notational clarity.

Computing (4) for large spatial datasets is computationally challenging, since
the determinant and linear solve operations have cubic time complexity and
quadratic space complexity using conventional dense linear algebra. Thus, for
large $n$ evaluating the log-likelihood directly becomes prohibitively
expensive. In response to the first of these computational challenges, a
number of approximations have been proposed, including Vecchia and “nearest-
neighbor Gaussian process” methods [37, 36, 20, 14], matrix tapering [11], and
Markov random field approximations [23], which each approximate $\bm{\Sigma}$
or $\bm{\Sigma}^{-1}$ using a sparse or block-sparse matrix. Low-rank updates
to a diagonal matrix have also been considered [9, 5, 10, 18, 32], although
they suffer from limitations for nonsmooth fields or when a nugget is not an
appropriate modeling assumption [34]. Approaches using hierarchical matrices
have also been proposed [6, 2, 8, 26, 24], which use rank-structured matrix
arithmetic to efficiently evaluate (4).

In spite of the development of these scalable methods, minimizing the negative
log-likelihood remains challenging depending on the covariance function and
its parameterization. Even for covariance functions with few parameters, for
example the stationary isotropic Matérn covariance

$\mathcal{M}_{\nu}(\left\lVert\bm{x}_{i}-\bm{x}_{j}\right\rVert)=\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\bigg{(}\frac{2\sqrt{\nu}\left\lVert\bm{x}_{i}-\bm{x}_{j}\right\rVert}{\rho}\bigg{)}\mathcal{K}_{\nu}\bigg{(}\frac{2\sqrt{\nu}\left\lVert\bm{x}_{i}-\bm{x}_{j}\right\rVert}{\rho}\bigg{)},$
(5)

where $\nu$ is a positive constant smoothness parameter and
$\mathcal{K}_{\nu}(\cdot)$ is the modified Bessel function of the second kind
of order $\nu$, simultaneously estimating the scale $\sigma^{2}$ and range
$\rho$ parameters poses significant challenge [39]. The likelihood surface is
far from convex, and even first-order methods may become trapped in nearly
flat nonellipsoidal regions of the likelihood surface and fail to make
meaningful progress. While a few relatively new articles attempt to treat the
optimization problem more seriously [14, 26, 12], the norm in practice is to
use derivative-free or first-order methods with finite difference derivatives.

For more complicated covariance functions that provide more flexibility and
thus require more parameters, the optimization problem becomes even harder,
and second-order optimization can mean the difference between an optimizer
stagnating and successfully reaching the MLE. The computational problem posed
in this setting is that efficient linear solves and log-determinants with
$\bm{\Sigma}$ are no longer sufficient. Derivatives
$\frac{\partial\bm{\Sigma}}{\partial\theta_{j}}$ need to be computed and
applied scalably; and if the gradient and Fisher matrix are computed directly,
matrix-matrix products of the form
$\bm{\Sigma}^{-1}\Big{(}\frac{\partial\bm{\Sigma}}{\partial\theta_{j}}\Big{)}$
need to be computed efficiently. Even for data sizes in which $\bm{\Sigma}$
does not require some form of approximation, with sufficiently many parameters
the computational burden of the matrix-matrix operations necessary for the
gradient alone can be problematic. In this work we address the problem of
second-order optimization for a covariance function with many parameters,
using as our motivating example a nonstationary spatial model, which we
introduce now.

For large spatial datasets it is often unrealistic to assume that process
parameters are constant over the entire domain, in other words that the
process is stationary. Therefore nonstationary covariance functions in which
the parameters vary in space become necessary in order to accurately capture
the dependence structure of the data. One such covariance function that we
will use here is derived by Paciorek and Schervish [27] as a modification of
the stationary Matérn covariance (5) and is frequently used in the
nonstationary Gaussian process literature [22, 28, 5, 29, 15]. We use the
following anisotropic version,

$k(\bm{x}_{i},\bm{x}_{j})=\sigma^{2}\frac{\left|\bm{\Lambda}(\bm{x}_{i})\right|^{\frac{1}{4}}\left|\bm{\Lambda}(\bm{x}_{j})\right|^{\frac{1}{4}}}{\left|\frac{\bm{\Lambda}(\bm{x}_{i})+\bm{\Lambda}(\bm{x}_{j})}{2}\right|^{\frac{1}{2}}}\mathcal{M}_{\nu}\Bigg{(}\sqrt{(\bm{x}_{i}-\bm{x}_{j})^{\top}\bigg{(}\frac{\bm{\Lambda}(\bm{x}_{i})+\bm{\Lambda}(\bm{x}_{j})}{2}\bigg{)}^{-1}(\bm{x}_{i}-\bm{x}_{j})}\Bigg{)},$
(6)

where $\bm{\Lambda}(\cdot)$ is a spatially varying function that assigns a
positive definite local anisotropy matrix at each location. In fact, Stein [4]
gives an extension of this nonstationary covariance that allows for spatially
varying scale $\sigma(\cdot)$ and smoothness $\nu(\cdot)$ functions. However,
jointly estimating range and scale parameters can be challenging even in the
simplest stationary settings [39], and robustly computing the derivatives of
$\mathcal{M}_{\nu}$ in the smoothness parameter $\nu$ is numerically
challenging, although progress on this topic has been made recently [13]. Thus
we restrict our investigations in this work to models that have only
nonstationary local ranges via the anisotropy parameters.

To estimate the spatially varying function $\bm{\Lambda}(\cdot)$ through
maximum likelihood, one can parameterize it using a set of parameter functions
$\theta(\cdot)$ such as its eigenvalues or elements of its Cholesky factors;
expand these spatially varying functions in a basis

$\displaystyle\theta(\bm{x})$
$\displaystyle=\sum_{i=1}^{m}\phi_{i}(\bm{x})\,\theta_{i},$ (7)

where $\phi_{i}$ are some set of basis functions; and estimate the
coefficients $\theta_{i}$. We adopt this approach here, using a radial basis
function (RBF) expansion (see Section 4). Other parameterizations and bases
can be used, however, and we describe here a framework for high-dimensional
parameter optimization with large data that is essentially agnostic of these
choices.

Since large spatial regions or highly nonstationary random fields require a
large number $m$ of basis functions, fitting the entire global model becomes a
difficult high-dimensional optimization problem for which derivative-free and
even first-order solvers are often ineffective. A number of past works using
this model circumvent this difficulty by fitting parameters only locally and
subsequently smoothing or regressing them into a global model using
$(\ref{eq:basis})$ where the $\phi_{i}$ are RBFs [22, 28, 15]. Without
considering all the basis coefficients simultaneously, however, relationships
between parameters in adjacent regions are ignored, and it can be difficult to
verify the quality of the global model fit. The primary goal of our work is to
present a computational framework in which this global model fitting is
tractable and to consider its benefits.

We emphasize here that the two computational challenges of large data size and
large parameter dimension are intimately linked. Large data often require
complex models in order to accurately characterize their dependence structure,
and, conversely, complex models often require large data in order to
accurately estimate their many parameters. Therefore, considering these two
problems together is imperative to developing effective practical methods.

We pursue here work in this vein, introducing methods to control the
computational complexity of second-order global optimization. Further, we
provide and discuss an application to high-resolution sea surface temperature
measurements, fitting a dataset of 107,600 measurements to a model with 192
parameters. The results of this application demonstrate that there may be
substantial gains to model fit by estimating the entire global model jointly.

We first introduce substantial improvements to the likelihood computations
provided by [12] by showing that the hierarchical matrix approximation therein
can be expressed more simply as the block full-scale approximation [31, 30].
We then outline the computations of exact gradients, Hessians, and expected
Fisher information matrices in linear time and storage complexity for the
block full-scale approximation with or without a nugget. In addition, we
discuss and apply an accurate stochastic trace estimation method to estimate
the gradient and full dense expected Fisher matrix using a single pass over
the derivative matrices, which offers significant performance gains in
practice. With these strategies, despite the computational burden of a highly
expensive covariance function, we demonstrate that fitting large datasets with
many-parameter models can be done effectively.

We assume that the mean function of the model is zero everywhere and focus
exclusively on covariance matrices. As in the present application, it is
fairly common to study anomaly fields of climatological variables and treat
these as having mean zero for a mean function specified up to some vector of
unknown linear parameters. The matrix calculus for derivatives with respect to
these parameters is much simpler and poses no computational concern.

##  2 Block Full-Scale Approximation

While a number of past works have approximated the covariance matrix as the
sum of a diagonal matrix and a low-rank matrix [9, 5, 10, 18, 32], these
models often fail to capture short-range behavior of the data. To improve the
model fit, one can add a banded or block diagonal matrix to the low-rank
approximation. We consider here a particular algebraic approximation of this
type referred to as the partially independent conditional approximation by
Snelson and Gharamani [31] or the block full-scale approximation by Sang et
al. [30]. This approximation allows one to modulate the block size and off-
diagonal rank independently to capture both smooth long-range dependence and
rough local dependence while still allowing computations that scale linearly
in the number of observations.

Let $N=\\{1,2,...,n\\}$ indicate the index set of all observations, and take
$I$ and $J$ to be subsets of $N$. Let $\bm{\Sigma}_{IJ}$ denote the submatrix
of $\bm{\Sigma}$ corresponding to rows $I$ and columns $J$. We start by using
a k-d tree to choose a set of $p$ landmark points
$\bm{X}_{P}=\\{\bm{x}_{i}\\}_{i\in P}$ indexed by $P\subset N$ that are
approximately equispaced over the spatial domain, and we construct the low-
rank Nyström approximation
$\bm{\Sigma}_{NP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{NP}^{\top}$.

Next, using the k-d tree, we partition the data into disjoint blocks of
observations at nearby spatial locations. Let block $\ell$ consist of
observations indexed by $B_{\ell}\subset N$ for $\ell=1,...,m$, where $m$ is
the total number of blocks, and denote the collection of these block index
sets as $B=\\{B_{1},...,B_{m}\\}$. For a matrix $\bm{A}\in\mathbb{R}^{n\times
n},$ define the block diagonalization operator with block structure $B$ as

$\Big{[}\operatorname{blkdiag}_{B}(\bm{A})\Big{]}_{ij}=\begin{cases}\bm{A}_{ij}&\text{for
$i,j\in B_{\ell}$ for some $\ell$}\\\ 0&\text{otherwise.}\end{cases}$ (8)

We use this operator to construct a block diagonal correction term that takes
the Nyström approximation to the exact values. This yields our approximate
covariance matrix

$\widetilde{\bm{\Sigma}}=\bm{\Sigma}_{NP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{NP}^{\top}+\operatorname{blkdiag}_{B}\Big{(}\bm{\Sigma}-\bm{\Sigma}_{NP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{NP}^{\top}\Big{)},$
(9)

where the block diagonal structure allows us to capture short-range
covariances exactly within disjoint local neighborhoods. Alternatively, one
can view this approximation as a two-level approximation to the covariance
function given by

$\tilde{k}(\bm{x}_{i},\bm{x}_{j})=\begin{cases}k(\bm{x}_{i},\bm{x}_{j})&\text{for
$i,j\in B_{\ell}$ for some $\ell$}\\\ k(\bm{X}_{P},\bm{x}_{i})^{\top}\
k(\bm{X}_{P},\bm{X}_{P})^{-1}\
k(\bm{X}_{P},\bm{x}_{j})&\text{otherwise.}\end{cases}$ (10)

This approximation is equivalent to assuming that observations in different
neighborhoods are conditionally independent given the observations at the
landmark points $\bm{X}_{P}$ [31]. Considering this conditional structure, one
can also see this approximation as a special case of the Vecchia approximation
[20].

To leverage recent advances in scalable hierarchical matrix operations and
factorizations, Geoga et al. [12] used the Nyström approximation to compress
off-diagonal blocks within a hierarchical off-diagonal low-rank (HODLR)
approximation to the covariance matrix that could be assembled and factorized
in quasilinear time and storage complexity. In that work, a set of $p$
landmark points $\bm{X}_{P}=\\{\bm{x}_{i}\\}_{i\in P}$ are selected from the
data and indexed by $P\subset N$, and each off-diagonal block is approximated
using the low-rank Nyström scheme
$\widetilde{\bm{\Sigma}}_{IJ}=\bm{\Sigma}_{IP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{JP}^{\top}$.
Crucially, the set $P$ of landmark points is fixed and used in every off-
diagonal block. Looking entrywise, this approach induces the two-level matrix
approximation

$\widetilde{\bm{\Sigma}}_{ij}=\begin{cases}\bm{\Sigma}_{ij}&\text{for $i,j\in
B_{\ell}$ for some $\ell$}\\\
\bm{\Sigma}_{Pi}^{\top}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{Pj}&\text{otherwise},\end{cases}$
(11)

where the block diagonal entries are exactly the entries of the full
covariance matrix $\bm{\Sigma}$ and the nonleaf entries are given by a low-
rank approximation. This is precisely the two-level approximation (9). The
simplification from a Nyström-based hierarchical covariance to a two-level
covariance is also discussed by Chen et al. [7].

##  3 Linear Complexity Computations

Here we discuss how to exploit the two-level covariance approximation
described above to obtain direct linear complexity computations of the
likelihood and its derivatives. We derive complexity estimates in terms of the
rank of the Nyström approximation and the block sizes in the block diagonal
term. We can tune these rank and block size parameters to trade off between
approximation accuracy and computational complexity.

###  3.1 Computing the log-likelihood

In order to perform maximum likelihood estimation, our first concern is to
efficiently compute the negative log-likelihood using our approximate
covariance matrix $\widetilde{\bm{\Sigma}}$, given by

$-\ell(\bm{\theta})=\frac{1}{2}\operatorname{logdet}\widetilde{\bm{\Sigma}}+\frac{1}{2}\bm{y}^{\top}\widetilde{\bm{\Sigma}}^{-1}\bm{y}.$
(12)

This requires the computation of the determinant of $\widetilde{\bm{\Sigma}}$
as well as the linear solve $\widetilde{\bm{\Sigma}}^{-1}\bm{y}$. For this
purpose, one would hope to use the matrix determinant lemma

$\det\big{(}\bm{A}+\bm{U}\bm{V}^{\top}\big{)}=\det\big{(}\bm{I}+\bm{V}^{\top}\bm{A}^{-1}\bm{U}\big{)}\det\big{(}\bm{A}\big{)}$
(13)

and the Sherman-Woodbury-Morrison formula

$(\bm{A}+\bm{U}\bm{V}^{\top})^{-1}=\bm{A}^{-1}-\bm{A}^{-1}\bm{U}(\bm{I}+\bm{V}^{\top}\bm{A}^{-1}\bm{U})^{-1}\bm{V}\bm{A}^{-1}$
(14)

which take advantage of the low-rank update structure for
$\bm{A}\in\mathbb{R}^{n\times n}$ and $\bm{U},\bm{V}\in\mathbb{R}^{n\times p}$
to reduce the complexity of these computations. Note, however, that for the
Nyström approximation (9) we have
$\widetilde{\bm{\Sigma}}_{PP}=\bm{\Sigma}_{PP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{PP}^{\top}=\bm{\Sigma}_{PP}$.
In other words, the Nyström approximation is exact on rows and columns
corresponding to landmark points, so
$\operatorname{blkdiag}_{B}\big{(}\bm{\Sigma}-\bm{\Sigma}_{NP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{NP}^{\top}\big{)}$
is zero on these rows and columns. In particular it is not invertible, and
thus we cannot use formulas (13) and (14). Previous works add a nugget
$\sigma^{2}I$ to $\widetilde{\bm{\Sigma}}$, circumventing this issue [31, 30].
However, this is not strictly necessary. We can remedy the rank deficiency
without a nugget by defining a matrix $\bm{\Pi}$ which permutes the landmark
indices $P$ to the last $p$ indices. Let $Q=N\setminus P$ denote the index set
of non-Nyström indices, and let
$B^{\prime}=[B_{1}^{\prime},...,B_{m}^{\prime}]$ be a new collection of block
index sets which partition $Q$ and thus contain no Nyström points. One method
of constructing $B^{\prime}$ is to simply remove any Nyström indices from $B$,
namely $B_{\ell}^{\prime}=B_{\ell}\cap Q$. Then we have

$\widetilde{\bm{\Sigma}}=\bm{\Pi}^{\top}\begin{bmatrix}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}+\operatorname{blkdiag}_{B^{\prime}}\Big{(}\bm{\Sigma}_{QQ}-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}\Big{)}&\bm{\Sigma}_{QP}\\\
\bm{\Sigma}_{QP}^{\top}&\bm{\Sigma}_{PP}\end{bmatrix}\bm{\Pi}.$ (15)

We note that the Nyström rank $p\ll n$, and thus the $(n-p)\times(n-p)$ upper
left block of the permuted matrix contains the vast majority of the covariance
information between observations, and the other blocks are small dense
matrices.

This permuted representation yields efficient and convenient computations,
since the block diagonal correction matrix in the upper left block is now full
rank. For ease of notation, we will write this matrix as

$\bm{D}=\operatorname{blkdiag}_{B^{\prime}}\Big{(}\bm{\Sigma}_{QQ}-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}\Big{)}.$
(16)

Returning to the determinant of $\widetilde{\bm{\Sigma}}$, we see that the
block diagonal correction matrix is the Schur complement of $\bm{\Sigma}_{PP}$
in the permuted matrix, and thus we can compute the determinant of the
approximate covariance matrix using the matrix determinant lemma as

$\det\big{(}\widetilde{\bm{\Sigma}}\big{)}=\det(\bm{D})\det(\bm{\Sigma}_{PP}).$
(17)

Here and for the remainder of this section we will assume for ease of analysis
that the block structure $B^{\prime}$ consists of equally sized blocks of size
$b=\left|B_{\ell}^{\prime}\right|$ for all $\ell=1,...,m$. We will also assume
we have factorized $\bm{\Sigma}_{PP}$ and $\bm{D}$ as a precomputation step
requiring $\mathcal{O}(p^{3}+nb^{2})$ work. Therefore, computing the
determinant of $\widetilde{\bm{\Sigma}}$ requires computing the determinant of
one factorized $p\times p$ matrix and $n/b$ determinants of factorized
$b\times b$ blocks, yielding $\mathcal{O}(n)$ complexity overall.

The linear solve $\widetilde{\bm{\Sigma}}^{-1}\bm{y}$ can also be computed in
a convenient way that leverages Schur complements. Defining the permuted
vector $\bm{z}=\bm{\Pi}\bm{y}$, we solve the permuted linear system

$\displaystyle\begin{bmatrix}\bm{w}_{1}\\\ \bm{w}_{2}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}+\bm{D}&\bm{\Sigma}_{QP}\\\
\bm{\Sigma}_{QP}^{\top}&\bm{\Sigma}_{PP}\end{bmatrix}^{-1}\begin{bmatrix}\bm{z}_{1}\\\
\bm{z}_{2}\end{bmatrix}$ (18) $\displaystyle\bm{w}_{1}$
$\displaystyle=\bm{D}^{-1}\Big{(}\bm{z}_{1}-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{z}_{2}\Big{)}$
$\displaystyle\bm{w}_{2}$
$\displaystyle=\bm{\Sigma}_{PP}^{-1}\Big{(}\bm{z}_{2}-\bm{\Sigma}_{QP}^{\top}\bm{w}_{1}\Big{)}$
$\displaystyle\widetilde{\bm{\Sigma}}^{-1}\bm{y}$
$\displaystyle=\bm{\Pi}^{\top}\bm{w},$

which requires $\mathcal{O}(np+nb)$ work to perform the matrix-vector products
and solve the linear systems. Continuing on to study derivatives of (12) in
much the same way, we show that even matrix-matrix products with this rank
structure can be worked with conveniently and efficiently, from which
scalability of all the required linear algebra follows.

###  3.2 Computing the gradient

To employ gradient-based optimization algorithms for maximum likelihood
estimation, we must compute the gradient of the negative log-likelihood (4).
Each component of the gradient is

$\big{[}-\nabla\ell(\bm{\theta})\big{]}_{j}=\frac{1}{2}\operatorname{tr}\Bigg{[}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bigg{)}\Bigg{]}-\frac{1}{2}\bm{y}^{\top}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bigg{)}\widetilde{\bm{\Sigma}}^{-1}\bm{y}.$
(19)

Alternative rank-structured approximations to $\bm{\Sigma}$ are often
constructed using early-terminating pivoted factorizations and are thus not
differentiable with respect to the kernel parameters. As a result, one must
introduce an additional approximation to the derivative matrices
$\frac{\partial\bm{\Sigma}(\bm{\theta})}{\partial\theta_{j}}$. For example,
the hierarchical matrix method by Minden et al. [26] computes an independent
hierarchical approximation to each
$\bm{\Sigma}(\bm{\theta})^{-1}\frac{\partial\bm{\Sigma}(\bm{\theta})}{\partial\theta_{j}}$
in order to compute the trace term in the gradient of the log-likelihood. In
contrast, Geoga et al. [12] use an algebraic hierarchical covariance matrix
approximation for which the derivatives matrices
$\frac{\partial\bm{\Sigma}(\bm{\theta})}{\partial\theta_{j}}$ can be computed
exactly in quasilinear complexity, but they must use stochastic estimators to
compute the trace term efficiently.

The key observation for efficient computation of the gradient for the block
full-scale approximation is that
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$ has a similar
rank structure to $\widetilde{\bm{\Sigma}}$. Since the Nyström factors
$\bm{\Sigma}_{QP}$ and $\bm{\Sigma}_{PP}$ are simply submatrices of the
covariance matrix $\bm{\Sigma}$, we can compute the derivatives of these
factors using the derivatives of the covariance function. Following basic
matrix differentiation, we have

$\displaystyle\frac{\partial}{\partial\theta_{j}}\Big{(}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}\Big{)}$
$\displaystyle=\bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}$
$\displaystyle-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}$
$\displaystyle+\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}\bigg{)}^{\top}.$
(20)

Since the second term on the right side shares a factor with each of the
others, the derivative of the rank-$p$ Nyström approximation has rank at most
$2p$.

The derivative of the approximate covariance matrix $\widetilde{\bm{\Sigma}}$
is then given by

$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}=\bm{\Pi}^{\top}\begin{bmatrix}\displaystyle\frac{\partial}{\partial\theta_{j}}\Big{(}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}\Big{)}+\frac{\partial\bm{D}}{\partial\theta_{j}}&\displaystyle\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}\\\
\displaystyle\frac{\partial\bm{\Sigma}_{QP}^{\top}}{\partial\theta_{j}}&\displaystyle\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}\end{bmatrix}\bm{\Pi},$
(21)

which has the same block rank structure as $\widetilde{\bm{\Sigma}}$ except
that the low-rank portion of the upper left block has double the rank.

Given this shared rank structure, we can compute the matrix-matrix linear
solve
$\widetilde{\bm{\Sigma}}^{-1}(\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}})$
with Schur complements as follows:

$\displaystyle\begin{bmatrix}\bm{W}_{1}&\bm{W}_{3}\\\
\bm{W}_{2}&\bm{W}_{4}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}+\bm{D}&\bm{\Sigma}_{QP}\\\
\bm{\Sigma}_{QP}^{\top}&\bm{\Sigma}_{PP}\end{bmatrix}^{-1}\begin{bmatrix}\displaystyle\frac{\partial}{\partial\theta_{j}}\Big{(}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}\Big{)}+\frac{\partial\bm{D}}{\partial\theta_{j}}&\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}\\\
\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}^{\top}&\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}\end{bmatrix}$
(22) $\displaystyle\bm{W}_{1}$
$\displaystyle=\bm{D}^{-1}\Bigg{(}\displaystyle\frac{\partial}{\partial\theta_{j}}\Big{(}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}\Big{)}+\frac{\partial\bm{D}}{\partial\theta_{j}}-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}^{\top}\Bigg{)}$
$\displaystyle\bm{W}_{2}$
$\displaystyle=\bm{\Sigma}_{PP}^{-1}\Bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}^{\top}-\bm{\Sigma}_{QP}^{\top}\bm{W}_{1}\Bigg{)}$
$\displaystyle\bm{W}_{3}$
$\displaystyle=\bm{D}^{-1}\Bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}\Bigg{)}$
$\displaystyle\bm{W}_{4}$
$\displaystyle=\bm{\Sigma}_{PP}^{-1}\Bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}-\bm{\Sigma}_{QP}^{\top}\bm{W}_{3}\Bigg{)}$
$\displaystyle\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bigg{)}$
$\displaystyle=\bm{\Pi}\bm{W}\bm{\Pi}^{\top}.$

This requires only linear solves and matrix-matrix products involving $n\times
n$ block diagonal matrices, $n\times p$ low-rank factors, and $p\times p$
matrices, resulting in $\mathcal{O}(np^{2}+nb^{2})$ complexity. In addition,
the resulting matrix
$\widetilde{\bm{\Sigma}}^{-1}(\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}})$
has the same permuted block diagonal plus low-rank structure as
$\widetilde{\bm{\Sigma}}$, where $\bm{W}_{1}\in\mathbb{R}^{(n-p)\times(n-p)}$
is a matrix of rank at most $2p$ plus a block diagonal term and
$\bm{W}_{2}\in\mathbb{R}^{(n-p)\times p}$,
$\bm{W}_{3}\in\mathbb{R}^{p\times(n-p)}$, and
$\bm{W}_{4}\in\mathbb{R}^{p\times p}$ are small dense matrices. The
preservation of this rank structure under linear solves will prove useful for
computing second-order information in the next section.

The remaining inner product term in the gradient of the negative log-
likelihood (19) can be computed in $\mathcal{O}(np+nb)$ time by using equation
(18) along with a straightforward block matrix-vector product. This allows us
to compute each entry of the gradient in linear complexity in $n$.

###  3.3 Computing the Fisher matrix

We can employ the linear solve method above to compute the entries of the
expected Fisher information matrix given by

$\mathcal{I}_{jk}=\frac{1}{2}\operatorname{tr}\Bigg{[}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bigg{)}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{k}}\bigg{)}\Bigg{]}.$
(23)

Computing the terms
$\widetilde{\bm{\Sigma}}^{-1}(\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}})$
and
$\widetilde{\bm{\Sigma}}^{-1}(\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{k}})$
using equation (22), applying a straightforward block matrix-matrix product,
and computing the trace of the resulting rank-structured matrix, we obtain
$\mathcal{O}(np^{2}+nb^{2})$ complexity per entry. Efficient methods for
computing $\mathcal{I}$ facilitate the use of Fisher scoring algorithms to
obtain the MLE and can be used to produce confidence intervals for estimated
parameters. Analogous methods can be used to compute the Hessian of the
negative log-likelihood for use in Newton-based optimization routines. See the
appendix.

###  3.4 Symmetrized trace estimation

Although the above computations of the gradient and the Fisher and Hessian
matrices have the desired linear scaling in the data size $n$, the matrix-
matrix solves and products in the trace terms require a large number of
intermediate allocations, which are computationally expensive in practice. To
provide fast unbiased estimates of these trace terms, we rely on a sample
average approximation based on the Hutchinson estimator [16]

$\operatorname{tr}(\bm{A})\approx\frac{1}{s}\sum_{\ell=1}^{s}\bm{u}_{\ell}^{\top}\bm{A}\bm{u}_{\ell}$
(24)

with $s$ samples, where $\bm{u}_{\ell}$ are independent symmetric Bernoulli
vectors, although alternatives exist (see [35]). In particular, factorizing
the approximate covariance as $\widetilde{\bm{\Sigma}}=\bm{W}\bm{W}^{\top}$,
one can use the symmetrized estimator presented by Stein et al. [35], which is
given by

$\operatorname{tr}\Bigg{[}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bigg{)}\Bigg{]}=\operatorname{tr}\Bigg{[}\bm{W}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bigg{)}\bm{W}^{-\top}\Bigg{]}\approx\frac{1}{s}\sum_{\ell=1}^{s}\bm{u}_{\ell}^{\top}\bm{W}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bigg{)}\bm{W}^{-\top}\bm{u}_{\ell}.$
(25)

Stein et al. [35] prove a bound on the variance of this estimator that is at
least as strong as the variance bound on the nonsymmetrized version, and Geoga
et al. [12] provide numerical results that indicate greatly improved accuracy
with the symmetrized estimator compared with the standard Hutchinson
procedure.

An analogous symmetrized estimator for the trace terms in the Fisher matrix
and Hessian that require only a small number of linear solves with $\bm{W}$
and matrix-vector products with
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$ is constructed in
[12]. The Fisher matrix estimator can be written as

$\mathcal{I}_{jk}\approx\frac{1}{4s}\sum_{\ell=1}^{s}\bm{u}_{\ell}^{\top}\bm{W}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}+\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{k}}\bigg{)}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}+\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{k}}\bigg{)}\bm{W}^{-\top}\bm{u}_{\ell}-\frac{1}{2}\mathcal{I}_{jj}-\frac{1}{2}\mathcal{I}_{kk},$
(26)

where the diagonal terms $\mathcal{I}_{jj}$ and $\mathcal{I}_{kk}$ can be
estimated in a trivially symmetric way. This gives fast stochastic estimators
of all quantities necessary for gradient-based and second-order optimization
solvers, which can be computed in $\mathcal{O}(snp+snb)$ complexity per entry
since matrix-vector products with
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$ and linear solves
with $\bm{W}$ require $\mathcal{O}(np+nb)$ work, which we now show.

###  3.5 Symmetric factor computation

The remaining concern is to obtain a rank-structured symmetric factor $\bm{W}$
for our two-level approximate covariance matrix $\widetilde{\bm{\Sigma}}$. The
principal challenge is the upper left block, which we hope to factorize as

$\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}+\bm{D}=\Big{(}\bm{X}\bm{Y}^{\top}+\bm{B}\Big{)}\Big{(}\bm{X}\bm{Y}^{\top}+\bm{B}\Big{)}^{\top}$
(27)

for $\bm{X},\bm{Y}\in\mathbb{R}^{(n-p)\times p}$ and $\bm{B}$ a block diagonal
matrix with the same block structure $B^{\prime}$ as $\bm{D}$. As discussed in
Section 2, our matrix structure is a two-level special case of the HODLR
format; thus we use a single step of the symmetric factorization algorithm of
Ambikasaran et al. [3], which can be written concisely by computing Cholesky
factors

$\displaystyle\bm{B}\bm{B}^{\top}$ $\displaystyle=\bm{D}$ (28)
$\displaystyle\bm{L}\bm{L}^{\top}$
$\displaystyle=\big{(}\bm{B}^{-1}\bm{\Sigma}_{QP}\big{)}^{\top}\big{(}\bm{B}^{-1}\bm{\Sigma}_{QP}\big{)}$
(29) $\displaystyle\bm{M}\bm{M}^{\top}$
$\displaystyle=\bm{I}+\bm{L}^{\top}\bm{\Sigma}_{PP}^{-1}\bm{L}$ (30)
$\displaystyle\bm{X}$ $\displaystyle=\bm{\Sigma}_{QP}$ (31)
$\displaystyle\bm{Y}$
$\displaystyle=\bm{L}^{-\top}\big{(}\bm{M}-\bm{I}\big{)}\bm{L}^{-1}\big{(}\bm{B}^{-1}\bm{\Sigma}_{QP}\big{)}^{\top}.$
(32)

This requires only Cholesky factorizations of block diagonal matrices and
$p\times p$ matrices, as well as linear solves and matrix-matrix products
involving $p\times p$ and $(n-p)\times p$ matrices, and thus can be computed
in $\mathcal{O}(np^{2}+nb^{2})$ time. The symmetric factor $\bm{W}$ is then
given by

$\bm{W}=\bm{\Pi}^{\top}\begin{bmatrix}\bm{X}\bm{Y}^{\top}+\bm{B}&0\\\
\bm{Z}^{\top}&\bm{G}\end{bmatrix}\bm{\Pi}$ (33)

where the remaining blocks are defined by

$\displaystyle\bm{Z}$
$\displaystyle=\Big{(}\bm{X}\bm{Y}^{\top}+\bm{B}\Big{)}^{-1}\bm{\Sigma}_{QP}$
(34) $\displaystyle\bm{G}\bm{G}^{\top}$
$\displaystyle=\bm{\Sigma}_{PP}-\bm{Z}^{\top}\bm{Z}$ (35)

with $\bm{Z}\in\mathbb{R}^{(n-p)\times p},\bm{G}\in\mathbb{R}^{p\times p}$.
Since $\bm{W}$ maintains the same permuted block diagonal plus low-rank
structure as $\widetilde{\bm{\Sigma}}$, we can compute matrix-vector products
and linear solves with $\bm{W}$ in $\mathcal{O}(np+nb)$ time using the
Sherman-Woodbury-Morrison formula. This facilitates the SAA approximations to
the gradient, Fisher matrix, and Hessian discussed above and provides a fast
sampling method for the process.

###  3.6 Prediction and conditional distributions

In addition to fast likelihood computations, the rank structure of
$\widetilde{\bm{\Sigma}}$ facilitates fast kriging and results in a
conditional covariance matrix with the same rank structure as
$\widetilde{\bm{\Sigma}}$. Given a set of locations
$\\{\bm{x}_{i}^{*}\\}_{i=1}^{n_{*}}$ indexed by $N_{*}$ with corresponding
process values $y_{i}^{*}=Z(\bm{x}_{i}^{*})$, we define the vector
$\bm{y}_{*}=[y_{1}^{*},...,y_{n_{*}}^{*}]$ and wish to compute the conditional
distribution

$\bm{y}_{*}|\bm{y}\sim
N(\widetilde{\bm{\Sigma}}_{*}^{\top}\widetilde{\bm{\Sigma}}^{-1}\bm{y},\widetilde{\bm{\Sigma}}_{**}-\widetilde{\bm{\Sigma}}_{*}^{\top}\widetilde{\bm{\Sigma}}^{-1}\widetilde{\bm{\Sigma}}_{*}).$
(36)

where
$\big{(}\widetilde{\bm{\Sigma}}_{*}\big{)}_{ij}=\tilde{k}(\bm{x}_{i},\bm{x}^{*}_{j})$
and
$\big{(}\widetilde{\bm{\Sigma}}_{**}\big{)}_{ij}=\tilde{k}(\bm{x}^{*}_{i},\bm{x}^{*}_{j})$.
We assign each point $\bm{x}_{i}^{*}$ to a block from the observed data, for
example by taking the block with centroid nearest to $\bm{x}_{i}^{*}$,
resulting in a collection of block index sets
$B_{*}=[B_{1}^{*},...,B_{m}^{*}]$ that partition $N_{*}$. To simplify the
discussion of computational complexity, we assume this results in blocks of
equal size $b_{*}=\left|B_{\ell}^{*}\right|$. Recalling the permutation
$\bm{\Pi}$ that orders the landmark points to the last indices, we obtain
covariance matrices of the form

$\displaystyle\widetilde{\bm{\Sigma}}_{**}$
$\displaystyle=\bm{\Sigma}_{N_{*}P}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{N_{*}P}^{\top}+\operatorname{blkdiag}_{B_{*}}\Big{(}\bm{\Sigma}_{N_{*}N_{*}}-\bm{\Sigma}_{N_{*}P}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{N_{*}P}^{\top}\Big{)}$
(37) $\displaystyle\widetilde{\bm{\Sigma}}_{*}$
$\displaystyle=\bm{\Pi}^{\top}\begin{bmatrix}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{N_{*}P}^{\top}+\operatorname{blkdiag}_{B^{\prime}B_{*}}\Big{(}\bm{\Sigma}_{QN_{*}}-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{N_{*}P}^{\top}\Big{)}\\\
\bm{\Sigma}_{N_{*}P}^{\top}\end{bmatrix},$ (38)

where the nonsymmetric block diagonalization operator is defined by

$\Big{[}\operatorname{blkdiag}_{B^{\prime}B_{*}}(\bm{A})\Big{]}_{ij}=\begin{cases}\bm{A}_{ij}&\text{for
$i\in B^{\prime}_{\ell}$ and $j\in B_{\ell}^{*}$ for some $\ell$}\\\
0&\text{otherwise.}\end{cases}$ (39)

We see that these covariance matrices have rank structures that are minor
variations on the block diagonal plus low-rank structure we have used in the
observed data covariance matrix, its derivatives, and its symmetric factor.

To compute the term $\widetilde{\bm{\Sigma}}^{-1}\bm{y}$ in the conditional
mean, we use the linear solve (18) followed by a straightforward block matrix-
vector product with $\widetilde{\bm{\Sigma}}_{*}$, which has complexity
$\mathcal{O}(nb_{*}+np+n_{*}p)$.

To compute the term $\widetilde{\bm{\Sigma}}^{-1}\widetilde{\bm{\Sigma}}_{*}$
in the conditional covariance, we use the first column of the structured
matrix-matrix solve (22). This yields a matrix with the same block structure
as (38). We then compute the block matrix-matrix product with
$\widetilde{\bm{\Sigma}}_{*}^{\top}$, which has complexity
$\mathcal{O}(n\,b_{*}^{2}+np^{2}+n_{*}p)$. Importantly, the resulting
conditional covariance matrix is a block diagonal plus a rank-$4p$ matrix.
Thus we can afford to compute and store it, facilitating further computations
such as symmetric factorization using (28) and thus yielding conditional
simulations in linear complexity.

###  3.7 Numerical verification of $\mathcal{O}(n)$ complexity

Before applying the methods developed above to a nonstationary process, we
demonstrate their linear scaling using a simple stationary process. We
consider fitting $n$ observations at locations selected uniformly at random in
$[0,1]^{2}$ using the stationary isotropic Matérn covariance (5) with fixed
$\nu=1$, where $\bm{\theta}=[\sigma^{2},\rho]$ are being estimated. We fix the
block size $b=128$ and the Nyström rank $p=32$. For various $n$ we then time
the computation of the approximate covariance matrix
$\widetilde{\bm{\Sigma}}$, the likelihood $\ell$, the symmetric factor
$\bm{W}$, the gradient $\nabla\ell$, and the Fisher information matrix
$\mathcal{I}$, as well as computation of the conditional distribution of
$\bm{y}_{*}|\bm{y}$ consisting of the mean and rank-structured conditional
covariance matrix at 512 sites also selected uniformly at random in
$[0,1]^{2}$. Figure 1 shows that all the aforementioned computations scale
linearly as expected and gives timings for our implementations on a single
core of an Intel Xeon CPU E5-2650 @ 2.00 GHz machine.

$10^{-4}$$10^{-2}$$10^{0}$$2^{9}$$2^{10}$$2^{11}$$2^{12}$$2^{13}$$2^{14}$$2^{15}$$2^{16}$$2^{17}$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$2^{9}$$2^{10}$$2^{11}$$2^{12}$$2^{13}$$2^{14}$$2^{15}$$2^{16}$$2^{17}$$\widetilde{\bm{\Sigma}}$$\ell$$\bm{W}$$\nabla\ell$$\mathcal{I}$$\bm{y}_{*}\,|\,\bm{y}$
Figure 1: Linear scaling of covariance matrix construction, symmetric
factorization, and log-likelihood evaluation (left); and of gradient, Fisher
matrix, and conditional distribution computation at 512 unobserved locations
(right). Dotted black lines show $\mathcal{O}(n)$ scaling for reference.

##  4 An Application at Scale

To test our computational framework in the relevant setting of nonstationary
modeling with many parameters, we study a large sea surface temperature
anomaly dataset from the NOAA Coral Reef Watch database consisting of a
40°$\times$ 40° domain in the central Pacific Ocean. These reanalysis data are
de-meaned sea surface temperature measurements generated by interpolating
polar-orbiting and geostationary satellite data from multiple sources to a
$0.05^{\circ}$ ($\approx 5$ km) grid using a Kalman filter-like approach [25,
21]. To study such a large spatial domain, we subsample this data to a 0.1°
grid. We then use local averaging on NASA’s MODIS Cloud Mask product [1] to
compute a holdout set on the 0.1°grid, which will serve as a testing set for
evaluating predictions and uncertainties. This leaves 107,600 non-cloudy
observations for maximum likelihood estimation. See Figure 2. These cloud-
masked data provide a realistic setting for interpolation in atmospheric
science applications.

-15-71917-15-8-161320-115-108-101-94-87-115-108-101-94-87-115-108-101-94-87

latitude

-5-4-3-2-1012

latitude

longitudelongitude-2-1.5-1-0.500.511.52longitude00.050.10.150.20.250.3

Figure 2: Artificially cloud-masked data (top left), interpolant (top center),
full data (top right), neighborhoods used for local fitting and approximation
blocks with RBF centers as white dots (bottom left), interpolant error (bottom
center), and interpolant standard deviations (bottom right).

Recall the nonstationary anisotropic model described in Section 1, which is
given by

$k(\bm{x}_{i},\bm{x}_{j})=\sigma^{2}\frac{\left|\bm{\Lambda}(\bm{x}_{i})\right|^{\frac{1}{4}}\left|\bm{\Lambda}(\bm{x}_{j})\right|^{\frac{1}{4}}}{\left|\frac{\bm{\Lambda}(\bm{x}_{i})+\bm{\Lambda}(\bm{x}_{j})}{2}\right|^{\frac{1}{2}}}\mathcal{M}_{\nu}\Bigg{(}\sqrt{(\bm{x}_{i}-\bm{x}_{j})^{\top}\bigg{(}\frac{\bm{\Lambda}(\bm{x}_{i})+\bm{\Lambda}(\bm{x}_{j})}{2}\bigg{)}^{-1}(\bm{x}_{i}-\bm{x}_{j})}\Bigg{)}.$
(6)

To estimate the spatially-varying anisotropy function $\bm{\Lambda}(\cdot)$,
we expand it in a normalized radial basis

$\displaystyle\bm{\Lambda}(\bm{x})$
$\displaystyle=\sum_{i=1}^{m}\phi(\left\lVert\bm{x}-\bm{a}_{i}\right\rVert)\
\bm{\Lambda}_{i}$ (40) $\displaystyle\phi(r)$
$\displaystyle=\frac{\psi(r)}{\sum_{i=1}^{m}\psi(r)}$ (41)
$\displaystyle\psi(r)$ $\displaystyle=e^{-(r/c)^{2}}$ (42)

using squared exponential bases $\psi$, and we estimate the positive definite
matrices $\bm{\Lambda}_{i}$. In particular, we express the anisotropy matrices
in terms of their Cholesky factors

$\displaystyle\bm{\Lambda}_{i}$ $\displaystyle=\bm{L}_{i}\bm{L}_{i}^{\top}$
(43) $\displaystyle\bm{L}_{i}$
$\displaystyle=\begin{bmatrix}\ell_{i}^{(1,1)}&0\\\
\ell_{i}^{(2,1)}&\ell_{i}^{(2,2)}\end{bmatrix},$ (44)

which guarantees positive definiteness. We find also that it is necessary to
parameterize the log of the diagonal elements in order to enforce uniqueness
of the Cholesky factors and avoid identifiability issues. This results in the
parameter vector

$\bm{\theta}=[\log\ell_{1}^{(1,1)},\ell_{1}^{(2,1)},\log\ell_{1}^{(2,2)},...,\log\ell_{m}^{(1,1)},\ell_{m}^{(2,1)},\log\ell_{m}^{(2,2)}]\in\mathbb{R}^{3m}.$
(45)

Since the nonnegative linear combination of positive definite matrices is
positive definite, we have obtained a parameterization in which
$\bm{\Lambda}(\bm{x})$ is positive definite for all $\bm{x}\in\Omega$. We fix
the smoothness parameter $\nu=1$ because it is difficult to estimate
simultaneously with the local anisotropies, and we are interested here
primarily in the nonstationary anisotropy parameters. We estimate the scale
parameter $\sigma^{2}$ using the profile likelihood by fixing $\sigma^{2}=1$,
estimating the anisotropy parameters $\bm{\theta}$ and then computing the
optimal scale parameter, which is given in closed form by
$\sigma^{2}=\frac{1}{n}\bm{y}^{\top}\widetilde{\bm{\Sigma}}(1,\bm{\theta})^{-1}\bm{y}$.

We start by partitioning our data into 64 disjoint subregions using a k-d
tree, and we use the centroid of each of these regions as the center of an RBF
within the full model, as shown in Figure 2. To determine a suitable initial
guess for our global model optimization and to compare against the current
state of the art in which locally fitted parameters are plugged into the RBF
expansion (40), we fit a local anisotropy $\bm{\Lambda}_{i}$ in each
subregion. Since each subregion contains only 1681 observations, we can afford
to use exact linear algebra and do not require any covariance matrix
approximation.

To fit these local models as well as to fit the global model we discuss below,
we use our own implementation of the trust-region algorithm adapted directly
from [38]. At iteration $k$ this algorithm minimizes a quadratic approximation
to the negative log-likelihood

$\min_{\bm{p}\in\mathbb{R}^{3m}}-\ell(\bm{\theta}^{(k)})-\nabla\ell(\bm{\theta}^{(k)})^{\top}\bm{p}-\frac{1}{2}\bm{p}^{\top}\bm{B}_{k}\bm{p}\quad\text{s.t.}\quad\left\lVert\bm{p}\right\rVert\leq\Delta_{k}$
(46)

where $\bm{\theta}^{(k)}$ are the current parameters, $\bm{B}_{k}$ is an
approximation to the Hessian $\nabla^{2}\ell(\bm{\theta}^{(k)})$, and
$\Delta_{k}$ is a radius parameter which indicates the size of the region in
that this quadratic objective is a good approximation to the true negative
log-likelihood. Trust-region algorithms are known to converge to stationary
points for various approximate solutions of the subproblem (46) and for
various Hessian approximations $\bm{B}_{k}$ as long as $\bm{B}_{k}$ is
bounded. In the spirit of Fisher scoring, we use the symmetrized stochastic
estimate of the full Fisher matrix $\bm{B}_{k}=\mathcal{I}(\bm{\theta}^{(k)})$
discussed in Section 3.4 and solve the subproblem (46) using a Newton-based
iterative method ([38], Section 4.3) which is inexpensive and effective for
problems of this size.

Following the local fitting of anisotropy parameters, we optimize the full set
of global model parameters $\bm{\theta}$. Here we consider the full data and
thus require the block full-scale approximation and accompanying computations
discussed in Sections 2 and 3. We use the local neighborhoods around each RBF
as the blocks in the approximation and use the k-d tree to select 72 Nyström
points for the low-rank portion that are approximately equispaced over the
spatial domain. We can therefore think of the global model as a combination of
the exact local neighborhoods on which the local parameters were estimated,
plus a low-rank term coupling these neighborhoods. Comparing the likelihood
under a disjoint neighborhood model using the locally estimated anisotropy
parameters (i.e. a block diagonal covariance) with the likelihood under the
Paciorek-Schervish model with the block diagonal plus low-rank covariance
approximation using the locally estimated anisotropy parameters as the
$\Lambda_{i}$ in equation (40), we observe a log-likelihood increase of
334,342 units. This large likelihood improvement indicates that observations
in adjacent neighborhoods are in fact highly correlated and that the low-rank
Nyström term captures some of this dependence.

Preparing for optimization of the global model parameters, we choose the width
parameter $c$ of the radial basis function $\psi$ to be half of the minimum
distance between RBF centers. This choice results in bases that are fairly
localized, helping avoid identifiability problems between adjacent RBFs caused
by more diffuse bases with larger $c$ values. In principle, estimating the
full parameter vector $\bm{\theta}$ allows one to also perform maximum
likelihood for the parameter $c$, which is impossible with local estimation.
We find, however, that if $c$ is included as a parameter alongside the locally
estimated parameters, the derivative in $c$ is much larger than the
derivatives in the anisotropy parameters in $\bm{\theta}$, and thus the solver
moves in the $c$ direction without materially improving the fit. If we instead
fix a localized $c$ and estimate $\bm{\theta}$ only, we find that the
likelihood surface in $c$ becomes uninformative because the basis width does
not have a large impact on the model fit if the local parameters are well
selected.

-15-8-161320-15-8-161320-115-108-101-94-87-115-108-101-94-87-115-108-101-94-87

latitude

latitude

longitudelongitudelongitude

Figure 3: Correlation functions at three points using parameters fit globally
(top row) and locally (bottom row), with contours at levels $0.1,0.2,...,0.9$.

After making the block full-scale approximation and selecting the basis width,
the key to fitting these 192 parameters in practice using the trust-region-
based Fisher scoring algorithm described above lies in the use of the SAA
methods discussed in Section 3 to compute the gradient and Fisher matrix in a
single pass over the derivative matrices
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$. The crucial
observation is that the symmetrized estimators for both the gradient and
Fisher entries (25, 26) require only inner products of terms of the form
$\bm{W}^{-\top}\bm{u}_{\ell}$ and
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bm{W}^{-\top}\bm{u}_{\ell}$.
Thus, we can compute each
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$ and matrix-vector
products with $\bm{W}^{-\top}\bm{u}_{\ell}$ efficiently and independently in
parallel. We can then compute the necessary inner products for each entry in
the gradient and Fisher matrix. We use 150 SAA vectors for this numerical
experiment, which leads to a dramatic reduction both in the computational
cost, because the inner products are less expensive than matrix-matrix
products, and in memory, because the vectors
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bm{W}^{-\top}\bm{u}_{\ell}$
are significantly smaller than the full derivative matrices
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$. While both exact
and SAA computations require $\mathcal{O}(n)$ time and storage, SAA provides a
dramatic reduction in prefactors. On an Intel Xeon CPU E5-2650 @ 2.00 GHz
machine, computing each Fisher matrix entry using the exact form takes about
40 seconds in our implementation after the matrices
$\widetilde{\bm{\Sigma}}^{-1}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$
have been computed, whereas each entry takes 0.05 seconds using SAA after the
vectors
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bm{W}^{-\top}\bm{u}_{\ell}$
have been computed. As a result, the full exact Fisher matrix requires about
206 hours, whereas the SAA approximation requires only 15 minutes. In our
implementation, parallelizing the construction of the derivatives of the
covariance matrix $\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$
on 16 cores leads to approximately 80 seconds for each of the 192 derivative
matrices, resulting in a total of about 4 hours. This easily dominates the
computational cost at each iteration, far more demanding than computing the
SAA gradient and Fisher matrix.

Comparing the log-likelihood of the global Paciorek-Schervish model using the
initial, locally fitted parameters with the log-likelihood of the same model
using the final, global parameters given by our Fisher scoring trust-region
algorithm, we observe an increase of 22,565 units. In total, parallelized over
16 cores, the time to optimize the local models in disjoint subregions was
about 7 hours, and the time to optimize the global model which stagnates after
16 iterations, was about 74 hours. The substantial likelihood difference
indicates that solving the high-dimensional optimization problem to fit all
parameters simultaneously can produce a significant improvement in model fit
when compared with purely local parameter estimation. The difference in
covariance structure can be seen in Figure 3, which shows the estimated
correlation function obtained using local parameter estimates versus optimized
global parameter estimates. In particular, we notice that the global model
fitting captures larger east-west correlations above and below the vortices
caused by equatorial currents and little correlation across these currents.

After maximum likelihood estimation is completed using the inexpensive SAA
gradients and Fisher matrices, we can efficiently compute interpolants and the
rank-structured conditional covariance matrix. Figure 2 shows the
interpolation results along with standard errors. In addition, we compute the
more expensive exact Fisher matrix at the MLE to evaluate parameter
uncertainties. We see that the inverse Fisher matrix contains non-negligible
terms away from the diagonal, indicating that some interaction exists between
spatially proximal parameters. See Figure 4, which shows the correlation
matrix given by normalizing the inverse Fisher matrix to have unit diagonal
entries. If one were to fit parameters only locally, the resulting Fisher
matrix would be block diagonal with $3\times 3$ blocks, which ignores these
off-diagonal contributions (also shown in Figure 4). Beyond the significantly
improved likelihoods, this additional second-order information about parameter
point estimates themselves is a material advantage of the global model
formulation used here. Interpolation can be nearly optimal even with a
completely misspecified covariance function [33] and may not be significantly
improved by a global model. We find negligible differences in the mean squared
error of interpolants computed in disjoint local neighborhoods with stationary
models and the interpolants shown in Figure 2 that are computed with the
global nonstationary model. However, uncertainties for estimated parameters
_can_ be significantly underestimated when global dependence is not accounted
for. Especially in an RBF model such as the one presented here, the MLE and
expected Fisher matrix from the global model may also serve as a reasonable
approximation to a Bayesian posterior, although one needs to be careful when
appealing to asymptotic results that depend on the consistency of parameter
estimates when working with spatial data [39]. The difference in parameter
correlation structure can be seen in Figure 5, which shows the correlation
between the log of the upper left Cholesky entry $\log\ell_{i}^{(1,1)}$ (which
serves as one of the nonstationary model parameters, see (45)) at various
spatial locations for the global Paciorek-Schervish model and the disjoint
local neighborhood model. Note the meaningful negative correlations between
adjacent parameters in some regions when using the global model, which cannot
be captured by disjoint local models.

-0.3-0.2-0.100.10.20.3-0.3-0.2-0.100.10.20.3-10-8-6-4-2024

Figure 4: Inverse Fisher matrix normalized to have unit diagonal entries for
the disjoint local neighborhood model (left) and for the global Paciorek-
Schervish model (center) shown with a reduced color range to emphasize off-
diagonal elements, as well as the logs of the magnitudes of the entries of the
latter Fisher matrix (right).

-15-8-161320-15-8-161320-115-108-101-94-87-115-108-101-94-87-115-108-101-94-87

latitude

-1-0.500.51

latitude

longitudelongitudelongitude-1-0.500.51

Figure 5: Spatially indexed parameter correlation for $\log\ell_{i}^{(1,1)}$
using the global Paciorek-Schervish model (top row) and the disjoint local
neighborhood model (bottom row).

To further identify where significant likelihood improvements have been made,
we investigate quality of fit at the lower tail of the spectrum of the
covariance matrix. If a covariance matrix $\bm{\Sigma}\in\mathbb{R}^{n\times
n}$ has orthonormal eigenvectors $\bm{q}_{j}$ and corresponding positive
eigenvalues $\lambda_{j}$, then a vector
$\bm{z}\sim\mathcal{N}(\bm{0},\bm{\Sigma})$ can be represented as
$\bm{z}\sim\sum_{j=1}^{n}\sqrt{\lambda_{j}}\varepsilon_{j}\bm{q}_{j},$ where
$\varepsilon_{j}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathcal{N}(0,1)$.
By virtue of the orthogonality of the eigenvectors, it follows that under this
model $\bm{q}_{j}^{T}\bm{z}\sim\mathcal{N}(0,\lambda_{j})$, and
$\bm{q}_{j}^{T}\bm{z}$ is independent of $\bm{q}_{k}^{T}\bm{z}$ for all $j\neq
k$.

Because $\bm{\Sigma}$ is inverted in the quadratic form that appears in the
likelihood, it is particularly revealing of log-likelihood improvements to
inspect these inner products for eigenvectors corresponding to the smallest
eigenvalues of $\bm{\Sigma}$. Figure 6 shows the quantity

$Z_{j}=\left|\frac{\bm{q}_{j}(\bm{\theta})^{T}\bm{z}}{\sqrt{\lambda_{j}(\bm{\theta})}}\right|,$
(47)

the absolute value of a standard $Z$-score, for the eigenpairs corresponding
to the 100 smallest eigenvalues using the model parameters $\bm{\theta}$
estimated locally in disjoint regions and simultaneously in the global model.

00.511.522.533.544.50.511.522.50.511.522.533.544.50510152025303540

Empirical quantiles

Half-normal quantiles

count

$|Z|$

Figure 6: Absolute values of standardized $Z$-scores for the eigenvector inner
products given by equation (47) corresponding to the 100 smallest eigenvalues
for the global model covariance matrix using plug-in local parameter estimates
(red) and globally optimized parameters (blue).

As can immediately be seen, even after taking into account small-sample
variability, the $Z$-scores using the local plug-in parameters are excessively
dispersed. This reflects a significant misfit, and the corresponding values
for the globally optimized parameters much more closely resemble the standard
half-normal density. Because the smallest eigenvectors for most standard
covariance models generally reflect high-frequency energy and are comparable
to some form of higher-order differencing, this indicates that the global
optimization leads to the model capturing fundamentally local behavior better
than what is achieved with plug-in local parameters. Since the model is
necessarily somewhat misspecified, this behavior is likely explained at least
in part as the global optimization sacrificing fit quality in some parts to
make significant enough improvement in other parts such that the total
likelihood is improved. Considering that almost every model used for serious
applications with environmental data will be to some degree misspecified, we
find this to be a valuable functionality that is not possible with plug-in
local parameters.

Overall, we see that the combination of the covariance matrix approximation
and second-order trust-region-based Fisher scoring algorithm presented here,
which scale favorably with the data size and the number of parameters, can
produce meaningfully different fits and richer parameter distributions for
many-parameter models. Therefore these methods represent an important
consideration when modeling nonstationary data at scale.

##  5 Discussion

In this work we present a complete set of linear cost computations for second-
order maximum likelihood estimation of Gaussian process parameters using the
block full-scale approximation with application to fitting many-parameter
nonstationary models. In particular, we give exact algorithms for computing
the Gaussian log-likelihood and its gradient, Fisher information matrix, and
Hessian, as well as methods for highly efficient stochastic approximation. The
ability to compute derivatives of the log-likelihood is crucial especially in
fitting expressive covariance models with many parameters, in which solvers
must navigate a complex, high-dimensional, nonconvex objective landscape. We
demonstrate using a large sea surface temperature dataset that our methods
facilitate parameter estimation for nonstationary models with a very large
number of parameters. Further, we demonstrate the value of complex global
parameterizations by the significantly improved likelihood and additional
inferred covariance structure between parameters, motivating the need for
methods that are designed to be scalable with respect to both data size and
parameter size.

The block diagonal plus low-rank approximation we use is well known in the
literature as the partially independent conditional [31] or block full-scale
approximation [30]. It is also a special case of the Vecchia approximation
[20] and is a two-level case of the more general hierarchical models of both
Katzfuss [17, 19] and Chen [8]. While these more sophisticated approaches can
achieve more accurate approximations to the log-likelihood [19, 8], they
provide no efficient methods for computing its derivatives, making high-
dimensional parameter estimation extremely challenging. Our principal
contributions are to demonstrate that the block full-scale approximation
yields efficient computations with and without a nugget using Schur
complements in a permuted covariance matrix and to apply these computations to
high-dimensional parameter estimation problems. The utility and efficiency of
computations with the block diagonal plus low-rank structure are due to the
fact that this structure is closed with respect to matrix-matrix addition,
multiplication, and inversion. Therefore, by making a single algebraic
approximation to the covariance matrix, we obtain rank-structured derivatives
matrices, symmetric factors, and conditional covariance matrices that can be
assembled using only $\mathcal{O}(n)$ computation and storage.

Although we find that this covariance matrix approximation in conjunction with
the SAA gradient and Fisher matrix methods makes parameter estimation possible
for hundreds of parameters, it has limitations in some circumstances. With
respect to the covariance approximation, Gong et al. [19] show that the
conditional mean can suffer from discontinuities at block boundaries, which
are noticeable upon close inspection of the interpolant in Figure 2. Although
these can in theory be alleviated by tapering at block boundaries with a
sufficiently smooth compactly supported function, this comes at the cost of
approximation accuracy. Regarding SAA and parameter dimension, the full Fisher
matrix contains $\mathcal{O}(m^{2})$ many entries, where $m$ is the number of
parameters. Due to the incredible efficiency of computing SAA inner products,
this is not a computational bottleneck in the application shown above.
However, if one were to need thousands of parameters with hundreds of
thousands of observations, or if exact Fisher entries were required for
extremely high-precision optimization, this quadratic scaling in parameter
dimension may become prohibitive. In this regime, some form of structured
quasi-Newton that approximates the Fisher matrix using fewer entries may
become necessary. Finally, we note that the cost of evaluating the basis
function expansion (40) to compute the covariance itself has cost
$\mathcal{O}(m)$, which may become a computational bottleneck for very large
datasets with an extremely large number of basis functions. This could
potentially be alleviated through the use of compactly supported basis
functions as in [15], although these authors use purely local parameter
estimation in place of the global model estimation we have developed here.

## Acknowledgment

This material was based upon work supported by the US Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research (ASCR)
under Contract DE-AC02-06CH11347.

## Appendix

### Computing the Hessian

To employ second-order Newton solvers instead of Fisher scoring to compute the
MLE, one must compute the Hessian, whose entries are given by

$\displaystyle\big{[}-\nabla^{2}\ell(\bm{\theta})\big{]}_{jk}$
$\displaystyle=-\mathcal{I}_{jk}+\frac{1}{2}\operatorname{tr}\Bigg{[}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial^{2}\widetilde{\bm{\Sigma}}}{\partial\theta_{j}\partial\theta_{k}}\bigg{)}\Bigg{]}$
$\displaystyle+\bm{y}^{\top}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}\bigg{)}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{k}}\bigg{)}\widetilde{\bm{\Sigma}}^{-1}\bm{y}$
(48)
$\displaystyle-\frac{1}{2}\bm{y}^{\top}\widetilde{\bm{\Sigma}}^{-1}\bigg{(}\frac{\partial^{2}\widetilde{\bm{\Sigma}}}{\partial\theta_{j}\partial\theta_{k}}\bigg{)}\widetilde{\bm{\Sigma}}^{-1}\bm{y}.$

We can again apply basic matrix differentiation rules to equation (3.2) to
obtain the second derivatives of the Nyström approximation, where the second
derivative of the rank-$p$ Nyström approximation

$\displaystyle\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}\Big{(}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}\Big{)}$
$\displaystyle=\bigg{(}\frac{\partial^{2}\bm{\Sigma}_{QP}}{\partial\theta_{j}\partial\theta_{k}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}$
(49)
$\displaystyle-\bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{k}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}$
$\displaystyle+\bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{QP}^{\top}}{\partial\theta_{k}}\bigg{)}$
$\displaystyle-\bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{k}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}$
$\displaystyle+\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{k}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}$
$\displaystyle-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial^{2}\bm{\Sigma}_{PP}}{\partial\theta_{j}\partial\theta_{k}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}$
$\displaystyle+\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{k}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}$
$\displaystyle-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{j}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{QP}^{\top}}{\partial\theta_{k}}\bigg{)}$
$\displaystyle+\bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{k}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}\bigg{)}^{\top}$
$\displaystyle-\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{PP}}{\partial\theta_{k}}\bigg{)}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial\bm{\Sigma}_{QP}}{\partial\theta_{j}}\bigg{)}^{\top}$
$\displaystyle+\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bigg{(}\frac{\partial^{2}\bm{\Sigma}_{QP}}{\partial\theta_{j}\partial\theta_{k}}\bigg{)}^{\top}$

has rank at most $4p$. The second derivative of the approximate covariance
matrix $\widetilde{\bm{\Sigma}}$ is then given by

$\frac{\partial^{2}\widetilde{\bm{\Sigma}}}{\partial\theta_{j}\partial\theta_{k}}=\bm{\Pi}^{\top}\begin{bmatrix}\displaystyle\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}\Big{(}\bm{\Sigma}_{QP}\bm{\Sigma}_{PP}^{-1}\bm{\Sigma}_{QP}^{\top}\Big{)}+\frac{\partial^{2}\bm{D}}{\partial\theta_{j}\partial\theta_{k}}&\displaystyle\frac{\partial^{2}\bm{\Sigma}_{QP}}{\partial\theta_{j}\partial\theta_{k}}\\\
\displaystyle\frac{\partial^{2}\bm{\Sigma}_{QP}^{\top}}{\partial\theta_{j}\partial\theta_{k}}&\displaystyle\frac{\partial^{2}\bm{\Sigma}_{PP}}{\partial\theta_{j}\partial\theta_{k}}\end{bmatrix}\bm{\Pi},$
(50)

which follows the same permuted block diagonal plus low-rank structure as
$\widetilde{\bm{\Sigma}}$ and
$\frac{\partial\widetilde{\bm{\Sigma}}}{\partial\theta_{j}}$, now with rank at
most $4p$ in the low-rank portion of the upper left block. Thus the trace term
in each Hessian entry can be computed using equation (22) exactly as was done
for the analogous term in the gradient. This results in a linear complexity
algorithm for computing entries of the Hessian.

## References

  * [1] S. Ackerman and R. Frey.  Modis atmosphere l2 cloud mask product.  NASA MODIS adaptive processing system, Goddard Space Flight Center, USA, 2015.
  * [2] S. Ambikasaran, D. Foreman-Mackey, L. Greengard, D. W. Hogg, and M. O’Neil.  Fast direct methods for Gaussian processes.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(2):252–265, 2015.
  * [3] S. Ambikasaran, M. O’Neil, and K. R. Singh.  Fast symmetric factorization of hierarchical matrices with applications.  arXiv preprint arXiv:1405.0223, 2014.
  * [4] E. B. Anderes and M. L. Stein.  Local likelihood estimation for nonstationary random fields.  Journal of Multivariate Analysis, 102(3):506–520, 2011.
  * [5] S. Banerjee, A. E. Gelfand, A. O. Finley, and H. Sang.  Gaussian predictive process models for large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(4):825–848, 2008.
  * [6] S. Börm and J. Garcke.  Approximating Gaussian processes with $\mathcal{H}^{2}$-matrices.  In European Conference on Machine Learning, pages 42–53. Springer, 2007.
  * [7] J. Chen, H. Avron, and V. Sindhwani.  Hierarchically compositional kernels for scalable nonparametric learning.  The Journal of Machine Learning Research, 18(1):2214–2255, 2017\.
  * [8] J. Chen and M. L. Stein.  Linear-cost covariance functions for Gaussian random fields.  Journal of the American Statistical Association, pages 1–18, 2021\.
  * [9] N. Cressie and G. Johannesson.  Fixed rank kriging for very large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(1):209–226, 2008.
  * [10] J. Eidsvik, A. O. Finley, S. Banerjee, and H. Rue.  Approximate Bayesian inference for large spatial datasets using predictive process models.  Computational Statistics & Data Analysis, 56(6):1362–1380, 2012\.
  * [11] R. Furrer, M. G. Genton, and D. Nychka.  Covariance tapering for interpolation of large spatial datasets.  Journal of Computational and Graphical Statistics, 15(3):502–523, 2006.
  * [12] C. J. Geoga, M. Anitescu, and M. L. Stein.  Scalable Gaussian process computations using hierarchical matrices.  Journal of Computational and Graphical Statistics, pages 1–11, 2019\.
  * [13] C. J. Geoga, O. Marin, M. Schanen, and M. L. Stein.  Fitting Matérn smoothness parameters using automatic differentiation.  arXiv preprint arXiv:2201.00090, 2022.
  * [14] J. Guinness.  Gaussian process learning via fisher scoring of Vecchia’s approximation.  Statistics and Computing, 31(3):1–8, 2021.
  * [15] H. Huang, L. R. Blake, M. Katzfuss, and D. M. Hammerling.  Nonstationary spatial modeling of massive global satellite data.  arXiv preprint arXiv:2111.13428, 2021.
  * [16] M. F. Hutchinson.  A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines.  Communications in Statistics-Simulation and Computation, 18(3):1059–1076, 1989.
  * [17] M. Katzfuss.  A multi-resolution approximation for massive spatial datasets.  Journal of the American Statistical Association, 112(517):201–214, 2017.
  * [18] M. Katzfuss and N. Cressie.  Bayesian hierarchical spatio-temporal smoothing for very large datasets.  Environmetrics, 23(1):94–107, 2012.
  * [19] M. Katzfuss and W. Gong.  A class of multi-resolution approximations for large spatial datasets.  Statistica Sinica, 30(4):2203–2226, 2020.
  * [20] M. Katzfuss and J. Guinness.  A general framework for Vecchia approximations of gaussian processes.  Statistical Science, 36(1):124–141, 2021.
  * [21] F. Khellah, P. Fieguth, M. J. Murray, and M. Allen.  Statistical processing of large image sequences.  IEEE Transactions on Image Processing, 14(1):80–93, 2004.
  * [22] Y. Li and Y. Sun.  Efficient estimation of nonstationary spatial covariance functions with application to high-resolution climate model emulation.  Statistica Sinica, 29(3):1209–1231, 2019.
  * [23] F. Lindgren, H. Rue, and J. Lindström.  An explicit link between Gaussian fields and Gaussian markov random fields: The stochastic partial differential equation approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(4):423–498, 2011.
  * [24] A. Litvinenko, Y. Sun, M. G. Genton, and D. E. Keyes.  Likelihood approximation with hierarchical matrices for large spatial datasets.  Computational Statistics & Data Analysis, 137:115–132, 2019.
  * [25] E. Maturi, A. Harris, J. Mittaz, J. Sapper, G. Wick, X. Zhu, P. Dash, and P. Koner.  A new high-resolution sea surface temperature blended analysis.  Bulletin of the American Meteorological Society, 98(5):1015–1026, 2017.
  * [26] V. Minden, A. Damle, K. L. Ho, and L. Ying.  Fast spatial Gaussian process maximum likelihood estimation via skeletonization factorizations.  Multiscale Modeling & Simulation, 15(4):1584–1611, 2017.
  * [27] C. J. Paciorek and M. J. Schervish.  Spatial modelling using a new class of nonstationary covariance functions.  Environmetrics: The official journal of the International Environmetrics Society, 17(5):483–506, 2006.
  * [28] M. D. Risser and C. A. Calder.  Local likelihood estimation for covariance functions with spatially-varying parameters: the convospat package for r.  arXiv preprint arXiv:1507.08613, 2015.
  * [29] H. Sang and J. Z. Huang.  A full scale approximation of covariance functions for large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(1):111–132, 2012.
  * [30] H. Sang, M. Jun, and J. Z. Huang.  Covariance approximation for large multivariate spatial data sets with an application to multiple climate model errors.  The Annals of Applied Statistics, pages 2519–2548, 2011.
  * [31] E. Snelson and Z. Ghahramani.  Local and global sparse Gaussian process approximations.  In Artificial Intelligence and Statistics, pages 524–531, 2007\.
  * [32] A. Solin and S. Särkkä.  Hilbert space methods for reduced-rank Gaussian process regression.  Statistics and Computing, 30(2):419–446, 2020.
  * [33] M. L. Stein.  Interpolation of Spatial Data: Some Theory for Kriging.  Springer, 1999.
  * [34] M. L. Stein.  Limitations on low rank approximations for covariance matrices of spatial data.  Spatial Statistics, 8:1–19, 2014.
  * [35] M. L. Stein, J. Chen, and M. Anitescu.  Stochastic approximation of score functions for Gaussian processes.  The Annals of Applied Statistics, pages 1162–1191, 2013.
  * [36] M. L. Stein, Z. Chi, and L. J. Welty.  Approximating likelihoods for large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(2):275–296, 2004.
  * [37] A.V. Vecchia.  Estimation and model identification for continuous spatial processes.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 50(2):297–312, 1988.
  * [38] S. Wright and J. Nocedal.  Numerical Optimization.  Springer Science, 35(67-68):7, 1999.
  * [39] H. Zhang.  Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics.  Journal of the American Statistical Association, 99(465):250–261, 2004.

Government License: The submitted manuscript has been created by UChicago
Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a
US Department of Energy Office of Science laboratory, is operated under
Contract No. DE-AC02-06CH11357. The US Government retains for itself, and
others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide
license in said article to reproduce, prepare derivative works, distribute
copies to the public, and perform publicly and display publicly, by or on
behalf of the Government. The Department of Energy will provide public access
to these results of federally sponsored research in accordance with the DOE
Public Access Plan. http://energy.gov/downloads/doe-public-access-plan.