  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.
