Papers

This paper proposes the certainty-equivalent first-order learning (CEFOL) algorithm, a deep learning algorithm for solving discrete-time dynamic programming problems with recursive utility. Dynamic programming with recursive utility is challenging because nonlinear certainty equivalent appears in the Bellman equation and the first-order optimality conditions but is difficult to evaluate. By introducing a separate neural network to represent the certainty equivalent, CEFOL enables the exploitation of the Bellman and model-specific first-order optimality conditions. In addition to certainty equivalent, CEFOL also uses neural networks to learn the value functions, policy functions, and Lagrange multipliers by using model-specific first-order conditions to construct residuals for minimization. By using first-order and KKT residuals to learn the policy, CEFOL directly accommodates general equality and inequality constraints on the controls, including occasionally binding constraints, without requiring penalty functions or problem-specific reformulations. We apply the algorithm to risk-sensitive and Epstein–Zin consumption-saving problems, a small-noise robust-control problem, and a DSGE model with recursive preferences and stochastic volatility. Across these applications, out-of-sample Bellman diagnostics and model-specific optimality residuals, including Euler or first-order residuals where applicable, are generally of order 1.0e-4 to 1.0e-3 over the relevant state regions, with larger values mainly near binding constraints, and the learned value and policy functions closely match VFI benchmarks when available. The CEFOL algorithm also works for dynamic programming problems with expected utility, as expected utility is a special case of recursive utility.

Maximizing a symmetric quadratic form under simultaneous L1 norm inequality and L2 norm equality constraints is a standard and widely used continuous relaxation for Sparse Principal Component Analysis (SPCA). This paper settles the computational complexity of this continuous formulation by proving it is NP-hard. Furthermore, the variant with both L1 and L2 norm inequalities is also shown to be NP-hard.

In this paper, we study the problem of optimizing Markov processes that interpolate between two prescribed probability distributions while minimizing a given cost. The main computational challenge is the curse of dimensionality: in high-dimensional state spaces, representing the full distribution is intractable. To address this, we reformulate the problem in terms of sequential couplings and develop convex relaxations based on local marginals and cluster moments. These relaxations exploit locality and sparse interaction structure, provide computable lower bounds, and recover low-order statistics of the intermediate laws. We identify dynamic optimal transport as a special case of our Markov process optimization problem and develop a procedure for recovering the underlying Benamou–Brenier dynamics from the relaxed solution. We also show that the procedure extends to more general Markov processes and illustrate it with a constrained process between Ising models.

Here we present PCPOP, a Julia package for polynomial optimization that supports non-commutative optimization, tracial polynomial optimization, trace polynomial optimization and state polynomial optimization. PCPOP fully supports exact arithmetic computations and incorporates convenient functionalities such as algebraic reductions based on Gröbner basis methods, automatized symmetrization via Wedderburn decompositions, and Jordan algebra reductions. As a distinguished feature, PCPOP implements a specialized framework for polynomial computations in partially commutative variables that provides significant computational advantages for problems appearing in quantum information.

In this paper, we address control design of Networks of nondeterministic and metric finite state Systems (NoS) where some systems are plants, others act as controllers, and some are uncontrolled processes. The control architecture of NoS extends decentralized control architectures. The problem addressed consists in designing local controllers for enforcing on the corresponding plants, local specifications given by regular languages over plants' state alphabet, up to desired accuracies. The approach used to solve the problem is termed incremental and consists in finding solutions on an increasing sequence of NoS extracted from the original NoS.

We present a polynomial-time algorithm for optimally coloring perfect graphs that is based entirely on graph-theoretic operations. At its core, the algorithm decides whether a perfect graph contains a clique of a given size by iteratively counting walks in the graph with certain weights assigned to its edges and nonedges. These weights are initialized according to a uniform scheme and then updated in each iteration based on the walk counts from the previous iteration.

We study matrix feasibility problems "find X in K intersect C" where K is a closed convex matrix set and C = {X : sigma(X) in S} is defined by a convex constraint S on the ordered singular values. Using the classical spectral transfer identity, projection onto C reduces to an SVD plus a small quadratic program. For seven natural polyhedral S we embed this projector in a plain alternating projection (AP) loop. We experiment with two concrete families of K - linear constraints (an affine subspace intersected with an entrywise box) and ellipsoidal constraints (a non-centered anisotropic Frobenius ellipsoid) - although the method applies equally to more general convex constraints. The experiments expose three regimes: rapid feasibility, slow tail convergence, or informative infeasibility plateaus.

