Geostatistical Modeling and Prediction Using Mixed Precision Tile Cholesky Factorization

Abstract

Geostatistics represents one of the most challenging classes of scientific applications due to the desire to incorporate an ever increasing number of geospatial locations to accurately model and predict environmental phenomena. For example, the evaluation of the Gaussian log-likelihood function, which constitutes the main computational phase, involves solving systems of linear equations with a large dense symmetric and positive definite covariance matrix. Cholesky, the standard algorithm, requires O(n 3 ) floating point operators and has an O(n 2 ) memory footprint, where n is the number of geographical locations. Here, we present a mixed-precision tile algorithm to accelerate the Cholesky factorization during the log-likelihood function evaluation. Under an appropriate ordering, it operates with double-precision arithmetic on tiles around the diagonal, while reducing to single-precision arithmetic for tiles sufficiently far off. This translates into an improvement of the performance without any deterioration of the numerical accuracy of the application. We rely on the StarPU dynamic runtime system to schedule the tasks and to overlap them with data movement. To assess the performance and the accuracy of the proposed mixedprecision algorithm, we use synthetic and real datasets on various shared and distributed-memory systems possibly equipped with hardware accelerators. We compare our mixed-precision Cholesky factorization against the double-precision reference implementation as well as an independent block approximation method. We obtain an average of 1.6X performance speedup on massively parallel architectures, while maintaining the accuracy necessary for modeling and prediction.

2D Static Resource Allocation for Compressed Linear Algebra and Communication Constraints

“…In the field of weather forecast, maximum-likelyhood estimation method requires the resolution of large scale linear systems. In [1] the authors tackle the problem using Cholesky decomposition of a matrix which tiles are compressed using BLR technique. Since full rank diagonal tiles are associated to longer work time than compressed off-diagonal ones, they are treated specifically to ensure better load balancing between nodes, in distributed platforms implementation.…”

“…It is difficult to get access to matrices corresponding to actual use cases. For our experiments we therefore relied on synthetic matrices, which nevertheless reproduce the main characteristics of actual matrices, as can be encountered in various large scale simulation problems ( [1,10]). Thoses synthetic matrices allowed us to validate the approaches over a larger range of parameters and to ensure the reproducibility of the experiments.…”

Mixed-Precision Tomographic Reconstructor Computations on Hardware Accelerators

“…This may actually create new opportunities for a wide range of scientific applications. For instance, a similar customized approach has already been successfully applied in the context of a climate and weather prediction application [22], though no hardware support for lower precision has been demonstrated.…”

“…As a result, users can instantiate functors and define a suitable heuristic to choose the precision of the GEMM operation on a given tile, without being intrusive into the underlying linear algebra routine. This implementation can, therefore, be used by other applications [22] and still exploit the matrix specific structure at hand, simply by customizing the heuristic. A similar mechanism can be enforced on the solve stage, i.e., the backward and forward substitution, but in our case, all GEMMs of the substitutions can be computed with HP1, without hindering the solution accuracy.…”

Accelerating Geostatistical Modeling and Prediction With Mixed-Precision Computations: A High-Productivity Approach With PaRSEC

“…We introduce ExaGeoStat_PaRSEC, i.e., ExaGeoStat powered by PaRSEC, extending the approach in [4] to accelerate the Cholesky factorization by mixing FP64 doubleprecision (DP), FP32 single-precision (SP) and FP16 halfprecision (HP) to take advantage of the tensor cores of modern GPUs, e.g., NVIDIA V100s. Precision adaptation inveighs against predictable load-balancing, which therefore requires reliance on a dynamic runtime system to schedule computationally rich tasks of tile-sized granularity and data exchanges.…”