Mixed-Precision GPU Acceleration for Confidence Regions of Spatial Excursion Sets

Mixed-Precision GPU Acceleration for Confidence Regions of Spatial Excursion Sets Xiran Zhang1 , Sameh Abdulah1 , Qinglei Cao2 , Joshua French3 , Rabab Alomairy1 , Marc G. Genton1 , David E. Keyes1 1

arXiv:2610.25039v1 [cs.DC] 18 Oct 2026

Computer, Electrical and Mathematical Sciences and Engineering (CEMSE) Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia 2 Saint Louis University, St. Louis, Missouri, USA 3 Department of Mathematical and Statistical Sciences, University of Colorado Denver, Denver, Colorado, USA Abstract—A spatial excursion set identifies the region where a field, such as temperature or pollutant concentration, exceeds a prescribed threshold. When the values are inferred from incomplete and noisy observations, inner and outer confidence regions quantify uncertainty in its spatial extent by providing lower and upper bounds on the excursion set at a specified confidence level. Constructing these regions requires repeated conditional simulations that account for joint spatial uncertainty, making dense covariance factorization and simulation major computational bottlenecks. We present a GPU-accelerated framework that combines adaptive tile-based mixed-precision arithmetic with distributed task execution to accelerate universal kriging, prediction-uncertainty estimation, conditional Gaussian simulation, and confidence region calibration. The central contribution is a parallel mixed-precision simulation workflow in which triangular matrix–matrix multiplication (TRMM) inherits the adaptive precision layout of the Cholesky factor and uses FP16 simulation operands, while covariance scaling and compensating operand rescaling accommodate the numerical range of reduced precision arithmetic. With appropriately chosen precision tolerances, the framework achieves empirical coverage close to the double-precision reference and the nominal 90% confidence level across weak, medium, and strong spatial correlation levels. Experiments on an ocean surface-salinity dataset show close agreement with the double-precision reference in spatial predictions and confidence regions. On four NVIDIA A100 GPUs (40 GB each), mixed-precision TRMM reaches 6.36× to 6.93× the theoretical FP64 peak throughput. On four NVIDIA GH200 Hopper GPUs, it reaches 10.20× to 11.03×. Each AMD MI250X node has eight GPU dies. Aggregate throughput rises by 9.1× to 11.5× from 1 to 16 nodes under weak scaling and by 1.6× to 1.9× from 4 to 16 nodes under strong scaling. These results show that adaptive mixed-precision accelerates large-scale excursion set analysis while preserving the statistical accuracy required for spatial uncertainty quantification. Index Terms—Geostatistical excursion sets, confidence regions, spatial uncertainty quantification, conditional Gaussian simulation, adaptive mixed-precision arithmetic, GPU acceleration, distributed task-based computing, PaRSEC.

I. I NTRODUCTION Spatial statistics enables the analysis of geographically distributed phenomena whose measurements exhibit spatial dependence, including temperature, precipitation, and pollutant concentration [4], [27], [31], [33]. Observations are typically available at a finite set of locations and may contain measurement error, so statistical models are needed to

infer an underlying latent process and predict its behavior at unobserved locations. In many applications, however, the scientific objective extends beyond predicting individual values to identifying regions where the measurements exceed a meaningful threshold. Examples include areas experiencing extreme temperatures or pollutant concentrations above a regulatory limit. Such regions are known as excursion sets [8], [39]. Let Z(s), s ∈ D, denote a latent spatial process over a domain D ⊂ R2 , and let u be a prescribed threshold. The positive excursion set consists of locations s where the process is at least u (Section III defines it formally). Because Z(s) is inferred from incomplete and noisy observations, the spatial extent and boundary of this set are uncertain. Neither predictions nor exceedance probabilities resolve this uncertainty. High confidence at individual locations does not imply simultaneous correctness across an entire region. Inner and outer confidence regions address this limitation by providing spatial lower and upper bounds on the unknown excursion set [8], [22], [24], [35]. The inner region identifies locations confidently within the excursion set, whereas the outer region provides a conservative envelope that contains it at the specified confidence level. The area between these bounds characterizes uncertainty in the excursion boundary, supporting decisions that depend on the location and extent of threshold exceedance. For Gaussian process models, constructing these confidence regions requires spatial prediction, uncertainty estimation, and repeated conditional simulations to calibrate spatially simultaneous inference [21]. These computations depend on dense covariance matrices and expensive linear algebra operations, including Cholesky factorization, triangular solves, and triangular matrix–matrix multiplication (TRMM). For a covariance matrix of order n, dense Cholesky factorization requires O(n3 ) arithmetic and O(n2 ) storage. Consequently, spatial resolution and simulation count can substantially increase computational cost. Supporting operations, including randomnumber generation, prediction-error calculations, and reductions over simulated fields, add data-movement and execution overhead. Efficient confidence region construction therefore requires coordinated optimization of the complete workflow.

GPU-accelerated architectures offer substantial computational throughput for these workloads, particularly when operations can exploit reduced precision arithmetic. Doubleprecision (FP64) remains necessary for numerically sensitive computations, while single-precision (FP32) and halfprecision (FP16) can increase throughput and reduce storage and data movement. Hardware features such as NVIDIA Tensor Cores [30] and AMD Matrix Core Units [34] further accelerate lower-precision operations. Exploiting these capabilities in spatial uncertainty quantification requires precision choices that preserve application-level statistical accuracy, rather than numerical performance alone. Mixed-precision algorithms combine arithmetic formats within a computation, assigning higher precision to sensitive operations and lower precision where reduced accuracy is acceptable. Precision selection may follow heuristic rules [2], [3] or criteria informed by numerical analysis [13], [25]. In tiled linear algebra, individual matrix tiles can use different precisions according to their numerical properties. Previous studies have demonstrated the effectiveness of such approaches for Cholesky factorization and large-scale scientific applications [1], [7], [13], [29], [32], [36]. However, incorporating these kernels into confidence region construction introduces additional challenges such as the precision layout must remain compatible with simulation operations, reduced precision operands must accommodate the numerical range of the covariance, and simulation errors must be assessed through their effects on statistical coverage. The precision selection must also adapt to the covariance structure, since a single prescribed tolerance may permit unsuitable FP16 assignments across different datasets and spatial-correlation levels. We address these challenges using a distributed GPUaccelerated framework built on PaRSEC [10]. The framework extends adaptive mixed-precision Cholesky factorization to conditional simulation through a parallel TRMM implementation. Automatic tolerance refinement and covariance scaling adapt reduced precision execution to the covariance structure, while empirical coverage experiments assess the statistical accuracy of the application. Beyond the dominant linear algebra operations, the framework accelerates supporting stages of confidence region construction. It computes pointwise mean squared prediction errors (MSPEs) without forming the full prediction-error covariance matrix, generates standard normal samples directly on GPUs, and fuses the minimum and maximum reductions required for confidence-region calibration. These optimizations reduce data movement and execution overhead throughout the statistical workflow. The main contributions of this paper are as follows: • We develop a distributed GPU-accelerated framework for excursion set confidence regions using PaRSEC. It integrates universal kriging, prediction-uncertainty estimation, conditional simulation, and confidence region calibration. • We implement parallel mixed-precision TRMM that inherits the adaptive precision layout of the Cholesky factor

