  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.
