Papers

This manuscript presents a performance evaluation of Fast Fourier Transform (FFT) implementations on emerging processors supporting the RISC-V Vector Extension (RVV 1.0). By introducing juFFTe, a light-weight high-performance library for discrete Fourier transforms, it is demonstrated how effective vectorization of performance-critical FFT kernels can be achieved on RVV-enabled hardware. Comprehensive benchmarks on three RVV 1.0-ready processors, the SiFive X280, the X100 core of the SpacemiT K3 and the C920v2 core of the Sophon SG2044, reveal substantial performance improvements of juFFTe (https://github.com/FZJ-JSC/juFFTe) over the widely used FFTW3 library. Although RVV-enabled platforms show promising results at this stage of development, a comparison with AMD's Zen 5 architecture indicates that RISC-V needs further maturing to reach the performance of established micro-architectures.

We introduce the wavemap: a spatial portrait of noise effects that assigns each site a per-noise-level arrival delay l_\gamma (v) and cross-entropy loss L_\gamma(v). These observables are exact at the lightcone frontier, where bond dimension \chi is small and the simulation is most faithful. Eigenvalue analysis of the composed gate-plus-noise Pauli transfer matrices confirms that the studied noise is pure amplitude damping: the spatial propagation pattern is entirely determined by the gate, making the wavemap a model-free noise diagnostic. We apply the multi-product formula (MPF) to recover the noiseless Pauli weight field from the noisy samples, subject to the Lieb-Robinson causal constraint nMPF <= nnl . Fitting time-adaptive coefficients \alpha(t) over the frontier recovers up to 55% of the information loss relative to the best noisy sample, exploiting the fact that the frontier is where truncation error is smallest. On an IBM heavy-hex lattice with heterogeneous hardware noise the method identifies an information-starved regime, pointing to calibrated synthetic noise as the next required experiment.

Quantum systems grow in complexity so rapidly that even modest models become difficult to simulate, creating a strong need for methods that can handle high-dimensional data, also known as tensors. In this work, we investigate a novel Jacobi-type tensor algorithm for tensor decomposition and develop a CUDA-based algorithm that supports tensors of arbitrary order on a single GPU. We test the implementation on NVIDIA H100 GPUs and show that the algorithm converges correctly for diagonalizable tensors up to nine dimensions, with runtime scaling in a predictable way as tensor order grows. Finally, our general algorithm outperforms the original MATLAB reference by more than two orders of magnitude.

Data assimilation (DA) and inverse modelling are indispensable tools for combining physics-based models with observations, yet the software that implements them is often fragmented: sequential (Kalman and particle) filters and variational (3D/4D-Var) estimators are typically developed as separate codebases, each tied to a specific class of forward model - for example, simple time-marching schemes for ordinary differential equations (ODEs), or full finite element (FE) discretisations of partial differential equations (PDEs). In this work, we present Opal.jl, a Julia package that overcomes this fragmentation by providing a unified environment for a wide range of DA methods and forward model backends, all accessible behind a single high-level interface. Starting from a handful of basic methods, complex inference tools can be built by composing them with advanced capabilities, such as covariance localisation/inflation, online noise covariance estimation, bias-aware correction, and kriging-based calibration of reduced-order surrogates. The package integrates natively with the SciML ecosystem for ODE-governed systems, and with Gridap/GridapROMs for both full-order and reduced-order discretisations of PDEs. Notably, the API is designed to allow a reduced-order surrogate to replace a full-order solver with no change to the code, even at the driver level. We demonstrate the library on a Lorenz-63 benchmark, on a Van der Pol oscillator problem with biased observations, on a turbulent Navier-Stokes flow past a square cavity, where the unknown Reynolds number is inferred from sparse velocity and pressure measurements, and finally on a heat equation using a reduced-order surrogate instead of a more standard full-order model.

Data assimilation (DA) is an essential tool for prediction and understanding in the geosciences. DA combines simulation programs representing scientific knowledge with observations that constrain system dynamics, resulting in analyses and forecasts that incorporate both knowledge and data. DA tasks can be addressed with a diverse toolset, including variational, ensemble and learning-based methods. In particular, many recent works have proposed using automatic differentiation tools for variational, learning-based or hybrid methods. However, comprehensive comparisons across algorithms and dynamical systems remain challenging, due to the incompatibility of simulation and assimilation codes, inflexible handling of spatial and temporal discretizations, specialization of DA methods to specific simulations, and limited support for automatic differentiation and parallel computation in simulations. To address this challenge, we introduce Automatic Differentiation for Data Assimilation (ADDA), a software framework for defining and working with system states, simulations, observation schemes and DA methods. ADDA provides a powerful and flexible set of base classes for representing dynamical systems and observation operators, with support for collocated and staggered grids, unstructured meshes, Lagrangian state variables and irregular or continuous-time observations. Parallel processing and differentiability are first-class features, with support for batch axes and automatic differentiation throughout. ADDA is implemented in PyTorch library, but supports DA for JAX-based computation of dynamics and their gradients. To demonstrate its features, we further provide differentiable, ADDA-compatible implementations of 10 dynamical systems of various dimensionalities and scales, from which we design multiple illustrative DA examples. All of our code is publicly available at https://github.com/m-dml/ADDA.

Superspace extends spacetime by anticommuting Grassmann coordinates, whose indices may transform under spin and flavour groups. Decomposing Grassmann monomials into independent invariant structures is a central step in constructing superfield expansions, but becomes increasingly difficult in extended superspace. We present a general representation-theoretic procedure for decomposing the space $\Lambda^n(\mathbb{C}^{d_S}\otimes \mathbb{C}^{d_F})$, where $d_S$ and $d_F$ are the dimensions of the spin and flavour representations. The branching multiplicities for $\mathrm{GL}(m)\downarrow{\mathrm{Sp}(m)}$ and $\mathrm{GL}(m)\downarrow{\mathrm{SO}(m)}$ are computed by exact Weyl-character comparison at deterministic rational sample points, avoiding Littlewood-Richardson modification rules. We establish completeness of the character ansatz and analyse the computational complexity, showing that it is polynomial in the number of partitions of $n$ and in the group rank, without dependence on Littlewood-Richardson combinatorics. We validate the resulting decompositions via dimension checks. We provide an accompanying standalone Julia script which implements this algorithm for groups of arbitrary rank, and produces the explicit contractions with invariant tensors.

We present "Alakazam", a Julia package for performing symbolic manipulations of tensor expressions arising in the study of field theories and quantum mechanics. This will serve as a powerful tool for studying spacetime and supersymmetric structures, including support for calculations in superspace. The package implements many key features such as raising and lowering indices, simplification of terms with dummy indices, replacement patterns, differential and variational calculus, bosonic and fermionic grading, gamma matrix identities, and Lie brackets and commutators. We present notable speed improvements over established systems. Some possible applications include the study of descendants and invariants.

Computing eigenvalues and eigenvectors is at the heart of the solution of many non-linear problems. For instance, finding the roots of polynomial systems reduces to computing joint eigenvectors of operators of multiplication. Similarly, tensor decomposition can be performed via the joint diagonalization of submatrices of the Catalecticant of the tensor. We describe and illustrate symbolic-numeric methods for computing the solutions of these algebraic problems from the computation of joint eigenvectors of commuting operators, and for analysing their multiplicity structure, as well as their implementation in the package AlgebraicSolvers.jl.

The titular object of this paper is an analogue of the Mandelbrot set over the 3-dimensional real algebra whose multiplication is the bilinear extension of the rock-paper-scissors operation. In order to study the dynamics of the mappings $x\mapsto x^2+c$ in this setting, a notion of holomorphy for functions on general finite-dimensional real algebras is introduced. Under a mild assumption it is shown that for such algebras holomorphy is always equivalent to solving a finite system of linear first-order PDEs generalizing the Cauchy-Riemann equations. Code is provided for generating animations of these fractals and for exploring them in a video game format.

We introduce a bilingual domain-specific language (DSL) for modelling compute kernels within a generic solver for hyperbolic partial differential equations (PDEs). Users express PDE terms, i.e.~the underlying physics, in a familiar native language such as C or C++, while the numerical scheme is specified in a Python-embedded DSL, DSLHyPE. DSLHyPE's compiler lowers the Python description to MLIR and introduces a translation pass that integrates it with native code likewise mapped to MLIR. Our approach keeps the numerical representation and the physics implementation separate for as long as possible, while delegating optimization to the compiler through existing MLIR optimization passes. This separation of concerns benefits researchers developing numerical schemes on top of existing PDE implementations or with applications involving nonlinear systems whose PDE terms must solve PDEs themselves. We demonstrate the feasibility of the approach using a gravitational-wave solver and a matter-evolution solver on x86 processors and H200 GPUs.

Anytime solvers return usable solutions before they terminate and improve them while time remains. Their progress is commonly summarized by a primal integral, a final gap, or convergence curves. Each summary leaves consequential choices open: how to score a run before its first feasible solution, whether poor incumbents are truncated, and whether the reference value is updated after the experiment or version-frozen beforehand. These choices become visible when runs produce no valid incumbent or improve a published best-known value. We study a normalized signed primal integral built from a bounded relative gap. It assigns an intrinsic worst value to an empty run, requires no acceptance threshold, and assigns negative instantaneous gaps to incumbents that beat a frozen reference. We compare the smooth difference-over-sum kernel with a signed version of Berthold's max-normalized gap and use the former as a working default. We also distinguish analysis-time from version-frozen reference policies and recommend reading the score with the raw final gap, mean convergence curve, and target-attainment curve. We evaluate these choices on a five-arm routing campaign and two model-fidelity ladders. The two signed kernels preserve every panel ordering in this study, whereas a common acceptance threshold compresses the distances between arms and reverses one panel ordering. Updating 54 of the 212 reference values changes score levels and removes all negative scores, while leaving the observed panel orderings unchanged. An exploratory screen also identifies 12 panel comparisons in which similar integral scores conceal materially different endpoints or attainment rates. The implementation, frozen inputs, and generators are openly released.

We present the Linear Algebra Aware Benchmarks (LAAB) framework for systematically assessing and reporting the performance of mathematical library installations on HPC systems. Mathematical libraries provide interfaces for operations that form the computational building blocks of scientific applications. Reporting their performance is important for assessing application efficiency, estimating compute-time requirements, and preparing resource-allocation requests. In this paper, we define four objectives for performance reporting: 1) traceability, linking each report to the exact library installation and execution settings; 2) compatibility, relating library-operation performance to higher-level scientific applications that use them; 3) reliability, supporting interpretation in the presence of measurement variability; and 4) accessibility, ensuring that reports, benchmark definitions, and relevant metadata are available for inspection and reproduction. We then present the design of LAAB and show how it addresses the challenges associated with these objectives.