and uses FP16 simulation operands. Task dependencies ensure correct in-place execution. • We automatically refine the precision-selection tolerance using FP16 representability checks and the relativemagnitude criterion. Covariance scaling and compensating operand rescaling support reduced precision execution, with statistical fidelity assessed empirically. • We accelerate prediction uncertainty estimation without forming the full prediction-error covariance matrix. GPU Gaussian-sample generation and fused calibration reductions further reduce workflow overhead. • We validate statistical accuracy across three spatialcorrelation levels and on the HYCOM surface-salinity dataset [26]. Results show coverage close to the nominal 90% level with appropriately chosen tolerances and close agreement with the double-precision reference. • We evaluate mixed-precision TRMM performance across different GPU architectures, demonstrating substantial single-node throughput gains and 9.1×–11.5× aggregatethroughput improvement from 1 to 16 AMD-based Frontier nodes. II. R ELATED W ORK This section reviews prior work on confidence region construction for excursion sets and tile-based mixed-precision linear algebra relevant to our proposed framework. Confidence Regions for Excursion Sets: Quantifying uncertainty in excursion set boundaries remains a significant challenge because threshold-based approaches to prediction surfaces do not produce statistically valid confidence regions. Lindgren and Rychlik [28] investigated confidence sets for level curves by leveraging crossing intensities to characterize the probabilistic behavior of contour crossings in stochastic processes. However, spatial smoothing of predictors often results in overly regular and biased boundaries [37]. Loss-based methods proposed by Craigmile and Cressie [18] and Zhang et al. [37] partially address this limitation but do not provide simultaneous confidence or credible regions. Recent advances have constructed regions that account for joint spatial uncertainty. French and Sain [24] introduced a framework based on kriging, conditional simulation, and hypothesis testing. This framework was later extended to contour uncertainty [22] and Bayesian contexts [23]. Concurrently, Bolin and Lindgren [8] developed a unified approach for excursion and contour uncertainty regions in latent Gaussian models, formulating the problem in terms of joint exceedance probabilities and addressing multiple-testing challenges through parametric set constructions and efficient numerical integration. The same framework has been implemented and scaled to large problem sizes in [39]. Although these methods provide strong statistical foundations, they introduce substantial computational challenges because they repeatedly evaluate high-dimensional probability distributions. This motivates scalable approaches. Tile-Based Mixed-Precision Parallel Linear Algebra: Task-based runtimes such as StarPU [6] and PaRSEC [9] support tiled linear algebra on heterogeneous systems. Previous

studies combine these runtimes with mixed-precision kernels to accelerate geostatistical modeling and large-scale scientific applications [1]–[3], [13], [14], [29], [38]. Early studies demonstrated that selectively applying lower precision within numerical kernels can substantially reduce computational costs without compromising solution quality. In geostatistical modeling, Abdulah et al. [3] introduced a tiled mixed-precision Cholesky factorization in which tiles near the diagonal are computed in high precision, while offdiagonal tiles associated with weaker spatial correlations are processed in lower precision. This strategy exploits the natural decay of spatial covariance structures to improve performance while incurring minimal accuracy loss. This method has since been extended to heterogeneous architectures using task-based runtimes that support adaptive precision selection and efficient data movement across distributed resources [2]. Moreover, Cao et al. [12] proposed automated precision-conversion strategies to reduce execution time, memory traffic, and energy consumption on modern GPU platforms. More recent work has explored multi-precision and out-of-core Cholesky factorizations that combine fine-grained task parallelism, adaptive precision techniques, and advanced interconnect technologies to achieve significant performance improvements in large-scale dense matrix computations [32]. III. M ATHEMATICAL F ORMULATION OF C ONFIDENCE R EGIONS OF E XCURSION S ETS We consider spatial data defined over a geographic domain D ⊂ R2 , where each location s ∈ D corresponds to a point in space (e.g., latitude and longitude coordinates). The observed process Y (s) represents measurements collected across the domain, such as temperature or air pollution. These observations typically exhibit spatial dependence, meaning that measurements at nearby locations tend to be more similar than those farther apart. The spatial data can be modeled as Y (s) = Z(s) + ε(s), s ∈ D, where Z(s) is the underlying (latent) spatial process we aim to estimate, and ε(s) represents measurement noise. For a threshold u ∈ R, the positive and negative excursion − sets are A+ u (Z) = {s ∈ D : Z(s) ≥ u} and Au (Z) = {s ∈ D : Z(s) ≤ u}, respectively. We focus on positive excursion sets. Negative excursions can be treated by reversing the sign of the process and threshold. Since excursion sets are inferred from noisy and incomplete observations, their spatial boundaries are uncertain. Inner and outer confidence regions provide simultaneous probabilistic bounds on the excursion set [35]. The outer confidence region O+ Eu,α (Z) contains the true excursion
set with confidence O+ level 1 − α: Pr A+ u (Z) ⊆ Eu,α (Z) = 1 − α. In contrast, I+ the inner confidence region Eu,α (Z) contains only those locations that exceed
the threshold with high confidence, I+ Pr Eu,α (Z) ⊆ A+ u (Z) = 1 − α. The inner and outer confidence regions provide lower and upper bounds on the unknown excursion set, characterizing the uncertainty in the spatial extent and boundary of threshold exceedance events.

Fig. 1: Positive excursion sets and confidence regions: (a) latent process Z(s) and threshold u; (b) excursion set A+ u (Z); (c) inner and outer bounds at confidence level 1 − α, with the intervening area indicating boundary uncertainty.

To identify these regions, we use predictions from a Gausb sian process model. Let Z(s) be the predicted value and σK (s) its uncertainty. We define the standardized score T (s) =

b −u Z(s) , σK (s)

(1)

which measures the evidence that a location exceeds the threshold, with larger values indicating stronger evidence. In practice, the domain D is discretized into a grid of m locations. Confidence regions are constructed by thresholding T (s) using a critical value estimated via conditional simulation. Specifically, multiple realizations of the spatial field are generated conditionally on the observed data, and the distribution of extreme values of T (s) is used to determine a threshold that guarantees the confidence regions have the desired confidence level. Specifically, the resulting regions are obtained as O+,m bu,α E (Z) = {gi : T (sgi ) ≥ TbαO+ },

(2)

where TbαO+ is the α-quantile of the distribution of the minimum values of T (s) for s ∈ A+ u (z) (see French and Sain [24] for more details). A similar approach is used for the inner region. Figure 1 illustrates the excursion set and its confidence regions. The following sections describe their computation using tiled GPU linear algebra. IV. T ILE -C ENTRIC M IXED -P RECISION L INEAR A LGEBRA Tile-based algorithms are a fundamental strategy in contemporary dense linear algebra libraries and runtime systems for addressing the scalability limitations of traditional panel-based factorizations, including those implemented in LAPACK [11]. By partitioning matrices into smaller tiles and representing computations as a directed acyclic graph (DAG) of finegrained tasks, tiled methods expose substantially more parallelism while preserving data dependencies. This execution model improves scalability and cache locality, and enables computation to overlap with communication across multicore CPUs, GPU-accelerated platforms, and distributed-memory systems. The tile abstraction also provides a natural granularity for applying different numerical representations or computational strategies to different regions of a matrix, including mixed-precision arithmetic [3] and low-rank compression [5], [15], [17]. These techniques exploit the numerical properties

