Papers

We introduce a framework for retraction-based convex optimization on Riemannian manifolds, which includes a notion of retraction-specific convex sets and retraction-based gradient projection algorithms. The standard theory of gradient projection algorithms generalizes easily to this framework. Within this framework, we establish convergence results for retraction-based gradient projection algorithms with various stepsize rules. As an application, we use our framework to study the weighted low-rank approximation. We also provide numerical validation of our convergence results on the image completion task.

We study whether the linear condition-number dependence in the stochastic complexity of SAPD+ is necessary for nonconvex-strongly-concave minimax optimization. For jointly $L$-smooth objectives with dual strong-concavity parameter $\mu$, we prove a lower bound that matches the SAPD+ upper bound under the same Moreau-envelope stationarity criterion and the same primal-dual initialization gap. Specifically, when $\sigma\ge\varepsilon$, the worst-case complexity of zero-respecting algorithms is $\Theta(\kappa LG\sigma^2\varepsilon^{-4})$ in the stated accuracy regime, where $\kappa=L/\mu$, $G$ bounds the initial primal-dual gap, and $\sigma^2$ bounds the variance of a general unbiased first-order oracle. The lower bound is realized on a smooth problem class with a bounded dual box. Our construction routes each link of a nonconvex zero-chain through a dual gradient of magnitude proportional to $\varepsilon/\sqrt{\kappa}$, while an undiscovered primal coordinate prevents stationarity. It also yields the primal-gradient lower bound $\Omega(L\Delta(\sqrt{\kappa}\varepsilon^{-2}+\kappa\sigma^2\varepsilon^{-4}))$ after combination with the known deterministic bound, where $\Delta$ bounds the initial primal function gap.

Quantum-inspired evolutionary optimization (QIEO) is a new class of population-based metaheuristic optimization algorithms which represents design variables as a set of qubits and searches a continuous, multi-dimensional landscape through rotation of the qubit's amplitude pair. Every generation rotates those amplitudes toward a single elite, which corresponds to that generation's best. The per-generation cost scales as $O(N_p N_g)$ for $N_p$ chromosomes and $N_g$ genes (decision variables). Production use of such solvers is rarely confined to a single machine class. Prototypes are run on laboratory servers- before moving to rented cloud workstations for more involved campaigns. The largest problems are reserved for leadership-class accelerators. This paper asks whether a single OpenMP~5 source of QIEO, offloaded with \#pragma omp target, is a viable production path in each of those settings. We report three independent, campaigns of the 0/1 knapsack problem against a same-source multi-core Intel CPU baseline. The study comprises approximately 3,000 runs spanning varying chromosome and gene counts, evaluated using both chromosome-level and gene-level offload strategies on the NVIDIA Tesla V100 SXM2, NVIDIA A100 80GB, and AMD Instinct MI300X GPUs. Deployment-specific nuances such as Volta's constant-memory cliffs, Ampere's L2 persistence and cp.async, CDNA~3's Infinity Cache and XCD occupancy are addressed to ensure high performance of these platforms. Results reveal gene-parallel offload achieved geometric-mean speedups of 90$\times$, 136$\times$, and 155$\times$ over a single CPU core on the V100, A100, and MI300X, respectively, and 12$\times$, 17$\times$, and 16.6$\times$ over 72 host threads. Furthermore DetermineElite, the $O(N_p)$ selection of the generation-best chromosome, is found to be better suited to the host than to the device.

Lee and Sidford [LS19] showed that the linear program $\min\{c^\top x\ :\ A^\top x=b,\ l\leq x\leq u\}$ over $x\in\mathbb{R}^m$, where $A\in\mathbb{R}^{m\times n}$, can be solved to accuracy $\epsilon$ using $O(\sqrt n\log^{13}m\cdot\log(mU/\epsilon))$ solves of linear systems in $A^\top\mathbf DA$ for positive diagonal matrices $\mathbf D$, where $U$ bounds the magnitudes of the input and of the initial point. Theirs is the first such bound governed by $\mathrm{rank}(A)=n$ rather than by the number of constraints $m$. Song [Son19] mentioned that reducing the $\log^{13}m$ factor is an interesting future direction. Given an interior point and a finite certified magnitude bound $U$, we give an algorithm that uses $O(\sqrt n\log^{4}m\cdot\log(mU/\epsilon))$ solves of such linear systems, improving the Lee–Sidford bound by a factor of $\log^{9}m$. We conjecture that $O(\sqrt n\log(mU/\epsilon))$ such linear-system solves suffice.

