Numerical evaluation of the Matérn covariance¶

The continuous-regularity functions gpmp.kernel.matern_kernel() and gpmp.kernel.matern_covariance() evaluate the Matérn correlation for a real regularity parameter \(\nu>0\). Their numerical implementation depends on the active backend. Use the default large_nu_method="auto" for model construction and parameter selection.

Primitive convention¶

The backend primitive evaluates

\[\rho_\nu(q) = \frac{2^{1-\nu}}{\Gamma(\nu)} z^\nu K_\nu(z), \qquad z=\sqrt{2\nu q}, \qquad \rho_\nu(0)=1.\]

Here \(q\) is a squared distance in the standard Matérn scaling. The public matern_kernel(nu, h) function passes \(q=2h^2\) to this primitive. Consequently its Bessel argument is \(2\sqrt{\nu}h\), which matches maternp_kernel at \(\nu=p+1/2\).

Backend policy¶

NumPy uses scipy.special.kve in the log domain for scalar \(\nu<70\). For scalar \(\nu\geq70\), "auto" evaluates the full correlation array with the uniform large-order expansion. This global switch avoids covariance matrices that mix scaled-Bessel values with individually repaired Gaussian-limit entries when SciPy becomes non-finite. Explicit "quadrature" forces the SciPy path, while "uniform_asymptotic" forces the expansion in its supported range. The nodes and t_max arguments are accepted for API consistency but do not change the NumPy computation.

Torch uses the following dispatch:

Torch evaluation policy¶
Method	Condition	Evaluation
`"auto"`	fixed scalar \(\nu\in\{1/2,3/2,5/2\}\) not differentiated with respect to \(\nu\)	exact closed form
`"auto"`	CPU float64, scalar \(0<\nu<70\), no gradients, sufficient input size under the active dispatch profile	value-only reduced-order recurrence
`"auto"`	CPU float64, scalar \(1/2\leq\nu<70\), gradients required, and the recurrence stability and calibrated-size conditions satisfied	reduced-order recurrence with an explicit first-order backward
`"auto"`	scalar \(\nu\geq70\)	uniform large-order expansion
`"auto"`	input below its calibrated threshold, calibration thread mismatch, or any other fallback case	log-space Gauss–Legendre quadrature
`"quadrature"`	any supported \(\nu\)	log-space Gauss–Legendre quadrature
`"uniform_asymptotic"`	scalar \(\nu\geq70\)	uniform large-order expansion

An explicit "uniform_asymptotic" request with nonscalar \(\nu\) or \(\nu<70\) raises ValueError. This guard prevents use of the expansion outside its validated domain.

Torch recurrence and first derivatives¶

For an eligible scalar order, the implementation writes

\[n=\left\lfloor\nu+\frac12\right\rfloor, \qquad \mu=\nu-n\in[-1/2,1/2).\]

It first computes the exponentially scaled values \(e^zK_\mu(z)\) and \(e^zK_{\mu+1}(z)\). A Temme series is used for \(z\leq2\), a Miller construction for \(2<z\leq17\), and a large-argument expansion for \(z>17\). The target order is then recovered with

\[K_{a+1}(z)=K_{a-1}(z)+\frac{2a}{z}K_a(z).\]

The lifting recurrence is evaluated in the log domain. Its order derivatives are propagated together with the values, and the resulting derivatives with respect to squared distance and \(\nu\) are saved by a custom first-order backward.

The recurrence has separate value-only and first-derivative implementations. The value-only path supports scalar \(0<\nu<70\). The derivative path supports scalar \(1/2\leq\nu<70\). Exact half-integer cases use their closed formulas when \(\nu\) is not differentiated. In a width \(10^{-7}\) around a reduced half-integer, excluding the exact point, the differentiated Miller construction is ill-conditioned, so "auto" retains quadrature. Gradient evaluations below \(1/2\) also retain quadrature because their stable recurrence requires a second reduced-order evaluation and is slower in the tested workloads. Float32, GPU arrays, vector orders, and higher derivatives retain quadrature.

Calibrated size dispatch¶

Recurrence has a larger fixed cost than quadrature for short vectors. Torch therefore compares q.numel() with the active process-wide dispatch profile. For a self-covariance matrix, GPmp evaluates only the strict upper triangle, so the primitive receives \(n(n-1)/2\) values rather than \(n^2\) values.

The shipped profile was calibrated with CPU float64 arithmetic and one Torch thread. Its global mode uses a threshold of 8192 values for both value and gradient evaluations. The optional "nu_buckets" mode uses:

