Papers

The facial distance of a polytope, equivalently its pyramidal width, is a geometric condition measure that appears in the convergence bounds for Frank–Wolfe algorithms. Peña and Wirth developed calculus rules for this quantity and conjectured that the ratio of the facial distance to the diameter multiplied by the square root of the dimension is bounded by a constant. We show that this quantity is uniquely maximized by a regular simplex in each dimension, which implies that the conjecture holds and that $\sqrt{2}$ is the smallest constant. We also compute the facial distance of the Birkhoff polytope.

Plug-and-Play (PnP) methods solve imaging inverse problems by incorporating deep denoisers into iterative optimization algorithms. Although practical implementations often decrease the denoiser noise level $\sigma$ along iterations, most existing convergence analyses assume a fixed denoiser. In this work, we establish convergence guarantees for a broad family of Plug-and-Play algorithms with annealed noise level, spanning deterministic methods (RED–GD and PnP–PGD) and stochastic methods (SNORE, equivariant RED, and a variant of PnP–Flow). For each method, we identify an explicit, nonconvex objective associated with the terminal denoising level and prove asymptotic stationarity of the iterates with respect to this objective. Our analysis does not prescribe any decay rate for the noise schedule, and our assumptions cover both learned gradient-step denoisers and exact MMSE denoisers. Overall, our theoretical results bridge the gap between existing PnP convergence theory and the decreasing-denoising practices used by state-of-the-art image restoration methods. We empirically demonstrate the benefits of such schedules and illustrate the predicted convergence behavior on several imaging inverse problems, including inpainting, super-resolution, demosaicing and tomography.

We show that every polytope $P\subseteq[0,1]^n$, and more generally every compact convex set, has Chvátal rank at most $12.22n^2+n\log_2 n+2n+4$. This improves the $O(n^2\log n)$ bound of Eisenbrand and Schulz and, together with the $\Omega(n^2)$ lower bound of Rothvoß and Sanità, shows that the maximum Chvátal rank of a polytope in $[0,1]^n$ is $\Theta(n^2)$. More precisely, if $P$ contains an integer point, then for every $c\in\mathbb{Z}^n\setminus\{0\}$ the inequality $cx\le\max\{cy: y\in P\cap\mathbb{Z}^n\}$ is valid for the $k$-th Chvátal closure of $P$ for some $k\le 12.22n^2+2n+2\log_2\|c\|_\infty+4$. Following Eisenbrand and Schulz, we derive this inequality along a chain of coarser and coarser integer vectors, but instead of halving the vector in each step, we round $\tau c$ for a scale $\tau$ chosen freely in each dyadic window $[2^{-t-1},2^{-t}]$. The main new ingredient is a multiscale version of Dirichlet's approximation theorem, proved by an elementary volume argument: for every $c\in\mathbb{R}^n$, the $\ell_1$-distances of $\tau c$ to $\mathbb{Z}^n$, minimized within each dyadic window and summed over all windows, total less than $1.222n^2$, independently of $\|c\|_\infty$.

We study corank-two distributions $D$ of odd rank $2N+1$ whose associated pencil of skew-symmetric forms $\{d\alpha|_{D}(q)\mid \alpha\in\Omega^1(M),\ \alpha|_{D}=0\}$ lies, at every point, in the generic orbit of the natural $\mathrm{GL}\bigl(D(q)\bigr)$-action; equivalently, $D$ has maximal first Kronecker index. In corank two, this class is the analogue of contact and even-contact distributions. We prove that the corresponding differential relation is ample and hence satisfies the multi-parametric $C^0$-close $h$-principle, and we deduce a complete topological characterization for a $(2N+3)$-dimensional manifold to admit such a distribution. The criterion is satisfied by non-parallelizable manifolds as well. To the best of our knowledge, ampleness has been verified so far for two classes of distributions defined by an open $\mathrm{Diff}$-invariant relation: even-contact distributions (McDuff, 1987) and distributions of rank greater than two with maximal small growth vector, in arbitrary ambient dimension (Martínez-Aguinaga, 2026); for $N\ge2$, our class is a proper open subclass of the latter, singled out by a further condition on the pencil. The verification of ampleness reduces, via Kronecker–Weierstrass normal forms of the restrictions of the pencil to hyperplanes, to computing the convex hulls of the connected components of the complement of the resultant hypersurface in the space of pairs of real binary forms of fixed degrees, components separated by a winding number, i.e., by the classical Cauchy index. This is in contrast with the even-contact case of McDuff, where the set to be removed is thin (i.e., of codimension at least two), so that its complement is automatically connected and ampleness is immediate.

