Software

We introduce stale synchronous parallel as a mode of execution in parallel sparse triangular linear system solve and present a general directed-acyclic-graph scheduler capable of producing such schedules. Stale-synchronous-parallel schedules allow the overlap of synchronisation and compute which results in a geometric-mean speed-up of 7-30% of our scheduler, ElasticDivide, over state-of-the-art synchronous scheduler GrowLocal on an ARM machine using 48 cores. On an x86 machine using 48 cores, we report geometric-mean speed-ups of 19-60% over SpMP.

We revisit the I/O complexity of attention in large language models. Given query-key-value matrices Q,K,V∈ℝ^(n× d), and a machine with fast memory size M, the goal is to compute the "attention matrix" A=softmax(Q K ^(T)/√(d)) V with the minimal number of data transfers between fast and slow memory. Existing methods in the literature, most notably FlashAttention and its variants, incur an I/O cost that depends quadratically on n, while a trivial lower bound only requires Ω(nd) I/O's to read the inputs and write the output. In this work, we present a technique for computing attention where the I/O cost only depends almost-linearly on n in most parameter regimes. This is achieved by developing I/O-efficient algorithms inspired by the recent approximate attention framework of Alman and Song. We also prove corresponding lower bounds in each parameter regime to show that our algorithms are indeed close to I/O-optimal.

Linear attention has emerged as a cornerstone for efficient long-context architectures, as evidenced by its integration into state-of-the-art open-source models including Qwen3.5/3.6, Kimi Linear, and RWKV-7. Models that incorporate linear attention layers with the so-called Delta-Rule involve the inversion of triangular matrices as a core sub-routine. This operation often forms a performance bottleneck, and, due to its high-sensitivity to numerical errors, it can significantly deteriorate end-to-end model accuracy if it is not carefully implemented. This work provides a systematic analysis of both direct and iterative triangular inversion algorithms, targeting methods that are rich in matrix products, and, therefore, have the potential to efficiently utilize modern hardware. To that end, our analysis covers a broad spectrum of mathematical and practical aspects, with a heavy focus on numerical stability, computational complexity, and, ultimately, hardware efficiency and practical considerations. We provide a rigorous experimental evaluation to verify these properties in practical scenarios, and in low-precision floating-point representations, highlighting the strengths and limitations of each method. Performance benchmarks on NPUs reveal up to 4.3× speed-up against the state-of-the-art implementations of SGLang for triangular matrix inversion, leading to significant performance improvements on the entire layer level, while maintaining full end-to-end model accuracy.

Property Directed Reachability (PDR) is a commonly used technique for automated hardware model checking, yet efficiently parallelizing it remains a significant challenge. Existing approaches, such as lemma sharing, often suffer from limited scalability as processor counts increase. In this work, we present two novel sharing-based parallelization strategies, preemptive propagation and ARPOS, and compare their performance with classical lemma sharing. To this end, we develop an asynchronous MPI-based message passing framework for the state-of-the-art rIC3 hardware model checker. Experimental results on the 2025 Hardware Model Checking Competition benchmark demonstrate that our preemptive propagation strategy yields a significant performance boost over classical lemma sharing.

We present a fast and scalable implementation of two general optimization algorithms: Simulated Annealing-Replica Exchange (SARE) and Simulated Bifurcation (SB). We particularly focus on the large-scale sparse Ising and Quadratic Unconstrained Binary Optimization (QUBO) problems, which are common formulations for many combinatorial optimization problems. The implementation of both algorithms is expressed as matrix and vector primitives using GraphBLAS. This enables nonblocking and multithreaded operations through the use of different backends. The SARE algorithm is also naturally expressed in an SPMD (Single Program Multiple Data) paradigm, so our implementation is expressed in this way. The multiprocess algorithm can run in a hybrid fashion by simply selecting a multithreaded backend.A big advantage of our approach is portability: it runs without change on different architectures. We compare our implementations with state-of-the-art general solvers from MQLib and dwave on Max-Cut graphs G22 and G81 from the Gset collection. In single-threaded operation, we observe runtime comparable to (SARE) or better than (SB) the reference, while achieving similar solution quality. Both algorithms show good performance and scalability across different architectures: we tested on both ARM and x86 nodes with similar results.We show that the SARE implementation scales well with consistent speedup using many concurrent processes and across multiple nodes.Also the hybrid approach achieves a measurable speedup compared to the single-threaded operation. Also SB shows almost linear speedup using the multithreaded backends, though its scalability is highly dependent on the size of the problem. In both cases the nonblocking backend achieves limited speedup.

