optimization engineers

Abstract: Automatic operator design (AOD) has increasingly adopted large language models (LLMs) to generate problem-specific search and variation operators for multi-objective combinatorial optimization (MOCO). Yet most methods evolve the operators while leaving the process that discovers them fixed. This limits adaptation as new Pareto trade-offs and operator interactions change where and how the search should proceed. We introduce HEMO, which makes the AOD harness itself a self-evolving optimization target. HEMO maintains a persistent harness state and updates one module at a time, revising the logic that determines which experience to retrieve, which operators to improve, which parents to select, and how to formulate the orchestration of search steps. Downstream solution quality and runtime guide these updates, while a discounted UCB scheduler directs adaptation toward promising modules as their utility changes. A persistent memory filesystem preserves candidate outcomes and diagnostic feedback across harness states, enabling successive states to build on accumulated optimization experience. Experiments across eight MOCO tasks and two structurally distinct evolutionary backbones demonstrate strong performance against conventional evolutionary, neural, and LLM-based AOD baselines, together with effective transfer to related tasks. HEMO extends automatic operator design beyond evolving individual operators to improving the search process that discovers them.

Abstract: In the study of constraint satisfaction problems (CSPs), the non-redundancy of a relation measures the maximum number of clauses in a CSP instance which can be logically independent. In the last few years, CSP non-redundancy has found a number of applications including in query complexity, sparsification, kernelization and streaming. In this paper, we present a novel technique for upper-bounding the non-redundancy of some relations using locally decodable codes (LDCs) and locally correctable codes (LCCs). More precisely, known limitations on the structure of 2-query LDCs and LCCs can be translated into limitations on the size of certain non-redundant CSP instances. This technique builds on a technique used by Brakensiek and Guruswami [STOC 2025] to bound non-redundancy via techniques from matching vector codes. Using LDCs and LCCs, we give improved upper bounds on non-redundancy of the BCK predicate pioneered by Bessiere, Carbonnel, and Katsirelos [AAAI 2020] as well as its well-studied Boolean variant. We also give the first example of a Boolean predicate whose non-redundancy is provably not of the form $n^{k\pm o(1)}$ for some positive integer $k$.

Abstract: L-extendable functions on the integer lattice are functions that admit a $\frac12$UJ-convex relaxation; that is, values can be assigned to the additional half-integral points, while keeping the original values at integer points, so that the resulting function is $\frac12$UJ-convex. Various examples of L-extendable functions and their applications are known. A general characterization of L-extendability for functions on the integer lattice, however, remains open. In this paper, we study the L-extendability of functions on the integer lattice. We first obtain a necessary condition: L-extendable functions are integrally convex. We then consider a special subclass of L-extendable functions, namely, those whose local convex extensions give L-convex relaxations. This subclass has a useful optimization property: the optimal value of the original function coincides with that of the $\frac12$UJ-convex relaxation given by its local convex extension, a property not guaranteed for general $\frac12$UJ-convex relaxations. We characterize this subclass by showing that a function is L-extendable through its local convex extension if and only if it is integrally convex and its restriction to every unit hypercube is D-convex, where D-convexity is a discrete convexity associated with reflected Freudenthal simplicial divisions. This subclass strictly contains the class of D-convex functions. We also give a combinatorial characterization of D-convex subsets of the vertices of a unit hypercube: a set $X\subseteq\{0,1\}^n$ is D-convex if and only if it is majority-closed and contains a relative opposite pair.

Abstract: Large language models (LLMs) have emerged as potential tools for numerical optimization. Despite promising results, their intrinsic optimization capabilities and limitations remain poorly understood, particularly in high dimensions. This paper investigates the ability of LLMs to exploit mathematical problem structures and generate feasible, near-optimal solutions. We develop an evaluation framework that distinguishes solution quality from execution failures, including numerical inaccuracies, constraint violations, and output-format errors. Using representative problems from convex optimization, fractional programming, and constrained projection, we assess LLM performance against deterministic reference solvers across different dimensions. Our results reveal that LLMs can capture useful optimization structures and produce near-optimal solutions, even when their outputs fail to satisfy numerical or formatting requirements. However, reliability deteriorates with increasing problem dimension and depends strongly on the underlying optimization structure. To further investigate the limitations, we deploy an iterative refinement mechanism that provides targeted feedback on feasibility and solution quality. The results demonstrate that iterative refinement substantially improves reliability and extends the range of problem dimensions over which LLMs produce valid solutions. Nevertheless, persistent numerical and structural errors limit their applicability as standalone optimizers. These findings provide insights into the capacity and limits of LLM-based optimization and motivate hybrid approaches that combine LLM reasoning with deterministic numerical verification.

Abstract: Reflective text optimization makes agentic systems better, but the optimization loop is itself a resource-hungry workload whose sequential structure leaves parallel evaluation capacity idle. We extend the open-source GEPA optimizer to propose and evaluate a batch of candidates on each optimization step instead of one candidate at a time. In our sweep on two tasks, most batched runs finished in half the wall-clock time or less, and the fastest in about a quarter to a third. Batched settings also achieved higher held-out test scores, from 68.9% to 72.1% on LiveBench-Math (with $2\times2$) and from 49.0% to 60.0% on HoVer (with $8\times1$). A simple runtime model explains the measured speedups under a finite worker pool, and batched settings also spend their metric-call and dollar budgets more efficiently. We frame the axes along which a reflective optimization step can scale in parallel, a first look at allocating a fixed optimization budget across them.

Abstract: Besides domination-based algorithms such as the NSGA-II, the decomposition-based MOEA/D is one of the most popular and successful optimization heuristic for multi-objective problems. It decomposes the problem into several, single-objective subproblems. By solving these in a co-evolutionary manner and archiving non-dominated solutions, it approximates or computes the Pareto front. Interestingly, so far no advantages from the co-evolutionary communication of subproblem solutions were proven. In this first mathematical runtime analysis addressing this research question, we investigate how the mutation-based MOEA/D solves the classic OneMinMax benchmark when solutions are exchanged among subproblems in different distances $r$. We show, for the first time, that the MOEA/D profits from an exchange of solutions, and that already very small communication radii $r$ are sufficient to maximally profit from this feature of the algorithm. More precisely, for all values $r \ge 1$ we prove a tight runtime guarantee of $Θ(n^2)$, where $n$ is the problem size. This is asymptotically faster than the known bound of $O(n^2 \log n)$ guarantee for $r=0$, which we complement with a matching bound of $Ω(n^2 \log n)$.

Abstract: The neighborhood topology of particle swarm optimization (PSO) determines information propagation and regulates the exploration-exploitation balance. However, existing evaluation methods based on fitness and spatial diversity fail to characterize the dynamic evolution internal interaction structure of the population. To address this gap, this paper focuses on small-scale PSO (20-50 particles) with regular topologies and models the iterative process as a dynamic weighted influence graph. We construct a multi-scale information-flow evaluation framework covering micro-, meso-, and macro-level metrics, and propose a freeze-point detection method based on effective nodes, a macro-level distributional metric. Experiments on CEC2017 benchmark functions show that the information flow structure freezes early in the search (3 percent - 22 percent for representative configurations and 15 percent-49 percent for global medians) and remains unchanged thereafter. The freeze time follows the order Gbest$<$ Von Neumann $<$Ring, consistent across standard, LDIW, and TVAC variants. Parameter-adaptive variants delay freezing without altering the topological order, and the freeze point is significantly negatively correlated with final solution quality. Intervention experiments under Gbest show that weakening hub nodes delays freezing and improves performance on multimodal functions, while strengthening hubs produces no additional deterioration due to saturation, revealing a decoupling between structural and performance effects. Among the three metric levels, only macro-level distributional metrics effectively capture freezing under regular topologies, because micro- and meso-level structures are predetermined by the topology. These findings reveal the freezing mechanism of information flow in small-scale PSO and provide structural prior signals for adaptive topology control.

