Solver Iterations

Follow this topic Ask about this topic

Momentum

7 papers in the last four weeks, against 1 the four weeks before. 0.1% of all new papers.

Jul 6Week of Sep 21

Latest papers 46

Sep 29, 2026cs.LG

The Advantages of Fresh Sketching for Ridge Regression

Over the past 25 years, sketching and sampling have become widely used tools for accelerating large-scale regression. In iterative randomized solvers, a basic design choice is whether to

\textit{reuse}

the same sketch or draw

\textit{fresh}

randomness at every step. For (under-constrained) iterative ridge regression with column sampling, whether fresh sketches offer provable advantages has remained open:

\textit{We show that they do.}

Fresh sketching lets us analyze error only along the current residual solution, rather than uniformly over the entire Gram matrix. This directional view yields sharper convergence guarantees for leverage score and ridge leverage score sampling and, more importantly, leads to residual-aware sampling rules. By minimizing the variance of the relevant sketched matrix-vector product, we derive an oracle distribution and practical approximations to the oracle distribution, including a mixture sampling distribution with (somewhat weaker) convergence guarantees. Experiments on synthetic and real data, including ridge probes on Qwen2.5 representations, support our theory, showing substantially faster convergence.

Linkai Ma, Qilin Li, Petros Drineas

Purdue University · University of Wisconsin-Madison

Kernel Ridge Regression Sketches Convergence Solver Iterations

Sep 28, 2026cs.DC

GPUPhysBench: Benchmarking Coding Agents for Correct and Efficient GPU Physics Simulation

Writing fast GPU code for physical simulation is difficult: implementations must preserve numerical accuracy while handling irregular data access, synchronization, and iterative solvers. We introduce GPUPhysBench, a benchmark of 50 tasks testing whether coding agents can meet these demands. Tasks cover fluids, deformable solids, and granular materials, from individual simulation operators to complete simulators. Agents write, compile, test, and optimize GPU code with access to a NVIDIA GPU under fixed time budgets. We report pass rates and runtime performance relative to expert-optimized reference implementations. In a single-attempt evaluation of six frontier model-harness pairs, the two strongest pass all 50 tasks, but even the fastest reaches at least 0.9 the reference speed on only 22% of them, and no submission is more than 5% faster than the reference. The largest gaps arise in collision detection, constraint solving, and iterative solvers. GPUPhysBench brings physical simulation workloads to coding-agent evaluation, testing both the ability to implement numerical methods correctly and the ability to make them run efficiently.

Yuchen Sun, Jinjin He, Sinan Wang +1

Georgia Institute of Technology

Physics Simulation Coding Agents Solver Iterations

Sep 27, 2026cs.LG

JIVE: Jacobian-Informed Volume Expansion for Diverse Generative Sampling

Generative models often suffer from mode collapse and limited sample diversity. While prior works attempt to mitigate this by jointly generating a batch of samples and repelling their trajectories, these heuristics do not explicitly maximize the diversity of the resulting endpoints. We introduce JIVE, a training-free framework that enhances generative diversity by injecting velocity perturbations aligned with the leading right singular subspace of the generator's endpoint Jacobian. By leveraging this local geometric structure, JIVE provably maximizes endpoint diversity while preserving sample quality. To maintain practical efficiency, we compute these perturbation directions via matrix-free iterations rooted in classical numerical linear algebra, requiring only a small computational overhead. Across different benchmarks, JIVE boosts both pixel and feature-level diversity in few-step and one-step generation.

Guangxun Zhang, Brian Cai, Boxuan Zhang +2

New York University · Stony Brook University · Rutgers University

Diversity Jacobian Solver Iterations Collapse

Sep 27, 2026math.OC

lapanda: A Matrix-Free Differentiable Solver for Nonconvex Constrained Optimization Layers

Differentiable optimization brings the structural guarantees of mathematical optimization to network pipelines, allowing them to be trained end-to-end. However, its application remains challenging for nonconvex constrained problems, as existing differentiable solvers often suffer from limited modeling expressiveness due to their reliance on specialized problem structures, while also incurring substantial computation time and memory overhead in both the forward and backward passes. To address these challenges, we propose lapanda, a matrix-free differentiable solver for nonconvex optimization with general constraints. It reformulates the problem to a sequence of augmented Lagrangian subproblems, each handled by a first-order inner solver through a proximal averaged quasi-Newton algorithm with adaptive linesearch, thus enabling efficient forward optimization. We establish local well-posedness of the solution map and convergence of the outer iterations, and further derive a sensitivity alignment between the original problem and the final subproblem in the backward pass, demonstrating that the subproblem sensitivity, which can be computed efficiently in a matrix-free manner, provides a principled approximation to the exact optimizer sensitivity. We evaluate lapanda on nonconvex constrained Rosenbrock benchmarks, imitation learning with several representative constrained optimal control problems, and embedded robotic obstacle-avoidance tasks. Compared with state-of-the-art differentiable solvers, lapanda delivers substantial reductions in computation time and memory footprint while maintaining reliable constraint satisfaction and learning performance.

Yuankun Chen, Zifei Nie, Kangyu Lin +2

School of Artificial Intelligence, Jilin University, China · Graduate School of Informatics, Kyoto University, JPN · Department of Civil and Systems Engineering, Johns Hopkins University, USA

Differentiable Optimization Convex Optimization Solver Iterations Optical Networks

Sep 23, 2026quant-ph

Repairability of Inexact Solvers in Recursive State Estimation with Machine Learning

Recursive state estimation often executes approximate numerical solutions inside a feedback loop, where highly accurate local steps do not guarantee better overall results. For a fixed linear Kalman model, we characterize when a correction within a prescribed subspace and norm budget can meet a local admissibility tolerance, and how the defects actually executed affect the finite-horizon covariance response. Centering each defect on the exact gain for the implemented covariance separates current solve error from inherited gain drift. Expanding the exact residual-drift identity reveals opposing quartic contributions beyond the quadratic response: innovation-covariance inflation enters positively, while local-gain reoptimization enters subtractively. Under matched initialization, an absolute sixth-order remainder bound, uniform over bounded defect sequences at fixed horizon, gives sufficient conditions for quadratic under- or overprediction. Machine learning proposes bounded corrections, while a learner-independent residual certificate and verified fallback govern execution of classical and quantum candidates without changing the reference estimator. In a power-grid tolerance study, learned correction lowers the minimum conjugate-gradient iteration count for deployment without fallback relative to uncorrected solves under the same residual certificate. Gains reconstructed from a variational quantum linear solver and from an annealing-based binary encoding, with small-scale terminal measurements on superconducting hardware and sampling on a quantum annealer, are executed through the same interface. By linking local repairability to nonlinear error propagation, the framework evaluates approximate solvers and learned corrections through independent certification and finite-horizon response, providing a practical basis for studying hybrid quantum--classical computation.

Yanjun Ji, Dennis Willsch, Orkun Şensebat +8

Jülich Supercomputing Centre, Forschungszentrum Jülich, 52425 Jülich, Germany · FH Aachen University of Applied Sciences, 52066 Aachen, Germany · Institute for Automation of Complex Power Systems, RWTH Aachen University, Aachen, Germany +5

State Estimation Solver Iterations Quantum Annealing Variational Quantum Algorithms Covariance Repair

Sep 23, 2026cs.LG

Probabilistic and Geometry Aware Neural Surrogate of Scrape Off Layer Plasma Simulations

Fast surrogates for tokamak boundary-plasma simulation are typically deterministic regressors mapping a global operating point to a flattened vector of cell values. Near the divertor detachment transition the steady state is not reliably single-valued. A point estimate must average over qualitatively different plasma states, and it arrives with no statement of confidence. Moreover, the flattened vector representation discards the geometric structure of the SOLPS-ITER mesh. This work addresses both problems. We unroll the curvilinear mesh into three fixed-size image tensors whose layout preserves cell adjacency and inverts exactly, letting a convolutional network act on the geometry without loss of information. A conditional flow matching model, well suited to highly sensitive systems, is then trained on this representation. The result is an efficient, scalable surrogate that captures multiple plausible outcomes even at sensitive operating points. Along a gas-puff scan, the predictive distribution splits into a hot and a cold mode across an early regime transition. A further check on synthetic data with an injected bifurcation of known size confirms the model recovers both branches rather than their average.

Gabriele Gianuzzo, Stefan Dasbach, Fleur Hendriks +2

Eindhoven University of Technology · DIFFER

Surrogate Models Surrogates Solver Iterations Geometry-Aware

Sep 21, 2026cs.RO

SPARSER: Sparse Variable Projection by Exploiting Separable Structure in Robotic Perception

