Papers

Affine-invariant ensemble samplers are widely used in Bayesian applications. However, their quantitative convergence theory, in particular geometric ergodicity, remains a basic open question. We study the affine invariant ensemble Langevin dynamics, an interacting particle system that uses the empirical covariance of the whole ensemble as a preconditioner. While effective in practice, theoretical understanding of this method is not available beyond plain qualitative convergence in total variation, with a central difficulty being that the empirical covariance can become degenerate. For potentials with bounded Hessian that are strongly convex outside a ball, we prove geometric ergodicity using a novel Lyapunov function that combines the inverse covariance with a coercive exponential energy. We then show that directly applying the Euler–Maruyama scheme can diverge with positive probability, even for a one-dimensional Gaussian target. This motivates a covariance-trace regularization in continuous time. We prove geometric ergodicity of the regularized diffusion and, for sufficiently small step size, of its unadjusted Euler–Maruyama discretization. We also show that the invariant distributions of the discretization converge weakly to the product target distribution as the step size tends to zero. Lastly, we provide an explicit discretization that is both affine-invariant and geometrically ergodic, based on an adaptive step size clamping.

We provide two complementary explanations of H-duality in smooth strongly convex optimization and contractive fixed-point problems. H-duality refers to the phenomenon where the worst-case performance of many fixed-step first-order methods (FSFOMs) for a given performance setup are exactly equal to the worst-case performance of their "time-reversed" counterparts and a paired performance setup. This concept gained interest as a mechanism for automatically converting an algorithm with guarantees on the final suboptimality to algorithms with guarantees on the final gradient norm. Since its initial discovery, additional H-dual performance setup pairings have been discovered empirically, yet a general explanation for these phenomena remained elusive. Our first explanation of H-duality shows that, on a class of instances with extremal curvature properties, FSFOMs produce a final iterate that is a scalar multiple of the initial iterate, and that this scalar is invariant under the time reversal of the underlying FSFOM. This, together with the empirical observation that many FSFOMs have such an instance with extremal curvature as their worst case, gives one explanation for why H-duality often occurs. The second explanation is complementary and uses the notion of a performance estimation certificate of the convergence rate of a FSFOM. We exhibit an explicit transformation of such certificates that preserves many of the properties that define such a certificate, and we show that in many numerical and analytical examples, this transformation in fact transforms certificates of the convergence rate of an algorithm into certificates for its time-reversal. As an application, we prove a convergence guarantee for the H-dual of ITEM matching the conjectured optimal gradient-to-gradient guarantee.

The dissertation considers various generalizations of the classical risk process describing the stochastic evolution of the capital of an insurance company (Cramer Lundberg model), including processes with variable deterministic or random premiums, non Poisson flows of premiums and claims, and a stochastic Markovian environment. Integral equations for the probability of nonruin as a function of initial capital are derived for these generalizations. For processes in a stochastic Markovian environment, systems of integral equations for nonruin probabilities corresponding to different initial states are obtained. General necessary and sufficient, as well as specific sufficient, conditions for the existence and uniqueness of solutions are established. A successive approximation method for numerical and analytical solution of the integral equations is theoretically and practically validated; its uniform convergence and rate of convergence are established. A technique for estimating the accuracy of approximate solutions is developed by constructing upper and lower approximations to the exact solution. The method is tested on numerical examples and compared with the Monte Carlo method and known solution approximations. The developed method increases the accuracy of actuarial calculations: it allows the probability of ruin to be calculated with any prescribed accuracy, empirical approximations to be verified and improved iteratively, and the accuracy and parameters of Monte Carlo simulations to be estimated and corrected.

We study a mean-field optimal stopping problem with an endogenous population-level shutdown. All remaining agents stop when the survival mass falls below a prescribed threshold. We recast the discontinuous objective as the singular, nonconvex constraint that the survival mass lie in $\{0\}\cup[\alpha,1]$. We prove the equivalence of strong and weak values via an approximation and the existence of an optimal rule via compactness and penalization. We also prove a dynamic programming principle. The value is continuous away from the critical boundary but may be discontinuous at the boundary itself. Under strict initial feasibility, finite-population values converge to the mean-field value. In the same regime, the laws of near-optimal empirical measures are tight and every mean-field optimizer admits a recovery sequence. At the threshold, however, finite-population convergence may fail.

