Non-Negative Conjugate Gradients

The conjugate gradient method (CG) solves a symmetric positive definite (SPD) system $Ax=b$, but does not respect $x\geq0$. We develop non-negative conjugate gradients, a solver for the bound-constrained quadratic $\min_x \tfrac12 x^\top A x - b^\top x$ subject to $x\geq0$, $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\'udice 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\mapsto M^\top(Mv)$ needs two products with $M$, no $O(n^2)$ storage, and retains CG's $O(\sqrt\kappa)$ rate, $\kappa=\kappa(A)$. Regularisation or preconditioning that compresses the spectrum lowers the iteration count; a ridge/Tikhonov split $A\mapsto(1-\alpha)A+\alpha R^\top R$ also secures the needed $P$-matrix property. A loop-free projected-gradient reference method, analysed alongside, converges at the slower $O(\kappa)$ 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.