Robotic perception often requires solving large nonlinear least-squares (NLS) problems. While sparsity has been widely exploited to scale solvers, a complementary and underused structure is \emph{separability}: some variables, such as visual landmarks, appear linearly in the residuals and admit a closed-form solution once the remaining variables, such as poses, are fixed. Variable projection (VarPro) exploits this structure by analytically eliminating the linear variables, yielding a reduced problem with favorable computational properties. However, its use in robotic perception has been limited by gauge symmetries, such as invariance to global translations and rotations, which introduce challenges for standard VarPro methods. We present SPARSER (\textbf{S}parsity \textbf{P}reserving \textbf{A}nalytic \textbf{R}eduction for \textbf{S}eparable \textbf{R}obotic \textbf{P}erception), a VarPro framework for gauge-symmetric problems that jointly exploits separability and sparsity. Our method constructs a \emph{matrix-free Schur complement operator} for efficient evaluation of reduced costs, gradients, and Hessian-vector products, enabling integration with iterative NLS solvers. We characterize the applicable problem class, identify common cases admitting further analytical simplifications, and show that IRLS-based robust costs preserve most of the exploitable structure. Across synthetic and real SLAM, SNL, and SfM benchmarks, SPARSER is on average

5\times

--

7\times

faster than state-of-the-art baselines on CPU and GPU, with gains exceeding

40\times

on individual datasets. On outlier-corrupted multi-robot SLAM data, the robust variant is

2\times

--

16\times

faster than a state-of-the-art GNC solver. We release open-source C++ code and all datasets.

Nikolas R. Sanderson, Andrew Fishberg, Haoyu Han +5

Department of Naval Architecture and Marine Engineering, University of Michigan, Ann Arbor, MI, USA · Institute of Experiential Robotics, Northeastern University, Boston, MA, USA · Aerospace Controls Lab, Massachusetts Institute of Technology, Cambridge, MA, USA +1

Robotic Perception Simultaneous Localization And Mapping Solver Iterations Sparsity Projection Separability

Sep 17, 2026cs.LG

PosteriorBench: From Point Estimates to Posterior Matching in Evaluating Generative Inverse Solvers

Generative models are increasingly used to solve scientific inverse problems, but existing evaluations still focus primarily on whether a method can produce a single plausible reconstruction. This is insufficient for ill-posed problems, where multiple solutions may be consistent with the same sparse or noisy observations. In these settings, a method can achieve strong pointwise accuracy while still failing to capture the true posterior through mode collapse, overconfident uncertainty, or averaging incompatible solutions. We introduce PosteriorBench, a benchmark for evaluating the distributional accuracy of generative inverse solvers. PosteriorBench evaluates four physics-based inverse problems: Darcy flow inversion, Poisson source recovery, carbon capture and storage, and light transport material inference. For each task, we construct high-fidelity reference posteriors using computationally heavy but established procedures such as rejection sampling and Markov chain Monte Carlo, enabling direct assessment of whether solvers recover the full set of solutions rather than the single best sample. We pair these references with a five-metric posterior evaluation suite: posterior-mean error, posterior-standard-deviation error, maximum mean discrepancy, sliced Wasserstein distance, and radially averaged power-spectrum error. These metrics assess pointwise accuracy, marginal uncertainty, distributional alignment, and global frequency fidelity. The benchmark spans sparse sensing, low-resolution observations, nonlinear forward models, varying noise levels, and multimodal priors, with a unified pipeline for distribution matching and uncertainty quantification. Our experiments reveal substantial distribution-matching gaps across current solvers, while showing that neural operators improve resolution robustness, and guidance weights and generation noise are key to posterior-variance calibration.

Jiachen Yao, Zi-Siang Hsu, Xi Deng +5

California Institute of Technology · National Taiwan University · Lawrence Berkeley National Laboratory +3

Bayesian Inverse Problems Posterior Sampling Posterior Solver Iterations Illumination Variations

Sep 3, 2026cs.LG

Time Without Timesteps: Simulating Coupled Dynamical Systems via Self-Consistency

Numerical simulation of dynamical systems is usually organized as a causal march through time: each state is computed from the previous one. We explore a different formulation for coupled systems. For each subsystem type we train a neural surrogate mapping a full driving trajectory and initial condition directly to a full output trajectory; following classical waveform relaxation, coupled systems are assembled by enforcing self-consistency among these trajectories: simulation becomes a fixed-point problem over complete trajectories rather than a stepwise rollout. On coupled van der Pol oscillators and Hodgkin-Huxley neuron networks, sequential depth becomes the number of solver iterations: 4-10 Newton iterations where the reference integrator takes 1500 steps. The gradient likewise loses its time recursion: it becomes a linear system solved by GMRES at memory independent of solver depth. A single scalar measured from the learned operator, the spectral radius of its Jacobian, predicts in advance where the coupled solve will converge; past that boundary, unrolled backpropagation diverges and a Neumann adjoint fails, while the implicit gradient remains correct to 0.04%. We report where the approach succeeds and where surrogate error degrades it.

Liyu Zerihun, Mark Shinyoung Lee

Dynamical Systems Slow Relaxation Modes Solver Iterations Self-Consistency Timestep

Aug 11, 2026cs.LG

DEFT: Data-Efficient Frequency-domain Top-k Sampling via Inverse Discrete Fourier Transform for Spatiotemporal Dynamical Systems Modeling

Modeling spatiotemporal dynamical systems governed by partial differential equations (PDEs) poses two major challenges: it either requires expensive physics-based simulators that entail iterative numerical solving at high computational cost, or it depends on abundant training data, yet purely data-driven models often generalize poorly to downstream dynamic operating conditions. We propose DEFT, a frequency-domain data sampling method that identifies the dominant Fourier modes of a physical system and systematically varies the corresponding amplitudes and phases to generate physically consistent training data via the inverse discrete Fourier transform. In addition, we derive a generalization bound of this method. We note that it also provides a theoretically principled criterion for selecting

K

. We evaluate the proposed method through three sets of experiments, each targeting a distinct aspect of its utility. First, we validate the framework on canonical PDEs solving demonstrating that it outperforms traditional methods when the system is dominated by a few prominent frequency components. Second, we employ DEFT as a data-value filter on the diffusion--sorption and Burgers equations of PDEBench, showing that it reduces data requirements by

40\%

while sacrificing less than

2\%

in predictive accuracy. Third, to evaluate DEFT for more challenging and practically relevant problems, we validate it in the battery degradation PDE system, achieving consistently high predictive accuracy across various test datasets with

R^2

values exceeding

0.99

. Moreover, the learned frequency-domain features transfer to other battery chemistries with only

20\%

of the fine-tuning data. These results demonstrate that DEFT is an effective data-sampling method for efficient operator learning.

Hengbo Xiao, Jiale Liu, Jiahao Song +1

Peking University

Frequency Domain Partial Differential Equations Fourier Spatiotemporal Solver Iterations Lithium-Ion Batteries

Aug 2, 2026math.OC

Rake-Compress Riccati Recursions for Parallel Scenario-Tree Model Predictive Control

Scenario-tree model predictive control (MPC) represents future information by a rooted tree and optimizes a nonanticipative policy over that tree. Numerical methods for solving the resulting nonlinear program typically compute their search directions through a sequence of branched linear-quadratic regulator (LQR) subproblems. The standard tree Riccati recursion requires linear work but has a dependency chain proportional to tree height. We present an algebraically exact parallel solver based on rake-compress tree contraction. After independent local control condensation, its two operations act on node and edge data that represent conditional quadratic functions. A rake eliminates a leaf and its parent edge, adding their reduced contribution to the parent-node data. A compress eliminates a unary node and replaces its two adjacent edges by one edge, using the same conditional-value composition as parallel Riccati methods on a chain. Together they contract an arbitrary rooted tree to its root; reversing the contraction recovers every Riccati coefficient, state, control, and multiplier. Given a reusable topology plan, a solve with

N

nodes and fixed state and control dimensions has

O(N)

arithmetic work and storage and

O(\log N)

span, independently of tree height, balance, and maximum out-degree. The formulation allows positive-semidefinite dual regularization, including the unregularized case, and an exact linear-size lifting covers the standard scenario-MPC convention of one control per information node. We prove the contraction identities and equivalence to the Karush-Kuhn-Tucker (KKT) system. Three MIT-licensed JAX packages implement the bidirectional contraction, the dual-regularized LQR solver, and a user-facing primal-dual interior-point solver for tree-structured optimal control.

João Sousa-Pinto

independent researcher

Model Predictive Control Linear Quadratic Regulator Optimal Control Solver Iterations Markov Decision Processes Data Compression Methods

Jul 30, 2026math.NA

Graph Neural Multilevel Preconditioners for Iterative Solvers

Solving large, sparse linear systems is a core task in scientific computing, and efficient iterative solvers rely critically on effective and robust preconditioning. While classical methods such as algebraic multigrid (AMG) are highly scalable, their robustness can degrade on indefinite or nonsymmetric systems where heuristics originally developed for elliptic PDEs are less reliable. Recently, Graph Neural Networks (GNNs) have emerged as data-driven preconditioners; yet, the practical impact of imposing an AMG-style hierarchy remains underexplored for general sparse matrices. In this work, we propose a Graph Neural Multilevel Preconditioner (GMP) that adopts an AMG hierarchy as a structural prior and learns smoothing, restriction, and interpolation operators in a unified framework. Our method targets general sparse systems and is instantiated as a drop-in preconditioner for standard Krylov solvers. On a benchmark of over 800 sparse matrices, we compare against classical AMG, single-level ILUT, and state-of-the-art GNN preconditioners, and characterize the regimes where multilevel graph neural preconditioning improves convergence or, conversely, introduces overhead relative to strong single-level baselines. These results highlight both the promise and the limitations of enforcing AMG-style multilevel structure in learned preconditioners for large-scale scientific simulations.

