Papers

We study the query complexity of finding approximate stationary points under higher-order smoothness using only function values or only gradients. Our results cover different dimensions, oracle models, randomization, and constraints. For constant dimension $d$ and $C^{p,\beta_p}$ objectives, with $\nu:=p+\beta_p$, we give deterministic function-value gradient trapping algorithms with query complexities $O(\varepsilon^{-(d-1)})$ on $\mathbb R^d$ and $O(\varepsilon^{-(d-1)/\nu})$ on $[0,1]^d$. High-order line integration enables gradient-only trapping. For two-dimensional $C^{1,1}$ objectives, randomized integration yields $O(\varepsilon^{-4/3})$ gradient queries on $\mathbb R^2$, breaking the deterministic gradient-only $\Theta(\varepsilon^{-2})$ barrier. For general dimension, we introduce a trapping-or-descent framework that lazily reuses boundary models. For $C^{p,\beta_p}$ objectives with $p\ge2$, it achieves deterministic query complexities $O(d^p+d^{p-1+\beta_p/(\nu-1)} \varepsilon^{-\nu/(\nu-1)})$ using function values and $O(d^{p-1}+d^{p-2+\beta_p/(\nu-1)} \varepsilon^{-\nu/(\nu-1)})$ using gradients. These bounds apply to approximate stationary points on $\mathbb R^d$ and approximate KKT points on $[0,1]^d$. Their dimension-independent accuracy exponent matches the classical AR$p$ exponent $(p+1)/p$ when $\beta_p=1$. For sufficiently high dimension, we prove randomized lower bounds for gradient-only methods. For $C^{1,1}\cap C^{p,1}$ objectives, our lower bound is $\Omega(\varepsilon^{-(3/2+1/(2p))})$, matching known upper bounds for $p=2$ and, up to logarithmic factors, for $p=3$. Without the $C^{1,1}$ assumption, we obtain an $\Omega(\varepsilon^{-5/3})$ lower bound for every fixed $p\ge3$.

Multiparametric optimal control in continuous time partitions the space of initial states into critical regions, each with its own control law involving a set of time-varying active constraints. The boundaries between two neighbouring regions are either hyperplanes, which are given in closed form by one linear equation, or curved boundaries, which must be computed. The sequence of time-varying active constraints plays a key role in defining the linear or nonlinear shape of a critical-region boundary. In this paper we show how to determine that shape prior to any computation. In particular, we show that the boundary is a hyperplane under two conditions. First, the new constraint becomes active at the same instant at every point of the boundary. Second, the durations for which the other constraints stay active are the same at every point of the boundary. Otherwise, the boundary is curved. Both conditions can be read from the sequences of active constraints of the two neighbouring regions, which gives five rules, each yielding a hyperplane or a curved boundary. We illustrate the rules on the Newell-Lee evaporator, where they classify the eight boundaries of a seven-region partition into four hyperplanes, given in closed form, and four curved boundaries, computed exactly.

Modern state space models (SSMs), such as Mamba2, provide a compelling alternative to transformers by combining linear-time sequence modeling with recurrent state-space dynamics. However, the behavior of SSMs in distributed learning settings remains poorly understood. In particular, the existing standard federated learning methods are largely architecture-agnostic, and do not account for the stability, selectivity, and state-space parameterization that characterize modern selective SSMs. To address this, we derive architecture-aware gradient and smoothness bounds for single- and multi-layer selective SSMs, and convergence bounds for FedAvg and FedProx, characterizing how recurrent stability, input-dependent discretization, and state projection norms affect federated optimization. We then numerically validate the single-layer bounds on sequences generated by a teacher SSM, using a learner that follows the analyzed recurrence. We use this analysis to formulate expectations about the effects of local training and client heterogeneity, and examine these expectations by comparing nine federated learning algorithms on Mamba2 language modeling across six text domains. These experiments illustrate how SSM-specific bounds can provide a basis for interpreting the behavior of practical federated learning algorithms.

