Papers

We provide a tutorial overview of receding horizon control across deterministic, stochastic, multi-agent, and game-theoretic settings. We contrast the minimizer of an optimal control problem with the equilibrium of a game in terms of cost and constraint handling to motivate the difference of MPC schemes involving multiple agents/players. Interestingly, dissipativity theory and the turnpike property have turned out to be a unifying thread and system-theoretic backbone across all these schemes. This paper is the first to give an overview of results, summarizing fundamental insights, drawing parallels, and highlighting technical differences in dissipativity and turnpike analysis across economic, game theoretic, and stochastic settings. We also explore recent trends in numerical methods for games. We conclude by pointing to a list of open challenges.

An exact Max-Cut semidefinite relaxation can have optima outside the convex hull of optimal cut matrices. We study when this hull equals the optimal face and how to recognize equality. For rational weights, the unique optimal dual slack is recoverable in polynomial time, but its kernel may contain directions no optimal matrix uses. Under strict complementarity, recognition and recovery of all optimal cuts on equality instances take deterministic polynomial time. At dual rank one, this hypothesis is unnecessary. Recognition is polynomial-time Turing equivalent to exact rational semidefinite feasibility, even on simple unweighted graphs. Testing whether a supplied cut matrix is the unique optimum is many-one equivalent, even on connected simple graphs with a unique optimal cut. NP-completeness of recognition would imply NP equals coNP. Replication realizes doubly nonnegative matrices with positive diagonal and rational kernel as sign-spanned optimal dual slacks, preserving rank, complete positivity, and cp-rank.

Polyak step size (PS) and Armijo line search (ALS) have received increasing attention in stochastic optimization, with encouraging empirical performance and theoretical guarantees. However, their convergence theory for stochastic heavy ball (SHB) methods remains limited. In this work, we develop a unified convergence analysis for SHB equipped with PS and ALS. To this end, we introduce a modified Armijo rule that closely parallels the Polyak step size, together with a decoupling analysis that isolates the historical dependence induced by momentum. For SHB with standard PS and ALS, we establish expected convergence for strongly convex, convex, and non-convex objectives without interpolation or restrictive conditions on the momentum parameter. Under interpolation or strong growth, we further strengthen the results to almost sure rates and last-iterate convergence. Moreover, for general settings beyond interpolation, we prove almost sure convergence to the exact optimum or to stationarity for SHB with diminishing variants of PS and ALS. These results provide a more comprehensive theoretical view of Polyak step size and Armijo line search for stochastic heavy ball methods.

We study federated learning (FL) with adversarial clients, where the goal is to minimize the average loss of the honest (non-adversarial) clients without knowing their identity. Under heterogeneity, a single shared model parameter is statistically inappropriate: it cannot capture the distinct data-generating processes across clients, incurring an irreducible model-heterogeneity bias and severely limiting robustness to adversarial clients (a.k.a. Byzantine-robustness). We address this problem through representation learning, where each client learns a personalized linear head, while collaboratively estimating a shared nonlinear representation through Byzantine-robust aggregation. We demonstrate that the heterogeneity among honest representation gradients is controlled by the representation error and statistical errors that decay either with the number of data samples per client ($\tau$) or the number of iterations ($T$). In particular, our non-asymptotic parameter recovery error bound reveals three terms: (i) an initialization-dependent error that goes away with $T$, (ii) finite-sample noise terms that decreases with $\tau$ and the number of honest clients, and (iii) a stochastic gradient variance term that also reduces with $T$. Importantly, with no irreducible model-heterogeneity bias in our bounds. We extend the regression analysis to multiclass classification, and empirically validate it on CIFAR-10, FEMNIST, and School Exam Score datasets.

Classical stochastic approximation methods rely on estimators of the first moment (mean) of a random regression function. We study methods that employ estimators of both the first and the second moments, which include modern deep-learning optimizers such as Adam and Muon as special cases. We derive second-moment stochastic approximation methods through the lens of optimal preconditioning for solving matrix equations, and develop a two-stage framework for their convergence analysis. The first stage focuses on the analysis of conceptual (impractical) methods that rely on the exact first and second moments. In the second stage, we replace the exact moments with their respective estimators, and invoke Dvoretzky's theorem to show that the resulting practical methods converge almost surely to a neighborhood of the target solution. The size of the neighborhood depends on the biases and variances of the first- and second-moment estimators. We derive concrete bounds for Muon and a spectral variant of Adam that determine the radius of their neighborhood of convergence.