Zechen Zhang, Rui Peng Li, Yousef Saad

University of Minnesota Minneapolis, MN, USA · CASA, Lawrence Livermore National Laboratory Livermore, CA, USA

Algebraic Multigrid Spectral Preconditioning Solver Iterations Graph Neural Networks Preconditioner Solvers

Jul 27, 2026cs.RO

Amortising Trajectory Optimisation for Residual MPC via Implicit Contact Differentiation

Differentiable simulation can accelerate contact-rich trajectory optimisation by exposing local sensitivities of task outcomes to controls. Existing approaches either use finite differences, which are expensive and step-size sensitive; differentiate iterative contact solvers by unrolling automatic differentiation (AD), which stores a growing computation trace; or require intricate, solver-specific KKT sensitivity derivations. We introduce an AD-assisted implicit derivative for regularised smooth contacts and apply it to Mujoco MJX, based on the Implicit Function Theorem (IFT). The method differentiates the stationarity residual at the tolerance-converged solution, avoiding both solver unrolling and hand-assembled KKT systems. IFT keeps compiled temporary memory nearly constant with solver effort, changing by less than 4

\%

from one to ten iterations versus 10.6

\times

growth for unrolled AD. IFT memory grows slower with active contacts and model dimension, using 20

\times

less memory at 256 contacts and 6

\times

less at 16 contacts and 96 DoF. We further introduce optimiser distillation for residual MPC, amortising batched full-horizon iLQR into a policy that guides short-horizon residual iLQR. Across Finger, Franka, and Unitree, this raises six-step success by 28-98 percentage points over standard iLQR.

Daniel Layeghi, Thomas Corbères, Calum Arnott +4

School of Informatics, University of Edinburgh, Scotland.

Robust Trajectory Implicit Differentiation Model Predictive Control Rigid-Body Dynamics Alternating Direction Method Of Multiplier Contact+2

Jul 24, 2026cs.LG

Complexity Bounds and Approaches to Learning Projected Gradient Descent Solver Iterates

Data scarcity poses a fundamental challenge in training generative models to produce initial guesses for parametric optimization problems that are otherwise numerically expensive to solve. We therefore study a

k

-neighborhood data collection strategy that augments datasets of converged solutions with intermediate solver iterates, increasing the amount of training data without additional solver runs. To understand the benefits of this approach, we derive a generalization bound based on Rademacher complexity that reveals the role of the

k

-neighborhoods and related parameters. To achieve this result, we focus on one-sided box-constrained quadratic programs solved by projected gradient descent. We illustrate the behavior of this solver on two examples. The approach proposed in this paper enables a more capable DDDAS paradigm by improving the efficiency of the data-model-optimization loop. We finish by discussing two views of learning solver-iterate data and connect our analysis with GLENS, a new data-efficient global search method.

Anjian Li, Ryne Beeson

Department of Electrical and Computer Engineering Princeton University, Princeton, NJ 08544, USA · Department of Mechanical and Aerospace Engineering Princeton University, Princeton, NJ 08544, USA

Gradient Descent Solver Iterations Quadratic Programs Generative Models

Jul 23, 2026physics.comp-ph

Cycle-Consistent and Uncertainty-Aware Neural Surrogates for Tokamak Edge Plasmas

The boundary and divertor plasma govern how a tokamak exhausts power and particles, setting heat fluxes, target conditions, and the onset of detachment. Predicting these quantities is essential for operating current and future devices, but edge simulations that resolve them are too slow for parameter scans, optimization, or real-time control. Machine-learning surrogates offer a fast alternative, yet most are forward-only: they cannot recover input parameters from observations or assess the reliability of their predictions. We introduce a cycle-consistent neural surrogate for edge plasmas, combining a conditional U-Net forward model with an optimization-based inverse method built on the frozen forward network. The forward model maps five control parameters to two-dimensional plasma-state fields on the SOLPS-ITER mesh; the inverse method enforces consistency between forward and inverse predictions, a self-supervised quality check needing no ground-truth labels at inference. An ensemble of multilayer perceptrons also predicts electron temperature and density profiles at the outboard midplane and divertor targets, with uncertainty estimates that flag where more simulations are needed. The forward model achieves normalized root-mean-square errors below 2.6% and Pearson correlations above 0.95 for all fields. Cycle-consistency regularization raises the average cyclical

R^2

from 0.59 to 0.99 without degrading forward accuracy and enables recovery of the core fueling rate; all five control parameters are recovered with Pearson

r\ge0.97

. A

k

-d tree warm start yields a database completion rate above 95%, versus roughly 30% outright failures when cold-started. With about

4\times10^6

parameters, the model produces full 2D predictions in milliseconds, five to six orders of magnitude faster than SOLPS-ITER, enabling real-time control, parameter scans, uncertainty analysis, and digital twins.

Abdourahmane Diaw, Sebastian De Pascuale, Jae-Sun Park +3

Fusion Energy Division, Oak Ridge National Laboratory, Oak Ridge, TN 37932, USA · DIFFER - Dutch Institute for Fundamental Energy Research, De Zaale 20, 5612 AJ Eindhoven, the Netherlands

Surrogate Models Nuclear Fusion Solver Iterations

Jul 8, 2026cs.LG

Eigenbasis-Independent Learnable Spectral Positional Encodings for Directed Graphs via Hermitian Block Krylov Subspaces

Spectral positional encodings (PEs) for \emph{directed} graphs face two obstacles: magnetic Laplacians require an

O(n^3)

Hermitian eigendecomposition per potential, and their complex eigenvectors are defined only up to unitary gauge, which prior work handles with basis-invariant architectures. We propose learnable spectral PEs of the form

h_θ(A_q)\,R

, where

A_q

is a normalized magnetic operator,

h_θ

a learnable scalar spectral response, and

R

a block of random probes. Because the PE is a \emph{matrix function} of the operator, it is gauge-invariant by construction. We compute it in a Hermitian block Krylov subspace from sparse matrix--vector products only, prove that

k = O(\log(1/\varepsilon))

block steps suffice uniformly over heat--resolvent response families, and give a covering-number argument for why low-dimensional structured families generalize where free per-eigenvalue weights overfit. On a directed SBM whose symmetrization is uninformative by construction, direction-blind PEs stay at chance while magnetic Krylov PEs converge to the exact-eigendecomposition oracle as the depth grows. The same probes yield gauge-invariant pairwise features with

1/\sqrt{s}

Monte-Carlo error, and the undirected

q{=}0

case improves heterophilous benchmarks over no-PE and polynomial baselines.

Jiaqing Xie, Yuxin Wang

Fudan University

Graph Edge Sparsification Graph Laplacians Spectral Decomposition Positional Encoding Solver Iterations Symmetry

Jul 6, 2026cs.MA

Collective Intelligence with Foundation Models

As foundation models grow in scale and diversity, coordinating multiple models into cooperative reasoning systems offers a path toward safer, more reliable AI. This chapter presents a multi-agent framework where solver models generate independent drafts, each undergoes structured critique and revision by a critic agent, and an aggregator agent synthesizes a final consensus solution. A scoring module provides semantic, numerical, and procedural evaluation across all agents. Through ablation studies on a benchmark spanning calculus, physics, chemistry, biology, economics, optimization, statistics, and mathematics, we isolate the contributions of framework architecture versus model diversity. We compare four configurations: (1) Individual Baseline, (2) Homogeneous Framework using one shared model, (3) Redundant Homogeneous Solvers using multiple instances of the same model, and (4) Heterogeneous Framework with diverse specialized models. Results show that while framework structure and redundant sampling yield modest gains, model heterogeneity is the critical factor driving substantial performance improvements. The heterogeneous configuration achieves superior step-wise accuracy (0.64 vs. 0.54 for individual models; 2.3x improvement over homogeneous configurations) with reduced variance across categories and difficulty levels. Step-wise reasoning quality (correctness of intermediate steps, not just final answers) improves dramatically only with model diversity, showing that heterogeneous agents provide complementary error detection and reasoning refinement essential for explainability and auditability. We discuss architectural principles, evaluation methodology, and implications for Global Applied AI, showing how heterogeneous multi-agent coordination supports transparent, auditable, high-confidence decision-making across scientific and industrial domains.

J. de Curtò, I. de Zarzà

Department of Computer Applications in Science & Engineering, BARCELONA Supercomputing Center, 08034 Barcelona, Spain · Escuela Técnica Superior de Ingeniería (ICAI), Universidad Pontificia Comillas, 28015 Madrid, Spain · Human Centered AI, Data & Software, LUXEMBOURG Institute of Science and Technology, L-4362 Esch-sur-Alzette, Luxembourg

