ABSTRACT Accurate wind speed nowcasting is crucial for optimizing wind energy production and grid stability, especially for countries such as Saudi Arabia with ambitious renewable energy targets. This article introduces a novel deep learning framework combining Deep Echo State Networks (DESNs) with spatial information through ‐nearest neighbors for high‐resolution wind speed prediction. Specifically, we develop a Spatio‐Temporal Attention Graph Autoencoder (STAGA) that effectively reduces spatial dimensionality while preserving critical spatio‐temporal patterns, enabling efficient processing of large‐scale meteorological data. Experiments using high‐resolution simulated wind data from Saudi Arabia demonstrate that our model consistently outperforms traditional methods. In addition, we provide uncertainty quantification through conformal prediction and demonstrate its practical value by assessing wind power estimates. The proposed methodology offers significant improvements for operational wind forecasting systems, supporting the efficient integration of wind energy into power grids.

In recent decades, statisticians have been increasingly encountering spatial data that exhibit non-Gaussian behaviors such as asymmetry and heavy-tailedness. As a result, the assumptions of symmetry and fixed tail weight in Gaussian processes have become restrictive and may fail to capture the intrinsic properties of the data. To address the limitations of the Gaussian models, a variety of skewed models has been proposed, of which the popularity has grown rapidly. These skewed models introduce parameters that govern skewness and tail weight. Among various proposals in the literature, unified skewed distributions, such as the Unified Skew-Normal (SUN), have received considerable attention. In this work, we revisit a more concise and intepretable re-parameterization of the SUN distribution and apply the distribution to random fields by constructing a generalized unified skew-normal (GSUN) spatial process. We demonstrate that the GSUN is a valid spatial process by showing its vanishing correlation at large distances and provide the corresponding spatial interpolation method. In addition, we develop an inference mechanism for the GSUN process using the concept of neural Bayes estimators with deep graphical attention networks (GATs) and encoder transformer. We show the superiority of our proposed estimator over conventional CNN-based architectures in terms of stability and accuracy by means of a simulation study. In addition, we demonstrate that the GSUN process offers enhanced flexibility compared to another model proposed in the literature through an application to Pb-contaminated soil data. Furthermore, we show that the GSUN process is different from the conventional Gaussian processes and Tukey g-and-h processes, through the probability integral transform (PIT).

Advancements in information technology have enabled the creation of massive spatial datasets, driving the need for scalable and efficient computational methodologies. While offering viable solutions, centralized frameworks are limited by vulnerabilities such as single-point failures and communication bottlenecks. This paper presents a decentralized framework tailored for parameter inference in spatial low-rank models to address these challenges. A key obstacle arises from the spatial dependence among observations, which prevents the log-likelihood from being expressed as a summation-a critical requirement for decentralized optimization approaches. To overcome this challenge, we propose a novel objective function leveraging the evidence lower bound, which facilitates the use of decentralized optimization techniques. Our approach employs a block descent method integrated with multi-consensus and dynamic consensus averaging for effective parameter optimization. We prove the convexity of the new objective function in the vicinity of the true parameters, ensuring the convergence of the proposed method. Additionally, we present the first theoretical results establishing the consistency and asymptotic normality of the estimator within the context of spatial low-rank models. Extensive simulations and real-world data experiments corroborate these theoretical findings, showcasing the robustness and scalability of the framework.

Gaussian Random Fields (GRFs) with Matérn covariance functions have emerged as a powerful framework for modeling spatial processes due to their flexibility in capturing different features of the spatial field. However, the smoothness parameter is challenging to estimate using maximum likelihood estimation (MLE), which involves evaluating the likelihood based on the full covariance matrix of the GRF, due to numerical instability. Moreover, MLE remains computationally prohibitive for large spatial datasets. To address this challenge, we propose the Fisher-BackTracking (Fisher-BT) method, which integrates the Fisher scoring algorithm with a backtracking line search strategy and adopts a series approximation for the modified Bessel function. This method enables an efficient MLE estimation for spatial datasets using the ExaGeoStat high-performance computing framework. Our proposed method not only reduces the number of iterations and accelerates convergence compared to derivative-free optimization methods but also improves the numerical stability of the smoothness parameter estimation. Through simulations and real-data analysis using a soil moisture dataset covering the Mississippi River Basin, we show that the proposed Fisher-BT method achieves accuracy comparable to existing approaches while significantly outperforming derivative-free algorithms such as BOBYQA and Nelder-Mead in terms of computational efficiency and numerical stability.

