Efficient Constant Optimization for Symbolic Regression with GPU-Accelerated Tree-Based Genetic Programming

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.

PyDoseRT Proton: A GPU Pencil-Beam Engine with a Convolutional Residual-Correction Network for Fast Proton Dose Calculation

Architecture category. Hybrid method: a physics-based analytical pencil-beam (PB) dose engine followed by a 3-D convolutional residual-correction network (RepVGG-U-Net). We addressed the DoseRAD2026 proton dose-prediction task with PyDoseRT Proton, a GPU-accelerated engine implemented in PyTorch and augmented by a learned residual toward Monte Carlo (MC) accuracy. A double-Gaussian PB kernel was calibrated to GATE/Geant4 integrated depth doses in water in two stages: a classical per-energy curve fit, then a gradient-based fit of the full 3-D dose through the PyTorch physics engine as it retains a differentiable execution path for gradient-based optimization of dose-dependent objectives. The engine computes each beamlet on a beam's-eye-view (BEV) lattice with variance-preserving Gaussian splitting, an analytic nuclear halo, and a Fermi-Eyges heterogeneity term, then rotates the result into the patient frame. Additionally, a compact residual U-Net predicts an additive correction in BEV space. It is conditioned on voxelwise material-label embeddings, a discrete energy embedding and spot size. The same model was used for all anatomical sites (thoracic and abdominal). It was trained with a patient-space L1 objective emphasizing the scored high-dose region and multi-scale BEV deep supervision. The submitted CT configuration obtained preliminary-test beamlet MAE 0.0066, image-z IDD distance 0.0025, plan MAE 0.0049, 98.30\% gamma pass rate (1\%/1 mm), and DVH error 0.460.

Learning to Difference: Adaptive Reversible Differencing (AdaRDiff) for Time Series Forecasting

Reliable long-horizon time series forecasting is an important yet difficult problem. Trends and seasonality introduce complex temporal structure that challenges learning-based forecasting models. Differencing, which subtracts nearby past values to remove such structure, is the classical remedy, but its reliance on hand-picked orders and periods has kept it largely absent from recent deep architectures. We propose \textbf{\underline{Ada}}ptive \textbf{\underline{R}}eversible \textbf{\underline{Diff}}erencing \textbf{(AdaRDiff)}, a generalized differencing approach that uses learnable weights to simplify the series through weighted differencing with previous time instants. This yields stabilized residuals on which forecasting is performed, after which the removed components are restored autoregressively to reconstruct the forecast, capturing trend and seasonality jointly through a single operator. This reconstruction admits a closed-form convolutional expression, which parallelizes on GPU and yields up to $33.7\times$ speedup over the naive recurrence. We furthermore rely on a two-phase training schedule that separates temporal structure discovery from reconstruction learning, as suggested by a theoretical analysis of the gradient when using a linear forecasting model. AdaRDiff attains state-of-the-art forecast accuracy across eight benchmarks spanning electricity, weather, traffic, and energy, at negligible parameter cost. Furthermore, it is designed as a plug-and-play module: integrating AdaRDiff improves eight diverse backbones, from linear models to Transformers, in the large majority of cases, by up to $25.9\%$ with a linear backbone and $18.3\%$ with iTransformer.

Tensor-Accelerated Eager Multi-Resolution Grids for Evolving Large-Scale Substrates

In neuroevolution, indirect encoding generates neural network connectivity from a compact genome rather than specifying each connection. ES-HyperNEAT automatically discovers where to place hidden nodes by examining CPPN output patterns: it recursively subdivides space using a quadtree, expanding regions where CPPN outputs show high variance. This adaptive approach discovers network topology without manual substrate specification, extending the fixed-grid HyperNEAT framework built on NEAT. However, the quadtree resists tensorization. Each depth level depends on the parent's variance, forcing sequential evaluation. Different CPPNs produce different subdivision patterns, preventing batching. And variable leaf counts are incompatible with JAX's static shape requirement for JIT compilation. Our prior work confirmed these limits at depths exceeding 5, and a JAX reimplementation of the quadtree yielded only marginal speedup despite batched optimizations, motivating the eager reformulation presented here. We present EMR-HyperNEAT, which evaluates all positions at all resolutions up front, then filters using the same variance criterion: ES-HyperNEAT's subdivide_if(var > $θ$) becomes eval_all(); filter(var > $θ$). This performs more CPPN queries than necessary, but all queries become independent and parallelizable across both cores and population members, reducing complexity from \BigO($4^D$) to \BigO($4^D/P$) across $P$ parallel cores. Recurrent substrate configurations become feasible through a connection type taxonomy. The experiments section validates 12-34$\times$ on-device GPU speedup on XOR at depths 5-7, and empirically higher solve rates across benchmarks.

