Reshaping Geostatistical Modeling and Prediction for Extreme-Scale Environmental Applications

Cao, Qinglei; Abdulah, Sameh; Alomairy, Rabab; Pei, Yu; Nag, Pratik; Bosilca, George; Dongarra, Jack; Genton, Marc G.; Keyes, David E.; Ltaief, Hatem; Sun, Ying

doi:10.1109/sc41404.2022.00007

Cited by 10 publications

(3 citation statements)

References 34 publications

Supporting

Mentioning

Contrasting

Order By: Relevance

“…This combination of MP with TLR-MVMs is the critical synergism of our approach to SRI. It has also recently been employed on covariance matrices in a geospatial statistics application Gordon Bell Finalist paper (Cao et al, 2022).…”

Section: Introductionmentioning

confidence: 99%

High performance computing seismic redatuming by inversion with algebraic compression and multiple precisions

Hong,

Ltaief,

Ravasi

et al. 2024

The International Journal of High Performance Computing Applica

Self Cite

View full text Add to dashboard Cite

We present a high-performance implementation of Seismic Redatuming by Inversion (SRI), which combines algebraic compression with mixed-precision (MP) computations. Seismic redatuming entails the repositioning of seismic data recorded at the surface of the Earth to a subsurface level closer to where reflections have originated. Marchenko-based redatuming is the best-in-class SRI method from a theoretical point of view because it can accurately handle full seismic wavefields with minimal data pre-processing. However, it requires intensive data movement, a large memory footprint, and extensive computations due to the need for manipulating large, formally dense, complex-valued frequency matrices in the bandwidth of the recorded seismic data. More precisely, at the heart of Marchenko modeling operator lies a Multi-Dimensional Convolution (MDC) operator, which requires the application of repeated complex-valued Matrix-Vector Multiplications (MVMs) with chosen set of seismic frequency matrices. SRI is deemed to be too time-consuming for realistic domain sizes and therefore industrially impractical. To alleviate these computational bottlenecks, we first improve the memory footprint of the MDC operator by using Tile Low-Rank Matrix-Vector Multiplications (TLR-MVMs); this exploits the low-rank nature of seismic recordings in the frequency domain, used as the kernel of the MDC modeling operator. MP computations are further introduced to increase the arithmetic intensity of the operators. We validate the numerical robustness of our synergistic approach on a benchmark 3D synthetic seismic dataset and demonstrate its use on various hardware architectures. We develop several FP32/FP16/INT8 MP TLR-MVM implementations with a limited and controllable impact on the overall accuracy of application. We achieve up to 63X memory footprint reduction and 12.5X performance speedup against the traditional dense MVM kernel on a single NVIDIA A100 GPU. Moreover, linear scalability is attained when deploying our implementation on eight A100 GPUs. Our implementations of the MDC operator integrate with the open-source PyLops framework for large-scale, matrix-free inverse problems.

show abstract

Section: Introductionmentioning

confidence: 99%

High performance computing seismic redatuming by inversion with algebraic compression and multiple precisions

Hong,

Ltaief,

Ravasi

et al. 2024

The International Journal of High Performance Computing Applica

Self Cite

View full text Add to dashboard Cite

show abstract

“…Other studies assume the sparsity of the original covariance matrix by partially including correlation between some spatial locations and ignoring others, i.e., sparse inverse covariance methods [4], [5]. Moreover, the emergence of modern hardware architectures that support low-precision computation, such as NVIDIA GPUs, has facilitated the optimization of sparse inverse covariance methods by applying different precisions to various parts of the dense covariance matrix to reduce the computational complexity instead of ignoring them [6], [7], [8]. For the latter, different types of low-rank approximations are exploited, which allows faster computation and less memory consumption compared to the original dense matrix [9], [10], [11], [12].…”

Section: Introductionmentioning

confidence: 99%

GPU-Accelerated Vecchia Approximations of Gaussian Processes for Geospatial Data using Batched Matrix Computations

Pan,

Abdulah,

Genton

et al. 2024

ISC High Performance 2024 Research Paper Proceedings (39th International Conference)

Self Cite

View full text Add to dashboard Cite