of individual tiles to reduce computational and data-movement costs while maintaining control over numerical accuracy. Mixed-precision arithmetic can be applied at the tile level to balance numerical accuracy and computational efficiency. Rather than using uniform precision for the entire matrix, individual tiles can be represented and processed using different arithmetic formats according to application-specific or numerical criteria. For example, in spatial statistics, tiles associated with strongly correlated regions may require double-precision, whereas tiles with weaker correlations can be processed using single-precision, half-precision, or other reduced precision formats [3]. This adaptive assignment can improve computational performance and reduce memory traffic while maintaining the statistical accuracy required by the application [13]. A. PaRSEC: A Task-Based Runtime for Scalable MixedPrecision Computing PaRSEC is a task-based runtime system that enables scalable execution of scientific applications across distributed, heterogeneous CPU and GPU architectures [9], [10]. In this work, we employ PaRSEC’s Parameterized Task Graph (PTG) model [20], which represents computations as parameterized DAGs of fine-grained tasks with explicit data dependencies. Based on these dependencies, PaRSEC asynchronously schedules tasks, manages distributed data movement and communication, and overlaps communication with computation. This dataflow execution model exposes parallelism, improves data locality, and reduces synchronization overhead at scale. For mixed-precision computing, PaRSEC provides the runtime infrastructure to execute tasks operating at different numerical precisions while preserving their data dependencies. In tile-based dense linear algebra, individual computational kernels can employ different precision formats, such as FP64, FP32, and FP16, with the necessary data transfers and precision conversions coordinated during execution. This separation of numerical precision selection from task execution enables mixed-precision algorithms to exploit heterogeneous hardware efficiently. Previous studies have demonstrated the scalability of PaRSEC-based mixed-precision methods for Gaussian process inference and environmental uncertainty quantification on distributed GPU systems [1], [15], [16]. B. Adaptive Mixed-Precision Cholesky Decomposition Our confidence region detection workflow employs the adaptive mixed-precision Cholesky decomposition introduced in [13]. The method exploits the correlation structure of spatial covariance matrices, where correlations typically decrease with increasing spatial separation. Consequently, tiles near the diagonal generally contain stronger correlations, while distant offdiagonal tiles often have smaller numerical contributions. This structure provides opportunities to reduce arithmetic precision selectively without significantly affecting numerical accuracy. The method partitions the covariance matrix into tiles and assigns computational precision according to their numerical significance [2], [3]. Diagonal tiles retain FP64 to maintain numerical stability, while less significant off-diagonal tiles

can employ FP32 or FP16. Compared with uniform FP64 computation, this tile-based mixed-precision strategy reduces arithmetic cost and memory traffic while controlling the impact of reduced precision on application-level accuracy. Its performance and accuracy have been extensively evaluated in previous studies of large-scale spatial modeling [1], [13]. In this work, we build upon the established mixed-precision Cholesky method and focus on optimizing its integration into the confidence region detection workflow. The following sections describe the proposed workflow optimizations and their impact on end-to-end performance. C. Adaptive MP Triangular Matrix–Matrix Multiplication The confidence region algorithm generates correlated simulations through triangular matrix–matrix multiplication (TRMM), leveraging the factor obtained from the adaptive mixed-precision Cholesky factorization (Further implementation details are provided in Section V). Specifically, we implement the left-sided, lower-triangular operation B ← αLs B, where Ls is the scaled Cholesky factor prepared in Algorithm 2 and preserves the tile-wise precision assignment determined during the Cholesky factorization. Let B (0) denote the input matrix before the TRMM operation. The update for each destination tile Bij is expressed as (0)

Bij ← α(Ls )ii Bij + α

i−1 X

(0)

(Ls )ik Bkj .

(3)

k=0

The diagonal contribution is computed by a TRMM task, while the off-diagonal contributions are accumulated through a sequence of dependent GEMM tasks. In our simulation pipeline, the TRMM coefficient is set to α = 1. The two operands use different precision policies. The factor Ls retains its adaptive tile precision map PL (i, k) ∈ {FP64, FP32, FP16}, whereas every right-hand-side tile is assigned FP16. The Frobenius norm of the rescaled righthand side is computed for diagnostic purposes only. It does not affect precision selection. Thus, TRMM inherits adaptivity from Cholesky rather than performing a new adaptive precision selection for the right-hand side. Algorithm 1 summarizes the implementation with FP16 storage for the right-hand side and destination tiles. For a diagonal update, the diagonal factor tile and the destination tile are represented in FP32, an FP32 TRMM kernel is applied, and the result is converted back to FP16. This conversion is necessary because the implementation lacks a native FP16 TRMM kernel. For each off-diagonal update, the contributing factor tile is converted from its stored precision to FP16 as needed, and an FP16 GEMM updates the destination with β = 1. These temporary conversions leave the stored precision map of Ls unchanged. After all updates finish, the result is converted to FP64 for the subsequent confidence region computations. Here, castp (X; P) converts a tile from its stored precision P to FPp, with no conversion required when the precisions coincide. The notation GEMM16 (A, X, C; α, β) denotes the FP16 update C ← αAX + βC, including storage of the updated destination in FP16.

Algorithm 1 Mixed-Precision Tile TRMM Input: Lower-triangular factor Ls with NT tile rows; random matrix Mrand with KT tile columns; factor precision map PL ; applied Cholesky scaling factor s; adaptive mixedprecision flag AMP Output: Simulation matrix newsim stored in FP64 Configuration: α = 1; FP16 destination storage in the mixed-precision path. 1: if AMP is 0 then 2: newsim ← TRMM64 (Ls , Mrand ) 3: return newsim 4: end if 5: B ← Mrand /s 6: Compute ∥B∥F 7: PB (i, j) ← FP16 for all tiles 8: B ← cast16 (B) 9: for i = NT − 1, NT − 2, . . . , 0 do 10: for j = 0, 1, . . . , KT − 1 do
11: D ← cast32 (Ls )ii ; PL (i, i) 12: X ← cast32 (Bij )
13: Bij ← cast16 TRMM32 (D, X) 14: for k = 0, 1, . . . , i − 1 do
15: A ← cast16 (Ls )ik ; PL (i, k) 16: Bij ← GEMM16 (A, Bkj , Bij ; α = 1, β = 1) 17: end for 18: end for 19: end for 20: newsim ← cast64 (B) 21: return newsim

The descending row order in Algorithm 1 provides a valid sequential ordering of the in-place computation: when row i is updated, every source tile Bkj with k < i still contains its original input value. PaRSEC expresses this requirement through task dependencies rather than a global sequential loop. Each diagonal task precedes the off-diagonal updates to its destination tile, and each source tile is protected from overwrite until all GEMM tasks that consume its original value have completed. Tasks whose dependencies are satisfied can therefore execute concurrently. V. PARALLEL A LGORITHM FOR C ONFIDENCE R EGION C ONSTRUCTION The algorithm for confidence region construction includes three stages: kriging prediction and uncertainty estimation, conditional simulation, and calibration of the inner and outer regions. Algorithm 2 summarizes its computational dependencies. The algorithm partitions matrices into tiles and distributes them over a P × Q process grid. PaRSEC manages task dependencies, data movement, and execution on the available hardware resources. np Let Do contain no observation locations, Dp = {si }i=1 contain np = m2 prediction locations on an m × m grid, and n = np + no . We let nsim denote the number of simulations.

With prediction locations ordered first, the joint simulation covariance V ∈ Rn×n is   ⊤ Vp Vop V = , (4) Vop Σy with Σy = Vo + τ 2 Ino . The blocks Vp and Vo are the latent-process covariances at prediction and observation locations, respectively, and Vop contains their cross-covariances. Measurement-error variance τ 2 is added only to the observation block Σy . We write Xo and Xp for the trend design matrices, W for the kriging weights, ẑ for the prediction vector, and v for the pointwise MSPE vector. Vector square roots, divisions, and comparisons are componentwise; ⊙ denotes the Hadamard product, and colsum returns the column sums as a column vector. The kriging stage uses FP64 arithmetic. Adaptive mixedprecision is applied to the joint covariance factorization and the subsequent triangular multiplication, while the correction and calibration stages operate on results stored in FP64. The following subsections describe these stages and the supporting precision-selection and tile kernels. A. Kriging Prediction and Uncertainty The first stage computes the universal kriging predictor and its MSPE at every prediction location. These quantities enter the standardized test statistic in Equation 1 and the confidence regions in Equation 2. For fixed data, locations, trend covariates, and covariance parameters, they are computed once before simulation. Let Xo ∈ Rno ×k and Xp ∈ Rnp ×k be the trend design matrices. The current two-dimensional implementation uses an intercept and two spatial coordinates, so k = 3. The universal kriging weights are h i −1 ⊤ W = Σ−1 Vop − Xo (Xo⊤ Σ−1 (Xo⊤ Σ−1 y y Xo ) y Vop − Xp ) ,

