Low-rank regularization is an effective technique for addressing ill-posed inverse problems when the solution tends to be rank-deficient. Global low-rank assumptions do not always hold for seismic wavefields; however, local low-rank features are observed in many practical situations. To leverage this feature, the unknown parameter is partitioned into tiles, each represented via low-rank factorization. This framework is applied to regularize multidimensional deconvolution in the frequency domain, considering two key factors. First, the unknown variable, referred to as Green's function, must maintain symmetry according to the reciprocity principle of wave propagation. To ensure symmetry within the tile-based low-rank framework, diagonal tiles are represented as the product of a low-rank factor and its transpose when numerically rank-deficient. Otherwise, they are represented by preconditioned dense forms. Symmetry in off-diagonal is achieved by parameterizing subdiagonal tiles as the product of two distinct low-rank factors, with the corresponding superdiagonal tiles assigned as their transposes. Second, the local rank of Green's function is allowed to vary with frequency. The numerical rank and optimal tile size at each frequency are estimated based on the solution from a benchmark solver. This approach is referred to as reciprocity-aware, factorization-based adaptive tile low-rank regularization. To enhance the applicability of the local low-rank approximation, receivers are reordered using a Hilbert space-filling curve, ensuring that receivers close in physical space are close in index in the unknown Green's function matrix. Our method is applied to a large-scale 3D deconvolution problem, demonstrating improved accuracy over the benchmark solver, especially in scenarios with added random noise and sparse source sampling.

Full waveform inversion (FWI) can recover high-resolution subsurface velocity models. Conventional waveform-difference objectives, however, are vulnerable to cycle skipping when the starting model is inaccurate. We introduce an FWI objective that compares features produced from modeled and observed seismic traces by SeisLM, a pretrained seismic foundation model. The SeisLM encoder remains frozen during inversion, and the feature discrepancy is differentiated with respect to the modeled traces to construct an adjoint source compatible with the standard adjoint-state framework. We also test a scheduled hybrid loss that combines the SeisLM feature loss with the conventional L_2 objective. Time-shift diagnostics show that the loss computed from features produced by the pretrained encoder has a broader and smoother basin around the correct alignment than either the waveform L_2 objective or the feature loss obtained from an encoder with the same architecture and randomly initialized parameters. In the Marmousi experiment, the SeisLM and hybrid objectives produce similar improvements during early-stage inversion and provide useful models for subsequent reflection-based L_2 refinement. In the 2D Overthrust experiment, which begins from a laterally invariant linear-gradient model, the SeisLM feature-loss workflow outperforms the conventional and hybrid workflows, indicating that introducing the L_2 contribution too early can reintroduce cycle-skipping sensitivity. In the 3D Overthrust experiment, conventional L_2 inversion stalls near the initial linear gradient, whereas the SeisLM feature loss guides the inversion toward a background model from which L_2 refinement recovers the principal structures. These results support using features produced by pretrained seismic networks to define early-stage FWI objectives rather than complete replacements for waveform-domain misfits.

Full waveform inversion (FWI) is a powerful technique for estimating high-resolution subsurface velocity models by minimizing the discrepancy between modelled and observed seismic data. However, the oscillatory nature of seismic waveforms makes point-wise discrepancy measures highly prone to cycle skipping, especially when the initial velocity model is inadequate. To address this challenge, various alternative misfit functions have been proposed in the literature, each with unique strengths and limitations. Dynamic time warping (DTW) is a popular technique in signal processing for aligning time series using dynamic programming. While a differentiable variant of DTW has been recently proposed, its use in FWI is hindered by high-frequency artefacts in the adjoint source and the substantial computational cost of gradient evaluations. In this study, we propose a neural network-based approach to learn the time shifts that align two time series in a supervised manner. The trained network is then utilized to compare traces from observed and modelled seismic data, offering a stable and computationally efficient alternative to DTW. Furthermore, the inherent differentiability of neural networks via backpropagation enables seamless integration into the FWI framework as a misfit function. We validate this approach on two synthetic datasets, namely the Marmousi model and the Chevron blind test dataset, demonstrating in both cases a similar convergence behaviour to that of Soft-DWT whilst drastically reducing the computational time of the adjoint source calculation.

