Papers

The key-value (KV) cache is the dominant memory bottleneck in long-context language model inference. Existing compression methods apply low-rank factorization or quantization independently, without jointly allocating rank and precision under a shared storage budget. We introduce JoLT, a training-free compressor that treats grouped prefill caches as fourth-order tensors and applies partial Tucker decomposition along the token and feature modes, the two axes that carry low-rank structure, while leaving the head and layer modes intact. A rotated low-bit quantizer captures the truncation residual, and a single Lagrangian dual allocates per-group Tucker ranks and residual bit-widths under a global byte constraint. FlashJoLT replaces the exact token-mode SVD with a randomized approximation that matches JoLT within the free zone at a fraction of the compression cost, and a fused Triton decode kernel evaluates attention directly over the stored factors without materializing dense KV tensors. Across five models from four architecture families, covering multi-head attention, grouped-query attention, and mixture-of-experts architecture, JoLT achieves 2 - 3x compression with less than 0.2% perplexity degradation, without retraining. On RULER at 64K context with LLaMA-3.1-8B, retrieval accuracy remains near-lossless through 3x and declines by only 0.90 and 2.40pp at 4x and 5x, respectively. JoLT demonstrates that tensor-aware low-rank decomposition and quantized residuals, unified under a single storage budget, achieve near-lossless KV-cache compression across diverse model architectures without retraining.

We provide sharp bounds on the dimension of the abnormal set for rank $2$ polarizations on metabelian Lie groups, establishing the Sard property for the end-point map of such groups. The proof is based on a novel approach that makes essential use of tools from tame geometry. We also obtain bounds for the dimension of the Goh-abnormal set for metabelian Lie groups where the codimension of the derived subgroup is at most $2$, with no assumption on the rank of the polarization. We thus infer that these polarized groups, equipped with sub-Riemannian structures, satisfy the minimizing Sard property.

In this paper, inspired by the idea of Adaptive Gradient Descent (AdGD) by Malitsky and Mishchenko, we propose a method that adaptively provides a stepsize interval, which represents the condition for stepsizes ensuring global convergence. By projecting the Barzilai Borwein stepsize onto this interval, we can ensure global convergence and expect further acceleration of convergence. Furthermore, we propose enlarging the stepsize interval obtained by AdGD and show that global convergence can be guaranteed within this enlarged stepsize interval. This extension of the interval enables us to adopt the full Barzilai-Borwein (BB) stepsize more frequently. We report the results of numerical experiments on the stepsize that combines the proposed interval with the BB stepsize.

Shor relaxations are standard tools for nonconvex quadratic optimization, but they may fail Slater's condition. The singularity degree quantifies this degeneracy by the number of facial-reduction steps needed to reach the minimal face. We determine how the constraints defining a binary feasible set affect this quantity. For every nonempty binary set defined by linear equations and inequalities, the Shor relaxation has singularity degree at most one. It is zero precisely when the linear programming (LP) relaxation contains a point in the open cube that strictly satisfies every inequality, and is one otherwise. As an algorithmic consequence of this characterization, Slater's condition can be restored by solving an LP rather than an auxiliary semidefinite program (SDP). In sharp contrast, for every $n\geq1$, we construct a binary feasible set defined by quadratic equalities whose Shor relaxation has a positive semidefinite matrix variable of order $n+1$ and singularity degree $n$. Thus, singularity degree is uniformly bounded in the linearly constrained case but can grow linearly with the matrix order when quadratic defining constraints are used.

This study extends variable-lattice-density optimization to high-Reynolds-number flows for the design of periodically arranged pin-fin heat sinks. Effective permeability and drag coefficient are identified from unit-cell Reynolds-averaged Navier–Stokes analyses of cylindrical pin-fin arrays and incorporated into a reduced model based on the Darcy–Forchheimer law for macroscopic design exploration. To enable stable optimization under high-Reynolds-number conditions, a dual-mesh framework is introduced, in which the flow field and sensitivities are evaluated on a fine mesh, whereas the design variables are updated on a coarse mesh corresponding to the unit-cell arrangement. For the base condition, the $L_2$ norm of the temperature deviation from the area-averaged temperature is decreased from 4.38 K to 1.00 K in the reduced model, and geometry-resolved analysis of the reconstructed design confirms a reduction from 5.98 K to 1.17 K. Additional calculations with higher inlet velocities and modified outlet locations show that the proposed method captures the dominant flow-redistribution trends under different operating and geometric conditions, although the thermal objective became less accurate when local solid-temperature variations became pronounced. These results indicate that the proposed approach is useful as practical design-exploration for identifying candidate pin-fin heat sink configurations under high-Reynolds-number conditions.

