Papers

Nested clustered optimization allocates within each cluster from the cluster's own covariance block and then across the resulting cluster portfolios. Block inversion says the unconstrained minimum-variance portfolio has the same two-tier shape, with each block replaced by its Schur complement against every other asset. Conditioning instead on one knot from each other cluster truncates the conditioning set in the manner of a Vecchia approximation, and we give the rank-one model of cross-cluster dependence under which it is exact. Damping the complement by $\gamma\in[0,1]$ then gives a bridge with nested clustered optimization at $\gamma=0$ and the global optimum at $\gamma=1$, with no linear solve larger than a cluster or the number of clusters. Under estimation error the optimal $\gamma$ can be strictly interior and full coupling can remain optimal, and we give the local theorem at the minimum-variance end with exact examples of both, including a symmetric family in which the optimum is a closed form.

We study online control of a scalar storage system with adversarial nonnegative arrivals, known retention coefficient, and convex costs depending on both state and action. Each action must respect current resource availability and is chosen before the current arrival and cost function are revealed. For the existing simplex disturbance-action policy class, we give an exact reparameterization by cumulative allocation fractions and a decay-weighted projected subgradient update. The resulting regret bound is independent of policy memory length. For fixed retention coefficient and cost constants, the controller achieves $O(\sqrt T)$ regret against the best fixed infinite-memory policy in this class, using $O(\log T)$ memory and arithmetic operations per round and one cost-subgradient query. A storage-specific block construction gives a matching lower bound against every causal feasible controller, including randomized controllers. Writing $\tau=(1-\alpha)^{-1}$, the minimax expected regret is $\Theta(\sqrt T\min\{T,\tau\}^{3/2})$ for every finite-memory simplex policy class and its infinite-memory extension, when $\alpha\in[1/2,1)$, $T\ge4$, and the positive cost constants are fixed. This identifies the joint horizon and retention-time dependence for these policy benchmarks.

Octopus crawling motivates soft robots that exploit redundancy, yet discovering and organizing diverse coordination modes for adaptation remains challenging. To address this, we introduce a Diffusion-based Uncertainty-aware Optimization (DUO) algorithm that learns demonstration-free crawling controllers for a simulated, muscle-actuated CyberOctopus. This work represents the first application of diffusion-based control to soft multi-arm robots in contact-rich simulations. By embedding a variety of locomotion behaviors within a shared control distribution, this approach enables the simulated octopus to navigate dynamic physical constraints, demonstrating that learned coordination diversity inherently facilitates robust adaptation. The main contributions include: (i) a symmetry-structured policy representation that folds radially equivalent controllers into a canonical directional sector, (ii) an online black-box optimization strategy, the DUO algorithm, that discovers and retains diverse coordination modes, and (iii) a control editing technique that adapts existing controllers to novel actuator constraints without retraining. These results show how learned coordination diversity makes motor abundance a practical resource for adaptation in soft multi-arm robots.

GPU implementations of the Muon optimizer and of first-order semidefinite programming (SDP) solvers share one computational pattern: a matrix factorization is replaced by a fixed-depth polynomial filter applied after normalization. When a few dominant spectral components carry most of the input's scale, normalization pushes the remaining spectrum toward zero, where the filter is least accurate. We propose spectral deflation: estimate the dominant components, remove them, filter the normalized residual with the unchanged filter, and restore them. Deflation preserves the target matrix function, and we prove that it strictly reduces the finite-step error of the classical Newton–Schulz family. Implemented with batched randomized SVD for Muon and warm-started subspace tracking for ADMM, deflation consistently improves GPT-2 pretraining over the corresponding Muon baselines with both the Newton–Schulz and Polar Express mappings, also in wall-clock time, and lowers the KKT residuals of factorization-free ADMM on large-scale SDPs within a similar projection time.