Deploying AI systems requires deciding when to delegate tasks and when humans should intervene to monitor and mitigate risk induced by AI operations. These decisions are challenging when failures cluster: a hallucination or harmful output can trigger further errors, creating periods of elevated risk. We introduce a continuous-time framework for learning adaptive human oversight under such turbulent AI risk. Existing oversight and delegation formulations condition on history but do not model incident clustering, or its suppression by supervision effort, jointly with the delegation decision and this study fixes this gap. The self-exciting dynamics capture how risk events increase the likelihood of subsequent events, making their timing and history central to decision-making. We formulate a stochastic control problem combining human actions, monitoring effort, and switching between human-AI-assisted operation and full AI delegation, balancing operational rewards against oversight costs, and cascading AI-failures and induced uncertainty. Human participation is an endogenous component of risk management: the policy determines both when oversight is needed and how much effort to allocate. We study a relaxed switching formulation and propose Hawkes-PPO, a policy-gradient method that uses a bank of exponential filters of observed incident times. In a synthetic environment it attains a higher risk-adjusted objective than either fixed regime and approaches an approximate full-information oracle. We illustrate our results with numerical simulations by examining how cascade risks influence intervention and delegation, connecting reinforcement learning with adaptive human oversight of AI systems. In particular, we illustrate the benefit of our switching strategy and Hawkes-PPO algorithm to monitor the project efficiently along time, reducing turbulent risks occurrences and costs.

This paper presents a constrained iterative Linear Quadratic Regulator (Box-iLQR) method for nonlinear aerospace optimal-control problems with bounded states and controls. Logarithmic barrier functions are used with a two-fold objective: the barrier parameter is reduced toward zero to recover the constrained optimal trajectory, while a finite barrier is retained to construct a smooth, constraint-aware feedback policy for closed-loop operation. The method is demonstrated on spacecraft attitude control, a minimum-fuel planar orbit transfer, and an L2 to L2 halo-orbit transfer. Box-iLQR respects the imposed constraints and produces nominal trajectories together with time-varying feedback policies. Closed-loop simulations under disturbance, measurement noise, and control uncertainty show substantial reductions in terminal error. These results demonstrate that Box-iLQR provides both high-accuracy constrained trajectory optimization and practical constraint-aware feedback within a single framework.

The classical method of Kaczmarz, introduced in 1937, is a textbook iterative method for solving linear systems $A x = b.$ Despite its widespread use, particularly in the context of solving inverse problems where it is commonly included in software packages within Matlab, Python, and Julia, the precise characterization of convergence has long been deemed difficult to obtain. While different bounds on convergence rates have been established, they are largely unsatisfying as they cannot explain the classical, cyclic-update method's efficient convergence in practice. In this work, we obtain a tight characterization of convergence of the classical Kaczmarz method, establishing both linear and sublinear convergence bounds. The worst-case tight (i.e., exactly attained by some instances in the considered family) bounds are expressed in terms of a fixed matrix that depends only on $A,$ but are not fully interpretable in terms of the matrix spectrum and row correlations, which had been observed to have an impact on convergence. We thus provide relaxations of these bounds that are fully expressible in terms of matrix row norms, row correlations, rank, and extremal positive singular values. The provided relaxed bounds explain one-cycle convergence in special cases where the matrix rows are all either parallel or orthogonal to each other. We further argue that the dependence on different parameters appearing in the bounds is necessary and within a small constant factor of the best attainable in the worst case. Finally, our bounds explain why the classical cyclic update is faster than the randomized one when matrix rows are weakly correlated, which is often observed in inverse problems where the cyclic method is used.

We consider the Lorentz positive cone of $n\times m$ matrices that map the Lorentz cone in $\Rbb^m$ into the Lorentz cone in $\Rbb^n$. The Lorentz positive cone and its dual, the cone of Lorentz separable matrices, are shown to provide polynomial-time algorithms for the problem of minimizing a bilinear objective over variables contained in ellipsoids in $\Rbb^n$ and $\Rbb^m$. We also demonstrate how these cones can be used to strengthen SDP relaxations of other nonconvex quadratic optimization problems, including the two-trust-region subproblem.

