Papers

The forward-reflected-backward splitting of Malitsky and Tam converges weakly for step sizes $\lambda\in(0,\tfrac{1}{2L})$, where $B$ is monotone and $L$-Lipschitz, and the question whether this bound is tight has been recorded as open \cite{GS}. We show that it is: for the matched-skew instance $A=LJ$, $B=LJ$ ($J$ the counterclockwise rotation by $\pi/2$ in $\R^2$), the iterates fail to converge for every $\lambda\ge\tfrac{1}{2L}$ and diverge for $\lambda>\tfrac{1}{2L}$. More generally, for the skew–rotation family $A=\gamma J$, $B=J$, the exact stability threshold is $\lambda^\star(\gamma)=1/\sqrt{(1+\gamma)(3-\gamma)}$, which attains its minimum $1/2$ at the matched skew $\gamma=1$ and recovers the reflected gradient constant $1/\sqrt3$ at $\gamma=0$. The same instance is a counterexample for the reflected–forward–backward method of Cevher and Vũ, which coincides with forward-reflected-backward on linear operators.

We analyze a one-evaluation-per-iteration method for $\mathrm{VI}(C,B)$ with monotone $L$-Lipschitz $B$: $x_{k+1}=P_C(x_k-\lambda B(u_k))$, $u_k=x_k+\theta_k(x_k-x_{k-1})+\beta_k(x_k-u_{k-1})$, $\theta_k+\beta_k=1$. For constant step and summable filter ($\sum_k\beta_k<\infty$) we prove weak convergence via a Lyapunov function with exact rational dissipation budgets, including Malitsky's reflected gradient method. It is robust to summable operator errors; if $C$ is bounded, the constant-step range reaches $\lambda<(\sqrt2-1)/L$. For unbounded $C$ this reduces to an a priori boundedness statement, certified to $\lambda L=0.387$ by dissipation trading; for affine $B$ on polyhedral $C$, $\dist(x_n,S)\to0$ with $\sum_n\dist^2(x_n,S)<\infty$, strong convergence for bounded $C$, and $R$-linear rates after face identification. The main result is a safeguarded adaptive step-size rule needing no $L$ and no extra evaluations, proved weakly convergent unconditionally: a data-driven Lyapunov weight removes $L$ from the dissipation budgets. For affine $B$ with $C=\mathcal H$, sharpness of $\lambda L=1/\sqrt3$ via a rotation lower bound; the same constant is sharp for unconstrained nonlinear $B$, with convergence for square-summable operator values. Unconditional convergence below $1/(\sqrt3\,L)$ reduces to a marginal-pole absolute-stability statement; the projected case remains open. Under strong monotonicity we prove $R$-linear convergence with explicit contraction, and certify $\Omega(L)$ speedups over constant steps. Numerics confirm the gains.

The classical perceptron cycling theorem of Block and Levin \cite{BlockLevin1970} bounds correction sequences whose selected updates come from a finite set and have nonpositive inner product with the current state. We prove a robust variant for additive trajectories $z_{k+1}=z_k+u_k$, for all integers $k\geq0$, with increments in a finite set $U\subset E$, where $E$ is a finite-dimensional real inner-product space and $0\in\conv U$. Let $A\colon E\to E$ have positive-definite symmetric part, without requiring symmetry, and let $B\geq0$. If each $u_k$ is a $B$-approximate minimizer of $u\mapsto\ip{Az_k}{u}$ over $U$, then $\sup_{k\geq0}\norm{z_k}\leq C(1+\norm{z_0}+B)$, with $C=C(E,U,A)$ independent of the initial state, $B$, and all admissible update choices. We derive two algorithmic consequences. The first is an $O(k^{-1})$ last-iterate norm bound for harmonic vertex-returning Frank–Wolfe for affine strongly monotone variational inequalities on polytopes with relatively interior solutions. It extends the quadratic Frank–Wolfe/herding guarantee of Bach, Lacoste-Julien, and Obozinski \cite[Section~4.2]{BachEtAl2012} to nonsymmetric affine operators. The second consequence concerns oblique relaxation for linear inequalities. Greedy corrections through a fixed matrix with positive-definite symmetric part remain bounded even for inconsistent systems. In particular, for a square matrix $G$ with positive-definite symmetric part, repeatedly increasing the coordinate corresponding to a most-violated inequality of $Gx\geq b$ terminates at an exactly feasible point after finitely many unit corrections. This remains true under bounded additive selection errors, provided the stopping test uses the true inequalities.