We present tight lower bounds for the size of any linear matrix inequality representation of the Euclidean ball, and on the semidefinite extension complexity of the ball. Specifically, we show that any linear matrix inequality representing the Euclidean ball must involve matrices of size at least $n$, and any spectrahedral extension of the Euclidean ball must involve matrices of size at least $\lceil 2 \sqrt{n-1}\rceil$. These match the upper bounds given by existing explicit constructions exactly. Our proofs rely on elementary facts about the dimension of the set of faces of spectrahedra, and indeed extend to any convex set whose extreme points form a semialgebraic set.

The permanent and mixed discriminant of positive matrices are classic problems for which we do not expect an efficient algorithm for exact computation. Thus much work has been done to understand how well we can bound and approximately compute these quantities. One line of research in this area begins with the results of the first author, where van der Waerden lower bounds of $\frac{n!}{n^n}$ are proven for doubly stochastic inputs for both problems, using a simple proof via stable polynomials. Along with the bound itself, the same techniques are used to show that the permanent and mixed discriminant are uniquely minimized at a certain natural symmetric input. In this paper, we generalize those results in two ways. First, we extend the unique minimization results beyond doubly stochastic inputs to other marginals which are near doubly stochastic. This yields the first such unique minimization results for the mixed discriminant beyond the doubly stochastic case. We also discuss why one cannot hope similar results to hold in general for all marginals. Second, we extend the unique minimization result for real stable polynomials to strongly log-concave (aka Lorentzian) polynomials in the doubly stochastic case. This captures an analogous previous result on unique minimization for the mixed volume. Finally, we discuss various open problems related to these results.

The iterative Linear Quadratic Regulator (iLQR) is a widely used algorithm for nonlinear trajectory optimization. At each iteration, it solves a local linear-quadratic approximation of the problem via dynamic programming, propagating a quadratic cost-to-go function. If the Hessian of the cost-to-go approximation is positive-semidefinite, one can derive a square root formulation of iLQR that propagates its Cholesky factor instead. This offers significant numerical advantages - much as square root Kalman filters improve upon their conventional counterparts - particularly when iLQR is used within an augmented Lagrangian framework for handling constraints, where large penalties degrade conditioning. Previous square root formulations of iLQR and related algorithms exist, but they are either numerically suboptimal, algorithmically complex, or both. In this paper, we show that the key to an effective square root formulation lies in the Gauss-Newton (weighted least-squares) structure of the cost function: this yields a positive semidefiniteness property that extends beyond the Hessian to the full augmented cost-to-go matrix, and enables a backward pass of remarkable simplicity in which each step reduces to a single QR-decomposition, from which the feedback gain and propagated Cholesky factor are extracted directly.

We report improved global-minimum configurations for the classical Thomson problem of $N=60$, $61$, $92$, and $99$ repulsive Coulomb charges confined to a disk. By combining the quenched molecular dynamics (QMD) method with the fixed-border heuristic introduced by Amore and Zarate, we systematically obtain configurations with energies $E_{\mathrm{QMD}}(60)=2159.3584240930$, $E_{\mathrm{QMD}}(61)=2237.19264190$, $E_{\mathrm{QMD}}(92)=5358.35353314$, and $E_{\mathrm{QMD}}(99)=6254.83029083$, which improve upon the previously best-known values. For $N=60$, the Voronoi diagram of our configuration differs from the one reported earlier. The symmetries for $N=61$ $(C_{2})$ and $N=99$ $(D_{1})$ are confirmed and are consistent with the known symmetry patterns for these configurations. For $N=92$, whereas the previously reported Voronoi diagram has lower symmetry, our solution exhibits a clear $D_{1}$ (axial) symmetry of defects. The results are highly reproducible across multiple independent runs, providing strong evidence that these configurations are robust global-minimum candidates.