This paper presents a new fast iterative solver for large systems involving kernel matrices. Advantageous aspects of 2 matrix approximations and the multigrid method are hybridized to create the 2-MG algorithm. This combination provides the time and memory efficiency of 2 operator representation along with the rapid convergence of a multilevel method. We describe how 2-MG works, show its linear complexity, and demonstrate its effectiveness on two standard kernels and on a single-layer potential boundary element discretization with complex geometry. The current zoo of 2 solvers, which includes a wide variety of iterative and direct solvers, so far lacks a method that exploits multiple levels of resolution, commonly referred to in the iterative methods literature as ``multigrid"" from its origins in a hierarchy of grids used to discretize differential equations. This makes 2-MG a valuable addition to the collection of 2 solvers. The algorithm has the potential for advancing various fields that require the solution of large, dense, symmetric positive definite matrices.

We present factorization and solution phases for a new linear complexity direct solver designed for concurrent batch operations on fine-grained parallel architectures, for matrices amenable to hierarchical representation. We focus on the strong-admissibility-based $\mathscr{H}<^>{2}$ format, where strong recursive skeletonization factorization compresses remote interactions. We build upon previous implementations of $\mathscr{H}<^>{2}$ matrix construction for efficient factorization and solution algorithm design, which are illustrated graphically in stepwise detail. The algorithms are 'blackbox' in the sense that the only inputs are the matrix and right-hand side, without analytical or geometrical information about the origin of the system. We demonstrate linear complexity scaling in both time and memory on four representative families of dense matrices up to one million in size. Parallel scaling up to 16 threads is enabled by a multi-level matrix graph coloring and avoidance of dynamic memory allocations thanks to prefix-sum memory management. An experimental backward error analysis is included. We break down the timings of different phases, identify phases that are memory-bandwidth limited, and discuss alternatives for phases that may be sensitive to the trend to employ lower precisions for performance.

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.

Long-context inference in large language models is increasingly bottlenecked by the memory and compute cost of the KV-Cache. Low-rank factorization compresses KV projections by writing W ≈ A * B, where A produces latent KV states and B can be absorbed into downstream weights. In modern RoPE-based LLMs, this absorption fails: RoPE forces latent KV states to be reconstructed to full dimension, reintroducing substantial memory and compute overhead. We propose RoPE-Aligned Pruning (RAP), which prunes entire RoPE-aligned column pairs to preserve RoPE's 2x2 rotation structure, restore B absorption, and eliminate reconstruction. Our evaluation on LLaMA-3-8B and Mistral-7B shows that RAP enables joint reduction of KV-Cache, attention parameters, and FLOPs by 20-30

Abstract Reconstructing the structural geology and mineral composition of the first few kilometers of the Earth's subsurface from sparse or indirect surface observations remains a long‐standing challenge with critical applications in mineral exploration, geohazard assessment, and geotechnical engineering. This inherently ill‐posed problem is often addressed by classical geophysical inversion methods, which typically yield a single maximum‐likelihood model that fails to capture the full range of plausible geology. The adoption of modern deep learning methods has been limited by the lack of large 3D training data sets. We address this gap with StructuralGeo, a geological simulation engine that mimics eons of tectonic, magmatic, and sedimentary processes to generate a virtually limitless supply of realistic synthetic 3D lithological models. Using this data set, we train both unconditional and conditional generative flow‐matching models with a 3D attention U‐Net architecture. The resulting foundation model can reconstruct multiple plausible 3D scenarios from surface topography and sparse borehole data, depicting structures such as layers, faults, folds, and dikes. By sampling many reconstructions from the same observations, we introduce a probabilistic framework for estimating the size and extent of subsurface features. While the realism of the output is bounded by the fidelity of the training data to true geology, this combination of simulation and generative AI functions offers a flexible prior for probabilistic modeling, regional fine‐tuning, and use as an AI‐based regularizer in traditional geophysical inversion workflows.

This paper introduces sTiles, a GPU-accelerated framework for factorizing sparse structured symmetric matrices. By leveraging tile algorithms for fine-grained computations, sTiles uses a structure-aware task execution flow to handle challenging arrowhead sparse matrices with variable bandwidths, common in scientific and engineering fields. It minimizes fill-in during Cholesky factorization using permutation techniques and employs a static scheduler to manage tasks on shared-memory systems with GPU accelerators. sTiles balances tile size and parallelism, where larger tiles enhance algorithmic intensity but increase floating-point operations and memory usage, while parallelism is constrained by the arrowhead structure. To expose more parallelism, a left-looking Cholesky variant breaks sequential dependencies in trailing submatrix updates via tree reductions. Evaluations show sTiles achieves speedups of up to 8.41X, 9.34X, 5.07X, and 11.08X compared to CHOLMOD, SymPACK, MUMPS, and PARDISO, respectively, and a 5X speedup compared to a 32-core AMD EPYC CPU on an NVIDIA A100 GPU. Our generic software framework imports well-established concepts from dense matrix computations but they all require customizations in their deployments on hybrid architectures to best handle factorizations of sparse matrices with arrowhead structures.