We present a general computational framework for macroscopic nonlocal traffic flow models on networks with multiple commodities. The model combines scalar conservation laws on edges with nonlocal velocity functions and couples them via buffer-based junction dynamics. Routing is defined in general as a prescription for distributing drivers across outgoing road segments. As an example, we implement dynamic k-shortest-path routing, where travel times along roads and waiting times at intersections are used to compute shortest paths to the commodities' destinations at each time step, and drivers are distributed accordingly. In another example, we optimize routing over a considered time horizon to minimize the total travel time. This framework naturally creates a feedback loop between traffic evolution and route choice. Numerical examples, ranging from small test cases to large grid-like networks, demonstrate the robustness of the approach and allow for a comparison of different routing strategies.

We consider $0/1$ packing problems $\max\{c^T x \colon Ax \leq 1, \, x \in \{0,1\}^n\}$, with $A \in \mathbb{R}_{\geq 0}^{m \times n}$. A way to solve such problems is via tightening the linear programming relaxation $P$ with Gomory cutting-planes. The Gomory-closure $P'$ of $P$ is the intersection of $P$ with all its cutting planes. The optimization problem over $P'$ is NP-hard. Mastrolilli (2020) has shown that for fixed ${\epsilon}>0$, the Lasserre hierarchy yields a polynomial-size convex but non-polyhedral extended formulation that approximates $P'$ up to a factor of $1+{\epsilon}$. Our main result is the construction of a polyhedral and polynomial extended formulation that approximates $P'$ with the same approximation guarantee. Our construction is based on first principles. Like Mastrolilli's approach, ours also applies to higher iterates $P^{(t)}$ for fixed $t$ and ${\epsilon}>0$. In contrast to an explicit construction, communication complexity provides an alternative way to describe extended formulations. Using this approach we obtain a quasi-polynomial polyhedral extended formulation for the above problem that is superior in some parameter regimes. To achieve this, we describe a communication protocol extending Yannakakis' protocol to decide whether the clique of Alice and the stable set of Bob intersect.

Schedule-Free methods have attracted growing interest for alleviating the burden of designing and tuning a learning rate scheduler, while matching and sometimes even outperforming optimizers with tuned schedulers. Despite their strong empirical results, their convergence theory in nonconvex optimization, where modern machine learning objectives typically arise, has remained largely unexplored. In this paper, we provide worst-case analyses of Schedule-Free gradient descent and Schedule-Free stochastic gradient descent, in their standard form and without auxiliary modifications or restrictive conditions, for smooth but possibly nonconvex objectives. Based on a Lyapunov analysis derived from the continuous-time limiting ordinary differential equation associated with these methods, we show that Schedule-Free gradient descent and Schedule-Free stochastic gradient descent achieve the optimal worst-case convergence rates attainable among first-order methods. We further formulate Schedule-Free gradient descent as a nonautonomous dynamical system and prove strict-saddle avoidance under an arbitrarily small one-time perturbation. These theoretical results provide a better understanding of the strong performance that Schedule-Free methods demonstrate.

We study the optimal transport of optimally controlled agents from a compactly supported absolutely continuous source to a discrete target measure. The ground cost for the transport is induced by the optimal cost of the agents' motion. When this ground cost satisfies the twist condition, the optimal transport map is given almost everywhere in terms of a Laguerre tessellation of the state space. We refer to this control-theoretic generalization of Laguerre tessellation as Control Laguerre Tessellation (CLT), and illustrate it for two ground costs induced by linear controlled agents with minimum energy and minimum time objectives.

We study conditions under which stability of the origin of stochastic differential equations is robust to small perturbations. We express robustness in two ways, firstly in the sense that stochastic stability is maintained under small parametric perturbations not exceeding a state-dependent bound vanishing at the origin but positive elsewhere, and secondly via stochastic input-to-state stability (ISS) which allows non-zero perturbations everywhere. We prove the former property assuming the existence of a Lyapunov function certifying stochastic stability of the nominal system. Under the same assumption, stochastic ISS holds under a suitable state-dependent perturbation scaling. Stochastic exponential stability is maintained under proportionally bounded perturbations and implies exponential ISS even without perturbation scaling. Finally, we propose a novel approach to stochastic integrator backstepping in pure-feedback form that uses the tools from our robustness analysis.

