2026/07/24 by Thomas Schmelzer, Martin Stoll
Computer Science · Mathematics · #math.OC #cs.NA #math.NA
The conjugate gradient method (CG) solves a symmetric positive definite (SPD) system Ax=b, but does not respect x≥0. We develop non-negative conjugate gradients, a solver for the bound-constrained quadratic minx \tfrac12 x^\top A x - b^\top x subject to x≥0, Bx=c, A\succ0, with Bx=c present only in the equality-augmented variant. The solver wraps CG in a primal-dual active-set loop: each step solves the unconstrained system over the free variables, drops variables that turn negative, re-admits any bound variable with negative reduced gradient, and restarts CG. These toggles are the principal pivots of the linear complementarity problem (A,-b), well defined since A is a P-matrix. Termination uses the block-principal-pivoting construction of Júdice and Pires, with a least-index Bland fallback, reaching the unique global minimiser in finitely many outer steps. Our contribution carries this guarantee to the inexact, matrix-free regime the Krylov inner solve requires: once CG's residual is tied to the decision margin, the approximate loop makes the same primal and dual sign decisions as the exact one, and a direct factorisation may replace CG unchanged. Each inner solve is matrix-free: for A=M^\top M, the operator v↦ M^\top(Mv) needs two products with M, no O(n2) storage, and retains CG's O(√κ) rate, κ=κ(A). Regularisation or preconditioning that compresses the spectrum lowers the iteration count; a ridge/Tikhonov split A↦(1-α)A+αR^\top R also secures the needed P-matrix property. A loop-free projected-gradient reference method, analysed alongside, converges at the slower O(κ) rate. On synthetic SPD problems the loop reaches the planted optimum in at most 6 outer steps, and with a direct free-set solve outpaces Lawson--Hanson and an interior-point solver by an order of magnitude on dense instances.