Papers

A smooth 3D function can be efficiently and accurately approximated by a 3D lookup table of significantly lower resolution combined with tetrahedral interpolation. The function is sampled once on a coarse grid in a preprocessing step. For each input, the enclosing grid cube is located, one of the 6 tetrahedra in this cube is selected, and the values at its 4 corners are combined by a weighted sum. For RGB to CMYK color conversion, a particularly frequent application of this, a table with $n = 16$ cells per axis ($17^3$ nodes, 16-bit values) requires 38\,KiB. A full 8-bit table requires 64\,MiB. As shown in an experimental evaluation, the mean error of the $n = 16$ table is 0.08\,% ink per channel. For 8-bit input, the lookup uses integer arithmetic only (Section~\ref{sec:integer}). It requires no floating point and no division, and all intermediate values fit in 32 bits. Floating point is needed only to build the table, once and offline. The integer result differs from tetrahedral interpolation with exact fractions by at most 1 unit of the 16-bit output, and white and black are reproduced exactly. This makes the approach suitable for low-power devices and real-time applications where computational resources are limited.

We present reference implementations of four published primal heuristics for mixed-integer programming inside the open-source solver HiGHS (Huangfu and Hall, 2018): Feasibility Jump (Luteberget and Sartor, 2023), fix-propagate-repair (FPR; Salvagnin, Roberti and Fischetti, 2025) with its LP-guided dive-time variant, LocalMIP (Lin, Zou and Cai, 2024), and the first-order-LP feasibility pump Scylla (Mexi et al., 2025). They share one integration interface with per-heuristic budgets and patiences, and each runs in parallel on independently seeded workers. The evaluation is on the mipfeas benchmark (Bussieck and Dirkse, 2026), 233 instances at 600 s, against a separately built unpatched HiGHS on the same machine. The primal integral improves by 23.7% over the full set and by 12.9% on the instances not used for tuning: the patched solver reaches a good solution sooner. At the time limit it is slightly ahead. It finds a solution on 4 instances where the unpatched solver finds none, against 1 the other way, and its final gap is smaller, significantly so only over the full set. Scylla and the LP-guided dive produced no accepted incumbent in any full-limit run that carried them and ship disabled; their implementations remain. The implementation is released as open source.

We propose GPU algorithms for high-fidelity simulation of incompressible turbulent flows. Discretization in space is performed with H(div)-conforming high-order Raviart-Thomas finite elements for the velocity and an $L^2$-conforming discontinuous Galerkin approximation for the pressure. In time, a consistent splitting scheme based on higher-order BDF time stepping is used, with convection treated explicitly. In this scheme, a pressure Poisson equation and a symmetric reaction-diffusion-type equation for the velocity need to be solved in each time step. We develop a solution framework with fast matrix-free operator evaluation for all ingredients, combined with multigrid solvers for the Poisson problem, and propose a robust mixed-precision algorithmic framework. A key to mixed-precision efficiency is a least-squares projection to generate accurate initial guesses for the iterative linear solvers, enabling us to work with relative residual tolerances of $10^{-3}$. In this regime, running the solvers entirely in single precision leads to almost no change in overall iteration counts and maintains the crucial turbulence statistics, while showing up to $1.7\times$ speedup over pure double-precision simulations.

Parallel numerical libraries such as PETSc are widely used in science and engineering applications where wrong results can have costly consequences. Despite this, numerical libraries are rarely formally verified. One of the challenges is the need for an expert to hand-write a specification and manually apply a verification tool, often requiring the development of a harness or driver, all of which can contain additional bugs that lead to false positives or false negatives during verification. With recent advancements in large language models (LLMs), it is tempting to generate such drivers and reference models automatically, but one-shot generation based on a simple prompt is brittle and leads to additional unverified code that needs to be audited. In this paper, we present an approach to use LLMs in a limited setting to generate a small, human-certifiable ACSL contract from the function's documentation, combined with a deterministic toolchain that supports a restricted ACSL profile and generates a driver that uses the CIVL verifier to check the implementation against the contract and, when available, an existing reference model. We demonstrate the pipeline end-to-end on three PETSc functions: MatAXPY (reference model already exists), MatAYPX (no reference model, so the certified contract is the sole oracle), and the non-compressing mode of MatFilter (no reference model, with more complex, conditional behavior). With this pipeline, we were able to discover a bug in PETSc's MatAYPX function that was previously undiscovered and had been present in the code since 1997.