A lattice Boltzmann scheme for the three-dimensional compressible Navier–Stokes–Fourier equations, derived automatically from the declared system by a symbolic compiler, is validated against exact solutions and published reference data. The declared system carries the viscous stress and the heat flux as transported state, and is discretized on a D3Q7 lattice in single precision. Against the exact Sod and Becker solutions the captured shock thickness converges at first order. On the supersonic Taylor-Green vortex at $M_0 = 1.25$ the scheme at $512^3$ matches the reference dilatational dissipation more closely than seven compared solvers, by thirty percent over the next best. Every operator in this solver, the generated collision and constitutive closure and the added shock sensor alike, reads only the cell it acts on, and data reaches a neighbor only by streaming along the lattice characteristics.

Thread scaling in Hexaly on the Time-Dependent Vehicle Routing Problem with Time Windows (TDVRPTW) depends on modeling, and on the language binding through which the model reaches the solver. We compare two encodings, both reaching the solver through its C++ binding: one evaluates continuous travel-time functions exactly through external callbacks, which that binding evaluates concurrently, while the other approximates them with time slices evaluated natively by the solver. Independent CPU accounting confirms that both use the cores they are allocated, but they respond differently to width. The external-function encoding stabilizes, reducing its seed dispersion by about 30% up to 8 threads before widening again, yet its pooled quality barely moves, because nine of the ten instances improve while the hardest one degrades by enough to cancel them. For the time-sliced encoding we first select a discretization on feasibility under the original travel-time functions rather than on approximation error, and the retained setting then improves monotonically with the thread count, gaining about a quarter of its final gap and reducing the seed dispersion by nearly half between 1 and 16 threads. At equal single-threaded budget the two encodings are close, and the small pooled edge of exact evaluation comes from that same hardest instance, so the practical advantage of the discretization is thread scaling rather than fidelity traded for speed. These descriptive results support multi-threading as a way to improve solution quality and stability when the model can exploit it, and they show that binding-level constraints on external evaluation must be measured, not assumed.