TorchMorph: CUDA-accelerated Morphological Transforms

Morphological transforms are long-standing tools for shape and mask processing, but the de facto reference implementation in the Python ecosystem, i.e. scipy.ndimage, is CPU-only, single-array, and therefore unusable inside a GPU training loop without an expensive device-to-host round trip. GPU vision libraries built on PyTorch cover a narrow subset of these operators, typically restricted to two spatial dimensions and flat structuring elements. We present TorchMorph, a lightweight PyTorch extension that closes this gap. TorchMorph exposes 22 public operators covering binary morphology, greyscale morphology, exact and approximate distance transforms, and entropy-regularised optimal transport, all implemented as fused CUDA kernels that operate directly on (B, C, Spatial...) CUDA tensors with up to eight spatial dimensions. The API deliberately mirrors scipy.ndimage argument-for-argument, including border modes, structuring-element origins and pre-allocated outputs, so that existing pipelines port with a change of import. We describe the layered architecture and the kernel designs behind each operator family. Against single-threaded CPU references, batched execution reaches up to 1.1e3 times the throughput of scipy.ndimage on greyscale morphology and up to 350x on exact Euclidean distance transforms, while the Sinkhorn solver runs up to 42x faster than POT. Binary and chamfer operators reproduce their SciPy counterparts exactly, and every float-valued operator agrees with the CPU reference to within 1.8e-6 absolute error. TorchMorph is released under the MIT licence at https://intcomp.github.io/tm.

A Feature-Major Codebook for Memory-Efficient Sparse-Binary Self-Organizing Maps: Scaling a MEDLINE Atlas to 1.05 Million Neurons on a Single Consumer GPU

A self-organising map turns a large corpus into a browsable two-dimensional atlas, but building one at MEDLINE scale has been impractical: the best-matching-unit (BMU) search that dominates training is bound by the bandwidth needed to read the codebook every epoch. I show that this bottleneck is largely an artefact of codebook layout. Storing it feature-major with each feature's weights contiguous, W[v.M+i], recasts the search as a tiled sparse-dense product in which every loaded weight column is reused across a tile of samples. Varying only the layout, with implementation, precision and update rule held fixed, accelerates the BMU search by 4.5-8.5x. Because an exact-argmin BMU is invariant to how the codebook is stored, this gain costs nothing: held-out quantisation error agrees with a cuSPARSE baseline to within 0.5% at every map size. Against that baseline the advantage is a crossover rather than a constant: cuSPARSE.SOM is faster at small maps, SparseBin.SOM is 1.5x faster at 128x128 and 2.6x at 256x256, and at 512x512 it is the only one that runs at all on 24 GB. Paired with a radius-independent box-blur update and a convergence-based stopping rule, it trains a converged map over 29.9 million MEDLINE articles in about 72 s at 64x64 on one 24 GB GPU, and accommodates 262,144 neurons (512x512 edges) where every alternative algorithm I tested exceeds memory constraints. On a 141 GB H200 it reaches 1,048,576 neurons (1024x1024 edges) - to my knowledge the largest self-organising map yet reported. Held-out error follows a smooth power law with no elbow across three decades of map size, so the limit on resolution is compute rather than any breakpoint in the data. At matched work the design is ~82x faster than MedSOM, the CUDA implementation behind our earlier MEDLINE atlases and, at 128x128, 621x faster than the best available multicore-CPU library.

On the Choice of Tensor Estimation for Corner Detection, Optical Flow and Denoising

Many image processing methods such as corner detection, optical flow and iterative enhancement make use of image tensors. Generally, these tensors are estimated using the structure tensor. In this work we show that the gradient energy tensor can be used as an alternative to the structure tensor in several cases. We apply the gradient energy tensor to common image problem applications such as corner detection, optical flow and image enhancement. Our experimental results suggest that the gradient energy tensor enables real-time tensor-based image enhancement using the graphical processing unit (GPU) and we obtain 40% increase of frame rate without loss of image quality.

FlashReg: GPU-Accelerated 3-Clique Point Cloud Registration for Real-Time Correspondence-to-Pose Estimation

Graph-based point cloud registration achieves high robustness by identifying geometrically consistent correspondence sets, but constructing second-order compatibility graphs and enumerating candidate cliques remain compute- and memory-intensive. This work presents FlashReg, a GPU-oriented correspondence-to-pose estimator that avoids materializing the dense scored second-order graph. Its Fast First- and Second-Order Graph (FFSOG) construction builds a capacity-bounded sparse second-order graph directly from the binary first-order graph. A dataflow-optimized three-node clique (3-clique) search then selects pivots from compact per-row candidate pools and enumerates triples through sorted sparse-neighborhood intersections. Across indoor and outdoor benchmarks, FlashReg reduces correspondence-to-pose latency by 2--3x relative to TurboReg at comparable registration recall, while using about 50% of its peak allocated tensor memory on an embedded GPU. These results make FlashReg suitable as a high-throughput registration backend within onboard perception pipelines.

