Papers

This paper introduces four groups of subspace methods for nonlinear monotone equations, with applications to large-scale machine learning problems. The methods use Jacobian-free subspace (JFS) directions of conjugate-gradient type, combined with either fixed step sizes or variable step sizes generated by the projected method of Solodov and Svaiter. To ensure convergence independently of the specific algebraic form of the subspace directions, we impose an angle condition together with an explicit scaling rule controlling the effective search directions. Under monotonicity and Lipschitz continuity of the operator, we establish global convergence for both the line-search and fixed-step frameworks, as well as a best-iterate residual rate $O(\ell^{-1/2})$. Under a local error bound, the distance to the solution set satisfies the sharper decay $o(\ell^{-1/2})$. If the operator is continuously differentiable and its Jacobian is nonsingular at a solution, the required local error bound and local isolation follow, yielding $R$-linear local convergence. The residual sequence then converges geometrically and hence satisfies the last-iterate rate $o(\ell^{-1})$, without strong monotonicity. We also derive iteration and residual-evaluation complexity bounds: $O(\varepsilon^{-2})$ for the baseline best-iterate guarantee and $O(\log(\varepsilon^{-1}))$ in the local linear regime, together with a uniform bound on backtracking residual evaluations. Under additional asymptotic assumptions, the proposed JFS directions and several classical update directions admit related optimistic gradient descent–ascent (OGDA)-type residual–memory representations. Numerical experiments illustrate the robustness and efficiency of the methods on representative min–max problems.

Intermodal freight networks are increasingly exposed to correlated, multi-mode disruptions, yet resilience assessments often rely on historical or uncorrelated scenarios that understate systemic risk. This paper develops a data-driven stress-testing framework integrating generative adversarial networks (GANs) with an intermodal optimization model to evaluate performance under realistic compound disruptions. The case study examines weather-related disruptions in the Tennessee Valley corridor. Each GAN-generated scenario is used as a simulation input, and the resulting routing problem is solved to obtain system costs. Aggregating outcomes enables estimation of expected costs and identification of major risk drivers. Results show that historical disruptions increase total cost by about 3%, whereas GAN-generated scenarios raise costs by over 25%, producing an expected annual cost of $5.11 million. Risk is concentrated in correlated multi-node failures and critical nodes such as the Port of Knoxville. The framework helps identify vulnerabilities and prioritize resilience investments.

Efficient resource use in packing and assembly-line applications requires decisions that jointly account for capacity and precedence constraints. The strongly NP-hard bin packing problem with generalized precedence constraints (BPP-GP) models such decisions by minimizing the number of ordered, capacitated bins required to pack weighted items, even when precedence requirements span multiple bins. Existing exact algorithms primarily focus on classical special cases, whereas general BPP-GP has been addressed only via compact integer models and heuristics, with no efficient open-source exact solver. We present PrecPack, a unified exact solver that extends branch-bound-and-remember (BBR) to arbitrary nonnegative precedence weights and naturally specializes to the classical cases. Generalized states capture restrictions that remain active across future bins, which are addressed through branching, dominance, and conflict-aware lower bounds. Root column generation uses fixed-point arithmetic to compute numerically valid dual bounds for pruning or to prove optimality. To support reuse and verification, we provide common programming and command-line interfaces, independent assignment checking, explicit termination statuses, and reproducible batch execution; the core procedures require no commercial software. In same-machine, single-threaded comparisons on classic assembly-line benchmarks, more instances are proven optimal, and average computing times are substantially reduced relative to leading source-available BBR implementations. Further comparisons with published benchmark results for bin packing with precedence constraints and BPP-GP also show that more instances were proved optimal and that reported average gaps were smaller on most benchmark sets. PrecPack is released under the MIT License at https://github.com/Sunkanghong-Wang/PrecPack.

Recently, dynamic state feedback controllers that are based on dynamic extensions have been presented for heterodirectional hyperbolic systems. In this paper, a similar concept for the control of coupled diffusion-reaction systems is suggested. The introduction of a specific controller dynamics leads to homogenized diffusion coefficients for the extended system. Then, a backstepping-based static state feedback for the dynamically extended system is designed, which, overall, results in a dynamic state feedback. Such a design allows stabilizing a more general class of parabolic systems as well as assigning arbitrary closed-loop dynamics. This can be used, e.g., to achieve a decoupled input-output behavior, which is, in general, not possible with a static state feedback. A simulation example illustrates the results.