We study stochastic fixed-point equations $\mathbf{T}(\mathbf{x}) = \mathbf{x}$ over normed spaces $(\mathcal{E}, \|\cdot\|)$, where the operator $\mathbf{T}$ is nonexpansive or contractive and is accessed only through unbiased stochastic evaluations with bounded second central moment. Given $\epsilon > 0, \delta \in (0, 1)$, the goal is to output $\mathbf{x} \in \mathcal{E}$ such that $\|\mathbf{T}(\mathbf{x}) - \mathbf{x}\| \leq \epsilon$ with probability at least $1-\delta$. We introduce VR-GHAL, a variance-reduced gradual Halpern method for quadratically smoothable Banach spaces. The key algorithmic ingredient is a recursive stochastic estimator based on clipped differences of oracle evaluations: instead of clipping $\tau(\mathbf{x}; \xi)$ itself, we clip stochastic differences at the Lipschitz scale $\gamma\|\mathbf{x} - \mathbf{y}\|$. This makes the estimator pathwise Lipschitz along the algorithmic trajectory while permitting martingale concentration under finite second moments in the native norm. Our main theorem gives an anytime high-probability residual bound: on a single event of probability at least $1 - \delta$, the residual decreases nearly geometrically across epochs, up to lower-order logarithmic factors. Under only bounded variance, displaying only the dependence on the target error $\epsilon$ and Lipschitz constant $\gamma \in (0, 1]$ of $\mathbf{T}$, the resulting oracle complexity is $\min\{\epsilon^{-5}, (1-\gamma)^{-3}\epsilon^{-2}\}$. Under a Lipschitz-in-expectation oracle, the dependence improves to the corresponding $\epsilon^{-3}$ nonexpansive rate (i.e., for $\gamma = 1$), and under samplewise nonexpansiveness to $\epsilon^{-2}$.

Interconnected financial systems are vulnerable to cascading failures arising from cross-holdings and nonlinear contagion, making the analysis and mitigation of systemic risk a challenging computational problem. In this work, we develop a unified optimization framework for financial network analysis based on Ising models and Quadratic Unconstrained Binary Optimization (QUBO). Starting from the Elliott Golub Jackson financial network model, we extend equilibrium valuation to incorporate threshold-induced failures, formulate the Maximum Cascade Failure Problem, and derive an equivalent QUBO representation. We then formulate the Optimal Bailout Allocation Problem as a controlled Ising model and transform the resulting bi-level optimization into a single joint QUBO that simultaneously determines equilibrium failures and optimal interventions under budget constraints. To characterize the influence of individual institutions, we introduce bailout susceptibility as a response-based measure of systemic importance and develop a susceptibility-driven greedy intervention strategy. Numerical simulations demonstrate equilibrium valuation, worst-case cascade identification, optimal bailout allocation, and susceptibility analysis on financial networks of varying sizes. Beyond optimization, the Ising representation provides a general statistical-mechanical framework for analyzing financial contagion, enabling the application of response theory, Monte Carlo methods, and other techniques developed for interacting spin systems. The proposed framework establishes a unified approach for systemic risk analysis that is compatible with classical annealing, quantum-inspired optimization, and emerging quantum annealing technologies

Soft metamaterials provide a promising platform for robotics, biomedical devices, and flexible electronics. The localized mechanical responses by nonuniform excitation are ubiquitous in soft materials, yet their controlled transmission across assemblies remains largely overlooked in metamaterial design, which critically constrains nontrivial functionalities with end-to-end and long-range deformation transmission. Here, we introduce an impedance-guided design framework that enables programmable transmission of localized deformation in modular soft metamaterials, achieving behaviors unattainable by intuitive design. By establishing a nonlinear model considering position-dependent interactions and integrating the concept of mechanical impedance within metamaterials, we regulate assembly-level transmission solely through unit-cell topology optimization. The resulting framework enables effective synthesis of module families, allowing both homogeneous and heterogeneous assemblies to be custom-built with markedly enhanced transmission characteristics. Leveraging the highly combinatorial and extensible design space, we physically realize diverse on-demand displacement manipulation architectures, including obstacle-bypassing modular soft-metamaterial assemblies, defect-tolerant soft gripping, and embodied signal processing. Beyond deformation programming, the reconfigurability and reassemblability of these soft modules can embed electric logic signals, enabling energy-efficient and low-latency information processing through compliant-switch-controlled mechanical LED displays and wearable finger-motion-sensing controllers. Our method provides fundamental insights into localized deformation transmission in modular soft metamaterials and establishes a scalable route toward embodied-intelligence material systems, particularly for soft-metamaterial-centric actuation, sensing, and collective computing.