Multi-Agent Reasoning Heterogeneous Modeling Tools Foundation Model Intelligence Diversity Solver Iterations

Jul 3, 2026cs.LG

A Near-Linear-Time Solver for Graph $p$ -Laplacian Semi-Supervised Learning via Continuation in $p$

Graph-based semi-supervised learning (SSL) propagates a few labels over a similarity graph by minimizing a Dirichlet-type energy. The standard quadratic (

p=2

) energy reduces to a single graph-Laplacian solve, but it degenerates exactly where SSL is most useful when labels are scarce: gathering more unlabeled data drives the

p=2

estimate to a near-constant function whenever

d\ge2

(Nadler-Srebro-Zhou). Well-posedness requires the nonlinear

p

-Laplacian energy with

p>d

. Existing solvers reduce this to a sequence of weighted Laplacian solves, but their reference implementations use a direct sparse factorization or ichol-preconditioned CG instead. Plugging a near-linear Laplacian solver is not straightforward: at large

p

the conductance weights degenerate near flat-gradient edges, making the system nearly singular and causing stagnation without a damped outer iteration. We close this gap. Recasting

p

-Laplacian SSL as a source-form nonlinear Laplacian flow

Bρ_p(B^\top x)=b

and solving by damped chord-Newton continuation in

p

, every linearized system stays well-conditioned and can be delegated to a near-linear Laplacian engine. On size-scaled graph families the wall-clock is empirically

m^{0.96}

-

m^{1.02}

per family (approximate Cholesky default), and a pooled fit across 228 SuiteSparse graphs gives

m^{1.19}

vs.\

m^{1.45}

for direct factorization; the solver handles a

6.8\times10^7

-edge social network in minutes. Memory is the binding constraint: Cholesky fill reaches

10

-

280\times

the graph nonzeros vs.\ our

O(m)

hierarchy. Against the released FCL solver we are

1.5

-

14\times

faster at matched accuracy. On MNIST

10

-NN,

p=3

scores

64\%

at one label per class vs.\

36\%

for

p=2

. Code: https://github.com/orenlivne/np.

Oren E. Livne

Pine Birch Analytics, Churchville, PA

Graph Laplacians Solver Iterations Semi-Supervised Learning Convex Relaxation Continuation

Jun 30, 2026cs.LG

SNAP-FM: Sparse Nonlinear Accelerated Projection for Physics-Constrained Generative Modeling

Generative models have emerged as scalable surrogates for physical simulation, yet they offer no guarantee that their outputs respect the conservation laws, boundary conditions, and nonlinear invariants that govern the underlying physics. Constrained sampling closes this gap, enforcing such constraints exactly at inference time without retraining, but at a computational cost: projection, correction and trajectory-optimization steps are repeated during sampling, with these steps becoming expensive for nonlinear constraints. Standard ML frameworks exacerbate this: their dense tensor algebra and limited sparse solver composability obscure the structure that physical constraints naturally induce, making efficient batched nonlinear optimization difficult to realize in practice. We address this bottleneck by exploiting the structure that sample-wise batching and local PDE couplings induce in the projection subproblems -- namely, block-sparse Jacobian and KKT systems -- exposing this structure using ExaModels.jl and solving the resulting sparse nonlinear programs with MadNLP.jl and GPU sparse factorization. Applied to Physics-Constrained Flow Matching (PCFM), on PDE benchmarks with linear, nonlinear, one-dimensional, and two-dimensional constraints, this approach accelerates nonlinear constraint projection while maintaining constraint satisfaction. These results show that sparse GPU nonlinear optimization is a practical foundation for constrained generative sampling in scientific machine learning.

Alaina Kolli, Theodoros Xenakis, Utkarsh Utkarsh +4

Massachusetts Institute of Technology, Cambridge, MA, USA · Norwegian University of Science and Technology, Trondheim, Norway

Inequality Constraints Solver Iterations Generative Models Conditional Flow Matching Physics Simulation Solvers+3

Jun 29, 2026stat.ML

Extrapolating from Regularised Solutions for Solving Ill-Conditioned Linear Systems in Machine Learning

Rapid prototyping of algorithms is a critical step in modern machine learning. Most algorithms exploit linear algebra, creating a need for lightweight numerical routines which -- while potentially sub-optimal for the task at hand -- can be rapidly implemented. For the numerical solution of ill-conditioned linear systems of equations, the standard solution for prototyping is Tikhonov-regularised inversion using a nugget. However, selection of the size of nugget is often difficult, and the use of data-adaptive procedures precludes automatic differentiation, introducing instabilities into end-to-end training. Further, while data-adaptive procedures perform multiple linear solves to select the size of nugget, only the result of one such solve is returned, which we argue is wasteful. This paper aims to circumvent the above difficulties, presenting autonugget; a Python package for automatic and stable numerical solution of linear systems suitable for rapid prototyping, and fully compatible with automatic differentiation using JAX. autonugget combines multiple linear solves using Richardson extrapolation to determine the solution of the ill-conditioned system, improving in accuracy over approximations based on a single nugget.

Disha Hegde, Jon Cockayne, Chris. J. Oates

University of Southampton, UK · Newcastle University, UK · The Alan Turing Institute, UK

Solver Iterations Regularization Extrapolation Full-Precision Performance

Jun 20, 2026math.NA

Spectrally Safe Neural Operator Warm-Starts for Large-Scale Newton Solvers

Neural operators are increasingly used to warm-start Newton solvers for nonlinear PDEs, on the premise that a low test error places the initial guess inside the basin of attraction. We show that this premise is unreliable. An operator trained to the relative

L^2

error

O(10^{-3})

can still produce an initial state in which the discrete Jacobian is indefinite, because the mean-squared training controls error on average while leaving localized pointwise violations of the underlying physics. For a nearly incompressible hyperelasticity problem, we trace this to the predicted volume change: the operator disperses

\mathrm{det} F

well away from one, and the resulting Jacobian acquires negative eigenvalues even when the predicted field is visually indistinguishable from the reference. At a small scale, this is a nuisance; at a multi-million degree-of-freedom scale, it is disqualifying, since the conjugate gradient and other Krylov solvers needed for memory-feasible Newton steps assume a definite spectrum. We then show that a short, label-free fine-tuning phase -- penalizing the operator against the discrete energy, with no additional solution data -- shifts the Jacobian spectrum back to positive definite. Combined with an inexact outer loop, this gives a warm-started Newton method that converges across the full loading range where the unregularized operator fails, reaching up to 5.4

\times

wall-clock speedup over incremental continuation on a 3D problem with 6.4 million degrees of freedom.

Jaemin Oh, Youngkyu Lee, Jerome Darbon +1

Division of Applied Mathematics, Brown University, Providence, RI 02912, USA

Solver Iterations Neural Operators Jacobian Partial Differential Equations

Jun 19, 2026cs.LG

Memory Is No Longer a Bottleneck: Memory-Efficient Graph Filtering for Scalable Collaborative Filtering

Graph convolutional networks (GCNs) have demonstrated significant success in capturing complex user-item relationships for collaborative filtering (CF). However, due to their reliance on extensive model training, training-free graph filtering (GF)-based CF methods have emerged as a promising alternative, offering computational efficiency by smoothing graph signals via matrix operations. In particular, polynomial GF-based approaches demonstrate improved accuracy through their ability to design more expressive and flexible filtering functions. Despite these advantages, existing GF methods suffer from a critical memory bottleneck: they necessitate storing the full item similarity graph, incurring prohibitive memory costs for large-scale datasets, which limits their practical applicability. To tackle this challenge, we propose Mem-GF (Memory-efficient GF), a new GF-based CF method that departs from conventional designs by principally leveraging the structure of Krylov subspaces as a core mechanism for approximating polynomial graph filters without explicitly storing the item similarity graph. We theoretically analyze the minimum Krylov subspace size that guarantees lossless approximation. Through extensive experiments, we demonstrate that Mem-GF achieves up to 5.74

\times

lower memory usage and 4.38

\times

speedup in runtime, while consistently exceeding the recommendation accuracy of state-of-the-art GF and GCN-based methods. Mem-GF robustly scales to datasets with tens of millions of interactions, establishing itself as a practically viable and theoretically grounded solution for efficient CF.

Jin-Duk Park, Won-Yong Shin

School of Mathematics and Computing (Computational Science and Engineering), Yonsei University, Seoul 03722, Republic of Korea

Collaborative Filtering Solver Iterations Similarity and Metrics Bottlenecks

Jun 18, 2026math.NA

A fast direct solver based neural network for solving PDEs

The matrices arising from large scale

N