We break the symmetry in classical spectral bi-partitioning in order to incentivise the alignment of directed cut edges. We use this to generate acyclic bi-partitions and furthermore topological orders of directed acyclic graphs with superb locality. The new approach outperforms the state-of-the-art Gorder algorithm by up to 17× on total reuse distance and minimum linear arrangement.

Large-scale solvers are critical in scientific and engineering domains, where achieving portable and near-optimal performance remains a challenge. Existing high-level parallel performance models provide coarse-grain communication cost estimates but fail to capture the hierarchical memory behavior and irregular access patterns prevalent in modern sparse solver workloads. Conversely, sparse modeling techniques focus on low-level kernel details and lack support for solver-level reasoning. To bridge this gap, we introduce PaHiS, a hierarchical cost model that represents hardware and algorithms through a multi-level abstraction of their communication, computation, and synchronization characteristics. We integrate this model into a GraphBLAS library and couple it with an automated microbenchmarking pipeline for hardware characterization, enabling solver cost estimation and thread autotuning for algorithms expressed in GraphBLAS. We evaluate the model on three CPU architectures for two solvers, demonstrating i) strong correlation between predicted and measured performance, ii) high-performance thread-level autotuning in ALP GraphBLAS, and iii) reliable guidance for high-level HW-SW co-design decisions such as hardware selection and algorithm comparison.

Approximating the k-th spectral gap Δ_k=|λ_k-λ_(k+1)| and the corresponding midpoint μ_k=(λ_k+λ_(k+1))/(2) of an N× N Hermitian matrix with eigenvalues λ_1≥λ_2≥…≥λ_N, is an important special case of the eigenproblem with numerous applications in science and engineering. In this work, we present a quantum algorithm which approximates these values up to additive error εΔ_k using a logarithmic number of qubits. Notably, in the QRAM model, its total complexity (queries and gates) is bounded by O( (N^2)/(ε^(2)Δ_k^2)polylog( N,(1)/(Δ_k),(1)/(ε),(1)/(δ))), where ε,δ∈(0,1) are the accuracy and the failure probability, respectively. For large gaps Δ_k, this provides a speed-up against the best-known complexities of classical algorithms, namely, O ( N^(ω)polylog ( N,(1)/(Δ_k),(1)/(ε))), where ω≲ 2.371 is the matrix multiplication exponent. A key technical step in the analysis is the preparation of a suitable random initial state, which ultimately allows us to efficiently count the number of eigenvalues that are smaller than a threshold, while maintaining a quadratic complexity in N. In the black-box access model, we also report an Ω(N^2) query lower bound for deciding the existence of a spectral gap in a binary (albeit non-symmetric) matrix.

Optimizing data movements during program executions is essential for achieving high performance in modern computing systems. This has been classically modeled with the Red-Blue Pebble Game and its variants. In existing models, it is typically assumed that the number of red pebbles, i.e., the size of the fast memory, is larger than the maximum in-degree in the computational directed acyclic graph (DAG). Graphs that do not satisfy this constraint need to be first transformed appropriately, which is not a trivial task for general graphs. In this work we propose a Pebble Game variant to model DAGs with arbitrary in-degrees, by allowing partial computations. In the new model, we show that it is NP-complete to decide whether there exists an optimal pebbling strategy with cost k, even for single-level DAGs and when only two words fit in the fast memory. Approximation algorithms for special cases are also outlined.