Asymptotic methods for hypothesis testing in high-dimensional data usually require the dimension of the observations to increase to infinity, often with an additional relationship between the dimension (say, p) and the sample size (say, n). On the other hand, multivariate asymptotic testing methods are valid for fixed dimension only and their implementations typically require the sample size to be large compared to the dimension to yield desirable results. In practical scenarios, it is usually not possible to determine whether the dimension of the data conform to the conditions required for the validity of the high-dimensional asymptotic methods for hypothesis testing, or whether the sample size is large enough compared to the dimension of the data. In this work, we first describe the notion of uniform-over-p convergences and subsequently, develop a uniform-over-dimension central limit theorem. An asymptotic test for the two-sample equality of locations is developed, which now holds uniformly over the dimension of the observations. Using simulated and real data, it is demonstrated that the proposed test exhibits better performance compared to several popular tests in the literature for high-dimensional data as well as the usual scaled two-sample tests for multivariate data, including the Hotelling's T^2 test for multivariate Gaussian data.

In this work, we introduce an integrated depth measure for functional data defined over complex multidimensional domains. We consider functional data whose discrete realizations are irregularly spaced, and may be available only over portions of the domain. To address this issue, we propose an integrated depth based on a Voronoi tessellation of the multidimensional domain. This approach ensures favorable statistical properties for the proposed depth, as well as computational efficiency, enabling the analysis of large-scale functional datasets. We validate our proposal with the study of air temperatures across the Earth surface, as provided by the CESM2 Large Ensemble Community Project. The proposed depth is able to capture the increase in global temperatures since the 1980s, coherently with global warming.

Voluminous highly resolved climate data requires high-performance statistical modeling frameworks that are both accurate and sustainable on modern supercomputing platforms. Spatio-temporal modeling methods, such as Maximum Likelihood Estimation (MLE), primarily rely on performing a dense Cholesky factorization of large covariance matrices, which becomes a significant bottleneck at scale. We present an enhanced version of ExaGeoStat, a high-performance geostatistical modeling framework, designed to address this challenge through tile-based matrix compression, mixed-precision arithmetic, and architecture-aware scheduling. Building upon our previous work on Fugaku, an Arm A64FX system, we extend our results to three supercomputers featuring diverse hardware architectures: Frontier (AMD MI250X GPUs), Alps (NVIDIA GH200 GPUs), and Shaheen III (AMD EPYC Genoa CPUs). We demonstrate that adaptive algorithms targeting low-rank and low-precision opportunities can deliver up to 4× speedup, 2× memory savings, and over 70% energy reduction while maintaining application-acceptable accuracy. Our work demonstrates the feasibility of running accurate climate-emulation workloads on next-generation AI hardware in a power-efficient manner, thereby advancing the sustainability of climate data science itself.

Wildfires pose an increasingly severe threat to air quality, yet quantifying their causal impact remains challenging due to unmeasured meteorological and geographic confounders. Moreover, wildfire impacts on air quality may exhibit heterogeneous effects across pollution levels, which conventional mean-based causal methods fail to capture. To address these challenges, we develop a Quantile-based Latent Spatial Confounder Model (QLSCM) that substitutes conditional expectations with conditional quantiles, enabling causal analysis across the entire outcome distribution. We establish the causal interpretation of QLSCM theoretically, prove the identifiability of causal effects, and demonstrate estimator consistency under mild conditions. Simulations confirm the bias correction capability and the advantage of quantile-based inference over mean-based approaches. Applying our method to contiguous US wildfire and air quality data, we uncover important heterogeneous effects: fire radiative power exerts significant positive causal effects on aerosol optical depth at high quantiles in Western states like California and Oregon, while insignificant at lower quantiles. This indicates that wildfire impacts on air quality primarily manifest during extreme pollution events. Regional analyses reveal that Western and Northwestern regions experience the strongest causal effects during such extremes. These findings provide critical insights for environmental policy by identifying where and when mitigation efforts would be most effective.