Scalable High-Fidelity Macromolecular Docking for GPU-Accelerated Supercomputers

Flexible macromolecular docking offers high-fidelity predictions of biomolecular interactions, but remains prohibitively expensive at scale. Among existing approaches, LightDock leverages Glowworm Swarm Optimization (GSO) for accuracy, yet suffers from limited parallelism, irregular computation, and severe load imbalance, preventing efficient execution on GPU supercomputers. We present SparkleDock, a scalable GSO-based docking framework enabling near-real-time flexible docking. We redesign GSO to expose massive fine-grained parallelism at the glowworm-agent level, and restructure the dominant energy scoring computation into a Tensor Core-compatible formulation, enabling efficient execution of irregular pairwise interactions through structured matrix operations. We further introduce a performance-model-driven scheduling for load balancing and out-of-core scaling across GPUs. SparkleDock achieves 9.7 $\times$ and 18.9 $\times$ speedups over LightDock on single A100 and H100 GPU, and delivers over two orders of magnitude acceleration at scale. On 512 GPUs, it reduces docking time from hours to seconds, enabling large-scale, high-fidelity virtual screening previously impractical with flexible docking.

SoRoMoX: Fast, Differentiable, and Parallelizable Soft Robot Models

Reduced-order models based on Cosserat-rod theory are now well established, and modeling theory is no longer the primary bottleneck in soft-robot control. Their implementations, however, do not support the differentiable, GPU-parallel, and control-oriented workflows that underpin advanced rigid-robotics applications. Here, we fill this gap with SoRoMoX (Soft Robot Models in JAX), a fully numerical, JIT-compilable Python/JAX framework. SoRoMoX implements articulated, Piecewise Constant Strain, and Variable Strain models through a unified, control-ready interface that provides inertia matrices, gravitational and elastic forces, Jacobians, and their derivatives. To our knowledge, it is the first rod/strain-based soft-robot modeling framework that runs directly on GPUs and is end-to-end differentiable with respect to states, inputs, and parameters. Sequential CPU rollouts are up to 18.1x faster than state-of-the-art alternatives, while GPU-parallel rollouts increase throughput by up to 234.6x. This performance enables workflows that were previously impractical or impossible: static-equilibrium system identification with 66% lower marker RMSE; residual-force learning with a further 64% reduction; computed-torque tracking with RMSE reduced by a factor of approximately 500 relative to model-free PD; control-gain optimization with up to 62% lower loss than untuned gains; safety-constrained control using high-order control barrier functions to keep the peak contact force within a prescribed 5 N bound, compared with 33.5 N without the safety constraint; and reinforcement-learning policy training up to 7x faster than a CPU PyElastica discrete-rod baseline through massively parallel rollouts.

Density-Functional Excited-State Gradients and Nonadiabatic Couplings on a Consumer GPU from a Contraction-DAG

Nonadiabatic dynamics needs an excited-state gradient and an interstate nonadiabatic coupling matrix element (NACME) at every nuclear geometry, and a double-hybrid functional's accuracy has been unavailable for the coupling. We report the first analytic derivative NACME for a double-hybrid excited state---deferred in the original hh-TDA method and supplied for hybrids only by Yu \emph{et al.}---derived, with the hole-hole and particle-particle Tamm--Dancoff (\hhTDA/\ppTDA) gradients and NACMEs, as a single reverse-mode transpose of one contraction graph closed under a non-symmetric atomic-orbital-direct $J/K$ kernel. Its double-hybrid excitation energy lowers the vertical-excitation mean absolute deviation from bare-\hhTDA\ $0.86$ to $0.47$~eV and removes the $+0.53\!\rightarrow\!+0.05$~eV over-excitation bias, improving seven of ten states while over-correcting the ionic $ππ^*$ states---the expected perturbative-doubles failure, reported not trimmed. Every coupling is validated to $\sim\!10^{-4}$ against an independent \emph{literal many-electron wavefunction-overlap} oracle that shares no code path with the method, and is physically meaningful at the ammonia $n\!\rightarrow\!σ^*$ \emph{covalent} conical intersection, where the \hhTDA/\ppTDA manifolds recover the $F\!-\!2$ seam and adiabatic linear-response TDDFT gives $τ\!\equiv\!0$ by construction. Gradients, NACMEs, and the double-hybrid coupling all run device-resident and AO-direct through one shared Cholesky-decomposed $J/K$ engine within the 8\,GB of a consumer RTX~4060 (a profile-guided $\sim\!10^2\times$ launch collapse preserving double-precision bit-identity)---placing on a commodity desktop card a correlated excited-state derivative capability that has until now required datacenter hardware.