Magnetic resonance elastography (MRE) is a noninvasive imaging modality for quantifying the viscoelastic properties of soft tissues from shear wave propagation. Recovering the complex-valued shear modulus from measured displacement fields leads to a severely ill-posed inverse problem, particularly in the presence of noise and limited boundary excitations. We investigate adjoint-based optimization, convolutional neural network (CNN) reconstruction, and a hybrid framework combining both approaches. The forward model is based on a scalar form of the modified stationary Stokes system with a complex shear modulus. We establish well-posedness of the forward problem, existence of minimizers, and first-order optimality conditions for the adjoint-based formulation, and implement a nonlinear conjugate-gradient method with Armijo line search. While PDE-constrained optimization can accurately refine coefficient reconstructions, its performance depends strongly on initialization. We therefore construct a two-dimensional CNN that maps complex-valued displacement measurements to spatially varying complex shear modulus fields and provides rapid, informative initial reconstructions. The proposed hybrid method uses the CNN reconstruction to initialize the adjoint-based optimization, yielding faster convergence and improved accuracy. The CNN is trained on coefficient fields containing individual perturbations and tested on both individual and previously unseen combined configurations. Numerical experiments demonstrate that the CNN generalizes to these more challenging configurations, while subsequent PDE-constrained optimization further refines the reconstructed coefficient. These results demonstrate the potential of combining data-driven initialization with physics-based optimization for efficient and accurate MRE reconstruction.

This work proposes a data-driven framework that synthesizes output-feedback control for linear time-invariant multi-agent systems (MAS). We first develop two algorithms that use offline data samples of state and input trajectories to determine a Nash equilibrium (NE) for MAS. Next, we formulate a semidefinite program that, given the NE solution, computes the static output-feedback controllers for the agents in MAS. Compared to prior work, virtues of our framework include weaker controllability assumptions, the absence of online information and data-driven estimation of states. Next, we use the proposed framework to design control schemes for voltage and frequency regulation for nonlinear multi-nodal electrical networks with grid-forming and grid-following inverters. Thus, compared to prior work that implements their data-driven frameworks in small-scale numerical examples, we validate the performance of our data-driven framework in high-fidelity nonlinear models of MAS. This implementation in realistic MAS provides insights that small-scale numerical examples may fail to capture. Results show successful grid stabilization and provision of regulation services from the inverters in the grid using our proposed framework. It can attain up to 78% less overshoot and 30% faster settling time than traditional strategies such as proportional-integral control.

This work proposes a semidefinite programming (SDP) model that synthesizes static output control feedback for networked linear time-invariant multi-agent systems (MAS). We theoretically prove that the output feedback control of MAS can be constructed using the Nash equilibrium previously derived for state-feedback control. The Nash equilibrium considers that agents seek to minimize their individual cost functions subject to the system dynamics of the MAS. The resulting control from the pro- posed SDP model ensures the stability of the closed-loop system and dissipativity when a disturbance is present in the MAS. A set of four case studies is implemented. The first three numerical cases serve to validate and compare cost effectiveness and time-domain performance of the proposed control framework with other state-of-the-art frameworks. The fourth case designs a frequency regulation scheme for a modified IEEE 24-bus power grid. Using the proposed framework, distributed controllers are synthesized for the inverters which are considered agents in the grid. Each inverter uses partial state information of the neighboring inverters. The resulting control scheme is implemented in a high-fidelity non-linear model of 261 states that describes the electromagnetic transient dynamics of the power grid. Hence, compared to prior studies that usually validate their frameworks in small-scale numerical models, this fourth case demonstrates a virtue of our pro- posed SDP: its scalability and practical utility in real-world networked MAS. We implement the proposed framework in an open-access toolbox in Python and MATLAB, enabling practitioners and researchers to readily implement it in their applications.

In this paper, we consider the control of Markov jump linear systems when the active system mode is not exactly known. We rather assume that this true mode is only known to belong to an observed cluster of modes, making both the dynamics and the associated transition probabilities uncertain. We construct a min-max optimal control problem in this setting, with the objective of robust regulation. Using results in regularized least mean-square optimization, we describe how the solution of this problem can be reduced to a Riccati recursion, such that the controller depends only on knowledge of which cluster the mode lies in. We subsequently establish convergence of the recursion and stability of the resulting closed-loop system. A numerical example illustrates the effectiveness of the proposed methodology.

