Fitting Matern smoothness parameters using automatic differentiation

The Matern covariance function is ubiquitous in the application of Gaussian processes to spatial statistics and beyond. Perhaps the most important reason for this is that the smoothness parameter nu gives complete control over the mean-square differentiability of the process, which has significant implications for the behavior of estimated quantities such as interpolants and forecasts. Unfortunately, derivatives of the Matern covariance function with respect to nu require derivatives of the modified second-kind Bessel function K-nu with respect to nu. While closed form expressions of these derivatives do exist, they are prohibitively difficult and expensive to compute. For this reason, many software packages require fixing nu as opposed to estimating it, and all existing software packages that attempt to offer the functionality of estimating nu use finite difference estimates for partial derivative K-nu(nu). In this work, we introduce a new implementation of K-nu that has been designed to provide derivatives via automatic differentiation (AD), and whose resulting derivatives are significantly faster and more accurate than those computed using finite differences. We provide comprehensive testing for both speed and accuracy and show that our AD solution can be used to build accurate Hessian matrices for second-order maximum likelihood estimation in settings where Hessians built with finite difference approximations completely fail.