This paper focuses on a class of nonconvex nonsmooth multicomposite optimization problems for training RNNs (Recurrent Neural Networks). We first establish easily verifiable conditions under which an approximate first-order d-stationary point of the problem is guaranteed to be an approximate second-order d-stationary point. Subsequently, we propose a sequential regularized piecewise affine algorithm, SRPA, that minimizes a regularized piecewise affine function at each iteration, where the idea of active-set strategy is incorporated to reduce the computational cost per iteration. Leveraging the aforementioned conditions for second-order d-stationarity, we establish the global convergence of SRPA to second-order d-stationary points and the complexity bound of $ \mathcal{O}(\epsilon^{-2}) $ for obtaining an $ \epsilon $-approximate second-order d-stationary point. Finally, numerical experiments for training RNNs on synthetic and real-world datasets demonstrate the promising performance of SRPA compared with state-of-the-art algorithms.

Loss spikes are recurrent instabilities in neural-network training and can arise from multiple mechanisms. For Adam in particular, macroscopic loss spikes have been linked to optimizer dynamics, yet how its two momentum timescales govern them remains unclear. We investigate this dependence by mapping training dynamics across the $(\beta_1,\beta_2)$ plane. Across a range of model–task settings, an approximately linear boundary, $1-\beta_2=C(1-\beta_1)$, separates spiky from non-spiky dynamics, whereas a one-dimensional quadratic loss produces approximately cubic slope. A one-dimensional superquadratic loss $L(x)\propto|x|^n$ recovers the near-linear scaling and links the boundary coefficient to the effective loss exponent $n$. We further show that confident cross-entropy losses develop a core–wall landscape comprising a narrow quadratic core followed by a steep wall, which produces effective superquadratic behavior at the scale of an optimizer update. Together, these results connect Adam loss spikes to both the mismatch between momentum timescales and finite-scale superquadratic loss geometry beyond the Hessian.

This paper studies predictive online control for linear dynamical systems with time-varying cost functions, motivated by applications such as real-time adaptive control of motorized LiDAR sensing systems, where the controller must adjust the sensing direction online to balance localization accuracy, scanning efficiency, and smooth actuation under changing environmental conditions. We propose a predictive online control (POC) algorithm that leverages short-term predictions of future cost functions while accounting for the memory effect induced by system dynamics. Using the disturbance-action controller (DAC) parameterization, the online control problem is transformed into an OCO-with-memory formulation over policy parameters. We develop a windowed receding-horizon update that incorporates short-term predictions and accommodates a total-variation regularizer to suppress abrupt policy variations. Theoretically, we prove that POC achieves a dynamic policy regret bound scaling with the path length of the comparator sequence, and provide sufficient conditions in terms of logarithmic prediction and memory horizons for the regret guarantee. The proposed method is evaluated on a motorized LiDAR sensing task, demonstrating improved localization accuracy together with a favorable trade-off between sensing accuracy and scanning completeness.

We perform a detailed study of a principal–agent problem in a continuous time version of the celebrated Holmström–Milgrom model (Econometrica 55 (2), 1987) where we add limited liability for the Agent. We develop a probabilistic methodology to prove that the Principal's value function is the unique bounded classical solution to a fully nonlinear and fully degenerate partial differential equation (PDE) with Cauchy-Dirichlet boundary conditions on $[0,T]\times[0,\infty)$. Indeed, we also prove infinite continuous differentiability of the solution in the interior of the domain. The strength of our regularity result is such that we can ensure existence of optimal controls in strong form—a rare occurrence in dynamic contracting—and we obtain fine properties of the optimal control map, including a characterisation via a further nonlinear degenerate PDE.

We give exact polynomial-time algorithms in fixed dimension for mixed-integer quadratic programming over arbitrary rational polyhedra and for integer cubic programming over bounded rational polyhedra. The quadratic objective may be indefinite; the algorithm detects infeasibility, certifies unboundedness, or returns an exact minimizer. For an instance with N variables and m constraints, the running time is 2O(Nlog N) times (m+ 1)O(N) times a factor in which the binary encoding lengths of the coefficients are raised to a power O(N); for pure integer problems, and for mixed-integer problems with integral right-hand side and linear term, only the constraint and quadratic matrices enter this factor, while the right-hand side and linear term enter polynomially. This replaces the dependence on coefficient magnitudes in earlier fixed-parameter results by dependence on encoding length, and the exponent O(N) is optimal under the exponential time hypothesis. For cubic objectives, integer minimization over a bounded polyhedron was previously known to be polynomial only in dimension two. Our algorithm extends to any polynomial whose second directional derivatives are concave in position, which includes cubics with an added concave quartic form. Both algorithms rest on a hidden convexity of such objectives on bounded lattice sets, exposed by a cover of the feasible integer points by reflecting cells. For quadratics, a negative-curvature integer displacement excludes a cell, and otherwise supporting inequalities on its integer points permit exact optimization through an integer-query separation oracle. For cubics and the broader class, negative curvature at a query point yields a linear cut valid for every global minimizer in the cell, and supporting inequalities among the remaining admissible points play the same role.

