Nearest-Neighbor Derivative Process (NNDP)

• NNDP is a scalable Gaussian-process framework that jointly infers a spatial process and its derivatives, lowering computational complexity from O(n^3) to O(n).
• It constructs an augmented state by bundling function values and derivatives through a nearest-neighbor factorization, ensuring coherent posterior inference without relying on finite difference approximations.
• Empirical results show accurate gradient estimation and improved efficiency in applications like spatial transcriptomics and air temperature mapping.

Searching arXiv for the cited NNDP paper and closely related nearest-neighbor GP context. The Nearest-Neighbor Derivative Process (NNDP) is a scalable Gaussian-process framework for joint inference on a spatial process $w(s)$ and its derivatives, particularly the gradient field $D(s)=\nabla w(s)$ and, when smoothness permits, second derivatives such as the Hessian. It is designed for settings in which full Gaussian-process inference on $[w,\nabla w]$ is computationally prohibitive because dense covariance matrices induce cubic complexity in the number of locations. NNDP replaces the dense joint density by a product of sparse nearest-neighbor conditionals on an augmented state that bundles function values and derivatives at each site, thereby reducing the computational time complexity from $O(n^3)$ to $O(n)$ when the neighborhood size and spatial dimension are treated as fixed (Chen et al., 2 Sep 2025).

1. Problem setting and motivation

Many tasks in spatial analysis require not only interpolation of a latent field $w(s)$ but also inference on spatial rates of change. The target object is often the gradient field $D(s)=\nabla w(s)$, and in some applications second derivatives are also relevant for curvature or Hessian analysis. A full Gaussian process for $[w,\nabla w]$ is statistically natural because any finite collection of function values and derivatives is multivariate Gaussian when the kernel is sufficiently differentiable. The computational obstacle is that the corresponding covariance matrices are dense and scale cubically with the number of spatial locations (Chen et al., 2 Sep 2025).

NNDP addresses this obstacle by extending Nearest-Neighbor Gaussian Process ideas to derivative inference. The central construction augments the state at each location $s_i$ as

$$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$

and then approximates the joint density of the stacked vector $D(s)=\nabla w(s)$0 by a directed acyclic nearest-neighbor factorization. This preserves a valid stochastic-process interpretation for the joint field of function values and derivatives while avoiding dense $D(s)=\nabla w(s)$1 algebra.

The method is positioned against a straightforward plug-in estimator that first fits a scalable Gaussian process for $D(s)=\nabla w(s)$2 and then approximates $D(s)=\nabla w(s)$3 by finite differences of predicted means. That approach is sensitive to the step size $D(s)=\nabla w(s)$4, to grid spacing, and to sampling irregularities; moreover, the approximate gradient may not be a valid Gaussian process, which compromises inference and uncertainty quantification. NNDP instead directly models $D(s)=\nabla w(s)$5 as a nearest-neighbor Gaussian process and therefore supplies coherent posterior inference for both fields without tuning $D(s)=\nabla w(s)$6 (Chen et al., 2 Sep 2025).

2. Gaussian-process and derivative-process foundation

NNDP begins from a parent Gaussian process

$D(s)=\nabla w(s)$7

The derivative process is defined in the mean-square sense, so kernel smoothness is decisive. First derivatives require $D(s)=\nabla w(s)$8 to be once differentiable in each argument, typically ensured for Matérn with $D(s)=\nabla w(s)$9 or for RBF. Second derivatives require $[w,\nabla w]$0 to be twice differentiable, for example Matérn with $[w,\nabla w]$1. Under these conditions, the joint Gaussian structure for function values and gradient components is determined by differentiating the kernel (Chen et al., 2 Sep 2025):

$[w,\nabla w]$2

$[w,\nabla w]$3

$[w,\nabla w]$4

$[w,\nabla w]$5

If second derivatives are modeled, analogous formulas hold for Hessian entries $[w,\nabla w]$6 and their cross-covariances.

For isotropic kernels with $[w,\nabla w]$7, $[w,\nabla w]$8, and $[w,\nabla w]$9, the derivative formulas take a compact vector form:

