Peer Revenue Publish or profit

Papers

Fresh from arXiv, typeset properly. Highlighted words trade on the market.

Computer Science Mathematical Software

cs.MS Watch Mathematical Software

Wednesday 14 October · still arriving

Mathematical Software · Cross-listed from cs.PF

How the Dirty Set Moves: Characterizing Schedule Reuse in Incremental Sparse Matrix Multiplication

Hunter Horne, Eric Spencer, Alex Henningfield et al. (+1 more)

Sparse matrices, in which most entries are zero, are central to scientific simulation, graph analytics and machine learning, and their multiplication with dense inputs is a core parallel kernel. In many applications the matrix stays fixed while only a few input rows change per step, so the result can be updated from those rows alone, which we call the dirty set. Sparse kernels often add a second optimization based on the inspector-executor principle, in which a tiled schedule is prepared once and reused across calls. However, when this reuse remains profitable as the dirty set moves is not well characterized. In this paper, we show that retained tiled execution rarely beats direct CSC in incremental sparse matrix multiplication, although incremental execution itself remains highly profitable. Across 21 SuiteSparse matrices and controlled patterns of motion, a direct column-wise incremental update is 2.1-8.2x faster than the best of eight full-recomputation codes, among them Intel oneMKL and PyTorch. Yet the same direct update outperforms a tiled executor that retains tile preparation and caches assembled schedules in 779 of 837 workloads. Separating reuse into per-tile preparation and cached schedules shows that the benefit comes from preparation; extra cache capacity seldom helps because a tile set rarely returns after a different one has intervened. Tiling is faster only when tiles merge many entries into each output write, a structural property that we predicted before measurement and confirmed on four held-out matrices of the same family. These results show that the structure of the matrix and the motion of the dirty set, rather than the sophistication of the schedule, determine whether reuse is profitable. We anticipate that our methodology and guidelines will help developers decide when to build schedule reuse into incremental sparse runtimes.
Read the full abstract

All authors Hunter Horne, Eric Spencer, Alex Henningfield, Arslan Bisharat

Tuesday 13 October · still arriving

Mathematical Software

TOMOS: Eliminating Data Duplication in Mixed-Precision Matrix Computations

Vasileios Georgiou and Petros Drineas

Mixed-precision algorithms exploit reduced-precision hardware, however computing with the same matrix at distinct precisions typically requires duplicate storage, repeated conversions, or operation-emulation. We introduce TOMOS, a composite-precision matrix format that provides implicit 32-bit and 16-bit access from a single representation. Each FP32 entry is split into a 16-bit primary word for low-precision computation and a complementary word that enables exact FP32 reconstruction. The encoding scheme improves numerical quality of the low-precision view through subnormal representations, and allows exact recovery of high-precision inputs. TOMOS uses FP32-sized storage, with data movement scaling with compute precision. Optimized GEMM and GEMV kernels achieve throughput comparable to vendor libraries on NVIDIA A100 GPUs, while a mixed-precision least-squares solver achieves up to 1.6x speedup over FP32 without redundant matrix allocations.
Read the full abstract

Mathematical Software · Cross-listed from quant-ph

PADO-Pauli: A GPU-Accelerated Implementation of Pauli Propagation

Hyunwoo Kim and Youngseok Lee

Variational quantum algorithms utilize parameterized circuits many times during training, at new parameter values, over batches of data, and with gradients. Pauli propagation approximates the expectation values of such circuits classically by evolving the observable in the Heisenberg picture as a truncated sum of Pauli strings. We present PADO-Pauli, a GPU engine for Pauli propagation that compiles a circuit and its observables once into sparse matrices in which each rotation angle contributes only a cosine and a sine factor, so every later evaluation reuses the compiled program. It evaluates batches of parameters and data inputs embedded as rotation angles. By default its gradients come from a vector-Jacobian product registered with PyTorch's automatic differentiation, so they flow through the angles to any classical layer that produces them. PADO-Pauli also provides truncation by Pauli weight and coefficient size, and noise channels. Zero-filtering removes terms that do not contribute to an expectation value in the all-zero reference state and makes the batched forward pass about 1.5 times faster. Compared to cuPauliProp on an NVIDIA A100, PADO-Pauli evaluates up to 10 times faster on 4 to 12 qubits and uses up to 3.1 times less memory on 16 qubits. We check its accuracy against exact statevector simulation up to 28 qubits. These features make PADO-Pauli a practical toolkit for training variational quantum algorithms, including quantum machine learning models.
Read the full abstract

Monday 12 October

Mathematical Software · Cross-listed from cs.LG

Learning to Design Convex Relaxations for Nonconvex Optimization

Jinrun Liu, Bai Yan, Tengfei Liu et al. (+5 more)

Convex relaxations make nonconvex optimization tractable by replacing hard nonconvex problems with convex surrogates that give valid bounds. Yet their design is still largely manual: a human expert must decide where to relax, which relaxation to use, and how much formulation growth is worth paying. We show that this design process can be learned. The method proposed in this paper, AutoConvexRelax, represents the symbolic formulation of a nonconvex problem as a tripartite graph linking variables, symbolic terms, and their objective or constraint positions. It then learns a policy that composes symbolic rewrites into problem-specific relaxation trajectories as the formulation evolves. We prove that finite action sequences generated by the learned policy preserve relaxation validity under explicit domain, occurrence-context, and representation-consistency conditions. We further establish computable error bounds for mixed relaxations and sufficient conditions for feasible recovery and global-optimality certification. As a first instantiation, we apply the method to nonconvex quadratic problems and fractional problems. On 241 benchmark quadratic instances, the learned policy strengthens McCormick lower bounds by 13.52%. Its mean lower bound is 3.70% below strengthened semidefinite relaxations on the normalized comparison metric, while it uses substantially smaller formulations. The learned policy retains this trade-off on fractional instances and reaches near-semidefinite bounds on beamforming and sensor-network localization instances. These results suggest that convex relaxation design can move from a fixed handcrafted choice to an automated learning process that adapts to problem structures. Code is available at https://github.com/auto4opt/AutoConvexRelax.
Read the full abstract