We present PDHCG-CQP, a GPU-accelerated first-order solver for large-scale conic convex quadratic programming. PDHCG-CQP supports affine constraints and Cartesian products of nonnegative, second-order, rotated second-order, exponential, and three-dimensional power cones. At its core is a restarted averaged primal-dual hybrid gradient (PDHG) method, whose primal update is computed inexactly by solving a conic quadratic proximal subproblem with projected gradient iterations. We establish local linear convergence of the restarted averaged scheme with both exact and inexact primal proximal evaluations under a uniform local quadratic-growth condition on the smoothed primal-dual gap. We further show that this condition holds under strict complementarity by exploiting a rotated second-order-cone lifting together with local primal and dual regularity conditions. Our C/CUDA implementation combines matrix-free linear algebra, batched cone projections, adaptive inner solves, reflected-Halpern acceleration, and fully device-resident KKT residual computations. It also supports multi-GPU execution through a two-dimensional partitioning of the problem data. Extensive experiments on standard and large-scale quadratic programming (QP), convex quadratically constrained quadratic programming (QCQP), second-order cone programming (SOCP), and quasilinear Fisher equilibrium benchmarks demonstrate that PDHCG-CQP achieves state-of-the-art robustness among first-order solvers while scaling efficiently to 8 GPUs and instances with up to $4.4\times10^8$ stored primal coordinates. PDHCG-CQP is open source and available at https://github.com/Lhongpei/PDHCG.