This paper investigates a regional optimal control problem for a nonlinear spatio-temporal epidemiological model involving fractional diffusion on a bounded domain. In this work, we consider a general form of disease transmission, and two types of control, vaccination and treatment. We establish the existence and uniqueness of global-in-time solutions to our proposed system. We also prove the existence of an optimal control that minimizes the infected individuals as well as the cost of vaccination and treatment. The necessary conditions for optimality are derived. We show that concentrating vaccination in a selected region of the domain can substantially mitigate the spread of the disease. Numerical simulations are given.

We classify the local maximizers of the entrywise fourth-power objective \[ Q\longmapsto \lVert Q\rVert_4^4 =\sum_{i,j=1}^r q_{ij}^4 \] over the real orthogonal group $\mathcal O(r)$. We prove that the signed permutation matrices are the only local maximizers, and hence the only global maximizers, in every dimension. More strongly, every other stationary point has an explicit rank-two tangent direction with strictly positive second variation. The proof is based on a maximum-entry pivot for the orthostochastic matrix $Q^{\circ2}$: the associated full Riemannian Hessian can be evaluated exactly and is positive at a largest nonunit squared entry. The argument is self-contained and handles zeros, repeated magnitudes, reducible support, and Hadamard-type stationary points.

We study the control and homogenization properties of a coupled hyperbolic system characterised by a rapidly oscillating coefficient present in the principal term. Due to the high-frequency microstructural variations of the density parameter $\varepsilon$, classical uniform controllability results across all initial data fail under boundary controls alone. To overcome this limitation, we divide the spectrum into low and high-frequency components. Utilizing Ingham-Beurling-type inequalities combined with spectral gap analysis, we demonstrate the uniform null controllability of the projected low-frequency and high-frequency states. Furthermore, we analyse the structural limits as $\varepsilon \to 0$ and prove that the associated sequence of boundary controls converge weakly to a null control for the underlying homogenized system. Finally, we show that by using an additional interior feedback control strategy, full space controllability can be achieved for all initial data in the energy space.

Interior-point methods (IPMs) are widely adopted due to their high practical efficiency in solving linear, convex, and nonconvex optimization problems. For convex optimization, this performance is theoretically well-supported: there are strong complexity guarantees for self-concordant barrier setups \cite{nesterov1994interior}, which cover linear and conic optimization. On the other hand, unconstrained convex optimization is well-studied. However, there is limited analysis of constrained convex optimization methods without the self-concordance assumption. We develop and analyze a regularized Newton method with line search applied to the log barrier function in the setting that the objective and constraints are thrice differentiable and have Lipschitz continuous first and second derivatives. Starting from a strictly feasible point, our method finds an $\epsilon$-approximately optimal solution in $\tilde{O}(\epsilon^{-2/3})$ iterations.

Nearest structured approximation and best structured preconditioning solve different matrix optimization problems. We determine their exact relation for Kronecker positive-definite matrices under the affine-invariant Riemannian metric. The Kronecker family is closed and geodesically convex, so every full matrix has a unique affine-invariant projection. Its logarithmic residual satisfies partial-trace normal equations and yields certified point and objective errors for an Armijo projection solver. Our central result shows that this unique projection is also a minimizer of the Hessian-relative condition number if and only if the extreme spectral states admit identical tensor marginals. A computable marginal-mismatch residual either vanishes at a condition-optimal projection or produces a strict descent direction. Two relative spectral levels always force projection optimality; more strongly, every $2\times 2$ Kronecker projection is condition-optimal. An explicit $2\times 3$ construction is therefore a dimension-minimal strict separation. Residual-calibrated bounds further bracket the best attainable Kronecker condition number and the suboptimality of the projection. Supporting results place classical diagonal and block Loewner sandwiches, fixed-basis primal–dual obstructions, and general log-spectral targets in the same certificate language. Given validated numerical enclosures and outward-rounded comparisons, an interval-safe corollary preserves the soundness of the full Kronecker tests. Deterministic small-matrix checks, including a multistart generic log-factor oracle independent of the partial-trace solver, verify the stated identities and bounds.