In flow matching (FM), a velocity model $v_{\theta}$ is trained using a predefined path $g_t$ that connects data and noise samples (e.g., $g_t(x_0, x_1) = (1 - t) x_0 + t x_1$). In this work, we study the choice of this path from an optimization perspective by analyzing the variance of stochastic gradients. We consider the class $G(p_t,v^\star_t)$ of paths that induce the same marginal distributions $p_t$ and marginal velocity field $v^\star_t$, and therefore the same FM objective. Our main finding is that the choice of path $g_t$ can fundamentally change the convergence rate of SGD, even when the FM objective remains exactly the same. (i) For a linear velocity model and one-dimensional Gaussian data, we derive a tight bound on the SGD iteration complexity up to logarithmic factors and find an analytically optimal path that minimizes this bound among linear paths inducing the same FM problem. (ii) We then extend the variance analysis to general FM problems and formulate path selection at a fixed $\theta$ as the variance-minimization problem PathOpt$_\theta$, constrained to $g_t\in G(p_t,v^\star_t)$. We show that this constraint is essential: reducing variance without it can lead to slower convergence. (iii) Since the constraint $g_t \in G(p_t,v^\star_t)$ cannot generally be verified directly, we derive an equivalent formulation with constraints that can be estimated from samples, allowing paths to be found numerically. Our theoretical results are supported by experiments with Gaussian data, Gaussian mixture models, and real datasets.

We study a state-constrained Bolza optimal control problem governed by a non-local continuity equation in the Wasserstein space $\sP_2(\R^d)$. We develop a Pontryagin maximum principle for these problems in which both the cost functionals and the constraints are only locally Lipschitz, and the initial measure does not necessarily have compact support. To derive the main result, we introduce a notion of a Clarke-type deterministic subdifferential over $\sP_2(\R^d)$. The proof of the main result combines a careful adaptation of the Lasry-Lions regularization to the Wasserstein space coupled with stability results of the Clarke subdifferential that allow us to recover the first-order information. We also show that, for functionals convex along plans, this subdifferential admits a global supporting characterization and provides a first-order characterization of global minimizers.

Drifting methods have recently been introduced as a promising new paradigm for training generative models, and have attracted a lot of attention so far. In this work, we focus on how the optimization procedure shapes their convergence properties, shedding light on the underlying machinery that may explain their empirical success. More specifically, when the drift field is a Wasserstein gradient, optimizing the drifting loss via plain Euclidean gradient descent seems to forfeit any hope of exploiting convexity and therefore of obtaining any convergence guarantees. Changing either the drift field or the optimization procedure, however, may lead to different dynamics: we show that keeping the Wasserstein gradient as drift field but replacing gradient descent with a Gauss-Newton scheme amounts to performing a natural gradient descent in parameter space. This connection allows us to provide global convergence guarantees for the optimization procedure under suitable convexity assumptions.

We study an optimal intervention problem for linear threshold models. This is a popular class of dynamical network systems whereby a number of agents, identified with the nodes of a graph, strategically change their binary action (0 or 1) according to a threshold rule. Specifically, an agent adopts action 1 if and only if the fraction of its neighbors in the interaction graph that do so is greater than or equal to a prescribed threshold. Assuming that a planner can modify the agents' thresholds at a cost equal to the aggregate threshold increase, we study the minimum intervention cost needed to ensure global convergence to the all-1 configuration. Our main contribution is the introduction of a new graph-theoretic quantity, called oriented path number, that is the minimum number of disjoint paths needed to cover the graph that can be oriented to form a directed acyclic graph. When thresholds are all equal to 1/2, the optimal cost is shown to coincide with the oriented path number, whereas, in the general case, it turns out to be the main ingredient of a bound on the optimal intervention cost.