Developing nonlinear models that are both expressive and computationally efficient remains a challenge in machine learning and nonlinear system identification. Tensor network kernel machines (TNKM) address this challenge by combining nonlinear feature representations with compact low-rank tensor-network parameterizations. However, practical and extensible software frameworks for developing TNKM models remain limited. In this work, we introduce "tnkm", an open-source Python library for constructing and training TNKM models using JAX. The library provides a unified interface for combining different feature maps, tensor-network architectures, and optimization strategies, including alternating least squares and gradient-based methods. We demonstrate the capabilities of "tnkm" on nonlinear benchmark problems, showing that the implemented models achieve competitive prediction accuracy while retaining compact parameterizations and efficient training. The proposed framework facilitates reproducible development and application of tensor-network-based learning methods.

Parameter-dependent linear matrix inequalities (PD-LMIs) require holding over a continuous domain. When the scheduling parameters vary with time, derivatives of parameter-dependent decisions may also enter the conditions. Since semidefinite programming solvers require finitely many constraints, we introduce GriD-LMIA, the Gridding-based Differentiable PD-LMI Assembler. It converts the conditions that need to hold on a continuous domain into finitely many sufficient LMIs in MATLAB. It first partitions the domain with a hyper-rectangular grid and represents known data and decisions on each cell with tensor Bernstein polynomials. The package exports direct Bernstein, Pólya-elevated, and sum-of-squares-based certificates through YALMIP. Examples are given to examine the balance among grid density, decision degree, and certificate choice.