Safety verification is a fundamental problem in control theory. Barrier certificates (BCs) provide a powerful formal mechanism, yet deriving BCs is computationally intensive. This paper introduces a generative framework that leverages large language models (LLMs) to synthesize BCs through reasoning. Based on the classical Sum-of-Squares (SOS) approach, we train a domain-specific LLM capable of generating high-quality BC candidates for nonlinear systems. Then, the LLM-generated BCs transform the intractable Bilinear Matrix Inequality (BMI) solving problems into convex Linear Matrix Inequality (LMI) feasibility test, significantly improving efficiency while preserving correctness. Experimental results show that our generative method achieves several orders of magnitude speedup over traditional numerical BC approaches and, perhaps surprisingly, surpasses the state-of-the-art dedicated neural BC model. These findings mark a substantive step toward integrating generative AI with formal safety verification for dynamical systems.

We revisit the classical topic in dynamic optimization: sequential vs simultaneous methods for solving DAE-constrained optimization problems, with a particular focus on how graphics processing unit (GPU) computing changes their effectiveness. Sequential methods offer key advantages through adaptive time stepping at the differential-algebraic equation (DAE) solver level, which is especially effective for handling stiff systems. However, long-time-horizon simulations remain a computational bottleneck, as time integration is inherently sequential and limits parallelization within the optimization algorithm. In contrast, simultaneous approaches are well-suited for parallel computing. They address these limitations by exploiting the highly repetitive structure of discretized DAE systems at the function evaluation level and leveraging sparse linear algebra routines that enable elimination tree-level parallelism. Although simultaneous methods typically lack adaptive time stepping, this limitation can often be mitigated by choosing a sufficiently fine initial mesh or iteratively adjusting mesh coarseness in an outer loop. In this work, we revisit the capabilities of the simultaneous approach in a GPU computing environment and assess its performance against a sequential method baseline. We employ a simultaneous approach based on orthogonal collocation within an open-source modeling framework and apply it to parameter estimation benchmarks from systems biology. We evaluate both GPU and CPU solvers on the simultaneous formulation. Our results show that, although less reliable, the simultaneous approach achieves up to 5.4x speedup compared to the sequential baseline among the largest instances where both methods solve successfully. The advantage of the simultaneous method with respect to problem size is more pronounced on GPUs than on CPUs.

Timely distribution of medical items (e.g., whole blood and vaccines) across regional healthcare networks often requires high-volume delivery operations from central facilities (e.g., blood banks) to intermediate regional hospitals, followed by rapid last-mile deliveries to points of injury. Traditional ground-based delivery systems often suffer from limited responsiveness (e.g., traffic congestion) and operational inefficiencies (e.g., blood waste). Leveraging aerial-drone-based delivery systems offers a promising solution for the fast and efficient delivery of time-sensitive medical items across regional healthcare networks. Therefore, we study a hierarchical hub-and-spoke drone-based network for delivering time-sensitive medical items with distinct release and due times to fixed and mobile delivery destinations. We consider a heterogeneous fleet of drones with distinct characteristics (e.g., cost, battery capacity, and speed) and different multi-trip delivery modes. We propose an efficient mixed-integer programming model for location/allocation of mobile delivery destinations, as well as routing and scheduling drones to minimize the total investment and operational costs of the drone delivery network while maintaining the delivery due times. We develop a customized exact solution method integrating problem-specific reformulations and dynamic cutting planes, as well as a fast heuristic algorithm by leveraging a simplified problem variant. Results based on a real-life case study of whole blood delivery data from Pendleton, Oregon, United States, and actual drone flight test data demonstrate that our exact and heuristic solution methods are 31 and 3,375 times faster, respectively, than the Gurobi solver. Results also show that allowing drones to perform multiple trips is 142.8% more cost-efficient than single trips.