Motivated by scenarios in control systems where decision making is informed by only partial information about the states, we introduce a novel formulation of the Schrödinger bridge problem for inertial interacting-agent systems. Within our setting, the agent ensemble is steered to match partial specifications of the collective configuration. Our analysis pertains to a mean-field framework in which the collective dynamics are described by a nonlinear Vlasov–Fokker–Planck integro-PDE with nonlocal interactions. We assume that only distributions over positions are specified at the endpoints of a finite time horizon, representing aggregate population snapshots, while the corresponding velocity degrees of freedom remain unspecified. We derive the associated optimality system, identifying the optimal controlled evolution as jointly informed by the initial and terminal specifications through forward and backward propagation in time. We also show that the partial endpoint specifications bring about distinct conditions on the optimal control and the consistent evolution. Subsequently, to compute solutions to the optimality system, we develop a nested fixed-point algorithm inspired by the classical Fortet–Sinkhorn iteration, but extended to handle the nonlinearities arising from nonlocal interactions and temporal coupling.

We study finite-horizon stochastic linear-quadratic (LQ) control for finite-activity jump diffusions with bounded random coefficients, indefinite weights and a linear recursive cost. The cost is defined by a backward stochastic differential equation (BSDE) whose generator has bounded deterministic coefficients multiplying the Brownian and jump integrands. A Brownian-Poisson change of measure and an integrating factor reduce the criterion to a nonrecursive LQ functional. On the intersection of the physical and transformed square-integrable control spaces, we establish recursive evaluation for L^1 quadratic data under the transformed measure and equality with the full transformed value. Under transformed uniform convexity and invertibility of the uncontrolled jump map, we identify the stochastic Riccati equation from the nonrecursive value kernel. Conditional energy estimates justify transferring its martingale coefficients to the physical measure. Completion of squares gives feedback representation, comparison and uniqueness in the feedback-admissible class. For each initial pair, the recursive value is attained exactly when the transformed optimizer has finite physical control energy. Analytical examples illustrate nonzero feedback and uniform convexity with a negative control weight.

We introduce an autoregressive differentiable method to solve 0-1 integer programs. We fix an arbitrary order of the binary variables and we train a transformer to predict the next bit while remaining in the feasible set. Our method is first trained on feasible incumbents provided by any solver, thus allowing us to initialize the transformer in the feasible set. Our procedure then implements a Lagrangian penalty to penalize infeasible solutions, and the transformer is further trained to explore the feasible set using Gumbel-softmax activations on the relaxed objective. We have tested our method on non-convex instances of quadratic knapsack problem and demonstrated consistent improvement upon state-of-the-art open-source solvers for dense problems up to 10,000 binary variables. In particular, we empirically demonstrate a phenomenon akin to a tunneling effect where the effective change of variables from binary variable to the continuous weights of the transformer that the method implements enables crossing barriers in the relaxed objective landscape.

This paper investigates the connection between mean-field systems and forward-backward stochastic systems arising from the mean-variance portfolio selection problem. By introducing the conditional expectation process of terminal wealth, the mean-field term is represented by a backward state, leading to a class of control problems for forward-backward stochastic systems with coupled initial-terminal values. The coupling between the terminal value of the forward state and the initial value of the backward state in the cost functional prevents classical methods for forward-backward stochastic optimal control from being directly applicable. To overcome this difficulty, we establish a stochastic maximum principle for control problems with running costs and coupled initial-terminal values, and provide sufficient conditions for optimality. Furthermore, we investigate a class of generalized stochastic linear-quadratic control problems with initial-terminal coupling, allowing the control weight in the running cost to be indefinite, and derive the optimal feedback representation via a Riccati equation. The classical mean-variance portfolio selection model is shown to be a special case, for which the corresponding efficient investment strategy and efficient frontier are obtained explicitly, and numerical simulations are presented to illustrate the results and examine parameter sensitivity and its financial implications.

This paper establishes sufficient conditions for finite-horizon solvability of fully coupled nonlinear mean-field forward–backward stochastic difference equations with dependence on unconditional first moments. The forward recursion uses conditional projections of the next backward state and its product with the innovation. Building on existing domination–monotonicity methods, we formulate deterministic matrix combinations in centered and mean coordinates and prove a continuation estimate uniform in the homotopy parameter. Global Lipschitz continuity and one active domination–coercivity direction yield a unique square-integrable adapted solution, an a priori bound, and coefficient stability. The active parameter can be normalized without imposing an additional smallness restriction on the original coefficients. A sign transformation handles the opposite monotonicity orientation. A nonlinear example with saturating state and mean interactions verifies the assumptions, including degenerate domination directions. A scalar mean-field LQ buffer-adjustment model illustrates the theorem: its Hamiltonian solution characterizes the unique open-loop optimizer. An exact finite scenario-tree calculation checks this characterization against direct quadratic optimization.