We study the problem of scheduling a general computational DAG on multiple processors in a 2-level memory hierarchy. This setting is a natural generalization of several prominent models in the literature, and it simultaneously captures workload balancing, communication, and data movement due to cache size limitations. We first analyze the fundamental properties of this problem from a theoretical perspective, such as its computational complexity. We also prove that optimizing parallelization and memory management separately, as done in many applications, can result in a solution that is a linear factor away from the optimum. On the algorithmic side, we discuss a natural technique to represent and solve the problem as an Integer Linear Program (ILP). We develop a holistic scheduling algorithm based on this approach, and we experimentally study its performance and properties on a small benchmark of computational tasks. Our results confirm that the ILP-based method can indeed find considerably better solutions than a baseline which combines classical scheduling algorithms and memory management policies.

We present several algorithms to generate tables for asymmetric numeral systems and prove that they are optimal in terms of discrepancy. In turn, this gives rise to the strongest proven bound on entropy loss. We further give improved theoretical bounds for the entropy loss in tabled asymmetric numeral systems and a brief empirical evaluation of the stream variant.

At the core of the Transformer, the softmax normalizes the attention matrix to be right stochastic. Previous research has shown that this often de-stabilizes training and that enforcing the attention matrix to be doubly stochastic (through Sinkhorn's algorithm) consistently improves performance across different tasks, domains and Transformer flavors. However, Sinkhorn's algorithm is iterative, approximative, non-parametric and thus inflexible w.r.t. the obtained doubly stochastic matrix (DSM). Recently, it has been proven that DSMs can be obtained with a parametric quantum circuit, yielding a novel quantum inductive bias for DSMs with no known classical analogue. Motivated by this, we demonstrate the feasibility of a hybrid classical-quantum doubly stochastic Transformer (QDSFormer) that replaces the softmax in the self-attention layer with a variational quantum circuit. We study the expressive power of the circuit and find that it yields more diverse DSMs that better preserve information than classical operators. Across multiple small-scale object recognition tasks, we find that our QDSFormer consistently surpasses both a standard ViT and other doubly stochastic Transformers. Beyond the Sinkformer, this comparison includes a novel quantum-inspired doubly stochastic Transformer (based on QR decomposition) that can be of independent interest. Our QDSFormer also shows improved training stability and lower performance variation suggesting that it may mitigate the notoriously unstable training of ViTs on small-scale data.

We develop and analyze new scheduling algorithms for solving sparse triangular linear systems (SpTRSV) in parallel. Our approach, which we call barrier list scheduling, produces highly efficient synchronous schedules for the forward- and backward-substitution algorithm. Compared to state-of-the-art baselines HDagg and SpMP, we achieve a 3.24× and 1.45× geometric-mean speed-up, respectively. We achieve this by obtaining an up to 11× geometric-mean reduction in the number of synchronization barriers over HDagg, whilst maintaining a balanced workload, and by applying a matrix reordering step for locality. We show that our improvements are consistent across a variety of input matrices and hardware architectures.

We study an extension of the well-known red-blue pebble game (RBP) with partial computation steps, inspired by the recent work of Sobczyk [23]. While the original RBP assumes that we need to have all the inputs of an operation in fast memory at the same time, in many concrete computations, the inputs can be aggregated one by one into the final output value. These partial computation steps can enable pebbling strategies with much smaller I/O cost, and in settings where such a step-by-step aggregation is possible, this extended red-blue pebble game offers a much more realistic cost model.

The well-studied red-blue pebble game models the execution of an arbitrary computational DAG by a single processor over a two-level memory hierarchy. We present a natural generalization to a multiprocessor setting where each processor has its own limited fast memory, and all processors share unlimited slow memory. To our knowledge, this is the first thorough study that combines pebbling and DAG scheduling problems, capturing the computation of general workloads on multiple processors with memory constraints and communication costs. Our pebbling model enables us to analyze trade-offs between workload balancing, communication and memory limitations, and it captures real-world factors such as superlinear speedups due to parallelization. Our results include upper and lower bounds on the pebbling cost, an analysis of a greedy pebbling strategy, and an extension of NP-hardness results for specific DAG classes from simpler models. For our main technical contribution, we show two inapproximability results that already hold for the long-standing problem of standard red-blue pebbling: (i) the optimal I/O cost cannot be approximated to any finite factor, and (ii) the optimal total cost (I/O+computation) can only be approximated to a limited constant factor, i.e., it does not allow for a polynomial-time approximation scheme. These results also carry over naturally to our multiprocessor pebbling model.

In this work we revisit the arithmetic and bit complexity of Hermitian eigenproblems. Recently, [BGVKS, FOCS 2020] proved that a (non-Hermitian) matrix A can be diagonalized with a randomized algorithm in O(n^(ω)log²(n/ε)) arithmetic operations, where ω≲ 2.371 is the square matrix multiplication exponent, and [Shah, SODA 2025] significantly improved the bit complexity for the Hermitian case. Our main goal is to obtain similar deterministic complexity bounds for various Hermitian eigenproblems. In the Real RAM model, we show that a Hermitian matrix can be diagonalized deterministically in O(n^(ω)log(n)+n²polylog(n/ε)) arithmetic operations, improving the classic deterministic Õ(n³) algorithms, and derandomizing the aforementioned state-of-the-art. The main technical step is a complete, detailed analysis of a well-known divide-and-conquer tridiagonal eigensolver of Gu and Eisenstat [GE95], when accelerated with the Fast Multipole Method, asserting that it can accurately diagonalize a symmetric tridiagonal matrix in nearly-O(n²) operations. In finite precision, we show that an algorithm by Schönhage [Sch72] to reduce a Hermitian matrix to tridiagonal form is stable in the floating point model, using O(log(n/ε)) bits of precision. This leads to a deterministic algorithm to compute all the eigenvalues of a Hermitian matrix in O(n^(ω)ℱ(log(n/ε)) + n²polylog(n/ε)) bit operations, where ℱ(b) ∈ Õ(b) is the bit complexity of a single floating point operation on b bits. This improves the best known Õ(n³) deterministic and O(n^(ω)log²(n/ε)ℱ(log(n/ε))) randomized complexities. We conclude with some other useful subroutines such as computing spectral gaps, condition numbers, and spectral projectors, and with some open problems.

We design and implement parallel prefix sum (scan) algorithms using Ascend AI accelerators. Ascend accelerators feature specialized computing units, the cube units for efficient matrix multiplication and the vector units for optimized vector operations. A key feature of the proposed scan algorithms is their extensive use of matrix multiplications and accumulations enabled by the cube unit.We present a multi-core scan algorithm that fully utilizes both the cube and vector units of Ascend, reaching up to 37.5% of the theoretical memory bandwidth. In addition, our radix sort implementation, which utilizes matrix multiplications for its parallel splits, highlights the potential of matrix engines in enhancing complex operations (3.3× speedup over the vector-only baseline).

We give sharp point-wise bounds in the weight-aspect on fourth moments of modular forms on arithmetic hyperbolic surfaces associated to Eichler orders. Thereby, we strengthen a result of Xia and extend it to co-compact lattices. We realize this fourth moment by constructing a holomorphic theta kernel on G × G × SL_(2), for G an indefinite inner form of SL_2 over ℚ, based on the Bergman kernel, and considering its L^2 -norm in the Weil variable. The constructed theta kernel further gives rise to new elementary theta series for integral quadratic forms of signature (2,2).

Let f be an L^2 -normalized holomorphic newform of weight k on Γ _0(N) ∖ ℍ with N squarefree or, more generally, on any hyperbolic surface Γ ∖ ℍ attached to an Eichler order of squarefree level in an indefinite quaternion algebra over ℚ. Denote by V the hyperbolic volume of said surface. We prove the sup-norm estimate ‖ Im(·)^((k)/(2)) f ‖_(∞) ≪_(ε) (k V)^((1)/(4)+ε) with absolute implied constant. For a cuspidal Maaß newform φ of eigenvalue λ on such a surface, we prove that ‖φ ‖_(∞) ≪_(λ,ε) V^((1)/(4)+ε). We establish analogous estimates in the setting of definite quaternion algebras.

The paper expresses parallel prefix sums using fixed-size matrix multiplications in the Tensor Core Unit computational model. It derives bounds on algorithm depth, execution time and the number of matrix operations as functions of input size, matrix size, latency and available tensor cores.

Recent years have seen stagnating improvements to branch predictor (BP) efficacy and a dearth of fresh ideas in branch predictor design, calling for fresh thinking in this area. This paper argues that looking at BP from the viewpoint of Reinforcement Learning (RL) facilitates systematic reasoning about, and exploration of, BP designs. We describe how to apply the RL formulation to branch predictors, show that existing predictors can be succinctly expressed in this formulation, and study two RL-based variants of conventional BPs.

This paper describes our winning solution for the ECML-PKDD ChAT Discovery Challenge 2020. We show that whether or not a Twitch user has subscribed to a channel can be well predicted by modeling user activity with boosting trees. We introduce the connection between target-encodings and boosting trees in the context of high cardinality categoricals and find that modeling user activity is more powerful then direct modeling of content when encoded properly and combined with a suitable optimization approach.

We present HiCR, a model to represent the semantics of distributed heterogeneous applications and runtime systems. The model describes a minimal set of abstract operations to enable hardware topology discovery, kernel execution, memory management, communication, and instance management, without prescribing any implementation decisions. The goal of the model is to enable execution in current and future systems without the need for significant refactoring, while also being able to serve any governing parallel programming paradigm. In terms of software abstraction, HiCR is naturally located between distributed heterogeneous systems and runtime systems. We coin the phrase Runtime Support Layer for this level of abstraction. We explain how the model's components and operations are realized by a plugin-based approach that takes care of device-specific implementation details, and present examples of HiCR-based applications that operate equally on a diversity of platforms.

In this work, we revisit prefix sums through the lens of linear algebra. We describe an identity that decomposes triangular all-ones matrices as a sum of two Kronecker products, and apply it to design recursive prefix sum algorithms and circuits. Notably, the proposed family of circuits is the first one that achieves the following three properties simultaneously: (i) zero-deficiency, (ii) constant fan-out per-level, and (iii) depth that is asymptotically strictly smaller than 2log(n) for input length n. As an application, we show how to use these circuits to design quantum adders with 1.893log(n)+O(1) Toffoli depth, O(n) Toffoli gates, and O(n) additional qubits, improving the Toffoli depth and/or Toffoli size of existing constructions.

Specialized computational units that perform small matrix multiplications as primitive operations are typically present in modern AI accelerators. However, these Matrix Multiplication Units (MMUs) are often underutilized for many fundamental deep learning operations besides dense matrix multiplications. Coincidentally, the lack of a rigorous theoretical model of computation for such architectures obstructs algorithmic design. In this work, we propose MMV-RAM, a computational model which judiciously extends the Vector-RAM model with an additional MMU. We provide a detailed theoretical analysis and carefully balance the computational power between the matrix and vector units, guided by the circuit complexity lower bound that parity is not in AC{[0]}. Given MMV-RAM, we proceed to algorithm design, starting with two fundamental parallel operations: segmented scan and sum. By expressing them as compositions of elementary parallel primitives (e.g., seg. sum reduces to: scan, compress, and vector differentiation), we can exploit MMUs to perform speculative blocked computations, ultimately leading to provable theoretical speed-ups against vector-only approaches. These results extend to other ubiquitous AI kernels, including dense matrix product, and sparse matrix-vector product. As a case study, we implemented the proposed algorithms on the Ascend 910B AI accelerator, which contains matrix and vector cores. We evaluate these implementations on synthetic and real-world datasets from various applications, including Large Language Models.

In the prequel, a sharp bound in the level aspect on the fourth moment of Hecke--Maaß forms with an inexplicit (in fact exponential) dependency on the eigenvalue was given. In this paper, we develop further the framework of explicit theta test functions in order to capture the eigenvalue more precisely. We use this to reduce a sharp hybrid fourth moment bound to an intricate counting problem. Unconditionally, we give a hybrid bound, which is sharp in the level aspect and with a slightly larger than convex dependency on the eigenvalue.

This work investigates how a small set of past branch outcomes can explain useful correlations for processor branch prediction. It uses sparse modeling to identify these relationships and proposes a compact prediction mechanism that stores the resulting models. Experiments examine prediction accuracy and the hardware cost of using sparse correlations.