Agentic reinforcement learning (ARL) with verifiable rewards improves the ability of large language models (LLMs) to tackle knowledge-intensive tasks by learning to interleave search and reasoning. However, most existing ARL methods optimize only LLM-generated tokens and treat retrieved evidence as environment observations. This creates an information-credit gap: failures caused by missing or misleading evidence are attributed to the LLM policy rather than to the retriever, which motivates training the LLM and the retriever jointly. In this paper, we show that retrieval and LLM policy learning are order-sensitive: adapting the retriever before optimizing the policy yields a larger reward gain than the reverse order. To preserve this hierarchy while allowing both components to co-adapt, we formulate retrieval-augmented agentic RL as a bilevel optimization problem. To solve it efficiently, we introduce BRIDGE, a memory-efficient first-order bilevel method motivated by a loss-landscape analysis of the RL and retrieval objectives. Across seven open-domain QA benchmarks, BRIDGE achieves the highest average accuracy with both 3B and 7B backbones, improving the multi-hop average over the strongest baseline by 9.6 and 3.4 EM points, respectively. It also achieves the best averaged answer accuracy and reasoning quality across medical QA benchmarks.

We study the realization problem for discrete-time input/output polynomial systems. These are formalized using tools from commutative algebra and algebraic geometry as systems whose state spaces are algebraic varieties, or more abstractly the set of k-points of an affine k-scheme, where k is an arbitrary infinite field. The input/output behaviors of such systems are described by "polynomial response maps" in which outputs are polynomial functions of past inputs. The main results show that every polynomial response map admits a canonical (quasi-reachable and algebraically observable) realization, which is unique up to isomorphism. The key to the approach is to linearize dynamics by considering a "dual" system in which states are functions defined on states. Finite dimensionality of the canonical realization, and its polynomiality, are characterized in terms of the space, algebra, and field of observables of the map, as well as in terms of algebraic input/output difference equations, and by a Jacobian rank criterion. A particular subclass consists of the maps that we call "bounded," defined by the property that their degree in the past inputs is uniformly bounded. Bounded maps are shown to be finitely realizable if and only if they are realizable by finite-dimensional state-affine systems, whose theory in turn reduces to that of rational formal power series. We also study the lattice of quasi-reachable realizations of a given map, including normal realizations. This work is an update of the PhD thesis written by the author in 1976; connections to recent work are briefly discussed in the last section.

This study investigates generalized boundary representations and transfer learning for physics-informed neural-network (PINN) estimation of spatially varying convective heat transfer coefficients (HTCs) in jet impingement cooling. Piecewise-linear (PWL) and neural-network boundary estimators are inferred from transient interior-temperature measurements and evaluated against CFD-CHT reference HTC profiles at Reynolds numbers of 5000, 10000 and 15000. A discrete Adjoint-FEM formulation with gradient regularization based on the $H^1$-seminorm provides an independent inverse-method comparison. Transfer learning is investigated under two source-target conditions: transfer between noise levels at the same Reynolds number, and transfer from a controlled Gaussian-profile problem to CFD-CHT-derived target problems at different Reynolds numbers. For the representative Re = 5000 target with 5% measurement noise, quadratic, PWL, and neural-network PINNs achieve relative HTC-profile errors of 5.08, 5.31 and 5.19%, respectively, while generalized representations reduce the error for the Re = 15000 target from 16.76% for the quadratic PINN to 14.41% and 14.46%. Transfer learning reduces target-training time by 57.86-81.71% for fixed-Reynolds-number transfer and by 50.59-84.94% for Gaussian-to-CHT transfer. Its effect on reconstruction accuracy is target dependent: fixed-profile PWL transfer remains within 0.30 percentage points of scratch training, whereas Gaussian-to-CHT transfer produces differences of 0.10-1.52 percentage points across targets and representations. These results demonstrate that generalized PINN representations can accommodate spatially varying HTC profiles without a prescribed quadratic form, while transfer learning can substantially reduce target-training cost for repeated inverse estimation, with the resulting reconstruction accuracy depending on source-target compatibility.