The choice of Riemannian metric can strongly influence the convergence of gradient-based optimization over covariance matrices. Euclidean, Bures-Wasserstein and affine-invariant metrics are common choices, but their relative effectiveness depends on the objective. We introduce a two-parameter family defined by $X^{p}LX^{q}+X^{q}LX^{p}=U$, solved for $L$ at each tangent vector $U$, that contains all three as exact members, at $(0,0)$, $(1,0)$ and $(1,1)$, and extends past them. We treat the choice of member as a particular way of preconditioning for a given problem. To this end, we analyze the conditioning of the Riemannian Hessian at the solution. We show that it obeys a lower bound that depends on $(p,q)$ only through the exponent $r=p+q$. When the Euclidean Hessian is a pure power that mixes no eigendirections, the member $p=q=r/2$ attains that bound, and a closed-form criterion identifies the other members that do. We discuss ways to tune $r$ for a given problem. Experiments on real covariance data confirm the predicted conditioning and the benefit of tuning $r$. A task covariance example shows a further gain from tuning the shape.

We study the problem of minimizing a weighted sum of squared Wasserstein distances with signed coefficients. In the case of a single positive coefficient, we derive two convex dual formulations: one expressed in terms of Brenier potentials, and one as a projection problem in convex order. Such dual formulations arise via Toland duality by exploiting different notions of convexity within the problem. They generalize those recently established for the metric extrapolation problem, which corresponds to signed barycenters between only two measures. For general signed barycenters, we identify necessary optimality conditions that reduce the problem to the case with one positive coefficient. On the numerical side, we propose a grid-free, particle-based algorithm built on entropic regularization and Sinkhorn iterations to compute signed barycenters with arbitrary coefficients.

Particle swarm optimization (PSO) is a widely used algorithm featured in many state-of-the-art optimization tool-kits. However, rigorous performance guarantees are still lacking. The standard PSO dynamics do not admit a natural mean-field description, which would provide an avenue for theoretical analysis. By modifying the PSO formulation, one can recover the consensus-based optimization (CBO) algorithm with memory, which admits a mean-field limit and facilitates rigorous convergence analysis. These theoretical guarantees rely heavily on the fact that for CBO, the drift and noise strengths can be chosen independently, whereas they are coupled for PSO. We analyze how the PSO parameter coupling affects existing convergence guarantees for CBO and its variant with memory effect. We show, by an explicit construction, that the coupling still leaves a non-empty set of admissible parameters for these convergence guarantees to hold. However, the admissible parameter ranges shrink in the limits used to recover PSO. The resulting convergence guarantees from CBO therefore do not directly extend to the classical PSO model. We provide numerical simulations illustrating the parameter tradeoffs shown in the theoretical analysis.

We introduce the Unified Separation Condition (USC) - a simple and general constraint qualification for finite-dimensional nonlinear programming. The USC requires that the origin be uniformly separated from the subdifferential of the constraint violation function in a neighborhood of the feasible set. The main result of this paper is two-fold. First, we show that the classical Mangasarian–Fromovitz constraint qualification (MFCQ) for inequalities and the linear independence constraint qualification (LICQ) for equalities imply the USC. Conversely, USC is strictly more general: it remains applicable in nonsmooth settings where classical conditions are not defined, and it may hold even when MFCQ fails. The proof that MFCQ implies USC is based on Gordan's theorem; the implication from LICQ to USC follows from the linear independence of the gradients of the active constraints. We also introduce a local version of USC and discuss its computational verification via convex quadratic programming. The USC provides a unified framework for constraint qualifications, bridging classical smooth theory and nonsmooth optimization.

Feature-based multi-product pricing uses customer characteristics to identify demand heterogeneity and tailor prices across products. Choice model trees segment customers through interpretable feature rules and fit a demand model within each leaf. Existing methods typically construct these trees greedily, selecting one myopic split at a time. We develop optimal choice model trees with multinomial logit leaves (OCMT-MNL), jointly optimizing the tree and leaf models within a prescribed depth. Our exact dynamic program derives closed-form Fenchel lower bounds during constrained Newton iterations and propagates them across nested and disjoint customer subsets, avoiding new fits and resuming unfinished fits without repeating completed work. In synthetic experiments, it reduces exact leaf fits by 99.98% and leaf evaluations by 86.13%, achieving up to 7.15-fold speedups over unpruned dynamic programming. One-dimensional lookup tables translate offline estimation into real-time pricing, with a revenue-loss bound quadratic in grid spacing under the fitted model. Compared with greedy trees, OCMT-MNL achieves lower revenue loss with fewer leaves on synthetic data and better predictive fit on real data. In a 23-week randomized experiment on ancillary seat pricing across 48 airline markets and 190,220 passengers, OCMT-MNL increases seat revenue per passenger by a statistically significant 11.3% over static pricing.