This work addresses the computation of extreme values of observables along trajectories of dynamical systems over an \ks{infinite time horizon}. The setting is intrinsically global: trajectories may exhibit large transient excursions, approach their largest values only asymptotically as time tends to infinity, or escape to infinity in finite time. The proposed framework treats these \mk{ scenarios} in a unified way, so that the resulting extreme values reflect both the transient and the asymptotic behavior of the dynamics. The main difficulty is that neither the time interval nor the state space is assumed to be compact. We address this by combining a radial compactification of the state space with a compactified clock variable and discounted occupation measures. This yields a description of trajectory ensembles on compact sets via occupation measures and the associated Liouville's equation. Extreme values are then characterized through a family of reachability problems for level sets of the observable, which can be solved by a bisection procedure. On the theoretical side, we prove a superposition principle for an infinite-horizon Liouville's equation with a time-distributed terminal measure. This result shows that the measure formulation represents ensembles of stopped admissible trajectories, including trajectories \ks{that never stop}, and establishes equivalence between the trajectory-level and measure-level descriptions.

In many safety-critical applications, control of uncertain dynamical systems relies on observers that estimate states and external disturbances. Neural network observers can improve estimation accuracy, but certifying their Lyapunov stability via Linear Matrix Inequality (LMI) constraints leads to large-scale semidefinite programs (SDPs) that are difficult to solve for large networks. To overcome this scalability bottleneck, we propose a novel two-stage training framework for provably stable neural network observers. Our approach decouples the optimization into a point-guided Lyapunov pre-training phase, which rapidly achieves high estimation accuracy and local stability over sampled states, followed by an LMI fine-tuning phase that efficiently satisfies a strict global Lyapunov stability certificate. We provide formal theoretical guarantees for local stability radii and probabilistic coverage over a prescribed compact error-state domain under specified regularity and sampling assumptions. Experiments on nonlinear control benchmarks and X-29 aircraft ablations show that our LMI-certified neural network observers train significantly faster than direct LMI-based methods and generalize robustly across diverse systems, achieving improved tracking accuracy over a range of observer baselines. The code is available at https://github.com/Berry-Myon/LearningNeuralNetworkObserver.

Modern service platforms often provide customers with real-time congestion information, such as anticipated waiting times, before they decide whether to use a service. This creates an intertemporal tradeoff in service recommendation: directing a customer to a service may generate immediate value, but the resulting congestion can make that service less attractive to future customers. We study this tradeoff through a finite-horizon stochastic optimization problem in which a platform dynamically recommends among multiple services, each represented by a queue, with customers' joining probabilities decreasing with congestion. We first formulate the problem as a Markov decision process and show that optimal decisions generally depend on the joint congestion state, giving rise to a prohibitively large state space. We develop two complementary approximation algorithms that overcome this curse of dimensionality. First, we provide a quasi-polynomial-time approximation scheme that uses truncation and rounding of the model primitives to compress the joint congestion state, while carefully accounting for how these approximations affect the stochastic evolution of the system to establish near-optimality. Second, we develop a polynomial-time LP-guided algorithm that replaces the joint-state representation with service-level marginal information and achieves a (1-1/e) approximation guarantee under substantially more general service and joining dynamics. Finally, we show that monotone joining behavior, whereby customers become less likely to join as a service becomes more congested, marks a fundamental tractability boundary: without monotonicity, the problem is NP-hard to approximate within any constant factor. This establishes monotonicity as a key behavioral property that enables tractable dynamic service recommendation despite the complexity introduced by congestion-dependent joining.