We study the deterministic oracle complexity of smooth convex optimization when the algorithm receives only exact function values. The objective is a globally $\beta$-smooth convex function, all queries and the final output are restricted to the Euclidean ball of radius $R$, and the unique minimizer lies in the ball of radius $R/2$. We establish an upper bound of $O(d\sqrt{\beta R^2/\epsilon})$ using coordinate finite differences together with an error-robust accelerated projected method. Our main contribution is a matching lower bound, up to the high-accuracy saturation of the construction: any deterministic adaptive value-oracle algorithm requires $\Omega\!\left(d\min\{\sqrt{\beta R^2/\epsilon},(d/\log(ed))^{1/3}\}\right)$ queries. Consequently, the minimax oracle complexity is $\Theta(d\sqrt{\beta R^2/\epsilon})$ throughout the moderate-accuracy regime $\beta R^2(\log(ed)/d)^{2/3}\leq\epsilon\leq c\beta R^2$ for a universal constant $c>0$. The lower bound must account for the fact that a single exact real value can encode arbitrarily much information. To overcome this difficulty, we construct a single fixed smooth convex hard instance using a Moreau-smoothed biased max chain, an exact prefix-shielding mechanism, and batched delayed rotations. These techniques preserve consistency with the full adaptive transcript and establish the optimality of the square-root complexity branch for deterministic bounded-query algorithms.

The article describes a technique for solving multiobjective stochastic optimization problems. As a generalized model of a stochastic system to be optimized a vector (input-random output) system is used. Ran- dom outputs are converted into a vector of deterministic performance and risk indicators. The problem is to find those inputs that correspond to a Pareto-optimal values of output indicators. The problem is approximated by a sequence of deterministic multi-criteria optimization problems, where, for example, the objective vector function is a sample average approximation of the original one and the feasible set is a discrete sample approximation of the feasible inputs. Approximate optimal solutions are defined as weakly Pareto-efficient ones within some vector tolerance. Convergence analysis includes establishing convergence of the general approx- imation scheme and establishing conditions of convergence with probability one under proper regulation of sampling parameters. The proposed technique is illustrated on a computer system for supporting multi-criteria optimization of insurance business.

Communication is a central bottleneck in distributed optimization, but its effect is often summarized by the spectral gap of a chosen mixing matrix. Since this gap depends on link weights as well as topology, it can obscure the intrinsic effect of the network. We clarify this distinction by relating graph diameter to the inverse spectral gap. We establish a universal one-sided bound that is tight up to constant factors, while constructions across several graph families show that no topology-only converse or universal two-sided scaling law exists. This motivates a class-based view of network optimality: diameter-dependent guarantees can be assessed uniformly over graph classes, whereas spectral-gap optimality requires both the graph–matrix pair class and the admissible algorithm class to be specified. Based on the resulting diameter-based minimax benchmark, we introduce Tree-Routed Gradient Tracking (Tree-RGT), a distributed method for stochastic nonconvex optimization. The method pipelines model dissemination and gradient aggregation over a rooted shortest-path spanning tree. Each iteration uses one round of one-hop communication and one stochastic-gradient evaluation per agent, without inner consensus or multi-gossip. Under smoothness and unbiased bounded-variance oracles, Tree-RGT achieves uniformly optimal network dependence over connected graphs of prescribed size and diameter. It recovers centralized stochastic scaling after $\mathcal{O}(nD_{\mathcal{G}}^2)$ transient iterations, where $n$ is the number of agents and $D_{\mathcal{G}}$ is the graph diameter. This network-dependent transient-iteration bound improves upon or matches those reported for representative distributed methods. Thus, topology-aware routing attains diameter-optimal guarantees without relying on a prescribed mixing matrix.