R has become a cornerstone of scientific and statistical computing due to its extensive package ecosystem, expressive syntax, and strong support for reproducible analysis. However, as data sizes and computational demands grow, native R parallelism support remains limited. This paper presents RCOMPSs, a scalable runtime system that enables efficient parallel execution of R applications on multicore and manycore systems. RCOMPSs adopts a dynamic, task-based programming model, allowing users to write code in a sequential style, while the runtime automatically handles asynchronous task execution, dependency tracking, and scheduling across available resources. We present RCOMPSs using three representative data analysis algorithms, i.e., K-nearest neighbors (KNN) classification, K-means clustering, and linear regression and evaluate their performance on two modern HPC systems: KAUST Shaheen-III and Barcelona Supercomputing Center (BSC) MareNostrum 5. Experimental results reveal that RCOMPSs demonstrates both strong and weak scalability on up to 128 cores per node and across 32 nodes. For KNN and K-means, parallel efficiency remains above 70 acceptable performance under shared and distributed memory configurations despite its deeper task dependencies. Overall, RCOMPSs significantly enhances the parallel capabilities of R with minimal, automated, and runtime-aware user intervention, making it a practical solution for large-scale data analytics in high-performance environments.

Parameter estimation with the maximum Lq-likelihood estimator (MLqE) is an alternative to the maximum likelihood estimator (MLE) that considers the q-th power of the likelihood values for some 0

Gaussian Processes (GPs) are vital for modeling and predicting irregularly-spaced, large geospatial datasets. However, their computations often pose significant challenges in large-scale applications. One popular method to approximate GPs is the Vecchia approximation, which approximates the full likelihood via a series of conditional probabilities. The classical Vecchia approximation uses univariate conditional distributions, which leads to redundant evaluations and memory burdens. To address this challenge, our study introduces block Vecchia, which evaluates each multivariate conditional distribution of a block of observations, with blocks formed using the K-means algorithm. The proposed GPU framework for the block Vecchia uses varying batched linear algebra operations to compute multivariate conditional distributions concurrently, notably diminishing the frequent likelihood evaluations. Diving into the factor affecting the accuracy of the block Vecchia, the neighbor selection criterion is investigated, where we found that the random ordering markedly enhances the approximated quality as the block count becomes large. To verify the scalability and efficiency of the algorithm, we conduct a series of numerical studies and simulations, demonstrating their practical utility and effectiveness compared to the exact GP. Moreover, we tackle large-scale real datasets using the block Vecchia method, i.e., high-resolution 3D profile wind speed with a million points.

Visualization and assessment of copula structures are crucial for accurately understanding and modeling the dependencies in multivariate data analysis. In this paper, we introduce an innovative method that employs functional boxplots and rank-based testing procedures to evaluate copula symmetry. This approach is specifically designed to assess key characteristics such as reflection symmetry, radial symmetry, and joint symmetry. We first construct test functions for each specific property and then investigate the asymptotic properties of their empirical estimators. We demonstrate that the functional boxplot of these sample test functions serves as an informative visualization tool of a given copula structure, effectively measuring the departure from zero of the test function. Furthermore, we introduce a nonparametric testing procedure to assess the significance of deviations from symmetry, ensuring the accuracy and reliability of our visualization method. Through extensive simulation studies involving various copula models, we demonstrate the effectiveness of our testing approach. Finally, we apply our visualization and testing techniques to two real-world datasets: a nutritional habits survey with five variables and wind speed data from three locations in Saudi Arabia.

Modified Bessel functions of the second kind are widely used in physics, engineering, spatial statistics, and machine learning. Since contemporary scientific applications, including machine learning, rely on GPUs for acceleration, providing robust GPU-hosted implementations of special functions, such as the modified Bessel function, is crucial for performance. Existing implementations of the modified Bessel function of the second kind rely on CPUs and have limited coverage of the full range of values needed in some applications. In this work, we present a robust implementation of the modified Bessel function of the second kind on GPUs, eliminating the dependence on the CPU host. We cover a range of values commonly used in real applications, providing high accuracy compared to common libraries like the GNU Scientific Library (GSL) when referenced to Mathematica as the authority. Our GPU-accelerated approach demonstrates a 2.68x performance improvement using a single A100 GPU compared to the GSL on 40-core Intel Cascade Lake CPUs. Our implementation is integrated into ExaGeoStat, the HPC framework for spatial data modeling, where the modified Bessel function of the second kind is required by the Matérn covariance function in generating covariance matrices. We accelerate the matrix generation process in ExaGeoStat by up to 12.62x with four A100 GPUs while maintaining almost the same accuracy for modeling and prediction operations using synthetic and real datasets.