We introduce an anytime, almost parameter-free algorithm for convex function-constrained optimization, in which the objective and constraint functions are unknown piecewise-smooth. Our algorithm, Restarted Penalty APEX (Penalty Accelerated Prox-level method for Exploring Piecewise Smoothness), is an accelerated bundle-level method based on the penalty approach. For problems satisfying quadratic growth, the proposed method is the first to achieve optimal oracle complexity without knowing the growth modulus or the exact penalty coefficient. Furthermore, Restarted Penalty APEX generates a verifiable certificate that bounds both the optimality gap and the constraint violation. This certificate also appears to be new to the literature.

We study an inverse optimal-control problem for identifying nonnegative weights in a convex-feature lower-level objective. The resulting model is an optimistic bilevel optimal-control problem whose lower level is a strongly convex, linearly constrained control problem with box constraints in \(L^2\). For the reduced response map \(x\mapsto (y_x,u_x)\), we prove global single-valuedness and an explicit Lipschitz estimate, without assuming differentiability of active sets. Exploiting the affine dependence of the lower-level functional on \(x\), we show that the value function is concave and Lipschitz and derive an asymptotically exact optimal-value relaxation. We further establish local identifiability and \(O(\delta)\) noise stability by a gradient-free argument based on second-order growth. Finally, a KKT reformulation yields a function-space MPCC; Scholtes relaxation gives C-stationarity under uniform multiplier boundedness. A decoupled scalar example provides an explicit strict-complementarity certificate via finiteness of threshold contacts.

J. Danskin proved in 1966 that the envelope of a family of continuous functions $\mathrm F : \mathbb R^n \times \mathrm K \longrightarrow \mathbb R$, parameterised by the points of a compact $\mathrm K \subseteq \mathbb R^m$, given by \[ f(u) := \max_{k \in \mathrm K} \mathrm F(u,k), \ \ \ \ f \, : \, \mathbb R^n \longrightarrow \mathbb R, \] is semi-differentiable on $\mathbb R^n$, if $\mathrm F$ is sufficiently regular. This result is of utmost importance in numerous applications. Extensions have been proved towards almost every direction, but none permits to replace “max over $\mathrm K$" with “essential sup over $\mathrm K$". We establish the semi-differentiability of generalised $\mathrm L^\infty$ envelopes when $ \mathrm F$ is defined on the product of a Banach space with a measure space. As an application, we obtain a new foundational regularity result in the Calculus of Variations in $\mathrm L^\infty$, asserting that supremal functionals are semi-differentiable everywhere, with an explicit formula for their semi-derivative, despite generally being non-differentiable. Moreover, we establish a variational characterisation of (absolute) minimisers of general level-convex $\mathrm L^\infty$ functionals via their semi-differentials, revealing the genuine $\mathrm L^\infty$ counterpart of the Euler-Lagrange equations. Conventional Aronsson equations and related divergence PDE involving measures, deducible from $\mathrm L^p$ approximations as $p\to \infty$, in general are not sufficient for minimality in $\mathrm L^\infty$. Our results dislodge the trademark necessity for extrinsic $\mathrm L^p$-approximations as $p\to \infty$ to derive and study PDEs, reinterpreting the existing $\mathrm L^\infty$ theory. Most crucially, they provide a new powerful intrinsic $\mathrm L^\infty$ framework, evocative of the straightforward variational methods for integral functionals.

The Job Sequencing and Tool Switching Problem (SSP) is a well-known combinatorial optimization problem arising in the context of flexible manufacturing. Since the seminal work of Tang and Denardo (1988), the SSP has received significant attention in the literature, leading to the development of numerous exact and heuristic approaches. Despite these efforts, several benchmark instances proposed decades ago and containing only 20 jobs have remained unsolved to proven optimality. In this work, we propose an exact algorithm for the SSP, namely the Combinatorial Branch-and-Bound (C-B&B) algorithm, which combines two distinct branch-and-bound algorithms, each introducing novel features compared with the existing literature. The former relies on a new branching scheme designed to reduce the size of the implicit enumeration tree, together with a collection of new bounding functions. The latter builds on the branching scheme introduced by Laporte et al. (2004) and strengthens it with a new bounding function and two dominance rules. Within C-B&B, these exact algorithms are complemented by a preprocessing phase that incorporates a new branch-and-bound-based heuristic capable of rapidly generating a high-quality initial incumbent solution. Extensive computational experiments show that C-B&B represents a strong breakthrough over previously published approaches, proving optimality for more instances with significantly less computational effort and closing several benchmark instances that have remained open for decades.