All authors Jinrun Liu, Bai Yan, Tengfei Liu, Qunfeng Liu, Jian Yang, Shi Cheng, Yuhui Shi, Qi Zhao

Mathematical Software · Cross-listed from cs.LG

GpFlash: Maximizing Speed and Minimizing Memory for Gaussian Process Training and Prediction

Andre Marchildon

This paper presents GpFlash, a Gaussian process library that is focused on fast training and posterior evaluations, minimal memory usage, and scalability to high-dimensional problems. GpFlash was compared to the popular GP libraries scikit-learn, GPy, GPflow, GPyTorch, and SMT on a MacBook Pro and the supercomputer Trillium. GpFlash was found to use less memory and to be faster for training and posterior evaluations relative to the other GP libraries for all cases considered. For example, when using $\nx = 10^4$ training points in a 10-dimensional space GpFlash required 9.6 to 35 times less memory and was 2.6 to 54 times faster than the other GP libraries for the evaluation of the marginal log likelihood and its gradient. The memory reduction was achieved by only using one matrix of size $\nx \times \nx$ to store different symmetric matrices in its lower and upper triangular entries. The computational speedup was the result of designing algorithms to access memory efficiently, using BLAS and LAPACK methods extensively, and factoring kernel equations so that expensive exponential calculations could be reused. GpFlash was more than 440 times faster than all of the other GP libraries considered when evaluating the posterior mean together with its gradient for the case with $10^4$ training points and 10 dimensions. GpFlash is available at https://github.com/marchildon/gpflash.
Read the full abstract

Sunday 11 October

Mathematical Software · Cross-listed from math.ST

Separating Statistical Inference from Experimental Runtime: An Architecture for Sequential Scientific Analysis

Alexander Stepanov and Andrej Novikov

Sequential statistical procedures are increasingly embedded in scientific systems in which observations arrive incrementally from heterogeneous sources. In such systems the statistical update is only one part of the computation: observations must also be acquired, aligned across inputs, ordered, delivered, persisted, and replayed, and results must be passed to systems that display them or act on them. Sequential methods are nevertheless missing from the general software toolkit. They run inside proprietary experimentation platforms and commercial statistical packages, or exist as offline functions over complete data, and the batch fit interface that dominates machine-learning libraries has no place for a procedure that listens, keeps state, and decides when to stop. We propose an architecture that separates statistical inference from data topology, execution runtime, run lifecycle, persistence, and external action. Its central abstraction is the Sequential Procedure: a state-transition function that maps an inferential state, a valid statistical input, and a statistical random state to a new state and a scientific output. The procedure does not define how observations are acquired, paired, scheduled, stored, or acted upon. We specify the architecture semantically rather than through an object model, derive seven conformance invariants that any implementation can be tested against (among them synchronous/asynchronous, checkpoint, replay, and batch equivalence), and describe a dependency-free reference implementation in which all seven hold bitwise. Lorden's 2-SPRT for the Kiefer-Weiss problem serves as the reference procedure; the same contract accommodates sequential estimators, Bayesian filters, change-point detectors, and stateful predictive models.
Read the full abstract

Saturday 10 October

Mathematical Software

manylp: reproducible, exact, and batched linear programming for genome-scale dynamic flux balance analysis

Fabio Cumbo

Dynamic flux balance analysis (dFBA) and its community, spatial and ensemble variants predict microbial metabolism by solving one linear program (LP) per species and time step, often 10^6 to 10^8 LPs per study. These LPs share their constraint matrix and objective and differ only in their bounds, yet general-purpose solvers treat each one independently, and GPU first-order solvers struggle with the degenerate, badly scaled matrices of genome-scale models. Moreover, metabolic optima are usually faces rather than vertices, so solvers return arbitrary flux vectors and dFBA predictions depend on the solver used. We present manylp, a batched solver for such bounds-parametric LP families that returns certified vertex solutions. Building on the bunching technique of stochastic programming, manylp caches optimal bases and certifies whole batches with one dense matrix product per basis on CPUs or GPUs; LPs outside the cached critical regions are re-solved by a warm-started lexicographic dual simplex, and new bases are propagated across the batch. Each LP is certified optimal and unique under a lexicographic flux rule, or infeasible, in independent double-precision arithmetic. On a genome-scale dFBA workload, manylp certified 47,889 LPs per second on a CPU, 12 to 33 times faster than HiGHS, Gurobi, FICO Xpress, GLPK and OR-Tools (32 processes), and reached 682,000 LPs per second on one GPU for community ensembles, needing 76 simplex solves for 94.6 million LPs. Integrated into the muODE simulator, it accelerated single trajectories 84-fold and 16-member ensembles 635-fold, and made adaptive Runge-Kutta integration affordable, reducing biomass errors 83-fold. Integrating 49 equally optimal flux policies in one batched run showed that alternative optima change predicted acetate 4.2-fold and shift community composition. manylp is open source: https://github.com/cumbof/manylp.
Read the full abstract

Thursday 8 October

Mathematical Software · Cross-listed from cs.LG

Coordinate-Free Automatic Differentiation

Alex Guerra, Joshua Aduol, Yaqian Tang et al. (+3 more)