Learning to Rank Tensor Network Contraction Plans for GPU-Accelerated Quantum Circuit Simulation

Classical simulation remains essential for developing and validating quantum algorithms, but its cost grows rapidly with circuit size. Tensor-network contraction can reduce this cost by exploiting circuit structure, although its efficiency depends strongly on the chosen contraction plan. On GPUs, plans with similar theoretical complexity may perform very differently because execution also depends on parallelism, reduction structure, memory traffic, and contraction geometry. We present a learning-to-rank framework for selecting efficient contraction plans before executing them. Each plan is represented by structural features derived directly from its sequence of pairwise contractions, and gradient-boosted rankers are trained from GPU measurements using listwise and pairwise objectives. We evaluate the resulting models on diverse circuit families, using separate in-distribution and circuit-family-shift test sets, and compare them with random and MinFill-based baselines. The learned rankers generally identify better plans, with the listwise model providing the strongest overall decision quality. We also study backend shift by comparing empirical plan orderings on two GPU architectures and evaluating the source-trained models on the second device without retraining. The rankings remain substantially, though not perfectly, stable across GPUs, and the models retain useful decision quality. These results support Learning to Rank as a practical way to reduce contraction-plan search, while showing that performance remains partly backend dependent.

Accelerating Dynamic Graph Clustering on GPU Architectures with cuGraph

This work addresses community detection in temporal networks through GPU-accelerated extensions of spectral clustering and modularity-based algorithms originally designed for static graphs. Built on the NVIDIA RAPIDS ecosystem, the framework enables the characterization and tracking of communities in snapshot-based dynamic graphs, either by Leiden greedy optimization with multi-GPU support via Dask-based workload distribution, or eigendecomposition of a symmetric Bethe-Hessian operator. Our multislice modularity backend achieves up to roughly three orders of magnitude speedup over the CPU reference under an equal-work budget, depending on graph density and snapshot count, while preserving compatibility with existing graph analytics pipelines. We demonstrate its applicability on real-world and synthetic datasets, facilitating exploratory analysis of structural network properties over time. Such capabilities are relevant across several application domains, such as epidemic spreading, financial systems, cybersecurity, and trajectory and mobility analysis. We release our implementation as free and open-source software, including Python bindings through the NetworkX-Temporal library for ease of use and zero-code acceleration with existing codebases.

TorchDCM: A Unified PyTorch-Native Package for Discrete Choice Modeling

Estimating large and simulation-intensive discrete choice models (DCMs) requires repeated evaluation of utilities, probabilities, derivatives, and simulated likelihoods over many observations, alternatives, and draws. Existing DCM software provides mature econometric workflows, while recent GPU-oriented tools accelerate selected models, leaving a gap between econometric coverage and scalable differentiable computation. We introduce TorchDCM, an open Python package for discrete choice modeling that compiles choice data and model specifications into a unified PyTorch-native likelihood engine for estimation, inference, prediction, and structured reporting on CPU or CUDA devices. The package covers the principal econometric functionality available across Biogeme and Apollo, including multinomial, nested, mixed, ordered, latent-variable, and panel likelihoods. It also supports ragged choice sets, constrained parameters, covariance estimation, willingness-to-pay analysis, elasticities, and extensible likelihood components. We evaluate TorchDCM against seven other estimation packages in aligned synthetic and real-data full-estimation experiments. TorchDCM completes all 45 synthetic cases, runs fastest in every comparable synthetic case, and satisfies the prespecified final-log-likelihood tolerance in every comparison with at least two comparable solutions. More precisely, it reduces median runtime by 89.1%-99.7% relative to Biogeme and Apollo across model-data settings. CUDA provides an additional 12.0-71.0x speedup over single-core TorchDCM. These results establish a scalable and reproducible foundation for econometric estimation and differentiable choice-model development. The open-source package and executed examples are available at https://github.com/mbc96325/torchdcm.

An automatic-differentiation framework for time-lapse electrical resistivity tomography inversion of hydrologic dynamics

Time-lapse electrical resistivity tomography (TL-ERT) provides spatially distributed information on subsurface hydrologic changes. However, inversion of long monitoring sequences is computationally demanding. Modifying the data misfit, regularization, model parameterization, or petrophysical transformation may also require new gradient derivations and separate implementations. Here, we present AD-TLERT, a unified, GPU-accelerated framework for time-lapse ERT inversion based on automatic differentiation. The framework integrates model parameterization, differentiable petrophysical transformations, forward modeling, data misfit, regularization and auxiliary constraints into a single computational chain. Alternative inversion formulations can therefore reuse the same PDE derivative implementation without re-deriving the complete ERT sensitivity for each case. Comparisons with pyGIMLi showed close agreement in the forward responses, gradients, and recovered resistivity models. Under the tested configuration, AD-TLERT achieved an approximately 51-fold speedup. Synthetic experiments showed that inversion choices affect the amplitude, geometry, and temporal behavior of recovered anomalies. By propagating gradients through the embedded petrophysical relationship, AD-TLERT enabled direct water-content inversion and yielded more accurate estimates than post-inversion conversion for the tested model. A field application further demonstrated how ERT, temperature, and soil-moisture observations can be combined to image snowmelt-driven hillslope wetting. AD-TLERT provides an efficient and flexible framework for time-lapse ERT inversion and hydrologic interpretation.