We present sparse dictionary differentiable predictive control (SD-DPC), a framework for learning sparse, interpretable feedback policies for nonlinear systems from data. A prediction model is first identified by rollout-based sparse identification of nonlinear dynamics (SINDy), building on gradient-based and multistep formulations. The policy is then parameterized as a sparse combination of dictionary functions and trained by differentiating a constrained finite-horizon predictive-control objective through this model, so that its terms are selected by closed-loop performance rather than by imitating a previously trained controller. The result is an explicit feedback law with only a handful of terms. Across three benchmark control problems, SD-DPC satisfies the constraints in all test scenarios, outperforms a policy distilled onto the same terms by up to an order of magnitude, and requires orders of magnitude less memory and online computation than an optimization benchmark, while admitting explicit sensitivity bounds.

This letter introduces an accelerated safe gradient flow (A-SGF) for strongly convex constrained optimization problems. The A-SGF filters a nominal second-order accelerated flow through a second-order control barrier function quadratic program. We establish conditions for strict feasibility and local Lipschitz continuity of the resulting vector field, prove forward invariance of a lifted safe set, and show that its unique equilibrium is the constrained minimizer. A variational Lyapunov analysis isolates the effect of the safety correction and gives explicit conditions on the barrier gains under which A-SGF retains the same convergence rate as the unconstrained Nesterov accelerated method. A nonlinear 2D example demonstrate anytime feasibility and faster convergence of the algorithm.

Since the 1980s, risk analysis for engineering systems has been defined based on a tuple of scenarios, $\mathcal{S}$, likelihoods, $L$, and consequences, $C$, where $Risk = (\mathcal{S},L,C)$. We argue that equating this triplet with risk alone is insufficient and potentially misleading. Instead, we define this tuple as a vulnerability space $\mathcal{V} = (\mathcal{S},L,C)$, which becomes risk only when one makes specific interpretations of likelihood and consequence. Here, an interpretation has mathematical meaning, and defines a map that assigns data elements about system vulnerability to a perspective. Formally, we define a vulnerability perspective as a specific choice of interpretations of likelihood and consequence data that allows vulnerability analysis measures to be defined and analyzed. We argue that risk is only one of four common perspectives on $\mathcal{V}$ and demonstrate mathematically why each perspective should remain a distinct concept. First, we show how interpretation of $C$ as a set-valued function produces reliability measures, not risk measures. Then, we show how interpretation of $L$ via possibility theory will produce adversarial and safety measures distinct from risk. This new theoretical framework clarifies the conceptual, mathematical, and practical differences between risk, reliability, adversarial, and safety perspectives—four common engineering system analyses often conflated in the literature. A consequence of our theory is that no single perspective will reveal all vulnerabilities of an engineering system, suggesting decision-making via risk analysis may incur unavoidable and underappreciated tradeoffs from other vulnerability perspectives.

Increasing the batch size during training is a common practice in large language model (LLM) pretraining, yet the theoretical justification behind its success is not well understood. Analyses of stochastic optimization often assume uniformly bounded stochastic gradient variance, yet recent evidence suggests that this assumption fails in many practical nonconvex problems. The Blum–Gladyshev (BG-$0$) noise model relaxes this assumption by allowing the variance to grow quadratically with the distance from initialization, suggesting that batch size schedulers can help by controlling the variance growth during training. However, this growth can be overly conservative in practice. We empirically investigate variance growth in LLM pretraining and observe that a generalized BG model with a tunable growth exponent provides a tighter description of practical noise behavior. Motivated by this observation, we introduce the generalized BG-$a$ noise model, which interpolates between bounded variance ($a=0$) and BG-$0$ noise ($a=2$). Under $L$-smoothness, we derive an information-theoretic lower bound with growth-dependent oracle complexity $\Omega(\epsilon^{-(4+a)})$ and establish a matching upper bound in $\epsilon$-dependence by increasing the batch size as the iterates move away from initialization. Finally, we propose an adaptive batch scheduler that controls variance growth through dynamic batch size adjustments during training. In pretraining OLMo2 models of up to 1B parameters on C4, our scheduler achieves a lower validation loss than both small and large batch training under matched token budgets, while using less than 10% of the iterations of small batch training.