Automatic differentiation (AD) is a cornerstone of modern machine learning. Since all AD frameworks implicitly fix a basis of Euclidean space to define gradients and Jacobians, they natively process mathematical objects that are arrays of a fixed set of coordinates. This makes it awkward to express objects such as rotations, volume forms, and curvatures that arise, for example, in applications of AI to science. Bespoke workarounds lack reusability. We show that (i) a single, basis-independent derivative operator can represent many geometric differential operations, and (ii) this operator can be expressed entirely in terms of existing AD primitives. This reusable transformation can be composed with JAX batching, compilation, and hardware acceleration. Geometry enters as algebraic data. We experimentally demonstrate its reuse across different manifold operators, its composition through a spacetime inverse-problem pipeline, and its use in optimization over different rotor subgroups.
Read the full abstract

All authors Alex Guerra, Joshua Aduol, Yaqian Tang, Daniel Williams, Ryan P. Adams, Peter Orbanz

Mathematical Software

asdex: Automatic Sparse Differentiation in JAX

Adrian Hill and Guillaume Dalle

Many tasks in scientific computing and machine learning require the Jacobian or Hessian matrix of a function. Automatic differentiation (AD) computes these derivatives to machine precision, but materializing a dense $m \times n$ Jacobian requires $n$ forward-mode or $m$ reverse-mode AD passes, one per column or row. For a large class of functions, each output depends on only a few inputs, making the derivative matrix sparse. Automatic sparse differentiation (ASD) exploits this structure in four steps: detection of the input-agnostic sparsity pattern, coloring of a graph to group columns or rows that can share an AD pass, compressed differentiation to compute a compressed derivative matrix with one AD pass per color, and finally decompression into the original sparsity pattern. The number of colors, and hence of AD passes, is often independent of the problem dimension: a banded Jacobian with $b$ contiguous bands, for instance, only ever requires $b$ colors, regardless of its size. asdex offers the first standalone ASD toolkit in the popular JAX ecosystem. With $\texttt{asdex.jacobian}$ and $\texttt{asdex.hessian}$, it provides sparse drop-in replacements for $\texttt{jax.jacobian}$ and $\texttt{jax.hessian}$.
Read the full abstract

Tuesday 6 October

Mathematical Software · Cross-listed from physics.comp-ph

acousticHelmholtzFoam: an MPI-parallel unstructured finite-volume solver for frequency-domain acoustics in heterogeneous two-phase media with perfectly matched layers

Chuanchao Xu, Jun Liu, Dieter Bothe et al. (+1 more)

acousticHelmholtzFoam is an OpenFOAM application that computes time-harmonic acoustic pressure fields in heterogeneous gas-liquid media on unstructured polyhedral meshes. The complex pressure is stored as real and imaginary volume fields, the material coefficients are constructed from a volume-of-fluid phase fraction with geometric face-fraction averaging, and open boundaries are truncated by Cartesian perfectly matched layers (PMLs) whose damping tensor is generated by the companion utility setPMLFields. The two coupled real-valued equations are discretized with the cell-centred finite-volume operators of OpenFOAM and assembled explicitly, including processor couplings and boundary contributions, into a distributed PETSc matrix that is solved by default with the MUMPS sparse direct solver. The solver writes pressure amplitude, complex velocity, acoustic radiation pressure and time-averaged momentum flux for acoustic-force evaluation. This paper documents the implemented equations, the case interface, the assembly and solution algorithm, the build procedure and the test infrastructure. A layered gas-liquid-PML case with an analytical solution serves as an end-to-end regression: on 8000 cells the whole-domain relative complex-pressure error is $1.197\times10^{-3}$, and serial and MPI solutions on up to eight ranks agree to within $10^{-12}$ in relative $L_2$ norm. A single-axis acoustic-levitator height sweep with constant-volume spheroidal drops demonstrates the two-phase material input and the force output, and strong- and problem-size-scaling measurements characterize the direct-solver workflow. The numerical formulation and its verification are presented in a companion paper (arXiv:2610.01593).
Read the full abstract

All authors Chuanchao Xu, Jun Liu, Dieter Bothe, Tomislav Maric

Mathematical Software · Cross-listed from math.NA

JAXIGA - A software framework for GPU-enabled and differentiable isogeometric analysis

Cosmin Anitescu and Yongjie Jessica Zhang

In this paper, we present JAXIGA, an open-source library for isogeometric analysis based on JAX. The library separates the construction of the mesh from the computations which are traced for compilation and differentiation. This allows geometry evaluation, assembly and the supported solvers to be combined in a just-in-time compiled program, for which first-order reverse-mode derivatives can be computed. We consider three numerical methods obtained from a single pointwise energy density. The Galerkin and energy-minimization methods use the same discrete potential and have the same stationary points, while collocation enforces the corresponding strong form. The sensitivities of a converged solution are obtained by solving one adjoint problem, and their accuracy depends on the tolerances of both the primal and adjoint problems. We also describe hierarchical extraction for local refinement of the variational methods and the use of smooth spline spaces for fourth-order problems. The numerical results are compared with analytical solutions and independent FEniCSx and JAX-FEM computations. For the reported GPU benchmark with 36,448 unknowns, assembly is twenty-four times faster than on the host CPU. When the conjugate-gradient solver is also run on the GPU, the total solution time is reduced from 31.8 to 1.38 seconds, and the result agrees with the direct solution to $5\times10^{-11}$. These results show that the choice of solver and the reuse of compilation are important for obtaining an overall reduction in computational cost.
Read the full abstract

Mathematical Software · Cross-listed from cs.LG

CNet: A Complex-Valued Deep Learning Framework with Wirtinger Autodifferentiation and FFT–Hadamard Convolution

Marcel Crasmaru