We initiate the study of prophet trading, an online trading model in which a trader interacts with sellers and buyers who arrive in a uniformly random order and whose prices are drawn independently from a common known distribution. The trader aims to maximize the expected gain from trade (GFT), with performance measured against an omniscient prophet that knows the realized arrival order and all prices in advance. Unlike the classical prophet inequality problem, the trader must make both purchasing and selling decisions while managing the inventory, which makes both the prophet benchmark and the analysis of online algorithms substantially more intricate. For arbitrary numbers of buyers and sellers, we propose a simple threshold-based algorithm and prove a distribution-free competitive ratio of~2, which is the best possible. Our main technical contribution is an exact characterization of the inventory process induced by threshold trading. By expressing the algorithm's expected GFT in terms of the cumulative holding probability at buyer arrivals, we derive an exact formula for this inventory term, which serves as the foundation of our analysis and yields sharper guarantees in several important special cases. For balanced trading (with $n$ sellers and $n$ buyers), we improve the competitive ratio to $(2n^2-n)/((n+4^{-n}-1)(n+1))$, which is strictly less than 2 for every $n>1$. For the single-seller case, we sharpen the analysis of the fixed-threshold algorithm to obtain an asymptotic competitive ratio of $2e^2/(e^2+1)\approx1.76$. Furthermore, by exploiting the additional structure of this setting, we design an adaptive-threshold algorithm whose competitive ratio approaches $e/(e-1)\approx 1.58$. Due to a symmetry property of our general algorithmic idea, the same $1.76$-competitive guarantee also holds for the single-buyer case.

We consider stochastic nonlinear optimization problems with deterministic equality constraints. While unconstrained stochastic optimization is well understood, the interplay between optimality and feasibility in the constrained setting poses significant challenges. Moreover, existing theoretical guarantees for constrained stochastic methods predominantly rely on bounded-variance assumptions, leaving the heavy-tailed noise regime largely unexplored. To address this gap, we propose a novel trust-region method within the stochastic sequential quadratic programming framework, termed TR-SSQP. Our method employs a normal-tangential decomposition in the step computation to balance optimality and feasibility. In addition, we incorporate a normalization mechanism in the design of the trust-region radius, together with Polyak momentum for gradient estimation, ensuring stable updates without gradient clipping. When the trust-region radius and the momentum parameter decay at appropriate rates, we establish global almost-sure convergence of the method. To the best of our knowledge, this is the first asymptotic convergence result for constrained stochastic optimization under heavy-tailed noise. We demonstrate the promising performance of the proposed method through extensive numerical experiments, including comparisons among its variants and with existing constrained stochastic optimization methods.

The integration of large-scale inverter-based resources (IBRs) into the grid presents two key challenges in terms of optimal power generation and transmission system stability. This work proposes concurrent transmission parameter learning and optimal power generation and current control for the grid-following IBRs. The proposed approach estimates transmission parameters online using the current measurements, enabling the IBRs to adapt their control actions to evolving grid conditions. The learned parameter estimates are incorporated into an optimization-based power distribution strategy to determine the optimal power setpoints. In parallel, a current control scheme is designed to regulate the grid-following (GFL) IBRs current output based on the estimates and received active and reactive power references from the high-level optimizer, thereby improving the dynamic response to the grid events. The effectiveness of the proposed framework is validated through a simulation in a MATLAB-Simulink environment. From the results, the impact of the concurrent estimator and current controller during the grid events can be seen.

This paper considers convex envelopes of disjunctive submodular functions—functions that are lattice family submodular over faces of a hypercube—and constructs the first strongly polynomial algorithm for their separation when there are two facial disjunctions. Submodular functions, whose convex envelopes are characterized by the Lovász extension, have occupied a fundamental role in constructing relaxations for combinatorial and nonlinear optimization problems. However, disjunctive submodular function envelopes have not been explored besides the use of ellipsoid algorithm, which remains practically intractable. Our algorithm is derived in three steps by expressing the disjunctive function as a minimum of two extended submodular functions, introducing a variable lifting technique, and constructing the sublinear envelope in the lifted space. The paper also makes several other contributions. First, we provide a disjunctive formulation for the case where each submodular function admits a linear programming formulation. Second, we derive the closed-form sublinear envelope characterization for intersecting submodular functions, yielding new structural insights into a multi-product inventory sales maximization problem. Third, we fully characterize the convex envelope of a bilinear function defined over a cycle graph in the original variable space. Finally, we show computationally that the cycle inequalities close approximately 60% of the gap for complete and Hadamard graphs, over 30% of the gap for complete bipartite graphs, and over 80% of the gap for sparse graphs such as cactus and Halin graphs. The resulting relaxations are also more efficient to solve than previous extended space formulations.