We investigate the sensitivity of solutions to multi-dimensional scalar balance laws with respect to perturbations of the initial data analytically and numerically. Reliable first-order sensitivity information is essential in gradient-based optimization and inverse problems constrained by hyperbolic balance laws. In hyperbolic problems, such perturbations affect both the smooth components of the solution and the locations of shocks. Consequently, classical difference quotients of the solution operator generally fail to converge in $L^1$, even in one space dimension. Although generalized tangent-vector techniques have been developed for one-dimensional problems, extending these ideas to multiple space dimensions is considerably more challenging because it requires a geometric description of hypersurfaces of discontinuity. We represent perturbations of hypersurfaces of discontinuity by normal displacements and account for the induced variation of the normal direction. This representation yields evolution equations for the components of a generalized tangent vector. Based on this calculus, we develop a numerical method for computing first-order variations with respect to the initial data. Numerical experiments in two space dimensions confirm the expected first-order accuracy of the resulting approximation.

We study efficient algorithms for realizing the first-order oracle complexity of optimization of $G$-Lipschitz convex functions with respect to the $\ell_{q}$-norm over an $\ell_{p}$-ball of radius $R$, where $1\leq p,q\leq \infty$. For $p<q$, we obtain error $\widetilde{O}_{p,q}(GR/T^{1/p-(1/q-1/2)_{+}})$ after $T$ oracle queries, efficiently realizing the nearly optimal rates of (MBG+26), thereby resolving the nonsmooth end of the COLT 2015 open problem (Guz15b). In particular, the rate is $\widetilde{O}(GR/T)$ for Euclidean Lipschitzness over an $\ell_1$-ball of radius $R$ ($p=1,q=2$). Our solution consists of reducing convex Lipschitz optimization to the chasing nested convex sets problem in sublevel sets of an evolving bundle (LNN95; BBE+20): at each query we either find a point with low function value or we produce a deep cut in the current sublevel of the bundle, that we chase. The dichotomy between stability of selectors and forced movement by deep cuts bounds the number of iterations of the algorithm near optimally. For nested subsets of $R B_{p}^{d}$, we introduce a novel notion of stable center whose movement is bounded by $\widetilde{O}_{p,q}(RT^{1-1/p+(1/q-1/2)_{+}})$ in the $\ell_{q}$-norm after $T$ steps, which we show is nearly optimal in high dimensions. A Monte Carlo average of the proposed selector achieves near-optimal rates with high probability and can be implemented in polynomial time for our optimization algorithm in the real-arithmetic model.

We study first-order black-box convex optimization over an $\ell_p$-ball for objectives Lipschitz in the $\ell_q$-norm, solving in the affirmative the nonsmooth version of the COLT open question (Guz15b) on whether the geometry of a smaller feasible set ($p < q$) can improve convergence rates in convex optimization, and matching prior lower bounds up to logarithmic factors. Our rates include \(\widetilde O(1/T)\) for convex Euclidean-Lipschitz optimization over the $\ell_1$-ball, improving on the $O(1/\sqrt{T})$ classical rate under general assumptions. The key technical device is a new online learning game, where the comparator is evaluated using the maximum of affine losses observed so far. We bound the value of this game above and below in terms of a combinatorial online learning quantity: the sequential fat-shattering dimension, which we characterize for the $\ell_p / \ell_q$ case. Our results generally apply when the feasible set $X$ and the set of possible subgradients $H$ are convex, centrally symmetric, and admit a type of minmax theorem, advancing on a fundamental question by Sridharan [Sri12, Section 10.1.2, Q3]. As a geometric consequence of our analysis, of independent interest, we obtain estimates for the expected distance of a convex hull of samples to their mean in several Banach geometries, a version of the celebrated Wendel's theorem (Wen62), but quantitative and for bounded general distributions as opposed to centrally symmetric ones.

We show that unless P = NP, there cannot be a polynomial-time (or even pseudo-polynomial-time) algorithm for output feedback stabilization of a linear dynamical system with a linear controller. This settles one of the best-known open problems in control theory. The result holds in both continuous and discrete time. We also present a family of stabilizable linear dynamical systems for which no polynomial-time algorithm can write down a stabilizing controller in its standard representation.