CNet is a C++/CUDA framework for building deep complex-valued neural networks (CVNNs) and optimizing complex functions by gradient descent with Wirtinger derivatives. Complex models are underexplored yet natural where data is intrinsically complex – RF/IQ communications, MRI k-space, radar/SAR, audio spectra – and phase carries information real networks discard. CNet is physics-native: a network is a cascade of complex (often unitary) operations on an amplitude vector, and classification is a Born-rule measurement p_k = |z_k|^2/||z||^2 rather than a softmax over real logits. Its library of complex layers makes conv(x,k) = IFFT(FFT(x).FFT(k)) learnable via signal-processing primitives – spectral-padding local kernels (Pad), the inverse DFT, a magnitude nonlinearity |z|^2 (CModulus2), holomorphic powers z^M (CPower), and mean pooling (MeanPool) – with GPU kernels for every layer, a GPU true-Adam optimizer, and a reduced-memory inference mode. Three studies. (1) A fully complex FNet-style causal language model, built on a new O(N log N) causal Fourier mixer, matches a param-matched real-valued causal FNet on tiny-shakespeare in under half the steps. We introduce Born-rule attention: a content-based causal mixer whose weights are quantum-measurement probabilities |<Q_k,K_j>|^2 of one token against another – the first attention mechanism built on the Born rule. With rotary position embeddings and per-token complex normalization, a fully complex Born-rule model reaches 1.52 nats/char on tiny-shakespeare, matching softmax attention and beating the Fourier mixer by 0.23. (2,3) Two bottleneck analyses – radio-modulation classification (RML2016.10a) and the Fourier phase problem of coherent-diffraction imaging – isolate where CVNNs need new operators. The complex machinery learns the correct structure; open problems are concrete and operator-level. Code: https://github.com/crasmarum/CNet
Read the full abstract

Mathematical Software · Cross-listed from math.QA

quantum-group: Exact symbolic verification of R-matrix and Yang–Baxter identities for $U_q(\mathfrak{sl}_2)$ and $GL_q(2|1)$

Taha Berk Terekli

Explicit R-matrices of quantum groups are routinely transcribed from the literature, yet whether they satisfy an identity depends on basis order, tensor placement, the coproduct and, for superalgebras, the Koszul sign rule. quantum-group is a Python package that builds such matrices exactly with SymPy and returns every identity as an exact residual matrix, whose nonzero entries locate the convention that disagrees. It covers U_q(sl_2) modules, tensor products and the fundamental R-matrix, and the graded Yang–Baxter equation for GL_q(2|1). Applied to the thesis it was written for, it exposed an R-matrix that satisfies the Yang–Baxter equation but intertwines the opposite coproduct.
Read the full abstract

Mathematical Software

Scen-Opt: A Scenario Optimization Toolbox for Data-Driven Convex Programming

Ben Wooding, Simone Garatti, Marco C. Campi et al. (+1 more)

The scenario approach is a well-established statistical framework for data-driven decision-making. In particular, in data-driven optimization, the scenario approach unveils how the problem structure governs out-of-sample generalization, and offers a principled basis for assessing and certifying the reliability of the optimal solution as per constraint satisfaction. Despite its strong theoretical development and wide applicability, no software toolbox has been available to date that enables user-friendly, data-driven convex optimization within the scenario-approach framework. In this paper, we introduce Scen-Opt, an open-source software tool that integrates convex programming with data samples while providing statistical guarantees grounded in scenario theory. Scen-Opt is implemented in Python, supporting data-driven linear, quadratic, and semidefinite programming, and offers a Python-based web application with an intuitive and reactive graphical user interface (GUI) built using modern web technologies. Scen-Opt can be used directly through its online interface or installed locally, accommodating both manual input and data-file uploads (CSV, JSON, TXT, TSV, MAT, Excel, NPY, NPZ, Parquet). Built on a Python backend with a modern JavaScript frontend, Scen-Opt offers a highly user-friendly experience and efficient usability across desktops, laptops, tablets, and mobile devices. In this paper, Scen-Opt is applied to a set of representative benchmarks, demonstrating its practical effectiveness for data-driven convex optimization with guaranteed performance.
Read the full abstract

All authors Ben Wooding, Simone Garatti, Marco C. Campi, Abolfazl Lavaei

Monday 5 October

Mathematical Software · Cross-listed from math.NT

CUDA-MPQS: A GPU-Resident Self-Initializing Quadratic Sieve, and the Factorization of RSA-155

Fabian Januszewski and Christoph Heinrichs

The quadratic sieve (QS) is an irregular, cache-hostile, branch-heavy integer factorization algorithm; prior GPU work accelerates individual stages of it. We present \swname{}, a self-initializing quadratic sieve in which every stage — polynomial initialization, sieving, relation post-processing, large-prime matching, \GFtwo{} matrix construction, Block Wiedemann linear algebra and the square root — executes on the GPU, the CPU restricted to orchestration, setup and I/O, the sieve window 99.9% GPU-busy with no device-wide or event synchronization and no synchronous copy. With it we factored the 512-bit RSA-155 challenge modulus end to end on GPUs: a 64$\times$H100 sieve feeding a single-H100 \GFtwo{} solve, 700.6~GPU-h and 242~kWh in total, 98.4% of it sieving. To our knowledge this is the largest integer factored by the quadratic sieve, exceeding the prior QS record RSA-150 (Buhrow, YAFU, 2025); all comparable general factorizations since 1996, RSA-155's own in 1999 included, used the number field sieve. On a single device we factor RSA-100 in 29.2~s on an H100 and 51.0~s on a consumer RTX~5070~Ti; in a controlled same-modulus measurement, with energy measured on both sides, the H100 is faster than the fastest CPU quadratic sieve, YAFU, on 96 cores of an EPYC~9655 (106–122~s), and than CADO-NFS's SIQS (292~s) and GNFS (263~s). Across four discrete GPUs (a 3.7$\times$ range in CUDA cores) RSA-100 throughput is close to proportional to core count, memory bandwidth being a weaker predictor, and measured cost over seven benchmark factorizations from RSA-100 to RSA-155 grows as $\LN{1/2}{1}^{1.04}$ ($R^2=0.993$). On RSA-150, the identical modulus on which the prior QS record was set, we measure 302.9~GPU-h against that run's 11{,}664 core-hours.
Read the full abstract