(5)

and the prediction vector is ẑ = W ⊤ yo . The weights satisfy Xo⊤ W = Xp⊤ to ensure unbiasedness for the trend model. The implementation evaluates these expressions without forming explicit inverses. A POSV call factors Σy and solves Σy F = Xo . Matrix products form G = Xo⊤ F and F ⊤ Vop − Xp⊤ , after which a second POSV call solves the small trend system. The retained factor of Σy is then reused in two TRSM calls to obtain W , avoiding another factorization of the observation covariance. A GEMM computes W ⊤ yo . For column wi and column ci of Vop , the pointwise MSPE 2 is vi = σK (si ) = (Vp )ii + wi⊤ Σy wi − 2wi⊤ ci .. The trendestimation contribution is included through the universal kriging weights. The prediction target is the latent process, so error variance enters Σy but is not added to (Vp )ii . To evaluate all MSPEs, the implementation restores the unfactored Σy and computes U = Σy W using CR_dsymm. Hadamard products and GEMM reductions against a vector of ones then give v = θ0 1np + colsum(W ⊙ U ) − 2 colsum(Vop ⊙ W ). (6)

Algorithm 2 Distributed Adaptive Mixed-Precision Confidence Region Construction Input: Locations Do , Dp , observations yo , designs Xo , Xp ; parameters θ, τ 2 ; simulation count nsim and seed; threshold u, confidence level c; precision settings and computing resources Output: Inner and outer regions Cin , Cout 1: Initialize PaRSEC and distributed tiles; generate Vp , Vop , Vo Kriging prediction and uncertainty 2: Σy ← Vo + τ 2 I; factor Σy = Lo L⊤ o and solve Σy F = Xo 3: G ← Xo⊤ F ; solve GH = F ⊤ Vop − Xp⊤ 4: Solve Lo Y = Vop − Xo H, then L⊤ o W = Y ▷ Reuse Lo 5: ẑ ← W ⊤ yo ; U ← Σy W 6: v ← θ0 1np + colsum(W ⊙ U − 2Vop ⊙ W ) Conditional simulation 7: Assemble V as in Equation 4; s ← 1 8: Generate Mrand ∈ Rn×nsim with independent N (0, 1) entries 9: if adaptive mixed-precision is enabled then 10: Select the tile precision map and applied scale s ▷ Refine tolerance in automatic mode 11: end if 12: Ls ← chol(s2 V ); B ← Mrand /s ▷ Use selected tile precisions 13: if adaptive mixed-precision is enabled then 14: Convert all tiles of B to FP16 15: end if 16: M ← Ls B; convert the result to FP64 if needed 17: Partition M = [Mp⊤ , Mo⊤ ]⊤ ; S ← Mp + W ⊤ (yo 1⊤ nsim − Mo ) Confidence region √ calibration 18: t ← (ẑ − u1np )/ v 19: For each j, compute ψj ← mini:Sij ≥u ti and ϕj ← maxi:Sij ≤u ti ▷ Fused scans 20: Gather and sort each sample, separately omitting emptyset results 21: qout ← Q1−c (ψ); qin ← Qc (ϕ) 22: return Cin = {si : ti > qin }, Cout = {si : ti ≥ qout }

Here, the constant marginal prediction variance is initialized from theta[0]. This calculation avoids forming the full np × np prediction-error covariance matrix. B. Conditional Simulation The second stage generates nsim realizations from the joint covariance in Equation 4 and applies the universal-kriging conditioning correction. An independent standard-normal matrix Mrand ∈ Rn×nsim is generated with Desc_rnorm. Its CUDA implementation follows the seeded LAPACK random-number stream and Box– Muller draw ordering used by the CPU path. When that mode is enabled, the joint covariance is factorized using the adaptive mixed-precision framework in Section IV. Otherwise, the

factorization uses FP64. Algorithm 2 gives the covariance factorization, input compensation, and triangular multiplication, with the applied scale s selected as described in Section V-E. Algorithm 1 describes the mixed-precision multiplication and returns the result in FP64 for the conditioning correction. Partition M = [Mp⊤ , Mo⊤ ]⊤ into prediction and observation rows. The conditional realizations are
S = Mp + W ⊤ yo 1⊤ (7) nsim − Mo . The residual is formed by negating Mo and broadcasting yo across its columns; a GEMM applies the correction to Mp . No separate observation-noise matrix is generated, since that noise is already represented by Σy in Equation 4. In exact arithmetic, the simulation mean is W ⊤ yo , and the diagonal of its covariance is the MSPE vector in Equation V-A. This provides consistency between the prediction uncertainty and the simulation-based calibration. C. Outer and Inner Confidence Regions The third stage calibrates the regions using the conditional simulations and the standardized statistic in Equation 1. On the prediction grid, we write ti = T (si ), where ẑi and vi are the prediction and pointwise MSPE defined in Equation V-A. Positive MSPEs are required. Let c = 1 − α denote the requested confidence level. For simulation j, define A+ j = {i : Sij ≥ u},

A− j = {i : Sij ≤ u}.

(8)

The two statistics used for calibration are ψj = mini∈A+ ti , ϕj = maxi∈A− ti . The implementation computes j j both extrema in one pass over the simulation tiles using Desc_col_summary_index_both. The resulting values are gathered and sorted separately. If a threshold set is empty, its sentinel value is excluded from the corresponding sample. Thus, the two samples may have different sizes. For a nonempty sorted sample a(1) ≤ · · · ≤ a(N ) , the implementation uses the empirical quantile Qp (a) = a(max{1,min{N,⌈pN ⌉}}) .

(9)

Both samples must contain at least one valid extremum for calibration. The critical values are qout = Qα (ψ), qin = Q1−α (ϕ), and the regions are Cout = {si ∈ Dp : ti ≥ qout }, Cin = {si ∈ Dp : ti > qin }. These correspond to the outer and inner regions in Equation 2. The nonstrict outer comparison and strict inner comparison follow the implementation. Descriptor kernels evaluate the comparisons to produce the final indicator vectors. D. Automatic fixedacc Selection Criteria The parameter fixedacc controls the tile-precision assignment used by the adaptive factorization. A smaller value generally keeps more tiles at higher precision, whereas a larger value allows more reduced precision tiles. In automatic mode, the implementation starts with fixedacc = 10−4 and recomputes the assignment, decreasing this value by a factor of 0.9 whenever either check fails. Iteration stops when both checks pass or no FP16 tiles remain. These checks

are inexpensive screening heuristics; they do not bound the complete factorization error or guarantee statistical coverage. The checks use entry extrema and tile Frobenius norms of the unscaled joint simulation covariance V . For the positive HP HP Matérn covariance entries considered here, let vmax and vmin denote the maximum and minimum entries in FP16-assigned tiles. The implementation reduces the stored entry values directly; it does not take absolute values or exclude zeros. The HP following range interpretation therefore assumes vmin > 0. 1) FP16 representability range: The first check is RHP =

HP xFP16 65504 vmax max < FP16 = −24 ≈ 1.10 × 1012 . HP 2 vmin xmin,sub

(10)

Uniform scaling preserves the ratio between the largest and smallest entries. A ratio exceeding this limit prevents both extrema from lying between the smallest positive FP16 subnormal and the largest finite FP16 value under a common scale. The implementation uses a strict inequality and reduces fixedacc when it fails. 2) Relative magnitude of FP16 tiles: The second check compares FP16-assigned entries and tile norms with their HP /vmax and ρ2 = global counterparts:  ρ1 = vmax max(i,j)∈T HP ∥Vij ∥F max(i,j) ∥Vij ∥F , where T HP is the set of FP16-assigned tiles and vmax is the maximum entry across all precision classes. The acceptance criterion is 10−4 max(ρ1 , ρ2 ) < , εFP16