This paper investigates the control of discrete-time linear time-invariant (LTI) systems subject to incomplete and corrupted measurements. Specifically, we focus on designing a Linear Quadratic Gaussian (LQG) controller without relying on explicit state estimation. By leveraging minimum variance duality, our approach allows the current control input to be represented as a linear function of available measurements and previously applied inputs, successfully reducing the task to a tractable deterministic optimization problem. We provide theoretical justification for this framework and demonstrate its practical effectiveness through numerical experiments.

A separable bilinear layer relaxation assigns a population cost to each layer and a nonpositive bilinear coupling to each pair of layers, then excludes a profile span whenever the relaxed minimum exceeds a reference value. We ask what a Kuznetsov-Sahinidis-type relaxation can certify at best. Let gamma = g(1) be the singleton-layer cost. For every valid relaxation with gamma > -1/2, every reference bounded per particle, and every underlying system satisfying two elementary low-energy witness conditions, we prove that the certified deficit N - rho(N) is at most of order sqrt(N), with an explicit bound depending only on the reference bound and gamma. Thus no relaxation in this class can certify a larger asymptotic deficit, independently of the detailed layer cost, interaction kernel, or reference. The threshold is sharp in the abstract witness class: at gamma = -1/2 we construct a valid relaxation with a bounded sound reference and deficit N - 1. The geometric Lennard-Jones system behaves differently. A three-atom chain restores the square-root ceiling at the endpoint, while sufficiently negative singleton costs allow valid relaxations with linear deficit. For still lower costs, beyond the stability constant, no sound relaxation excludes any span. These regimes define a certification spectrum for layer relaxations and reveal a reentrant transition between no exclusion, linear exclusion, and the universal square-root ceiling.

This article describes a novel numerical algorithm for tracking the motion of a collection of open curves in two space dimensions. The proposed strategy combines two complementary representations of these curves at each iteration of the evolution process: on the one hand, they are meshed explicitly, as a sub-collection of the entities of a mesh of the total computational domain. Concurrently, using a variant of the Level Set Method, they are captured implicitly as algebraic combinations of the negative, zero and positive subsets of two auxiliary scalar functions, defined on the whole ambient space. This coupling of representations allows to perform accurate geometric or mechanical computations on these curves, while leaving room for large evolution of their shape. After the description of its main numerical ingredients, several applications examples of this methodology are proposed, where it is used to simulate the motion of open physical discontinuities, such as a vortex sheet roll-up, and to optimize the shape of open-ended curves, e.g. with respect to their anisotropic length, with the aim to improve the trajectory of a laser acting on a powder bed in the context of additive manufacturing, or to reconstruct fracture sets.

This article establishes a new numerically viable technique for solving a class of constrained, nonconvex, continuous-time optimal control problems (OCPs) for linear systems that commonly arise in aerial and aerospace applications. The lossless convexification technique is employed to translate the original nonconvex OCP with annular control magnitude constraints into a convex problem, and then by finitely parametrizing the control space with piecewise constant functions, an efficient numerical approach is established that guarantees exact solutions while ensuring the satisfaction of an uncountable family of constraints over a compact time interval. The effectiveness of the approach is demonstrated on a spacecraft landing problem involving three degrees of freedom (DoF), underscoring its potential for real-world aerospace guidance and control tasks.

We study deterministic first-order minimization of a convex function without prior knowledge of the objective's growth, smoothness regime, or associated parameters. We develop anytime, parameter-free bundle-level methods that adapt simultaneously to these unknown properties and attain best-known oracle complexities. For nonsmooth Lipschitz objectives satisfying quadratic growth, the proposed bundle-level W-certificate method (BLW) achieves the optimal complexity without requiring the growth modulus or target accuracy as input. We then introduce an accelerated variant, A-BLW. Without knowing the Hölder smoothness parameters, the quadratic-growth modulus, or the target accuracy, A-BLW attains the optimal rates in the nonsmooth, weakly smooth, and smooth regimes. Central to both methods is an affine W-certificate, a condition based on the descent-slowness of an affine minorant that converts the geometry of a bundle model into an optimality-gap guarantee under quadratic growth. A stopping-time analysis further shows that the same A-BLW algorithm, without modification, achieves the corresponding best-known rates for general convex objectives and for objectives satisfying Hölder growth of order at least two. Numerical experiments illustrate the practical performance of the proposed methods.

