We introduce the QuadratiK package that incorporates innovative data analysis methodologies. The presented software, implemented in both R and Python , offers a comprehensive set of goodness-of-fit tests and clustering techniques using kernel-based quadratic distances, thereby bridging the gap between the statistical and machine learning literature. Our software implements one, two, and \( k \) -sample tests for goodness-of-fit, providing an efficient and mathematically sound way to assess the fit of probability distributions. Expanded capabilities of our software include supporting tests for uniformity on the \( d \) -dimensional sphere based on Poisson kernel densities. Particularly noteworthy is the incorporation of a unique clustering algorithm specifically tailored for spherical data that leverages a mixture of Poisson kernel-based densities on the \( d \) -dimensional sphere. Alongside this, our software includes additional graphical functions, aiding the users in validating, as well as visualizing and representing clustering results. This enhances the interpretability and usability of the analysis. In summary, our R and Python packages provide a powerful suite of tools, offering researchers and practitioners the means to delve deeper into their data, draw robust inferences, and conduct potentially impactful analyses and inferences across a wide array of disciplines.
Hardware trends favor algorithm designs that maximize data reuse per FLOP. We develop and benchmark high-performance Multipole-to-Local (M2L) translation operators for the kernel-independent Fast Multipole Method (kiFMM), a widely adopted FMM variant that supports a broad class of kernels and has been favored by recent implementations for its simple specification. Naively implemented, M2L is bandwidth-limited and therefore a key bottleneck in the FMM. State-of-the-art FFT-based M2L implementations, though elegant and with a fast setup time, suffer from low operational intensity and require architecture-specific optimizations. We demonstrate that a Basic Linear Algebra Subprograms (BLAS)-based M2L, combined with randomized low-rank compression, achieves competitive performance with greater portability and a simpler implementation leveraging existing BLAS infrastructure, at the cost of higher setup times—especially for high-accuracy settings in double precision. Our Rust-based implementation enables seamless switching between strategies for fair benchmarking. Results on CPUs show that FFT-based M2L is favorable in low-accuracy settings or dynamic particle simulations, while BLAS-based M2L is favored for high-accuracy settings for static particle distributions, where its higher setup costs are amortized in many practical applications of the FMM.
This review evaluates the reproducibility of the results presented in “Exploiting (min, +)/(max, +) Isomorphism to Speed up Convolutions.” The focus is exclusively on the verification of the authors’ data, code, and experimental procedures. Using the provided materials, all reported results were successfully reproduced on an independent system, confirming the reliability and transparency of the computational workflow.
Time series distance measures are fundamental in numerous domains, including finance, healthcare, and signal processing, enabling crucial tasks such as pattern recognition, anomaly detection, and predictive modeling. However, many applications require computing distances between all pairs of time series in large datasets, a computationally intensive task that can become a significant bottleneck in analysis pipelines. The tsdistances library is a high-performance Python package designed for computing distances between time series, with GPU support for accelerated processing. This article introduces tsdistances and its key features, focusing on the implementation of elastic distance algorithms and their optimizations. We present both CPU and GPU implementations, highlighting the use of dynamic programming techniques and GPU-specific optimizations such as warp-based parallelization. The performance of tsdistances is compared with existing alternatives in the literature, demonstrating significant speed improvements, especially for large-scale time series analysis tasks.
A Lipschitz continuous function \( f \) that has a high-precision, albeit slow reference implementation can be evaluated more rapidly by piecewise polynomials, obtained by Chebyshev approximation. The partition into subdomains should be made such that the computation of subdomain index and reduced function argument is fast and introduces no rounding errors. A rigorous analysis of truncation and rounding errors shows that the maximum relative error can be chosen arbitrarily close to \(2\epsilon\) , except around zeroes of \( f \) . Alignment of coefficients minimizes cache loading times. Code generation and target algorithms are implemented in the open source project ppapp , maintained at https://jugit.fz-juelich.de/mlz/ppapp , with a snapshot archived in the Collected Algorithms of the ACM. It has been successfully applied to three real-valued functions in the open source complex error function library libcerf .
(min,+) and (max,+) algebra lie at the core of theories for the analysis of worst-case performance bounds, such as Deterministic Network Calculus and Real-Time Calculus. In both algebras, computing expressions can be computationally expensive. In particular, the convolution operation can be very time-consuming, sometimes taking hours or not completing at all. In fact, its operands are represented as pseudo-periodic sequences of segments and points (henceforth elements for short), which may have different periods. As already observed in literature, a convolution requires that every couple of elements belonging to different operands be elaborated, up to the least common multiple (lcm) of their periods. In this article, leveraging the isomorphism between (min,+) and (max,+) algebras, we prove formally that there is a much smaller bound than said lcm, allowing to greatly reduce the elementary operations required. We devise a new algorithm for (min,+) and (max,+) convolution, called super-isospeed , which avoids unnecessary computations and is considerably faster than the ones known so far, reducing the computation times by orders of magnitude. Unlike other works on the same problem, our method is both exact (i.e., does not introduce any approximation) and is not limited to operands of particular shapes (e.g., concave/convex, sub-/superadditive).
A high-order finite element method based on Hermite interpolation polynomials is described and a software implementation in the deal.II finite element library (https://dealii.org) is proposed. This method can be used to strongly enforce continuity of spatial derivatives as well as solution value over element boundaries, which is necessary for solving fourth order partial differential equations (PDEs) in primal form and useful for some second order PDEs. The implementation is verified by solving the wave equation based on a modification of the PDE with a Lax-Wendroff-like procedure, allowing a high-order solution method with negligible additional computing time.
Mixed precision “memory accessor” approaches decouple storage and compute precisions (data are stored and accessed in low precision, but computations are kept in higher precision) to reduce data accesses, improve accuracy, and simplify programming. In this work, we develop such a memory accessor aimed at accelerating sparse direct solvers and propose several new improvements. In particular, we propose a BLAS-based block approach that can directly rely on BLAS libraries for efficiency and portability. When considering BLAS-2 memory-bound operations like triangular solves, we observe that the performance adequately matches the storage cost, in multiple parallel settings, provided that the conversion from storage to compute precision is efficient and that the block size is suitably chosen. For the storage precision, we leverage custom floating-point types unsupported by hardware, and we take advantage of the recent AVX512-VBMI instruction set to reach an improved efficiency. We also consider rank-structured matrix representations such as the Block Low-Rank (BLR) format and explain how to optimize the memory accessor for such matrices. We present preliminary performance experiments using the sparse direct solver MUMPS with adaptive precision BLR approximations. Our results confirm the potential of these memory accessor approaches to achieve efficiency while optimizing storage.
Many real-world optimization problems exhibit dynamic characteristics, posing significant challenges for traditional optimization methods. Evolutionary Dynamic Optimization Algorithms (EDOAs) have been developed to address these challenges by adapting to changing environments over time. However, the reproducibility and consistency of experimental results in the literature remain limited due to the lack of publicly available source codes and the complexity of accurately re-implementing algorithms and performance evaluation protocols. To support the community, we introduce Evolutionary Dynamic Optimization LABoratory (EDOLAB), an open source MATLAB platform designed for both research and educational purposes. EDOLAB includes 27 EDOAs, four highly configurable benchmark generators, and a growing suite of performance indicators. The platform supports full parameter tuning, batch experiment management, parallel execution, and automated statistical comparisons-including rankings, significance testing, box plots, and performance trend visualizations over time. An educational application allows users to observe: (a) dynamic changes in a 2D problem landscape, (b) the movement of individuals in response to these changes, and (c) the ability of an algorithm to track moving optima. By providing an integrated environment for experimentation, benchmarking, and instructional use, EDOLAB promotes reproducibility, comparative analysis, and a deeper understanding of EDOAs in dynamic environments.
Modern GPUs feature specialized hardware for low-precision floating-point arithmetic to accelerate compute-intensive workloads that do not require high numerical accuracy, such as those from artificial intelligence. However, despite the significant gains in computational throughput, memory bandwidth utilization, and energy efficiency, integrating low-precision formats into scientific applications remains difficult. We introduce Kernel Float , a header-only C++ library that simplifies the development of portable mixed-precision GPU kernels. Kernel Float provides a generic vector type, a unified interface for common mathematical operations, and fast approximations for low-precision transcendental functions that lack native hardware support. To demonstrate the potential of mixed-precision computing unlocked by our library, we integrated Kernel Float into nine GPU kernels from various domains. Our evaluation on Nvidia A100 and AMD MI250X GPUs shows performance improvements of up to \(12\times\) over double precision, while reducing source code length by up to 50% compared to handwritten kernels and having negligible runtime overhead. Our results further show that mixed-precision performance depends not only on choosing appropriate data types, but also on tuning traditional optimization parameters (e.g., block size and vector width) and, when relevant, even domain-specific parameters.
SUNDIALS is a well-established numerical library that provides robust and efficient time integrators and nonlinear solvers. This paper overviews several significant improvements and new features added over the last three years to support scientific simulations run on high-performance computing systems. Notably, three new classes of one-step methods have been implemented: low storage Runge-Kutta, symplectic partitioned Runge-Kutta, and operator splitting. In addition, we describe new time step adaptivity support for multirate methods, adjoint sensitivity analysis capabilities for explicit Runge-Kutta methods, additional options for Anderson acceleration in nonlinear solvers, and improved error handling and logging.
We present PyRigi, a novel Python package designed to study the rigidity properties of graphs and frameworks. Among many other capabilities, PyRigi can determine whether a graph admits only finitely many ways, up to isometries, of being drawn in the plane once the edge lengths are fixed, whether it has a unique embedding, or whether it satisfied such properties even after the removal of any of its edges. By implementing algorithms from the scientific literature, PyRigi enables the exploration of rigidity properties of structures that would be out of reach for computations by hand. With reliable and robust algorithms, as well as clear, well-documented methods that are closely connected to the underlying mathematical definitions and results, PyRigi aims to be a practical and powerful general-purpose tool for the working mathematician interested in rigidity theory. PyRigi is open source and easy to use, and awaits researchers to benefit from its computational potential.
Trilinos is a community-developed, open-source software framework that facilitates building large-scale, complex, multiscale, multiphysics simulation code bases for scientific and engineering problems. Since the Trilinos framework has undergone substantial changes to support new applications and new hardware architectures, this document is an update to “An Overview of the Trilinos project” by Heroux et al. (ACM Transactions on Mathematical Software, 31(3):397-423, 2005). It describes the design of Trilinos, introduces its new organization in product areas, and highlights established and new features available in Trilinos. Particular focus is put on the modernized software stack based on the Kokkos ecosystem to deliver performance portability across heterogeneous hardware architectures. This paper also outlines the organization of the Trilinos community and the contribution model to help onboard interested users and contributors.
Many classical and modern finite element spaces are derived by dividing each computational cell into finer pieces. Such macroelements frequently enable the enforcement of mathematically desirable properties such as divergence-free conditions or C^1 continuity in a simpler or more efficient manner than elements without the subdivision. Although a few modern software projects provide one-off support for particular macroelements, a general approach facilitating broad-based support has, until now, been lacking. In this work, we describe a major addition to the FIAT project to support a wide range of different macroelements. These enhancements have been integrated into the Firedrake code stack. We provide numerical evaluation of the new macroelement facility.
We present the GPU implementation of the general-purpose interior-point solver Clarabel for convex optimization problems with conic constraints. We introduce a mixed parallel computing strategy that processes linear constraints first, then handles other conic constraints in parallel. The GPU solver currently supports linear equality and inequality constraints, second-order cones, exponential cones, power cones and positive semidefinite cones of the same dimensionality. We demonstrate that integrating a mixed parallel computing strategy with GPU-based direct linear system solvers enhances the performance of GPU-based conic solvers, surpassing their CPU-based counterparts across a wide range of conic optimization problems. We also show that employing mixed-precision linear system solvers can potentially achieve additional acceleration without compromising solution accuracy.
A numerical algorithm to decompose an exact low-rank skew-symmetric tensor into a sum of elementary (rank-1) skew-symmetric tensors is introduced. The algorithm uncovers this Grassmann decomposition based on linear relations that are encoded by the kernel of the differential of the natural action of the general linear group on the tensor, following the ideas of [Brooksbank, Kassabov, and Wilson, Detecting null patterns in tensor data, arXiv:2408.17425v2, 2025]. The Grassmann decomposition can be recovered, up to scale, from a diagonalization of a generic element in this kernel. Numerical experiments illustrate that the algorithm is computationally efficient and quite accurate for mathematically low-rank tensors.
C++ leans towards a memory-inefficient storage of structs: The compiler inserts padding bits, while it is not able to exploit knowledge about the range of integers, enums or bitsets. Furthermore, the language provides no support for arbitrary floating-point precisions. We propose a language extension based upon attributes through which developers can guide the compiler what memory arrangements would be beneficial: Can multiple booleans or integers with limited range be squeezed into one bit field, do floating-point numbers hold fewer significant bits than in the IEEE standard, and is a programmer willing to trade attribute ordering guarantees for a more compact object representation? The extension offers the opportunity to fall back to normal alignment and native C++ floating point representations via plain C++ assignments, no dependencies upon external libraries are introduced, and the resulting code remains (syntactically) standard C++. As MPI remains the de-facto standard for distributed memory calculations in C++, we furthermore propose additional attributes which streamline the MPI datatype modelling in combination with our memory optimisation extensions. Our work implements the language annotations within LLVM and demonstrates their potential impact through smoothed particle hydrodynamics benchmarks. They uncover the potential gains in terms of performance and development productivity.
We introduce the shape module of the Python package Geomstats to analyze shapes of objects represented as landmarks, curves and surfaces across fields of natural sciences and engineering. The shape module first implements widely used shape spaces, such as the Kendall shape space, as well as elastic spaces of discrete curves and surfaces. The shape module further implements the abstract mathematical structures of group actions, fiber bundles, quotient spaces and associated Riemannian metrics which allow users to build their own shape spaces. The Riemannian geometry tools enable users to compare, average, interpolate between shapes inside a given shape space. These essential operations can then be leveraged to perform statistics and machine learning on shape data. We present the object-oriented implementation of the shape module along with illustrative examples and show how it can be used to perform statistics and machine learning on shape spaces.
Efficiently solving nonlinear equations underpins numerous scientific and engineering disciplines, yet scaling these solutions for challenging system models remains a challenge. This article presents NonlinearSolve. jl-a suite of high-performance open source nonlinear equation solvers implemented natively in the Julia programming language. NonlinearSolve.jl distinguishes itself by offering a unified API that accommo-dates a diverse range of solver specifications alongside features such as automatic algorithm selection based on runtime analysis, support for static array kernels for improved GPU computation on smaller problems, and the utilization of sparse automatic differentiation and Jacobian-free Krylov methods for large-scale problem-solving. Through rigorous comparison with established tools such as PETSc SNES, Sundials KINSOL, and MINPACK, NonlinearSolve.jl demonstrates robustness and efficiency, achieving significant ad-vancements in solving nonlinear equations while being implemented in a high-level programming language. The capabilities of NonlinearSolve.jl unlock new potentials in modeling and simulation across various domains, making it a valuable addition to the computational toolkit of researchers and practitioners alike.
The finite precision implementation of mathematical functions frequently depends on polynomial approx-imations. A key characteristic of this approach is that rounding errors occur both when representing the coefficients of the polynomial on a finite number of bits, and when evaluating it in finite precision arithmetic. Hence, to find a best polynomial, for a given fixed degree, norm, and interval, it is necessary to account for both the approximation error and the floating-point evaluation error. While efficient algorithms were already developed for taking into account the approximation error, the evaluation part is usually a posteriori handled, in an ad hoc manner. Here, we formulate a semi-infinite linear optimization problem whose solution is a best polynomial with respect to the supremum norm of the sum of both errors. This problem is then solved with an iterative exchange algorithm, which can be seen as an extension of the well-known Remez exchange algorithm. An open source C implementation using the Sollya library is presented and tested on several examples, which are then analyzed and compared against state-of-the-art Sollya routines.