Abstract: We revisit the variational characterization of Nesterov's celebrated method and study when Nesterov trajectories are eventually locally minimizing in a variational sense, a notion that we define precisely. We show that this question is equivalent to the nonnegativity of an associated Schrödinger operator acting on variations of the trajectory. For power objectives of the form \[ f(x)=\frac{1}{p}\|x\|^p \] (where $\| \cdot\|$ denotes the Euclidean norm), we give a complete characterization: Among all exponents $p\ge 2$, the Nesterov trajectory with non-zero initial position is eventually locally minimizing if and only if $p=4$. We further interpret this exceptional exponent through the stability theory of semilinear elliptic equations.

Abstract: Classical optimization and learning-guided search offer complementary approaches to open problems in discrete mathematics, but their outputs are difficult to compare when search spaces, validity checks, and stopping rules differ. We present \emph{RadonBench}, a benchmark for constructing maximum-cardinality Radon-free lattice configurations: finite subsets of $\mathbb{Z}^d$ admitting no integral Radon partition. It remains open whether Radon-free configurations with $11$ points in $\mathbb{Z}^3$ or $21$ points in $\mathbb{Z}^4$ are possible. If they exist, they would improve the corresponding lower bounds on the integer Radon numbers. We compare exhaustive extension, greedy local search, PatternBoost, AlphaZero, and AlphaEvolve. All outputs use one success predicate: a convex target-sized configuration succeeds exactly when its number of integral Radon points is zero. PatternBoost and AlphaZero construct verified Radon-free configurations of sizes $10$ and $20$, one point short of the respective targets; AlphaZero's best convex target terminals have two integral Radon points in both dimensions. Exhaustive extension over fixed seeds and declared finite search boxes likewise reaches a minimum of two integral Radon points. Greedy local search and AlphaEvolve each produce target-sized near misses with one integral Radon point. No audited run reaches zero. These results establish reproducible near-miss baselines over explicitly scoped search domains.

Abstract: We give an algorithm for computing an $\eps$-approximate fixed point of any $d$-dimensional nonexpansion with respect to the $\ell_\infty$ norm in time $\log(1/\eps)^{O(\log^2 d)}$. At its core, our method uses a boosting argument to geometrically shrink a search region. Among other consequences, our algorithm implies that simple stochastic games can be solved in quasipolynomial time.

Abstract: We prove that, for every \(\varepsilon>0\), a \((\frac{k}{2k-1}+\varepsilon)\)-approximation algorithm for maximizing monotone \(k\)-submodular functions requires exponentially many value-oracle queries. This improves the \(\frac{k+1}{2k}\) oracle bound of Iwata, Tanigawa, and Yoshida. Their \(\frac{k}{2k-1}\)-approximation algorithm therefore gives the exact threshold in the monotone setting. For \(k\ge5\), our bound also matches, up to an arbitrarily small loss, Aggarwal's recent approximation guarantee for non-monotone \(k\)-submodular maximization.

Abstract: Heterogeneous optimization problems require an optimizer to decide not only which candidate solutions to evaluate, but also which objectives and constraints should be evaluated under different costs and fidelities. Existing SOTA methods such as MFE-NSGA-III address this setting with manually designed selection metrics, but their behavior is largely determined by predefined heuristic rules and hyperparameters. In this work, we propose REP-HO, a reasoning-based evaluation planner that replaces the rule-based mixed-fidelity selection step with a large language model (LLM) planner. The LLM receives the same primitive optimization-state information as the heuristic planner, such as surrogate predictions, uncertainties, evaluation costs, and budget status. We first show that REP-HO achieves competitive performance compared to the human-designed metric-based method (MFE-NSGA-III). We then further enhance each decision with an adaptive nearest-neighbor summary, leading to better and more stable results. We evaluate REP-HO across five heterogeneous regimes covering cost imbalance, constraint difficulty, surrogate fidelity, budget level, and objective conflict. Across these regimes, one of the two variants obtains the best mean ranking score in every regime, REP-HO(Deepseek) is first or second throughout, and the planner shows adaptive allocation behavior, such as prioritizing constraint evaluations under tighter budgets and reallocating constraint evaluation budget when surrogate reliability changes. These results suggest that LLMs can serve as flexible evaluation planners for heterogeneous optimization, offering a general alternative to fixed hand-designed metrics.

Abstract: We give faster algorithms for sampling from the uniform distribution over a polytope. Our main result is an algorithm that, when applied to a polytope with $m$ constraints in $n$ dimensions, uses $m \cdot \mathrm{poly}\log(\frac m δ)$ iterations, each running in time $O(mn^{ω- 1})$, to return a sample within total variation distance $δ$ to the uniform distribution. When $m =O(n)$ and $δ= \mathrm{poly}(\frac 1 n)$, our runtime improves upon all prior works by an $Ω(n^{2/3})$ factor. Our sampler is based on the recent \emph{functional stochastic localization} framework of [GST26], and achieves an iteration count nearly matching the Poincaré constant of the uniform distribution with respect to the logarithmic barrier.

Abstract: We prove in the affirmative a conjecture of Brlek and Reutenauer (2003), which relates unambiguous rational sets to a well known result due to Khovanskii and Pukhlikov (1992) and independently to Lawrence (1990) on the valuations of the algebra spanned by the indicator functions of polyhedra.

Abstract: Acyclic Unique Sink Orientations (AUSOs) of the binary hypercube form an elegant abstraction of various optimization problems, ranging from linear programs whose feasibility region is a combinatorial hypercube, to 2-player turn-based stochastic games (2TBSGs) in which one of the players has two actions from each one of its states. The currently fastest known randomized algorithm for finding the sink of such AUSOs is the Random-Facet algorithm of Kalai (STOC 1992), and of Matoušek, Sharir and Welzl (SoCG 1992), adapted to the binary hypercube by Gärtner (RSA 2002). The expected number of vertex queries made by the algorithm is at most $e^{2\sqrt{n}}$, where $n$ is the dimension of the cube. We design an improved randomized sink-finding algorithm that performs an expected number of only $e^{\sqrt{(2+o(1))n}}$ queries. This is the first improved bound for the problem in over two decades. The improved algorithm is obtained by adapting improved versions of the Random-Facet algorithm for solving general linear programs, and LP-type problems, developed by Kalai (STOC 1992), Gärtner (FOCS 1992), and by Hansen and Zwick (STOC 2015). The improved bound is not implied by these much more general algorithms. A key feature of the improved algorithm and its analysis is its simplicity. The structure of binary AUSOs enables a simpler description and a much simpler analysis of the adapted algorithm. In particular, the recurrence relations bounding the complexity of the algorithm depend on only one variable, and not two, and they can be solved using fairly straightforward techniques.