Many numerical algorithms for optimal control leverage an elimination of the state via the control-to-state map such as condensed approaches or preconditioned conjugate gradient methods for the optimality system. As such, the norm of the control-to-state map directly enters the convergence estimates for these methods, e.g., via the condition number of the associated linear system. In this work we show that using feedback transformations one may reformulate the optimal control problem to decrease the norm of the (feedbacked) control-to-state map, leading to a drastic improvement of the involved condition numbers. We illustrate the abstract approach for ordinary and partial differential equations such as parabolic, hyperbolic or elliptic equations. For each of these problem classes we provide a constructive method to improve solution operator norms via feedbacks. Further, we showcase the efficacy of the method by means of various numerical examples with elliptic, parabolic and hyperbolic partial differential equations.

This article presents a novel, numerically viable algorithm for solving sparse robust optimal control problems in continuous time. We consider a constrained linear noisy system governed by an ordinary differential equation (ODE), with an $L^1$-type objective function in line with the sparse optimal control literature. The resulting optimal control problem is shown to admit a semi-infinite programming (SIP) formulation. Building upon this insight, we develop a new framework that enables the computation of exact solutions – to our knowledge, the first such achievement in the context of sparse optimal control. We demonstrate that a finite and computationally viable convex optimization problem can be solved to recover, in a lossless manner, both the optimal value and the corresponding optimizers of the original SIP, while also guaranteeing satisfaction of uncountably many constraints. We also show that the parameter-dependent noisy systems and the minimum attention problem fall into our framework and can be solved efficiently via our algorithm. The efficacy of our algorithm is illustrated through a benchmark numerical example.

In multi-objective combinatorial optimization, unsupported non-dominated points typically outnumber supported points and are often significantly more challenging to compute. Recent studies show that extreme supported non-dominated points provide high-quality representations of the non-dominated set for certain binary problems. We demonstrate that this observation does not generalize to capacitated network optimization problems: representation quality decreases with increasing arc capacities, whereas supported non-dominated points consistently provide high-quality representations with respect to several quality indicators. However, supported point sets may still be too large in practical applications, where only a small, fixed number of alternatives is typically desired. Selecting fixed-size representations from the non-dominated set requires its computationally expensive generation and thus diminishes the computational advantages that representations are intended to provide. We therefore suggest the (extreme) supported points as alternative candidate sets in subset selection problems. Our numerical results show that restricting the candidate set to supported non-dominated points yields fixed-size representations of nearly the same quality as those selected from the complete non-dominated set. Overall, supported non-dominated points serve both as high-quality representations and as reasonable candidate sets for subset selection.

For decades, static solution concepts (Nash, Correlated, and Coarse Correlated Equilibria) and the Price of Anarchy (PoA) have formed the bedrock of algorithmic game theory, with no-regret learning proving fast convergence to such game-theoretic equilibria. We show that reducing multi-agent learning to static equilibrium and black-box regret analysis obscures underlying dynamic disequilibrium and game theoretic bounds. First, interior Nash equilibria lack $C^1$ vector field information, meaning agents cannot distinguish aligned from strictly opposing incentives. Inheriting this geometry, the worst-case pure Nash equilibria dictating robust PoA bounds manifest as topologically unstable strict saddles, and in canonical congestion games, as global repellers supported on almost everywhere strictly dominated strategies. Anchoring efficiency guarantees to these unstable states causes algebraic sensitivity; we prove that accommodating all strictly positive affine costs renders the PoA unbounded. Furthermore, projecting learning trajectories onto the discrete simplex of correlated play systematically accommodates non-rationalizable behavior. Evaluating dynamics via Coarse Correlated Equilibria or proximal refinements fails to preclude strictly dominated strategies. Moreover, optimal $O(1/T)$ swap-regret minimization does not preclude macroscopic turbulence, manifesting as chaotic limit sets even in minimal games. Finally, we examine the non-atomic limit of congestion games. Though considered highly stable with tight sub-linear $\Theta(p/\ln p)$ PoA bounds (where $p$ is the polynomial degree), we prove that under discrete-time learning, the unique equilibrium destabilizes into Li-Yorke chaos and global attractors whose time-averaged inefficiency degrades exponentially as $2^p$. These results necessitate re-evaluating worst-case equilibrium frameworks for dynamically grounded metrics.