Background and Objectives: Spatial mapping retains the spatial distribution of radiomic features, but computational cost and fragmented software limit its use. We developed Flash-Radiomics with scalar extraction and spatial mapping, a central processing unit (CPU) backend, a hybrid Compute Unified Device Architecture (CUDA) backend, consistent feature names, and Hierarchical Data Format version 5 (HDF5) storage. Methods: We evaluated Image Biomarker Standardisation Initiative (IBSI) compliance, CPU-CUDA concordance, and end-to-end processing time. Compliance testing included 825 chapter 1 (IBSI-1) tests covering 165 high-consensus features and 323 chapter 2 (IBSI-2) tests with numerical references. Concordance testing included 1,148 scalar pairs and 93 spatial-map pairs. End-to-end processing time was measured five times per input volume of interest (VOI) size. Comparisons included the Medical Image Radiomics Processor (MIRP) and PyRadiomics for 102 shared scalar features and PyRadiomics for 93 shared spatial maps. Results: Both backends passed all 1,148 IBSI tests, and all paired results were concordant. At the largest scalar input, CPU required 76.343 s and hybrid CUDA 81.915 s; CPU was 4.7 times faster than MIRP and 190.7 times faster than PyRadiomics. At the largest spatial input completed by both backends, hybrid CUDA reduced processing time by 68.6% relative to CPU (79.280 versus 252.791 s). At PyRadiomics' largest completed spatial input, hybrid CUDA was 84.9 times faster. Conclusions: Flash-Radiomics unified standardized scalar extraction, spatial mapping, concordant CPU-CUDA results, and HDF5 storage. CPU processing time was similar or shorter for scalar extraction, whereas hybrid CUDA was faster for spatial mapping under the tested conditions.

Reaching compressible flow has usually forced lattice Boltzmann methods to abandon the compact stencil or the strict locality that make them efficient. We give up neither, solving the compressible Navier–Stokes–Fourier equations with a scheme that transports only the conserved mass, momentum and energy. The viscous stress and heat flux depend on gradients of the velocity and temperature. To capture them, existing compressible schemes go beyond a single lattice of the conserved fields. They enlarge the velocity set, carry the stress and heat flux as extra transported fields, place the energy on a separate grid, or give up exact streaming for an off-lattice advection. We instead recover the gradients from the non-equilibrium part of the distributions already present, inside the collision. No field beyond the conserved state is transported, no neighbour is read, streaming stays exact, and every computation runs in single precision. The method carries five fields where a transported-flux scheme carries fourteen, at a fraction of the memory and $4.9$ times the throughput. Verified against the exact Sod and Becker shock solutions and a supersonic Taylor–Green vortex, it tracks the direct-numerical-simulation reference on the kinetic energy and stays within the reference-solver spread on both dissipation rates.

We introduce OpenQARP, the Open Quantum Application Research Package: an open-source Python framework for quantum application research, built on a compiled C++ core. In OpenQARP, an application is assembled from an algorithm library that orchestrates three interchangeable layers: blocks describe circuits, primitives describe what to extract from them, and engines describe how they run. The release spans near-term through fault-tolerant methods, targeting problems from electronic structure to combinatorial optimization, with device-aware compilation, circuit cutting, and resource estimation over a single set of numerical conventions that the code must follow. Every numerical feature is checked against an independent oracle, and every published timing carries a correctness check. Against OpenFermion, Qiskit/Aer, PennyLane/Lightning, Qulacs, qsim, and pytket, OpenQARP is one to two orders of magnitude faster on operator algebra and among the fastest on simulation. It matches mature compilers on two-qubit gate count, and expresses a complete algorithm in the lines a framework needs rather than the plumbing a primitive stack demands.

The array programming paradigm applies operations to entire arrays rather than individual elements, eliminating explicit loops and abstracting many low-level details of computation. Two influential approaches have shaped array programming: NumPy's broadcasting and Kenneth Iverson's array-oriented notation, with the latter first introduced in APL and later developed in the J language. Although NumPy's array programming model was influenced by Iverson's work, the two systems formalize array operations differently: NumPy primarily through array shapes and broadcasting, and J through verb rank. This paper establishes a formal correspondence between NumPy broadcasting and J rank and presents an implementation of NumPy broadcasting in J. The results provide a formal connection between two influential approaches to array-oriented computation and show how they address the common problem of applying operations to arrays with differing shapes without explicit loops.