We study an elliptic optimal control problem governed by Darcy's equation in heterogeneous porous media, with pointwise box constraints on the control. The objective functional is formulated via a Radon measure, which allows the desired pressure state to be tracked on observation sets of varying dimension, including points, curves, and subdomains, within a single formulation. The state and adjoint equations are discretized by a symmetric interior penalty discontinuous Galerkin method, yielding a locally mass-conservative approximation that is robust across strong permeability discontinuities, while the control is approximated by piecewise constants. The state, adjoint, and control are retained as primary unknowns in a single coupled optimality system and solved monolithically by a primal-dual active set strategy. We establish stability and well-posedness of the discretization and derive a priori $L^2$ error estimates for both variables. The principal difficulty is the reduced regularity of the adjoint state induced by the measure-valued tracking data. The analysis controls the resulting adjoint-control coupling through an intermediate adjoint driven by the continuous optimal state. Numerical experiments confirm the predicted convergence rates for point, curve, and subdomain observations, exhibit mesh-independent primal-dual active set iteration counts, and demonstrate local mass conservation.

This paper studies constraint certification in synchronous distributed model predictive control (DMPC) when neighboring predictions change between sampling instants. Before the parallel local solves, each agent communicates a shifted prediction and an announced update budget. A hard trajectory trust region makes that budget enforceable, while an edge-wise feasibility cap computed from the shifted packets keeps the fallback feasible without using any current optimizer output. Distance and relative-attitude constraints are tightened with explicit Lipschitz constants and two budget layers: one accounts for the simultaneous neighbor update and the other retains a checkable shift reserve. We prove hard pairwise constraint satisfaction and recursive feasibility under stated nominal-execution and terminal assumptions, give the additional residual caused by execution error, and derive a local practical value-decrease bound. A spacecraft formation example uses hard terminal and pairwise constraints, a geodesic relative- attitude constraint on $\SO$, and reproducible terminal-set checks. Comparisons with fixed, trajectory-only, and windowed online margins show that the proposed budget reduces conservatism while preserving a positive shifted-feasibility margin.

We determine the minimax dimension dependence of noisy zeroth-order optimization for globally strongly convex functions with higher-order smoothness in multilinear operator norm. Each query returns one function value with fresh independent Gaussian noise. For dimension $d$, query budget $T$, and any fixed finite smoothness order $\beta>2$, the minimax expected objective error satisfies $\mathcal{E}_{\beta}(T,d)=\Theta_{\beta}\left(\min\left\{1,(d^2/T)^{(\beta-1)/\beta}\right\}\right)$, with dimension-independent strong-convexity, gradient-smoothness, higher-order-smoothness, noise, and minimizer-radius bounds. Thus, for sufficiently small target error $\varepsilon$, the necessary and sufficient number of noisy values is $\Theta_{\beta}(d^2\varepsilon^{-\beta/(\beta-1)})$, with no logarithmic loss. The new lower bound matches the dimension and budget dependence of the upper bound of Akhavan et al. (2024), and holds for arbitrary randomized adaptive algorithms with unrestricted query points. Its key step localizes a normalized hypercube with one radial cutoff: a query can emphasize one hidden coordinate, but the aggregate information about all coordinates remains small. A self-contained multiscale spherical estimator attains the matching upper bound by cancelling lower-order bias terms. Constants in the rates depend only on the fixed smoothness order under the stated normalization.