Mathematical Software · Cross-listed from cs.NE

A Flexible and Generic Approach for Explainable Landscape Analysis and the pyXla Toolbox

Tony Ombaso, Anna S. Bosman, Arnaud Liefooghe et al. (+2 more)

Landscape analysis has been successfully applied to understand complex optimisation problems, gain insights into algorithm behaviour, and automate algorithm selection and configuration. Although many landscape analysis techniques have been developed over the last decades, it remains difficult for researchers and practitioners to decide which approaches are appropriate and to implement them in practice. Some tools are available, but these are either restricted to particular problem domains (e.g., unconstrained black-box continuous optimisation), or are limited in what they model and measure. In addition, output from landscape analysis is often not easily interpretable, especially when computed landscape features do not correspond with aspects of problems that practitioners are familiar with. In this paper, we introduce a principled approach for explainable landscape analysis (XLA) with an associated Python package called pyXla. The approach is generic in that it applies to problems with different representations (continuous or combinatorial), with single or multiple objectives, with or without constraints. The extent of analysis provided by the XLA framework depends on the data available, with richer analysis offered as additional information is provided by the user. We demonstrate the explainable output produced by pyXla on a selection of hand-crafted problems with diverse landscape characteristics.
Read the full abstract

All authors Tony Ombaso, Anna S. Bosman, Arnaud Liefooghe, Katherine M. Malan, Sébastien Verel

Mathematical Software · Cross-listed from cs.DC

Acceleration of Data Analytics on Heterogeneous Supercloud Systems

Georgios Zacharopoulos, Ilias Bournias and Lukas Cavigelli

Heterogeneous Supercloud systems are transforming data analytics by enabling scalable and efficient task distribution across diverse resources. This paper presents computational models and performance estimation techniques tailored for accelerating Dense Cholesky (CH) Decomposition and Singular Value Decomposition (SVD)-based analytics on Heterogeneous Superclouds. Our ReDSEa tool-chain automates mapping, load balancing, scheduling, parallelism, and computation overlap. Implemented on a heterogeneous system with a Huawei Kunpeng 920 ARM CPU and an Ascend 910 AI accelerator, our LLVM compiler tool-chain employs novel performance models for recursive, iterative, and blocked computations, achieving up to 17x speedup for CH and 88x for SVD over fully optimized 48-core CPU implementations.
Read the full abstract

Mathematical Software · Cross-listed from cs.CV

Robust Local Optimization Done Right

James Pritts and Kevin Köser

RANSAC scoring and local optimization (LO) impose different robustness requirements, motivating the separation of hypothesis selection from refinement. We systematically isolate the effects of robust-loss shape, incorrectly specified inlier scales, and optimization strategy on essential matrix, fundamental matrix, and homography estimation. A profile-marginal score marginalizes the nuisance inlier scale and selects an inlier partition, from which we estimate the scale that sets the LO loss width; this makes LO robust to an inlier scale specified too large, whereas one specified too small degrades selection itself. Refinement needs gradient from correspondences the seed currently rejects: optimizers that reweight from current residuals stay pinned to their seed, whereas methods with broad basins recover strongly perturbed seeds yet degrade accurate score-selected hypotheses, so basin size alone is insufficient to assess RANSAC LO. Joint half-quadratic optimization balances the two and is the most consistent strategy across model classes. An optimizer matched to the profile-marginal score, which never decreases it, does not reach the best accuracy, challenging the prescription that scoring and refinement objectives should match. Composed from these findings, our RANSAC reduces the median essential-matrix pose error of a state-of-the-art RANSAC on PhotoTourism from 2.23 degrees to 1.58 degrees with a correctly specified inlier scale and from 38 degrees to 6.2 degrees when it is grossly misspecified (128x too large).
Read the full abstract

Thursday 1 October

Mathematical Software · Cross-listed from eess.SY

Interactive Power Flow in the Browser

Samuel Talkington, Frederik Geth, Qian Zhang et al. (+2 more)

This paper introduces tellegen, an open source framework for interactive power flow (PF) and optimal power flow (OPF) studies that run in the browser. This provides intuitive and democratized access to power system analysis tools compiled to WebAssembly. A user can drag and drop a case file, click and drag to change a nodal demand or line rating, preview the impacts via sensitivity analysis, and obtain an exact re-solve on release. User case files and results stay entirely on the device: tellegen transmits zero Critical Energy/Electric Infrastructure Information (CEII). The framework comprises a core numerical engine for PF and OPF, reusable browser components, saved studies, and structured WebMCP tools for agentic interaction. We evaluate the transmission OPF solver by comparing objectives with PGLib baselines; the distribution PF solver by comparing voltages and currents with OpenDSS; and the WebAssembly execution times by comparing with PowerModels.jl. On realistic synthetic grids, WebAssembly OPF solves take only 25-43% longer than native binary solves. The implementation shows how an engineer can distribute an executable numerical study as a URL, reducing installation and hosting requirements while keeping case data local.
Read the full abstract

All authors Samuel Talkington, Frederik Geth, Qian Zhang, Le Xie, Skyler Liu

Tuesday 29 September

Mathematical Software · Cross-listed from math.NA

Running error bounds in finite element kernels

Michal Habera, Paul T. Kühner, Matteo Croci et al. (+1 more)