Matrix multipliers available on recent GPUs do not conform with the IEEE 754 floating point standard. Features of matrix multipliers differ across vendors and architectures of the same vendor, such as accumulator width, rounding behaviour, normalisation points, intermediate underflow and overflow logic, the handling of subnormals, and the treatment of special inputs. As a result, reproducibility of small matrix multiplier results across devices is not possible and cannot be achieved by software control. Implementation details of matrix multipliers are not documented, making it difficult to interpret discrepancies in the computed results. We characterise the numerical behaviour of matrix multipliers across three AMD GPU architectures: CDNA 1, CDNA 2, and CDNA 3, using the MI100, MI210/250, and MI300A/300X GPUs, respectively. We design test vectors to target numerical features for all supported input formats and provide the derivation and the reasoning for why each vector allows to determine a particular numerical feature based on the outputs of the devices. MATLAB-based software models of the matrix multipliers are then developed for each architecture and validated for bit-level reproducibility against hardware using a randomized test suite consisting of 10 million sets of random input vectors. To achieve this, we applied a previously developed technique to iteratively refine the accuracy of the models in a loop, by randomized testing followed by test-refinement until the model matches the hardware for every test case. Finally, as a proof of concept for what experimental research can be done with the models, we have utilised them in two demonstrative numerical applications, quantifying application-level accuracy differences between AMD matrix cores and the NVIDIA tensor cores.

We introduce cvxgenrust, an open-source tool for generating custom Rust code that solves families of parameterized convex optimization problems modeled in CVXPY. cvxgenrust canonicalizes a problem family, extracts affine maps to Clarabel cone-program data, and generates a specialized Rust crate that updates parameters and calls Clarabel natively at runtime. The generated solver can also be exposed to Python and registered as a custom CVXPY solver. Our code generator supports a wide range of convex optimization problems up to semidefinite programs and exponential-cone problems. Numerical experiments show reduced runtime relative to direct CVXPY solves and performance comparable to CVXPYgen on shared problem classes.

LP-Newton solves a linear program through a sequence of nearest-point problems. Given an interior feasible point for an inequality-form LP, we construct the compact hull formed by the origin and its normalized constraint rows. Polar LP-Newton (P-LPN) follows the objective ray to the boundary of this row-polar hull, where it recovers a primal-dual optimum or a recession direction proving unboundedness. For rational row-polar data, we prove an outer-iteration bound quadratic in dimension and linear in binary input length, excluding inner Wolfe work. Restarting Wolfe at every outer iteration discards its terminal corral even though the hull is unchanged and the next target lies on the same ray. Corral-repair P-LPN (CR-P-LPN) instead repairs that corral and uses it to start the next projection. Verification over the full hull preserves P-LPN's projections, outer targets, and LP conclusion in exact arithmetic. Experiments in Julia and MATLAB show where repair helps. After common initialization, median P-LPN times across 180 single-LP tests are 1.62 times the corresponding CR-P-LPN times in Julia and 1.55 times in MATLAB. On three problems with 5,000 to 100,000 rows, the corresponding geometric-mean factors are 2.03 and 1.46. A separate end-to-end study includes initialization on 19 application-derived LPs. CR-P-LPN has lower median total time than the faster cold HiGHS mode, either dual simplex or the interior-point method with crossover, on 18 of the 19 workloads in each language. The geometric means of the ratios of HiGHS time to CR-P-LPN time are 4.93 and 3.06. Component comparisons identify terminal-corral reuse as the useful component, whereas candidate-first ordering alone shows no consistent gain. Repair is not always beneficial: it is slower on Klee-Minty tests, and an author-implemented simplex is faster on the retained Netlib models.

Recent developments in the US knowledge economy have created a significant growth in datacenter loads, with two major consequences for power generation. First, to compensate for growth in load, datacenters are encouraged to bring their own generating units. Second, in search of the remaining pools of dispatchable generation, utilities are increasingly turning to subtransmission and distribution level generating assets. Coupled with increasing loads and the resulting tighter grid conditions, both trends are likely to create a need to commit a large number of localized generating units under grid constraints. We propose an extension of Relax-and- Round Unit Commitment (RRUC) that is capable of committing generating units for larger problems faster than conventional methods, while staying within intertemporal and spatial MVA constraints. We demonstrate the performance of RRUC using synthetic congestible test systems ranging from 100 to 20,000 buses. RRUC consistently finds low cost solutions, independent of the problem size, and its run time increases sub-quadratically in the number of buses. RRUC can solve the 100 bus system in less than a second and the 20,000 bus system 7 minutes. In contrast, a leading state of the art solver cannot find a feasible solution to the 100 bus system in 15 minutes.