We study an optimal control problem for a two-equations Keller-Segel system on a bounded domain with Neumann boundary conditions, where the chemotactic sensitivity is modulated by a control $f(x,t)$ acting on the taxis flux $\nabla\!\cdot(f\,u\,\nabla v)$. We prove the well-posedness of the state equation under low regularity assumptions. The analysis is based on the study of an $\varepsilon$-regularized problem, the application of Schauder's fixed point theorem, $\varepsilon$-independent, positivity and a priori estimates, the passage to the limit as $\varepsilon \to 0$, and a uniqueness argument. For a cost functional combining trajectory tracking, taxis penalization, control regularization, and terminal objectives, we establish the existence of optimal controls and derive the first-order optimality conditions through a Lagrange multiplier framework, leading to an adjoint system and a variational inequality for the optimal control. A numerical scheme based on spline collocation in space and Runge-Kutta time-stepping is set and implemented. Numerical experiments are performed to simulate the free and controlled dynamics in conservative and non-conservative cases.

Selective state space models (SSMs) have recently emerged as a compelling alternative to transformers, combining competitive performance with substantially improved inference efficiency. At each SSM layer, a sequence of hidden states are propagated by a recurrence, mixing information of different tokens. Despite using a different mechanism, this mixing plays a role analogous to attention in transformers. In fact, recent works have shown that the two architectures may be closer than they first appear, as this recurrence admits a formulation akin to linear attention. In transformers, attention is known to drive the tokens to cluster, i.e., to reach consensus, collapsing in the limit to a single direction. Thus, we ask: does the recurrence at the core of SSMs drive the tokens to consensus, as attention does in transformers? To answer this question, we take a dynamical systems perspective on SSMs, modeling the evolution of tokens across layers as an ordinary differential equation. By exploiting input-to-state stability arguments, we establish local exponential stability of the consensus equilibria and characterize their domain of attraction for time-varying weight matrices, a setting not addressed by previous results. We thereby show that the resemblance between SSMs and transformers does run deeper: the recurrence at the core of SSMs aggregates tokens just as attention does. Numerical experiments on a pretrained Mamba-2 model point to the output gate as the component that regulates the extent of this consensus, preventing the tokens from reaching it in full.

We present a solution of the limiting five-expert prediction problem with geometric stopping, based on stochastic calculus and partial differential equations. If $\delta$ is the stopping probability, the minimax expected regret from tied initial scores is $45\pi^2/(512\sqrt{2\delta})+o(\delta^{-1/2})$ as $\delta\downarrow0$. The adversarial rank control selecting the best and third-best experts maximizes the limiting Hamiltonian at every state. Following the four-expert construction of Bayraktar, Ekren, and Zhang (2020), we represent the value correction as a discounted boundary-local-time expectation for a degenerate obliquely reflected Brownian motion. The hyperbolic systems governing its boundary traces are derived and solved explicitly. To verify the nonlinear Hamilton–Jacobi–Bellman equation, we combine heat-equation minimum principles in equality directions, ODE minimum principles applied to convex combinations of multiple controls, and comparisons across controls. These arguments reduce the 48 regional control inequalities to a few lower-dimensional boundary problems, with the remaining signs established by explicit monotonicity arguments and by factorization of multivariate integral kernels that are positive polynomials in some of their variables. We prove global $C^2$ regularity across both regional interfaces and changes in rank ordering. We also show that the alternating COMB control attains the limiting Hamiltonian maximum exactly when the two highest scores coincide and the third- and fourth-highest scores coincide.

The rapid expansion of data centers presents both economic opportunities and systemic risks for Caribbean Small Island Developing States~(SIDS). This paper develops a PyPSA-based capacity expansion framework for optimal co-location of data centers with renewable energy infrastructure, incorporating behind-the-meter battery storage, workload flexibility (70%/30% inflexible/flexible split), and multi-bus transmission topology. Two case studies are examined: Trinidad~&~Tobago (single-bus, 92\,MW Brechin Castle solar, 5\,MW data center) and Jamaica (4-bus 138\,kV JPS network, 15\,MW data center at Kingston). Jamaica's optimum deploys 405\,MW solar, 278\,MW wind, 800\,MWh grid-scale BESS, and 84\,MWh behind-the-meter BESS, achieving 42.5% renewable penetration; Trinidad's gas-anchored optimum builds no new RE, and a sensitivity analysis shows solar entering at gas costs above \$65--70/MWh, a carbon price of \$25–50/tCO$_2$, or solar CAPEX of \$800/kW. A post-hoc distributionally robust stress test evaluates hurricane resilience, and a novel 10-factor siting framework addresses the regulatory gap in SIDS that lack data-center-specific zoning. Results show that energy-optimal locations may score poorly on societal criteria—noise, freshwater stress, zoning readiness. This demonstrates the need for integrated planning tools in small island contexts.