εFP16 = 2−10 .

(11)

This heuristic limits the relative magnitude of FP16-assigned data but does not bound accumulated rounding error. With a user-specified tolerance, failed checks are reported without automatic refinement. E. Automatic Scaling of the Covariance Matrix The joint covariance can contain entries spanning many orders of magnitude. For the Matérn model, off-diagonal covariances decay with distance, whereas diagonal entries contain the marginal variance and, in the observed block, measurement-error variance. Small entries in reduced precision tiles may therefore become subnormal or round to zero. Scaling increases their magnitudes subject to the overflow limits of the assigned formats. For each present reduced precision format f ∈ (f ) (f ) {FP32, FP16}, let vmin > 0 and vmax be the entry extrema in its assigned tiles. To place the scaled entries within the normal range, the factor s would need to satisfy q q (f ) (f ) (f ) (f ) ℓf = xmin,norm /vmin ≤ s ≤ xmax /vmax = uf . (12) The routine combines the present formats using ℓ = maxf ℓf and us = min{65504, minf uf }. For the positive upper bounds used here, it chooses an integer candidate approximately one unit below the upper bound, scand = ⌊us ⌋ − 1. If the normal-range lower bound exceeds the upper bound, it retains the upper-bound-based choice and reports diagnostics using the subnormal limits. It does not guarantee that every

scaled entry remains normal, nor do these input bounds control all intermediate values produced by Cholesky. The driver applies the candidate when its integer value exceeds one; otherwise, no scaling is performed. We use s for this factor, so s = 1 in the unscaled case. The compensation in Algorithm 2 preserves the simulation covariance in exact arithmetic and is applied before FP16 conversion. F. Half-Precision Casting of the Random Matrix The compensated TRMM input prepared in Algorithm 2 has entry variance s−2 before rounding. Standard-normal values are concentrated far below the FP16 overflow limit, and input compensation further reduces their magnitudes. However, it can increase the number of subnormal or zero values after conversion. The suitability of this representation therefore depends on the resulting simulation errors, not solely on the small magnitudes of the discarded entries. The numerical and statistical experiments in Section VI assess its effect alongside finite-simulation variability. The experiments use the default FP16 destination-storage path with optional FP32 accumulation disabled. VI. E XPERIMENTAL E VALUATION A. Experimental Testbeds We run experiments on four GPU-accelerated platforms: single-node of Polaris, Guyot, and JUPITER, and Frontier, using 1, 4, and 16 nodes. Polaris at the Argonne Leadership Computing Facility. Each compute node contains one 32-core AMD EPYC Milan CPU, 512 GB of DDR4 memory, and four NVIDIA A100 GPUs, each has with 40 GB of HBM2 memory. • Guyot at Innovative Computing Laboratory. A compute node contains two EPYC 7742 64-Core CPUs at 2.25GHz, 2063 GB of main memory, and eight NVIDIA A100-SXM4-80GB GPUs. • JUPITER at the Jülich Supercomputing Centre. Each node contains four NVIDIA GH200 Grace Hopper Superchips. Each GH200 combines a 72-core NVIDIA Grace CPU with 120 GB of LPDDR5X memory and an NVIDIA Hopper GPU with 96 GB of HBM3 memory. • Frontier at the Oak Ridge Leadership Computing Facility (OLCF) is an HPE Cray EX system. Each node features a 64-core AMD EPYC CPU with 512 GB of DDR4 memory and four AMD MI250X accelerators, providing eight GPU dies with 64 GB of HBM2E each. •

The implementation is built on HiCMA-X (https://github. com/ecrc/hicma-x) and the PaRSEC task-based runtime system. The software stack further includes the GNU Scientific Library (GSL) and optimized BLAS and LAPACK implementations available on each platform. Each experiment includes an initial warm-up execution followed by multiple measured runs. We report the median execution time together with the observed run-to-run variability.

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

(a) Weak (1 : 0.033 : 0.5)

1

0

0.2

0.4

0.6

0.8

(b) Medium (1 : 0.1 : 0.5)

1

0

0.91

0.90

0.90

0.89

0.89

0.88 0.87 0.86 0.85 0.84 0.83

target = 0.9 strong (0.234) medium (0.1) weak (0.033)

4 5 6 7 8 9 10 11 12 13 14 DP

log10 (fixedacc) (DP at rightmost tick)

(a) Inner coverage rate

0.88 0.87 0.86 0.85 0.84 0.83

target = 0.9 strong (0.234) medium (0.1) weak (0.033)

4 5 6 7 8 9 10 11 12 13 14 DP

log10 (fixedacc) (DP at rightmost tick)

(b) Outer coverage rate

Fig. 3: Empirical coverage rates of the inner and outer confidence regions obtained with adaptive mixed-precision and full doubleprecision across spatial-correlation levels. The horizontal axis shows the adaptive precision parameter fixedacc, and the dashed horizontal line marks the nominal confidence level of 0.9.

For the representative cases shown in Figure 2, we compare which cells are included in the confidence regions when using DP versus fixedacc = 10−4 . Under weak correlation (β = 0.033), the inner and outer confidence regions differ in their inclusion for only 0.54% and 0.47% of cells, respectively. The corresponding percentages are 0.08% and 0.12% under medium correlation (β = 0.1), and 0.03% and 0.06% under strong correlation (β = 0.234). Although the approximate setting does not identically preserve the statistical properties of the confidence regions constructed using DP, there is little practical difference in the resulting confidence regions. C. HYCOM Surface-Salinity Real Dataset Analysis

0.2

0

0.91

Outer coverage rate

1) Synthetic Data Generation: We generate synthetic datasets over the two-dimensional unit square D = [0, 1]2 . For a prescribed number of observations no , we sample the observation locations Do independently from a uniform distribution over D using a fixed random seed. Prediction np locations Dp = {si }i=1 are placed on a regular m × m grid, 2 with np = m . The grid points are located at the centers of equally sized cells, with spacing h = 1/m and coordinates ((2j + 1)h/2, (2i + 1)h/2) and i, j = 0, . . . , m − 1. The prediction and observation locations are combined into a single set of size n = np + no , with the prediction locations ordered first. A covariance matrix V is then constructed jointly over Dp ∪ Do using the specified spatial covariance parameters. This joint construction preserves the spatial dependence between the latent field values at prediction and observation locations. The deterministic mean is modeled using the design vector X(s) = (1, x, y)⊤ and regression coefficient vector b = (0, 15, 10)⊤ , yielding m(s) = X(s)⊤ b = 15x + 10y. A realization of the latent spatial process is generated as Z = Xb + Lη, η ∼ N (0, I), and LL⊤ = V, where X is the design matrix over all n locations and L is the Cholesky factor of the joint covariance matrix. 2) Simulation Studies: We evaluate the statistical accuracy of the double-precision (DP) and adaptive mixedprecision confidence region implementations using the synthetic data-generation procedure in Section VI-B. We use a stationary power-exponential model with variance parameter σ 2 = 1, smoothness parameter ν = 0.5, and settings β ∈ {0.033, 0.1, 0.234} for weak, medium, and strong spatial correlation, respectively. Each replicate combines no = 3,200 observation locations with a 120 × 120 prediction grid, giving np = 14,400 prediction locations and a total problem size of n = 17,600. We set the threshold to u = 10. The generated field combines spatially correlated residual variation with the increasing linear trend defined in Section VI-B. Under the deterministic trend, the u = 10 contour is y = 1 − 1.5x. The Gaussian residual field perturbs the trend-induced boundary. Figure 2 shows representative confidence regions at fixedacc = 10−4 across the three correlation regimes.

