Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation

TL;DR

This paper introduces an automatic differentiation-based implementation of the Matérn covariance function's derivatives with respect to its smoothness parameter, enabling faster and more accurate estimation in Gaussian process models.

Contribution

We develop a new AD-compatible implementation of the modified Bessel function for the Matérn covariance, improving derivative computation for parameter estimation.

Findings

AD derivatives are faster and more accurate than finite differences.

Our method enables reliable second-order maximum likelihood estimation.

Hessian matrices built with AD derivatives outperform finite difference approximations.

Abstract

The Mat\'ern 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 Mat\'ern covariance function with respect to $\nu$ require derivatives of the modified second-kind Bessel function $\mathcal{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_\nu \mathcal{K}_\nu$. In this work, we introduce a new implementation of $\mathcal{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.