Counterfactual regret minimization (CFR) is one of the few large numerical workloads that still runs faster on CPUs than on GPUs. Each iteration sweeps a game tree with up to billions of states in millions of small, interdependent gather and scatter steps issued through a generic tree interface. On a GPU every kernel finishes in microseconds, so kernel launches and framework dispatch dominate the run time, and prior GPU implementations have lost to optimized CPU code. We observe that for a fixed game, everything about a CFR iteration except the numerical values is known before the first iteration runs. We propose GPU-CFR, a compiler and runtime built on this observation. It compiles any game once into static dataflow: flat edge and information-set arrays, precomputed indices, and depth-level batched passes fix the entire operation sequence, and only solver state changes between iterations. Static chance folding, depth-level execution blocks, and a dual-lane reach buffer cut the number of framework operations by up to 18.1x. Because shapes, indices, and buffer addresses never change, CUDA Graph Replay records the iteration once and replays it with a single graph launch. On one A100, across an eight-game suite that spans card games, dice games, and board games, GPU-CFR runs 29.8–80.4x faster than the fastest prior GPU CFR on the same accelerator, and 14–258x faster than LiteEFG, one of the fastest open-source CPU implementations, on the four largest games. The compiled representation carries most of that margin: on eight CPU threads with no accelerator it is already 2.2–51.1x faster than the GPU baseline. On the CPU the optimized path reproduces the reference iterates bitwise, and tree construction and graph capture pay for themselves within the first solve. GPU-CFR beats every CPU and GPU baseline on the mid-to-large games of the suite without changing the update rule.

Agentic search and automated optimization of GPU kernels are powerful tools for large language model inference. Their effectiveness, however, depends not on the sophistication of the search itself, but on the clarity of the optimization problem being solved. We provide an application-first approach that drives a hierarchical code generation tool from operator specifi cation down to GPU instructions, and show that a clearly defined computational model makes the optimization problem tractable.

Fluorescence lifetime imaging (FLI) measures the temporal decay of fluorescence after excitation and provides quantitative information about a fluorophore's local environment and molecular interactions. Depending on the fluorophore and experimental design, lifetime can report changes associated with pH, oxygenation, cellular metabolism, and Forster resonance energy transfer. These properties make FLI useful across microscopy, biophysics, biomedical optics, and preclinical imaging, where the same type of molecular contrast can be studied across different biological scales. FLI measurements, however, are acquired with instruments that record fluorescence in different ways. Intensified charge-coupled device (ICCD) cameras, single-photon avalanche diode (SPAD) arrays, and time-correlated single-photon counting (TCSPC) systems differ in temporal sampling, data organization, detector noise, instrument response, and file format. Lifetime-estimation methods also make different assumptions about the recorded decay, while learning-based approaches require realistic training and validation data for which the underlying parameters are known. PyFLI is an open-source framework for FLI processing and standardized data simulation. It imports measurements from acquisition systems, simulates labeled data under configurable acquisition and noise conditions, and provides complementary approaches for parameter estimation. These include nonlinear least-squares fitting (NLSF), maximum-likelihood estimation (MLE), phasor analysis, rapid lifetime determination (RLD), Laguerre-based estimation, and optional Bayesian and deep-learning inference. CPU and GPU processing support image-scale analysis, while reconstruction, visualization, statistical analysis, and cross-software comparison provide tools for evaluating results. PyFLI includes compressed-sensing reconstruction for single-pixel hyperspectral FLI.