We consider the task of optimizing smooth, possibly non-convex functions of a matrix variable given access only to directional derivatives rather than full gradients. This setting arises when fine-tuning large neural networks on consumer-grade hardware: the network's weights are matrices, memory constraints rule out backward-mode automatic differentiation, but directional derivatives remain available through forward mode. We frame gradient estimation in this setting as a structured recovery problem, in the spirit of signal processing. From this perspective we provide three contributions. First, we introduce an alternative to the standard random gradient estimator; the difference corresponds to replacing the adjoint of the sampling operator with its pseudoinverse. Second, when the gradient satisfies an approximate low-rank condition, techniques from matrix sensing yield a family of highly accurate gradient estimators that drop into any first-order method. Third, we note that while the computational cost of such estimators is high, this can be amortized by combining them with a matrix-aware optimizer such as spectral descent. Specifically, the gradient estimator computes a factorization of the gradient, allowing for the projection step of spectral descent to be done at no extra cost. We demonstrate our findings with two careful numerical experiments on synthetic functions with approximately low-rank gradients. We show that by exploiting this low-rank property one obtains much faster convergence to good approximate solutions.

Energy-based analysis of the three-dimensional incompressible Navier–Stokes equations with spatially localized dissipation is developed at the level of kinetic energy and Leray–Hopf weak solutions. A localized coercivity property linking viscosity and damping to global $L^2$ control is introduced and related to observability and spectral inequalities for the associated Stokes operator. Under this condition, exponential energy decay is established without regularity or linearization assumptions, and time-dependent damping is treated within a nonautonomous framework. A finite-dimensional spectral quantity is proposed as a computable proxy for coercivity, yielding explicit decay estimates for reduced-order models.

Modern feedback design for distributed parameter systems presupposes that the closed-loop dynamics define a well-posed evolution problem. This presupposition becomes nontrivial for temporally degenerate parabolic systems, where temporal degeneracy affects not only the analytical properties of the evolution equation but also the mathematical formulation of the feedback interconnection itself. It is shown that admissible feedback interconnections for temporally degenerate parabolic systems are completely characterized by an operator compatibility condition linking the singular reaction operator with the actuator and observation operators. This characterization removes the singular component of the closed-loop dynamics and reduces the degenerate evolution equation to a regular evolution equation. Building upon this regularized formulation, a critical–residual decomposition yields a uniform exponential stability certificate, which is subsequently extended to the original infinite-dimensional evolution through a finite-to-infinite lifting theorem. A constructive static output feedback synthesis is finally obtained as a consequence of these results. Numerical experiments illustrate the regularization mechanism, validate the stability certificate, and confirm the finite-to-infinite lifting.

The Gromov-Wasserstein (GW) distance provides a principled framework for aligning metric measure (mm) spaces based solely on their intrinsic structure. Its ability to identify isomorphic representations of distributions across spaces renders it valuable for comparing data where equality up to isomorphism occurs naturally such as in graphs or, more generally, distributions on graphs. Recently, a type of dual form for the GW distance between Euclidean distributions with the squared Euclidean or inner product costs was derived, spurring the development of new statistical and algorithmic results for this setting. This work furnishes a novel duality result for GW distances with and without entropic regularization that is applicable to all finitely supported mm spaces. Leveraging this result, we derive the sample complexity of empirical GW distances between finite mm spaces, as well as limit distributions under proper centering and scaling. Furthermore, we propose new algorithms for solving the regularized GW problem which are subject to formal convergence guarantees. These statistical and algorithmic advancements give rise to a principled and efficient framework for testing whether two distributions on the set of graphs with a fixed number of nodes are isomorphic based on samples.