This paper presents a novel factorization-based, low-rank regularization method for solving multidimensional deconvolution problems in the frequency domain. In this approach, each frequency component of the unknown wavefield is represented as a complex-valued square matrix and approximated using the product of one rectangular matrix and its transpose. The benefit of such a parametrization is two-fold: first, the size of the unknown matrix is greatly reduced compared to that of the original wavefield of interest (and halved compared to conventional factorization-based, low-rank approximations); second, the retrieved wavefield is implicitly guaranteed to comply with the reciprocity principle, as expected from theory. We further show that the proposed objective function can be successfully optimized using the accelerated proximal gradient algorithm and discuss a robust strategy to define the initial guess of the solution. Numerical examples on synthetic and field data demonstrate the effectiveness of the proposed method in compressing the retrieved Green's function while preserving its accuracy.

This paper explores the performance optimization of out-of-core (OOC) Cholesky factorization on shared-memory systems equipped with multiple GPUs. We employ fine-grained computational tasks to expose concurrency while creating opportunities to overlap data movement asynchronously with computations, especially when dealing with matrices that cannot fit on the GPU memory. We leverage the directed acyclic graph of the task-based Cholesky factorization and map it onto a static scheduler that promotes data reuse while supporting strategies for reducing data movement with the CPU host when the GPU memory is exhausted. The CPU-GPU interconnect may become the main performance bottleneck as the gap between the GPU execution rate and the traditional PCIe bandwidth continues to widen. While the surface-to-volume effect of compute-bound kernels partially mitigates the overhead of data motion, deploying mixed-precision (MxP) computations exacerbates the throughput discrepancy. Using static task scheduling, we evaluate the performance capabilities of the new ultra-fast NVIDIA chip interconnect technology, codenamed NVLink-C2C, that constitutes the backbone of the NVIDIA Grace Hopper Superchip (GH200), against a new four-precision (FP64/FP32/FP16/FP8) left-looking Cholesky factorization. We report the performance results of a benchmarking campaign on various NVIDIA GPU generations and interconnects. We highlight 20 cuSOLVER on a single GH200 with FP64 while hiding the cost of OOC task-based Cholesky factorization, and we scale almost linearly on four GH200 superships. With MxP enabled, our statically scheduled four-precision tile-based Cholesky factorization scores a 3X performance speedup against its FP64-only counterpart, delivering application-worthy FP64 accuracy when modeling a large-scale geospatial statistical application.

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.

We propose a novel method for solving linear systems with large dense operators, right-hand sides, and unknown matrices. We extend the idea of performing a low-rank or ℋ^2 approximation of the operator to simultaneously approximating of all three matrices, using a shared bases assumption. We showcase the effectiveness of our method through its application to a particularly challenging problem in the field of seismic imaging, specifically the use of Multidimensional Deconvolution (MDD) for redatuming applications. While offering more accuracy than conventional correlation-based redatuming methods, MDD faces challenges due to the ill-posed nature of the underlying inverse problem and the requirement to handle large, dense, complex-valued matrices. These obstacles have long limited the adoption of MDD in the geophysical community. Recent interest in this technology has spurred the development of new strategies to enhance the robustness of the inversion process and reduce its computational overhead. Our proposed approach can greatly alleviate the data-heavy nature of MDD. Moreover, since in 3d applications the matrices do not lend themselves to global low rank approximations, we introduce a novel ℋ^2 -like approximation. With this work, we aim to streamline MDD implementations, fostering efficiency and controlling accuracy in wavefield reconstruction. This innovation holds potential for broader applications in the geophysical domain and beyond.

Radial Basis Function (RBF) interpolation is a widely used technique for maintaining high-quality 3D unstructured meshes in fluid-structure interaction, ensuring smooth deformation while preserving mesh integrity. This paper studies the rank distribution of global and compact support RBF kernels using Tile Low-Rank (TLR) Cholesky solver, focusing on their impact on computational efficiency and scalability. We analyze the arithmetic intensity and performance of these RBFs through the hierarchical Roofline model to understand computational trade-offs. Large-scale experiments are conducted to evaluate the scalability of TLR Cholesky solver for two kernels: Gaussian (global support) and Wendland (compact support). As a real-world application, we examine the SARS-CoV-2 viral mesh, where RBF-based deformation is employed to model viral populations in dynamic fluid environments. We adopt a sphere-packing approach to arrange viruses in a dense, non-overlapping configuration, modeling real-world constraints. Our findings show that Gaussian's lower-rank intensity hinders scalability, whereas Wendland's compact support results in higher computational cost. The results demonstrate that the TLR approach reveals performance characteristics that dense factorization cannot, offering insights into kernel selection, scalability, and computational trade-offs for large-scale RBF-based mesh deformation.

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.