We prove a joint stochastic zeroth-order lower bound for Goldstein stationarity on a Euclidean query ball, even when the ball is guaranteed to contain a stationary point. In dimension $d$, let $f=\mathbb{E}[F(\cdot;\xi)]$, assume $\mathbb{E}[\operatorname{Lip}(F(\cdot;\xi))^2]\le L_0^2$, and bound the initial objective gap over the ball by $\Delta$. For neighborhood radius $\delta>0$ and residual tolerance $\varepsilon>0$, our smooth hard family requires $\Omega(dL_0^2\Delta/(\delta\varepsilon^3))$ scalar evaluations for success probability $1/2$, against randomized adaptive algorithms that may retain and repeatedly query each sampled function. The result holds for $\varepsilon\le cL_0$, $\Delta\ge C\delta\varepsilon$, and $d\ge C[1+\log(2+\Delta L_0^2/(\delta\varepsilon^3))]$, on a constructed ball of radius $\Theta(\Delta/\varepsilon)$. A sequence of localized regions forces repeated direction estimation, and an adaptive Gaussian posterior argument controls sample reuse. The result establishes a joint dimension and accuracy obstruction on solvable bounded-domain instances. The gap is local to the query ball; a matching minimax characterization at a common radius and the corresponding unrestricted global-gap lower bound remain open here.

We establish sharp high-probability statistical rates for nonsmooth stochastic convex optimization under heavy-tailed noise, closing the confidence-dependence gap in both general convex and strongly convex settings. We consider composite objectives with a Lipschitz convex loss and a known proximal regularizer, accessed through fresh unbiased stochastic subgradients. The centered noise has conditional norm and directional $p$th moments bounded by $L^p$ and $s^p$, respectively, where $1<p\le2$ and $0<s\le L$. For $T$ oracle calls and failure probability $0<\delta<1/10$, set $q=1-1/p$, $\ell=\log(3/\delta)$, and $A=s(L^2/s^2)^q$. With initial distance bound $D$, the minimax objective error is $\Theta(D(A+s\ell^q)/T^q)$ in the nonsaturated statistical regime. For a $\mu$-strongly convex regularizer, the corresponding quadratic-regime rate is $\Theta((A^2+s^2\ell^{2q})/(\mu T^{2q}))$. These rates require domination of the deterministic and finite-variance terms; matching lower bounds require $d\ge\lceil L^2/s^2\rceil$, while the upper bounds hold in every finite dimension. The new upper bounds attain known statistical lower-bound orders, replacing the previous confidence coefficient $s^{1/p}L^{1-1/p}$ by $s$, and its square by $s^2$ under strong convexity. The key is an online inequality that separates directional fluctuations from the total second moment created by rejection sampling. Localization preserves this separation under strong convexity without additional horizon logarithms. The results characterize oracle complexity with unrestricted internal computation.

We establish matching polynomial query bounds for low-rank approximation from exact matrix–vector products. Given an unknown matrix $A\in\mathbb{R}^{m\times n}$, at each step a randomized algorithm chooses either $v\in\mathbb{R}^n$ and receives $Av$, or $u\in\mathbb{R}^m$ and receives $A^\top u$. The choice may depend measurably on all previous queries and replies and on the algorithm's private randomness; each vector product costs one query. The output is a rank-$k$ right projector with Schatten-$p$ residual at most $1+\varepsilon$ times optimal. Write $N=\min\{m,n\}$ and let $Q_p^*$ denote the worst-case query budget for success probability $2/3$ on every input. For every $1\le k<N$ and sufficiently small $\varepsilon$, our lower bounds, combined with existing Krylov upper bounds, give $Q_p^*=\widetilde{\Theta}\!\left(\min\{N,k\min\{p^{1/6}\varepsilon^{-1/3},\varepsilon^{-1/2}\}\}\right)$ $(2\le p<\infty)$, $Q_\infty^*=\widetilde{\Theta}\!\left(\min\{N,k\varepsilon^{-1/2}\}\right)$. These bounds have universal constants and allow $p$ to vary with the problem parameters, identifying the transition at $p\varepsilon\asymp1$. A complementary result for each fixed $1\le p<2$ gives $\widetilde{\Theta}_p(\min\{N,k\varepsilon^{-1/3}\})$, with constants and an accuracy threshold that may depend on $p$. Together, the results recover this fixed-norm rate for every fixed finite $p$, supplying the multiplicative rank dependence missing from previous lower bounds. Tildes suppress logarithmic factors. The proof extends adaptive Wishart deferred decisions to a rectangular factor with a $k$-dimensional nullspace. Posterior overlap gives a short fixed-norm argument, while persistence of small compression eigenvalues controls growing $p$ and the spectral endpoint. Exact range recovery handles target costs of order $k$; the Wishart family covers the remaining regimes.