We introduce a new single-loop algorithmic framework for smooth nonconvex–concave minimax optimization. The resulting projected damped extragradient method combines projected extragradient updates, dual momentum, and a moving proximal center. Under both the optimization-stationarity and game-stationarity criteria, our method achieves the best-known complexity among single-loop first-order methods. For optimization stationarity, our method achieves a gradient complexity of $O(L^2D_Y\bar\Delta_0\varepsilon^{-3})$, where $L$ is the gradient Lipschitz constant, $D_Y$ bounds the diameter of the dual feasible set, and $\bar\Delta_0$ is an initialization quantity involving the value-function gap and the initial gradients. Moreover, by incorporating a fixed-center warm-up phase, the complexity can be improved to $O(L^2D_Y\Delta_\phi\varepsilon^{-3})$, up to an additive lower-order cost, where $\Delta_\phi:=\phi(x_0)-\inf_x\phi(x)$. We further establish a lower bound of $\Omega(L^2D_Y\Delta_\phi\varepsilon^{-3})$ for optimization stationarity over projected zero-respecting first-order methods. This lower bound proves that the warm-started version of our algorithm is optimal up to a constant factor for optimization stationarity within this oracle class. For game stationarity, our method achieves $\mathcal{O}\!(L^{3/2}D_Y^{1/2}\Delta_\phi\varepsilon^{-5/2})$ gradient complexity. This matches the best-known complexity of multi-loop first-order methods, thereby establishing the same complexity with a single-loop algorithmic structure. Under dual strong concavity, the proposed framework achieves $O\!(\sqrt{\kappa}\,L\Delta_\phi\varepsilon^{-2})$ leading complexity for both stationarity criteria, where $\kappa=L/\mu$ is the dual condition number, up to an additive initialization cost. The $\varepsilon^{-2}$ accuracy dependence is optimal under fixed regularity and initialization bounds.

We reformulate average-cost (ergodic) control of a Markov jump process as a finite-horizon stochastic optimal transport (SOT) problem that jointly optimizes over the controlled evolution and the marginal law from which it starts and returns to (i.e., a self-transport). The constraint relating the marginal flow to the controlled generator is the basic adjoint relationship (BAR), so we call the resulting problem BAR-SOT. For any time horizon T>0, its optimal value is a constant scaling of the long-run average-cost rate. We give three equivalent formulations (through controlled processes, a Fokker–Planck constraint, and relaxed marginal measures) and show that the optimal dual is a stationary potential plus a term linear in time, with slope equal to the rate. A relative-entropy penalty on the control yields a cost-tilted Schrodinger bridge problem, computable by a Sinkhorn-type iteration when every transition rate is controlled, whose value converges to the unregularized optimum as the penalty vanishes. We develop the theory for finite Markov decision processes and then for general controlled Markov jump processes. A neural parametrization of the dual, trained as a physics-informed neural network (PINN), that encodes Harrison's equivalent-workload formulation (see Harrison 2000) matches the strong reinforcement-learning baseline of Dai and Gluzman (2022); numerically conditioning the fit of a sub-dominant transverse correction then improves on both that baseline and the best priority heuristic. On the input-queued switch, a graph-attention parametrization of the dual improves on the strongest matching heuristic we are aware of, with a parameter count that does not grow with the switch size.

Coupling-matrix synthesis of cross-coupled resonator filters is difficult when the desired topology is noncanonical. Automatic differentiation removes the need for problem-specific sensitivity expressions, but does not by itself overcome the nonconvex optimization landscape. This work develops a staged differentiable synthesis framework and identifies the importance of its initial objective through a controlled comparison that varies only the objective of the first optimization phase. For tenth- and sixteenth-order benchmarks, directly matching a frequency-sampled response under the prescribed topology reaches the target from $0$ of $24$ and $0$ of $12$ random starts, respectively. Replacing the first-stage objective with characteristic-polynomial matching increases the success rate to $18$ of $24$ and $12$ of $12$ (exact paired $p=7.6\times10^{-6}$ and $4.9\times10^{-4}$). The observed difference is not explained by local conditioning and is consistent with saturation of the bounded response residual far from the solution. A differentiable front end converts symmetric multiband specifications into target polynomials and, when applied end to end to a fabricated dual-band filter specification, finds a realization with one fewer cross-coupling than the published design. Finally, differentiating through Monte Carlo tolerance simulations enables yield-driven design centering, increasing the model-predicted yield from $50.9\%$ to $67.0\pm0.5\%$ under an additive coupling-error model, an improvement of $16.1$ percentage points.