We study finite element approximations for an optimal control problem with pointwise tracking governed by a semilinear elliptic PDE. The control variable enters the state equation as a reaction coefficient, resulting in a bilinear control, and the cost functional includes point evaluations of the state variable. We propose two discretization schemes: a semidiscrete scheme, where the control variable is not discretized, and a fully discrete scheme, where the control variable is discretized using piecewise constant functions. In two and three dimensions, for both schemes, we establish convergence of discrete solutions and prove a priori error bounds that behave as $\mathcal{O}(h|\log h|)$ in the $L^2(\Omega)$-norm for approximating a locally optimal control. In two dimensions, these bounds are improved to $\mathcal{O}(h)$ for the fully discrete scheme and $\mathcal{O}(h^{2}|\log h|^2)$ for the semidiscrete scheme. We conclude with numerical experiments that demonstrate the performance of the proposed schemes.

The phenomenon of the top subspace stabilization of the Hessian matrix is an surprising and critical aspect in study of the second-order information of neural network training. Prior work argues that the top subspace of the Hessian stabilizes by measuring the overlap between the top subspaces of the step-wise Hessian, and explains this stabilization with diminishing parameter change in the late phase of training. In this paper, we define a new instability metric for the subspace evolution, and use it to detect subspace stabilization that is independent of the magnitude of parameter change. In the meantime, we observe that the gradient covariance matrix has a similar property of its top subspace to the Hessian. By using a between-class and within-class decomposition of the gradient covariance matrix, we identify an explicit form that gives a near-perfect approximation of the top-$(C-1)$ subspace of the Hessian and the gradient covariance matrix. In the gradient flow set-up, we show that the slow evolution of the idenfied approximation is due to the separation between the outlier and the bulk eigenvalues of the Hessian matrix, thus providing an explanation to the phenomenon of the top subspace stabilization of the Hessian matrix.

Multi-objective Bayesian optimization (MOBO) is commonly approached through specialized acquisition functions or scalarization schemes designed to explicitly account for trade-offs among non-preferential objectives. In this work, we show that such complexity might be unnecessary. We propose a framework that extends standard single-objective acquisition functions directly to the multi-objective setting through a hypervolume-based transformation. We further extend hedge strategies for acquisition functions, which are typically used only in single-objective optimization, to the multi-objective regime. Our approach requires minimal modification to existing Bayesian optimization pipelines and avoids the need for bespoke multi-objective formulations. We demonstrate how a broad class of commonly used single-objective acquisition functions and hedge strategies can be adapted in a principled manner to handle multiple objectives, while preserving their intuitive interpretation and computational efficiency. Empirically, we evaluate the proposed methods across a range of synthetic and real-world multi-objective benchmarks. Despite their simplicity, our extensions consistently match or outperform more complex state-of-the-art MOBO methods in terms of optimization performance and sample efficiency. These results suggest that effective multi-objective Bayesian optimization can be achieved by reusing and carefully extending well-established single-objective acquisition strategies, offering a simpler and more flexible alternative to existing approaches.

The Gromov-Wasserstein (GW) distance measures the discrepancy between metric measure (mm) spaces and identifies optimal alignments between them based solely on their intrinsic structure. Since it identifies isomorphic mm spaces, it provides a natural notion of distance for heterogeneous datasets which may admit isomorphic representations. In order to accelerate computation of GW distances, many practitioners employ entropic regularization to obtain an Entropic GW (EGW) problem. The most popular EGW solver is the Mirror Descent (MD) algorithm, which reduces EGW computations to an iterative process where an entropic optimal transport (EOT) problem is solved at each iteration. Despite its widespread use, the convergence of MD for this problem has only been established for restricted classes of costs. On the other hand, a recently proposed dual gradient method is available for general costs, but requires a choice of step size which depends on the regularization parameter. To address these two issues, we introduce Averaged Mirror Descent (AMD), which averages consecutive MD steps, and prove its convergence for arbitrary costs. Then, we establish that the dual gradient method with a fixed step size also converges for arbitrary costs at the cost of a more complicated iteration. In both cases, we also account for inexact iterations which are inescapable in practice. We compare the empirical performance of these methods across various settings and, in particular, show that AMD and the dual gradient method both converge on an example where classical MD fails.