-body problems can be efficiently represented using hierarchical matrices, whose key idea is that the admissible off-diagonal sub-matrices can be well approximated by low-rank matrices across a hierarchy of matrix partitions. HODLR (Hierarchical Off-Diagonal Low-Rank) matrices are a subclass of hierarchical matrices in which all off-diagonal submatrices at every level of a recursive binary partition are low-rank. In this article, we present a neural network that learns the inverse operation of HODLR matrices based on the fast direct solver for HODLR matrices developed by Ambikasaran and Darve (2013). We further extend the architecture to learn nonlinear solution operators associated with PDEs by replacing some of the linear layers with deep sub-networks. We demonstrate the performance of the proposed architecture by performing a comprehensive set of experiments that include (i) solving a linear problem such as the Fredholm integral equation of the second kind, (ii) solving PDEs such as the nonlinear Schrödinger equation, Burgers' equation, and the steady-state Darcy's flow equation, (iii) generalization study across varying parameter values, (iv) comparing the inference time of the proposed network with the run time of a classical numerical solver, and (v) comparing the proposed network with some of the existing neural operator learning networks.

Jashwanth Reddy Kadaru, Vaishnavi Gujjula

Department of Computer Science & Engineering, International Institute of Information Technology Bangalore (IIIT-B), India · Department of Data Science and Artificial Intelligence, International Institute of Information Technology Bangalore (IIIT-B), India

Neural Partial Differential Equation Solvers Partial Differential Equations Classical Numerical Solver Neural Operators Solver Iterations Neural Network+2

Jun 6, 2026math.NA

Stable and Scalable Probabilistic Numerical Solvers for Stiff and High-Dimensional ODEs

Filtering-based probabilistic numerical solvers for ordinary differential equations (ODEs) have been established as a flexible and efficient simulation framework with built-in numerical uncertainty quantification. However, problems that are both stiff and high-dimensional remain a challenge, as current methods are either stable and have cubic cost in the ODE dimension, or scale linearly at the expense of stability. In this paper, we close this gap and develop probabilistic ODE solvers that are both stable and scalable. We propose two complementary strategies. First, we develop a matrix-free update step that uses Jacobian-vector products, iterative linear solvers, and stochastic covariance estimation to enable linear scaling, all while retaining stability. Second, we propose iterative re-linearization to further improve stability without sacrificing scalability, turning probabilistic ODE solvers into fully implicit methods. We evaluate the proposed approaches on a range of stiff and high-dimensional problems and demonstrate improved stability and scalability over established probabilistic solvers.

Nathanael Bosch

EPFL, Lausanne, Switzerland

Solver Iterations Ordinary Differential Equations Uncertainty Quantification

May 29, 2026cs.LG

GLENS: Global Search via Learning from Solver Iterates with Diffusion Models

We consider the problem of generating a large collection of initial guesses for local minima of multimodal non-convex continuous optimization problems. The goal is for these initial guesses to be high-quality (i.e., a numerical solver converges quickly) and diverse (i.e., represent many different local minima). Identifying multiple locally optimal solutions enables flexible downstream decision-making, but typically requires expensive global search. Existing data-driven methods predict initial guesses using only the final converged optima from offline solver runs, which discards information about the local neighborhoods of solutions and limits the available training data. We propose GLENS (Global Search via Learning from Solver Iterates), a data-efficient global search method that leverages intermediate solver iterates as free data augmentation. GLENS consists of two components: a neighborhood structure model that uses diffusion models to learn the local geometry around optima conditioned on problem parameters, and a solver behavior model that learns refinement directions to further guide samples towards nearby optima during diffusion sampling. Experiments on modified non-convex benchmark problems and a two-robot obstacle-avoidance navigation problem show that GLENS generates high-quality initial guesses while preserving the multimodal distribution of diverse local optima. The resulting initial guesses lead to faster solver convergence across different problem settings and solvers. We also analyze how key hyperparameter choices affect the performance.

Anjian Li, Bartolomeo Stellato, Ryne Beeson

Department of Electrical and Computer Engineering, Princeton University, Princeton, 08544, NJ, USA. · Department of Operations Research and Financial Engineering, Princeton University, Princeton, 08544, NJ, USA. · Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, 08544, NJ, USA.

Large Neighborhood Search Convex Optimization Solver Iterations Diffusion Models Obstacle Avoidance

May 27, 2026cs.LG

LoRe: Adaptive Interaction-Evaluation Routing with Per-Step Interaction Budgets for Iterative Graph Solvers

Diffusion-based neural solvers for combinatorial optimization repeatedly re-evaluate dense edge/factor interactions, making inference expensive in wall-clock time and often memory-bound at scale. Inspired by the computational methodologies of many-body physics, we introduce LoRe, a training-free, inference-time drop-in wrapper that enforces per-step interaction-evaluation budgeting: at each iteration, it evaluates only a fixed fraction of interactions by dynamically routing computation to high-conflict or high-uncertainty interactions, instead of using a fixed sparsification (e.g., static kNN graphs or static masks). Under fully inclusive end-to-end wall-clock accounting, LoRe substantially improves scalability on the Maximum Independent Set (MIS) problem, extending feasible inference more than

3\times

beyond the baseline's out-of-memory limit, delivering a

\sim 8\times

speedup and a

\sim 12\times

peak-memory reduction, with solution quality preserved in this regime. Demonstrating cross-task generality on the large-scale Traveling Salesperson Problem (TSP) and zero-shot robustness to topology shifts, LoRe achieves a

\sim 15\times

speedup at

n=1000

with a

44\times

memory reduction and competitive tour quality.

Jintao Li, Yong-Yi Wang, Zheng-An Wang +1

Beijing Key Laboratory of Fault-Tolerant Quantum Computing, Beijing Academy of Quantum Information Sciences, Beijing, China · Beijing National Laboratory for Condensed Matter Physics, Institute of Physics, CAS, Beijing, China · Beijing Key Laboratory of Advanced Quantum Technology, Beijing, China +1

Neural Solvers Solver Iterations Inference-Time Token Budget Allocation Peak Memory Many-Body Physics

May 26, 2026cs.DC

Accelerating Birkhoff Projection for Manifold-Constrained Hyper-Connections

Manifold-constrained hyper-connections (mHCs) have recently been proposed as a principled extension of hyper-connections, where the residual mixing matrices are constrained to be doubly stochastic via projection onto the Birkhoff polytope. In practical mHC implementations, this constraint is enforced by Sinkhorn-Knopp iterations, and the backward pass relies on unrolling the iterative solver. This design introduces substantial computation and memory overhead, and may also yield inaccurate projections when the algorithm converges slowly on challenging inputs, undermining the intended norm-control and stability guarantees of mHCs. In this work, we focus on the practically important 4x4 Birkhoff projection setting and develop an end-to-end acceleration framework. By leveraging the dual formulation, we reduce the problem to a three-dimensional unconstrained convex problem and solve it with Newton's method, achieving fast convergence and high accuracy. For the backward pass, we replace the unrolled differentiation with implicit differentiation, yielding exact gradients without storing intermediate states. To exploit massive parallelism, we design a warp-level CUDA kernel that uses only register-level primitives, avoiding global and shared memory I/O. Extensive experiments against representative open-source baselines demonstrate that the proposed solver yields substantially more reliable doubly stochastic projections -- especially when the input magnitude is large -- and achieves significant end-to-end speedups (including the backward pass), reaching over 20x acceleration at large batch sizes while maintaining orders of magnitude smaller marginal errors.

Chenrui Wang, Yixuan Qiu

RenminSchoolUniversityof Statisticsof China · School of Statistics and Data Science & Institute of Big Data Research · Shanghai University of Finance and Economics

Solver Iterations Projection Matrix Multiplication Convergence Cuda

May 18, 2026cs.LG

Randomized Advantage Transformation (RAT): Computing Natural Policy Gradients via Direct Backpropagation

Natural policy gradients improve optimization by accounting for the geometry of distribution space, but their practical use is limited by the cost of estimating and inverting the Fisher matrix. We present Randomized Advantage Transformation (RAT), a method for estimating Tikhonov-regularized natural policy gradients via direct backpropagation. By applying the Woodbury formula, we reformulate the regularized natural policy gradients as vanilla policy gradients with a transformed advantage. RAT computes this transformation efficiently via randomized block Kaczmarz iterations on on-policy mini-batches, avoiding explicit Fisher construction, conjugate-gradient solvers, and architecture-specific approximations. We provide convergence guarantees for RAT and demonstrate empirically that it matches or exceeds established natural-gradient methods across continuous and visual control benchmarks, while remaining simple to implement and compatible with various architectures.

Mingfei Sun

The University of Manchester, United Kingdom.

Policy Gradient Backpropagation Solver Iterations

May 18, 2026cs.LG

RL4RLA: Teaching ML to Discover Randomized Linear Algebra Algorithms Through Curriculum Design and Graph-Based Search