This paper proposes a deterministic bilevel optimization framework for the coordinated operation of Virtual Power Plants (VPPs) embedded in an active distribution network. The Distribution System Operator (DSO) acts as the upper-level leader, minimizing a weighted combination of expenditure, active losses and voltage deviation subject to nonlinear AC power flow constraints, while each VPP operates as a lower-level follower that maximizes its profit under the uniform price signal issued by the DSO. Unlike most existing formulations, which linearize the lower-level subproblem to obtain a Mixed-Integer Linear Program (MILP), the proposed model retains the full AC Optimal Power Flow (AC-OPF) equations, producing a bilevel Mixed-Integer Nonlinear Program (MINLP). The lower-level problem is replaced by its Karush-Kuhn-Tucker (KKT) optimality conditions and the Strong Duality Theorem, yielding a single-level Mathematical Program with Equilibrium Constraints (MPEC). Complementarity conditions are then linearized via the Fortuny-Amat big-M transformation. The framework is validated on the IEEE 33-bus feeder over a 24-hour horizon, with four distributed resources aggregated into a single VPP. Compared with individual dispatch against a regulated time-of-use tariff, aggregation reduces active losses by 10.5 %, the accumulated voltage deviation by 18.3 %, and the bus-hours below 0.95 p.u. from 132 to 29. These gains cost 0.31 % in social cost and 0.26 % in DSO expenditure, while the rent of the aggregator is preserved.

We give an exact algorithm for minimizing an arbitrary rational quadratic polynomial $x^T Q x + c^T x + \gamma$ over the integer points of a bounded rational polyhedron $\{x : Ax \le b\}$. For $n$ variables and $m$ inequalities, the running time is $$2^{O(n \log(n+1))} (m+1)^{O(n)} \beta^{O(n)} (1+L)^{O(1)},$$ where $\beta$ is one plus the maximum binary encoding length of an entry of $A$ or $Q$, and $L$ is the full input encoding length. The polynomial degree in $L$ is absolute; in particular, the encoding of $b$, $c$, and $\gamma$ enters only this fixed-degree factor. Our algorithm refines the integral parallelepiped cover of Goemans and Rothvoss into cells with controlled displacement sets. On each cell, either an integer displacement of negative quadratic value rules out every point as a global minimizer, or the quadratic identity supplies a separation oracle for the convex hull of each integer objective sublevel set. The integer-query feasibility theorem of Hildebrand and Goess then finds the cell minimum. A refinement through scaled lattices solves the displacement search, and a localization argument separates the encoding dependence of the constraint matrix and quadratic part from that of the remaining data. Combining the bounded algorithm with the unboundedness test of Ari and Hildebrand extends the method to arbitrary rational polyhedra.

This paper considers feedback optimization for optimal steady-state operation of nonlinear discrete-time systems when the steady-state input-output map and its sensitivity are expensive to compute. We propose a hybrid extension of sequential feedback optimization (SFO) that augments the model-based SFO gradient with correction terms through a convex combination with summable diminishing weights. Two variants are studied: one based on recursive least-squares (RLS) sensitivity estimation, and another on extremum seeking control (ESC) gradient estimation. Under contractivity and smoothness assumptions, we show that both hybrid schemes preserve the convergence of SFO to a neighborhood of the optimal steady state. The proposed methods are validated through a wind farm power maximization problem using a medium-fidelity model, demonstrating improved early transient performance compared to pure SFO.

We consider the problem of finding stationary points of stochastic convex functions and related variational inequalities. For each, we show that regularized empirical risk minimization, coupled with a random tilting perturbation, obtains stationarity residual order $\sqrt{d/n}$ for $d$-dimensional problems given $n$ observations. We present a few complementary results that show that some dimension dependence is necessary, in distinction from standard stochastic optimization and empirical risk minimization, by providing minimax lower bounds scaling as $\sqrt{\log d / n}$.