Abstract: Multi-objective optimization involves diverse preferences over conflicting objectives. Rather than returning a finite set of trade-off solutions, preference-conditioned Pareto set learning seeks to map a specified preference directly to a corresponding Pareto optimal solution. Yet in many real-world scenarios, optimization must proceed offline without further objective evaluations, leaving only a limited static dataset that sparsely covers the decision and objective spaces. Existing methods either optimize separately for each preference or train a problem-specific Pareto solution model, requiring repeated search across preferences or repeated training across problems. We introduce in-context preference-conditioned Pareto set learning (ICPSL), a framework that amortizes both preference-wise solution search and problem-specific model construction. ICPSL pretrains a foundation-style hypernetwork on a large and diverse set of synthetic prior data. Given the offline observations of an unseen problem, the pretrained hypernetwork learns to generate a problem-specific, preference-conditioned Pareto solution model in a single forward pass. The generated model can then be reused to predict solutions for arbitrary preferences without problem-specific training or preference-wise search. Experiments on synthetic and real-world problems show that ICPSL produces high-quality Pareto solutions while providing substantial speedups in problem-specific solution-model construction.

Abstract: Extremal constructions play a central role in the design and analysis of algorithms. Integrality-gap instances expose the limitations of relaxations, while hard configurations reveal structural sources of computational difficulty. Such constructions, however, are still typically discovered through a slow, manual process. This paper investigates how machine learning can accelerate their discovery to support algorithm design and analysis. We build on PatternBoost, a framework alternating local search with generation from a trained transformer, using its Axplorer implementation to investigate three problems: integrality gaps for maximum independent set of rectangles, integrality gaps for stabbing unit squares with axis-parallel lines, and hard instances for recursive guillotine cutting. By exploring candidate representations, local search operations, and surrogate models used to rank candidates before exact verification, we carefully adapt the framework to these combinatorial geometry problems.While naive representations and search strategies fail to find non-trivial instances, our tailored adaptation discovers numerous non-trivial constructions. For stabbing unit squares, the method reaches a gap of $20/13\approx 1.54$, improving on the literature's best known lower bound of $3/2$. Purely automated search converges rapidly and then plateaus, prompting the question of what is required to push the bound higher still. We find that initial instance diversity plays an important role: with modest expert input via a structured staircase construction, the system reaches an improved integrality gap of $14/9\approx 1.555$. These results demonstrate that machine learning and domain expertise can be combined to discover new extremal constructions that advance algorithmic lower bounds.

Abstract: Expected improvement (EI) based on stochastic kriging (SK) surrogate is widely used in continuous simulation optimization, but unknown heteroscedastic noise on simulation output makes three basic design choices consequential: the incumbent used to define improvement, the replication level at each sampled point, and the terminal recommendation. Existing methods typically either assume the noise levels are known or choose these components separately through heuristic rules or empirical calibration. We instead derive a high-probability upper bound on the simple regret and systematically choose the three components within general rule classes to tighten this bound. We first establish a finite-sample confidence bound for the SK predictor with unknown and estimated noise levels as a foundation for the three design choices. We then analyze general posterior-quantile families for both the incumbent and the terminal recommendation. For the incumbent, minimizing the regret bound selects a lower quantile coefficient that matches the coefficient in the SK predictor confidence bound. For the terminal recommendation, minimizing the same bound instead selects the sampled point with the best upper posterior quantile under the same coefficient. These opposite quantile choices reflect the different roles of search and implementation in simulation optimization. Intermediate simulations are exploratory with an optimistic incumbent, while the terminal recommendation is implemented in the real system and is therefore relatively conservative about the remaining posterior uncertainty. Under mild conditions satisfied by common kernels, minimizing a tractable kernel-specific bound yields $Θ(\log T)$ replications. To our knowledge, this is the first simple-regret guarantee for a continuous-domain EI-type simulation optimization method that learns unknown heteroscedastic variances through replications.

Abstract: Test-time scaling has enabled LLMs to solve increasingly complex reasoning, discovery, and optimization tasks. However, existing approaches often gather observations while seeking evidence for current hypotheses, then generate new hypotheses conditioned on those hypotheses and their associated evidence. This couples observation discovery to existing hypotheses, favoring local refinements over exploration of unknown hypotheses regions. Inspired by the hypothetico-deductive method (Hempel, 1945) and dual-space models of scientific discovery (Klahr & Dunbar, 1988), we introduce SHELDON, a scientific reasoning harness that separately represents and jointly explores hypothesis and observation spaces. By actively seeking new observations along with new hypotheses, SHELDON can discover observations that its current hypotheses do not explain. These unexplained observations expose gaps in its understanding, motivate new investigations, and guide the generation of hypotheses beyond refinements of existing ones. Across diagnosis, abstract reasoning, and optimization benchmarks, SHELDON outperforms evaluated baselines by up to 26.7 points. Additional analysis shows that SHELDON generates up to 4.5x more novel solutions.

Abstract: Convex relaxations make nonconvex optimization tractable by replacing hard nonconvex problems with convex surrogates that give valid bounds. Yet their design is still largely manual: a human expert must decide where to relax, which relaxation to use, and how much formulation growth is worth paying. We show that this design process can be learned. The method proposed in this paper, AutoConvexRelax, represents the symbolic formulation of a nonconvex problem as a tripartite graph linking variables, symbolic terms, and their objective or constraint positions. It then learns a policy that composes symbolic rewrites into problem-specific relaxation trajectories as the formulation evolves. We prove that finite action sequences generated by the learned policy preserve relaxation validity under explicit domain, occurrence-context, and representation-consistency conditions. We further establish computable error bounds for mixed relaxations and sufficient conditions for feasible recovery and global-optimality certification. As a first instantiation, we apply the method to nonconvex quadratic problems and fractional problems. On 241 benchmark quadratic instances, the learned policy strengthens McCormick lower bounds by 13.52%. Its mean lower bound is 3.70% below strengthened semidefinite relaxations on the normalized comparison metric, while it uses substantially smaller formulations. The learned policy retains this trade-off on fractional instances and reaches near-semidefinite bounds on beamforming and sensor-network localization instances. These results suggest that convex relaxation design can move from a fixed handcrafted choice to an automated learning process that adapts to problem structures. Code is available at https://github.com/auto4opt/AutoConvexRelax.

Abstract: Mixed-integer programming (MIP) solvers are widely used, yet their performance can vary substantially even under seemingly neutral perturbations. Such variability complicates benchmarking, and the interpretation of computational studies. We introduce the Weighted Variability Score (WVS), a robust multi-component metric that aggregates dispersion in normalized runtime, final optimality gap, LP iterations, and branch-and-bound nodes across repeated runs. On a large subset of MIPLIB 2017 instances, each solved in 100 randomized runs under a four-hour limit, WVS reveals a right-skewed distribution with a substantial upper tail: many instances are essentially stable, but some exhibit large variability that is not explained by basic properties alone. Notably, this experiment proved optimality for the open instance 'neos-5045105-creuse' in two SCIP runs. Using static descriptors and root-level features, a tree-based predictor with a gated tail uplift achieves a mean absolute error of 0.058 and Spearman rank correlation of 0.65 on a sealed test split. Feature-importance analysis highlights root degeneracy, matrix structure, coefficient ranges, and presolve reductions as the main drivers of high WVS. Because variability is systematic and partially predictable from inexpensive early-run features, our findings offer a practical tool for preemptive tuning and lay a foundation for fairer evaluation of new algorithms.