Randomized linear algebra (RLA) algorithms are a modern class of numerical linear algebra techniques that play an essential role in scientific computing and machine learning, with broad and growing adoption. However, their discovery remains mostly a manual process that requires deep expert knowledge and inspiration. While Reinforcement Learning (RL) offers a pathway to automation, standard approaches struggle with sparse reward landscapes and vast search spaces inherent to high-performing RLA algorithms. In this paper, we present RL4RLA, a general RL framework that automates the discovery of interpretable, symbolic RLA algorithms. Unlike black-box approaches, our method builds explicit algorithms from basic linear algebra primitives, ensuring verifiable and implementable representations. To enable efficient discovery, we introduce: (1) a numerical curriculum that progressively increments problem difficulty to encode inductive bias specific to the RLA domain; (2) Monte Carlo Graph Search, which optimizes exploration by identifying and merging equivalent partial algorithms. We demonstrate that RL4RLA rediscovers state-of-the-art methods, including sketch-and-precondition solvers, Randomized Kaczmarz, and Newton Sketch, and can be targeted to produce algorithms optimized for specific trade-offs between accuracy, speed, and stability. Code is available at https://github.com/Tim-Xiong/RL4RLA.

Jinglong Xiong, Xiaotian Liu, Ruoxin Wang +4

Pratt School of Engineering, Duke University, Durham, NC, USA · Department of Computer Science, Dartmouth College, Hanover, NH, USA.

Solver Iterations Low-Rank Structure Search Algorithms Offline Reinforcement Learning Riemannian Optimization Curriculum

May 18, 2026cs.LG

Lightweight Gaussian Process Inference in C++ on Metal and CUDA

Gaussian process (GP) inference in Python is dominated by libraries such as GPyTorch and GPflow, which are built on deep-learning frameworks and inherit their dispatch overhead and dependency footprint. We present LightGP, a dependency-free C++17 library for GP regression with Python bindings, supporting Apple Metal and NVIDIA CUDA backends alongside tuned CPU paths via Apple Accelerate and OpenBLAS. LightGP provides four inference paths -- exact Cholesky, matrix-free conjugate gradients, sparse variational free energy, and structured kernel interpolation with FFT -- covering problems from

N{=}100

to

N{=}500{,}000

. On an Apple M4, LightGP CPU is 2.6--8.7

\times

faster than GPyTorch CPU for exact GP and

{\sim}1.5\times

faster for sparse GP at every scale tested. On an NVIDIA RTX~3060, LightGP CUDA is 2.3--6.7

\times

faster than GPyTorch CUDA for exact GP up to

N{=}2{,}048

, with GPyTorch closing the gap at

N{=}4{,}096

. A fused matrix-free kernel-vector product on Metal achieves 32

\times

over the explicit path at

N{=}20{,}000

with

O(N)

memory, and an FFT-accelerated SKI matvec via Accelerate vDSP runs in sub-millisecond time at

N{=}200{,}000

