Kriging for large datasets via penalized neighbor selection

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.

💡 Research Summary

The paper addresses the prohibitive O(N³) computational cost of classical kriging when dealing with large spatial datasets. While local kriging based on a fixed number K of nearest neighbors reduces the cost, the choice of K is typically ad‑hoc and ignores the underlying spatial correlation and redundancy of information. The authors propose a penalized kriging framework that embeds an ℓ₁ (LASSO‑type) penalty directly into the kriging equations, thereby turning neighbor selection into a data‑driven sparsity problem.

Key methodological steps:

1. The observation locations are ordered by Euclidean distance to the prediction site. The first p nearest observations are forced to remain in the model by exploiting the unbiasedness constraint (Xᵀλ = x₀) to express these p weights as linear functions of the remaining N‑p weights. The ℓ₁ penalty is then applied only to the latter group, guaranteeing that the closest p points always receive non‑zero weights while more distant points can be shrunk to zero.
2. The resulting optimization problem is
    min λ λᵀΣλ − 2λᵀc₀ + σ₀² + η‖λ_{‑p}‖₁ subject to Xᵀλ = x₀,
   where η≥0 controls sparsity. The authors solve this convex problem with a combination of coordinate‑descent and ADMM, ensuring existence and uniqueness for any η.
3. An adaptive LASSO extension is introduced: after an initial solution λ̂, penalty weights w_j = 1/|λ̂_j|^γ (γ>0) are computed, and the penalty term becomes η∑w_j|λ_j|. This weighting improves variable‑selection consistency and better respects heterogeneous spatial correlation.
4. Selecting η is critical. Instead of costly cross‑validation, the authors develop a novel criterion based on the effective sample size (ESS). ESS quantifies the amount of independent information remaining after accounting for spatial correlation and the ℓ₁ shrinkage; as η grows, ESS declines. The optimal η̂ minimizes a loss that balances prediction mean‑squared error against ESS, automatically trading off accuracy and redundancy.

Experimental evaluation:

• Simulated data with various covariance models (exponential, Gaussian, Matérn) and differing smoothness/variance levels show that the method automatically selects fewer neighbors for smooth fields and more for highly variable regions. Prediction mean‑squared error matches that of full‑sample kriging, while computational time is reduced by an order of magnitude compared with K‑nearest‑neighbor kriging and by two orders of magnitude versus global kriging.
• Real‑world case studies: (a) the Jura heavy‑metal dataset (≈300 points) where the method uses K≈30–40, achieving RMSE within 0.02 of the global kriging benchmark; (b) the COBE sea‑surface‑temperature dataset (≈43 000 points) where the algorithm adapts K locally (≈120 near the equator, ≈30 in polar regions), preserving predictive skill while completing the analysis in roughly two hours—a substantial speed‑up over traditional approaches.

The paper’s contributions are threefold: (i) a principled ℓ₁‑penalized kriging formulation that yields automatic, spatially adaptive neighbor selection; (ii) an ESS‑based, computationally cheap tuning‑parameter rule that respects spatial dependence; (iii) an adaptive LASSO weighting scheme that enhances selection stability. Limitations include reliance on Gaussianity and linear mean structures, and the need to pre‑specify the small set p of forced neighbors. Future work is suggested on extending to non‑linear mean models, multi‑scale neighbor sets, and high‑performance parallel implementations.

Overall, the study provides a practical, theoretically grounded tool for scalable spatial prediction, bridging the gap between computational feasibility and statistical optimality in the era of massive geospatial data.