Abstract: Presolve strongly affects mixed-integer programming (MIP) performance, yet learning-based methods only optimize parameter configurations and cannot express the non-commutative temporal dependencies among actions, whose default order is nearly unique on most domains, yet functionally necessary: artificially shuffling the order of the same sequence inflates the tail of the solve-time distribution by up to several-fold. We recast presolve planning as autoregressive sequence generation over a unified atomic action space, moving the decision object to action sequences; we call this framework ORDO---Operation-level Round-aware Dynamic Ordering for MIP Presolve. Its payoff is cross-domain generalization: on multiple unseen domains it attains end-to-end zero-shot speedup---to our knowledge the first for presolve action sequences---varying by domain and not explained by corpus richness, the strongest domain reaching the largest speedup once racing is added. Deployment uses sequence racing, in which candidate sequences run concurrently and the winner is kept, enabled by an execution-and-observation facility, added by modifying the SCIP source, that injects sequences along the native path and records which actions actually execute and in which round.

Abstract: The existential theory of the reals asks whether polynomial constraints with integer coefficients have a real solution. We give a proof placing this problem in the counting hierarchy. The first argument is intended to expose the essential steps, and a separate analysis lowers the bound to $\exists \mathbb{R}\subseteq\textsf{BPP}^{\textsf C_3\textsf P}\subseteq\textsf C_4\textsf P$, the fourth level of the hierarchy. For each fixed $w$, sentences with $w$ alternating real quantifier blocks lie in $\textsf C_{9w+17}\textsf P$. We then treat exact semidefinite feasibility, PosSLP, square-root sum, geometric real counting, Euler characteristic, and complex feasibility in separate applications. The corresponding bounds include $\textsf{BPP}^{\textsf C_2\textsf P}$ for general SDP, $\textsf{BPP}^{\textsf{PP}}\cap\textsf{P}^{\textsf{NP}^{\textsf{PP}}}$ for PosSLP and square-root sum, and $\textsf{FP}^{\textsf C_4\textsf P}$ for total geometric real counting. Note: These proofs were discovered by ChatGPT after a series of conversations ending on September 29th 2026. A group of researchers has been working to digest the proof, and while the most essential arguments appear correct, we are endeavoring to give this result the treatment it deserves and a proper exposition and development to benefit of the community. However, on October 6th, OpenAI released a very similar result, with a slightly weaker bound. While we work to improve our exposition of this proof, the current version has been uploaded as a service to the community to compare the different proof techniques. While the listed author takes responsibility that the proofs appear to be correct, he has not played a nontrivial role in developing them, and believes the human value will be in good exposition and canonicalization of the results.

Abstract: We determine the hundredths digit of the real Grothendieck constant by proving $1.773 \leq K_G \leq 1.7799$. Furthermore, numerical heuristics suggest the estimate $K_G \approx 1.779$. We start by simplifying the Gaussian duality theory of $K_G$, based on Hermite projection games and Krivine rounding schemes. The bounds are obtained by a heuristic primal/dual search for these objects, followed by rigorous {\em certification} of their values. The search and certification are performed using AI tools, and in particular the rigorous certificates are very large (analytical reductions to thousands of numerical inequalities).

Abstract: Maximizing a submodular function subject to a matroid constraint is a cornerstone problem in combinatorial optimization. The state-of-the-art algorithm for this problem obtains $0.401$-approximation~\cite{buchbinder2024constrained}, and the state-of-the-art hardness result shows that no polynomial time algorithm can obtain better than $0.478$-approximation for this problem~\cite{oveisgharan2011submodular}. In this work, we present the first improvement in $15$ years for the hardness result, showing that no polynomial time algorithm in the value-oracle model can obtain better than $8/17 \approx 0.471$-approximation, even for the special case of a cardinality or (simplified) partition matroid constraint.

Abstract: Landscape analysis has been successfully applied to understand complex optimisation problems, gain insights into algorithm behaviour, and automate algorithm selection and configuration. Although many landscape analysis techniques have been developed over the last decades, it remains difficult for researchers and practitioners to decide which approaches are appropriate and to implement them in practice. Some tools are available, but these are either restricted to particular problem domains (e.g., unconstrained black-box continuous optimisation), or are limited in what they model and measure. In addition, output from landscape analysis is often not easily interpretable, especially when computed landscape features do not correspond with aspects of problems that practitioners are familiar with. In this paper, we introduce a principled approach for explainable landscape analysis (XLA) with an associated Python package called pyXla. The approach is generic in that it applies to problems with different representations (continuous or combinatorial), with single or multiple objectives, with or without constraints. The extent of analysis provided by the XLA framework depends on the data available, with richer analysis offered as additional information is provided by the user. We demonstrate the explainable output produced by pyXla on a selection of hand-crafted problems with diverse landscape characteristics.

Abstract: We present an oracle-optimal method for finite-sum monotone inclusions under mean-square Lipschitz continuity. Our switching regularization method finds a point $y$ and a certificate $g\in G(y)$ with $(\mathbb{E}\|F(y)+g\|^2)^{1/2}\le\varepsilon$ using $\mathcal{O}(n+\sqrt{n}LR/\varepsilon)$ expected component evaluations and resolvent evaluations. It removes the additive $n\log n$ cost of restarting a variance-reduced solver at every regularization stage by switching to a centered stochastic proximal iteration at regularization strength $L/\sqrt{n}$. Carrying an operator estimate between the remaining stages limits their total cost to $\mathcal{O}(n)$. A matching $Ω(n+\sqrt{n}LR/\varepsilon)$ lower bound holds for randomized linear-span component-oracle algorithms with adaptive stopping and expected query budgets. Thus, for $0<\varepsilon\le LR/2$, our method attains the optimal worst-case expected component complexity in this oracle model, up to universal constants.

Abstract: Ising machines have emerged as promising hardware accelerators for intractable optimization and sampling problems, yet their practical impact increasingly hinges on the co-design of algorithms and hardware, where algorithmic demands shape new architectures and new hardware capabilities inspire entirely new algorithms. In this Review, we survey probabilistic algorithms designed for portability across diverse Ising platforms, advocating a top-down perspective that prioritizes principled methods with provable guarantees. We cover foundational methods such as simulated annealing and parallel tempering, including two-dimensional extensions that natively encode hard constraints, and examine approaches that expand the scale of solvable problems from cluster mean-field methods to variational samplers. We highlight the Probabilistic Approximate Optimization Algorithm (PAOA), a classical analog of QAOA that emerged directly from probabilistic hardware development, and explore how generative AI and Ising machines might reinforce each other: learned models propose global moves to accelerate optimization, while probabilistic techniques improve inference in large language models. Much as quantum computing has seen algorithms co-evolve with hardware, probabilistic and Ising computing stand at a similar inflection point. We outline a co-design framework for accelerating the capabilities and adoption of next-generation Ising machines.