We introduce the T-Adaptive Segment Routing, a multi-period optimization problem that emerges in the traffic engineering of core IP/MPLS networks when scheduled maintenance operations are involved. The problem seeks a sequence of Segment Routing (SR) configurations that can adapt to a multi-period scheduled maintenance, while allowing limited number of path reconfigurations between successive time steps. Instead of only minimizing the classical Maximum Link Utilization (MLU) criteria, we propose a more refined lexicographic objective that minimizes the sorted vector of link loads in its entirety, thereby enabling a more efficient resource utilization across all links. To tackle this problem, we develop three exact formulations for solving the resulting lexicographic optimization problem. Among these, the Stela and Carla formulations provide the most favorable computational basis for trajectory-based optimization. We thus develop a trajectory column-generation scheme for these formulations, using pricing oracles that vary from heuristic to exact. Exact pricing gives an exact solution of the LP relaxation of the trajectory master, but it does not itself give a certificate of integer optimality: solving the integer master over the columns generated at the root node only certifies optimality within the generated column pool. To overcome this gap, we incorporate trajectory column generation into a Branch-and-Price framework. The branching rule operates on aggregate variables representing segment usage, which correspond to the original compact routing variables. Thus, Branch-and-Price guarantees global integer optimality for every Stela and Carla rank, and when all ranks are solved to exactness, it guarantees lexicographic optimality of the final solution.

While nonlinear systems can be generally described as infinite-dimensional linear systems via their Koopman semigroups, the direct synthesis of state observers on the premise of learned Koopman operator models, in an analogous manner to finite-dimensional linear systems, remains an open problem. In this paper, the concepts of Koopman observability (exact and approximate), Koopman detectability (exponential and strong), and Koopman–Luenberger observer of nonlinear systems are provided. By defining the Koopman semigroup on a reproducing kernel Hilbert space (RKHS), the conditions for Koopman exponential detectability are established, which involves solving an operator Lyapunov inequality that only implicitly involves the the infinitesimal generator of the Koopman semigroup. Through the use of a finite-horizon Yosida approximation of the Koopman generator and an empirical sampling operator, the data-driven solution of the observer gain can be reduced to a finite-dimensional subspace specified by the sample and hence fulfilled by convex semidefinite programming. The state observation error is guaranteed to be bounded in mean squares, under closed-loop dissipation conditions on the sampling operator and its empirical operator.

Motivated by networked systems in which equilibrium conditions apply to a regulated state rather than directly to the decision variable, we study transformed-state general variational inequalities. Using a projection-residual reformulation, we develop residual feedback methods that operate in the decision space without inverting the state mapping $H$. For a one-step method, we establish an $\mathcal{O}(T^{-1/2})$ best-iterate residual rate under solution-restricted cocoercivity and linear convergence under solution-restricted strong monotonicity and residual Lipschitz continuity. We also obtain an $\mathcal{O}(T^{-1/2})$ best-iterate residual rate for a predictor-corrector method applied to the induced GVI residual under monotonicity and Lipschitz continuity. The assumptions are imposed directly on the computable residual and are illustrated by affine network and rank-deficient examples. Finally, we extend the residual framework to generalized quasi-variational inequalities with decision-dependent feasible sets.

This study develops a projection-based descent method for minimizing weakly semismooth functions over polyhedral sets, taking an initial step toward extending projection-based descent algorithms to nonsmooth, nonconvex optimization problems. At each iteration, the method considers an inner approximation of the Clarke $\varepsilon$-subdifferential at the current iterate over the feasible region and uses its least-norm element to generate a search direction. An exponential limited backtracking line search is then developed to assess the quality of the generated direction, while a projected variant of Mifflin's line search is proposed to identify nonredundant subgradients that enrich the Clarke $\varepsilon$-subdifferential approximation when necessary. To quantify stationarity, a characterization of stationary points based on the projection operator yields a computable optimality measure for the proposed method. The convergence properties of the proposed method are established under mild assumptions. If the algorithm generates infinitely many serious steps, we identify two subsequences such that every cluster point of each subsequence is stationary. If the number of serious steps is finite, we show that the final serious step generates a stationary point. Numerical experiments demonstrate the efficiency and broad applicability of the proposed method across a wide range of test problems, including applications in image denoising, multiobjective optimization, and data clustering. To extend the applicability of the proposed method to problems with general smooth nonlinear constraints, we propose a heuristic approach based on sequential linearization.