$O(n^3)$0

and

$O(n^3)$1

These identities furnish the $O(n^3)$2 cross-covariance matrix for gradient components.

Two kernel families are explicit in the NNDP formulation. For the RBF kernel

$O(n^3)$3

with $O(n^3)$4 and $O(n^3)$5,

$O(n^3)$6

$O(n^3)$7

RBF is infinitely differentiable, so all derivative processes exist. For the Matérn kernel with smoothness $O(n^3)$8,

$O(n^3)$9

the radial derivatives are

$O(n)$0

$O(n)$1

Plugging these into the isotropic formulas yields explicit cross-covariances for gradients and Hessians. The smoothness thresholds are explicit: first derivatives exist in mean square when $O(n)$2, and second derivatives exist when $O(n)$3 (Chen et al., 2 Sep 2025).

3. Nearest-neighbor factorization and process construction

The defining step in NNDP is the nearest-neighbor factorization of the augmented state. After choosing a directed acyclic ordering of locations $O(n)$4, one selects for each site $O(n)$5 a neighbor set $O(n)$6 consisting of the indices of the $O(n)$7 nearest neighbors among $O(n)$8. NNDP then defines the approximation

$O(n)$9

where each conditional factor is multivariate normal with mean and covariance derived from the corresponding blocks of the full covariance matrix of the stacked augmented state (Chen et al., 2 Sep 2025):

$w(s)$0

$w(s)$1

Here $w(s)$2 collects the mean components

$w(s)$3

$w(s)$4 is the on-site covariance block of $w(s)$5 and its derivatives at $w(s)$6, $w(s)$7 comprises cross-covariances between $w(s)$8 at $w(s)$9 and all stacked $D(s)=\nabla w(s)$0 at neighbors $D(s)=\nabla w(s)$1, and $D(s)=\nabla w(s)$2 is the block covariance over neighbors. These blocks are assembled from $D(s)=\nabla w(s)$3, its first derivatives, and its mixed second derivatives evaluated at the relevant pairs of locations.

This construction extends NNGP from $D(s)=\nabla w(s)$4 to $D(s)=\nabla w(s)$5 by retaining the same sparse conditional structure while replacing scalar or single-field covariance entries with block covariances formed by kernel differentiation. The result is a scalable joint model for the parent field and its derivative field. The recommended neighborhood size is fixed $D(s)=\nabla w(s)$6, with values in $D(s)=\nabla w(s)$7–$D(s)=\nabla w(s)$8 typically balancing accuracy and speed. This suggests that the practical operating regime is one in which local conditioning captures the dominant dependence structure while keeping local linear algebra modest (Chen et al., 2 Sep 2025).

4. Likelihood, inference, prediction, and computational scaling

Under the NNDP approximation, the log-likelihood of the augmented state is

$D(s)=\nabla w(s)$9

In practice, derivatives need not be observed. The likelihood may instead be written on observed $[w,\nabla w]$0 and latent $[w,\nabla w]$1, with NNDP supplying the conditional structure for either maximum-likelihood or Bayesian inference. The observation model for noisy function evaluations is

$[w,\nabla w]$2

and if derivative observations are available,

$[w,\nabla w]$3

independently. The joint model remains Gaussian after adding the appropriate block-diagonal noise components (Chen et al., 2 Sep 2025).

Parameter inference can proceed by maximum likelihood or by Bayesian computation. Maximum likelihood uses optimization of $[w,\nabla w]$4 or the corresponding latent-variable objective via quasi-Newton methods such as L-BFGS. In the Bayesian formulation, priors are placed on $[w,\nabla w]$5 and $[w,\nabla w]$6; the specified examples are IG priors on $[w,\nabla w]$7, uniform on $[w,\nabla w]$8, fixed $[w,\nabla w]$9 such as $s_i$0, and IG on $s_i$1, with posterior sampling through Gibbs/Metropolis within the NNDP conditional structure.