Abstract: LLM-guided evolutionary methods, such as AlphaEvolve, have emerged as powerful approaches for challenging computational optimization problems, such as circle packing. However, prior work typically optimizes performance gain over a fixed number of iterations. We argue that practical optimization should maximize gain per unit cost. To this end, we propose FrugalEvo, a cost-aware evolutionary framework where a stronger, higher-cost LLM explores solution strategies, and a cheaper LLM implements them and iteratively refines the resulting code. We also design a cache-efficient evolution process, where our harness and prompts maximize the sharing of prefixes across different evolution steps, to improve cache reuse. To measure solution quality throughout a fixed cost budget, we introduce Budget-Aware Area Under the Curve (BA-AUC), defined as the area under the best-so-far evaluation score curve over cumulative LLM cost, up to the budget. Across 10 mathematical and systems optimization tasks, FrugalEvo matches or surpasses state-of-the-art baselines, including OpenEvolve, ShinkaEvolve, AdaEvolve, and EvoX, in final solution quality and achieves higher BA-AUC on 9 tasks. It also achieves higher average performance than these baselines on 10 algorithmic optimization tasks from ALE-Bench-Lite. Notably, on circle packing, FrugalEvo achieves new state-of-the-art performance with GPT-5.6 Terra and Luna for only 1.68 USD and with GLM-5.3 and its Flash variant for only 0.55 USD, matching or surpassing all baselines, including multi-agent methods such as CORAL and SwarmResearch, which cost approximately 50 USD on average.

Abstract: This paper studies the complexity of convex optimization using lazy second-order oracles (Doikov, Chayti, and Jaggi, ICML 2023), where an algorithm queries gradients every iteration and Hessians once per $m$ iterations. Under this setting, we show a lower bound of $Ω(m+ m^{1/7} ε^{-2/7})$ on the number of total iterations to find an $ε$-solution using a novel block zero-chain construction. Then we propose a novel method that achieves a new upper bound of $\tilde{\mathcal{O}}(m+ m^{1/7} ε^{-2/7})$, which significantly improves the prior one (Chen, Liu, Luo, and Zhang, COLT 2026) of $\tilde{\mathcal{O}}(m+ m^{13/21} ε^{-2/7})$ and is tight up to logarithmic factors.

Abstract: Can programmable radio propagation serve as computational depth rather than only as a communication channel or one-shot analog transform? We introduce the Radio Alternating Operator Ansatz (RAOA), a recurrent computing architecture that alternates an energy-derived problem update with a mixing update over a persistent latent state. Recomputing the problem field after each mix makes repeated passes compositional even when the same learned controls are reused across depth. We evaluate this idea through exact discrete optimization, constrained programmable-propagation simulation, and pretrained-model adaptation. On discrete objectives, repeated execution can improve solution quality without increasing the learned-control count, and the same formulation handles higher-order interactions directly. A passive phase-only free-space model further shows that the required operators can be approximated by programmable propagation while retaining useful downstream behavior despite realization error. When inserted as a zero-initialized residual adapter, RAOA adapts pretrained language models with WikiText performance close to a matched shallow MLP across three model families, while reasoning-task transfer remains model-dependent. Together, these results connect alternating-operator computation, programmable radio propagation, and neural adaptation within one recurrent framework. The RF realization evidence is simulation-based rather than a hardware demonstration.

Abstract: Decision-focused learning for linear optimization is complicated by the discontinuity of the optimizer, where small cost errors may leave the decision unchanged or move it to a different vertex. We show that this non-smooth pointwise behavior becomes locally quadratic after averaging over the data distribution, and we derive the curvature in closed form, specifically, a matrix-valued measure supported on the walls of the normal fan. This measure depends only on the feasible set, with the data distribution entering only as a weight. We then offer a tractable approximation for this curvature, computable with just one projection to the feasible set. We prove that the approximation weakly converges to the true population curvature. We offer one application of our findings, a decision-aware scenario generation method for expected-cost linear optimization. Our experiments test the quadratic and weak convergence laws and show a 30.8% regret improvement over uniform allocation on battery arbitrage.

Abstract: Mixed-integer linear programs (MILP) model many real-world decision problems, motivating machine-learning methods that exploit recurring structure to accelerate MILP solving. MILPs can admit many equivalent formulations: integrality-preserving changes of variables and the addition of redundant constraints can alter their formulations while preserving the optimization problem. We leverage these reformulations as a source of self-supervision for learning general-purpose representations of MILP variables and constraints. We characterize the affine reformulations that are valid for every input instance, and distinguish re-descriptions, which leave variables unchanged, from substitutions, which transform them predictably. Building on equivariant self-supervised learning, we introduce ReMILP (reformulation-contrastive MILP representation learning), which jointly trains a graph neural network and a hypernetwork to predict how variable embeddings transform under changes of variables. Without solver-derived labels, ReMILP learns representations that exhibit the intended invariance and equivariance on unseen problem classes. Across binary solution, constraint activity and integrality gap prediction, these representations carry task-relevant information when frozen and provide a useful initialization for fine-tuning.

Abstract: We give the first sublinear-time classical solvers for sparse semidefinite programs in the bounded-radius regime, without low-rank assumptions or Frobenius norm dependence on the constraint matrices. For constant precision and bounded primal and dual radii, prior quantum algorithms of Brandão et al. (2019) and van Apeldoorn and Gilyén (2019) achieved $\widetilde{O}(\sqrt{n}+\sqrt{m})$ dependence on matrix dimension $n$ and constraint number $m$. Compared with the $\widetilde{O}(mn)$ runtime of existing classical methods, this suggests a quartic quantum speedup when $m \approx n$. Beyond a usual Grover speedup, this separation relies on the Quantum OR lemma, whose sample-reuse mechanism decouples the cost of Gibbs-state preparation from constraint search. We show that this reuse mechanism is classically realizable for sparse SDPs. Our main technical contribution is a classical procedure for simultaneously estimating many expectation values with respect to a sparse Hamiltonian's Gibbs state. This combines randomized Lánczos filtering with an efficient sampling-based estimator. We also introduce a stochastic online-learning framework for SDP solving, substantially improving accuracy-dependence over standard oracle-based MMWU approaches. Let $s$ denote the the input matrix sparsity and $γ:=Rr/\varepsilon$ capture dependence on the primal $(R)$ and dual $(r)$ radii as well as target accuracy $(\varepsilon)$. When $γ^2\leq\min\{m,n/s\}$, our solver runs in time $\widetilde{O}\left(nsγ^{4.5}+msγ^2\right)$. For $γ=O(1)$, this is $\widetilde{O}\left((n+m)s\right)$ and sublinear in the $O(mns)$ input size. Similar to the quantum algorithms, this matches known lower bounds with respect to $m$ and $n$, up to logarithmic factors. This implies that, with respect to dimensions $m$ and $n$, there is no super-quadratic quantum advantage for generic sparse SDP solving.

Abstract: We introduce a framework for combinatorial optimization using sum-constrained continuous quadratic programs solvable by QCi's Dirac-3S photonic entropy computer. This is enabled by the Motzkin-Straus theorem which provides a powerful bridge between discrete clique problems and optimization over the probability simplex. We demonstrate this framework's versatility by solving constraint satisfaction problems, providing extensive benchmarks on the DIMACS suite. The Dirac-3S platform matches or outright leads two independently implemented classical baselines on more than four-fifths of the benchmark instances, reaching the best known solution on nearly all structured graph families, even outperforming both classical solvers on several of the largest instances tested. On the other hand, well-tuned classical continuous optimizers retain an edge only on the hardest planted-clique instances. This work establishes a viable pathway for solving combinatorial optimization problems using natively analog unconventional computing platforms, while positioning entropy computing as a competitive approach for navigating non-convex landscapes and providing rigorous baselines for an emerging computational paradigm.