We propose an algorithmic framework for solving adaptive robust optimization with provable guarantees on both tractability and solution accuracy. The framework introduces soft separation, a probabilistic mechanism for identifying worst-case uncertainty realizations via a time-inhomogeneous Markov chain. Rather than solving an exact separation problem in each iteration, which is intractable in general, the chain carries adversarial information across iterations, co-evolves with the optimization iterates, and recovers exact separation at terminal iterates with high probability. For continuous first-stage (here-and-now) decisions, we design a first-order method that uses soft separation to produce adaptive gradient estimates. Notably, we prove polynomial-time convergence in expectation to the global optimum. For mixed-integer here-and-now decisions, we embed soft separation within a branch-and-cut framework to generate valid cuts for the robust objective and obtain a high-probability certificate of global optimality. Numerical experiments demonstrate that the proposed methods scale favorably with problem dimension and scenario size relative to state-of-the-art approaches. The instances and code are available online at https://github.com/xuqy2002/SoftSeparationARO.

We exhibit two vector fields on $S^2$, which generate the entire Lie algebra of algebraic vector fields. Compositions of the flows associated with these two vector fields enable us to approximate every diffeomorphism of $S^2$ isotopic to the identity in the $\mathrm{C}^r$-sense for any $r\geq 0$.

Does representation learning stop when the training loss stops improving? We study this question for matrix Muon, whose polar-normalised updates have a step length set by the gradient's rank rather than its norm. Near the edge of stability, full-batch Muon on teacher-student problems enters approximately period-2 loss oscillations that persist for thousands of steps: the cycle-mean loss stays flat or rises, yet the weights keep moving and the learned features continue to align with the teacher subspace. For linear teacher-student learning toys, we derive explicit cycle and alignment formulas and conditional plateau and decay bounds. For a population mean-field ReLU model, we prove that, under stated dimension, initialisation and small-head conditions, the leading eigenspace of the average gradient outer product (AGOP) recovers the teacher subspace exactly during a loss plateau, before the loss later drops. In all 33 ReLU, GELU and SiLU teacher configurations we study, direction-only alignment metrics show the student AGOP aligned with, or still aligning to, the teacher subspace during the period-2 oscillations; projected head refitting on selected configurations shows that the learned directions are useful for prediction, and further measurements distinguish AGOP alignment from weight-mass concentration. In deep residual ReLU students, freezing the downstream layers while the first layer trains with full-batch exact polar updates recreates a nearly flat cycle-mean loss with improving input-AGOP alignment; freezing and unfreezing switch between this plateau and loss decrease, and the effect is sensitive to momentum and to the choice of orthogonaliser.

Inverse optimal control provides optimality and robustness guarantees without solving the Hamilton–Jacobi–Bellman (HJB) equation, but its application to 3D nonholonomic vehicles has been precluded by the absence of strict control Lyapunov functions (CLFs). Recently, such strict CLFs were constructed for the 3D nonholonomic vehicle in spherical coordinates, which circumvent Brockett's obstruction, on the largest possible domain. Building on this construction, this paper develops a general inverse optimal stabilizing controller for the 3D nonholonomic vehicle actuated by surge velocity, pitch rate, and yaw rate. The design is a damping $L_gV$ feedback that minimizes a meaningful cost and inherits stability margins. Its flexibility is demonstrated through two choices of penalty function: a quadratic running cost on the state, recovering a near-classical optimal control formulation, and stabilization under user-defined input constraints. Numerical simulations illustrate stabilization performance, control-effort tradeoffs, and the effect of input constraints.