Rounding errors in finite element computations can lead to a complete loss of accuracy, stalled convergence, and incorrect results. Moreover, the effects of rounding errors accumulated within automatically generated and compiled kernels are difficult to analyze a priori. We present the first software framework for automated rounding error estimation within finite element kernels. The proposed methodology is based on an a posteriori technique called Running Error Analysis (REA), where a forward error estimate is automatically computed concurrently with the value. An open-source implementation is provided for the FEniCS Form Compiler (FFCx), based on a C++ backend for generating type-generic templated kernels over a custom arithmetic type that tracks both the value and its error estimate. We demonstrate REA on two examples. First, we use it to detect catastrophic cancellation in the assembly of a Neo-Hooke hyperelastic model in the small deformation regime. A series expansion circumvents the cancellation problem and the computed error estimates show this. Second, we study the assembly of the Laplace operator on a near-degenerate mesh. We demonstrate that our REA implementation typically incurs only 2-4x performance overhead. Applications of this work include robust reduced-precision computations in embedded systems, numerical debugging of new, possibly ill-conditioned or unstable PDE formulations, and guiding the design of mixed-precision kernels.
Read the full abstract

All authors Michal Habera, Paul T. Kühner, Matteo Croci, Andreas Zilian

Mathematical Software · Cross-listed from math.NA

Mixed-Precision Computing for Scientific Discovery: Formats, Co-Design, and Responsible Approximation

Emmanuel Agullo, Hartwig Anzt, Daniel Bauer et al. (+37 more)

Reduced and mixed precision have moved from a niche optimization to a central design axis in scientific computing and engineering, driven by energy constraints, heterogeneous accelerators, and the convergence of simulation and machine learning. This paper organizes the landscape around seven coupled themes—number formats, floating-point emulation, emerging architectures, hardware/software co-design, relation to other approximations, software design, and precision as a multilevel resource —and, for each theme, synthesizes the state of the art, future directions, and open questions. We emphasize energy per trusted solution as the core objective, and we frame \say{recklessly responsible} computing as a pragmatic doctrine: exploit low precision aggressively, but with systematic detection, escalation, and certification pathways.
Read the full abstract

All authors Emmanuel Agullo, Hartwig Anzt, Daniel Bauer, David Bindel, Alfredo Buttari, Alexandru Calotoiu, Erin Claire Carson, Pasqua D'Ambra, Ieva Daužickaitė, James W. Demmel, Jack Dongarra, Iain Duff, Massimiliano Fasi, Dominik Göddeke, Stef Graillat, Laslo Hunhold, Roman Iakymchuk, Fabienne Jézéquel, Nils Kohl, Harald Köstler, Jakub Kružík, Julien Langou, Xiaoye Sherry Li, Hatem Ltaief, Piotr Luszczek, Yuxin Ma, Theo Mary, Mantas Mikaitis, Hiroyuki Ootomo, Daniel Osei-Kuffuor, Enrique S. Quintana-Ortí, Ulrich Rüde, Jennifer Scott, John Shalf, Linda Stals, Rasmus Tamstorf, Stefan Turek, Petr Vacek, Bastien Vieublé, Rio Yokota

Monday 28 September

Mathematical Software · Cross-listed from cs.LG

E3J: An Efficient and Open-Source Backend for Euclidean Equivariant Operations on GPU and TPU

Olivier Peltre, Armand Picard, Adrien Pichard et al. (+6 more)

We present e3j, a fast Euclid-equivariance backend for geometric deep learning applications with JAX bindings for GPU and TPU. Leveraging both optimized CUDA and Pallas kernels and algorithmic improvements, the library achieves state-of-the-art throughput and runtime on both forward and backward paths. On a machine learning interatomic potential (MLIP) use case, it outperforms established backends, measuring up to 34% speed-up over cuEquivariance on water box NPT simulation using MACE, while remaining fully open source. E3j achieves over 80% efficiency over the H100 maximum memory bandwidth on tensor product operations, and in many cases more than doubles throughput of message passing convolutions forward compared to previously available backends. In addition, with the release of dedicated Pallas TPU kernel, e3j opens the possibility of large scale equivariant deep learning workloads on TPU architectures, which has so far been difficult to achieve. Our benchmarks show that e3j also achieves over 80% of a TPUv6e memory bandwidth, up to one order of magnitude more than e3nn-jax. The library is available on GitHub, PyPI and is released under an open source Apache 2.0 license.
Read the full abstract

All authors Olivier Peltre, Armand Picard, Adrien Pichard, Miguel Bragança, Luca Giacomoni, Valentin Heyraud, Zachary Weller-Davies, Christoph Brunken, Jules Tilly

Mathematical Software

A performance enhancement of the Payne-Hanek range reduction algorithm

Tue Ly

Range reduction plays a crucial role in the accuracy and performance of evaluating trigonometric functions, and is often the primary bottleneck for large floating-point inputs. While fast algorithms such as Cody–Waite work efficiently over narrow intervals, the Payne–Hanek algorithm remains the standard technique for accurate reduction across large floating-point inputs. However, many existing implementations of Payne–Hanek suffer from high latency due to heavy branching, conversion overheads, and the use of multi-word integer arithmetic, which hinders SIMD vectorization. In this paper, we analyze and present a branch-free variation of the Payne–Hanek algorithm using only floating-point arithmetic. Our method operates directly over large double-precision inputs ($|x| \ge 2^{16}$) and is well suited to hardware with FMA instructions. We formulate the precision constraints in terms of a truncation error budget, construct a compact lookup table indexed by the input exponent, and prove that the scaled reduced argument has absolute error below $2^{-110}$ and relative error below $2^{-60}$ for every input, including worst cases. The same routine can serve both as the complete range reduction of a single-stage implementation and as the fast path of a correctly rounded one, achieving higher throughput than existing implementations and lower latency than those returning a double-double reduced argument. The algorithm is currently implemented in the LLVM libc project.
Read the full abstract