Chambolle-Dossal acceleration guarantees $F(x_n)-F^*=o(n^{-2})$ for every fixed smooth convex loss with a minimizer. We show that this qualitative improvement admits no universal quantitative rate. For every damping parameter $\alpha>3$ and positive nondecreasing gain $G(n)\to\infty$, we construct a fixed one-dimensional smooth convex loss whose exact CD orbit satisfies $\sup_{n\ge1} n^2G(n)\bigl(F(x_n)-F^*\bigr)=\infty$. Thus no divergent gain improves the $n^{-2}$ scale for all fixed losses, even with instance-dependent constants. The construction prescribes queried gradients and realizes infinitely many slow blocks within one smooth convex objective. Under local $p$-power growth with $p>2$ and sufficiently strong damping, we also construct a fixed loss whose exact CD orbit satisfies $F(x_n)-F^*\sim Dn^{-2p/(p-2)}$, $D>0$, establishing the sharpness of the known convergence rate. Both main results are formally verified in Lean 4.

On a network with a fixed set of verified sources, discounted averaging induces an equilibrium support $h_i^S=\mathbb{E}_i[\lambda^{T_S}]$, the discounted probability that a random walk reaches $S$ before attenuation. We define the discounted hitting domination number $\delta_{\lambda,\tau}(G)$ as the minimum number of sources required to guarantee $h_i^S\ge\tau$ at every vertex. Although this potential is known through penalized and group hitting probabilities, the associated minimum-cardinality uniform-coverage problem appears to be new. Aggregate support is monotone submodular, while the uniform-floor problem is an exact submodular-cover problem. Moreover, if $\lambda^{r+1}<\tau\le\left(\frac{\lambda}{\Delta}\right)^r$, then $\delta_{\lambda,\tau}(G)$ equals the distance-$r$ domination number. This yields NP-completeness and APX-completeness at $(\lambda,\tau)=(1/4,1/14)$ on graphs of maximum degree three. For spiders, we obtain an exact finite-state characterization and a polynomial-time algorithm for every fixed rational pair $(\lambda,\tau)$, and show that the branching vertex need not belong to a minimum source set. Finally, an exact mixed-integer linear formulation certifies optimal placements on a real network and a synthetic graph and demonstrates substantial differences from degree, closeness, and classical domination.

We investigate a class of backward stochastic differential equations (BSDEs) with at most quadratic growth which explode at a possibly unbounded random horizon, defined through the first hitting of zero of an adapted Ito process. In contrast with the classical theory of singular BSDEs, the explosion is generated by the nonlinear dependence of the generator on the martingale integrand, rather than by a superlinear coercivity condition in the solution component. We construct a minimal singular solution and derive two-sided estimates on its explosion, together with weighted BMO estimates for the martingale integrand. We also obtain the uniqueness and exact explosion rates under additional structural assumptions. For Hamiltonian generators, we establish a verification theorem for an infinite-horizon stochastic optimal control problem with possible non-Markovian state constraints, showing that the BSDE feedback induces the unique optimal constrained law. We finally specialize the theory to exit times of uniformly elliptic Markov diffusions and recover the connection with large solutions of viscous Hamilton-Jacobi equations.

Outliers are essential for evaluating and improving the robustness of machine learning systems, especially when future distributions may differ significantly from historical training data. In high-stakes applications, robustness often depends on rare cases that finite datasets fail to capture, making simple resampling or perturbation insufficient for stress scenario generation. Existing outlier synthesis methods typically rely on sparse neighborhoods, low support latent regions, or classifier boundary crossings, which can be heuristic, unstable, and tied to specific modalities or architectures. We therefore propose Sinkhorn Boundary Outlier Generation (SBOG), a structured framework for latent-space outlier generation that couples Sinkhorn optimal transport geometry with distributionally robust boundary modeling. The resulting Sinkhorn-induced support cost guides the sampler toward weakly supported boundary regions, while semantic constraints prevent uncontrolled drift from the intended context, yielding controlled deviations from the in-distribution reference measure rather than arbitrary sparse-region samples. Experiments on time series anomaly generation and image outlier synthesis show that our framework produces informative, semantically controlled outliers and improves downstream robustness evaluation across modalities, providing a foundation for stress scenario generation beyond empirical support.