. LightGP compiles as a single static library with zero external dependencies and is installable via \texttt{pip install lightgp

Yu-Hsueh Fang

Department of Information Management, National Taiwan University · 2H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology

Gaussian Process Cpu-Gpu Hybrid Designs Cuda Pytorch Solver Iterations

May 17, 2026cs.LG

Learning Fill-in Reduction Ordering via Graph Policy Optimization for Sparse Matrices

Matrix reordering in large sparse solvers seeks a permutation that minimizes factorization fill-in to reduce memory and computation. Because the minimum fill-in ordering problem is NP-complete and fill-in is implicit in the sparsity pattern, graph-theoretic heuristics are used. Existing reinforcement learning methods either ignore sparsity patterns--missing the global fill-in--or lack local exact fill-in feedback. We propose a graph policy optimization method, modeling fill-ins from global and local views: both the policy and value networks use a multi-hop graph neural backbone to embed global fill-in; the policy further interacts with symbolic factorization over graphs to extract local, step-level fill-ins, and the resulting feedback is aligned with the value network via an adaptive saturation function to improve convergence. On the SuiteSparse Matrix Collection, our method achieves mean reductions of 29.3 in fill-ins and 31.3 in peak memory usage over state-of-the-art baselines.

Ziwei Li, Shuzi Niu, Huiyuan Li +2

Institute of Software, Chinese Academy of Sciences, Beijing, China

Tensor Completion Solver Iterations Convex Relaxation Offline Reinforcement Learning Feedback

May 13, 2026cs.LG

Local Inverse Geometry Can Be Amortized

Nonlinear inverse problems often trade inexpensive but fragile first-order updates against curvature-aware methods such as Gauss-Newton and Levenberg-Marquardt, which obtain stronger directions by repeatedly solving Jacobian-based linearized systems. We propose a learned alternative: amortize local inverse geometry into a reusable reverse operator. Our framework learns a bidirectional surrogate, Deceptron, and deploys it through D-IPG (Deceptron Inverse-Preconditioned Gradient), an iterative solver that pulls residual-corrected measurement-space proposals back to latent space. The key mechanism is a Jacobian Composition Penalty (JCP), which trains the reverse Jacobian to act as a local left inverse of the forward Jacobian; its runtime counterpart, RJCP, measures the same inverse-consistency error along optimization trajectories. We prove that D-IPG is first-order equivalent to damped Gauss-Newton under local pseudoinverse consistency, with deviation controlled by composition error and conditioning. Across seven PDE inverse-problem benchmarks, D-IPG outperforms standard baselines, achieves 94.8% mean success across the six-problem reliability suite, and reaches comparable or better recovery quality at up to 77x lower inference-time solve cost on the main benchmarks.

Aaditya L. Kachhadiya

Independent Researcher Surat, India

Bayesian Inverse Problems Solver Iterations Jacobian Shape Optimization Curvature-Aware Spectral Framework

May 11, 2026eess.SP

Adaptive Diagonal Loading using Krylov Subspaces for Robust Beamforming

Reliable adaptive beamforming is critical for large microphone arrays operating in highly dynamic acoustic environments. In scenarios characterized by fast-moving talkers and interferers, the available sample support for estimating the spatial correlation matrix is often snapshot-deficient. This deficiency degrades the White Noise Gain (WNG), leading to severe target signal cancellation. To ensure stable and robust beamforming, we previously proposed an adaptive diagonal loading method that leverages the Kantorovich inequality to guarantee the WNG remains strictly within specified bounds. However, accurately determining the smallest necessary loading level requires calculating the extreme eigenvalues of the spatial correlation matrix, a computationally expensive

\mathcal{O}(M^3)

operation for large arrays. In this paper, we introduce a highly efficient

\mathcal{O}(kM^2)

estimation technique using Lanczos iterations to build a small Krylov subspace. By projecting the correlation matrix onto a tridiagonal matrix of dimension

k \ll M

, we extract Ritz values that rapidly converge to the exact extreme eigenvalues. Our evaluations demonstrate that this Lanczos-accelerated approach achieves performance identical to exact Eigenvalue Decomposition (EVD), ensuring optimal interference suppression and strict WNG adherence at a fraction of the computational cost.

Manan Mittal, Ryan M. Corey, John R. Buck +1

Stony Brook University1 · University of Illinois Chicago2 · University of Massachusetts Dartmouth3

Beamforming Microphone Array Spectral Decomposition Solver Iterations Interference

May 11, 2026stat.ML

Affine Tracing: A New Paradigm for Probabilistic Linear Solvers

Probabilistic linear solvers (PLSs) return probability distributions that quantify uncertainty due to limited computation in the solution of linear systems. The literature has traditionally distinguished between Bayesian PLSs, which condition a prior on information obtained from projections of the linear system, and probabilistic iterative methods (PIMs), which lift classical iterative solvers to probability space. In this work we show this dichotomy to be false: Bayesian PLSs are a special case of non-stationary affine PIMs. In addition, we prove that any realistic affine PIM is calibrated. These results motivate a focus on (non-stationary) affine PIMs, but their practical adoption has been limited by the significant manual effort required to implement them. To address this, we introduce affine tracing, an algorithmic framework that automatically constructs a PIM from a standard implementation of an affine iterative method by passing symbolic tracers through the computation to build an affine computational graph. We show how this graph can be transformed to compute posterior covariances, and how equality saturation can be used to perform algebraic simplifications required for computation under specific prior choices. We demonstrate the framework by automatically generating a probabilistic multigrid solver and evaluate its performance in the context of Gaussian process approximation.

Disha Hegde, Marvin Pförtner, Jon Cockayne

University of Southampton · AI Center, University of Tübingen

Solver Iterations Affine Algebraic Multigrid Least Squares

May 9, 2026cs.MS

cuRegOT: A GPU-Accelerated Solver for Entropic-Regularized Optimal Transport

Optimal transport (OT) has emerged as a fundamental tool in modern machine learning, yet its computational cost remains a significant bottleneck for large-scale applications. While harnessing the massive parallelism of modern GPU hardware is critical for efficiency, the de facto standard Sinkhorn algorithm, despite its ease of parallelization, often suffers from slow convergence in challenging problems. More recently, the sparse-plus-low-rank quasi-Newton method offers a balance between convergence rate and per-iteration complexity; however, its efficiency on GPUs is severely hindered by the serial nature of sparse matrix symbolic analysis and irregular memory access patterns. To bridge this gap, we present cuRegOT, a high-performance GPU solver tailored for entropic-regularized OT. We introduce a suite of algorithmic and architectural optimizations, including an amortized symbolic analysis strategy to mitigate CPU bottlenecks, an asynchronous Sinkhorn iterates generation mechanism, and a fused kernel for bandwidth-efficient gradient evaluation. These strategies are backed by rigorous theoretical guarantees ensuring algorithmic convergence. Extensive numerical experiments demonstrate that cuRegOT achieves significant speedups over state-of-the-art GPU-based solvers across a variety of benchmark tasks.

Yixuan Qiu

School of Statistics and Data Science & Institute of Big Data Research Shanghai University of Finance and Economics

Fused Sinkhorn-Localized Similarity Cpu-Gpu Hybrid Designs Solver Iterations Newton Solvers

May 8, 2026math.NA

NSPOD: Accelerating Krylov solvers via DeepONet-learned POD subspaces

The convergence of Krylov-based linear iterative solvers applied to parametric partial differential equations (PDEs) is often highly sensitive to the domain, its discretization, the location/values of the applied Dirichlet/Neumann boundary conditions, body forces and material properties, among others. We have previously introduced hybridization of classical linear iterative solvers with neural operators for specific geometries, but they tend to not perform well on geometries not previously seen during training. We partially addressed this challenge by introducing the deep operator network Geo-DeepONet and hybridizing it with Krylov-based iterative linear solvers, which, despite learning effectively across arbitrary unstructured meshes without requiring retraining, led to only modest reductions in iterations compared to state-of-the-art preconditioners. In this study we introduce Neural Subspace Proper Orthogonal Decomposition (NSPOD), a multigrid-like deep operator network-based preconditioner which can dramatically reduce the number of iterations needed for convergence in Krylov-based linear iterative solvers, even when compared to state-of-the-art methods such as algebraic multigrid preconditioners. We demonstrate its efficiency via numerical experiments on a linearized version of solid mechanics PDEs applied to unstructured domains obtained from complex CAD geometries. We expect that the findings in this study lead to more efficient hybrid preconditioners that can match, or possibly even surpass, the convergence properties of the current gold standard preconditioning methods for solid mechanics PDEs.

Francesc Levrero-Florencio, Youngkyu Lee, Jay Pathak +1

Synopsys UK Ltd, UK · Division of Applied Mathematics, Brown University, Rhode Island, USA · Synopsys Inc, California, USA

Algebraic Multigrid Solver Iterations Spectral Preconditioning Subspace Mesh Generation

May 7, 2026cs.CV

Solving Minimal Problems Without Matrix Inversion Using FFT-Based Interpolation

Estimating camera geometry typically involves solving minimal problems formulated as systems of multivariate polynomial equations, which often pose computational challenges when using existing Gröbner-basis or resultant-based methods due to matrix inversion needed in the online solver. Here we propose a sampling-based, matrix inversion-free method that constructs the solvers using sparse hidden-variable resultants. The determinant polynomial in the hidden variable is efficiently reconstructed via inverse fast Fourier transform interpolation from sampled evaluations, avoiding symbolic expansion. Solving this polynomial yields the hidden variable, and the remaining unknowns are recovered by identifying rank-1 deficient submatrices and applying Cramer's rule. A greatest common divisor-based criterion ensures robust submatrix identification under noise. Experiments on diverse minimal problems demonstrate that the proposed solver achieves strong numerical stability and competitive runtime, particularly for small-scale problems, providing a practical alternative to traditional Gröbner-basis and resultant-based solvers.

Haidong Wu, Snehal Bhayani, Janne Heikkilä

Center for Machine Vision and Signal Analysis University of Oulu, Oulu, Finland

Solver Iterations Matrix Polynomials Interpolation

May 6, 2026math.ST

Direct Estimation of Schrödinger Bridge Time-Series Drifts: Finite-Sample, Asymptotic, and Adaptive Guarantees

We study nonparametric estimation of Schrödinger bridge (SB) drifts from i.i.d.\ data observed on a single time interval. Starting from the conditional-ratio form of the Schrödinger bridge time-series (SBTS) drift formula, we analyze a direct Nadaraya--Watson plug-in estimator built from kernelized numerator and denominator terms. Unlike recent SB analyses based on entropic-OT potentials, Sinkhorn iterations, or iterative bridge solvers, our approach works directly at the drift level and isolates \emph{statistical error} from optimization, approximation, and discretization error. Under Hölder regularity, a marginal-density floor, and bounded support, we prove a uniform non-asymptotic bound for admissible bandwidth pairs, a pointwise CLT under genuine undersmoothing, and an adaptive bandwidth selector satisfying an oracle inequality. We also prove a pivot-local minimax lower bound which, through an explicit uniform pivot, yields a global minimax lower bound under transparent compatibility conditions; hence the adaptive selector is minimax-rate optimal up to logarithmic factors. Synthetic experiments provide theorem-targeted diagnostics for finite-sample scaling, Gaussian approximation, and adaptive behavior.

Othmane Mazhar, Huyên Pham

LPSM, Sorbonne Université and Université Paris Cité · École Polytechnique, CMAP

Schrödinger Bridges Finite-Sample Solver Iterations Entropy

May 2, 2026cs.LG

Decision-Focused Learning via Tangent-Space Projection of Prediction Error

Decision-Focused Learning (DFL) trains predictors to improve downstream decision quality, but computing regret gradients typically requires differentiating through solvers or relying on surrogate losses, which can be computationally expensive or deviate from the true objective. We show that, under standard regularity with locally stable active constraints, the regret gradient admits a closed-form geometric characterization, equivalent to the prediction error projected onto the tangent space of active constraints, scaled by local curvature. This reveals that regret gradients can be obtained by filtering decision-irrelevant components from the MSE gradient, providing a simpler and more direct alternative to existing approaches. Based on this, we propose PEAR (Projected Error As Regret-gradient), which computes regret gradients via a reduced linear system over active constraints, avoiding differentiation through solver iterations or additional optimization solves. Experiments on LP benchmarks and a real-world QP task show that PEAR achieves the best decision quality among all baselines while being the most computationally efficient, with gains that persist under constraint shifts.

Junhyeong Lee, Sangjin Jin, Yongjae Lee

Department of Industrial Engineering, Ulsan National Institute of Science and Technology, Ulsan, South Korea.

Plan Regret Prediction Error Structural Constraints Solver Iterations

Apr 21, 2026cs.GR

An Efficient Multilevel Preconditioned Nonlinear Conjugate Gradient Method for Incremental Potential Contact

Incremental Potential Contact (IPC) guarantees intersection-free simulation but suffers from high computational costs due to the expensive Hessian assembly and linear solves required by Newton's method. While Preconditioned Nonlinear Conjugate Gradient (PNCG) avoids Hessian assembly, it has historically struggled with poor convergence in stiff, contact-rich scenarios due to the lack of effective preconditioners; simple Jacobi preconditioners fail to capture the global coupling, while advanced hierarchy-based preconditioners like Multilevel Additive Schwarz (MAS) are computationally prohibitive to rebuild at every nonlinear iteration. We present MAS-PNCG, a method that unlocks the power of hierarchical preconditioning for nonlinear optimization. Our key technical innovation is a Sparse-Input Woodbury update algorithm that incrementally adapts the fine-level MAS components to rapidly evolving contact sets. This bypasses the need for full preconditioner rebuilds, reducing maintenance cost to near-zero while capturing the complex spectral properties of the contact system. Furthermore, we replace heuristic PNCG search directions with a Hessian-aware 2D subspace minimization that optimally combines the preconditioned gradient and previous direction. We also apply a fast per-subdomain conservative CCD method that ensures penetration-free trajectories while avoiding overly restrictive global step sizes. Experiments demonstrate that our MAS-PNCG outperforms state-of-the-art Newton-PCG solvers, GIPC and StiffGIPC, both preconditioned with MAS up to 5.66

\times

and 2.07

\times

respectively.

Yu Zhang, Xing Shen, Kemeng Huang +5

S-Lab, Nanyang Technological University, Singapore · Shanghai AI Laboratory, China · University of Hong Kong, China +3

Solver Iterations Spectral Preconditioning Gradient

Apr 17, 2026math.NA

Towards Universal Convergence of Backward Error in Linear System Solvers

The quest for an algorithm that solves an

n\times n

linear system in

O(n^2)

time complexity, or

O(n^2 \text{poly}(1/ε))

when solving up to

ε

relative error, is a long-standing open problem in numerical linear algebra and theoretical computer science. There are two predominant paradigms for measuring relative error: forward error (i.e., distance from the output to the optimum solution) and backward error (i.e., distance to the nearest problem solved by the output). In most prior studies, convergence of iterative linear system solvers is measured via various notions of forward error, and as a result, depends heavily on the conditioning of the input. Yet, the numerical analysis literature has long advocated for backward error as the more practically relevant notion of approximation. In this work, we show that -- surprisingly -- the classical and simple Richardson iteration incurs at most

1/k

(relative) backward error after

k

iterations on any positive semidefinite (PSD) linear system, irrespective of its condition number. This universal convergence rate implies an

O(n^2/ε)

complexity algorithm for solving a PSD linear system to

ε

backward error, and we establish similar or better complexity when using a variety of Krylov solvers beyond Richardson. Then, by directly minimizing backward error over a Krylov subspace, we attain an even faster

O(1/k^2)

universal rate, and we turn this into an efficient algorithm, MINBERR, with complexity

O(n^2/\sqrtε)

. Finally, we extend this approach via normal equations to solving general linear systems in

O(n^2\log(n)/ε)

time complexity. We report strong numerical performance of our algorithms on benchmark problems.

Michał Dereziński, Yuji Nakatsukasa, Elizaveta Rebrova

University of Michigan, Ann Arbor, MI · University of Oxford, Oxford, UK · Princeton University, Princeton, NJ

Solver Iterations Convergence Backward Passes Optimal Sample Complexity

Mar 23, 2026math.NA

Stability and Bifurcation Analysis of Nonlinear PDEs via Random Projection-based PINNs: A Krylov-Arnoldi Approach

We present a numerical framework for the stability and bifurcation analysis of nonlinear partial differential equations (PDEs) in which the solutions are sought in the function space spanned by physics-informed random projection neural networks (PI-RPNNs), and discretized via a collocation approach. These are shallow, single-hidden-layer networks with randomly sampled and fixed a priori hidden-layer weights; only the linear output layer weights are optimized, reducing training to a single least-squares solve. This linear output structure enables the explicit generalized eigenvalue formulation governing the linear stability of stationary solutions, separating the physical domain interior dynamics from the algebraic constraints imposed by boundary conditions (BCs), at no additional training cost and without requiring additional PDE solves. However, the random projection collocation matrix is inherently numerically rank-deficient, rendering naive eigenvalue computation unreliable and contaminating the true eigenvalue spectrum with spurious near-zero modes. To overcome this limitation, we introduce a matrix-free shift-invert Krylov-Arnoldi method that operates directly in weight space, avoiding explicit inversion of the numerically rank-deficient collocation matrix and enabling the reliable computation of several leading eigenpairs of the physical Jacobian - the discretized Frechet derivative of the PDE operator with respect to the solution field, whose eigenvalue spectrum determines linear stability. We further prove that the PI-RPNN-based generalized eigenvalue problem is almost surely regular, guaranteeing solvability with standard eigensolvers, and that the singular values of the random projection collocation matrix decay exponentially for analytic activation functions, rigorously explaining the observed ill-conditioning.

Gianluca Fabiani, Michail E. Kavousanakis, Constantinos Siettos +1

Hopkins Extreme Materials Institute, Johns Hopkins University, Baltimore, 21218, MD, USA · Dept. of Chemical and Biomolecular Engineering, Johns Hopkins University, Baltimore, 21218, MD, USA · School of Chemical Engineering, National Technical University of Athens, 9 Iroon Polytechniou Street, Athens, 15772, Greece +2

Partial Differential Equations Parametric Physics-Informed Neural Network Solver Iterations Bifurcation Physics-Informed Learning Random

Mar 16, 2026cs.LG

Unbiased and Biased Variance-Reduced Forward-Reflected-Backward Splitting Methods for Stochastic Composite Inclusions

This paper develops new variance-reduction techniques for the forward-reflected-backward splitting (FRBS) method to solve a class of possibly nonmonotone stochastic composite inclusions. Unlike unbiased estimators such as mini-batching, developing stochastic biased variants faces a fundamental technical challenge and has not been utilized before for inclusions and fixed-point problems. We fill this gap by designing a new framework that can handle both unbiased and biased estimators. Our main idea is to construct stochastic variance-reduced estimators for the forward-reflected direction and use them to perform iterate updates. First, we propose a class of unbiased variance-reduced estimators and show that increasing mini-batch SGD, loopless-SVRG, and SAGA estimators fall within this class. For these unbiased estimators, we establish a

\mathcal{O}(1/k)

best-iterate convergence rate for the expected squared residual norm, together with almost-sure convergence of the iterate sequence to a solution. Consequently, we prove that the best oracle complexities for the

n

-finite-sum and expectation settings are

\mathcal{O}(n^{2/3}ε^{-2})

and

\mathcal{O}(ε^{-10/3})

, respectively, when employing loopless-SVRG or SAGA, where

ε

is a desired accuracy. Second, we introduce a new class of biased variance-reduced estimators for the forward-reflected direction, which includes SARAH, Hybrid SGD, and Hybrid SVRG as special instances. While the convergence rates remain valid for these biased estimators, the resulting oracle complexities are

\mathcal{O}(n^{3/4}ε^{-2})

and

\mathcal{O}(ε^{-5})

for the

n

-finite-sum and expectation settings, respectively. Finally, we conduct two numerical experiments on AUC optimization for imbalanced classification and policy evaluation in reinforcement learning.

Quoc Tran-Dinh, Nghia Nguyen-Trung

Department of Statistics and Operations Research The University of North Carolina at Chapel Hill 318 Hanes Hall, UNC-Chapel Hill, NC 27599-3260.

Variance Reduction Solver Iterations Convergence

Jan 2, 2026cs.LG

Precision autotuning for linear solvers via contextual bandit-based RL

We propose a reinforcement learning (RL) framework for \xy{responsive} precision tuning for linear solvers, which can be extended to general algorithms. The framework is formulated as a contextual bandit problem and solved using incremental action-value estimation with a discretized state space to select optimal precision configurations for computational steps, \xy{retaining} precision and computational efficiency. To verify its effectiveness, we apply the framework to iterative refinement for solving linear systems

Ax = b

. In this application, our approach dynamically chooses precisions based on calculated features from the system while maintaining acceptable accuracy and convergence. In detail, an action-value estimator takes discretized features (e.g., approximate condition number and matrix norm) as input and outputs estimated action values, from which a policy selects the actions (chosen precision configurations for specific steps), optimized via an

ε

-greedy strategy to maximize a multi-objective reward to balance accuracy and computational cost. Empirical results demonstrate effective precision selection, \xy{increasing the use of lower-precision arithmetic} while maintaining accuracy comparable to double-precision baselines. \xy{We further evaluate the learned policies in a compiled CPU GMRES-IR implementation using FP16, FP32, and FP64 arithmetic for solver-level native validation.} The framework generalizes to diverse out-of-sample data and provides insights into applying RL precision selection to other numerical algorithms, advancing mixed-precision numerical methods in scientific computing. To the best of our knowledge, this is the first work on precision autotuning with RL with verification on unseen datasets.

Erin Carson, Xinye Chen

Charles University Prague, Czech · Sorbonne Université, CNRS, LIP6 Paris, France

Solver Iterations Full-Precision Performance Precision Offline Reinforcement Learning Greedy Algorithm Value Functions

Dec 1, 2023stat.ME

Multiple Testing of Linear Forms for Noisy Matrix Completion

Many important tasks of large-scale recommender systems can be naturally cast as testing multiple linear forms for noisy matrix completion. These problems, however, present unique challenges because of the subtle bias-and-variance tradeoff of and an intricate dependence among the estimated entries induced by the low-rank structure. In this paper, we develop a general approach to overcome these difficulties by introducing new statistics for individual tests with sharp asymptotics both marginally and jointly, and utilizing them to control the false discovery rate (FDR) via a data splitting and symmetric aggregation scheme. We show that valid FDR control can be achieved with guaranteed power under nearly optimal sample size requirements using the proposed methodology. Extensive numerical simulations and real data examples are also presented to further illustrate its practical merits.

Wanteng Ma, Lilun Du, Dong Xia +1

Department of Statistics and Data Science, University of Pennsylvania, USA · Department of Decision Analytics and Operations, City University of Hong Kong, Hong Kong · Department of Mathematics, Hong Kong University of Science and Technology, Hong Kong +1

Tensor Completion False Discovery Rate Sample Size Noisy Data Solver Iterations Bias-Variance Decomposition

Date pendingstat.ME

Scalable Krylov Subspace Methods for Generalized Mixed-Effects Models with Crossed Random Effects

Mixed-effects models are widely used to model data with complex grouping structures and high-cardinality categorical predictor variables. However, for high-dimensional crossed random effects, current standard computations relying on Cholesky decompositions can become prohibitively slow. In this work, we present Krylov subspace-based methods that address existing computational bottlenecks, and we analyze them both theoretically and empirically. In particular, we derive new results on the convergence and accuracy of the preconditioned stochastic Lanczos quadrature and conjugate gradient methods for mixed-effects models, and we develop scalable methods for calculating predictive variances. In experiments with simulated and real-world data, the proposed methods yield speedups of several orders of magnitude and are more computationally robust than Cholesky-based computations, while maintaining essentially the same accuracy.

Pascal Kündig, Fabio Sigrist

Lucerne University of Applied Sciences and Arts · University of Basel · Seminar for Statistics, ETH Zurich

Generalized Linear Model Solver Iterations Singular Value Decomposition Bayesian Quadrature Newton