We study the convergence rate of the moment-SOS (sum-of-squares) hierarchy for polynomial optimization problems (POPs) on a bounded subset of the real line described by arbitrary polynomial inequalities. We prove that, for every fixed univariate POP, the relaxation error is bounded by $O(1/r^2)$, where $r$ is the relaxation order. In particular, boundary degeneracies in the polynomial description of the feasible set affect the constant but not the convergence exponent. The proof combines the structure of finitely generated univariate quadratic modules with a Chebyshev polynomial construction that approximately recovers the natural generators of the feasible set while controlling the degree of the certificate. We also give an elementary degree-four example for which the relaxation error is exactly $1/(2r(r-1))$, showing that the quadratic rate is optimal. Equivalent reformulations connect this example to a cubic univariate problem and to a bivariate POP whose feasible set has a cusp singularity at the minimizer.

Federated learning trains a shared model from private client data. In practice, data-generating distributions may differ, and the true mixture across clients is often unknown, making the underlying group distribution difficult to specify. Existing approaches address cross-client mixture uncertainty by optimizing against the worst-case mixture, yet assume accurate client-wise distribution estimates. However, these estimates can be unreliable when based on finite samples. To handle both cross-client mixture uncertainty and within-client distributional ambiguity, we construct a global ambiguity set as the union of admissible mixtures of local ambiguity sets. The construction allows client-specific ambiguity radii and admits a client-wise separable reformulation. Leveraging this structure, we establish a high-probability out-of-sample performance guarantee. We further develop a federated algorithm for a penalty-based reformulation and prove its convergence under milder regularity conditions. Simulations validate the algorithm's effectiveness.

Two-timescale stochastic approximation (TTSA) is a fundamental tool for analyzing coupled iterative algorithms in reinforcement learning, optimization, and stochastic control. However, finite-time guarantees for nonlinear two-timescale schemes remain difficult to obtain, especially under constant step-sizes. In this paper, we study nonlinear TTSA with step-sizes $\alpha\gg\beta$. Under standard stability, regularity, and Markovian noise assumptions, we upper bound the mean-squared error and the bias of both iterates around their limiting equilibria. Our bounds scale as $O(\alpha+\beta^2/\alpha^2)$, which we prove to be tight when $\beta\le\alpha^{3/2}$. The analysis separates the contributions of initial conditions, fast-timescale tracking error, Markovian dependence, and timescale coupling, thereby clarifying the origin of the $\beta^2/\alpha^2$ term. Our results reveal qualitative differences from the linear TTSA setting previously studied, showing that nonlinear dynamics introduce additional finite-time effects that are absent in the linear case.

We establish sharp convergence rates for quadratically regularized optimal transport with quadratic cost in the regime of small regularization $\varepsilon$. In particular, we quantify the sparsity of the support of the regularized optimizer. For smooth marginal densities in $\mathbb{R}^d$, this support lies within a distance of order $\varepsilon^{1/(d+2)}$ from the Brenier graph, and every section of the support is sandwiched between balls with radii of that order. The geometry of the support is closely linked to the dual potentials. We show that the potentials satisfy uniform two-sided Hessian bounds and converge uniformly at rate $\varepsilon^{2/(d+2)}$, while their gradients converge at rate $\varepsilon^{1/(d+2)}$. The rates are sharp, and we further identify the regularity threshold where convergence breaks down.

Mathematical Programs with Blocks of Vanishing Constraints (MP-BVCs) are a challenging class of optimization problems due to the loss of standard constraint qualifications. Based on existing theory for classical MPVCs, i.e., MPBVC with block size one, we extend the constraint qualification theory for MPBVCs and show that the known LICQ-type condition for MPBVCs implies the Guignard Constraint Qualification (GCQ). A novel piecewise gradient flow approach converges under GCQ. It directly computes strongly stationary points and is an extension of the piecewise gradient flow approach for classic MPVC problems. Numerical experiments on several truss topology design instances with two and even three applied load forces demonstrate the robustness and effectiveness of the method.