Machine-learning policies have attracted increasing interest for ride-hailing fleet control in recent years. Reinforcement learning, in particular, requires a structured simulation environment that specifies observations, actions, rewards, and decision epochs for training and evaluation. For electric fleets, this environment must also capture the interaction among stochastic demand, vehicle operations, and capacitated charging infrastructure. We present OpenHail, an open-source Gymnasium environment for joint control of electric ride-hailing fleets. Its fixed-size observation–action interface exposes request assignment, repositioning, and charging to a single policy. The event-driven simulator represents requests with pickup deadlines, vehicle job queues, battery dynamics, and finite-capacity charging facilities with first-in–first-out queues. A configurable decision-epoch mechanism separates internal simulator events from policy interactions, supporting event-driven, periodic, hybrid, and policy-requested control within the same operational model. The software provides seeded instances, feasible-action utilities, evaluation tools, operational metrics, and baseline policies. The source code is available at https://github.com/tommaso-schettini/openhail.

In this paper, we establish the first asymptotic convergence guarantees for the Muon algorithm through a more accurate proxy for the Newton-Schultz iteration than the typical matrix sign function. We prove that, for appropriate choices of hyperparameters, the iterates satisfy $\lim_{k\to\infty}\|\nabla f(x_k)\|=0$, and, under a global Polyak-Łojasiewicz condition, that the sequence of function values converges linearly. The key insight is that the regularization, implicit in Muon's Newton-Schulz implementation, induces a bounded preconditioner, exposing Muon as a preconditioned Polyak heavy-ball method and enabling a classical Lyapunov analysis. This observation naturally motivates applying the same preconditioning structure to the Nesterov gradient evaluation. We formalize this idea by introducing Muesterov, a Nesterov-based variant of Muon, and prove that it enjoys the same convergence guarantees, extending the theoretical framework beyond the heavy-ball setting. Numerical experiments on a scalar cross-entropy problem corroborate the theory and illuminate the joint role of the learning rate and the Newton-Schulz regularizer in controlling convergence. Preliminary numerical simulations training the nanoGPT dataset provide intuition regarding the relevance of the observations in this paper to practical applications.

Parameter estimation for ordinary differential equation (ODE) models is a fundamental task that is often complicated by the limitations of conventional optimization-based methods. In theory, differential-algebraic approaches offer an appealing alternative: they reduce the problem to polynomial system solving and do not require user-supplied initial guesses for parameter values. In practice, however, algebraic methods have been limited by their sensitivity to measurement noise, because they require accurate derivatives of observed outputs. In this work, we integrate Gaussian Process Regression (GPR) into the differential-algebraic method and derive a first-order error analysis in terms of noise level and algebraic sensitivity. We evaluate the method across several noise levels on a benchmark of 25 dynamical systems arising in applications including mechanical engineering and systems biology. The proposed method achieves the highest aggregate performance among the methods considered, recovering all sought parameter values and initial conditions to within 10% relative error in 88.5% of runs. These results demonstrate that robust derivative estimation can make differential-algebraic parameter estimation practical for dense, noisy synthetic data while retaining key advantages of the algebraic formulation.

In this paper, we consider linear stochastic systems with control-linear noise and periodically sampled measurements. We address the problem of feedback stabilization in the mean-square sense. The main contribution of the paper is to provide a necessary and sufficient condition for feedback stabilization. In particular, we relate feedback stabilizability of the stochastic system to the existence of a fixed point of a Riccati-type algebraic equation and, further, to the existence of a solution to an infinite-horizon optimal control problem.