We prove that the maximal cp-rank $p_n$ satisfies $p_n\ge n(n-3)/2$ for odd $n\ge5$ and $p_n\ge n(n-3)/2-1$ for even $n\ge6$. Together with the known upper bound $p_n\le\binom{n+1}{2}-4$, this gives $p_n=n^2/2+O(n)$.

Full-parameter fine-tuning of large language models (LLMs) incurs substantial optimizer state memory overhead, limiting the model sizes that fit on modern GPUs. Existing approaches either compress optimizer state, abandon first-order gradients, or change the update geometry while retaining dense state. The recently introduced Muon optimizer reduces optimizer memory through matrix-valued updates. Still, its geometry differs from AdamW and can lead to performance degradation when fine-tuning AdamW-pretrained models. To reduce optimizer memory without sacrificing accuracy or computational efficiency in LLM fine-tuning, we propose Ternary Absolute-max Column-wise One-sparse optimizer, or TACO, which follows Muon's operator-norm steepest-descent view but takes the geometric route further. TACO computes the exact steepest-descent direction under a dimension-normalized $1\to1$ operator norm by selecting the sign of the largest magnitude entry in each column of two-dimensional weight matrices. This retains first-order gradients while making optimizer state memory nearly negligible. Our practical TACO optimizer maintains only a small set of low precision gradient components per column, reducing persistent optimizer state by $174\times$ relative to AdamW8bit (from 27.7 GB to 0.16 GB) and peak training memory by $2.9\times$ (from 80.6 GB to 27.5 GB) on OPT-13B, while achieving comparable accuracy and runtime. TACO further enables full-parameter fine-tuning of 30-32B-parameter models on a single 80 GB H100 GPU across multiple model families and tasks.

Step-size selection remains a central challenge in large-scale neural network optimization; conservative steps slow convergence, while aggressive steps can destabilize it. We combine Zero-and-First-Order optimization~(ZFO) and propose a lightweight framework that decouples direction selection from step-size. ZFO uses a trusted first-order optimizer to determine the direction and performs zeroth-order evaluations only along this one-dimensional subspace to choose how far to move. Using the current {gradient information} and two additional objective function evaluations, ZFO instances construct a local model of the objective function along the proposed direction and select a curvature-aware step within a bounded search interval. This yields an adaptive step-selection mechanism that costs less than a full line search. We provide theoretical guarantees to show that shared-sample evaluations produce reliable finite-difference curvature estimates, that the induced local model selects a near-optimal step along the search interval, and that ZFO converges to a neighborhood of a stationary point. Across the evaluated settings, language models and datasets, ZFO frequently improves optimization and final performance relative to fixed-step first-order baselines, with the magnitude and preferred local model depending on the objective. Our code is publicly available at: https://github.com/nizswan/Zeroth-First-Order-Framework.

On broad classes of linear systems, the shortest experiments are almost as good as the best possible ones. For $n$ states and $m$ inputs, the shortest input sequences that support robust data-driven stabilization of every controllable plant have $mn+1$ steps with exact states and $m(n+1)$ with noisy states. We show that, when the spectral radius is bounded and the spectrum is well separated near the unit circle, these sequences tolerate a fixed fraction of the error level achievable by any experiment, even one designed with full plant knowledge and allowed to use any finite duration. This constant-factor comparison can fail for slowly actuated systems. For $A=I+hG$ with controllability depth $\nu\ge2$, short experiments lose a factor of order $h^{\nu-1}$, and duration of order $1/h$ is both necessary and sufficient to recover a fixed fraction of the optimal tolerance.