Manufacturing firms face uncertainty in both supply and demand. They implement process flexibility, enabling plants to produce multiple products, to hedge against supply–demand mismatches. Assigning more products to a plant results in efficiency losses at the plant, whereas limited flexibility may improve efficiency because of the specialization effect. Furthermore, the flexibility design determines the sourcing locations of products, thereby affecting demand due to the country-of-origin effect and exposure to tariffs. We incorporate these mechanisms into the supply network flexibility design problem by modeling supply and demand uncertainties as endogenous, that is, as functions of the first-stage design decisions. We propose a two-stage stochastic program and solve it via a decomposition scheme that relies on distribution-specific optimality cuts. To avoid exhaustive enumeration of all distributions, we derive distribution-free formulations that are valid for all distributions and for subsets of supply and demand distributions. Through extensive computational experiments, we demonstrate the superior performance of our solution approach and distribution-free formulations relative to exhaustive distribution enumeration and state-of-the-art inequalities, respectively. Moreover, we explain how the problem structure affects the effectiveness of the proposed formulations. From a managerial perspective, we show the trade-offs between endogenous effects and how they affect the optimal flexibility design.

Optimization of mathematical models involving integer decisions and neural networks with ReLU activation (ReLU ANNs) is a challenging task. Nevertheless, such models are an enabling technology in many application domains. A prominent example is superstructure optimization in chemical engineering, where ReLU ANNs are frequently employed as surrogate models for complex nonlinear processes. We survey recent developments in this area. We argue that in addition to network size and training options of the ANNs, the ReLU activation geometry and the number of linear regions on the domain of interest have a strong impact on computational optimization performance. While standard model compression approaches such as structured pruning reduce network size, they do not explicitly address geometric considerations. Therefore, we propose a novel region-reduced model compression approach that combines the stabilization of unstable neurons and the merging of redundant neurons to reduce the number of linear regions while maintaining predictive accuracy through error compensation. We evaluate our method against standard compression approaches on multiple optimization use cases. First, the two-dimensional peaks function for which we can visualize the activation geometry. Second, on optimization over individual surrogate ReLU ANNs for three chemical processes, and third, on a hybrid superstructure optimization problem that involves the three ReLU ANNs, additional process submodels, and binary variables. The results for the superstructure problem demonstrate the large potential of region reduction with a decrease of 40% to 50% in computational time, yielding solutions closer than 1% to the reference at negligible effort of obtaining the compressed model.

We consider the problem of controlling a three-dimensional falling liquid film using same-fluid blowing and suction through the substrate. Building upon recent results for the corresponding two-dimensional problem, we employ a reduced-dimensional weighted-residual model for the film height and downstream flux as the basis for controller design. Linear stability analysis is used to estimate the number of actuators required to stabilise the film, while a linear-quadratic regulator (LQR) framework is used to construct feedback controls for the linearised system. The resulting gain matrix is then combined with observations from direct numerical simulation (DNS) of the three-dimensional Navier-Stokes equations, allowing controls derived from the reduced-dimensional model to be applied directly to the full flow. Through systematic numerical tests, we demonstrate that these controls successfully stabilise the flat-film solution across a physically relevant range of Reynolds numbers and domain configurations. We further show that both the number and placement of actuators play a crucial role in control performance, with spanwise instabilities requiring additional actuator coverage in the cross-stream direction. Although the extension from two to three dimensions substantially increases the cost of controller construction, this cost is incurred entirely offline during the computation of the gain matrix, making the overall framework computationally tractable. Our results demonstrate that LQR feedback controls derived from a three-dimensional weighted-residual model retain predictive value when applied to fully resolved DNS, extending previous two-dimensional control methodologies and helping bridge the gap between theoretical control design and experimentally relevant falling-film flows.