Prefabricated prefinished volumetric construction moves most building work into module factories, whose production floor operates as a flexible job shop. A major complication is decisive: long post-operation time-lags caused by concrete curing, watertightness ponding tests, and paint drying, during which a module is blocked while its workstation stays free. On benchmark instances grounded in an official national prefabrication guidebook, these lags inflate even the optimal reference makespan by about 67% on average, and ignoring them at decision time, then repairing to feasibility, is worse than every dispatching rule. We adapt a state-of-the-art dual-attention deep reinforcement learning solver through three minimally invasive, individually ablatable extensions: lag-aware dynamics with an admissible reward bound, two anticipatory lag feature channels, and liveness-masked operation- and station-type embeddings. With every extension disabled the implementation reproduces the original solver exactly, so all gains are attributable to the adaptations. We release a public, guidebook-grounded benchmark generator. On held-out instances the learned policy is the strongest solver-free scheduler: it reaches within about 4% of a constraint-programming reference and beats every dispatching rule and a genetic-algorithm metaheuristic, with its advantage widening under capacity contention, and a single size-mixed policy carries this lead across the trained range of factory sizes. It needs no solver, model, or license in the loop and re-plans within seconds of a disruption; where an exact solver can be deployed, that solver remains the quality ceiling, a boundary we map explicitly.

We study a general decision-dependent contextual stochastic program (DD-CSP) in which uncertainty depends on both exogenous contextual information and endogenous decisions. To learn the potentially complex dependence of uncertainty on decisions and contextual information, we employ several nonparametric regression models, including k nearest neighbors (kNN), classification and regression trees (CART), and ReLU neural networks. To account for estimation errors in predicting the uncertainty, we adopt an empirical residuals-based decision-dependent sample average approximation (ER-DD-SAA) framework, which adds empirical residuals to the point predictions from the learned regression models. For each nonparametric regression model, we develop exact mixed-integer programming (MIP) representations that can be seamlessly embedded within the ER-DD-SAA framework. For two-stage ER-DD-SAA problems with kNN, we further propose a tailored decomposition algorithm, named BD-CG, that combines Bender's decomposition with constraint generation. Under suitable assumptions, we prove that the proposed BD-CG converges to a global optimum within a finite number of iterations. From a statistical perspective, we establish the consistency and asymptotic optimality of ER-DD-SAA with all three nonparametric regression models under mild regularity conditions. Numerical experiments on a newsvendor problem with pricing and a two-stage facility location problem demonstrate that the ER-DD-SAA model with nonparametric learning consistently outperforms a parametric benchmark in out-of-sample performance and the proposed reformulations and algorithm substantially improve computational tractability.

This paper investigates the dependence of the control cost for a wave equation with respect to perturbation by a time-independent potential $\lmbd V$ scaled by a large parameter $\lmbd$ on a compact Riemannian manifold. We introduce the geometric control condition~\eqref{GCC+}, a variant of the geometric control condition of Bardos–Lebeau–Rauch–Taylor, tailored to accommodate the influence of the potential $V$. We show that~\eqref{GCC+} is necessary and sufficient for the existence of a uniform observability cost with respect to the large parameter $\lmbd$. We provide geometric examples satisfying~\eqref{GCC+} and estimate the blow-up rate of the observability cost in situations where it fails. The proofs rely on semiclassical and second microlocal defect measures.

Lecture notes for a course on convex optimization taught by Andreas Habring at Graz University of Technology in 2026. The notes cover mathematical preliminaries, fundamental results about existence and of solutions for minimization problems, projected subgradient descent, proximal-gradient methods, heavy ball gradient descent, Nesterov accelerated gradient descent and FISTA, primal-dual methods, and a short intro to optimal transport.