Gaussian processes (GPs) are commonly used for geospatial analysis, but they suffer from high computational complexity when dealing with massive data. For instance, the log-likelihood function required in estimating the statistical model parameters for geospatial data is a computationally intensive procedure that involves computing the inverse of a covariance matrix with size n × n, where n represents the number of geographical locations in the simplest case. As a result, in the literature, studies have shifted towards approximation methods to handle larger values of n effectively while maintaining high accuracy. These methods encompass a range of techniques, including low-rank and sparse approximations. Among these techniques, Vecchia approximation is one of the most promising methods to speed up evaluating the log-likelihood function. This study presents a parallel implementation of the Vecchia approximation technique, utilizing batched matrix computations on contemporary GPUs. The proposed implementation relies on batched linear algebra routines to efficiently execute individual conditional distributions in the Vecchia algorithm. We rely on the KBLAS linear algebra library to perform batched linear algebra operations, reducing the time to solution compared to the state-of-the-art parallel implementation of the likelihood estimation operation in the ExaGeoStat software by up to 700X, 833X, 1380X on 32GB GV100, 80GB A100, and 80GB H100 GPUs, respectively, with the largest matrix dimension that can fully fit into the GPU memory in the dense Maximum Likelihood Estimation (MLE) case. We also successfully manage larger problem sizes on a single NVIDIA GPU, accommodating up to 1 million locations with 80GB A100 and H100 GPUs while maintaining the necessary application accuracy. We further assess the accuracy performance of the implemented algorithm, identifying the optimal settings for the Vecchia approximation algorithm to preserve accuracy on two real geospatial datasets: soil moisture data in the Mississippi Basin area and wind speed data in the Middle East.

show abstract

“…Although having a suboptimal O(N2) factorization complexity and O(N1.5) in memory consumption, it is still much better than the dense O(N3) complexity. HiCMA in Akbudak et al (2017); Cao et al (2019) and LORAPO in Cao et al (2020, 2022) provide a highly optimized implementation in the BLR format that delivers excellent performance on CPUs and distributed memory environment. Many of the BLR implementations including the ones mentioned above use the 2-D block-cyclic process distribution used by ScaLAPACK (Blackford et al 1997).…”

Section: Introductionmentioning

confidence: 99%

An inherently parallel ℋ²-ULV factorization for solving dense linear systems on GPUs

Ma,

Yokota

2024

The International Journal of High Performance Computing Applica

View full text Add to dashboard Cite

Hierarchical low-rank approximation of dense matrices can reduce the complexity of their factorization from [Formula: see text] to [Formula: see text]. However, the complex structure of such hierarchical matrices makes them difficult to parallelize. The block size and ranks can vary between the sub-blocks, which creates load imbalance. The dependency between the sub-blocks during factorization results in serialization. Since many sub-blocks are low-rank, their small computational load exposes the overhead of runtime systems. The combination of these factors makes it challenging to implement these methods on GPUs. In this work, we show that dense matrices can be factorized with linear complexity, while extracting the potential parallelism of GPUs. This is made possible through the [Formula: see text]-ULV factorization, which removes the dependency on trailing sub-matrices.

show abstract

Reshaping Geostatistical Modeling and Prediction for Extreme-Scale Environmental Applications

Cited by 10 publications

References 34 publications

High performance computing seismic redatuming by inversion with algebraic compression and multiple precisions

High performance computing seismic redatuming by inversion with algebraic compression and multiple precisions

GPU-Accelerated Vecchia Approximations of Gaussian Processes for Geospatial Data using Batched Matrix Computations

An inherently parallel ℋ²-ULV factorization for solving dense linear systems on GPUs

Contact Info

Product

Resources

About

Reshaping Geostatistical Modeling and Prediction for Extreme-Scale Environmental Applications

Cited by 10 publications

References 34 publications

High performance computing seismic redatuming by inversion with algebraic compression and multiple precisions

High performance computing seismic redatuming by inversion with algebraic compression and multiple precisions

GPU-Accelerated Vecchia Approximations of Gaussian Processes for Geospatial Data using Batched Matrix Computations

An inherently parallel ℋ2-ULV factorization for solving dense linear systems on GPUs

Contact Info

Product

Resources

About

An inherently parallel ℋ²-ULV factorization for solving dense linear systems on GPUs