We consider the problem of sampling from Gibbs distributions on matrix spaces whose potential energies are neither convex nor globally gradient-Lipschitz. We introduce a family of non-quadratic kinetic energies that lead to a new underdamped Langevin system with momentum preconditioning, in which the gradient of the kinetic energy acts as a smooth spectral taming of the momentum. We prove that, under these relaxed assumptions on the potential, the resulting dynamics leaves the target Gibbs measure invariant, and we establish exponential convergence to equilibrium in a weighted total variation distance. Finally, we show that the corresponding Euler-Maruyama discretization admits moment bounds that are uniform in time, without any modification of the potential gradient, which ensures the stability of the resulting sampling algorithm.

We study linear programming (LP) representations and strongly polynomial algorithms for robust Markov decision processes (RMDPs) with rational polyhedral state-action rectangular uncertainty in rewards and transitions. By encoding a finite sequence of robust policy-iteration steps, we construct a single LP whose optimal solutions recover the robust optimal value and all optimal stationary randomized policies. At fixed discount, the LP has polynomial dimension and encoding length and can be constructed in strongly polynomial time. We also develop a general complexity analysis of robust policy iteration that combines the cost of minimizing over uncertainty sets with the number of iterations needed to evaluate a policy. For a fixed discount factor, we use this analysis to improve the known complexity bounds for $\ell_1$ and $\ell_\infty$ RMDPs and establish new strongly polynomial bounds for general interval, weighted $\ell_1$, and Wasserstein RMDPs, as well as turn-based stochastic games with these uncertainty sets.

We prove randomized matrix-vector query lower bounds for two normalized matrix-game geometries: a Euclidean unit ball against a probability simplex, with row norms at most one, and two probability simplices, with entries of absolute value at most one. Each query returns $(Ax,A^\top y)$ for arbitrary real vectors. The algorithm must return a feasible pair with full saddle-point gap at most $\varepsilon$, with probability at least $2/3$ on every admissible matrix. For sufficiently small $\varepsilon$, the worst-case query complexities are $\Omega(\varepsilon^{-2/3}/(\log^2(1/\varepsilon)\log\log(1/\varepsilon)))$ for ball-simplex games and $\Omega(\varepsilon^{-2/3}/(\log^{7/3}(1/\varepsilon)\log\log(1/\varepsilon)))$ for simplex-simplex games. The hard instances have dimensions of order $\varepsilon^{-2/3}$ and $\varepsilon^{-2/3}/\log^{1/3}(1/\varepsilon)$, respectively, and the bounds extend to larger dimensions. These lower bounds match the deterministic upper bounds of Karmarkar, O'Carroll, and Sidford up to logarithmic factors. The proof extracts a fresh Gaussian core after adaptive two-sided queries and uses uncertainty in its smallest singular value to establish linear-system solve hardness. Two reductions transfer this hardness to matrix games by converting a small full gap into a small residual, with an additional logarithmic normalization loss only for simplex-simplex games.

We study quadratically regularized optimal transport with quadratic cost in the regime of small regularization $\varepsilon$, where the support of the optimal coupling is sparse. For smooth marginal densities in $\mathbb{R}^d$, we show that the support contains the graph of the Brenier map. After centering by the Brenier image and rescaling by $\ell=\varepsilon^{1/(d+2)}$, the sections of the support have explicit ellipsoidal limits for $\varepsilon\to0$, except at boundary points, where the limits are half-space profiles. The optimal dual potentials admit expansions at order $\ell^2$, uniformly up to the boundary. We identify the leading coefficients, which consist of a common local profile and opposite global corrections determined by a linear Neumann problem. Finally, we analyze two approximations to the Brenier map, namely the gradient of the dual potential and the conditional mean of the coupling. We obtain the sharp $L^p$ rate $\ell^{1+1/p}$ with exact leading constants and further identify the leading interior and boundary biases. Taken together, our results illustrate that quadratic regularization induces an accurate sparse approximation of classical optimal transport.

Given a subvariety of the Grassmannian and a data point, we seek to find a point on the subvariety minimizing the Euclidean distance to the data point. The number of complex critical points of this optimization problem is the Euclidean distance (ED) degree. We show that the ED discriminant of the whole Grassmannian, that is, the set of data points with a number of critical points different than the ED degree, is the discriminant of the characteristic polynomial of the data as a projection matrix. Another closely connected algebraic complexity measure for a subvariety in the Grassmannian is the Grassmann distance (GD) degree, which is the number of complex critical points of the distance optimization problem when the data point itself is in the Grassmannian. We give formulae for ED and GD degrees of geometrically meaningful subvarieties of the Grassmannian, namely, products of smaller Grassmannians, matroid realization varieties, and Schubert varieties.