confidence region procedure uses nsim = 80,000 conditional simulations to estimate the critical values of the lower and upper test statistics, with nominal confidence level 1−α = 0.9. We compare the DP implementation with the adaptive mixedprecision implementation across a range of fixedacc values. For each run, we compute the true positive excursion set on the prediction grid, A+ u (Z)∩Dp , from the generated full field. A replicate is successful when the outer confidence region contains this set and the inner confidence region is contained within it Cin ⊆ A+ u (Z) ∩ Dp ⊆ Cout . Figure 3 reports empirical coverage rates across the common set of valid seeds, comparing DP and adaptive mixed precision with the nominal confidence level.

Inner coverage rate

B. Synthetic Data Generation and Simulation Studies

0

0.2

0.4

0.6

0.8

1

(c) Strong (1 : 0.234 : 0.5)

Fig. 2: Confidence regions for excursion sets using adaptive mixedprecision. Orange and blue indicate the inner and outer confidence regions. Black dots represent observations exceeding the threshold.

For each correlation setting, we repeat the experiment 1,000 times using independent random seeds. In each replicate, the

For the real-data evaluation, we use surface salinity from the HYCOM + NCODA Gulf of Mexico 1/25◦ Analysis product, GOMu0.04 experiment 90.1m000 [26] (Figure 4a). Produced by the Naval Research Laboratory, the product provides hourly fields on a uniform longitude–latitude grid. We convert the native HYCOM output to NetCDF and interpolate it to 40 standard z-levels. Available variables include sea-surface elevation, temperature, salinity, and the two horizontal velocity components. We use sea_water_salinity at the surface level (0 m) and model it as a Gaussian spatial process with Matérn covariance. The resulting confidence regions quantify uncertainty in low-salinity threshold exceedances across the domain. HYCOM uses a hybrid vertical coordinate, and its

Salinity [psu]

Latitude

28 26 24 22 20 97.5

92.5

87.5 Longitude

82.5

77.5

(a) Surface salinity field. The black square marks the selected study region.

27.5 Latitude

30

36.75 36.50 36.25 36.00 35.75 35.50 35.25 35.00

26.1

6

6

11

11

11

16 21

16 21

26

26

31 34

31 34

1

6

11 16 21 26 31 34

Tile column

1

1

double-precision

23.1 -92.5 -90.1 Longitude

-87.6

16 21 26

(a) Weak

24.6

21.6 -95.0

1

6

Tile row

29.0

1

Tile row

Tile row

analyses use the Navy Coupled Ocean Data Assimilation (NCODA) system [19].

6

11 16 21 26 31 34

Tile column

(b) Medium single-precision

31 34

1

6

11 16 21 26 31 34

Tile column

(c) Strong half-precision

Fig. 5: Lower-triangular tile-wise precision assignments resulting from adaptive mixed-precision Cholesky factorization under three spatial-correlation levels.

(b) Confidence regions for low-salinity excursions.

Fig. 4: Surface salinity and confidence regions for the HYCOM dataset. In (b), blue and orange indicate outer and inner regions; black points mark observations below 36 psu.

From the full surface-salinity field, we select the square subdomain λ ∈ [−95.00◦ , −87.64◦ ] and ϕ ∈ [21.64◦ , 29.00◦ ], where λ and ϕ denote longitude and latitude. The native grid supplies 185 × 185 = 34,225 valid surface values in this window. Both coordinate spans are 7.36◦ . Let Z(s) denote the original salinity and let Z̄ = 36.2215 psu be its mean over the selected grid. Since the implementation constructs positive excursion sets, we model the centered, sign-reversed response Y (s) = Z̄ − Z(s). The target is the lower-salinity region Z(s) ≤ 36 psu. In the transformed field, this corresponds to the upper threshold uY = Z̄ − 36 = 0.2215 psu. The implementation’s strict comparison Y > uY represents Z < 36, with equality having probability zero under the continuous Gaussian-process model. We discretize the normalized prediction domain on a 200×200 grid and model the transformed field with a Matérn model. The selected window spans the western and central Gulf of Mexico, from offshore Louisiana and Mississippi to the Mexican margin and Yucatán–Campeche region. In the north, Mississippi and Atchafalaya River inputs create low-salinity surface plumes, sharp horizontal gradients, and near-coastal variability. In the south, freshwater enters through groundwater discharge from the peninsula’s low-lying carbonate karst rather than surface rivers. The domain thus combines coastal and open-Gulf waters, broad mesoscale variability, and sharper salinity gradients near freshwater sources and coastlines. Figure 4b shows the resulting confidence regions. Both are concentrated near the northern boundary, where the field has its strongest coastal gradient. The inner region is narrow and lies along the upper boundary. The outer region is wider and includes additional scattered locations in the northern and northwestern parts of the domain. This pattern is consistent with the low-salinity surface plume near the northern Coast. D. TRMM Performance Analysis Figure 5 shows the precision assigned to the covariancematrix tiles before Cholesky factorization for the three spatialcorrelation levels. The colors represent the selected arithmetic

precision. Changes in the covariance structure alter the tile norms the adaptive precision rule uses. Each map shows the lower triangle of a 34 × 34 tile grid, comprising 595 tiles. In automatic mode, we selected half-precision for 36 tiles (6.05%) under weak correlation (β = 0.033), 66 tiles (11.09%) under medium correlation (β = 0.1), and 15 tiles (2.52%) under strong correlation (β = 0.234). The share therefore peaked in the medium case rather than increasing monotonically with correlation strength. The assignment rule compares the norm of each tile with the global matrix norm, while validation can adjust the fixedacc between cases. The observed shares reflect both the tile-norm distributions and the final accuracy settings, so they cannot be attributed to correlation strength alone. We evaluate single-node TRMM throughput on Polaris, Guyot, and JUPITER, followed by weak- and strong-scaling experiments on Frontier. Figure 6 presents the single-node TRMM performance as a function of matrix size. The experiments show increasing computational throughput with problem size, benefiting from reduced precision arithmetic and GPU matrix acceleration. Across the three correlation regimes, the maximum performance ratios relative to the theoretical FP64 peak reach 6.36×–6.93× on Polaris, 7.00×–7.24× on Guyot, and 10.20×–11.03× on JUPITER. These ratios compare mixedprecision throughput to the theoretical FP64 hardware limit. Figure 7 presents distributed-memory weak-scaling results for mixed-precision TRMM on Frontier using 1, 4, and 16 nodes under weak, medium, and strong spatial correlation. The global problem size grows with node count, with singlenode matrix sizes of N = 32,768, 65,536, and 98,304. From 1 to 16 nodes, aggregate throughput increases by 9.1– 10.0×, 10.8–11.5×, and 10.0–11.0× under weak, medium, and strong correlation, respectively. These results demonstrate substantial, though sublinear, throughput scaling for largescale conditional Gaussian simulation. Figure 8 presents strong-scaling results for mixed-precision TRMM on Frontier using 4, 8, and 16 nodes, with fixed matrix sizes of N = 163,840, 196,608, and 262,144. Increasing the node count from 4 to 16 improves throughput by 1.6–1.7×, 1.7–1.8×, and 1.7–1.9× under weak, medium, and strong correlation, respectively. Larger matrices generally achieve

200 100 0

Max Speedup: Weak 7.00× | Medium 7.24× | Strong 7.14× 1200 Weak 1000 Medium Strong 800 600 400 200 0

4 8 2 6 0 4 8 2 1638 3276 4915 6553 8192 9830 11468 13107 Matrix Size

8 3276

(a) 1 node on Polaris

6 6553