Abstract: Population-based optimization methods often use previous search information through successful solutions, parameter adaptation, or operator performance, but they rarely retain the context in which a search behavior succeeded or failed. We introduce Cognitive Memory-Driven Optimization (CMDO), a derivative-free population-based optimizer that represents experience as the relationship between search context, search behavior, and observed outcome. CMDO organizes these experiences across working, episodic, and consolidated memory, retrieves them according to similarity with the current search state, and uses both positive and negative evidence to guide subsequent search. Retrieved experience does not replay previous candidate locations; instead, it selects search recipes that are reconstructed from the current population through exploratory, directed, and local search behaviors with adaptive search geometry. We evaluate CMDO on selected Blackbox Optimization Benchmarking test suite on COCO (BBOB/COCO) and Congress on Evolutionary Computation 2017 (CEC2017) problems against DE, CMA-ES, SHADE, GWO, HHO, and ORCA, and further study its application to seven-parameter photovoltaic model estimation using measured current--voltage data. The results show problem-dependent but competitive optimization performance, including the lowest median error among the compared methods on CEC2017 F10. More importantly, analysis of the search traces shows that context-dependent recall changes the distribution of executed search behaviors, while unsuccessful experiences remain available as negative evidence for later decisions, showing that accumulated experience directly influences subsequent search behavior. These results support the use of explicit context--behavior--outcome memory as an active mechanism for controlling population-based search.

Abstract: Recovering probability measures from their moments has numerous applications, esp. in connection with the method of moments in statistics and optimization. In the setting where measure need not be finitely atomic, but its support is known to be compact and semialgebraic with codimension at least one, the problem is still open. We combine moment-matrix kernel information with the Christoffel--Darboux kernel to provide a discrete approximation of the support. To validate the proposed approach, we test our algorithm on analytically computed moments and pseudo-moments arising from polynomial optimization problems without unique global minimizers. This complements well-known recent work on recovery of measures with algebraic support, where the kernel of a moment matrix can reveal polynomials vanishing on the support, and on recovery of sufficiently regular full-dimensional supports, where estimators constructed by thresholding the Christoffel--Darboux kernel are known to converge asymptotically to the support.

Abstract: One of the most important problems in computational algebra is the computation of real solutions of a polynomial system. Lasserre, Laurent, and Rostalski introduced a numerical method to compute the real radical of an ideal defining a finite set of real points. Their method is based on moment matrices and reduces the problem to solving a semidefinite program. In this work, we build on this method and on the work of Laurent and Mourrain to exploit the sparsity of the input polynomials. For this purpose, we extend the theory of moment matrices to the setting of affine toric varieties. We prove that, for sparse polynomial systems, this approach allows us to recover real points on affine toric varieties. Moreover, we show that a projection of the spectrahedron associated with the semidefinite program is, in fact, a polytope, and we propose a variation of the method to compute only one real solution of the system.

Abstract: We show that computing a Nash equilibrium is CLS-complete for linear network congestion and network coordination games. As a result, finding a KKT point of a bilinear polynomial is CLS-complete.

Abstract: We give an exposition of the recent polynomial time algorithm of Guo, Fang, and Lu for the Komlos problem. We simplify various arguments, and highlight the key new spectral potential idea and how the algorithm follows naturally from it.

Abstract: We prove nearly quadratic lower bounds for randomized algorithms for linear optimization and uniform sampling over convex bodies in the membership oracle model. For linear optimization, this matches the known nearly quadratic upper bound up to a polylog factor in the dimension. For uniform sampling, this improves on the previous linear lower bound. Our construction also implies the same lower bound for volume estimation.

Abstract: The Komlós conjecture, now a theorem, asserts that whenever the columns of a matrix $A\in\mathbb{R}^{m\times n}$ have Euclidean norm at most one, some signs $\varepsilon\in\{-1,1\}^n$ make every coordinate of $A\varepsilon$ bounded by an absolute constant. Guo, Fang, and Lu gave the first polynomial-time algorithm for finding such signs, a deterministic spectral signing procedure with discrepancy $8272$ and running time $O((mn^9+n^{10})\log(m+n))$. We give a deterministic algorithm that finds signs with $\|A\varepsilon\|_\infty<99$ using $O(mn+n^{ω+2}\log^3 n)$ arithmetic operations, where $ω>2$ is any fixed attainable matrix-multiplication exponent; with the current bounds on $ω$ this is $\widetilde O(mn+n^{4.372})$. Our algorithm uses the same framework: it rounds a single fractional coloring and watches all rows through the top eigenvalue of a Gram matrix of energy-corrected barriers. Steps follow flat directions, rescaled so that no barrier near its threshold moves faster than a constant, and the regularizer grows as coordinates freeze; together these bound the number of updates by $O(n^2\log n)$. A weak quadratic charge on the tracked row sums leaves $O(n\log^2 n)$ rows to evaluate at any time, and a motion clock bounds when any other row could approach its barrier. Each update is a short sequence of matrix products. Its direction is read off by conditional expectations from a polynomial soft projector, its small constraint residual is repaired in affine row representations, and one identity accounts for every change of representation. For rational input the algorithm has polynomial bit complexity.

Abstract: Large language models (LLMs) can improve solutions to verifiable scientific and algorithmic problems by spending additional computation at test time. Recent systems achieve strong results with increasingly elaborate evolutionary search harnesses or by updating model parameters during test-time training. We ask how much of this machinery is necessary. We introduce Hill Sampling, a simple procedure that repeatedly samples candidate program edits from a frozen LLM, retains the best program found so far, and conditions all subsequent samples on that program. We evaluate the method on circle packing, sums/differences of sets, and Erdos' minimum-overlap problem using three open-weight models. Hill Sampling sets a new state of the art on circle packing among published methods, improves over the AlphaEvolve reference on Erdos' minimum-overlap problem, and achieves strong results on sums and differences of finite sets. The circle-packing and Erdos results require only hours of wall-clock time on eight NVIDIA H100 GPUs. To our knowledge, we also conduct, the largest study, by parameter count, of evolution strategies (ES) applied directly to LLM weights at test time. Surprisingly, learning the weights is worse than setting the ES learning rate to zero: at zero learning rate, the method is still searching in weight space through fixed random perturbations. Those perturbations can help exploration, but randomness from token sampling is stronger still, and repeated sampling remains substantially weaker than Hill Sampling. These results suggest a simple test-time compute allocation strategy: repeatedly sample edits to the best verified solution found so far, before introducing additional complexity such as adding archives, diversity mechanisms, evolutionary scaffolds, or test-time parameter learning.

Abstract: We show that $1.4<K_G^{\mathbb{C}}<1.404898554746$, where $K_G^{\mathbb{C}}$ is the complex Grothendieck constant. The upper bound improves on Haagerup's bound of $1.40490913\ldots$ from 1987, and the lower bound improves on Davie's bound of $1.33807$ from 1984. The upper bound combines ideas from the real and noncommutative complex cases of Grothendieck's inequality.