We study nonasymptotic convergence of primal-dual methods for a class of nonconvex constrained optimization problems under local multiplier regularity. Both the objective and functional inequality constraints are compositions of convex Lipschitz outer functions with smooth nonlinear mappings. Unlike existing analyses requiring regularity globally or along the entire iterate sequence, we impose multiplier regularity only near the feasible set. Using an auxiliary minimax problem with a compact dual set, we derive Lyapunov estimates that identify nearly feasible iterates under feasible initialization and a sufficiently large penalty, without assuming multiplier boundedness. Local regularity then bounds the associated multipliers and makes dual truncation inactive, yielding approximate KKT certificates for the original problem. The smoothed prox-linear augmented Lagrangian method with dual regularization yields an $O(K^{-1/3})$ KKT residual rate. Without regularization, piecewise linear outer functions and additional local structural assumptions yield a local dual error bound and the sharper $O(K^{-1/2})$ rate.

We study Model Predictive Path Integral (MPPI) control for nonlinear systems with additive process disturbances whose covariance is unknown, spatially varying, and slowly time-varying. A mismatched disturbance covariance produces a persistent penalty in closed-loop stability certificates, while online estimation can reduce this penalty as data are collected. We propose a cell-wise recursive covariance estimator with spatial diffusion and prove a finite-horizon error bound that separates stochastic-approximation error, spatial-smoothing bias, and temporal-drift effects. The diffusion kernel is chosen to be reversible with respect to the stationary visitation measure, making the diffusion operator dissipative in the weighted Lyapunov analysis. We then substitute the resulting covariance estimate into the MPPI sampling distribution and derive an adaptive stability certificate with an explicit learning penalty. The main result is a payoff theorem: after a computable crossover time, the adaptive controller achieves a strictly tighter certified stability bound than any fixed covariance choice whose mismatch exceeds the residual smoothing and drift allowance. Numerical experiments illustrate the estimator convergence and the resulting stability-tightening effect.

We introduce a framework for tackling questions in discrete optimization associated with parametric constraint matrices. More precisely, the constraint matrices have entries that are polynomials in one variable and all subdeterminants of these matrices are polynomials in a given prescribed set $S$. Two key problems arise in this context. The first is the recognition problem: can a matrix of this form be recognized in polynomial time? The second is the optimization problem: given an integer program whose constraint matrix is of this form, can it be solved in polynomial time? We answer both questions affirmatively for a particular set $S$ consisting of nine linear forms. The matrices we consider are of themselves independent interest; they arise as matrix projections of certain bimodular matrices that admit two distinct unimodular projections.

This work surveys the different parallel variants of the mesh adaptive direct search (MADS) algorithm for constrained blackbox optimization. These problems can inherently imply high computational costs due to the possible large number of variables and multi-modality of the search space. In addition, the potential time-intensive nature and time heterogeneity of the blackboxes defining the problem prompts the need for efficient implementations. Parallelism emerges as an actionable solution to mitigate computation time, as modern computer systems rely on multi-core architecture. The reviewed methods employ diverse levels of parallelism and distinct parallel strategies to effectively tackle each aspect outlined above. The manuscript details the practical implementations, provides computational results, and offers insights into the advantages and limitations of each MADS parallel method.

In this article we investigate the theoretical convergence properties of the fictitious-play approximation procedure applied to coupled FBSDE systems for finite-player non-zero-sum stochastic differential games. Under one set of assumptions, the convergence is shown to be geometric. Under an additional structural assumption, the geometric convergence rate further improves to a super-exponential rate in a special class of games. To the best of our knowledge, this provides the first convergence analysis of fictitious play for fully coupled FBSDEs. A numerical experiment with a linear-quadratic interbank borrowing and lending problem confirms the geometric convergence.