A Fortran General-Purpose Transpiler: Proof of Concept

Fortran has been the cornerstone of high-performance computing for decades and remains unmatched in many domains. Yet the language faces an expertise gap: a new generation of scientists is barely familiar with it, while many experienced Fortran developers are only now transitioning to modern ecosystems such as JAX. This gap often results in "Fython" - Python code written with a Fortran mindset - that fails to leverage modern frameworks. We present FGPT, a transpiler framework designed to bridge this gap. It provides a systematic pipeline that transpiles Fortran into GPU-adapted Fortran, auto-differentiable Fortran via Tapenade, or NumPy and JAX scripts. Its architecture comprises three stages: (i) a frontend that parses Fortran and extracts target procedures along with all their dependencies; (ii) a middle-end that lowers the code into an intermediate representation, then into GPU-adapted or auto-differentiable Fortran, or a NumPy class; and (iii) a backend that transforms NumPy scripts into JAX modules optimized for GPU acceleration and automatic differentiation. Large language models fail when applied to the scale of community scientific codes-often spanning hundreds of thousands of lines-where consistent transformations, strict numerical fidelity, and validation against production tests are non-negotiable. FGPT addresses these challenges by preserving program semantics throughout the entire translation. We verified the framework on representative climate modeling kernels and demonstrated that it produces correct, differentiable Python implementations without requiring manual intervention. By combining rigorous compiler techniques with modern accelerator support, FGPT offers a scalable, trustworthy path for modernizing legacy Fortran code.

Quasi-SVD: Learning a Lie-constrained matrix factorisation for real-time imaging

Singular Value Decomposition (SVD) underlies matrix factorisation tasks across computational imaging, with medical applications increasingly demanding real-time processing. Yet SVD algorithms are inherently sequential, constraining real-time GPU throughput and limit online deployment in clinical pipelines. This study introduces Quasi-SVD, a differentiable, fully parallelized matrix factorization framework for GPUs. Rather than enforcing orthogonality on both factors, it guarantees exact orthogonality for a single Lie-parameterized factor while recovering the remaining components through soft constraints, enabling efficient parallel decomposition without iterative singular-vector optimization. This asymmetric design, provably sufficient for valid factorisation, achieves reconstruction fidelity of SSIM = 0.89-0.94 and accelerates computation by 3-20x relative to cuSOLVER and randomised SVD, enabling throughput above 25 FPS. Performance is evaluated on two medical imaging tasks spanning complementary computational regimes: (1) spatio-temporal background subtraction for ultrasound localisation microscopy, requiring high-dimensional matrix separation, and (2) Mueller matrix polarimetry for neurosurgical tissue characterisation, requiring massive batch processing of small matrices. Across both regimes and multiple imaging instruments, the proposed framework demonstrates robust domain transfer and throughput exceeding 25 FPS at clinical matrix scales, a rate sufficient for live image-guided workflows that classical solvers cannot currently support in these settings. By prioritising downstream reconstruction fidelity over exact spectral recovery, Quasi-SVD makes structured matrix factorisation practical for real-time imaging.

Fast, accurate, and differentiable: a neural-network surrogate for NRSur7dq4 precessing binary black hole waveforms

We present a neural network surrogate model that emulates the NRSur7dq4 gravitational waveform model for precessing binary black hole mergers. The surrogate decomposes the waveform into constituent quantities and trains an independent multilayer perceptron (MLP) for each. We validate the surrogate against NRSur7dq4 on 10,000 waveforms spanning its full parameter space ($1 \leq q \leq 4$, $|χ_{A,B}| \leq 0.8$). For representative total masses between 60 and 300 $M_\odot$, median sky-averaged frequency-domain mismatches range from $8.0 \times 10^{-5}$ to $1.7 \times 10^{-4}$, with 95th percentiles below $10^{-3}$. On an NVIDIA L40S GPU the JAX surrogate evaluates a single waveform in about 1 ms end-to-end, roughly 10 times faster than the LALSimulation C implementation of NRSur7dq4, and sustains about 140 times the LALSimulation throughput at batch size 64, making it well suited for both low-latency parameter-estimation samplers and large-scale waveform generation. The full NRSur7dq4 NN waveform-to-likelihood pipeline is implemented in JAX and is differentiable. This is the first neural-network surrogate of a precessing numerical-relativity waveform model to combine validated NR-faithful accuracy with a fully differentiable, GPU-accelerated inference pipeline, enabling gradient-based inference approaches via automatic differentiation including Fisher information matrices, GPU-accelerated nested sampling, gradient-based MCMC and importance sampling.