Gaussian-process Bayesian optimization (GP-BO) excels at black-box optimization of costly functions, e.g., hyperparameter optimization (HPO) and multi-agent system (MAS) design. Convergence-rate guarantees exist for select methods, notably GP upper confidence bound (GP-UCB), but require a fixed kernel. Critically, the kernel encodes how input proximity affects objective value similarity. When raw coordinates poorly match this geometry - as with log-scaled hyperparameters or localized peaks - input warping can greatly improve sample efficiency, yet known GP-UCB proofs require a fixed kernel. We propose Finite-Library Input-Warped Bayesian Optimization (FLIWBO), which selects warps from a finite library of smooth input maps by any history-dependent rule. It adapts the input geometry to accelerate learning while retaining high-probability convergence guarantees under mild hypotheses, with an explicit $\sqrt(N_\varepsilon)$ library-size cost. Controlled diagnostics show that finite-library warping repairs planted geometry mismatches and identify FLIWBO failure cases. Across four repeated benchmarks - warped synthetic objectives, a confidence-fence trap, and Fashion-MNIST HPO - FLIWBO-UCB beats raw-coordinate GP-UCB under misspecified geometry, escapes traps that defeat even oracle-warp expected improvement, and recovers much of the gain from manual log scaling, while leading the tested methods that admit a matching regret guarantee. A 20-dimensional MAS design study further shows feasibility under costly noisy evaluations. Code for experiments is available: https://github.com/edvin-ketabati/bogp-paper-experiments.

We study Gaussian-restricted barycenters for quadratic two-sided Kullback–Leibler unbalanced optimal transport with independent marginal penalties and no coupling entropy. Exact profiling of the barycenter mass reduces the problem to a smooth Gaussian shape functional with endogenous Gibbs weights. We establish global attainment, derive the stationary moment equations, and construct a reverse-KL majorization–minimization (MM) iteration whose full sequence converges from every nondegenerate Gaussian initialization to a stationary fixed point. The diagonal second variation induces a parallel-sum tensor coupling the Bures–Wasserstein and Fisher–Rao metrics; its finite-mass extension admits a radial cone representation. Under common penalty scaling, global minimizers converge to a Gaussian Wasserstein barycenter with effective weights; for sufficiently large penalties, the minimizer is unique and admits a first-order analytic expansion. In the small-penalty regime, the distance of every global minimizer to the compact maximizer set of a weighted Chernoff affinity functional vanishes with respect to the mean–covariance parameter distance. Numerical experiments illustrate MM descent, local contraction, and the two penalty limits.

We study how far gradient descent (GD) can be accelerated by predetermined stepsizes in smooth convex optimization. Going beyond the classical $\Omega(n^{-2})$ first-order oracle lower bound of Nemirovsky and Yudin (1983), we prove an $\Omega(n^{-1.6342})$ non-anytime lower bound and an $\Omega(n^{-1.2408})$ anytime lower bound. These improve the recent $\Omega(n^{-1.932})$ non-anytime lower bound of Ma and Chen (2026) and the $\Omega(n^{-4/3})$ anytime lower bound of Tsai et al. (2026), respectively. Both results continue to hold when the stepsizes may be negative. Our anytime lower bound also shows that the $O(n^{-\log_2(1+\sqrt{2})})$ rate of non-anytime silver schedules (Altschuler and Parrilo, 2025; Grimmer et al., 2025) is unattainable in the anytime setting. This establishes a strict separation between the two settings.

Data-consistent inversion (DCI) constructs probability measures whose push-forward distributions agree with observed data, while iterative data-consistent inversion (iDCI) extends this framework to generalized stochastic inverse problems by enforcing multiple push-forward constraints sequentially. Although iDCI avoids the direct approximation of high-dimensional joint densities, its relationship to the original joint DCI solution has remained unclear. In this work, we establish this relationship through copula theory. Using Sklar's theorem, we derive a factorization of the DCI update into separate marginal and dependence transformations and show that the discrepancy remaining after convergence of the iDCI algorithm is entirely characterized by the copulas associated with the observed and predicted joint distributions. This characterization motivates a copula-transformed iDCI solution, and we prove that an exact copula transformation recovers the original DCI solution. We further establish convergence results for approximate copula transformations under converging sequences of reference measures and progressively enriched feasible sets. Numerical examples demonstrate how the geometry induced by the quantity-of-interest map governs the importance of the copula transformation, illustrate an adaptive reference-measure refinement strategy for improving computational accuracy under a fixed sampling budget, and demonstrate the progressive refinement of generalized stochastic inverse problems through heterogeneous, asynchronously acquired experiments.