Abstract: Large language models (LLMs) enable automated algorithm design (AAD) through reasoning and code synthesis. However, most existing LLM-based AAD methods separate algorithm design from target optimization, deploying a fixed design even as the optimization state evolves. Conventional adaptive optimizers can respond to such changes, but their adjustments remain confined to predefined parameters, operators, or strategies. To address these limitations, we introduce online LLM-based AAD, a novel optimization paradigm that treats the algorithm itself as a state-dependent decision variable. At each stage, LLM agents synthesize an algorithm with new behavior logic from the current optimization state. Executing the generated algorithm advances the search and provides feedback for subsequent designs, coupling algorithm design with target optimization without requiring a separate offline algorithm pretraining stage. To implement this paradigm, we propose OnDesign, a multi-agent framework that reconciles competing design perspectives to synthesize executable algorithms and uses execution feedback to refine how runtime evidence is interpreted for subsequent designs. We evaluate OnDesign across two mainstream black-box optimization paradigms on three scenarios: Bayesian optimization, evolutionary continuous optimization, and evolutionary mixed-variable optimization. Extensive experiments on six benchmark suites and one engineering problem across multiple problem dimensions demonstrate superior overall performance over conventional optimizers and offline LLM-based AAD methods.

Abstract: We study the complexity of finding $(δ,ε)$-Goldstein stationary points of nonsmooth nonconvex Lipschitz functions. By now, it is known that randomized first-order algorithms can solve this task with a dimension-free oracle complexity [Zhang et al., 2020], whereas deterministic algorithms cannot, as their complexity must scale at least linearly with the dimension $d$ [Jordan et al., 2023, Tian and So, 2024]. This leaves open whether deterministic algorithms can nevertheless solve the problem with oracle complexity polynomial in $d$. We answer this question negatively by proving a lower bound of order $(1/ε)^{Ω(d)}$ for deterministic algorithm, closing the exponential gap between the previously known lower and upper bounds and resolving an open problem posed by Jordan et al. [2023]. We further discuss several extensions and implications of this result to weaker stationarity notions, finding a descent direction and deterministic smoothing. Overall, our results establish an exponential computational advantage in nonsmooth nonconvex optimization offered by randomization.

Abstract: Scaling analog and digital Ising machines to larger problems requires overcoming finite device capacity and the cost of communication between devices. We present cluster mean-field theory (CMFT), a framework that partitions a large interaction graph into clusters sized to fit available hardware. Each cluster performs local updates independently, while interactions across cluster boundaries enter through periodically updated effective biases computed from boundary-spin averages. To reduce the error introduced by fixed cluster boundaries, we introduce dynamic partitioning, which cycles through multiple partitions so that interactions approximated by mean fields at one stage can act through instantaneous spins at another. On three-dimensional spin glasses and planted Pegasus instances, dynamic CMFT exhibits power-law decay of residual energy over sweep budgets. This shows that solution quality continues to improve with computational effort despite the mean-field approximation. An automated graph partitioner combined with a weighted recovery ratio provides a practical heuristic for selecting partition combinations on graphs without natural cut directions. We demonstrate CMFT on four GPUs with approximately four million p-bits, reaching comparable energy densities up to 15 times faster than a single-GPU implementation of the full graph. By coupling locally evolving clusters through programmable effective biases, CMFT provides a route to extreme-scale Ising machines on both analog and digital hardware.

Abstract: We compute the convex envelope of (nonconvex) bivariate piecewise linear-quadratic (PLQ) functions (bivariate quadratic functions defined on a polyhedral subdivision). Our algorithm is composed of the following steps: (1) compute the convex envelope of each quadratic piece obtaining piecewise rational functions (quadratic divided by linear function) defined over a polyhedral subdivision; (2) compute the (Legendre-Fenchel) conjugate of each resulting piece to obtain piecewise quadratic functions defined over a parabolic subdivision; (3) compute the maximum of all those functions to obtain the conjugate of the original PLQ function as a piecewise quadratic function defined on a parabolic subdivision; (4) compute the conjugate of each resulting piece; and finally (5) compute the maximum over all those functions to obtain the convex envelope (biconjugate) as rational functions (quadratic divided by linear function) defined over a polyhedral subdivision. Our contribution includes a practical algorithm running in linear time, and proving that the convex envelope is a piecewise rational function.

Abstract: The Matrix Spencer conjecture asks whether any $n$ real symmetric matrices A_1,...,A_n \in \mathbb{R}^{m \times m} of operator norm at most one admit a signing $x\in\{-1,1\}^n$ such that the operator norm of the signed sum is at most O(\sqrt{n \log(2m/n)}) We give a randomized algorithm establishing this bound with polynomial runtime in the real-arithmetic model. We first prove the $O(\sqrt n)$ bound for $m\le n$, resolving the square case, and then obtain the rectangular bound by changing the regularizer. As in earlier algorithmic discrepancy methods \cite{lovettmeka2012,bansalLaddhaVempala2022,pesentivladu2026}, we run a covariance-controlled random walk from the origin of the hypercube, rounding coordinates near its faces and keeping them fixed. Our potential measures a soft spectral edge of the evolving discrepancy matrix perturbed by an operator-valued free semicircular element. Inspired by the free interpolation approach of \cite{bbvh2023}, we combine Lehner's variational formula for the free edge \cite{lehner1999} with spectral Tsallis regularization \cite{allenZhuLiaoOrecchia2015,pesentivladu2026}. This puts the discrepancy and remaining covariance in a single smooth optimization problem. The potential has a finite-dimensional semidefinite formulation. Stability of its optimizer, governed by equations related to the matrix Dyson equation \cite{erdos2019}, lets us find a large subspace in which to move while controlling discrepancy. The square case uses the Tsallis--$1/2$ regularizer; the rectangular case uses a suitable generalized Tsallis power regularizer. Our companion paper \cite{kathuria2026ks} applies these ideas to give an algorithmic proof of Weaver's discrepancy theorem, whose existence proof by [MSS15] resolved the Kadison--Singer conjecture \cite{mss2015}.Lean formalizations of our main discrepancy theorems have been completed and will be released shortly.

Abstract: We give an approximately uniform sampler for satisfying assignments of constraint satisfaction problems that satisfy $4\mathrm e p(Δ+1)^2\le1$, where $p$ is the largest constraint-violation probability under the uniform product distribution, and $Δ$ is the maximum degree of the dependency graph. The algorithm invokes the recent efficient approximate counting algorithm of Liu, Wang, Yin, Zhang, and Zhou as a subroutine and returns a satisfying assignment sampled within total-variation distance $\varepsilon$ of the uniform distribution in $(n+m/\varepsilon)^{O(kΔ\log D)}$ time, where $n$ and $m$ are the numbers of variables and constraints, $D$ is the common domain size, and $k$ bounds the constraint arity.

Abstract: This work studies the stability and convergence of augmented primal-dual dynamics when constraint values are estimated from samples. Unbiased constraint observations can produce a biased augmented multiplier signal, shifting the equilibria of the mean dynamics. For componentwise inequalities, we give a necessary and sufficient condition for preserving the Karush-Kuhn-Tucker (KKT) equilibria and construct a convex example with a locally exponentially stable equilibrium that violates complementarity. To address this bias, constraint values are estimated recursively before forming the augmented multiplier signal. For smooth convex conic problems, a joint energy analysis establishes boundedness of the primal, dual, and estimation states, vanishing estimation error, and almost sure convergence of the primal-dual iterates to a single KKT point under global regularity and bounded conditional second moments. The result allows nonunique solutions and multipliers while keeping the number of samples per iteration fixed. Numerical studies illustrate the predicted equilibrium bias and examine convergence with nonunique KKT points and nonlinear constraints.