Certified Parallel-in-Time Sinkhorn for Dynamic Entropic Optimal Transport

Dynamic applications, including optimal-transport Flow Matching, repeatedly solve related entropic optimal transport problems, yet conventional distributed Sinkhorn processes frames sequentially and synchronizes after every iteration. We present TemporalSinkhorn, a parallel-in-time executor that batches future candidates and their repairs without making output accuracy speculative. A centered, row-sharded certificate accepts only a deterministic safe prefix. The remaining candidates share packed Sinkhorn updates; an online projective forgetting rate places audit milestones, while a posteriori residual checks recover from every depth underestimate. Prediction can therefore change work placement but cannot authorize an inaccurate output. On 4 A100 GPUs, a 60-run, five-seed grid at n = 2048 shows that forgetting-guided milestones reduce wall time by 1.15x-1.47x relative to auditing every packed iteration in five statistically resolved regime cells. Against a sequential soft c-transform warm start, temporal execution is 1.42x-3.55x faster across six synthetic streams, with zero marginal-tolerance violations. On Flow Matching minibatch streams, temporal execution is 3.054x-3.632x faster than sequential carry at n = 2048, with no tolerance violations. A separate fixed-kernel test on an RTX 4060 Laptop GPU gives a 4.315x geometric-mean speedup. These are complementary deployment studies rather than a controlled hardware comparison. End-to-end Flow Matching integration, optimized-solver comparisons, and multi-node validation remain open.

Singular value soft-thresholding via the polar decomposition

Singular value soft-thresholding can be computed via a reduction to the matrix polar decomposition, which allows one to exploit GPU-friendly algorithms for computing the polar decomposition. Empirically, there is a significant speed-up on GPUs compared to the standard approach using the SVD. We leave the investigation of robustness to future work, but note that due to the discontinuous nature of the sign function, the reduction to the polar decomposition is likely only suitable for low-accuracy applications.

A Hardware-oriented Approach for Efficient Bayesian Inference Computation and Deployment

Bayesian inference provides a principled foundation for reasoning under uncertainty, but its computational cost hinders deployment on resource-constrained edge devices. In this paper, we present a hardware-oriented methodology for accelerating discrete Bayesian inference on commercial off-the-shelf embedded GPUs. We identify that the latency of a broad class of variational message-passing algorithms is dominated by tensor contractions. Our approach restructures the memory layout of these operations using two complementary merging strategies that produce compact, regularly-shaped primitives better suited for efficient GPU execution. We then introduce optional sparse array representations and a tensor-clustering scheme to reduce the memory footprint. We instantiate the methodology and produce optimized variants of three message-passing algorithms for Hidden Markov Models (HMMs), namely variational filtering, variational message passing, and marginal message passing. Furthermore, we complement this with a machine-learning-based autotuner that automatically selects the best-performing algorithmic variant for a given generative model specification. Benchmarked on an NVIDIA Jetson Orin AGX across 770 randomly sampled realistic Partially Observable Markov Decision Process (POMDP) configurations, our implementations achieve speedups of up to 5x, with typical gains of 2-2.5x, while producing numerically identical outputs to the baseline implementations.

CayleyR: Solving the TopSpin puzzle via cycle intersection

We present cayleyR, an R package for solving permutation puzzles by detecting cycle intersections in Cayley graphs. The core algorithm performs an iterative bidirectional search: from both the initial and target permutation states, random operation sequences generate cycles in the Cayley graph of the symmetric group Sn; their intersection yields a connecting path. When no direct intersection is found, a distance-guided bridge selection narrows the gap, and the process repeats. The package targets the TopSpin(n,k) puzzle, whose state space is a Cayley graph of Sn generated by a cyclic shift and a prefix reversal. We describe the mathematical framework, the algorithm, and its implementation, which combines a C++ hash-indexed state store with optional Vulkan GPU acceleration. The software is publicly available on CRAN.

VQCSim: When Does Compile-Once Statevector Simulation Beat Generic Quantum Frameworks?