Distributionally robust optimization (DRO) provides a principled framework for learning under distribution shift, but its practical use is hindered by the difficulty of evaluating worst-case risks for nonconvex loss functions. We study a penalized DRO formulation in which the adversary may choose any distribution but incurs a Wasserstein penalty for deviating from the empirical distribution. We show that the adversary's problem can be reformulated as an optimization problem over transport maps that push empirical samples to adversarial ones, and we prove that optimal maps are cyclically monotone. We also show that standard adversarial training—based on per-sample local optimization—violates cyclical monotonicity and wastes transport costs unless the adversary is severely restricted. We propose two remedies. First, we introduce multi-start particle ascent, which alternates parallel gradient ascent with reassignment to enforce cyclical monotonicity across samples. Second, we parameterize adversarial maps as gradients of input-convex neural networks, which guarantees cyclical monotonicity by construction. Experiments on robust regression, image classification, and robust control show that our methods consistently outperform standard adversarial training and state-of-the-art baselines, achieving improved robustness and better generalization under distribution shift.

Gradient-based optimization (GBO) of problems constrained by partial differential equations is central to parameter identification and structural design in solid mechanics. Lattice Boltzmann method (LBM) schemes for linear elastic solids have only recently been derived. Combined with automatically generated discrete adjoint collision kernels, they make GBO possible with a purely LBM-based primal and adjoint solver. We propose, to our knowledge, the first optimization benchmarks for solid LBM: three inverse problems, each with a manufactured solution of the forward problem and a known reference control. These are the identification of two scalar amplitudes of a body force field on a periodic domain, the distributed control of the same field with up to $12800$ degrees of freedom, and the identification of Young's modulus on an elliptic plate with Dirichlet and Neumann boundaries. The primal scheme converges with second order on the periodic domain and with first order on the elliptic plate. The recovered controls follow these rates, with experimental orders of convergence of $1.92$ and $1.06$, respectively. Gradients obtained via forward-mode automatic differentiation (f-AD) and the discrete adjoint LBM (ALBM) require the same number of optimization steps and yield control errors agreeing to at least eleven significant digits. The ALBM gradient is further validated node-wise against a finite difference quotient, with deviations decreasing under mesh refinement. The benchmarks are implemented in the open-source LBM framework OpenLB and provide a verification basis for optimization in LBM-based solid mechanics and its extension to fluid-structure interaction.

In this paper we prove a version of Shnirelman's inequality in the two-dimensional square $[0,1]^2$. More precisely, given a target fluid configuration $f$ (a volume-preserving diffeomorphism) we construct a divergence-free velocity field $v\in L^q_tL^p_x$ whose flow connects $f$ to the identity and whose $L^q_{t} L^p_x$-norm can effectively be bounded by the $L^p$-difference of $f$ and $id$. However, higher regularity of the vector field might be affected. The main difference with the higher dimensional case is that, because of the topological obstructions of dimension $\nu=2$, we have to allow that different fluid trajectories 'intersect' in space. We also observe that this choice leads to the emergence of irreversible behaviors in the Eulerian dynamics.

Convergence in parameter estimation classically requires “persistency of excitation," which is a non-degeneracy condition on a trajectory-dependent signal. This paper develops, to the best of our knowledge, the first such excitation theory for a nonlinear partial differential equation. In the context of identifying the unknown viscosity from spectral observations in the two-dimensional Navier–Stokes equations for incompressible fluids, our excitation condition is computable, verifiable a priori, and sharp with respect to scaling. Our approach employs a data assimilation methodology to account for an unknown initial state and select candidate viscosities by minimizing the loss between the low-mode observations of the fluid velocity and low-mode projection of a nudging-based filter that assimilates these observations. The main novelty of our framework is to distinguish a new functional, $\mathsf{W}$, representing the work done by the filter's associated sensitivity variable on the enstrophy, around which our entire analysis is centered. We show that $\mathsf{W}$ is explicitly comparable to the observability Gramian of the associated Gauss–Newton iteration, and subsequently establish the following dichotomy at every critical point of the observational loss: either $\mathsf{W}$ exceeds a certain threshold, in which case the candidate viscosity obeys an explicit error estimate that depends inversely on $\mathsf{W}$ and the observational density $N$, but directly on the error between initial conditions, or else $\mathsf{W}$ is below the threshold and the observations are quantitatively insensitive to parameter updates in a way that is detectable to the user. Notably, our approach is energy-based, and therefore expected to be adaptable to many other nonlinear dissipative systems.