Abstract: We introduce cvxgenrust, an open-source tool for generating custom Rust code that solves families of parameterized convex optimization problems modeled in CVXPY. cvxgenrust canonicalizes a problem family, extracts affine maps to Clarabel cone-program data, and generates a specialized Rust crate that updates parameters and calls Clarabel natively at runtime. The generated solver can also be exposed to Python and registered as a custom CVXPY solver. Our code generator supports a wide range of convex optimization problems up to semidefinite programs and exponential-cone problems. Numerical experiments show reduced runtime relative to direct CVXPY solves and performance comparable to CVXPYgen on shared problem classes.

Abstract: Let $A_1,\ldots,A_n$ be independent $d \times d$ real symmetric Gaussian random matrices, and consider the linear operator $A(x) = n^{-1/2}\sum_{i=1}^n x_i A_i$, $x\in \mathbb{R}^n$. We construct an iterative algorithm in the Approximate Message Passing family which iterates over $A$ and its adjoint $A^*$, and establish a state evolution result which characterizes its behavior in the limit $d\rightarrow \infty, 2n/d^2 \rightarrow α$ in terms of a correlated Gaussian-semicircular process in a free probability space, in the sense of strong convergence of operators. We then apply this iteration to the random matrix discrepancy problem which asks for a binary vector $x \in \{-1,+1\}^n$ such that $A(x)$ has a small operator norm. Our algorithm achieves an operator norm $2σ(α)$, for an explicit expression of the standard deviation $σ(α)<1$ for all $0<α<α_* \simeq 5.74$. This resolves the algorithmic question of Kunisky-Zhang (2023) and Maillard (2025) in this interval.

Abstract: We construct counterexamples to Rockafellar's sum conjecture in which two maximally monotone operators satisfy the interior-domain condition but their sum is not maximally monotone. We give one counterexample on $c_0$ and another on $\ell^1$ with its usual norm. We establish a general construction theorem that computes the entire monotone polar of a class of graphs, gives a necessary and sufficient condition for their maximal monotonicity, and shows how a positive rank-one perturbation yields a nonmaximal sum under this condition. We verify the theorem's hypotheses and its maximality criterion on $c_0$, thereby obtaining a counterexample to the conjecture. Furthermore, we construct a bounded linear surjection from $\ell^1$ onto $c_0$ and use it to obtain the counterexample on $\ell^1$.

Abstract: We study the oracle complexity of computing a point with small fixed-point residual $\|T(x)-x\| \leq ε$, for a general norm $\|\cdot\|$ and a self-map $T$ of a compact convex set. We study this problem in the setting where $T$ is nonexpansive with respect to the same norm $\|\cdot\|$ and accessed via an unbiased stochastic oracle with bounded variance $σ^2$. We provide an algorithm that solves such instances for any norm with a weak Rademacher type $q > 1$, with high probability. The algorithm is based on a recursive anchoring technique. For type-$2$ spaces, such as $\ell_p$-spaces for $p \in [2, \infty]$, our algorithm attains stochastic oracle complexity $\tilde O(σ^2 ε^{-3} + ε^{-1})$. We further prove a near-matching lower bound (i.e., matching up to poly-log factors) for such $\ell_{\infty}$-norm instances in high dimensions. Our lower bound holds against any randomized algorithm that succeeds with constant probability. It further extends to settings with ``sparse'' noise, where variance measured with respect to any $\ell_p$ norm is of the same order, ruling out the possibility of improving oracle complexity as a function of $\varepsilon$ by measuring variance in a non-matching $\ell_p$ norm.

Abstract: We present a quantum-inspired Walsh/PCE solver for MaxCut based on sparse Pauli-correlation encodings. Instead of assigning one qubit or one variable to each graph vertex directly, the method represents relaxed binary variables through expectation values of diagonal Pauli/Walsh observables. These correlators are computed classically from sparse Walsh autocorrelations, producing a compact differentiable relaxation of the MaxCut objective. We evaluate the method on selected Gset instances, G1, G6, G12, and G18, and compare it with random search and tabu search over 10 independent seeds. The proposed model uses $801$ active parameters, corresponding to only $0.306\%$ of the full Walsh space over $18$ qubits. After a final bitflip local search, Walsh/PCE achieves approximation ratios of $0.99033 \pm 0.00226$ on G1, $0.95647 \pm 0.01604$ on G6, $0.96007 \pm 0.00951$ on G12, and $0.92964 \pm 0.02202$ on G18, outperforming both baselines on all tested instances. The method also yields the lowest average runtime in all cases. These results suggest that sparse Walsh/PCE representations provide an efficient quantum-inspired route for MaxCut and may be further extended to hardware-based estimation of Pauli/Walsh correlators.

Abstract: Digraphs with constant cut balance admit nearly linear directed cut sparsifiers. This condition requires the total arc weights in the two directions of every cut to be within a constant factor of each other. We ask whether this condition also permits nearly linear spectral sparsification with respect to the energy of the nonlinear directed Laplacian. For a weighted digraph $G=(V,E,w)$, let \[ Q_G^+(x)=\sum_{(u,v)\in E}w_{uv}(x_u-x_v)_+^2, \qquad (t)_+:=\max\{t,0\}. \] This energy agrees with the outgoing-cut function on binary vectors. A spectral sparsifier is a nonnegatively reweighted subgraph that preserves $Q_G^+(x)$ within a factor of $1\pm\varepsilon$ simultaneously for all $x\in\mathbb R^V$. We show that cut balance alone does not yield nearly linear spectral sparsifiers: for constant error, the worst-case support size for simple unweighted Eulerian digraphs is $\widetildeΘ(n^{3/2})$, although Eulerian digraphs are perfectly cut-balanced and admit nearly linear directed cut sparsifiers. In contrast, we prove that every $n$-vertex tournament has a spectral sparsifier with $\widetilde O(n/\varepsilon^3)$ arcs, without any assumption on its cut balance. This includes the transitive tournament, whose cut balance is unbounded. Thus perfect balance does not guarantee nearly linear spectral sparsification, while unbounded imbalance does not preclude it. Finally, we use convex duality to show that preserving $Q_G^+$ also preserves, for every feasible demand vector, the optimum quadratic cost of a nonnegative flow. Hence the guarantee contains information beyond directed cut values.

Abstract: Fourier-Motzkin elimination computes an inequality description of the projection of a polyhedron onto a subset of its coordinates by eliminating one variable at a time. It is used in several areas of optimisation and computer science, and it is a standard way of obtaining the entropic constraints of a causal structure, where the marginalisation over the latent variables produces such a projection. Its limitation is the growth of the intermediate systems of inequalities, which can be doubly exponential in the number of eliminated variables even though the projection itself grows only as a single exponential. In practice the computational overload of the method therefore depends on two choices: how the redundant inequalities are removed after each step, and the order in which the variables are eliminated. We consider both. We first show, by an explicit example, that Imbert's redundancy test cannot be interleaved with redundancy removal by linear programming. We show that the two methods, however, can be combined soundly if the derivation records used by Imbert's test are re-initialised after every step at which linear programming is used. We then propose a rule for choosing the elimination order of the variables that gives a significant computational advantage, however, at the cost of increased resource usage. We demonstrate this advantage on some random polytopes, where the rule reduces the running time by factors of between 6 and 25 compared with the same elimination under a fixed order. For entropic descriptions of causal structures, with more than 250 inequalities and more than 100 variables to eliminate, our rule keeps the number of inequalities handled at each step one to two orders of magnitude lower than a fixed order.