Kriging and Spatial Design Accelerated by Orders of Magnitude: Combining Low-Rank Covariance Approximations with FFT-Techniques

Computational power poses heavy limitations to the achievable problem size for Kriging. In separate research lines, Kriging algorithms based on FFT, the separability of certain covariance functions, and low-rank representations of covariance functions have been investigated, all three leading to drastic speedup factors. The current study combines these ideas, and so combines the individual speedup factors of all ideas. This way, we reduce the mathematics behind Kriging to a computational complexity of only $\mathcal{O}(dL^{*} \log L^{*})$, where L∗ is the number of points along the longest edge of the involved lattice of estimation points, and d is the physical dimensionality of the lattice. For separable (factorized) covariance functions, the results are exact, and nonseparable covariance functions can be approximated well through sums of separable components.Only outputting the final estimate as an explicit map causes computational costs of $\mathcal{O}(n)$, where n is the number of estimation points. In illustrative numerical test cases, we achieve speedup factors up to 108 (eight orders of magnitude), and we can treat problem sizes of up to 15 trillion and two quadrillion estimation points for Kriging and spatial design, respectively, within seconds on a contemporary desktop computer. The current study assumes second-order stationarity and simple Kriging on a regular, equispaced lattice, without working with restricted neighborhoods. Extensions to many other cases are straightforward.

[1]  Wolfgang Nowak,et al.  Efficient Computation of Linearized Cross-Covariance and Auto-Covariance Matrices of Interdependent Quantities , 2003 .

[2]  Noel A Cressie,et al.  Statistics for Spatial Data. , 1992 .

[3]  Joseph R. Kasprzyk,et al.  A new epsilon-dominance hierarchical Bayesian optimization algorithm for large multiobjective monitoring network design problems , 2008 .

[4]  C. Loan Computational Frameworks for the Fast Fourier Transform , 1992 .

[5]  C. Loan,et al.  Approximation with Kronecker Products , 1992 .

[6]  Wolfgang Nowak,et al.  Measures of Parameter Uncertainty in Geostatistical Estimation and Geostatistical Optimal Design , 2010 .

[7]  Boris N. Khoromskij,et al.  Low-rank Kronecker-product Approximation to Multi-dimensional Nonlocal Operators. Part I. Separable Approximation of Multi-variate Functions , 2005, Computing.

[8]  Ana Cortés,et al.  Parallel ordinary kriging interpolation incorporating automatic variogram fitting , 2011, Comput. Geosci..

[9]  Hermann G. Matthies,et al.  Low-rank data format for uncertainty quantification , 2010 .

[10]  J. A. Vargas-Guzmán,et al.  Sequential kriging and cokriging: Two powerful geostatistical approaches , 1999 .

[11]  Richard S. Varga,et al.  Eigenvalues of circulant matrices. , 1954 .

[12]  Timothy C. Coburn,et al.  Geostatistics for Natural Resources Evaluation , 2000, Technometrics.

[13]  Patrick Bogaert,et al.  Comparison of kriging techniques in a space-time context , 1996 .

[14]  Pejman Tahmasebi,et al.  Accelerating geostatistical simulations using graphics processing units (GPU) , 2012, Comput. Geosci..

[15]  Jürgen Pilz,et al.  Spatial sampling design and covariance-robust minimax prediction based on convex design ideas , 2010 .

[16]  H. Omre Bayesian kriging—Merging observations and qualified guesses in kriging , 1987 .

[17]  Steven G. Johnson,et al.  FFTW: an adaptive software architecture for the FFT , 1998, Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP '98 (Cat. No.98CH36181).

[18]  B. Minsker,et al.  Cost‐effective long‐term groundwater monitoring design using a genetic algorithm and global mass interpolation , 2000 .

[19]  Fred C. Schweppe,et al.  Uncertain dynamic systems , 1973 .

[20]  W. Nowak,et al.  Application of FFT-based Algorithms for Large-Scale Universal Kriging Problems , 2009 .

[21]  Ali H. Sayed,et al.  Displacement Structure: Theory and Applications , 1995, SIAM Rev..

[22]  C. R. Dietrich,et al.  A stability analysis of the geostatistical approach to aquifer transmissivity identification , 1989 .

[23]  E. G. Vomvoris,et al.  A geostatistical approach to the inverse problem in groundwater modeling (steady state) and one‐dimensional simulations , 1983 .

[24]  Gene H. Golub,et al.  Linear algebra for large scale and real-time applications , 1993 .