For linear-quadratic $N$-dimensional Nash systems with mean-field-like scaling and strongly displacement monotone data, we prove uniform-in-time a priori derivative estimates exhibiting the expected scaling as $N \to \infty$ in order to pass the system to the limit. As a consequence, in a heterogeneous (i.e., non-symmetric) mean-field setting, we obtain both uniform-in-time convergence of the Nash system on any arbitrarily long horizon $[0,T]$ (as $N \to \infty$), and convergence of both the $N$-dimensional Nash system and the corresponding master equation to their respective ergodic counterparts (as $T \to \infty$).

We study fundamental information-theoretic limits of robust stochastic optimization when the distribution is known only through its exact moment sequence. We develop a unified framework that produces families of distinct distributions sharing all moments yet inducing radically different optimal decisions, thereby establishing strong impossibility results for a range of decision problems under moment ambiguity. Our approach gives two explicit constructions: a binary and an $N$-way construction showing that distributions with identical moment sequences can nevertheless exhibit arbitrarily different quantiles, order-statistics and threshold regions, forcing incompatible optimal actions. These families of distributions yield, in fact, strong impossibility results across several stochastic optimization problems. First, for the newsvendor problem, moment equivalence causes quantile ambiguity, inducing any fixed or randomized order quantity to fail arbitrarily badly. Second, for revenue maximization, no deterministic or randomized posted pricing scheme can secure a nontrivial approximation relative to the full-information benchmark. Third, for the secretary with cardinal observations setting, the worst-case robust value over all exact moment disclosures is exactly the classical $1/e$ finite-horizon value as opposed to the celebrated result of $0.58$ success probability with full-information by Gilbert and Mosteller (J. Am. Stat. Assoc., 1966). We also recover and expand upon the recent impossibility result of Correa et al. (STOC, 2026) for prophet inequalities with moment knowledge. Indeed, exact moment knowledge can yield at best a $\Theta(1/\log n)$ competitive ratio, even when competing against relaxed benchmarks based on expected $r$-th order statistics or when the algorithm is allowed to select $r$ items.

This paper studies constrained optimization problems through the lens of feedback control. Building on the interpretation of Lagrange multipliers as feedback controllers, we propose the saddle-point PID (SPPID) dynamics: a unified proportional–integral–derivative (PID) framework for continuous-time saddle-point dynamics. The proposed dynamics employ PID control for equality constraints and anti-windup PI control for inequality constraints. We show that SPPID is equivalent to a preconditioned primal–dual gradient flow of the augmented Lagrangian. This equivalence reveals the distinct role of each feedback component: integral action enforces constraint satisfaction, proportional action induces the augmented Lagrangian structure, and derivative action modifies the geometry of the primal dynamics via a state-dependent Riemannian metric. For convex problems, we establish convergence to the KKT set together with conditions for global asymptotic stability. For equality-constrained problems, we establish local exponential convergence for nonlinear equality constraints under standard assumptions, and show that projected gradient flow emerges as the infinite derivative-gain limit. For strongly convex problems with affine equality or inequality constraints, we establish global exponential convergence by leveraging contraction theory. Finally, we provide various numerical examples to illustrate the utility of SPPID.

The empirical success of pretraining large language models has inspired a deeper investigation into the underlying loss landscapes and the optimization dynamics. Recent empirical and theoretical study suggest that the training loss landscape often exhibits a "river-valley" structure, which features a low-loss manifold (river) flanked by sharp orthogonal directions with higher loss (mountains). In the long term, the optimization progress is determined primarily by the progress along the river. Within such a landscape, gradient descent with large learning rates can move faster along the river despite high apparent loss due to vertical oscillations, while a subsequent sharp decay in the learning rate suppresses these oscillations, revealing genuine optimization progress. This explains the recent success of warmup-stable-decay (WSD) learning rate scheduler which, unlike cosine scheduling, keeps stable high learning rate and decays before producing intermediate checkpoints. Building on this foundation, in this work we take a step further and study the role of momentum within such a loss landscape. We establish theoretical analysis that characterizes how momentum accelerates optimization by stabilizing large learning rates that can not be tolerated by vanilla GD without deviating significantly from the river. The enabled large learning rate in-turn gives greater speed along the river and makes faster essential progress in the long run. Another intriguing observation from theory is that for a river-valley landscape with very flat and slow-spinning river, the momentum itself does not contribute directly to acceleration in terms of the speed of tracking the river, while the main acceleration comes from the admissible larger learning rate.