The connected p-median problem is defined as a variant of the classical p-median problem when the facility nodes induce a connected subgraph. In this paper, we introduce the weighted version of the above problem when the weight of the facility connection in the objective function is defined by the minimum weight spanning tree of the facility nodes. This approach is motivated by the sink node selection in distributed sensor networks, in which the collected information is shared among the sink nodes through the minimum spanning tree. The weights of the graph determining the network topology of the candidate sink nodes as connection costs are distinguished from the standard access costs of the p-median problem. The fixed deployment costs for the setup of facilities are also considered. The objective is to minimize the overall cost as the sum of deployment cost, access cost and connection cost. We show that the problem is NP-hard and propose three mixed-integer linear programming (MILP) formulations adapted from the traveling salesperson problem literature. Since these formulations are poorly scalable with respect to network size, we develop a four-phase matheuristic method based on linear programming rounding. We conduct an extensive computational study to evaluate the performance of the MILP formulations and 22 variants of the matheuristic under different parameter settings. The results indicate that the MILP models perform effectively on small instances but struggle to solve medium- and large-scale instances within a two-hour time limit. In contrast, several matheuristic variants consistently produce high-quality solutions within minutes. Finally, we analyze the impact of network structure, size, density, and the parameter $p$ on solution quality, providing further insights for network design.

Semidefinite programming hierarchies for commutative and non-commutative polynomial optimization represent a powerful computational tool with many applications in quantum information. In such applications, a given variable is typically not either commuting or non-commuting with all other variables, but instead commutes with some variables and does not commute with others, i.e., the variables satisfy some partial commutation relations. While such partial commutation relations can always be incorporated in a fully non-commutative setting through suitable linear constraints in the semidefinite programming relaxations, exploiting their algebraic properties from the onset can result in more compact relaxations. This leads us to introduce partially-commutative polynomial optimization, a framework that encompasses commutative and non-commutative polynomial optimization, allowing for arbitrary commutation relations among the variables. We point out that the underlying algebraic structure is that of a partially-commutative monoid. We present and review several key aspects of such monoids and show how they can be used to build SDP relaxations for partially-commutative polynomial optimization problems in which the partial commutations are natively implemented in the monomial structure, without the need of additional linear constraints.

Progress in mathematics often requires more than a certificate of truth: it requires proof structures that are transparent, checkable, and reusable. Automated systems can increasingly certify that a result is true; what they typically return, however, is a dense certificate rather than an interpretable, reusable proof structure. Recent work on performance estimation problems has shown that performance bounds and complexity analyses of first-order optimization methods can be discovered by searching over a structured space of Lagrangian dual certificates. We cast the search for simpler proof structures as a second-stage optimization problem over these certificates. Starting from dual certificates, we develop post-processing procedures using tools from sparse optimization and statistical learning. We measure complexity through features such as active hypotheses and residual structure, and introduce methods based on exhaustive sparsification, weighted $\ell_1$-type heuristics, and semidefinite programming (SDP) formulations for discovering simple proofs and intermediate lemmas. Examples on gradient descent, proximal methods, and fast-gradient methods show that these procedures can autonomously prune redundant inequalities, reveal structured proof patterns, and, in the proximal setting, recover Lyapunov functions as intermediate lemmas that lead to simple, streamlined proofs. By distilling dense machine-generated certificates into compact proof structures, this workflow acts as a pre-processing step for the final proof, reducing the complexity that must be managed during human interpretation, reuse, and formalization.

We study stochastic barycentric estimators for proximal points and metric projections obtained by exponentially reweighting Gaussian perturbations. Our main result is an abstract comparison theorem for probability measures with densities proportional to an exponential weight, under a radial dominance condition relative to a prescribed profile. This yields an explicit upper bound for the norm of the associated barycenter in terms of a one-dimensional comparison measure. We also provide tractable sufficient conditions for radial dominance, including strong convexity, addition of nonnegative convex terms, and star-shaped constraints. As a consequence, we obtain a refined convergence rate for stochastic proximal estimators of weakly convex functions, together with asymptotic sharpness of the constant. The same framework yields a corresponding rate for stochastic projection estimators onto closed convex sets. We further establish basic structural properties of the barycentric approximation operator, such as smoothness and cocoercivity. Numerical experiments illustrate the predicted rate, the dimensional scaling of the constant, and its asymptotic sharpness.