Hybrid quantum-classical machine learning workflows repeatedly evaluate many small parametrized circuits during training and model exploration. In this regime, framework dispatch and orchestration overhead often dominate runtime. Prior simulators accelerate execution but leave open the question of when compile-once specialization is the right choice for static variational circuits. We answer this question with VQCSim, a compile-once, PyTorch-native statevector execution path with native autograd. In a systematic MQT Bench study, VQCSim compiles all tested static circuits and provides 87.7% end-to-end semantic validation. Across a five-GPU evaluation set, VQCSim delivers pooled median speedups of 4.49x for native inference and 26.78x for native training, while retaining a 3.31x advantage under matched finite-difference training. Ablation identifies native autograd as the dominant source of acceleration (27.6x), with compile-once caching and batch vectorization contributing additional gains. The speedup trades higher GPU memory (VQCSim is memory-limited at the high end) for lower runtime. We derive a hardware-aware regime map and release vqcsim-oracle, an open-source backend selector with 91.1%-97.7% top-1 agreement (including cross-GPU transfers), enabling automatic simulator selection in QML design loops.

Scalable Optimal Transport Algorithm for Network Alignment

Network alignment identifies node correspondences across different networks and is a fundamental primitive in many data science applications, including social network analysis, fraud detection, and knowledge graph integration. However, state-of-the-art network alignment methods often achieve high accuracy by repeatedly constructing and updating dense matrices, sacrificing scalability in the process. To address this scalability limitation without compromising alignment accuracy, we present FastAlign, a scalable, sparsity-aware framework for optimal transport-based network alignment. Rather than introducing a new alignment model, FastAlign preserves the original OT formulation and reinterprets its computation as a set of recurring mixed sparse-dense operations. FastAlign combines sparsity-aware graph computation with domain-specific kernel fusion, including a custom SpMM kernel. Our results show that FastAlign achieves alignment quality comparable to state-of-the-art OT-based methods while substantially reducing end-to-end runtime up to 3.89x-9.45x on CPU and 2.24x-32.54x on GPU.

FlashTrie: A GPU-Accelerated Constrained Beam Search for Generative Retrieval

Constrained decoding is essential in generative retrieval, where document identifiers generated directly from a query must exactly match a predefined library of valid IDs. At scale, decoding is often constrained using a trie with beam search but most implementations run on CPU. Limited parallelism then makes trie traversal and candidate validation a serving bottleneck as beam width grows. We present FlashTrie, which addresses this limitation by optimizing constrained beam search on GPUs. It introduces an integer-aware succinct trie layout that uses bit compression to reduce memory footprint while keeping the full index in GPU high-bandwidth memory reducing memory stalls, and a cooperative CUDA kernel that performs beam expansion, validation, and pruning entirely on-device without per-step host orchestration. It further replaces CPU-style irregular lookup and heap maintenance with GPU-aware parallel primitives, improving warp utilization and reducing divergence. Together, these designs significantly reduce decoding latency and increase throughput while preserving retrieval quality. On a library of 800M keywords with beam widths up to 1000, FlashTrie reduces trie-search latency to under 3 ms, achieving up to 24x speedup over a highly optimized multi-threaded CPU baseline. These improvements enable FlashTrie to scale beam sizes by up to 5x in latency-critical applications such as sponsored search. In a large-scale online A/B experiment on a popular commercial search engine, it delivers a statistically significant +0.71% revenue lift, enabling real-time constrained decoding at a scale previously feasible only offline. The FlashTrie code will be publicly released after the review process.

Scaling Weisfeiler-Leman Expressiveness Analysis to Massive Graphs with GPUs

The stable coloring of the Weisfeiler-Leman (1-WL) test is a cornerstone of Graph Neural Networks because it provides an upper bound to the expressive power of message-passing architectures. Unfortunately, computing it presents two fundamental bottlenecks. First, classic algorithms are inherently sequential and cannot exploit modern massively parallel hardware. Second, these are \emph{global} algorithms, i.e., they require availability in memory of the full graph, severely limiting applicability to real-world instances. We leverage a linear-algebraic interpretation of 1-WL stable coloring and introduce two key contributions: (i)~a randomized refinement algorithm with tight probabilistic guarantees and (ii)~a correctness-preserving batching scheme that decomposes the graph into independently processable subgraphs while provably returning a stable coloring of the original graph. This approach maps directly to GPU-efficient primitives. In numerical experiments, our CUDA implementation delivers speedups up to two orders of magnitude over classical CPU-based partition refinement and, for the first time, successfully computes stable colorings on web-scale graphs with over 30 billion edges, where CPU baselines time out or fail.

GPU-Accelerated Inverse Structural Anastylosis from Block Collapse Dynamics