The linear-time scaling is obtained because each site only conditions on $s_i$2 neighbors. With fixed $s_i$3 and $s_i$4, the per-site Cholesky factorization of $s_i$5 costs $s_i$6, and the remaining local solves and block operations cost $s_i$7 per site. Consequently, the total time is

$s_i$8

which is linear in $s_i$9 when $$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$0 and $$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$1 are small constants, and the memory cost is

$$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$2

without forming dense $$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$3 matrices. Numerical stability is handled through Cholesky decompositions of $$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$4 with jitter, such as adding $$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$5 with $$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$6 if needed, together with kernel parameter constraints, reasonable nugget $$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$7, and coordinate standardization (Chen et al., 2 Sep 2025).

Prediction at a new location $$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$8 follows the same local-conditioning principle. Defining

$$z_i=\left[w(s_i),\frac{\partial w}{\partial s_1}(s_i),\ldots,\frac{\partial w}{\partial s_d}(s_i)\right]^T \in \mathbb{R}^{1+d},$$9

and selecting a neighbor set $D(s)=\nabla w(s)$00 of size $D(s)=\nabla w(s)$01 from the training locations, the predictive distribution is

$D(s)=\nabla w(s)$02

$D(s)=\nabla w(s)$03

The requisite blocks use exactly the same kernel derivatives as the training-time construction.

5. Theoretical properties and methodological contrasts

NNDP is presented as a valid stochastic process for

$D(s)=\nabla w(s)$04

over $D(s)=\nabla w(s)$05, where $D(s)=\nabla w(s)$06 is a set of Lebesgue measure zero associated with neighbor-set discontinuities, such as ties in nearest-neighbor selection. The validity claim follows the NNGP construction together with derivative-process existence, and yields a proper joint Gaussian density over any finite collection. Smoothness is inherited in a corresponding sense: if the parent Gaussian process $D(s)=\nabla w(s)$07 is mean-square differentiable of order $D(s)=\nabla w(s)$08 or $D(s)=\nabla w(s)$09, then the NNGP approximation $D(s)=\nabla w(s)$10 is also mean-square differentiable of order $D(s)=\nabla w(s)$11 or $D(s)=\nabla w(s)$12 outside a measure-zero set, and the derivative and cross-covariance functions for NNDP are continuous away from $D(s)=\nabla w(s)$13. For Matérn kernels, mean-square differentiability of order $D(s)=\nabla w(s)$14 requires $D(s)=\nabla w(s)$15 (Chen et al., 2 Sep 2025).

As a Vecchia-type approximation, NNDP is also described as consistent in the limit of increasing neighborhood size. Increasing $D(s)=\nabla w(s)$16 shrinks approximation error and recovers the full Gaussian-process derivative process as $D(s)=\nabla w(s)$17 with appropriate orderings. Error bounds are inherited from Vecchia approximations for Gaussian-process likelihoods, while derivative-process blocks improve because they are constructed by differentiating the same parent kernel. Under correct kernel specification, the NNDP conditional means deliver approximately unbiased gradient estimates, and uncertainty quantification is valid in the NNDP model.

The principal methodological contrast is with plug-in finite differences. In that approach, a directional derivative is approximated as

$D(s)=\nabla w(s)$18

where $D(s)=\nabla w(s)$19 is a Gaussian-process predictor. The stated deficiencies are the need to tune $D(s)=\nabla w(s)$20, sensitivity to grid spacing and sampling irregularity, numerical errors for small $D(s)=\nabla w(s)$21, bias for large $D(s)=\nabla w(s)$22, absence of a joint Gaussian process for gradients, and increasing instability for higher-order derivatives. NNDP avoids these issues by modeling $D(s)=\nabla w(s)$23 jointly with valid cross-covariances and by enabling principled likelihood-based or posterior-based inference. The comparison to other scalable Gaussian-process approaches is narrower but equally important: NNGP, low-rank, and tapering methods target $D(s)=\nabla w(s)$24 only and do not natively supply derivative-process covariances. NNDP differs by augmenting the state with derivatives and building blockwise covariances using kernel derivatives, thereby preserving valid uncertainty for spatial rates of change (Chen et al., 2 Sep 2025).