A basic question in adaptive control is how rapidly a discrete-time nonlinear plant with unknown parameters may grow while remaining stabilizable. When the nonlinear growth exponent is known and only a scalar coefficient is unknown, existing theory identifies 4 as the exact critical exponent for stabilizability. We determine how this critical exponent changes when the growth exponent is also unknown. For any positive disturbance bound, one feedback law stabilizes all plants in some neighborhood of a nominal parameter pair if and only if the nominal exponent is below $3\sqrt{3}/2$. At equality, every compact parameter set whose exponents do not exceed this value remains stabilizable. If the exponent belongs to a known finite set, the critical exponent remains 4 for every nondegenerate compact coefficient interval. Hence finite exponent sets and arbitrarily short exponent intervals can have different critical exponents. Uncertainty in the growth rate therefore changes the range of growth that feedback can stabilize, not merely the size of the uncertainty set. The critical exponent remains $3\sqrt{3}/2$ under a known positive state-dependent multiplier bounded above and away from zero. For a system whose coefficient and exponent are known functions of an unknown parameter, a compact parameter family is stabilizable when its exponents do not exceed this value. At an interior parameter point where the Jacobian of these functions has rank two and the exponent is at least this value, no compact neighborhood is stabilizable.

We study when and how momentum improves large-batch training in the one-pass regime, using power-law kernel regression as a tractable setting. We first characterize risk stability through the critical learning rate, defined as the largest learning rate for stable training, and obtain $\eta_{\mathrm{SGD}}^{\mathrm{crit}}\eqsim 1$, $\eta_{\mathrm{Polyak}}^{\mathrm{crit}}\eqsim \min\{1,B(1-\rho)\}$, and $\eta_{\mathrm{Nesterov}}^{\mathrm{crit}}\eqsim \min\{1,B^\beta(1-\rho)\}$, where $B$ is the batch size, $\rho$ is the momentum factor, and $\beta>1$ is the capacity exponent. Within this admissible region, we derive scaling laws for the full risk dynamics, capturing the progression from an early transient, through power-law decay, to a noise floor. We then minimize the final-step risk over the admissible learning rates and momentum factors under a fixed data budget, yielding a three-regime batch-size phase diagram that reveals how the role of momentum changes with batch size. Notably, Polyak enlarges the critical batch size, the largest batch size preserving the best small-batch data-scaling exponent, thereby enabling greater parallelism without sacrificing data efficiency. In contrast, Nesterov achieves better data efficiency in the large-batch regime because its look-ahead mechanism suppresses noise accumulation. Numerical experiments validate the predicted stability boundaries, risk dynamics, and batch-size phase diagram.

We consider stochastic multi-objective optimization over a nonempty closed convex set, where every objective is an expectation and only sample-gradient information is available. We develop a line-search-free and function-value-free adaptive projected-gradient algorithm for the sample-average approximation (SAA) problem. Each iteration computes a feasible regularized multi-gradient step and updates the regularization parameter from the projected step length. A normal-cone-based certificate yields descent estimates and an explicit complexity bound for the Pareto-stationarity residual of the SAA problem. The consistency of SAA gradients then transfers vanishing SAA residuals to Pareto stationarity for the population problem, while an additional concentration argument gives a finite-sample residual bound on compact sets. Experiments on synthetic problems, classification, portfolio selection, multi-task learning, and robot control illustrate the practical performance of our algorithm.

We study insensitizing control problems for a class of coupled linear stochastic parabolic systems. We establish the existence of controls such that a sentinel functional, involving localized observations of the state variables and their spatial gradients, is insensitive to small perturbations of the null initial data. We first reformulate the insensitizing control problem as a null controllability problem for a coupled forward–backward stochastic parabolic system, in which the observation terms induce both zeroth- and second-order coupling terms. By duality, the analysis is reduced to an observability inequality for the corresponding adjoint system. The main analytical contribution is the derivation of new global Carleman estimates for coupled stochastic parabolic systems with zeroth- and second-order coupling terms, under suitable geometric assumptions on the control and observation regions. These estimates yield the required observability inequalities and, consequently, the existence of insensitizing controls. Furthermore, depending on the value of a weighting parameter $\beta\in[0,1]$, which determines the relative contributions of the two state components to the sentinel functional, we consider two cases. If $\beta\in\{0,1\}$, the sentinel functional depends on only one state component, and a single localized control acting in the drift of the first equation is sufficient. In contrast, if $\beta\in(0,1)$, both state components contribute to the sentinel functional, and two localized controls acting in the drift terms of the two equations are sufficient. Moreover, the control strategy in this paper involves two additional controls acting throughout the diffusion terms.