The physical anastylosis of collapsed architectural monuments -- the meticulous reassembly of fallen stone elements into their original structural configuration -- represents one of the most intellectually demanding challenges in conservation science. Traditional approaches depend heavily on expert archaeologist judgement and manual block-by-block correspondence, a process that is both labour-intensive and inherently subjective. Inspired by the combinatorial complexity of this problem as manifested in the game of Jenga, we present Jenga Inverse Predictor , a GPU-accelerated deep learning framework that addresses structural anastylosis as an inverse prediction task. Given an image of a collapsed block assembly, JIP-2 reconstructs the most probable prior tower configuration by: (1) implementing a complete rigid-body physics engine with OBB/SAT collision detection and a Projected Gauss-Seidel (PGS) contact solver accelerated with Numba JIT and CuPy CUDA; (2) applying the analytical force thresholds of Ziglar (CMU, 2006) -- F_app = 3*mu_s*m*g (Y-axis, torque-free) and F_app = 4*mu_s*m*g (X-axis, torque risk) -- over three friction levels (mu_s in {0.25, 0.40, 0.60}) across 450 simulated episodes; (3) training a dual-stream ResNet-18 that injects a friction one-hot vector and jointly predicts block removal count, per-position removal probabilities, centre-of-mass imbalance, and Ziglar torque risk; and (4) generating a smooth 3-D video of the block-by-block reverse reconstruction. We discuss implications for computer-assisted anastylosis at the UNESCO Maya site of Uxmal, Yucatan, and provide a detailed technical description of the full pipeline, architecture, and loss formulation.

TurboMPC: Fast, Scalable, and Differentiable Model Predictive Control on the GPU

Robotics increasingly relies on GPUs for parallel simulation, large-scale learning, and neural-network inference. For model predictive control (MPC) to scale with this paradigm, solvers must run efficiently on this hardware while remaining fast, differentiable, and compatible with expressive MPC formulations used in robotics. We present TurboMPC, a differentiable MPC solver that runs entirely on the GPU and supports state and control inequality constraints, implicit integrators, cross-time-coupled costs, and slack variables. TurboMPC combines sequential quadratic programming (SQP), an alternating direction method of multipliers (ADMM) inner solver, implicit differentiation, and a co-designed JAX-CUDA implementation for efficiency and ease of use. In simulation, we validate TurboMPC on constrained planning, humanoid imitation learning, and reinforcement learning with neural-network cost function tasks, achieving up to $15\times$ and $58\times$ speedups over state-of-the-art CPU and GPU differentiable solvers, respectively. We deploy TurboMPC on a full-scale car for minimum-time racing and find that batched, GPU-accelerated tuning of MPC parameters via Bayesian optimization yields significantly faster driving than a hand-tuned baseline. TurboMPC also scales to planning horizons of over $8000$ knot points while maintaining control of the vehicle. We open-source TurboMPC at: https://github.com/ToyotaResearchInstitute/turbompc

Mojo: A Promising Tool for Scalable Financial AI Efficiency

For thirty years, quantitative finance has paid a costly two-language tax: models researched in Python are rewritten in C++ for production, often introducing numerical discrepancies. GPU-accelerated deep learning exacerbates this problem, as nondeterministic floating-point reductions can produce drift in long backtests, challenging regulatory reproducibility and auditability expectations. This article surveys Mojo, Modular's 2026 Python-like systems language, as a structural response for capital markets engineering. While closing the Python-to-C++ performance gap, Mojo uniquely combines native interoperability with the low-level systems control required to construct bit-exact deterministic kernels. Its MLIR compilation infrastructure further allows a single codebase to target scalar, SIMD, multicore, and GPU execution, reducing the translation bottleneck between research and production. We benchmark four core financial AI workloads: Monte Carlo option pricing, LLM sentiment inference, multi-asset backtesting, and portfolio Value at Risk. On Apple Silicon, Mojo demonstrates 20x to 180x speedups over pure Python on directly measured kernels; larger-scale GPU workload results are projections calibrated from published benchmarks. Alongside transparent performance data, we introduce mojo-deterministic, an open-source library of reproducible reduction kernels, and provide a candid assessment of the problems Mojo does and does not yet solve.

GraphGP: Scalable Gaussian Processes with Vecchia's Approximation

Gaussian processes are a powerful tool for modeling continuous fields, but their naive $\mathcal{O}(N^3)$ computational cost and $\mathcal{O}(N^2)$ memory requirement often limit their practical use. Vecchia's approximation is a sparse precision matrix approximation for stationary, decaying kernels that conditions each point only on its $k$ nearest neighbors. We present GraphGP, a GPU algorithm for Vecchia's approximation that scales to nearly a billion parameters with linear time and memory requirements, handling arbitrary point distributions over a large dynamic range. Our key contributions are (1) a bit-reversed k-d tree ordering that allows efficient neighbor searches while also maximizing batch parallelism, and (2) a differentiable CUDA implementation, which is substantially faster and more memory efficient than our pure JAX baseline. GraphGP provides the building blocks for inference, including forward generation, inverse application, log-determinant, and kernel parameter derivatives.