A complete characterization of the stability-preserving exponent set for fractional Hadamard powers of a monic Hurwitz-stable polynomial f of degree five is presented. The stability problem is reduced to the analysis of a single scalar function depending only on three parameters formed from the coefficients of the polynomial f. This reduction leads to a unique stability threshold $p_*(f) \in (0,1)$ such that the $p$-th Hadamard power of $f$ is Hurwitz stable if and only if $p > p_*(f)$. Consequently, the stability-preserving exponent set is precisely $(p_*(f),\infty)$. This threshold depends smoothly on the parameters and provides a global coordinate on the admissible parameter region. Finally, smooth dependence is illustrated by a one-parameter family of Hurwitz-stable polynomials whose complex-conjugate zeros approach the imaginary axis.

We prove that unrestricted static output-feedback stabilization is NP-hard, already in the single-input multiple-output case. More precisely, we give a polynomial-time many-one reduction from the NP-complete Betweenness problem. Given an instance with $n$ elements and $m$ betweenness constraints, the reduction constructs integer matrices $A\in\mathbb Z^{D\times D}$, $B\in\mathbb Z^{D\times1}$, and $C\in\mathbb Z^{n\times D}$ such that the instance is satisfiable if and only if there exists an unrestricted real row vector $K\in\mathbb{R}^{1\times n}$ for which $A+BKC$ is Hurwitz. The state dimension is $D=6(m+n)$, and every matrix entry has logarithmic bit length in the input size. The reduction uses one universal cubic source polynomial whose Routh determinant enforces the sign condition encoding betweenness, while fixed anchor values provide a mechanism for excluding spurious stabilizing gains of arbitrarily large magnitude. Frequency separation then compresses the simultaneous source conditions into one common affine polynomial, which is realized directly as the characteristic polynomial of $A+BKC$. As a by-product, deciding whether an affinely parameterized rational polynomial family contains a Hurwitz member on a rational box is NP-hard.

Increasingly severe and unpredictable extreme weather events-such as high-intensity forest firesthreaten long-term ecosystem service provision, driving a critical reassessment of current forest management planning. While planning frameworks exist, few incorporate optimization methods that account for the uncertainty and variability introduced by such extremes. Meanwhile, process-based models have advanced our understanding of fire impacts on forest dynamics and ecosystem services, yet their potential to inform strategic management decisions remains largely untapped. Focusing on intensive plantation forests, we propose a combined approach that uses fire simulation within a growth model to optimize forested landscape management under increasing fire risk. We demonstrate that incorporating fire behavior and its uncertainties significantly alters management outcomes, favoring diversified portfolios that balance high-yield but fire-prone regimes with lowerrisk alternatives. While diversification improves robustness for objectives like timber provision, it does not guarantee complementary benefits for carbon storage. Our method exposes the limitations of current management regimes and identifies key areas for improved planning to better mitigate fire risk and strengthen long-term resilience.

Penalizing a state constraint creates a singular localization problem: whole-space values may grow like the inverse penalty parameter, while numerical diffusion allows even an inward feedback to cross the boundary. For deterministic discounted optimal control, we show how the penalty itself supplies confinement that offsets this growth. With mesh size $h$ and penalty parameter $\varepsilon$, an inward barrier bounds the additional cost of numerical leakage by $O(h/\varepsilon)$ for a monotone centered-difference scheme with vanishing viscosity. Under $h\le\varepsilon$, an occupation estimate then yields $O(h)$ localization error with a sufficient box margin logarithmic in $1/h$ and independent of $\varepsilon$. We extend the result to bounded squared-distance penalties and combine it with residual-based evaluation bounds and discounted policy-error propagation. Under coupled refinement with vanishing penalization and discretization errors, explicit conditions on evaluation errors and greedy gaps ensure convergence of the neural value approximations to the constrained value, allowing measurable, nonunique greedy selectors. Reference calculations isolate leakage and localization across mesh and penalty scales, and distinguish evaluation, iteration, and approximation errors. An explicit cylindrical state-constraint solution provides a benchmark in arbitrary dimension, tested up to dimension twenty. Paired obstacle-navigation experiments illustrate why policy-value accuracy must be assessed alongside sampled residuals when comparing raw-residual and finite-grid-assisted neural evaluation.