The Boltzmann equation plays an important role in modeling mesoscopic behavior in a wide range of scientific and engineering applications. However, its numerical solution is computationally expensive due to the high dimensionality of the model and the nonlinear nonlocal collision operator, especially for steady-state problems that require iterative solvers. This cost becomes prohibitive for inverse problems, where the induced optimization problem requires repeated forward solves. In this work, we propose a reduced-order model (ROM) for the parametric Boltzmann equation to address this computational challenge. The ROM constructs a low-dimensional approximation space for the parameter-induced solution manifold through a residual-based greedy strategy, and the reduced solution is then obtained via residual minimization over the reduced space, subject to mass conservation. The overall efficiency of the ROM is achieved by exploiting the quadratic structure of the collision operator and a precomputed separable approximation of the collision kernel. The resulting ROM is further applied to a thermally-driven inverse problem for reconstructing collision parameters from the observed macroscopic temperature data. This is accomplished either by directly replacing the PDE constraint with the ROM, leading to a bilevel optimization formulation, or by reformulating the task as a single-level optimization problem through the Karush–Kuhn–Tucker (KKT) conditions. Numerical experiments in both collision-dominated and transport-dominated cases are performed to demonstrate the efficiency and accuracy of the proposed ROM and its effectiveness in inverse problems. In particular, the resulting inverse problem is computationally much more tractable, achieving speedups of several orders of magnitude over that based on the full-order model while maintaining comparable accuracy.

We present a combinatorial optimization problem for the strategic deployment of UAVs equipped with 5G antennas to assist rescue operations in regions hit by natural disasters. Our goal is to optimize the placement of UAVs to provide coverage in flying ad-hoc networks among given candidate sites. Our formulation aims to maximize signal coverage and minimize interference while ensuring network connectivity. To mitigate interference effects, we incorporate the use of multiple frequencies. We formulate this problem as an integer quadratic program (IQP). We present numerical solutions obtained via the CPLEX solver and conduct a preliminary analysis of the problem's scalability in realistic network configurations. Our findings reveal a significant exponential increase in Time-to-Solution (TTS) as the number of sites grows, which poses a critical challenge in urgent, time-sensitive scenarios. To address this issue, approximate suboptimal solutions can be produced by enforcing a time limit on the solver. Although these solutions are not optimal, they preserve connectivity in most cases, providing a practical trade-off between solution quality and computational times that remain within feasible limits for real-time UAV redeployment. Recognizing the limitations of classical solvers in these contexts, we explore quantum computing as a promising alternative. Specifically, we reformulate the problem as a quadratic unconstrained binary optimization (QUBO) problem, suitable for most quantum algorithms. Through high-performance computing emulation, we show that the quantum adiabatic algorithm (QAA) can accurately solve small-scale instances, paving the way for future application of quantum computing to large-scale, time-critical optimization problems in disaster response.

We investigate the optimization problem of minimizing a nonsmooth function that satisfies a nonsmooth version of the descent lemma over a nonempty and closed but not necessarily convex set. The objective function belongs to the class of upper-$\mathcal{C}^2$ functions, whereas the constraints may promote a sparse or low-rank structure. We propose a projected subgradient method with two different globalization strategies: (a) a nonmonotone linesearch and, under additional assumptions, (b) an auto-conditioned method, where the stepsize is given by a formula depending on data from past iterations. We show that both methods converge to solutions that satisfy a stronger stationarity concept than one would expect from the subdifferential sum-rule, which is particularly important since the optimization problems of interest are inherently nonconvex. Finally, we present promising numerical results when applying the algorithm to an MPEC-style problem as well as the matrix optimization problems MAXCUT and Robust PCA.

We study the problem of learning from data a parametric Nash equilibrium (NE) problem that is monotone (or strongly monotone) for all parameter values. In the presence of local and shared convex constraints, monotonicity enables efficient computation of generalized Nash equilibria of the learned game. We consider two learning scenarios: (i) direct learning of the NE problem from samples of the agents' costs, and (ii) inverse learning of surrogate agents' costs from samples of their best responses. We propose two methods to solve these tasks. The first is a penalty-based approach that promotes monotonicity of the learned game during training. The second, based on a representation theorem we introduce for a broad class of monotone games, parameterizes the agents' costs so that monotonicity of the NE problem is guaranteed by construction, for all parameter values, regardless of the training data used. We illustrate the applicability of the proposed methods on several numerical examples. A Python library and the examples reported in the paper are available at https://github.com/bemporad/learn_monotone_games.