We prove unique continuation for the real periodic Benjamin–Bona–Mahony equation in the energy space: a solution that equals a constant on a nonempty space–time strip is that constant everywhere. In particular, Rosier's unique continuation conjecture [23] holds without any smallness or mean assumption. It even suffices to observe, at one point $x_0$ and for $t$ in an interval, the solution and the potential $\psi=(1-\partial_x^2)^{-1}(u+u^2/2)$ of its flux: $u(x_0,t)=C$ and $\psi(x_0,t)=C+C^2/2$ imply $u\equiv C$, and the second condition cannot be dropped when $C=-1$. The proof uses an exponentially weighted Hamiltonian of indefinite sign with a coercive dissipation identity, and it applies to every nonconstant real polynomial flux. Consequently, Rosier's stabilization results for two localized feedback laws and for a boundary feedback hold unconditionally, for arbitrary data. Stationary step profiles show that the theorem fails in $H^s$ for $0\leq s<1/2$. We also give a generating-series proof of single-point unique continuation for the linearized equation and discuss related control results.

When a free energy admits several local minimisers, the associated mean-field particle system exhibits a metastable behavior, undergoing random transitions between these minimisers. However, when the number of particles is large, these transitions are rare events and typically occur at a time-scale which is out of reach for simulations. In this work, we introduce a system of particles coupled with a slow macroscopic external variable which approximately follows a Langevin dynamics associated to a coarse-grained free energy. The empirical distribution of the particle system thus samples all local minimisers of the initial free energy, with transition times which are now independent from the number of particles. The construction is inspired by a theoretical tool: the two-scale approach for proving log-Sobolev inequalities for Gibbs measures.

We study feasibility correction for linear programs. Suppose an initial point satisfies most of the constraints. Can it be moved into the feasible region without changing the objective value? Our main deterministic result says that this is possible if the constraint matrix satisfies a restricted isometry property (RIP). We apply the correction theorem to obtain lower bounds for random linear programs with finite-moment or sub-Weibull entries. We complement these with matching-order upper bounds under a uniform $p$th-moment bound for some $p>2$ and a delocalization condition on the objective direction. As an application of our techniques, we also derive lower bounds for the Gaussian width of feasible regions under RIP assumptions.

In this work, we develop a second-order optimization framework for physics-informed neural networks (PINNs) applied to high-dimensional parametric partial differential equations (PDEs). The framework is built on the Gauss–Newton pullback metric, which provides an operator-informed notion of curvature in parameter space and connects the method to the broader family of natural gradient schemes. We show that, for coordinate-separable neural architectures and linear differential operators (or linearized operators in the nonlinear case) admitting a finite separable representation, the residual Jacobian inherits a structured separable factorization. This yields an exact compressed formulation of the Gauss–Newton step in a reduced space, without assembling the full residual Jacobian on the exponentially large tensor-product collocation grid. The dimension of the reduced space (the effective compressed dimension) is determined by the local collocation grid sizes, the separable operator structure, and the contraction pattern of the architecture, thereby avoiding dependence on the full tensor-product grid size and replacing dense linear algebra in parameter space by a substantially smaller structured problem. Within our framework, we investigate canonical polyadic and tensor-train parametrizations and derive their full algebraic characterization relevant to the Gauss–Newton method, including the structure of the residual Jacobian, the resulting compressed system, and its effective compressed dimension. Numerical experiments on high-dimensional PDEs, including parametric problems, demonstrate the high efficiency of the proposed compressed Gauss–Newton method, which achieves substantially lower errors than tensor-compressed first-order baselines with orders of magnitude fewer iterations and only a fraction of the computing time.

Given a bounded open set $V\subset \mathbb{R}^2$, we consider the problem of slicing this set by a family of curves, whose average length is as small as possible. The first sections of the paper prove the existence of an optimal slicing and establish a connection with the problem of eradicating an invasive biological species from an island, considered in \cite{BMS}. Namely, an optimal slicing can be obtained as a limit of minimum time eradication strategies. The second part of the paper derives various necessary conditions for optimality, describing the behavior of an optimal slicing strategy in the interior of $V$ and near the boundary $\partial V$.

The paper considers a class of optimal set motion problems, where the control determines the velocity of boundary points in the normal direction. Results and counterexamples are given, on the continuous dependence of the minimum cost from the initial configuration. The last sections of the paper study the slicing cost of bounded open sets. For strongly Lipschitz domains, it is proved that the slicing cost varies continuously w.r.t. the Hausdorff distance of the boundaries. We conjecture that the slicing cost for a disc is maximal among all sets with the same area.