This paper presents a stability analysis of the generalized multi-source Weber problems under the influence of data perturbation within the framework of the Minkowski function. First, we establish explicit Lipschitz constants for the objective functions and optimal value functions. Then, we prove several properties of the global optimal solution sets. Furthermore, we provide sufficient conditions for the upper semicontinuity of global solution mappings and the inner semicontinuity of local solution mappings. Three illustrative examples are constructed. Our results give indirect answers to several open questions regarding the stability of local optimal solution mappings of the optimization problems discussed.

In this manuscript, we study the master equation for mean field game of controls, assuming only that the data are Lipschitz continuous in the measure variable. Accordingly, we propose a weaker notion of solution called weak solutions for the master equation. We establish the global well-posedness of the master equation for such solution under the Lasry-Lions monotonicity and displacement $\lambda$-monotonicity conditions, respectively. The arguments rely on the analysis of the associated Pontryagin forward-backward stochastic differential equation system, especially the stability in the measure variable. Finally, we review several existing notions of non-smooth solution from the literature, namely good solutions, weak viscosity solutions, Lipschitz solutions, monotone solutions, and examine their relationship with our weak solution. We show that under appropriate assumptions all of them are equivalent by their definitions, and thus our manuscript provides a unification of weak solution notions for the master equations derived from mean field game of controls.

We solve the four-dimensional Hamilton-Jacobi-Bellman (HJB) equation for two diffusively coupled Stuart-Landau-like oscillators to obtain full-state optimal feedback control. A sweep over coupling strength reveals a sharp change in the numerically optimal strategy: below a threshold coupling value, the controller steers the phase difference toward anti-phase while keeping both oscillators near the limit cycle, as reduced-order models would suggest. Above this threshold, the HJB solution changes qualitatively; the controller transiently collapses one oscillator's amplitude to near zero, thereby enabling large phase repositioning near the origin before rebuilding its amplitude. Direct gradient-based and stochastic optimization do not recover this lower-cost collapse trajectory from the initializations considered, suggesting that it occupies a region of the control landscape that is difficult to access by direct search. A joint sweep over nonisochronicity and coupling shows that collapse can occur even for an isochronous oscillator: phase repositioning near the origin can favor an off-cycle strategy. Nonisochronicity provides an additional energetic benefit through a phase-velocity surplus at small amplitude, quantitatively accounting for the observed reduction in control cost. Comparisons with uncoupled and coupled phase-reduced baselines show that phase models become increasingly inaccurate and cost significantly more energy for strong coupling. Results for coupled Van der Pol oscillators further demonstrate that exploitation of off-cycle dynamics is not specific to the Stuart-Landau-like oscillators.

We study reachability and exact controllability for systems driven by a nondecreasing left-continuous Stieltjes derivator whose measure may have continuous, singular-continuous and atomic parts, while flat intervals represent inactivity. We construct an evolution family and post-jump variation-of-constants formula, and characterize the reachable subspace as the essential span of input directions. For constant coefficients, the Kalman condition is necessary for every derivator and, under regressivity, sufficient when the derivator is continuous and nonconstant on a subinterval. If the measure is purely atomic with exactly $N$ positive atoms, the reachable subspace is the truncated Krylov space generated by $B,AB,\ldots,A^{N-1}B$; for a Kalman-controllable pair, the minimum atom count among purely atomic derivators equals the controllability index. At the threshold count, one vanishing mass with the others converging to positive finite limits causes Gramian degeneration and worst-case energy blow-up; at or above threshold, compact positive mass ranges yield uniform coercivity. A mixed double-integrator example illustrates the results.

Using a novel effective non-asymptotic concentration inequality, we establish a distribution-uniform generalization of Sturm's strong law of large numbers for inductive means of $L^1$-sequences of i.i.d. random variables on (separable) Hadamard spaces. Building on that, we establish a strong law of large numbers for integrable random monotone vector fields on (separable) Hilbert-Hadamard spaces, that is Hadamard spaces where all tangent cones isometrically embed into Hilbert spaces, extending a previous result of Salim set in Hilbert spaces. We use this latter strong law to establish a probabilistic Lie-Trotter-Kato formula for the gradient flow generated by an integral function over (separable) Hilbert-Hadamard spaces where all tangent cones are actually full Hilbert spaces. This application leverages the previous Lie-Trotter-Kato formula established for gradient flows of sums of convex functions by Stojković together with a result relating resolvent and gradient flow convergence established by Bačák, as well as a new result on the interchangeability of the subdifferential and the integral, which we establish for $L^2$-Lipschitz integrands. The last ingredient provides, to our knowledge, the first nonlinear version of (a particular case of) a previous result due variously to Ioffe and Tikhomirov, Levin, Hiriart-Urruty, Thibault as well as Rockafellar and Wets. Throughout the paper, we highlight various remaining open problems.