This paper introduces a novel second-order dynamical system driven by a forward-backward splitting operator for solving the structured inclusion $0\in(A+B)(x)$ in a real Hilbert space, where $A$ is a maximal $\rho$-comonotone operator and $B$ is a $\nu$-cocoercive operator. The well-posedness of the system is established, and Lyapunov analysis yields accelerated convergence rates of order $o(\frac{1}{t})$ for the velocity and $o(\frac{1}{t^2})$ for the forward-backward residual, together with weak convergence of the trajectories to $\operatorname{zer}(A+B)$. Temporal discretization further leads to a class of double inertial Halpern forward-backward splitting algorithms that encompasses the classical forward-backward splitting method and its inertial variants as special cases. Numerical experiments on split feasibility, sparse signal recovery, and image deblurring illustrate the effectiveness of the proposed algorithm.

In this paper, we study the emergent collective dynamics of multi-agent systems driven by infinite distributed fading memory of Volterra type. We establish a unified theoretical framework covering both first-order opinion consensus dynamics and second-order velocity alignment flocking kinematics. By introducing Dafermos past-history transformations, the governing integro-differential systems are reformulated into dynamical systems on an extended product Hilbert spaces. For first-order dynamics, we prove that fading memory inherently provides a hidden dissipative mechanism, guaranteeing unconditional global exponential consensus with or without instantaneous communication forces. For second-order dynamics, we obtain unconditional exponential flocking for the pure fading memory system and we give a sufficient condition for flocking when an instantaneous interaction is also present. In particular, this condition is always satisfied, namely the flocking occurs unconditionally, when the influence function has a divergent tail.

Feature and sample screening reduce the cost of machine learning by eliminating irrelevant features and noninformative samples, respectively. Although recognized as primal-dual counterparts, their relationship remains informal and model-dependent. Viewing screening and duality as transformations of objective functions, we introduce Fenchel-Rockafellar (FR) representations, a class of convex problems encompassing the Lasso and SVM that is closed under both transformations. We then prove that feature and sample screening form an equivariant pair: dualization followed by feature screening is equal to sample screening followed by dualization.

This paper shows that the optimal gain of the risk-sensitive linear quadratic regulator (LQR) problem is invariant under input randomization, i.e., when the controller deliberately injects noise into the nominal control input. This appears counterintuitive at first glance because certainty equivalence does not hold for risk-sensitive LQR and input randomization inflates the effective process noise. Nonetheless, the gain is preserved because the input noise enters not only the system dynamics but also the cost functional, and its total effect on the gain eventually vanishes. Consequently, the optimal gain and its associated Riccati recursion need not be recomputed, and the increment in the optimal cost can be readily evaluated in closed form. This result facilitates the use of risk-sensitive LQR in applications that employ input randomization for privacy or exploration, such as watermarking for replay attack detection, differential privacy, and path integral control.

In this lecture note, we present a concise and self-contained derivation of the discrete-time Kalman filter equations that requires only a basic understanding of least squares estimation. The treatment is designed to minimize mathematical overhead while preserving both rigor and generality.

The cycle space of a directed graph is generated by a cycle basis, where, in general, cycles are allowed to have both forward and backward arcs. In a forward cycle, all arcs must follow the given direction. Several open questions remain regarding the complexity of the minimum cycle basis problem, in particular the minimum-weight integral cycle basis problem, and the minimum-weight weakly and strictly fundamental forward cycle basis problems. In this paper, we address these open questions. First, we study the existence, structure, and computational complexity of minimum-weight forward cycle bases. We give a complete structural characterization of digraphs that admit weakly fundamental (and hence integral) forward cycle bases. We further provide a characterization when a strongly connected digraph admits a forward fundamental cycle basis, proving that such a basis exists if and only if the set of directed cycles has cardinality equal to the cycle rank; in this case, the basis is unique. Lastly, we show that while minimum-weight forward fundamental cycle bases can be found in polynomial time whenever they exist, the minimum-weight forward weakly fundamental cycle basis problem is APX-hard via an L-reduction from the minimum-weight weakly fundamental cycle basis problem on digraphs with metric weights. Second, we introduce opt-in graphs, i.e., the family of graphs for which minimum cycle bases are integral for any weight function. We show that this family is minor-closed and hence, by the Robertson-Seymour theorem, is characterized by a finite set of forbidden minors, so that the opt-in recognition problem is solvable in polynomial time. Lastly, we present an algorithm to check whether a graph is opt-in, and if not, to identify which of its minors belong to the set of forbidden minors. Applying this algorithm, we show that the complete graph $K_n$ is opt-in if and only if $n \leq 7$.