6. Empirical behavior, applications, and limitations

The empirical evaluation reported for NNDP comprises simulations and real-data analyses. In simulations, two setups are emphasized. Pattern 1 is a two-dimensional sinusoidal surface,

$D(s)=\nabla w(s)$25

with $D(s)=\nabla w(s)$26, regular grids, and Matérn $D(s)=\nabla w(s)$27, compared across exact GP, NNDP, and plug-in finite differences. Pattern 2 is a one-dimensional highly oscillatory function,

$D(s)=\nabla w(s)$28

used to assess sensitivity to grid spacing and irregular sampling. The evaluation metrics are correlation with truth, mean squared error for gradients, and runtime (Chen et al., 2 Sep 2025).

The reported quantitative results are specific. For Pattern 1 on a regular grid, NNDP correlations are approximately $D(s)=\nabla w(s)$29–$D(s)=\nabla w(s)$30 with low MSEs; finite differences are substantially worse and sensitive to $D(s)=\nabla w(s)$31, and exact GP matches NNDP when feasible for small $D(s)=\nabla w(s)$32. For Pattern 2, NNDP consistently outperforms finite differences, while the optimal finite-difference scale depends on grid size and irregular sampling. Runtime results are used to illustrate the asymptotic scaling: for $D(s)=\nabla w(s)$33 up to $D(s)=\nabla w(s)$34 with $D(s)=\nabla w(s)$35, NNDP times scale linearly, with an example of approximately $D(s)=\nabla w(s)$36 seconds, whereas exact GP exhibits cubic behavior, with an example of $D(s)=\nabla w(s)$37 seconds for $D(s)=\nabla w(s)$38. Finite differences are described as slower than NNDP in comparable settings because they require separate predictions at perturbed locations $D(s)=\nabla w(s)$39.

Two real-data applications illustrate the intended use cases. In VisiumHD spatial transcriptomics on mouse brain data, involving $D(s)=\nabla w(s)$40 genes and $D(s)=\nabla w(s)$41 spots, NNDP with Matérn $D(s)=\nabla w(s)$42 produced high-resolution gradient maps aligned with anatomical boundaries, including the hippocampal formation, CA fields, and dentate gyrus; representative gene gradients were computed in approximately $D(s)=\nabla w(s)$43 seconds. In NARR air temperature data over the continental United States with $D(s)=\nabla w(s)$44 grid cells, NNDP identified large temperature gradients along topographic features such as the Sierra Nevada, Cascade Range, and Rockies (Chen et al., 2 Sep 2025).

Practical guidance in the NNDP description focuses on smoothness, neighborhood size, and model assumptions. Values of $D(s)=\nabla w(s)$45 in $D(s)=\nabla w(s)$46–$D(s)=\nabla w(s)$47 typically balance accuracy and speed, while larger $D(s)=\nabla w(s)$48 improve approximation at the cost of increased local cubic work in $D(s)=\nabla w(s)$49. Matérn kernels are recommended for controlled smoothness, with $D(s)=\nabla w(s)$50 for first derivatives in practice and $D(s)=\nabla w(s)$51 for second derivatives; fixing $D(s)=\nabla w(s)$52 to values such as $D(s)=\nabla w(s)$53 or $D(s)=\nabla w(s)$54 is described as stabilizing inference. Anisotropy can be modeled by separate length scales per coordinate. The limitations are equally explicit: derivative estimation near domain boundaries is harder because there are fewer neighbors and more extrapolation, nonstationary gradients require kernel extensions because NNDP assumes stationarity via $D(s)=\nabla w(s)$55, and computational cost scales with $D(s)=\nabla w(s)$56 per site, making the method practical for $D(s)=\nabla w(s)$57–$D(s)=\nabla w(s)$58 rather than very high-dimensional spatial domains. This suggests that NNDP is best viewed as a scalable derivative-process model for low-dimensional spatial settings with smooth kernels and locality-compatible dependence structure (Chen et al., 2 Sep 2025).