Instead of the sequential plant–control–estimator design pipeline, which is in general not optimal, we propose the Control and Estimation Co-Design (ContEst) framework, which poses the entire system design as a single two-stage problem solved by first-order methods. A key enabler is that the gradient of each inner control/estimation optimal cost with respect to the design parameters is available from the dual variables the inner solver already returns, by the envelope theorem: exactly under mild regularity (a unique inner minimizer), and as a subgradient under more relaxed conditions. That is, no differentiation through the optimizer or its optimality conditions is needed. We first present the framework in a linear quadratic Gaussian (LQG) regime, where the control and estimation costs become two semidefinite programs (SDPs) whose design gradients are read from the duals of the Lyapunov constraints. We extend the results to robust ($H_\infty$) formulation to cover worst-case control and filtering problems. We also present an extension to constrained and nonlinear systems, using an extended Kalman filter (eKF) and a local linearization of the information-state dynamics to yield an approximate but convex model predictive control (MPC) problem whose design gradients are read directly from the costates. We apply our methods to a variety of realistic design examples from different fields and report significant design improvements.

Given a functional dependent on the spectrum of a differential operator, we address the problem of finding a domain which optimizes this functional. PDE solvers might be used to tackle this optimization. It is however computationally expensive. We propose two neural network models which learn the spectrum directly from the geometry of the domain and can be used to optimize the domain from one or more eigenvalues. We investigate two representations. The first encodes the domain through Fourier coefficients and a light MLP, which is efficient on star-shaped geometries, achieving a precision of 0.2%. Through a rescaling of the coefficients the designed models satisfy the scaling law of the eigenvalues. Additionally, averaging the outputs of the trained surrogates over rotations and reflections induces invariance for these transformations. The second is a model that takes the landscape function, the indicator function and the gradient of the landscape function. A Gram-Schmidt process produces orthogonal eigenfunctions as output of the model along with the associated eigenvalues. The landscape model reaches 1% mean relative error on the first ten eigenvalues, compared with 4% for an FNO model. Replacing the landscape by an SDF worsened both prediction and optimization errors. The trained model also generalizes from synthetic shapes to domains given as classical image dataset. The resulting surrogates of both approaches recover classical spectral optima such as the disk for the first eigenvalue or the conjectured minima of higher eigenvalues. This confirms that our models produce accurate differentiable estimates of eigenvalues, which can be used in shape optimization problems involving spectral quantities.

Among distinct optimal acceleration mechanisms for deterministic monotone root-finding problems and fixed-point problems, dual-anchoring has recently been shown to admit a more robust direct stochastic extension than standard anchor acceleration. However, without additional strong monotonicity, the existing stochastic dual-anchoring guarantee has two limitations: first, it requires cocoercivity in expectation, and second, it attains only $O(\epsilon^{-3})$ oracle complexity, leaving a gap to the near-optimal $\tilde{O}(\epsilon^{-2})$ complexity achieved by other methods. In this work, we address both of these limitations by combining dual-anchoring with stochastic resolvent approximation and optimized variance control. For unbiased stochastic oracles with variance bounded by $\sigma^2$, where sample operators are monotone and uniformly $L$-Lipschitz, our algorithm finds a point with $\epsilon$-residual with a near-optimal oracle complexity of $O ( (LD / \epsilon) \ell + (\sigma^2 / \epsilon^2) \ell^2)$, where $\ell = \log (1 + LD / \epsilon)$ and $D$ is the initial distance to a solution. This result improves the best known oracle complexity in the noise-dominated regime under these samplewise assumptions, reducing the poly-logarithmic factor from cubic to quadratic.

The lack of rigorous safety and performance certificates remains a key bottleneck to the deployment of modern learning-based methods. Sample compression has recently emerged as a powerful tool for deriving such certificates, with particularly sharp bounds available for algorithms satisfying a so-called preference property – also known as stability in the learning theory literature. These bounds find direct application across domains as different as the Scenario Approach, Pick-to-Learn, and Support Vector methods. However, whether they are tight has remained an open problem. In this paper we resolve this question affirmatively and show that the state-of-the-art bound for preferent compressions is provably tight. We establish this by exhibiting an explicit construction based on the uniform distribution and order statistics that attains the bound in the limit. Along the way, we also provide a considerably shorter and more accessible proof of this bound, requiring only elementary counting arguments and no infinite-dimensional duality.