FP8 Ozaki II emulates FP64 matrix multiplication by tensor-core products over a CRT residue system; converting the operands into residue planes (the deconstruction term in the Tensor-Memory Equilibrium model of the companion paper "FP8 is All You Need, Part 1") costs integer-pipe and memory resources before tensor instructions issue. This paper engineers that path; every result is a model projection pending measurement. First, a deconstruction-aware model: on the NVIDIA Rubin GPU the emulated rate reaches the arithmetic roof $P_{\rm FP8}/(3r+1)$ ($\approx 473$ TFLOPS at $r=12$) only within one thread-block cluster; larger outputs are re-split on the fly and held at a floor of $\approx 235$ TFLOPS (half the roof, a ratio of three design integers, not a fit), while real solvers' tall/skinny shapes stay near the crossover, $1.6$-$1.9\times$ over simple deconstruction today. Second, the method: convert-once residue workspaces, an exact two-limb constant-reduction GEMM on integer tensor pipes (or pure-SIMT dp4a), and conversion pipelined behind the MMAs, moving the crossover from $\approx 1211$ to $\approx 480$-$730$. Third, modulus co-design: all-byte and hybrid sets, two supply bounds and a carry-corrected E4M3 split of tail moduli. Fourth and central, the closed-form floor names its hardware escape, and the prize is Rubin's: a stream-side residue-conversion mode on the asynchronous copy path (Option C), a narrow fixed-function block sized as a bill of materials, takes plane formation off the arithmetic pipes and lifts the floor from 235 TFLOPS to the full 473-TFLOPS roof at unchanged cluster reach, about doubling HPL-class FP64 per Rubin GPU, and unbinds conversion-bound sparse kernels. The NVIDIA GB300 GPU, whose 135-TFLOPS roof sits at its own floor, gains little; floor and remedy are Rubin-scale. Application traces ground the analysis; constants are script-checked.

We describe geometric-function-atlas, our open-source Python package for the sharp extremal problems of geometric function theory. We organise it around a catalogue of thirty-nine Ma–Minda starlike generators. From this catalogue we compute exact Taylor coefficients, closed-form Fekete–Szegő constants, exact coefficients of the Ma–Minda extremal function, and admissibility screens. Our verifier answers membership questions for normalised polynomials at three levels of evidence: a floating-point grid screen, an exact sufficient condition decided in rational arithmetic, and a certified interval enclosure at the worst screened point. Every answer names the level at which we obtained it. We ship a checksummed artifact snapshot with three hundred and six coefficient certificates and seven hundred and two directed inclusion radii. Eight reviewed radius lanes carry certificates whose proof chains we replay symbolically, and we re-execute every coefficient certificate through our exact Schur-parameter machinery on request. We emit all results through one versioned envelope that records the method, the evidence status, the assumptions, and the artifact identifiers. Two optional laboratories apply the same discipline to cryptographic S-box metrics and to image-quality metrics. We present our design, state as propositions what each tier establishes, follow one radius lane from screen to replayed certificate, report measured timings, and place our package among symbolic-algebra, rigorous-numerics, and mathematical-database software. We release geometric-function-atlas under the MIT licence on the Python Package Index and at https://github.com/Prasanna28Devadiga/geometric-function-atlas.

The automatic derivation turns a declared system of conservation laws into a lattice Boltzmann scheme, giving each conserved physical quantity a set of q populations whose linear equilibrium embeds the physical flux in their first moment. When the flux depends on gradients of the conserved state, those gradients are supplied by tracking them as additional transported fields. Since lattice Boltzmann is commonly memory-bound, these additional degrees of freedom reduce the achievable throughput. To reclaim it, we reconstruct the gradients from the moment structure of the equilibrium instead of transporting them. To leading order, the first moment of a conserved quantity's non-equilibrium part carries its gradient. Because the reconstructed flux enters its own equilibrium reference, that moment is the image of the gradient under a linear operator built from the diffusive-flux Jacobian. The reconstruction is that operator's algebraic inverse, generic across gradient-form constitutive closures and generated automatically from each declared flux. The resulting scheme carries the conserved quantities alone, reconstructing the required gradients from the populations and forming the fluxes locally. It converges at second order in double precision across advection-diffusion-reaction, Allen-Cahn, Navier-Stokes, resistive magnetohydrodynamics and homogenized compressible Navier-Stokes-Fourier systems, matching the accuracy of gradient tracking at equal resolution. On an NVIDIA RTX A5000 it is up to 3.7 times faster in single precision, the memory-bound kernels reaching up to 97% of peak memory bandwidth.