We generalize the theory of Wasserstein barycenters between N given measures from the Wasserstein distance to the weak optimal transport problem with barycentric cost. We find a multi-marginal formulation and its dual problem. The multi-marginal problem is new in itself; contrary to the classical multi-marginal optimal transport problem, the first M marginals of its competitors are not fixed but instead only need dominate the first M given measures in convex order, while the last N-M marginals of its competitors are dominated in convex order by the last N-M given measures.

We prove that first-order algorithms require $\Omega(\delta^{-1}\epsilon^{-3})$ gradient queries (in the worst case) to find a $(\delta,\epsilon)$-Goldstein stationary point of a Lipschitz function, at which there is a convex combination of gradients within distance $\delta$ whose norm is at most $\epsilon$. This lower bound is tight, matching known algorithms up to absolute constants, therefore resolving the complexity of convergence to stationarity in nonsmooth nonconvex optimization. We further prove a tight lower bound of $\Omega(\lambda^{1/2}\epsilon^{-7/2})$ for finding points satisfying the recently proposed relaxed notion of $(\lambda,\epsilon)$-stationarity, which allows combining further-away gradients. Our results reveal that convergence rates to nonsmooth stationarity are not affected by gradient stochasticity, in sharp contrast to smooth optimization.

Muon updates matrix-valued neural-network parameters by orthogonalizing a gradient-based momentum matrix. Its reliance on derivatives limits its use when gradients are unavailable or unreliable. We develop a derivative-free framework that constructs Muon-style updates from structured finite differences. Four variants are considered: full entrywise recovery, random low-rank surrogates, basis-aligned rank-one probing, and direct structured search. Exhaustive basis-aligned probing is equivalent, up to positive scaling before ideal polar orthogonalization, to coordinate finite differences. Matrix-regression experiments show that random rank-one probing can reduce the number of function evaluations substantially, at the cost of less accurate updates. Controlled noisy-gradient experiments on regression and a neural network illustrate when accurate function values can compensate for an unreliable gradient oracle. A small CartPole study further examines orthogonal rank-one probes under a fixed episode budget. These results support structured probing as a practical option for selected black-box problems; they do not establish a general convergence guarantee or an advantage over accurate, inexpensive gradients.

We study an infinite-horizon stochastic control problem for the optimal exploitation of an exhaustible resource with unknown total reserves. Information is generated both endogenously through continued extraction without depletion and exogenously through an external information flow. This interaction makes the natural problem non-Markovian and time-inconsistent. We show that it nevertheless admits an equivalent time-consistent formulation with the same optimal controls. The associated value function is characterized as the unique viscosity solution of a stochastic Hamilton–Jacobi equation with random coefficients. We prove comparison on the infinite horizon through a Snell-envelope-based strict-contact argument and establish uniqueness by an independent Brownian regularization and a BSDE correction, avoiding piecewise Markovian approximations. Finally, we identify the deterministic benchmark and show that, under persistent reserve uncertainty, the rescaled stochastic value function converges to a pathwise deterministic control problem, with the long-run optimal extraction rate determined by the asymptotic hazard rate of the limiting reserve distribution.

The $\textit{Colored Knapsack Problem}$ (ColKP) generalizes the classical Knapsack Problem by partitioning the items into color classes and requiring the selected items to admit an ordering in which consecutive items have different colors. The problem is weakly $\mathcal{NP}$-hard and admits two pseudo-polynomial dynamic programming (DP) algorithms proposed in the literature. These two DP algorithms have worst-case running times $O(b \, n^4)$ and $O(b^2 \, n^3)$, respectively, where $b$ is the knapsack capacity and $n$ is the number of items. We develop the first approximation algorithm for the ColKP. By rounding an optimal basic solution of the linear programming relaxation of its natural integer programming formulation and repairing color feasibility, we obtain a linear-time approximation-algorithm whose worst-case performance ratio is $1/3$. We then reformulate both DP algorithms so that profit, rather than knapsack capacity, indexes their pseudo-polynomial dimension, and combine them with profit scaling to obtain two fully polynomial-time approximation schemes (FPTASs). The first FPTAS runs in $O(n^5/\varepsilon)$ time for nonnegative profits and in $O(n^6/\varepsilon)$ time for arbitrary integer profits. The second FPTAS runs instead in $O(n^5/\varepsilon^2)$ and $O(n^7/\varepsilon^2)$ time, respectively. The approximation guarantee, along with new structural insights, provides the bounds needed to control the scaled profit range and establish these running times.