We develop accelerated gradient-free methods for stochastic convex optimization with constraints defined by expectations. Our batched primal-dual sliding method uses two-point evaluations sharing a random sample and guarantees expected objective error and expected maximum constraint violation at most $\varepsilon$. It achieves $O(\varepsilon^{-1/2})$ sequential oracle rounds for smooth data and $O(d^{1/4}/\varepsilon)$ for nonsmooth data in dimension $d$, recovering the accuracy and dimension dependence of the corresponding unconstrained accelerated methods. The smooth rate is optimal in accuracy. Each round collects all samples needed for its inner primal-dual updates. Total evaluations retain quadratic dependence on inverse accuracy, with only polylogarithmic dependence on the number of constraints under vector feedback. Complementary lower bounds distinguish the dimension cost of gradient estimation from the unavoidable logarithmic cost of estimating noisy constraint levels. For smooth strongly convex objectives, restarts give logarithmic round complexity and objective evaluation cost linear in inverse accuracy, while constraint-level estimation necessarily remains quadratic. Known affine constraints require no constraint-oracle calls.

This paper presents Peppy, an AI-assisted workflow for discovering tight, analytic convergence proofs for first-order optimization algorithms. Generic approaches to using LLMs to conduct mathematical research target an unspecified, broad spectrum of problems and sometimes use the Lean 4 proof assistant for formalization. On the other hand, Peppy leverages domain-specific knowledge more heavily and is thereby capable of constructing the proofs in a more structured manner, which allow a minimal and accessible verification through SymPy. We experimentally demonstrate through examples that Peppy provides a rigorous, practical, and reproducible paradigm for AI-assisted theorem synthesis in optimization. We further highlight its capability of closing several open problems on tight convergence analysis of first-order optimization algorithms, including conjectures for Nesterov's FGM. Overall, Peppy is designed to turn the art of optimization algorithm analysis into a science.

We study exponential output-to-state stability (eOSS) of linear infinite-dimensional systems in Banach spaces with bounded output operators. It is shown that eOSS is equivalent to the existence of a coercive eOSS Lyapunov function in implication form. Also exponential detectability guarantees the existence of a coercive eOSS Lyapunov function in dissipation form and therefore eOSS. Counterexamples demonstrate that the converse implications fail in general: exponential zero-detectability does not imply eOSS, eOSS does not imply the existence of an eOSS Lyapunov function in dissipation form, which, in turn, does not imply exponential detectability. If the unstable subspace is finite-dimensional, zero-detectability implies exponential detectability and yields equivalent characterizations of eOSS. The results are illustrated with a parabolic equation.

Interior-point methods (IPMs) are among the most widely used algorithms for constrained optimization, yet their Newton-based search directions require costly second-order information and large linear-system solves. Learning to optimize offers cheaper updates learned from data, but the singular behavior of logarithmic barriers near constraint boundaries makes IPMs highly sensitive to perturbations, complicating both warm starting and learning reliable updates. We introduce pdLIP, an IPM for smooth nonlinear programs that integrates learned preconditioning with pdProj, an all-shifted primal-dual projected-search IPM. A shared coordinate-wise recurrent network predicts a positive diagonal preconditioner that scales the right-hand side of the reduced Newton system for the primal step, and the remaining slack and multiplier directions are recovered analytically. The learned iterations avoid Hessian evaluations and Newton-system solves, using only first-order and coordinate-wise operations amenable to GPU parallelization. Training is self-supervised, with a loss based on a penalty-barrier merit function and the residual of perturbed optimality conditions, requiring neither target directions nor precomputed solutions. Primal and dual shifts mitigate the barrier's sensitivity to perturbations near constraint boundaries, enabling effective warm starting. Across four classes of 200-dimensional convex and nonconvex constrained problems, pdLIP warm starts reduce pdProj refinement iterations by 63-67% compared with cold starts at the same KKT residual tolerance of $10^{-8}$, with negligible warm-start generation cost relative to the subsequent pdProj solve. Improvements persist on box-constrained QPs with 1000 variables and extend to applications including portfolio optimization, support vector machines, and a nonlinear control example.