In this paper, we derive a closed-form for the proximity operator of the weighted squared $\ell_{2,\infty}$ norm. Moreover, we derive a closed-form for the proximity operator of a weakly convex version of this function. The proof involves using known results for compute the proximity operator of a supremum function, which requires finding a solution of an auxiliary problem. We find the explicit solution of this auxiliary problem by finding the KKT multipliers and the critical point associated to the first-order optimality conditions. Additionally, we present applications to denoising problems and weighted max-min dispersion problems.

This paper investigates the numerical resolution of the Beurling-LASSO (BLASSO), a convex optimization framework that promotes sparsity in the space of measures. We consider its application to the estimation of Gaussian mixture models (GMMs) with an unknown number of components and unknown diagonal covariance matrices. Our approach combines the Conic Particle Gradient Descent (CPGD) principle with Riemannian gradient descent, to account for the underlying Fisher-Rao geometry of Gaussian distributions. Our contributions are twofold. First, we provide theoretical guarantees for the convergence of our algorithm. In particular, we establish exponential local convergence under a non-degeneracy condition on the solution and relate this assumption to a separation condition on the underlying statistical target. Second, we address practical implementation aspects of CPGD and present numerical experiments illustrating its performance. On the test cases considered, these experiments suggest that CPGD is more robust to overspecification of the number of components than the EM algorithm. We also investigate the impact of component separation on recovery accuracy.

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.

We study the deterministic oracle complexity of finding approximate solutions to composite monotone inclusion problems, formed by the sum of a smooth single-valued monotone operator and a maximally monotone set-valued operator, under the tangent-residual criterion. We introduce the Anchored Extra-Proximal (AEP) framework, which combines an anchored extrapolation step with an inexact anchored proximal update satisfying a relative-error condition. The framework recovers the composite Fast Extragradient method in the first-order setting and yields natural second- and higher-order extensions by replacing the operator in the implicit update with its Taylor approximation at the extrapolated point. For every $p\geq 2$, assuming that the $(p-1)$th derivative of the single-valued operator is Lipschitz continuous, we combine this construction with a bisection line search to obtain a $p$th-order method that finds a point with tangent residual at most $\varepsilon$ in $\widetilde{O}(\varepsilon^{-2/(3p-1)})$ oracle calls. This improves all prior upper bounds for $p$th-order methods: in particular, it improves the previous best-known $\widetilde{O}(\varepsilon^{-1/p})$ tangent-residual complexity as well as the classical $O(\varepsilon^{-2/(p+1)})$ bound of higher-order hybrid proximal extragradient methods under the weaker duality-gap criterion. We complement this result with a worst-case lower bound of $\Omega(\varepsilon^{-2/(3p-1)})$ for every deterministic algorithm in the $p$th-order oracle model, without restricting the algorithm to tensor steps or any other prescribed update structure. Thus, the proposed method attains the optimal dependence on $\varepsilon$, up to logarithmic factors, for all $p\geq2$.

We introduce GridSFM, a framework that combines a pretrained foundation model across grid topologies with physics-informed fine-tuning for solving AC Optimal Power Flow (AC-OPF) at scale. It is a $15$ million parameter physics-inspired graph neural network pretrained across $54$ topologies of $500$ to $4{,}000$ buses. Our model attains a $2.45\%$ zero-shot generation-cost error on a $10{,}000$ bus case held-out operating conditions with no degradation as system size grows. Building on this, we pair the pretrained backbone with a physics-informed fine-tuning design based on Newton's method for power flow. With only $100$ solved instances, GridSFM adapts to unseen grids up to $10{,}000$ buses. We show it out performs single topology, dedicated neural network models that are trained more data, both in terms of cost and solver iterations when deployed as warm starting points. In designing this foundation model, we overcome the fact that the feasible set for AC-OPF can be disconnected. This is an obstruction that prevents any continuous neural network from approximating the solution map. To do so, we lift the problem and relax its constraints with logarithmically penalized slacks. We prove that the resulting elastic feasible set is contractible, that the AC-OPF minimizers remain minimizers of the elastic problem above an explicit penalty threshold, and that projecting an approximate solution back onto the AC-OPF feasible set is well posed. We release all models, data, and code so that the community can build on a shared starting point for AC-OPF.