Operator splitting methods are a fundamental class of algorithms for solving structured monotone inclusion problems arising in optimization, control, and game theory. We consider the case, common in stochastic regimes, where the forward evaluation of one of the constituent operators is either unavailable in closed form or computationally expensive to evaluate, and is therefore approximated using a finite number of noisy oracle samples. We establish distribution-free, finite-sample certificates for the quality of the output produced by data-driven Davis-Yin splitting algorithms. Unlike previous works, our analysis directly controls the residual error via martingale inequalities instead of relying on algorithmic stability arguments for a tailored surrogate loss, yielding the first a priori certificates whose statistical excess provably vanishes with the sample size. We further show that, under linear convergence of the Davis-Yin splitting algorithm, the dependence of our bounds on the iteration count improves from linear growth to exponential decay. We validate our theoretical results on a stochastic portfolio optimization problem with uncertain asset returns.

This paper is concerned with multivariable stochastic Hamiltonian systems governed by an ordinary differential equation for the position and an Ito stochastic differential equation for the momentum. The momentum dynamics involve both conservative and nonconservative forcing including Langevin viscous damping and the Wiener process as a random force. The setting allows for rotational degrees of freedom, a non-quadratic potential energy function and position-dependent mass, diffusion and damping matrices. In addition to modelling nonlinear physical dynamics subject to random forcing, the system is motivated by a stochastic optimisation viewpoint, where the search for global minima of the potential energy (for example, using the gradient descent) is carried out through convergence to a Boltzmann equilibrium measure which favours smaller values of this function. Under certain regularity conditions, we show that the $\chi^2$-divergence or, equivalently, Renyi's second-order relative entropy with respect to the invariant measure has a nonpositive time derivative and thus provides a Lyapunov functional candidate for this convergence. However, this derivative vanishes at some moments of time, which is related to the behaviour of the probability distribution of the momentum conditioned on the position and has a link to a family of auxiliary quantum harmonic oscillator Hamiltonians parameterised by the position space of the underlying classical system. This position-momentum conditioning is used in order to prove that such pre-equilibrium break times in the entropy dissipation are isolated, thus making the $\chi^2$-divergence a strictly decreasing Lyapunov functional for the Fokker-Planck-Kolmogorov equation which governs the joint position-momentum distribution in such systems.

A fulfillment center (FC) is a specialized logistics facility where e-commerce inventory is received, stored, and processed. Its internal organization directly impacts space utilization and the performance of put-away and picking operations. Items are stored in bins following strict operational guidelines to ensure productivity in both activities, so bin dimensions must be tailored to the product profile of each facility to maximize space utilization. % We address the bin dimensioning problem arising in e-commerce fulfillment centers: given a set of candidate bin types and a large, heterogeneous inventory, determine the number of bins of each type and assign all items so as to minimize total bin volume, subject to operational constraints governing how products may be combined within a single bin. We develop a dimensional decomposition that exploits these constraints to reduce the three-dimensional packing problem to a one-dimensional block-positioning problem. Building on this decomposition, we formulate the problem as a mixed-integer linear program and propose two solution methods: a Best-Fit-Decreasing (BFD) heuristic and a column generation scheme that produces tight LP lower bounds. Computational experiments on four synthetic datasets, calibrated on proprietary data from an American e-commerce fintech and ranging from 1.5 thousand to 1.5 million SKUs, show that the BFD heuristic scales to instances with up to 1.7 million blocks, and that column generation certifies a BFD–LP gap of 2.26% on the smallest dataset. The experiments further reveal a structural divide in solution quality across datasets, driven by the interaction between operational constraints and each facility's product profile.