Shipped one-thread bucket thresholds¶
Smoothness interval	Value threshold	Gradient threshold
\(0<\nu<2\)	8192	8192, with quadrature below \(\nu=1/2\)
\(2\leq\nu<10\)	8192	8192
\(10\leq\nu<30\)	8192	4096
\(30\leq\nu<50\)	4096	4096
\(50\leq\nu<70\)	4096	4096

The profile records the Torch thread count used for calibration. If the active thread count differs, "auto" uses quadrature. This prevents thresholds measured under one parallel configuration from being reused under another.

Local calibration is explicit:

import torch
from gpmp.config import get_config
from gpmp.num import matern_torch

torch.set_num_threads(4)
result = matern_torch.calibrate_matern_dispatch()

dispatch = get_config().matern_torch_dispatch
dispatch.mode = "nu_buckets"  # or "global"

The function benchmarks paired recurrence and quadrature calls, verifies their values and gradients, updates both profiles for the current process, and returns the complete measurements. Threshold selection uses the lower quartile of paired speedups and, by default, requires a 50 percent gain at every larger measured size. A None threshold means that recurrence did not reach the margin in the tested range. Calibration does not run during import or model fitting and does not write a machine-specific cache.

Torch quadrature and derivatives¶

For \(z>0\), the scaled modified Bessel function has the integral representation

\[e^z K_\nu(z) = \int_0^\infty \exp\{-z(\cosh t-1)\}\cosh(\nu t)\,\mathrm{d}t.\]

Torch applies a Gauss–Legendre rule on a truncated interval. nodes=160 and t_max=20 define the base quadrature density and interval. The implementation locates the rightmost estimated integrand saddle, keeps a guarded tail beyond it, and shortens or extends the interval accordingly. The node count is scaled with the interval so that the base quadrature density is preserved. It can therefore be smaller or larger than nodes.

The quadrature custom backward does not differentiate through the full quadrature graph. For scalar \(\nu\), one normalized quadrature table supplies the Matérn value, the adjacent-order ratio used for the distance derivative, and the order derivative. The resulting derivatives with respect to squared distance and \(\nu\) are saved for the backward pass. This avoids repeating the costly quadrature normalization for each derivative.

Large-order expansion¶

For scalar \(\nu\geq70\), both backends evaluate the exponentially scaled Bessel function with a ten-term uniform large-order expansion. Its Debye-polynomial recurrence uses the coefficient organization from the SLATEC DASYIK routine. NumPy evaluates the correlation values. Torch also propagates derivatives with respect to the Bessel argument and order through the same recurrence [3, 14, 15].

The switch at \(\nu=70\) is a numerical dispatch boundary, not a change in the covariance definition. Tests compare values from both methods around the boundary and also compare Torch derivatives. large_nu_method="quadrature" remains available as the backend reference computation near the boundary.

Diagonal and covariance checks¶

The primitive returns exactly one at \(q=0\). Its numerical diagonal convention sets the derivatives with respect to \(q\) and \(\nu\) to zero. Self-covariance functions also set diagonal squared distances to zero before evaluating the correlation. They evaluate the strict upper triangle once and reflect it to form the symmetric covariance matrix.

The test suite checks recurrence values and derivatives against quadrature, regional-boundary gradients, exact and neighboring half-integer orders, agreement around \(\nu=70\), covariance matrices with near-duplicate points, likelihood values, and restricted-likelihood gradients. These checks control numerical error. They do not replace the positive-semidefiniteness theory of the exact Matérn covariance.

API controls¶

large_nu_method: Use "auto" in model code. Use "quadrature" to force SciPy with the NumPy backend or quadrature with the Torch backend. Use "uniform_asymptotic" only for a deliberate scalar large-order evaluation with \(\nu\geq70\).
nodes: Base Gauss–Legendre node count for Torch quadrature. It has no effect with NumPy.
t_max: Base upper integration limit for Torch quadrature. It has no effect with NumPy. Torch can shorten it when the guarded quadrature tail is covered or increase it for difficult small-argument evaluations.
gpmp.num.matern_torch.calibrate_matern_dispatch: Benchmark and install process-local Torch recurrence thresholds. This Torch-specific function is intentionally not re-exported by gpmp.num.

See gpmp.kernel covariance functions for the public covariance functions and Higher-dimensional interpolation with Matérn regularity selection for parameter selection with \(\nu\) in the covariance-parameter vector.