Consider a parametric ODE control model. A function of the states and parameters is called observable if its value can in principle be reconstructed from input-output data. The observable functions form a field, called the observation field, represented naturally by a set of generators. Even when the model is not fully observable, this field captures the information still accessible from input-output data. We present an algorithm for computing a concise generating set for the observation field of a model with rational dynamics. The algorithm relies on two new results: one allows observable functions to be extracted from the coefficients of repeated Lie derivatives of the outputs, while the other reduces the required orders of differentiation by exploiting identifiable parameter combinations. We implement the resulting algorithm in StructuralIdentifiability$.$jl (https://github.com/SciML/StructuralIdentifiability.jl). For computational efficiency, we employ recent techniques for differential elimination and rational function field simplification. Using models from epidemiology, chemical kinetics, and cancer modeling, we show that the algorithm produces generators with domain-specific interpretations that can inform model analysis and development.

Ideal point estimation is widely used to analyze and visualize political data. However, selecting the corresponding spatial model involves various trade-offs: while model-based approaches such as Item Response Theory (IRT) are based on utility functions rather than optimized for predictive accuracy, most Machine Learning (ML) alternatives struggle to generalize beyond training data when embedding sparse test responses. We introduce IXPLORE, a bounded ideal point estimation algorithm that combines a predictive fit objective with a sparsity-aware likelihood function. On five benchmark datasets spanning surveys, roll calls, and deliberation, this approach surpasses model-based and ML-based algorithms on reconstruction and imputation error - especially for users with sparse responses. Furthermore, we show that non-linear feature transforms can further reduce the reconstruction error while remaining visually interpretable. To quantify uncertainty, IXPLORE applies grid-based posterior inference on a bounded 2D latent space. Available as a Python package on PyPI, IXPLORE offers a flexible framework for constructing bounded, interpretable political maps with fast inference and strong imputation performance.

Modern CPUs increasingly integrate high-throughput matrix engines optimized for low-precision AI workloads, while many scientific computing applications still rely on FP32 and FP64 GEMM to meet their numerical accuracy requirements. This mismatch motivates an algorithmic bridge that uses low-precision matrix products to emulate higher-precision GEMM. This paper presents a CPU-oriented method based on Intel Advanced Matrix Extensions (AMX) and BF16 matrix products. For FP32, each operand is decomposed into three BF16 components and six selected component products are evaluated, targeting FP32-level accuracy relative to oneMKL SGEMM without claiming elementwise or bitwise identity. For FP64 inputs within the supported BF16 exponent range, the method uses a simplified fixed six-slice Ozaki decomposition. Each retained BF16 product is first produced in FP32, then widened and accumulated in FP64. Four product-count settings retain 6, 10, 15, or 21 component products, exposing the accuracy–performance tradeoff relative to oneMKL DGEMM. The implementation combines precomputed packed component buffers, VNNI-packed $B$ panels, and an FP32 tile-resident operand-reuse schedule. On the tested square matrices, AMX-FP32 exceeds oneMKL SGEMM throughput. For AMX-FP64, low-product-count variants can exceed DGEMM at sufficiently large orders, while retaining more products improves accuracy at additional cost.

Constant optimization refines the numerical coefficients of candidate expressions in tree-based genetic programming for symbolic regression. But its per-generation cost has led modern GPU-accelerated frameworks to omit it or restrict it to lightweight forms. We present a GPU-resident, batched Levenberg–Marquardt solver that optimizes constants across a structurally heterogeneous population of expression trees using a fixed number of population-wide CUDA launches per iteration. Reverse-mode automatic differentiation assembles the per-tree Jacobian in one backward sweep, making the dominant per-iteration cost independent of the number of constants per tree, and a double-precision delivery guard guarantees that returned constants are never worse than their initial values. On early-generation populations, the solver sustains up to $5.1{\times}10^{5}$ trees per second on an NVIDIA A100; at a GPU-saturated benchmark configuration it delivers roughly $9.9{\times}$ the throughput of Operon running on a 64-core EPYC 7763, while matching fp64-reference quality. Integrated in-process into EvoGP, the solver enables end-to-end search to recover governing equations on $10$ of $18$ constructed problems versus 0 for stock EvoGP. Our code is at https://github.com/TensorConv/CuSR.

Ordinary differential equations (ODEs) underlie models in science and engineering, and many applications need derivatives of their solutions with respect to parameters. Ensembles of independent trajectories suit graphics processing units (GPUs), but current GPU software forces a trade-off: the fastest ensemble solvers cannot be differentiated in reverse mode at the speed they solve, and the solvers built for differentiation solve more slowly. No single tool has yet offered a reverse-mode gradient at the speed of a fused-kernel solve. We present GRADSOLVE, an open-source JAX library for solving and reverse-mode differentiating low-dimensional ODE ensembles on NVIDIA GPUs. It records the steps an adaptive solver accepts and differentiates a fixed-step replay of them; the returned gradient is the exact discrete adjoint of those steps, the same derivative Diffrax returns by default, obtained more cheaply from a fixed-length chain than from an adaptive loop. It targets ensembles differentiated many times against one recorded mesh, keeps Diffrax as a fallback, and supports explicit and Rosenbrock integrators. Used as a solver, GRADSOLVE's forward-only kernel ran 2.8x faster than DiffEqGPU.jl; used for gradients, once a record exists, it computed them 5.6-14.1x faster than Diffrax's checkpointed adjoint at matched forward-state accuracy across three GPU generations, the advantage narrowing on large ensembles and, on stiff systems, down to parity at tight accuracy. GRADSOLVE is released at https://github.com/ECLIPSE-AI4Science/gradsolve.

Scientific benchmarks are commonly built by domain experts who write tasks and cross-check one another's work, or who adapt existing material from textbooks, published papers, and online resources. These routes can produce strong evaluations, but they require substantial per-item labor. Language models can reduce this repeated work by proposing candidates quickly. The remaining problem is acceptance. We target scientific questions whose answers require computations with specialist software rather than unaided reasoning alone. A candidate is invalid if its script fails or returns a different answer, or trivial if a model answers it without the software. We present ToolGate, which treats every generated item as a proposal and keeps it only if three gates pass. First, an executable solution script must reproduce the proposed answer when run with the scientific software. Second, randomized no-tool screening rejects candidates that models can already solve from the prompt alone. Third, a tool-using agent must solve each survivor within a fixed time limit. We instantiate ToolGate in FEniCSx with 500 generation attempts. The local-verification gate retains 478 candidates. For final reporting, we rescreen this pool after generation: two randomized no-tool screens exclude 222 from the reported pool, and direct GPT-5.5 API calls at medium reasoning (the API default) exclude another 121. Of the remaining 135, a GPT-5.5 Codex CLI agent with access to FEniCSx solves 130; exact deduplication leaves 128 unique protocol survivors. ToolGate turns repeated answer checking and difficulty screening into an auditable process while leaving domain design and final review to experts.

The Time-Dependent Vehicle Routing Problem with Time Windows (TDVRPTW) captures a central difficulty of urban logistics: travel times vary with the time of day, so the cost of a route depends on when it is driven. This paper introduces KAYROS, an open-source solver for the duration-minimization TDVRPTW that is both anytime and exact. It returns improving valid solutions throughout its time budget and can close instances with conditional computational optimality certificates through an integrated branch-price-and-cut component. We extend the composition of continuous arrival-time functions from the literature, previously described only at proof or proposition level, to the left-continuous functions that stepwise benchmarks require. We prove that the composition routine shipped in the checker computes that operation exactly. We analyze when its evaluation and breakpoint normalization are exact in IEEE-754 double-precision arithmetic, and argue that a canonical, epsilon-free checker must define the objective. The anytime layer combines greedy construction, granular time-dependent local search over balanced route trees, and iterated local search, with a fleet-aware descent for a fleet-cost objective. We further propose an anytime evaluation methodology based on a normalized signed primal integral under a pre-committed statistical design. On 212 instances from five families at a one-hour single-threaded budget, KAYROS achieves a pooled anytime score 63.3% lower than the strongest of three available contenders (Timefold, Hexaly, jsprit). All three contrasts are significant under Holm correction. KAYROS improves all 30 high-effort references of Blauth et al. and 10 literature best-known solutions, and publishes 704 such certificates. The solver, benchmarks and campaign data are openly released.

We present a cohesive framework for simulating seismic wave propagation utilizing quantum computing paradigms and their classical tensor network equivalents. We detail a quantum circuit-based formulation for the explicit finite-difference time-domain (FDTD) solution of the two-dimensional acoustic wave equation and map this quantum architecture onto a tensor train representation, namely for Matrix Product State (MPS). The MPS solver enables deterministic simulation of large-scale wavefield dynamics on classical high-performance computing systems. We demonstrate the MPS representation by computing 2D seismic wavefields on the Marmousi model. Our results indicate that the MPS representation is a viable direction for computing and scaling wavefield propagation.

Bilevel optimization provides a natural modeling language for hierarchical decision problems. However, applying existing numerical solvers usually requires substantial manual analysis and reformulation. In this paper, we introduce disciplined bilevel programming (DBLP), a symbolic framework that allows users to specify and solve optimistic bilevel problems in a high-level, human-readable way that is close to the mathematical formulation. For problems with a disciplined nonlinear upper problem and a convex lower problem satisfying the disciplined parameterized programming rules, DBLP automatically canonicalizes the lower problem into conic form and constructs an equivalent single-level reformulation using the conic Karush-Kuhn-Tucker conditions. We relax the resulting complementarity constraint and use a gap continuation procedure to approximately solve a sequence of smooth nonlinear problems. We implement DBLP in the open-source Python package BLVPY, an extension of CVXPY for bilevel programming. We demonstrate the modeling and solution capabilities of BLVPY on a range of bilevel optimization problems from several application domains. The proposed framework and implementation allow users to specify and solve bilevel optimization problems within a few lines of code, without prior expertise in bilevel modeling and numerical optimization.

The complexity of time-domain simulation of modern power systems has increased significantly because converter-based resources introduce control dynamics that must be simulated alongside slower system-level and fast electromagnetic dynamics. The resulting wide range of timescales may force classical time-domain solvers to use small timesteps, complicate the solution of nonlinear equations at each timestep, and reduce solver reliability under strongly nonlinear and multi-timescale transient conditions. This paper introduces RAPTOR, a first-of-its-kind time-domain simulation framework for power system dynamic simulations. At its core, RAPTOR introduces a new integration technique that represents the unknown trajectory of hybrid differential-algebraic equations (DAEs) over a time interval using a physics-informed random-projection neural network (PIRPNN) built from fixed Gaussian radial basis functions (RBFs). A nonlinear solver, e.g., Newton-Raphson, determines how to combine the fixed RBFs so that the resulting trajectory satisfies the governing equations. By combining RBFs with broad and localized shapes, RAPTOR can represent complex multi-timescale behavior over comparatively long time intervals. This enables larger effective simulation advances, fewer overall nonlinear solver iterations, and faster integration over long time horizons. Numerical studies on stiff RMS models, IEEE 9-, 14-, 39-, 57-, and 118-bus RMS simulation benchmarks, and EMT test cases show that RAPTOR accurately finds solutions with substantially fewer sequential simulation advances while offering superior accuracy-runtime performance. With speedups of over 10x in certain cases, these results showcase RAPTOR's potential to outperform long-standing and widely used integration methods, such as the Radau method and the trapezoidal rule.

Forecasting a correlation or covariance matrix is common in risk management and portfolio construction, but evaluating such a forecast correctly is not routine: naive matrix-comparison metrics are not proper scoring rules, walk-forward evaluation windows are easy to overlap with the estimation window in ways that silently leak information, and significance testing on serially dependent forecast-error sequences needs machinery few analysts implement from scratch. corrscore is a Python package that provides matrix-aware implementations of two established proper scoring rules for this setting – the energy score and the variogram score – dispatched across a closed-form tractability spectrum (point, discrete-mixture, and isotropic-Gaussian-mixture forecasts are scored exactly; a general Monte Carlo ensemble falls back to sampling), a geometry-aware variant of the variogram score built from the affine-invariant distance on the correlation manifold, a zero-overlap-by-construction walk-forward backtest harness, and a bundled significance-testing suite (circular block bootstrap, the Diebold-Mariano test, and the Model Confidence Set). We describe the package's design, its point of departure from the existing scoringRules and properscoring packages, and walk through a complete worked example.

Low-rank methods can reduce the memory and computational requirements of deep neural network training in approaches such as GaLore. Randomized range finders offer an attractive alternative to singular value decompositions, particularly when the required rank is determined adaptively from a prescribed approximation tolerance. We introduce a blocked adaptive randomized range finder based on an implicit Householder QR decomposition and an optimized hybrid CPU–GPU implementation. The proposed method avoids explicit reorthogonalization. Numerical experiments show that it preserves orthogonality and approximation accuracy in regimes where block Gram–Schmidt without reorthogonalization becomes unstable. The blocked formulation exposes matrix–matrix operations and enables overlap of CPU panel factorization with GPU updates. On an NVIDIA GH200, the overlapped implementation reduces the runtime for the largest tested matrix from 9.91 seconds on the CPU to 0.407 seconds. The method provides a stable and efficient building block for low-rank approximation on heterogeneous systems with applications in computational science and engineering.