We study smooth strongly convex–strongly concave minimax optimization in the deterministic unconstrained setting, without assuming a bilinear or separable structure. Although existing multi-loop methods attain near-optimal condition-number dependence, standard single-loop methods generally exhibit a substantial complexity gap. To close this gap, we propose the Single-Loop Predictor–Corrector Extragradient Method with Damped Momentum (PCE-DM), which combines an extragradient prediction–correction scheme with a novel auxiliary feedback recursion for the weaker-curvature variable. PCE-DM uses fixed parameters and two new full-gradient evaluations per iteration after one initialization query, while requiring no inner solves, accuracy schedules, or staged restarts. We develop a Lyapunov analysis that controls the predictor–corrector mismatch through corrected-gradient increments and establish last-iterate linear convergence. Specifically, PCE-DM computes an $\varepsilon$-accurate relative solution, measured by the squared Euclidean distance to the saddle point, within $\mathcal{O}\!\left(\sqrt{\kappa_x\kappa_y} \log(2\kappa_x\kappa_y/\varepsilon)\right)$ full-gradient queries. This result closes the condition-number complexity gap between standard single-loop methods and near-optimal multi-loop methods, matching the known lower-bound order up to logarithmic factors while retaining fixed, explicit single-loop updates. Numerical experiments on regularized linear regression and AUC maximization demonstrate the computational efficiency of PCE-DM.

We study continuous-time distributed optimization on the three-dimensional rotation group over undirected tree graphs, where agents with heterogeneous local costs seek a common attitude minimizing the geodesically strongly convex sum of their costs on a geodesically convex operating ball. We propose a nonsmooth protocol with fixed gains that combines local Riemannian gradient descent with signum-based consensus feedback computed from the Lie-group logarithms of relative rotations between neighboring agents. Since geodesic strong convexity alone does not guarantee invariance, as demonstrated by a certified counterexample, we impose an inward-pointing boundary condition and prove that every Filippov solution with all agents initialized in the ball remains there. We construct a nonsmooth intrinsic disagreement Lyapunov function and derive a uniform dissipation bound through a cluster-contraction argument exploiting the tree structure. With bounded local gradients, this bound gives an explicit gain-ratio condition under which we establish exact finite-time consensus for every such solution and obtain a settling-time estimate requiring no structural graph constant beyond the number of agents. After consensus, we identify the sliding dynamics as a scaled Riemannian gradient flow of the summed objective. When the minimizer lies in the interior of the ball, the dynamics converge to it exponentially in geodesic distance. Simulations illustrate our theoretical results.

Gromov–Wasserstein optimal transport (GW-OT) extends classical optimal transport by introducing structure-aware transport cost. This is particularly relevant for spatial transcriptomics, where dynamical reconstruction should preserve tissue structure in addition to matching expression patterns. While static formulations have been widely used for such structure-aware alignment, a general dynamic formulation for reconstructing continuous trajectories is still missing. We introduce Travelling Pair Dynamical Alignment and Trajectory Estimation (TP-DATE), a theoretical and computational framework to generalize GW-OT dynamically in a simulation-free manner. We formulate a broad class of static and dynamic Quadratic-form OT (QOT) through path actions and prove the static dynamic equivalence. We further develop travelling-pair flow matching, which allows interacting conditional paths and marginalizes their interactions into a single vector field. On synthetic and real spatial transcriptomics data, TP-DATE better preserves spatial structure and improves continuous 3D dynamics reconstruction.