Quantum-inspired evolutionary optimization (QIEO) represents design variables as a set of qubits and searches a continuous, multi-dimensional landscape through rotation of the qubit's amplitude pair. Every generation rotates those amplitudes toward a single elite, which corresponds to that generation's best. The update is cheap, almost parameter-free, and well-suited for massive parallel implementation, which has encouraged its adoption in engineering, design, and planning applications. However, there are critical issues with this formulation, principally, the treatment of design variables as independent probability components which make it incapable of exploiting local curvature, anisotropy, or variable coupling. Despite this, QIEO is believed to hold promise, and has been used extensively to solve real-world problems, with significant qualitative and computational advantage over its classical counterpart, Genetic Algorithm (GA). A collection of 256 (actually 508; 256 unshifted + 252 shifted, 4 could not be shifted) continuous function are selected from the prior works, in such a way that they represent eleven landscape characteristics, namely continuity, differentiability, separability, scalability, modality, convexity, conditioning, symmetry, maximum dimensionality, dimension dependency, and the coupling pattern of the design variables. These functions are then solved by three QIEO variants, two GA encodings and Hansen's Covariance Matrix Adaptation Evolution Strategy (CMA-ES). The results are evaluated in terms of computational cost, solution precision, and specialization across landscape characteristics. They identify the conditions under which QIEO provides competitive performance, clarify where its independent-variable representation becomes limiting, and establish whether particular QIEO variants offer advantages for specific landscape characteristics.

Quantum-inspired algorithms emulate quantum mechanical principles, such as, superposition, interference, and probabilistic amplitude evolution, on classical hardware by representing candidate solutions as qubit vectors and evolving them through rotation-gate operators. This approach offers higher optimization performance without physical qubits, and has been shown to achieve order-of-magnitude speedups (10–80$\times$) over traditional solvers on combinatorial, high-dimensional NP-hard problems. A critical barrier to adoption, however, is the lack of a unified execution framework that delivers both algorithmic performance and hardware portability. We present Cross-Backend Quantum Inspired Evolutionary Optimizer (QIEO), the runtime core of BQP's BQPhy solver, which addresses this gap through a single-source-of-truth architecture. One C++ implementation of the QIEO algorithm is compiled once per hardware target and exposed to multiple high-level languages via thin binding layers. The framework dispatches to CPU (sequential), OpenMP~5 (multi-core), CUDA (NVIDIA), and HIP (AMD) backends at runtime, adapting kernels to each device's memory hierarchy and warp/wavefront execution model. The framework's real-world utility is validated through binding demonstrations that share the identical C++ runtime. BQPhy's Python library is demonstrated on a neural network hyperparameter optimisation achieving 88.60% test accuracy on MNIST. BQPhy's MATLAB's Toolkit is tested on wind farm layout optimisation attaining $365\,399 \pm 4\,552$~MWh/yr, which is statistically indistinguishable from particle swarm optimisation and $+7.6\%$ above genetic algorithms on a 32-variable constrained engineering problem. The Julia package tackles the Lotka–Volterra parameter estimation where BQPhy replaces native Julia solvers on the same residual, cutting mean SSE by $2.1\times$.

In this paper, we consider bilevel variational inequality problems, where the feasible set is the solution set of another variational inequality. We propose a Bregman inertial regularized extragradient method for solving these problems. By performing the inertial extrapolation in the dual space, the proposed method aligns the inertial step with the structure of the Bregman three-point identity. We also derive explicit non-asymptotic bounds for the generalized gap function. More precisely, using diminishing regularization, we establish convergence rates of $\mathcal{O}(1/k^{1-b})$ for the optimality gap and $\mathcal{O}(1/k^b)$ for the feasibility gap with $ 0 < b< 1$. For constant regularization parameter $\eta>0$, we establish an $\mathcal{O}(1/k)$ rate for the optimality gap and an $\mathcal{O}(1/k)+\mathcal{O}(\eta)$ bound for the feasibility gap. Furthermore, numerical experiments on traffic networks and multi-portfolio Nash equilibrium problems demonstrate the practical effectiveness of the proposed framework.