We recognize the emergence of a statistical computing community focused on working with large computing platforms and producing software and applications that exemplify high-performance statistical computing (HPSC). The statistical computing (SC) community develops software that is widely used across disciplines. However, it remains largely absent from the high-performance computing (HPC) landscape, particularly on platforms such as those featured on the www.top500.org or Green500 lists. Many disciplines already participate in HPC, mostly centered around simulation science, although data-focused efforts under the artificial intelligence (AI) label are gaining popularity. Bridging this gap requires both community adaptation and technical innovation to align statistical methods with modern HPC technologies. We can accelerate progress in fast and scalable statistical applications by building strong connections between the SC and HPC communities. We present a brief history of SC, a vision for how its strengths can contribute to statistical science in the HPC environment (such as HPSC), the challenges that remain, and the opportunities currently available, culminating in a possible roadmap toward a thriving HPSC community. This article is categorized under:

Emulating computationally intensive scientific simulations is crucial for enabling uncertainty quantification, optimization, and informed decision-making at scale. Gaussian Processes (GPs) offer a flexible and data-efficient foundation for statistical emulation, but their poor scalability limits applicability to large datasets. We introduce the Scaled Block Vecchia (SBV) algorithm for distributed GPU-based systems. SBV integrates the Scaled Vecchia approach for anisotropic input scaling with the Block Vecchia (BV) method to reduce computational and memory complexity while leveraging GPU acceleration techniques for efficient linear algebra operations. To the best of our knowledge, this is the first distributed implementation of any Vecchia-based GP variant. Our implementation employs MPI for inter-node parallelism and the MAGMA library for GPU-accelerated batched matrix computations. We demonstrate the scalability and efficiency of the proposed algorithm through experiments on synthetic and real-world workloads, including a 50M point simulation from a respiratory disease model. SBV achieves near-linear scalability on up to 512 A100 and GH200 GPUs, handles 2.56B points, and reduces energy use relative to exact GP solvers, establishing SBV as a scalable and energy-efficient framework for emulating large-scale scientific models on GPU-based distributed systems.

We develop a new method based on the Lq-likelihood to perform robust covariance estimation on spatial data. We propose the maximum composite Lq-likelihood estimator (MCLqE), which is obtained from a combination of the maximum Lq-likelihood and the composite likelihood. In order to down-weight a single value at a spatial location, we make use of the composite likelihood of the spatial data, which is in the form of the product of the pairwise likelihood of all pairs of the observations. We apply the maximum Lq-likelihood on each pair to form the MCLqE, which effectively down-weights abnormal observations in a spatial dataset and gives more robust covariance parameter estimation results when the data are contaminated by outliers. We adapt the hyperparameter tuning procedure developed for the MLqE to the new MCLqE framework. The robustness of the MCLqE is further verified by both simulation studies and experiments on a precipitation dataset from the US.

We present a new paradigm, called functional multiple-point simulation, in which multiple-point geostatistical simulation can be performed when functions or curves are observed at each location of a random field. Multiple-point simulation is a non-parametric method used for conditional geostatistical simulation of complex spatial patterns by inferring multiple-point statistics from a training image, rather than from a two-point variogram or covariance model. When the observable at each spatial location is a functional random variable, such multiple-point simulation must take into account not only the spatial correlation among locations but also the similarity of functions or curves observed at each location. The data events to be compared in this case are now functional, in the sense that they consist of spatial arrangements of functions. Consequently, we propose four distances, inspired by the functional data analysis literature, for measuring similarities between functional data events and use these to extend the direct sampling method to perform multiple-function geostatistical simulation with functional fields. We coin the new method Functional Direct Sampling and carry out extensive qualitative and quantitative performance comparison between the four proposed distances using simulation techniques on two well-known applications of multiple-point simulation: simulating copies of a functional random field and gap-filling of locations in a functional random field. We apply the proposed method to a gap-filling task of simulated wind profiles spatial functions over the Arabian Peninsula.

In the past decades, clean and renewable energy has gained increasing attention due to a global effort on carbon footprint reduction. In particular, Saudi Arabia is gradually shifting its energy portfolio from an exclusive use of oil to a reliance on renewable energy, and, in particular, wind. Modelling wind for assessing potential energy output in a country as large, geographically diverse and understudied as Saudi Arabia is a challenge which implies highly non-linear dynamic structures in both space and time. To address this, we propose a spatio-temporal model whose spatial information is first reduced via an energy distance-based approach and then its dynamical behaviour is informed by a sparse and stochastic recurrent neural network (Echo State Network). Finally, the full spatial data is reconstructed by means of a non-stationary stochastic partial differential equation-based approach. Our model can capture the fine scale wind structure from a high-resolution simulation and produce more accurate forecasts of both wind speed and energy in lead times of interest for energy grid management.