Parallel Approximation of the Maximum Likelihood Estimation for the Prediction of Large-Scale Geostatistics Simulations

Abstract

Maximum likelihood estimation is an important statistical technique for estimating missing data, for example in climate and environmental applications, which are usually large and feature data points that are irregularly spaced. In particular, the Gaussian log-likelihood function is the de facto model, which operates on the resulting sizable dense covariance matrix. The advent of high performance systems with advanced computing power and memory capacity have enabled full simulations only for rather small dimensional climate problems, solved at the machine precision accuracy. The challenge for high dimensional problems lies in the computation requirements of the log-likelihood function, which necessitates O(n 2 ) storage and O(n 3 ) operations, where n represents the number of given spatial locations. This prohibitive computational cost may be reduced by using approximation techniques that not only enable large-scale simulations otherwise intractable, but also maintain the accuracy and the fidelity of the spatial statistics model. In this paper, we extend the Exascale GeoStatistics software framework (i.e., ExaGeoStat 1 ) to support the Tile Low-Rank (TLR) approximation technique, which exploits the data sparsity of the dense covariance matrix by compressing the off-diagonal tiles up to a user-defined accuracy threshold. The underlying linear algebra operations may then be carried out on this data compression format, which may ultimately reduce the arithmetic complexity of the maximum likelihood estimation and the corresponding memory footprint. Performance results of TLR-based computations on shared and distributed-memory systems attain up to 13X and 5X speedups, respectively, compared to full accuracy simulations using synthetic and real datasets (up to 2M), while ensuring adequate prediction accuracy.

Leveraging PaRSEC Runtime Support to Tackle Challenging 3D Data-Sparse Matrix Problems

“…We compress off-band tiles and retain their most significant singular values (and associated vectors) above the accuracy threshold of 10 −8 (except in Section VIII-G), which ultimately yields an absolute numerical error of order 10 −9 in the solution of the linear system in Equation 1 to make it consistent as in [22]. This 10 −9 tolerance is sufficient to satisfy the prediction accuracy requirements of the 3D climate and weather prediction applications, as described in [8]. We employ the "band distribution" and a 2DBCDD for tiles off-band with a process grid P × Q (as square as possible) where P ≤ Q.…”

“…They highlight the impact of batch executions on GPUs to increase the hardware occupancy [2], [19], [20] for iterative solvers. They employ the HiCMA taskbased numerical library with the StarPU dynamic runtime system [8], [21] for attenuating load imbalance effects on distributed-memory systems in the context of TLR Cholesky factorizations. Further performance improvement of HiCMA has been obtained using the hybrid data distribution implemented in the PaRSEC dynamic runtime system [22], [23].…”

“…We extend PaRSEC with new functionality that takes into account the rank information in the task dataflow: (1) a dynamic data structure management driven by a performance model to reduce extra floating-point operations, (2) a dynamic memory allocator to further optimize memory footprint, (3) a rank-aware data distribution to cope with the workload imbalance, and (4) a recursive formulation of computational kernels to expose concurrency during the critical path. We use these features to leverage the performance of the TLR Cholesky factorization at the heart of the Maximum Likelihood Estimation (MLE) [8].…”

Performance Analysis of Tile Low-Rank Cholesky Factorization Using PaRSEC Instrumentation Tools

“…In particular, we compress off-diagonal tiles and retain their most significant singular values (and associated vectors) above the accuracy threshold of 10 −8 , which ultimately yields absolute numerical error of order 10 −9 in the solution of linear system in Equation 2in [5]. This 10 −9 tolerance is sufficient to satisfy the prediction accuracy requirements of the 3D climate and weather prediction applications, as described in [30]. We employ a process grid P × Q across computational nodes and make it as square as possible, with P < Q when this square is not possible.…”

Comments on: Data science, big data and statistics

“…Recent developments include multi-resolution methods (Nychka et al, 2015), nearest-neighbor approaches (Datta et al, 2016), and hierarchical low rank approximations (Huang and Sun, 2018). There is also a rich literature on taking advantage of computational tools in numerical analysis to approximate the exact covariance matrix (Abdulah et al, 2018a;Baugh and Stein, 2018). Recently, HPC has been firstly introduced to calculate the maximum likelihood estimator exactly and perform kriging prediction for exascale spatial data (Abdulah et al, 2018b).…”