In 1992, Billera and Sturmfels introduced coherent monotone paths on polytopes as part of their description of the fiber polytope construction, and later in 1994 showed with Kapranov that these coherent monotone paths capture the topology of the space of all monotone paths, paths from a minimum to a maximum, in the directed graph of a polytope with orientation induced by a linear function. Those results motivate the following analog of the polynomial Hirsch conjecture: Does there always exist a coherent monotone path of polynomial length on a polytope for any choice of orientation induced by a linear function? We show this is not the case by exhibiting a family of polytopes and corresponding linear functions for which every coherent monotone path is exponentially long. As applications, we strengthen longstanding results pertaining to lower bounds for the shadow simplex method, geometric transversals in discrete geometry, and parametric linear optimization.

We study the policy convergence of unregularized policy mirror descent (PMD) with arbitrary constant step sizes for finite discounted Markov decision processes. We focus on decomposable mirror maps of the form $h(p)=\sum_a \psi(p(a))$, where $\psi$ satisfies standard Legendre-type assumptions. Under these conditions, we prove that the policy sequence generated by PMD converges in the policy domain to a limiting optimal policy, even when the optimal policy set is not a singleton. This result covers a broad class of commonly used mirror maps, including the squared Euclidean mirror map underlying projected Q-ascent, the negative Shannon entropy underlying softmax natural policy gradient, Tsallis entropy, the Hellinger mapping, and the Fermi-Dirac entropy. Although policy convergence has been established previously for specific PMD instances or for regularized variants such as homotopic PMD, to the best of our knowledge, this is the first systematic and unified policy convergence theory for unregularized PMD under general decomposable mirror maps and arbitrary constant step sizes. Our analysis further reveals that the convergence behavior is governed by the differentiability of $\psi$ at $0$ and $1$, leading to different behaviors, including finite-time termination, asymptotic convergence, and an MDP-dependent dichotomy. When $\psi$ is twice continuously differentiable with strictly positive finite curvature, we further establish local policy convergence rates for the asymptotic convergence cases, covering the standard mirror maps mentioned above.

This paper introduces numerical optimizations for maximizing throughput on GPU when solving large batches (10,000 to over 100,000) of sequential quadratic programming (SQP) iterations, where all problems have the same structure. The optimizations are implemented in a toolbox WarpMPC for model-predictive control (MPC) in JAX and Warp. Based on the insight that all MPC problem instances in a batch share the same sparsity in time, cost, and constraints, we propose unrolling sparse linear factorizations and solves, which dominate alternating direction method of multipliers (ADMM) solver runtime. We avoid memory access bottlenecks and wasting computations via optimized memory layout, padding-reducing segmentation of the unrolled factorization, and dependency level scheduled backsolves, additionally accelerating sensitivity computation. We achieve throughputs of 8,000 to 250,000 SQP iterations per second on nonlinear cartpole, quadrotor, and humanoid robot benchmarks, outperforming baselines by 3$\times$ to 25$\times$. We illustrate practical usefulness by synthesizing a dataset and training a neural network approximation of an MPC in under 4 minutes that stabilizes a nano quadrotor in hardware experiments.

Two-stage stochastic mixed-integer programs are important tools for decision-making under uncertainty. Representing the uncertainty with many scenarios, however, can make them challenging to solve. Scenario reduction addresses this by finding a distribution supported on fewer scenarios that still yields similar optimal first-stage decisions. In this paper, we revisit the classical scenario reduction theory based on distances between probability distributions and the optimal mass transportation problem. The transportation problem's cost function captures scenario similarity and is central to the effectiveness of scenario reduction. We then review and compare various transportation cost functions from the literature and propose a new one. Using the Forward Selection Algorithm, we prove that our proposed cost function selects the best possible scenario from a given sample on the first draw with respect to the relative approximation error. To reduce the computational cost of evaluating this cost function, we further propose a hybrid algorithm with a scenario pre-selection phase. We assess solution quality and computational complexity on the two-stage stochastic unit commitment problem for small 24-bus and large 300-bus case studies. With only around five scenarios, the proposed cost function approximates the full-distribution optimum to within roughly 2.1% and 0.4% error for the small and large cases, respectively. In contrast, prevalent cost functions often need 25 scenarios or more to achieve that solution quality. The hybrid algorithm achieves similar solution quality while reducing wall-clock time by a factor of 18 and work (per Gurobi solver) by a factor of 66 on the large case study.