We study computable error bounds and certified early stopping for regularized inverse problems, where a data-fidelity term is traded against a regularizer. The analysis relies on an exact duality-gap identity that splits the total gap of $F(\Phi\mu)+\lambda R(\mu)$ into a data-fidelity Fenchel–Young loss and a regularizer Fenchel–Young loss, $ \Delta(\mu,h)=L_F(\Phi\mu\parallel h)+\lambda L_R(\mu\parallel\eta),\qquad \eta=-\Phi^\star h/\lambda, $ valid for any primal point $\mu$ and any dual point $h$. The data-fidelity term $F$ is strictly convex, so wherever $F^\star$ is differentiable the loss $L_F(\Phi\mu\parallel h)$ is the Bregman divergence of~$F$ between the prediction $\Phi\mu$ and $\nabla F^\star(h)$, and it vanishes exactly at Mirror Alignment $h=\nabla F(\Phi\mu)$. Evaluated at a dual-feasible point $\tilde h$, the gap~$\Delta(\mu,\tilde h)$ is computable and oracle-free, meaning that it uses no knowledge of the solution, and it bounds the suboptimality of $\mu$. Under the source condition, the same Fenchel–Young losses give a priori bounds on the estimation and prediction errors. Their scale is the irreducible model and noise error $L_F(\Phi\mu^\star\parallel h^\star)$, which vanishes exactly when Mirror Alignment holds at the certificate. This gives an early-stopping rule: run the algorithm until the regularizer Fenchel–Young loss falls below a tolerance $\epsilon$. A constructive version of the Brøndsted–Rockafellar theorem then turns the current pair into an exact dual-feasible one, and this proxy lifts to an exact primal certificate. We build the proxy by a proximal step in the geometry of the fidelity, with Bregman kernel $F^\star$ and tilted by the prediction $\Phi\mu$: it recovers the Euclidean step of Carlier when $F$ is the squared error, and it reduces the duality gap by the regularizer Fenchel–Young loss, up to a second-order remainder that vanishes in the quadratic case. Our running example is the Generalized Beurling–Lasso (GBL), where $R$ is the total-variation norm on signed measures. It contains the classical Beurling–Lasso, obtained with the squared error, and also covers robust, logistic, entropic and inverse-optimal-transport losses. The same duality gap certifies deep-learning optimizers such as Lion-K and Muon, in their proximal form, as solvers of the regularized program. A companion paper by the same authors builds on these error bounds to establish exact support recovery for the GBL under a non-degenerate source condition.

Modern large-scale machine learning tasks often require multiple workers, devices, CPUs, or GPUs to compute stochastic gradients in parallel and asynchronously to train model weights. Theoretical results typically distinguish between two settings: (i) the homogeneous setting, where all workers have access to the same data distribution, and (ii) the heterogeneous setting, where each worker operates on different data distributions. Known optimal time complexities in these settings reveal a significant gap, with far more pessimistic guarantees in the heterogeneous case. In this work, we investigate whether these pessimistic optimal time complexities can be overcome under different assumptions. Surprisingly, we show that improvement is provably impossible under widely used first- and second-order similarity assumptions for any randomized algorithm. We then turn to the interpolation regime and demonstrate that the weak interpolation assumption alone is also insufficient. Finally, we introduce a minimal combination of irreducible assumptions, strong interpolation and the local Polyak-Lojasiewicz condition, to derive a new time complexity bound that matches the dependence on worker computation times in the best-known result in the homogeneous setting, without requiring identical data distributions.

Optimizing industrial process flowsheets is often computationally prohibitive due to the high cost of rigorous simulations and the curse of dimensionality inherent in complex design spaces. To address these challenges, we present a reduced-space multi-fidelity Bayesian optimization (RS-MFBO) framework designed for high-dimensional, expensive black-box functions. The approach integrates Global Sensitivity Analysis (GSA) for dimensionality reduction with a fidelity-augmented Gaussian process that captures correlations between low-cost approximations and expensive high-fidelity evaluations. A cost-aware acquisition strategy, augmented with cooldown and promotion mechanisms, adaptively guides the allocation of samples across fidelities. The framework is validated on two distinct industrial process simulators: a plasmid DNA bioprocess in SuperPro Designer and a green fuel synthesis plant in Aspen HYSYS. Results across diverse economic and physical objectives demonstrate that the proposed method substantially reduces the number of high-fidelity simulator evaluations while maintaining competitive optimization performance compared to single-fidelity baselines. These results highlight RS-MFBO as a scalable, simulator-agnostic approach for cost-constrained black-box optimization.