Saturday 26 September

Mathematical Software · Cross-listed from cs.NI

Analyzing 10 Petabit/s Network Data with Accelerated Associative (Token) Arrays

Jeremy Kepner, Hayden Jananthan, LaToya Anderson et al. (+19 more)

As networks expand and become an ever more critical infrastructure to modern society the need to analyze these networks with the highest regard for privacy is essential to ensure their proper function. Depending on the level of the network layer to be analyzed, sources and destinations can be any combination of physical, logical, or persona/agentic endpoints, which requires the ability to handle diverse data. Invaluable to these analyses are mathematical tools that enable sophisticated mathematical algorithms to be expressed succinctly while achieving scalable vertical (within a compute node), horizontal (across compute nodes), and temporal (over different generations of hardware) performance. Associative (token) array mathematics and corresponding libraries is one approach that can meet these requirements. Accelerating these libraries with GPUs enables the analysis of the largest networks. The MIT/IEEE/Amazon Anonymized Network Sensing Graph Challenge provides a venue for highlighting the applicability of accelerated associative arrays for these types of problems. The D4M associative library has been implemented in a number of languages. This work benchmarks a prototype Matlab D4M GPU accelerated implementation of the Anonymized Network Sensing challenge across a wide range of CPU and GPU hardware. Scalable performance is demonstrated within and across CPU cores, CPU nodes, and GPU nodes. Horizontal scaling across multiple nodes was linear. Running on hundreds of GPU nodes simultaneously achieved a sustained processing rate sufficient to potentially analyze a 10 Petabit/s network.
Read the full abstract

All authors Jeremy Kepner, Hayden Jananthan, LaToya Anderson, William Arcand, David Bestor, William Bergeron, Chansup Byun, Alex Bonn, Daniel Burrill, Vijay Gadepally, Michael Houle, Matthew Hubbell, Michael Jones, Piotr Luszczek, Peter Michaleas, Lauren Milechin, Julie Mullen, Andrew Prout, Albert Reuther, Antonio Rosa, Charles Yee, Alex Pentland

Friday 25 September

Mathematical Software

DistributedDesignOptimizer: A modular Python framework for Setup, Execution and Processing of Distributed Design Optimization

Sebastian Ellmaier, Albert J. de Wit, Akilesh Raveendran et al. (+3 more)

Multidisciplinary Design Optimization (MDO) enables the application of optimization algorithms to complex engineered systems, ranging from aerospace and automotive to robotics and microelectronics. Such multi-component systems typically consist of coupled subsystems, each characterized by its own design variables, constraints and objectives. Without adequate coordination of these couplings, subsystems may be optimal in isolation while the overall system remains suboptimal or even infeasible. Distributed design optimization coordinates coupled subsystem optimization problems while allowing them to retain control over their local design variables. Although promising coordination methods exist, their application is hindered by the lack of a comprehensive software framework which supports intuitive problem definition, provides suitable distributed coordination algorithms, is modular and extensible, enables interactive (post-)processing, and supports the distributed computation of subsystem executions. This work derives requirements for such software, assesses existing frameworks against them, and introduces DistributedDesignOptimizer, an open-source Python framework to fill this gap.
Read the full abstract

All authors Sebastian Ellmaier, Albert J. de Wit, Akilesh Raveendran, Marc-Eric Vogt, Thomas Bäck, Anna V. Kononova

Mathematical Software

KnottedGraph: Scalable knotted-graph topology for scientific and mathematical discovery

Hakan Akgün, Xianquan Yan, Kehan Liu et al. (+2 more)

Scientific data span heterogeneous structures, including coordinates, networks, surfaces, volumes and fields, yet their topology can be quantified within a common framework through graph connectivity, cycle structure, genus and spatial embedding. Graph- and homology-based summaries do not determine spatial embedding, while standard knot and link polynomials require extensions to accommodate branching graphs. Here, we introduce KnottedGraph, a computational framework that converts such scientific representations to knotted graphs that retain graph connectivity and spatial embedding together. It constructs projected diagrams and PD codes, enabling various topological analyses, including Yamada-polynomial evaluation for topological classification. For scalable exact evaluation, it combines partial resolutions that leave the same unresolved connections and optimizes their processing order; the resulting algorithm is verified against published topological invariants of knotted graphs with up to 500 crossings. This scalability enables us to introduce an LLM-assisted mathematical-discovery methodology, in which computational topological data generated across knotted-graph families are used to identify candidate closed-form formulas. With this approach, we identify analytical Yamada-polynomials for generic graph motif families exhibiting Abelian and non-Abelian word sequences. Together, these scalable capabilities make knotted-graph topology computationally accessible across scientific domains, enabling large-scale classification and introducing a route from topological data to LLM-assisted AI4Math discovery.
Read the full abstract

All authors Hakan Akgün, Xianquan Yan, Kehan Liu, Zhaoyun Chen, Ching Hua Lee

Thursday 24 September

Mathematical Software · Cross-listed from eess.SY

Practical Algebraic Parameter Estimation for Noisy Data via Gaussian Process Regression

Oren Bassik, Alexander Demin and Alexey Ovchinnikov

