Kriging for large datasets via penalized neighbor selection

TL;DR

This work tackles the cubic computational burden of classical kriging by embedding $\ell_1$ penalties into kriging equations to automatically select informative neighbors. The authors introduce penalized kriging (and adaptive LASSO) with an effective sample size–based criterion to tune the penalty, balancing sparsity and prediction variance without cross-validation. The method adapts the neighborhood to the underlying spatial correlation, yielding accuracy close to global kriging while dramatically reducing computation, as demonstrated on simulated data and real-scale Jura and COBE SST datasets. The approach provides a scalable, principled alternative to fixed $K$-NN or covariance tapering, with potential extensions to non-Gaussian, non-stationary, and spatio-temporal settings.

Abstract

Kriging is a fundamental tool for spatial prediction, but its computational complexity of $O(N^3)$ becomes prohibitive for large datasets. While local kriging using $K$-nearest neighbors addresses this issue, the selection of $K$ typically relies on ad-hoc criteria that fail to account for spatial correlation structure. We propose a penalized kriging framework that incorporates LASSO-type penalties directly into the kriging equations to achieve automatic, data-driven neighbor selection. We further extend this to adaptive LASSO, using data-driven penalty weights that account for the spatial correlation structure. Our method determines which observations contribute non-zero weights through $\ell_1$ regularization, with the penalty parameter selected via a novel criterion based on effective sample size that balances prediction accuracy against information redundancy. Numerical experiments demonstrate that penalized kriging automatically adapts neighborhood structure to the underlying spatial correlation, selecting fewer neighbors for smoother processes and more for highly variable fields, while maintaining prediction accuracy comparable to global kriging at substantially reduced computational cost.

Figures (16)

• Figure 1: Left: surface $\lambda_{1} = 1 - \lambda_{2} - \lambda_{3}$ obtained from Problem \ref{['eq:proposed_lasso_problem_mod']} under the assumption that the mean is unknown (ordinary kriging setup). The green and the blue curves are the intersection of the LASSO restrictions $|\bm{\lambda}_{-1}| = 0.291$ (OK with two locations) and $|\bm{\lambda}_{-1}| = 0.615$ (exact OK solution) with the restriction plane respectively. The red curve is the solution path $\bm{\lambda}(\eta)$, that connects the solution $\bm{\lambda}(M) = [1,0,0]^{\top}$ and the ordinary kriging solution $\widehat{\bm{\lambda}}(0) = [0.683, 0.466, -0.149]^{\top}$. Right: Trace plot created with the solution path $\widehat{\bm{\lambda}}(\eta)$, showing how the method assigns zero weight to the farthest location.
• Figure 2: Behavior of the tuning parameter selection criterion for exponential covariance models with practical range (PR) of $0.05$ (left), $0.15$ (center), and $0.2$ (right). The variance ratio $v_{\eta}$ (red line) increases monotonically with $\eta$, while the sparsity measure $s_{\eta}$ (black line) shows non-monotonic behavior. The harmonic mean criterion (blue line) achieves its maximum at intermediate values, with optimal $\eta^{\star}$ increasing for smoother processes (larger PR).
• Figure 3: Locations used in the simulation experiment. The gray points are used as the training set, while the red points are the prediction locations used to explore the neighborhood structure selected with penalized kriging. The labels stand for farthest neighbor (FN), average neighbor (AN), densest neighbor (DN), side neighbor (SN), and corner neighbor (CN).
• Figure 4: Number of neighbors selected by penalized kriging as a function of the practical range for three covariance models. Each curve represents one of the five prediction locations (FN: black line; AN: red lines; DN: green line; CN: blue line; SN: magenta line). The number of selected neighbors decreases as the practical range increases for all models, with the rate of decrease depending on the process smoothness and the geometry of the prediction locations.
• Figure 5: Relative variance increase (%) compared to global kriging (GK) using all 500 observations. Red lines correspond to penalized kriging (PK), green lines to local kriging (LK) using the same number of $K$-nearest neighbors as PK, and the black horizontal line indicates the global kriging baseline (0%). Top row shows the effect of the covariance model with prediction location fixed at AN, while the bottom row shows the effect of prediction location under exponential covariance.