In [H. Lavenant, J. Luckhardt, G. Mordant, B. Schmitzer, L. Tamanini, The Riemannian geometry of Sinkhorn divergences. Ann. Inst. H. Poincaré Anal. Non Linéaire 43 (2026)] we introduced a Riemannian metric $\mathsf{d}_S$ on the space of probability distributions obtained from entropic optimal transport, specifically from the Sinkhorn divergence $S_\varepsilon$. In the present work we discuss how to approximate and compute $\mathsf{d}_S$. Spatially, we prove Gromov–Hausdorff convergence of the metric and convergence of geodesics for increasingly fine Eulerian discretization of the base space. Temporally, we show $\Gamma$-convergence of the chain discretization $N \sum_{k=0}^{N-1} S_\varepsilon(\mu_k, \mu_{k+1})$ to the energy functional defining $\mathsf{d}_S$. We deduce and implement numerical schemes to compute approximations of $\mathsf{d}_S$.

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.

The design of satellite constellations for Earth observation requires balancing spatial coverage, revisit time, cost, and operational complexity. This paper considers the problem of designing the orbits of a given number of Low Earth Orbit (LEO) or Very Low Earth Orbit (VLEO) satellites to maximize the spatial and temporal resolution achieved over a prescribed set of ground targets. This kind of problem is inherently nonconvex and possibly combinatorial, making its solution computationally demanding for large constellations and target sets. To address this challenge, we propose a two-stage optimization strategy that separates spatial-coverage design from temporal-resolution optimization, thereby reducing the complexity of the overall problem. Two variants of the method are developed. The first employs continuous decision variables during the spatial-optimization stage, whereas the second discretizes these variables and reformulates the problem as a Quadratic Unconstrained Binary Optimization (QUBO) problem. The latter formulation enables the use of efficient classical QUBO solvers and is directly compatible with quantum-annealing hardware. The proposed framework provides a scalable approach to the design of heterogeneous LEO and VLEO Earth-observation constellations and establishes a pathway for exploiting emerging quantum-optimization technologies in satellite mission design. Preliminary simulation results are presented to demonstrate the effectiveness of the strategy.

We study nonconvex–strongly-convex bilevel optimization under a stochastic first-order oracle. We introduce MRT-FD, a single-loop first-order method that simultaneously tracks the upper-level variable, the lower-level solution, and the auxiliary response arising from implicit differentiation of the hyperobjective. MRT-FD performs one update of each variable per iteration and approximates the second-order derivative actions using order-$p$ finite differences. For any fixed finite smoothness order $p\ge1$ in the lower-level variable, MRT-FD finds an $\varepsilon$-stationary point using $\mathcal{O}(\varepsilon^{-4-2/p})$ stochastic gradient queries. We also prove a matching $\Omega(\varepsilon^{-4-2/p})$ oracle lower bound. The lower-bound construction starts from a hard nonconvex minimization chain with a stronger stochastic oracle, and lifts it to a bilevel problem through a sinusoidal coupling with a scalar lower-level variable. Consequently, the dependence on $\varepsilon$ is optimal for every fixed finite $p$, closing the upper–lower complexity gap in this stochastic first-order oracle setting.

Signomials—finite sums of generalized monomials with real exponents—are a natural machine learning model, combining parsimonious representations of power laws, inverse relationships, and multiplicative interactions with universal approximation and interpretable parameters. In fact, 45 of the 100 equations in the AI Feynman benchmark admit signomial representations. Moreover, signomial optimization is widely used in engineering design, communications, economics, and machine learning, but is computationally intractable in general. Convexity, the gold standard for efficient optimization, offers three routes to tractability for signomials: a signomial may be convex in its original variables, become convex after a logarithmic change of variables, or, when positive, become convex after additionally taking the logarithm of its value—the structure underlying disciplined geometric programming. We show that recognizing each form is strongly NP-hard, both on compact domains and globally. Our proofs start from gap-promise variants of polynomial convexity and develop polynomial-to-signomial reductions that preserve curvature gaps under the relevant logarithmic transformations.