Recent research has focused on developing GPU-accelerated algorithms for solving linear programs (LPs), with results that are nothing short of extraordinary. But a complete solver pipeline consists of more than the core algorithm, and the elephant in the room is that stages such as presolve have remained on the CPU, where they have become an increasingly large bottleneck. Presolve has stayed on the CPU because of the longstanding view, held by academic and industrial developers alike, that it is inherently sequential and irregular, and therefore difficult to parallelize efficiently. In this paper we challenge this view and show that presolving can benefit substantially from GPU acceleration. We describe simple design principles that avoid severe load imbalance and expose parallelism in reductions that appear inherently sequential. We implement these principles in cuPSLP, a GPU-accelerated version of PSLP, a CPU-based presolver that is itself several times faster than a state-of-the-art commercial presolver. In our experiments, cuPSLP reduces the shifted geometric mean presolve time of PSLP by a factor of 11 on the Mittelmann LP benchmark set and 42 on the GAMS large-scale LP benchmark set, with similar reduction quality.

This paper extends the fixed-point bundle framework for finite-dimensional variational inequalities (VIs) from the simplex domain to the product-of-simplex domain, which is directly applicable to solving Nash equilibria. The fixed-point bundle for VIs on the product-of-simplex domain reveals a composite fiber bundle structure. The key innovation is to construct an equivalent VI on the simplex domain and establish the equivalence between the two fixed-point bundle frameworks via a fiber bundle isomorphism. Exploiting this geometric equivalence, the predictor-corrector path-following algorithm for the VI on the product-of-simplex domain is shown to inherit the convergence guarantee of the simplex-domain framework, namely, global convergence with linear gap reduction near solutions. Numerical experiments on 5600 randomly generated instances with dimensions ranging from 2-player 128-action to 128-player 2-action demonstrate robust performance. The algorithm converges in every tested instance.

Reinforcement learning with verifiable rewards (RLVR) improves the reasoning capabilities of large language models (LLMs). We introduce Bellman Policy Optimization (BPO), a critic-free method derived from Policy Mirror Descent (PMD). For autoregressive generation with terminal rewards, BPO uses the Bellman equations to reformulate PMD as a trajectory-level objective. The reformulation avoids estimating state values at intermediate states. We prove that it has the same unique optimal solution as the original PMD objective. We derive the practical BPO loss by approximating this objective. Its mismatch-correction weight is a smoothed ratio of complementary token probabilities. Experiments on mathematical reasoning benchmarks demonstrate the effectiveness of BPO.

The forward-reflected-backward (FRB) splitting solves inclusion problems involving the sum of a maximally monotone operator and a monotone Lipschitz continuous operator. Each iteration performs one resolvent step and one evaluation of the Lipschitz operator, plus a reflection term with coefficient one that reuses the previous operator value. We consider the setting where the maximally monotone operator is the subdifferential of a proper, lower semicontinuous, convex function. We identify the tight admissible range of constant reflection coefficients, namely all values larger than one half for a sufficiently small stepsize. We then propose a linesearch-free adaptive variant of FRB for locally Lipschitz continuous operators, in which both the stepsize and the reflection coefficient change across iterations. The stepsize is computed in closed form from a simple local Lipschitz estimate, without requiring the existence or knowledge of a global Lipschitz constant, while maintaining the iteration cost of FRB. The analysis of both methods relies on a new Lyapunov function that combines the distance to a solution, successive differences of iterates, and a gap-type term. The individual terms need not decrease along the iterates, but the stepsize conditions ensure that their weighted combination does. Moreover, alongside nonasymptotic rates, we derive explicit lower bounds on the generated stepsizes, and showcase the effectiveness of the adaptive method through numerical experiments on minimax problems, equilibrium models, and regularized regression.

We study strategic trading around index reconstitution in a continuous-time, multiasset game with transient cross-asset price impact and heterogeneous beliefs about future index membership. Opportunistic traders position before a public announcement, adjust to the revealed composition, and trade around an indexer following a prescribed execution schedule. Under a no-price-manipulation condition, we construct a subgame-perfect Nash equilibrium on every finite horizon. The affine feedback policies are computed from a non-standard matrix Riccati equation and linear ordinary differential equations whose number and dimensions do not grow with the number of traders. Aggregate inventories and price impact depend on beliefs only through the population-average belief, while differences in beliefs affect individual inventory positions. Under mean-field scaling, we construct a mean-field equilibrium and obtain quantitative convergence and approximate-Nash bounds. Numerical illustrations show how competition and impact decay determine the balance between adverse price displacement from anticipatory trading and savings in the indexer's execution costs when opportunists trade against its orders during implementation.