We investigate a special class of nonlinear infinite dimensional systems. These are obtained by subtracting a nonlinear maximal monotone (possibly multi-valued) operator M from the semigroup generator of a scattering passive linear system. While the linear system may have unbounded linear damping (for instance, boundary damping) which is only densely defined, the nonlinear damping operator M is assumed to be defined on the whole state space. We show that this new class of nonlinear infinite dimensional systems is well-posed and incrementally scattering passive. Our approach uses the theory of maximal monotone operators and the Crandall-Pazy theorem about nonlinear contraction semigroups, which we apply to a Lax-Phillips type nonlinear semigroup that represents the whole system.

We study a family of preconditioned primal dual algorithms for convex-concave saddle point problems by the dynamics introduced in \cite{apidopoulos2026preconditioned}. The proposed framework exploits the possible smooth + nonsmooth structure of the saddle point formulation. It includes, but is not limited to, linearly constrained convex optimization problems. The proposed antisymmetric preconditioners allow us to establish non ergodic convergence rates, accounting for possible computational errors in the implementation of the method. Finally, we present numerical experiments to indicate our well performed preconditioned primal dual algorithms.

This paper proposes a robust optimization formulation to calculate dynamic operating envelopes (DOEs) to safely operate unbalanced three-phase distribution systems. Unlike conventional formulations that satisfy network constraints only at the envelope bound, the robust formulation covers the entire envelope range. We formulate a robust non-linear programming (NLP) problem with the full AC power flow equations, as well as an approximate linear programming (LP) model. Numerical simulations are run with real-world data from Belgium and two different distribution test feeders. The paper compares the conventional approaches with their robust counterparts and examines the trade-off between constraint violation and envelope size as well as accuracy and solve time aspects.

We investigate a special class of nonlinear infinite dimensional systems. These systems are obtained by modifying the second order differential equation that is part of the description of conservative linear systems out of thin air introduced by M. Tucsnak and G. Weiss in 2003. The modified differential equation contains a new nonlinear damping term, that is maximal monotone and possibly set-valued. We show that this new class of nonlinear infinite dimensional systems is incrementally scattering passive (hence well-posed). Our approach uses the theory of maximal monotone operators and the Crandall-Pazy theorem about nonlinear contraction semigroups, which we apply to a Lax-Phillips type nonlinear semigroup that represents the whole system. We illustrate our result on the n-dimensional wave equation.

This article concerns interacting particle systems with singular kernels, driven either by degenerate deterministic controls or by degenerate decomposable noise. In the deterministic setting, we establish global exact controllability and a topologically robust property called solid controllability. Moreover, we prove a result that guarantees global approximate controllability with prescribed trajectories. For stochastic dynamics, we obtain ergodicity and exponential mixing by utilizing coupling and recurrence mechanisms based on controllability. Our approach exploits the singularity and applies to a broad class of models, including Biot-Savart, Coulomb, Riesz, and Yukawa interactions, as well as heterogeneous multi-species systems.

An important quantity in the theory of gradient descent (GD) is the sharpness, defined as the largest eigenvalue of the objective Hessian. Classical analyses typically require the step size to be uniformly smaller than twice the reciprocal of the sharpness, but this condition is frequently violated in the training of deep neural networks. Recent work bridges this gap in the setting of overparametrised least-squares with a single scalar output, providing a normal form for large-step GD in a neighbourhood of an isolated flat minimum and establishing three corresponding convergence results. In this paper, we extend this theory in two directions: (1) to overparametrised least-squares with vector-valued outputs (including regression with arbitrarily many observations), and (2) to a neighbourhood of a manifold of flat minima (which we show is essential for applications such as matrix factorisation). We generalise both the normal form and all three convergence theorems of \cite{macdonaldeos} to this broader setting, overcoming several technical challenges. We further show that our framework applies to deep matrix factorisation under mild assumptions, yielding several new structural results. In particular, we prove that the set of flat minima forms a fibre bundle over a product of spheres, and that the sharpness is Morse-Bott along this manifold.