We introduce a new model for a regulated day-ahead type auction-based electrical market in which the system operator can measure and limit the producer's emissions when choosing its optimal dispatch. We prove properties of this market model that describe it as a generalization to other previous works in the electrical market literature. Furthermore, we use standard sensitivity analysis tools to measure the impact of these pollution variables on market equilibria and present numerical examples of this effect.

Teams of mobile robots rely on continuous communication with their neighbors for coordination, yet most distributed model predictive control (DMPC) schemes assume the communication network stays connected rather than actively enforcing it. Adding such a guarantee is hard since the usual mathematical condition for connectivity is nonconvex and links every agent to every other, which is incompatible with a scalable distributed real-time controller. We propose a DMPC framework in which each agent is assigned a connectivity contract: a local region prescribing where its predicted positions may lie over the prediction horizon. The contracts are designed so that, as long as every agent stays within its own contract, the team is guaranteed to remain connected. Given the maintained contract graph, an agent builds its contract from a single exchange with its immediate neighbors, after which every agent solves its own optimization problem independently. We prove that the resulting closed-loop system maintains connectivity, avoids collisions, and respects local state and input constraints. Simulation and hardware experiments on miniature autonomous car-like robots demonstrate the approach.

We study Percolation Markov Decision Processes (PMDPs), in which the decision maker repeatedly moves a token through the d-dimensional integer lattice with deterministic transitions and random payoffs assigned to the edges. Payoffs are revealed before the beginning of the decision process and the decision maker aims to maximize the average accumulated payoff over a fixed horizon. We establish the existence of the uniform value (as the horizon tends to infinity) and 0-optimal strategies. In the particular case of Bernoulli payoffs, we establish continuity results for the uniform value. We also present several one-dimensional examples to illustrate that optimal strategies may be very complex. Then, we introduce a dimensional lifting property that allows PMDPs to be approximated by PMDPs of lower dimension and could be used to approximate critical probability thresholds in percolation theory.

In this paper, we reinterpret pairwise comparison theory through the lens of optimal transport. The main contribution is the connection established between pairwise comparison matrices, doubly stochastic scaling, and entropy-regularized transport problems. Using Sinkhorn scaling, we show how consistency, priority derivation, and inconsistency assessment can be studied within a common optimal-transport framework. As methodological consequences of this connection, we obtain a new characterization of consistency, a Sinkhorn-based prioritization method, and an entropy-based inconsistency index. We study these proposals both analytically and numerically and compare them with established methods from the literature

The research studies the extended mean-field anticipated forward-backward stochastic delayed differential equations with common noise (extended MF-AFBSDDEswCN) by expanding the domination-monotonicity conditions. We generalize the conditions from a linear setting to a nonlinear one by including nonlinear adjoint functions which are key to guaranteeing the well-posedness of extended MF-AFBSDDEswCN. Utilizing this broader well-posedness framework for extended MF-AFBSDDEswCN, combined with other refined analytical tools, we examine two classes of stochastic optimal control problems. These include a linear-convex problem and a linear-quadratic problem with input constraints that can be both time-dependent and random. Moreover, the control elements regarding the initial values are path-dependent, which will cause some essential difficulties. For each problem, we establish the existence and uniqueness of optimal controls and provide their explicit, closed-form representations.

We study the convergence of the vanilla stochastic policy gradient method applied to the linear quadratic regulator (LQR) problem. The method is cheap in the following sense: (1) at each iteration only $\tilde{O}(1)$ interactions with the environment are needed, therefore allowing frequent policy improvement steps, and (2) to ensure stability throughout and convergence to an $\epsilon$-optimal policy with probability $1-\delta$, only $O(\mathtt{Polylog}(1/\delta)/\epsilon)$ interactions are needed. To the best of our knowledge, this appears to be the first time that a stochastic model-free policy optimization method for LQR converges with high probability with $\tilde{O}(1)$ per-iteration computation and polylogarithmic dependence on the confidence level. The convergence analysis presented here is agnostic to LQR specifics and hence could be potentially generalized to a broader class of problems.