Max Speedup: Weak 10.20× | Medium 11.03× | Strong 10.93× 1500 Performance (Tflop/s)

300

Performance (Tflop/s)

Performance (Tflop/s)

Max Speedup: Weak 6.36× | Medium 6.93× | Strong 6.81× 600 Weak 500 Medium Strong 400

Weak Medium Strong

1250 1000 750 500 250 0

4 72 40 08 9830 1310 1638 1966 Matrix Size

8 3276

(b) 1 node on Guyot

6 6553

4 72 08 44 9830 1310 1966 2621 Matrix Size

(c) 1 node on Jupiter

N: 32768 N: 65536 N: 98304

10.0x 9.4x 9.1x

100

1

4 Number of Nodes

N: 32768 N: 65536 N: 98304

1000

11.4x 10.8x 11.5x

100

16

1

(a) Weak correlation

4 Number of Nodes

Performance (Tflop/s)

1000

Performance (Tflop/s)

Performance (Tflop/s)

Fig. 6: Shared-memory performance under different correlation levels. The dashed line indicates the theoretical FP64 peak performance. For each correlation setting, the maximum speedup relative to the theoretical FP64 peak is annotated at the top of the corresponding panel. N: 32768 N: 65536 N: 98304

1000

11.0x 10.6x 10.0x

100

16

1

(b) Medium correlation

4 Number of Nodes

16

(c) Strong correlation

N: 262144 1.7x

500

1.7x 1.6x

300 4

8 Number of Nodes

16

(a) Weak correlation

N: 163840 N: 196608

1000 800

N: 262144 1.8x 1.7x

600

1.7x 400 4

8 Number of Nodes

(b) Medium correlation

16

Performance (Tflop/s)

N: 163840 N: 196608

700

Performance (Tflop/s)

Performance (Tflop/s)

Fig. 7: Mixed-precision TRMM weak-scaling performance on Frontier under different correlation levels. Here, N denotes the matrix size on a single node, while the global problem size increases with the number of nodes according to the weak-scaling configuration. N: 163840 N: 196608

1000 800

N: 262144 1.9x 1.7x

600

1.7x

400 4

8 Number of Nodes

16

(c) Strong correlation

Fig. 8: Mixed-precision TRMM strong-scaling performance on Frontier under different correlation levels. Here, N denotes the matrix size across different node configurations.

higher throughput and greater scaling gains, although the improvements remain below the ideal fourfold increase. E. Application-Level Time to Solution We measure full application pipeline execution time on JUPITER using 4, 8, and 16 nodes, with four GH200 GPUs per node. Each configuration uses no = 10,240 observations. Prediction grids of 448 × 448, 512 × 512, and 640 × 640 give joint covariance dimensions of n = 210,944, 272,384, and 419,840, respectively. The corresponding simulation counts are 204,800, 327,680, and 417,792. Both the covariance dimension and simulation count increase with the available resources. For a consistent comparison across correlation levels, each configuration directly uses the final precision tolerance obtained from a preceding automatic-selection execution. Figure 9 shows that execution times are similar across correlation levels. Under weak correlation, the median times are approximately 80.04, 107.21, and 146.88 seconds on 4, 8, and 16 nodes. The corresponding medium correlation times are 80.94, 90.70, and 136.98 seconds. Strong correlation requires

approximately 89.13, 106.12, and 159.59 seconds. Thus, the larger configurations accommodate substantially larger covariance matrices and simulation ensembles, although elapsed time increases with node count. Because the workload dimensions vary between configurations. VII. C ONCLUSIONS We presented a distributed, GPU-accelerated framework for geostatistical excursion set confidence region applications that integrates kriging, uncertainty estimation, conditional simulation, and calibration through adaptive, tile-based mixed precision and PaRSEC. A parallel mixed-precision TRMM implementation, automatic precision-tolerance refinement, and covariance and operand scaling enable efficient reduced precision execution across covariance structures. Workflow optimizations further reduce memory and computational costs. Synthetic experiments and ocean surface-salinity results closely match double-precision while preserving empirical coverage near the nominal confidence level. Performance evaluations demonstrate high throughput across NVIDIA and AMD GPUs

Time (seconds)

160 140 120 100 80 60 40 20 0

159.587 s 146.881 s 136.983 s

Correlation

Weak Medium Strong

89.130 s 80.943 s 80.041 s

210,944 4 nodes

107.213 s 106.120 s 90.698 s

272,384 8 nodes

N

419,840 16 nodes

Fig. 9: Application execution time on JUPITER for different spatial correlation levels using 4, 8, and 16 nodes (four GH200 GPUs/node).

and increasing aggregate throughput across multiple nodes. These findings support adaptive mixed-precision as an effective approach to accelerating large-scale excursion set analysis while maintaining statistical accuracy. R EFERENCES [1] S. Abdulah, A. H. Baker, G. Bosilca, Q. Cao, S. Castruccio, M. G. Genton, D. E. Keyes, Z. Khalid, H. Ltaief, Y. Song et al., “Boosting earth system model outputs and saving petabytes in their storage using exascale climate emulators,” in SC24: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2024, pp. 1–12. [2] S. Abdulah, Q. Cao, Y. Pei, G. Bosilca, J. Dongarra, M. G. Genton, D. E. Keyes, H. Ltaief, and Y. Sun, “Accelerating geostatistical modeling and prediction with mixed-precision computations: A high-productivity approach with parsec,” IEEE Transactions on Parallel and Distributed Systems, vol. 33, no. 4, pp. 964–976, 2021. [3] S. Abdulah, H. Ltaief, Y. Sun, M. G. Genton, and D. E. Keyes, “Geostatistical modeling and prediction using mixed precision tile cholesky factorization,” in 2019 IEEE 26th international conference on high performance computing, data, and analytics (HiPC). IEEE, 2019, pp. 152–162. [4] T. Ahola, K. Virrantaus, J. M. Krisp, and G. J. Hunter, “A spatiotemporal population model to support risk assessment and damage analysis for decision-making,” International Journal of Geographical Information Science, vol. 21, no. 8, pp. 935–953, 2007. [5] K. Akbudak, H. Ltaief, A. Mikhalev, A. Charara, A. Esposito, and D. Keyes, “Exploiting data sparsity for large-scale matrix computations,” in European Conference on Parallel Processing. Springer, 2018, pp. 721–734. [6] C. Augonnet, S. Thibault, R. Namyst, and P.-A. Wacrenier, “Starpu: a unified platform for task scheduling on heterogeneous multicore architectures,” in European Conference on Parallel Processing. Springer, 2009, pp. 863–874. [7] M. Baboulin, S. Donfack, O. Kaya, T. Mary, and M. Robeyns, “Mixed precision randomized low-rank approximation with gpu tensor cores,” in European Conference on Parallel Processing. Springer, 2024, pp. 31–44. [8] D. Bolin and F. Lindgren, “Excursion and contour uncertainty regions for latent gaussian models,” Journal of the Royal Statistical Society Series B: Statistical Methodology, vol. 77, no. 1, pp. 85–106, 2015. [9] G. Bosilca, A. Bouteiller, A. Danalis, M. Faverge, T. Hérault, and J. J. Dongarra, “PaRSEC: Exploiting heterogeneity to enhance scalability,” Computing in Science & Engineering, vol. 15, no. 6, pp. 36–45, 2013. [10] A. Bouteiller, T. Herault, Q. Cao, J. Schuchart, and G. Bosilca, “Parsec: scalability, flexibility, and hybrid architecture support for task-based applications in ecp,” The International Journal of High Performance Computing Applications, vol. 39, no. 1, pp. 147–166, 2025. [11] A. Buttari, J. Langou, J. Kurzak, and J. Dongarra, “A class of parallel tiled linear algebra algorithms for multicore architectures,” Parallel computing, vol. 35, no. 1, pp. 38–53, 2009.