Selected inversion is essential for applications such as Bayesian inference, electronic structure calculations, and inverse covariance estimation, where computing only specific elements of large sparse matrix inverses significantly reduces computational and memory overhead. We present an efficient implementation of a two-phase parallel algorithm for computing selected elements of the inverse of a sparse symmetric matrix A, which can be expressed as A = LL^T through sparse Cholesky factorization. Our approach leverages a tile-based structure, focusing on selected dense tiles to optimize computational efficiency and parallelism. While the focus is on arrowhead matrices, the method can be extended to handle general structured matrices. Performance evaluations on a dual-socket 26-core Intel Xeon CPU server demonstrate that sTiles outperforms state-of-the-art direct solvers such as Panua-PARDISO, achieving up to 13X speedup on large-scale structured matrices. Additionally, our GPU implementation using an NVIDIA A100 GPU demonstrates substantial acceleration over its CPU counterpart, achieving up to 5X speedup for large, high-bandwidth matrices with high computational intensity. These results underscore the robustness and versatility of sTiles, validating its effectiveness across various densities and problem configurations.

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.

The Hamming distance, a fundamental measure of dissimilarity between data points, plays a crucial role in various fields, including error detection, machine learning, and genomic sequence alignment, where it is commonly used for identifying mismatches in nucleotide or protein sequences. This work introduces two implementations for computing Hamming distances for sequence alignment: synchronous and asynchronous matrix-based approaches. While most existing implementations rely on vector-based methods due to their simplicity and ease of use, they are not efficient for large-scale data. Our work focuses on enhancing performance by introducing matrix-based implementations that significantly improve computational efficiency and scalability. Our asynchronous implementation showcases Julia for sequential task flow and PaRSEC for parameterized task graph execution models on homogeneous and heterogeneous architectures. CPU computations use INT8 GEMM from oneMKL, while GPU implementations employ Tensor/Matrix Core INT8 GEMM from cuBLAS/hipBLAS and 1-bit TensorOps GEMM capabilities from CUTLASS. For constructing bitmask matrices on GPUs, we develop both a naive CUDA implementation using global memory and an optimized implementation utilizing shared memory at the warp level, with the optimized version achieving a 5X speedup over the naive approach. The results demonstrate significant performance improvements, with the asynchronous matrix-based implementation achieving up to 284X speedup over the vector-based approach on CPUs, while the asynchronous GPU-enabled implementation on A100 GPUs delivers a 15X speedup compared to the CPU matrix-based approach and a three orders of magnitude improvement over the CPU vector-based approach. Furthermore, the asynchronous implementation of PaRSEC scales well on up to 256 nodes of Summit and Frontier. These advancements highlight the scalability and efficiency of matrix-based Hamming distance computation, leveraging GPU acceleration and advanced asynchronous execution, paving the way forward for large-scale genomic sequence alignment and data analysis.

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:

In this work, we explore the data sparsity characteristics of different Radial Basis Function (RBF) kernels in the context of 3D unstructured mesh deformation. While RBF interpolation is a powerful method for generating high-quality adaptive meshes, solving the resulting boundary problems leads to large, dense linear systems that are computationally expensive and memory-intensive due to their cubic computational complexity and quadratic memory requirement. To address these challenges, we exploit the rank structure of the matrix operators by employing a Tile Low-Rank (TLR) Cholesky-based solver, which approximates off-diagonal matrix tiles up to an application-specific accuracy threshold. Our study compares global support RBFs and compact support RBFs, focusing on their effects on rank distribution and numerical accuracy. Using realistic 3D geometries of SARS-CoV-2 viruses from the Protein Data Bank (PDB), we evaluate various RBF kernels, analyze the corresponding matrix rank structures, and assess the backward error resulting from low-rank approximations for different kernel types. We conduct experiments on various shared-memory systems, demonstrating the performance scalability on massively parallel architectures. Leveraging the Hierarchical Computations on Manycore Architectures (HiCMA) library and the PaRSEC runtime system, we show how data sparsity accelerates large-scale mesh adaptation, providing valuable insights into the balance between computational efficiency and numerical accuracy.