We study a Stackelberg model where two players compete in a sequential contest, and the player that successfully invests the most effort wins. Our formulation considers a probability that each player fails to invest their chosen level of effort, as captured by a failure function that increases with the amount of effort. The players each have their own cost of failure, as well as a maximum level of effort (capturing different abilities of the players). We prove three key results. First, we provide the optimal strategy for Player 1 (the follower), to maximize their utility given an observed effort exerted by Player 2 (the leader). Next, we characterize the optimal strategy to maximize the expected utility of Player 2, anticipating the optimal response by Player 1. Finally, we analyze the expected utilities of each player when both play their optimal strategy, and provide key insights into how heterogeneous player abilities and failure costs affect the outcomes of the game.

Epidemic mitigation comes with a price in practice, and when budget is limited, its effective allocation becomes an important question on public policy. In this paper, we study epidemic control through a Stackelberg mean-field game in which a principal, representing the government, chooses social distancing guidelines and vaccination levels for a large population of rational minor agents subject to a budget-spending process. The principal anticipates the minor agents mean-field Nash equilibrium (MFNE) and selects intervention policies to minimize her own objective, inducing a bi-level optimal control problem. We define the Stackelberg mean-field game equilibrium (SMFE) and give its necessary condition by a forward-backward ordinary differential equation (FBODE) system, and provide existence and uniqueness results. We further propose an iterative fixed-point algorithm to solve the FBODE numerically and present computational observations that illustrate the resulting equilibrium behavior.

The Kemeny rule aggregates multiple strict rankings into a single strict ranking that minimizes the sum of its distances from the input rankings. The resulting optimization problem, called the Kemeny problem (KP), is a special case of the linear ordering problem (LOP). The Kemeny rule satisfies several desirable properties in social choice theory, including the extended Condorcet criterion (XCC). Ando et al. strengthened this result by introducing the strong Condorcet criterion (SCC) and showing that it holds for every optimal solution to an arbitrary LOP instance. Yoo and Escobedo extended the Kemeny rule to rankings with ties and showed that the resulting rule satisfies the non-strict extended Condorcet criterion (NXCC). This criterion gives a condition under which one candidate must be ranked strictly above another in every optimal solution. In this paper, we introduce the non-strict strong Condorcet criterion (NSCC), a counterpart of the SCC for rankings with ties, and show that it holds for every optimal solution to an arbitrary instance of the linear ordering problem with ties (LOPT). We also establish a complementary structural property that gives conditions under which two candidates must be tied in every optimal solution to an arbitrary LOPT instance.

Large-scale nonconvex quadratic programming remains challenging because the sparse matrix factorizations used by many existing solvers can incur substantial computational and memory costs and are difficult to parallelize efficiently on GPUs. We present PDNQP, a factorization-free first-order solver for finding stationary points of large-scale nonconvex quadratic programs. The method combines a proximal augmented-Lagrangian framework with restarted accelerated primal–dual hybrid gradient to solve the resulting strongly convex QP subproblems. Its key ingredients are a residual-form reformulation that avoids penalty-dependent normal matrices and preserves a fixed sparse constraint operator across outer iterations, and an adaptive strategy that coordinates inner-solve accuracy with outer progress. The resulting computations rely on sparse matrix–vector products, projections, and vector operations suited to GPU execution. On $107$ nonconvex CUTEst instances, PDNQP achieves the highest success rate among the tested solvers at both accuracy levels. With absolute and relative tolerances set to $10^{-4}$, it achieves speedups of $2.6\times$–$6.0\times$ in shifted geometric mean runtime; at $10^{-6}$, its aggregate runtime remains comparable to that of the fastest tested solvers. On $30$ large-scale synthetic nonconvex QPs with approximately one million variables each, PDNQP returns points satisfying the common external termination criteria on $27$ instances, while none of the other tested solvers succeeds under the same experimental protocol.