Parameter estimation for ordinary differential equation (ODE) models is a fundamental task that is often complicated by the limitations of conventional optimization-based methods. In theory, differential-algebraic approaches offer an appealing alternative: they reduce the problem to polynomial system solving and do not require user-supplied initial guesses for parameter values. In practice, however, algebraic methods have been limited by their sensitivity to measurement noise, because they require accurate derivatives of observed outputs. In this work, we integrate Gaussian Process Regression (GPR) into the differential-algebraic method and derive a first-order error analysis in terms of noise level and algebraic sensitivity. We evaluate the method across several noise levels on a benchmark of 25 dynamical systems arising in applications including mechanical engineering and systems biology. The proposed method achieves the highest aggregate performance among the methods considered, recovering all sought parameter values and initial conditions to within 10% relative error in 88.5% of runs. These results demonstrate that robust derivative estimation can make differential-algebraic parameter estimation practical for dense, noisy synthetic data while retaining key advantages of the algebraic formulation.
Read the full abstract

Wednesday 23 September

Mathematical Software · Cross-listed from cs.LG

hyperbolix: Hyperbolic Deep Learning in JAX

Timo Klein, Thomas Lang, Yllka Velaj et al. (+1 more)

We present hyperbolix, an open-source library for hyperbolic deep learning in JAX, built on Flax NNX. To our knowledge, it is the first comprehensive, general-purpose hyperbolic deep learning library in JAX. It includes six manifolds with a common interface: Euclidean space, the Poincaré ball, the hyperboloid, the $\kappa$-stereographic model, mixed-curvature product spaces, and the proper velocity space. We implement layer families that cover linear layers, convolutions, attention, normalization, positional encoding, regression, and vector quantization. These building blocks span methods ranging from Ganea's original hyperbolic neural networks to recent fully hyperbolic architectures such as Hypformer and Lorentzian ResNet. Additionally, hyperbolix contains Riemannian optimizers implemented as optax transformations, wrapped distributions, and hyperbolic dimensionality-reduction techniques. Its API uses idiomatic JAX: Manifolds are stateless, with curvature being passed at call time, while manifold operations act on single points, with vmap enabling batch operations. The precision of every checked operation is tested against a closed-form NumPy/SciPy transcription from the source paper or a finite difference, for both float32 and float64. On the hyperboloid, standard formulas for two-point operations, such as the distance, lose precision far from the origin, because they subtract two large, nearly equal terms. hyperbolix replaces these subtractions with cancellation-free formulas that stay accurate in float32 at distances where prior implementations return NaN. hyperbolix is available under the MIT license at https://github.com/timoklein/hyperbolix .
Read the full abstract

All authors Timo Klein, Thomas Lang, Yllka Velaj, Sebastian Tschiatschek

Mathematical Software · Cross-listed from cs.RO

GLASS: Architecture-Tuned, Composable, Device-Side Linear Algebra for Edge Robotics and Beyond

Brian Plancher

GPU robotics lacks the reusable numerical infrastructure of mature CPU stacks, instead relying on compiler frameworks that introduce overhead or repeatedly reimplementing numerical libraries. To address this, we introduce GLASS (GPU Linear Algebra Simple Subroutines), a header-only CUDA C++ library that provides thread-, warp-, block-, and NVIDIA-backed implementations of robotics-scale linear algebra and geometric computations under one composable device API. GLASS treats implementation choice, execution scope, and launch packing as architecture-specific placement decisions determined by offline measurement and resolved statically at compile time. This is critical as the best and worst placements differ by a median of 4.9x (max 81x), with 145 of 396 recommended placements changing between a Jetson AGX Orin and an RTX 5090, and 162 of 396 versus an AGX Xavier. These stakes are highest at the edge as GLASS's advantage over the best of PyTorch and JAX is as much as 73x on the Orin versus 12x on the RTX 5090. GLASS is released open source with independent numerical oracles and source-bound local-GPU test attestation. Finally, integrating GLASS with published robotics systems both exposed a pre-existing numerical bug and improved embedded runtimes by up to 1.5x.
Read the full abstract

Monday 21 September

Mathematical Software

AWE: Adaptive Weight Encoding for Exact Integer Matrix Products with Fewer GEMMs on FP4 Tensor Cores

Shun-ichiro Hayashi, Daichi Mukunoki, Tetsuya Hoshino et al. (+1 more)

Emulation of high-accuracy floating-point matrix multiplication, as in the Ozaki scheme, splits the inputs into low-precision components and multiplies them pairwise. These products must be error-free, and each is an integer matrix product times a scale factor. FP4 Tensor Cores are the fastest on the NVIDIA B200 and B300 but cannot hold INT8 operands. The FP4 values scaled by 2 form the set $S = \{0, \pm1, \pm2, \pm3, \pm4, \pm6, \pm8, \pm12\}$, which contains every residue modulo 13, so with carries any integer splits into base-13 digits that FP4 can store. Prior work splits each INT8 operand into 3 such digits (limbs) with weights $(1, 13, 169)$ and multiplies them pairwise, 9 FP4 matrix multiplications (GEMMs) for INT8$\times$INT8. The classical ways to reduce products, such as the Karatsuba and Toom–Cook methods, do not apply as they stand: sums of limbs reach $\pm 24$ and leave $S$. This paper asks how many FP4 GEMMs are needed for one integer matrix product. We propose Adaptive Weight Encoding (AWE): the limbs take freely chosen integer weights, the stored planes are linear combinations of limbs, and the exact product is the sum of the FP4 GEMMs scaled by reconstruction coefficients. For each input range, we searched these choices for encodings with fewer products and found INT8$\times$INT8 in 6 products and INT4$\times$INT8 in 4. The formulation also holds modulo $m$, which covers the residue number systems of Ozaki scheme II: for the FP64 significand, the 75 products of prior work are reduced to 59. The boundary in product count between encodings with and without residues lies near input width 15. We release the encodings found.
Read the full abstract

All authors Shun-ichiro Hayashi, Daichi Mukunoki, Tetsuya Hoshino, Takahiro Katagiri