Publications | Bartolomeo Stellato

Papers listed by type and year of publication or year of submission if unpublished.
We present FOMTools, an open-source Python package for developing, analyzing, and deploying first-order algorithms for parametric linear and convex quadratic optimization. The user specifies an iteration in Python, with syntax close to its mathematical description. The package represents the iteration as a computation graph of affine and piecewise affine operations. To analyze worst-case performance, FOMTools formulates a mixed-integer verification problem that maximizes a convergence residual over a specified set of problem instances and initial iterates. For probabilistic performance analysis, a batched GPU simulator collects algorithm trajectories and computes statistical bounds on the probability of exceeding a specified convergence tolerance. From the same algorithm description, FOMTools generates specialized solver code for embedded hardware, multicore CPUs, and GPUs. The package first generates the numerical data, factorizations, and matrix-vector routines without any agent. A coding agent then translates the iteration into C or CUDA code and iteratively optimizes the implementation. After each phase, the package runs the generated code on sampled inputs and compares it with a reference execution. FOMTools also serves as a CVXPYgen backend, generating and optimizing solvers for parametric CVXPY problems. We demonstrate FOMTools on sparse coding, network optimization, optimal control, and portfolio rebalancing. The verification bounds reveal worst-case behavior that sampling misses. The generated solvers run up to two orders of magnitude faster than NumPy implementations, and they outperform hand-written GPU code and the OSQP solvers generated by CVXPYgen.
We study Nesterov's accelerated gradient method with its standard constant parameters for smooth strongly convex optimization. We characterize the exact worst-case asymptotic factor for the squared distance to the minimizer. To this end, we construct candidate worst-case trajectories by combining planar spirals with a common contraction factor that satisfy a self-similarity condition. Each component rotates by a fixed angle at every iteration. The exact factor is the optimal value of a linear optimization problem over probability measures on a fixed disk, with constraints derived from smooth strongly convex interpolation. Every optimal measure has a common contraction factor, which is also the largest factor for which the angular interpolation conditions are feasible. We also show that this worst-case asymptotic rate can be realized by a \emph{single} function in a finite-dimensional space whose iterates attain this factor. A matching upper bound holds uniformly over all functions in the class and all finite dimensions. Consequently, allowing the function, dimension, and initial point to vary with the iteration horizon does not increase the worst-case asymptotic factor. We also show that the exact factor is 1−4/(3κ)+o(κ−1/2)1-4/(3\sqrt\kappa)+o(\kappa^{-1/2})1−4/(3κ​)+o(κ−1/2) as κ→∞\kappa\to\inftyκ→∞, where κ=L/μ\kappa=L/\muκ=L/μ is the ratio of the gradient Lipschitz constant to the strong convexity parameter.
We consider the problem of generating a large collection of initial guesses for local minima of multimodal non-convex continuous optimization problems. The goal is for these initial guesses to be high-quality (i.e., a numerical solver converges quickly) and diverse (i.e., represent many different local minima). Identifying multiple locally optimal solutions enables flexible downstream decision-making, but typically requires expensive global search. Existing data-driven methods predict initial guesses using only the final converged optima from offline solver runs, which discards information about the local neighborhoods of solutions and limits the available training data. We propose GLENS (Global Search via Learning from Solver Iterates), a data-efficient global search method that leverages intermediate solver iterates as free data augmentation. GLENS consists of two components: a neighborhood structure model that uses diffusion models to learn the local geometry around optima conditioned on problem parameters, and a solver behavior model that learns refinement directions to further guide samples towards nearby optima during diffusion sampling. Experiments on modified non-convex benchmark problems and a two-robot obstacle-avoidance navigation problem show that GLENS generates high-quality initial guesses while preserving the multimodal distribution of diverse local optima. The resulting initial guesses lead to faster solver convergence across different problem settings and solvers. We also analyze how key hyperparameter choices affect the performance.
We introduce the concept of disjunctive sum of squares for certifying nonnegativity of polynomials. Unlike the popular sum of squares approach where nonnegativity is certified by a single algebraic identity, the disjunctive sum of squares approach certifies nonnegativity with multiple algebraic identities which can be found in parallel. Our main result is a disjunctive Positivstellensatz proving that we can keep the degree of each algebraic identity as low as the degree of the polynomial whose nonnegativity is in question. Based on this result, we construct a semidefinite programming based converging hierarchy of lower bounds for the problem of minimizing a polynomial over a compact basic semialgebraic set, where the size of the largest semidefinite constraint is fixed throughout the hierarchy. We further prove a second disjunctive Positivstellensatz which leads to an optimization-free hierarchy for polynomial optimization. We specialize this result to the problem of proving copositivity of matrices. Finally, we describe how the disjunctive sum of squares approach can be combined with a branch-and-bound algorithm and we present numerical experiments on polynomial, copositive, and combinatorial