[12] Q. Cao, S. Abdulah, R. Alomairy, Y. Pei, G. Bosilca, J. Dongarra, M. G. Genton, D. E. Keyes, and H. Ltaief, “Reducing data motion and energy consumption of geospatial modeling applications using automated precision conversion,” in 2023 IEEE International Conference on Cluster Computing (CLUSTER). IEEE, 2023. [13] Q. Cao, S. Abdulah, R. Alomairy, Y. Pei, P. Nag, G. Bosilca, J. Dongarra, M. G. Genton, D. E. Keyes, H. Ltaief et al., “Reshaping geostatistical modeling and prediction for extreme-scale environmental applications,” in SC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2022, pp. 1–12. [14] Q. Cao, S. Abdulah et al., “Reducing data motion and energy consumption of geospatial modeling applications using automated precision conversion,” IEEE International Conference on Cluster Computing, 2023. [15] Q. Cao, R. Alomairy, Y. Pei, G. Bosilca, H. Ltaief, D. Keyes, and J. Dongarra, “A framework to exploit data sparsity in tile low-rank Cholesky factorization,” in 2022 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2022, pp. 414–424. [16] Q. Cao, Y. Pei, K. Akbudak, G. Bosilca, H. Ltaief, D. Keyes, and J. Dongarra, “Leveraging PaRSEC runtime support to tackle challenging 3d data-sparse matrix problems,” in 2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2021, pp. 79– 89. [17] Q. Cao, Y. Pei, K. Akbudak, A. Mikhalev, G. Bosilca, H. Ltaief, D. Keyes, and J. Dongarra, “Extreme-scale task-based cholesky factorization toward climate and weather prediction applications,” in Proceedings of the Platform for Advanced Scientific Computing Conference, 2020, pp. 1–11. [18] P. F. Craigmile and N. Cressie, “Spatial statistics and environmental science,” Technical Report, 2005. [19] J. A. Cummings and O. M. Smedstad, “Variational data assimilation for the global ocean,” in Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications, S. K. Park and L. Xu, Eds. Berlin, Heidelberg: Springer, 2013, vol. II, pp. 303–343. [20] A. Danalis, G. Bosilca, A. Bouteiller, T. Herault, and J. Dongarra, “PTG: An Abstraction for Unhindered Parallelism,” Proceedings of WOLFHPC: 4th International Workshop on Domain-Specific Languages and High-Level Frameworks for High Performance Computing, 2014. [21] J. P. French, “Confidence regions for the level curves of spatial data,” Environmetrics, vol. 25, no. 7, pp. 498–512, 2014. [22] ——, “Spatio-temporal exceedance locations and confidence regions,” The Annals of Applied Statistics, vol. 8, no. 3, pp. 1294–1314, 2014. [23] J. P. French and J. A. Hoeting, “Credible regions for exceedance sets of geostatistical data,” Environmetrics, vol. 27, no. 1, pp. 4–14, 2016. [24] J. P. French and S. R. Sain, “Confidence regions for the level curves of spatial data,” Environmetrics, vol. 24, no. 1, pp. 61–73, 2013. [25] N. J. Higham and T. Mary, “Mixed precision algorithms in numerical linear algebra,” Acta Numerica, vol. 31, pp. 347–414, 2022. [26] HYCOM Consortium, “HYCOM + NCODA Gulf of Mexico 1/25◦ Analysis (GOMu0.04/expt 90.1m000),” https://www.hycom.org/data/ gomu0pt04/expt-90pt1m000, accessed: Oct. 6, 2026. [27] S. K. Jha, G. Mariethoz, J. P. Evans, and M. F. McCabe, “Demonstration of a geostatistical approach to physically consistent downscaling of climate modeling simulations,” Water Resources Research, vol. 49, no. 1, pp. 245–259, 2013. [28] G. Lindgren and I. Rychlik, “Stationary stochastic processes and their representations,” International Statistical Review, vol. 63, no. 3, pp. 261– 272, 1995. [29] H. Ltaief, R. Alomairy, Q. Cao, J. Ren, L. Slim, T. Kurth, B. Dorschner, S. Bougouffa, R. Abdelkhalak, and D. E. Keyes, “Toward capturing genetic epistasis from multivariate genome-wide association studies using mixed-precision kernel ridge regression,” in SC24: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2024, pp. 1–12. [30] S. Markidis, S. W. Der Chien, E. Laure, I. B. Peng, and J. S. Vetter, “Nvidia tensor core programmability, performance & precision,” in 2018 IEEE international parallel and distributed processing symposium workshops (IPDPSW). IEEE, 2018, pp. 522–531. [31] Q. Pan, S. Abdulah, M. G. Genton, D. E. Keyes, H. Ltaief, and Y. Sun, “Gpu-accelerated vecchia approximations of gaussian processes for geospatial data using batched matrix computations,” in ISC High performance 2024 research paper proceedings (39th international conference). Prometeus GmbH, 2024, pp. 1–12.

[32] J. Ren, H. Ltaief, S. Abdulah, and D. E. Keyes, “Accelerating mixedprecision out-of-core cholesky factorization with static task scheduling,” in ISC High Performance 2025 Research Paper Proceedings (40th International Conference). Prometeus GmbH, 2025, pp. 1–12. [33] M. L. O. Salvaña, S. Abdulah, H. Ltaief, Y. Sun, M. G. Genton, and D. E. Keyes, “Parallel space-time likelihood optimization for air pollution prediction on large-scale systems,” in Proceedings of the Platform for Advanced Scientific Computing Conference, 2022, pp. 1–11. [34] G. Schieffer, D. A. De Medeiros, J. Faj, A. Marathe, and I. Peng, “On the rise of amd matrix cores: Performance, power efficiency, and programmability,” in 2024 IEEE International Symposium on Performance Analysis of Systems and Software (ISPASS). IEEE, 2024, pp. 132–143. [35] M. Sommerfeld, S. Sain, and A. Schwartzman, “Confidence regions for spatial excursion sets from repeated random field observations, with an application to climate,” Journal of the American Statistical Association, vol. 113, no. 523, pp. 1327–1340, 2018.

[36] I. Yamazaki, S. Tomov, and J. Dongarra, “Mixed-precision cholesky qr factorization and its case studies on multicore cpu with multiple gpus,” SIAM Journal on Scientific Computing, vol. 37, no. 3, pp. C307–C330, 2015. [37] H. Zhang, N. Cressie, and P. F. Craigmile, “Loss function approaches to predict a spatial quantile and its exceedance region,” Technometrics, vol. 50, no. 2, pp. 216–227, 2008. [38] Q. Zhang, R. Alomairy, D. Wang, Z. Gu, and Q. Cao, “Highperformance mixed-precision matrix multiplication via tile-centric design on modern architectures,” SN Computer Science, vol. 7, no. 1, p. 24, 2026. [39] X. Zhang, S. Abdulah, J. Cao, H. Ltaief, Y. Sun, M. G. Genton, and D. E. Keyes, “Parallel approximations for high-dimensional multivariate normal probability computation in confidence region detection applications,” in 2024 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2024, pp. 265–276.

Related documents

• Mixed-Precision GPU Acceleration for Confidence Regions of Spatial Excursion Sets#1389985
• LLMQ: Efficient lower-precision LLM training for consumer GPUs#254429
• OmniPilot: An Uncertainty-Aware LLM Inference Advisor for Heterogeneous GPU Clusters#332499
• Do Deep Ensembles Actually Capture Uncertainty in Graph Neural Networks?#216853
• Sensitivity Fidelity Reveals Physical Dependencies Beyond Forecast Accuracy in Deep-Learning Earth-System Models#1326790