[25]  Peter K. Kitanidis,et al.  Efficient methods for large‐scale linear inversion using a geostatistical approach , 2012 .

[26]  A. Journel,et al.  Geostatistics for natural resources characterization , 1984 .

[27]  S. De Iaco,et al.  Towards an automatic procedure for modeling multivariate space-time data , 2012, Comput. Geosci..

[28]  S. Goreinov,et al.  A Theory of Pseudoskeleton Approximations , 1997 .

[29]  Gene H. Golub,et al.  Matrix computations (3rd ed.) , 1996 .

[30]  W. Nowak,et al.  A modified Levenberg-Marquardt algorithm for quasi-linear geostatistical inversing , 2004 .

[31]  Wolfgang Nowak,et al.  First‐order variance of travel time in nonstationary formations , 2004 .

[32]  Dale L. Zimmerman,et al.  Computationally exploitable structure of covariance matrices and generalized convariance matrices in spatial models , 1989 .

[33]  Peter K. Kitanidis,et al.  Introduction to geostatistics , 1997 .

[34]  Hans Wackernagel,et al.  Multivariate Geostatistics: An Introduction with Applications , 1996 .

[35]  P. K. Kitanidis Introduction to Geostatistics: Generation of realizations , 1997 .

[36]  Boris N. Khoromskij,et al.  Low-rank Kronecker-product Approximation to Multi-dimensional Nonlocal Operators. Part II. HKT Representation of Certain Operators , 2005, Computing.

[37]  P. A. Finkea,et al.  Mapping groundwater dynamics using multiple sources of exhaustive high resolution data , 2004 .

[38]  Alexander Litvinenko,et al.  Data Sparse Computation of the Karhunen-Loeve Expansion , 2008 .

[39]  Stephen M. Wesson,et al.  Radar rainfall image repair techniques , 2004 .

[40]  M. Genton Separable approximations of space‐time covariance matrices , 2007 .

[41]  W. Hackbusch Tensor Spaces and Numerical Tensor Calculus , 2012, Springer Series in Computational Mathematics.

[42]  Ying Sun,et al.  Geostatistics for Large Datasets , 2012 .

[43]  Raymond H. Chan,et al.  Conjugate Gradient Methods for Toeplitz Systems , 1996, SIAM Rev..

[44]  Phaedon C. Kyriakidis,et al.  Geostatistical Space–Time Models: A Review , 1999 .

[45]  C. R. Dietrich,et al.  Bounds on the size of nonnegative definite circulant embeddings of positive definite Toeplitz matrices , 1994, IEEE Trans. Inf. Theory.

[46]  Gene H. Golub,et al.  Matrix computations , 1983 .

[47]  Michael W.D. Davis,et al.  Kriging in a global neighborhood , 1984 .

[48]  Hermann G. Matthies,et al.  Sparse Data Representation of Random Fields , 2009 .

[49]  Werner G. Müller,et al.  Collecting Spatial Data: Optimum Design of Experiments for Random Fields , 1998 .

[50]  B. Kozintsev,et al.  Computations With Gaussian Random Fields , 1999 .

[51]  P. Kitanidis Parameter Uncertainty in Estimation of Spatial Functions: Bayesian Analysis , 1986 .

[52]  Patrick Reed,et al.  Comparative analysis of multiobjective evolutionary algorithms for random and correlated instances of multiobjective d-dimensional knapsack problems , 2011, Eur. J. Oper. Res..

[53]  M. W. Davis,et al.  Contouring Very Large Datasets Using Kriging , 1984 .

[54]  B. Khoromskij Tensors-structured Numerical Methods in Scientific Computing: Survey on Recent Advances , 2012 .

[55]  A. Zinober Matrices: Methods and Applications , 1992 .

[56]  C. R. Dietrich,et al.  Fast and Exact Simulation of Stationary Gaussian Processes through Circulant Embedding of the Covariance Matrix , 1997, SIAM J. Sci. Comput..

[57]  P. Jones,et al.  A European daily high-resolution gridded data set of surface temperature and precipitation for 1950-2006 , 2008 .

[58]  Analytical expressions of conditional mean, covariance, and sample functions in geostatistics , 1996 .

[59]  Steffen Börm,et al.  Hybrid cross approximation of integral operators , 2005, Numerische Mathematik.

[60]  James G. Nagy,et al.  Optimal Kronecker Product Approximation of Block Toeplitz Matrices , 2000, SIAM J. Matrix Anal. Appl..