High-dimensional extreme eigenvalue problems often arise from molecular vibrational models, electronic structure calculations, and quantum mechanics. Directly solving these problems suffers from the curse of dimensionality. Instead of tackling the problem in high-dimensional ambient space, we reformulate the problem through low-rank tensor formats, which can significantly reduce the computational cost and storage. Specifically, we consider the Rayleigh–Ritz problem on bounded-rank tensors in the tensor train format, which enables a more flexible choice of rank parameters. Moreover, we consider a smooth parametrization for bounded-rank tensors by introducing slack variables, leading to a smooth manifold structure. The original Rayleigh–Ritz problem is therefore transferred to an optimization problem on the manifold. We develop the Riemannian geometry and propose optimization methods for solving extreme eigenvalue problems. In practice, the Kronecker-product structure in Hamiltonian and PDE operators is employed to simplify the computation. Numerical experiments on the harmonic oscillator, the Laplace operator, the Schrödinger equation, the layered cluster problem, and the Bose–Einstein model demonstrate that the proposed method achieves accuracy comparable to full-space eigensolvers with reduced computation time and avoids storing full vectors. The results also show numerical rank reduction during iteration and robustness when the rank parameter is over-estimated. In addition, the fixed-point iterations illustrate that the proposed methods are able to serve as an inner eigensolver for solving nonlinear eigenvalue problems.

The main purpose of this paper is to establish the local quadratic convergence of Newton's method for nonlinear optimal control problems. Under suitable smoothness assumptions, the cost functional is first shown to be twice Gâteaux differentiable with respect to the control. An explicit operator representation of its second derivative is derived, and its continuity in different control spaces is examined. The second-order coercivity condition in $L^2(0,T;\maathbb{R}^m)$ is then shown to be equivalent to the existence of a strongly regular solution to an associated matrix Riccati differential equation, thereby expressing an infinite-dimensional quadratic-form condition in terms of the solvability of a matrix differential equation. On this basis, local quadratic convergence of Newton's method in $L^\infty(0,T;\maathbb{R}^m)$ is established using linear-quadratic optimal control theory. Under suitable structural assumptions, local quadratic convergence in $L^2(0,T;\maathbb{R}^m)$ is also established. Explicit estimates of the corresponding convergence neighborhoods are derived, and the theoretical results are illustrated by numerical examples.

Graph-structured optimization with linear constraints is fundamental to critical infrastructure but faces scalability limits due to massive strict hard constraints and high dimensionality. While recent projection-based methods such as Trainable Sampling Kaczmarz-Motzkin Net (T-SKM-Net) guarantee feasibility, they face high computational costs in dynamic environments by processing the entire constraint set and requiring expensive matrix factorizations. To bridge this gap, we propose the Accelerated Trainable-SKM (AT-SKM) Net framework. To concentrate computation on the active constraints and eliminate redundant calculations, we introduce a hybrid sampling strategy guided by a topology-aware heterogeneous GNN model. To efficiently handle topological shifts in graph-based constraints, we employ a Cholesky Update mechanism that theoretically reduces the equality projection complexity from O(N^3) to O(N^2) under low-rank perturbations. Experiments on random geometric graphs, N-1 Security-Constrained DC-OPF, and minimum-cost gas transport problem demonstrate that AT-SKM reduces iteration counts by up to 85% and achieves 2.95x-7.29x SKM layer speedups, while maintaining zero constraint violations.

We study the computational complexity of peak-oriented rationalizability, a survey based revealed-preference test introduced by Seror (2026). We provide a polynomial-time algorithm for testing rationalizability and recovering a utility function, establishing that peak-oriented preference elicitation is computationally tractable. In contrast, we show that computing the peak-oriented Houtman-Maks index is NP-hard. These results delineate the precise computational boundaries of peak-oriented revealed-preference analysis.