Finite-budget primal-dual Newton methods for LP, convex QP, and symmetric cones

This procedure has two admissible modes. All input values are finite stored floating-point numbers. Audits interpret them as exact dyadic rationals, and rounded operations obey Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add.

LP/QP mode. Let $n\ge 1$, $p\ge 0$, $Q\in {\mathbb{R}}^{n\times n}$ be symmetric positive semidefinite, and $A\in {\mathbb{R}}^{p\times n}$ have full row rank. The supplied data are $A,b,c,Q$ for $\min \frac{1}{2}{x}^{T}Qx+c^{T}x$ subject to $Ax=b$, $x\ge 0$; LP is $Q=0$. The iterate inputs are $x>0$, $s>0$, and free $y$. Also supply $0<\sigma <1$, $0<\eta <1$, finite tolerances ${ϵ}_p,{ϵ}_d,{ϵ}_g\ge 0$, and an integer update cap $N\ge 0$.

Symmetric-cone mode. Let $K$ be a nonempty finite product of nonnegative orthant, Lorentz, or positive-semidefinite cones, let $A$ have full row rank, and use the primal-dual pair of Conic, LP, and convex QP primal-dual models. Supply $A,b,c,G,h$, a point $x$ with $s=h-Gx\in \mathrm{int}K$, a dual point $z\in \mathrm{int}{K}^{\ast }$, and free $y$. Equality feasibility $Ax=b$ is checked as a residual, not required of the initial point. Also supply $0<\sigma <1$, $0<\eta <1$, finite tolerances ${ϵ}_p,{ϵ}_d,{ϵ}_g\ge 0$, and integer cap $N\ge 0$. The positive barrier parameter $\nu$ is the sum of the cone-factor parameters.

The input audit checks all dimensions, exact symmetry of $Q$, and full row rank of $A$ by rational elimination. It checks $Q⪰0$ by finite exact symmetric Schur elimination: reject a negative diagonal; if a positive diagonal exists, permute it to the first position and recurse on its Schur complement; if every diagonal is zero, accept only when every off-diagonal entry is zero. This test is equivalent by completing the square at each positive pivot; with zero diagonal, a nonzero off-diagonal makes one of the vectors $e_i+e_j$ or $e_i-e_j$ have negative quadratic form. In symmetric-cone mode it also checks that each listed factor has a positive finite dimension, that the data match the declared product spaces, and that $A$ has full row rank. All these checks use exact rational arithmetic on the stored dyadic inputs.

In LP/QP mode form $r_d=Qx+c-A^{T}y-s,r_p=Ax-b,\mu =x^{T}s\mathrm{/}n,r_c=Xs-\sigma \mu \mathbf{1.}$ In symmetric-cone mode form $r_d=A^{T}y-G^{\ast }z-c,r_p=Ax-b,\mu =⟨s,z⟩\mathrm{/}\nu ,r_c=s\circ z-\sigma \mu e\mathrm{.}$ The cone $r_c$ is the stated diagnostic Jordan residual; for Lorentz factors its target normalization differs from that of the selected log-homogeneous barrier. In either mode, ${ϵ}_p$ bounds $\mathrm{\parallel }{r}_p{\mathrm{\parallel }}_2$, ${ϵ}_d$ bounds both $\mathrm{\parallel }{r}_d{\mathrm{\parallel }}_2$ and $\mathrm{\parallel }{r}_c\mathrm{\parallel }$, and ${ϵ}_g$ bounds the surrogate complementarity gap $n\mu =x^{T}s$ in LP/QP mode or $\nu \mu =⟨s,z⟩$ in symmetric-cone mode. On product cones use Euclidean norms on orthant and Lorentz factors and Frobenius norms on matrix factors. Norm tests compare exact squared values with squared tolerances. The output is a record containing status, the last accepted iterate when one exists, update count, audited residuals, and audited gap; it carries the exact audited gap as an objective-suboptimality bound only for status certified-gap. The possible statuses are invalid-input, strict-start-unavailable, nonpositive-gap, certified-gap, residual-candidate, singular-system, lost-range, uncertified-newton-solve, boundary-status-unavailable, and budget-exhausted. invalid-input has no iterate. A residual candidate is diagnostic, never a certificate.

Interior audits. For an orthant factor, require every coordinate to be positive. For a Lorentz factor $(t,u)$, require $t>0$ and $t^2-\mathrm{\parallel }u{\mathrm{\parallel }}_2^2>0$. For a symmetric matrix factor $Y$, test positive definiteness by exact rational Schur elimination: require $Y_{11}>0$, form $S=Y_{2:n,2:n}-Y_{2:n,1}{Y}_{1,2:n}\mathrm{/}{Y}_{11}$, and repeat on $S$ until the matrix is empty. Completing the square gives $v^{T}Yv=Y_{11}(v_1+Y_{1,2:n}{v}_{2:n}\mathrm{/}{Y}_{11}{)}^2+v_{2:n}^{T}S{v}_{2:n}$, so this finite test succeeds exactly for $Y\succ 0$.

Instructions. First reject invalid dimensions, nonfinite data, parameters outside their stated domains, or a negative cap with invalid-input. Reject an LP/QP start unless $x,s>0$, or a cone start unless $s,z$ pass the interior audits, with strict-start-unavailable; this is not an infeasibility claim. Set the accepted-update count to zero and repeat the following audit and update loop.

Compute the residuals and complementarity gap exactly as rational expressions in the stored inputs; compare exact squared norms with squared tolerances. A nonpositive or nonfinite gap returns nonpositive-gap. Return certified-gap in LP/QP mode only when $Ax=b$, $x,s\ge 0$, $Qx+c-A^{T}y-s=0$, and $x^{T}s\le {ϵ}_g$; in cone mode require $Ax=b$, $s=h-Gx\in K$, $z\in K^{\ast }$, $A^{T}y-G^{\ast }z=c$, and $⟨s,z⟩\le {ϵ}_g$. Otherwise, if $\mathrm{\parallel }{r}_p{\mathrm{\parallel }}_2\le {ϵ}_p$, $\mathrm{\parallel }{r}_d{\mathrm{\parallel }}_2\le {ϵ}_d$, $\mathrm{\parallel }{r}_c\mathrm{\parallel }\le {ϵ}_d$, and the surrogate gap is at most ${ϵ}_g$, return residual-candidate. If the update count equals $N$, return budget-exhausted.

Assemble the Newton derivative specified in Log-barrier central paths and primal-dual Newton residuals. In LP/QP mode solve its displayed block system. In symmetric-cone mode solve $A\mathrm{\Delta }x=-r_p,A^T\mathrm{\Delta }y-G^{\ast }\mathrm{\Delta }z=-r_d,\mathrm{\Delta }s=-G\mathrm{\Delta }x,\mathrm{\Delta }s\circ z+s\circ \mathrm{\Delta }z=-r_c\mathrm{.}$ The exact equations form a square system after eliminating $\mathrm{\Delta }s$. Form each system coefficient as an exact dyadic expression in the stored inputs, then round that coefficient once to the nearest normal representable value, preserving exact zeros; any failed conversion returns lost-range. Use the row order of the displayed system: $(-r_d,-r_p,-r_c)$ in LP/QP mode and, after eliminating $\mathrm{\Delta }s$, $(-r_p,-r_d,-r_c)$ in symmetric-cone mode. Convert that exact rational right-hand side componentwise by the same nearest-normal rule, preserving exact zeros; a failed conversion also returns lost-range. In LP/QP mode, The LP/QP primal-dual Newton matrix is nonsingular from a strict positive start proves nonsingularity of the exact Newton matrix under the stated rank and strict-start hypotheses; in either mode an exact rank audit detects a singular stored system.

Stored floating-point numbers are treated as exact dyadic rationals (rationals with denominator a power of two) for audits; the unit roundoff $u$ and failure behavior of elementary operations are those of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. For a matrix $M$, $\mathrm{\parallel }M{\mathrm{\parallel }}_F^2$ is the sum of the squares of its entries. Check exact rank by rational row elimination: scan columns from left to right, find the first nonzero entry at or below the next pivot row, swap it into that row, subtract exact multiples to clear entries below, and count pivots. If fewer pivots than rows are found, return singular-system. For a full-rank matrix, use a reduced QR solve: for each active column apply Scaled Householder reflector formation and application to the column and trailing block, store the reflector, apply the stored reflectors in order to the right-hand side using the same sequential dot, scale, multiply, and subtraction order as that procedure, and solve the resulting upper-triangular system by descending substitution $t_i=\mathrm{fl}(q_i-\mathrm{fl}({\sum }_{j>i}{\widehat{R}}_{ij}{\widehat{d}}_j))$, ${\widehat{d}}_i=\mathrm{fl}(t_i\mathrm{/}{\widehat{R}}_{ii})$. A range exception returns lost-range; a zero pivot returns singular-system. Let $d$ be the system dimension. If $(d+3)u\ge 1$, return uncertified-newton-solve. Otherwise define ${\gamma }_{d+3}=(d+3)u\mathrm{/}(1-(d+3)u)$. Require ${\gamma }_{d+3}<1\mathrm{/}8$, and, for the computed triangular factor $\widehat{R}$, verify with exact dyadic arithmetic that ${\gamma }_{d+3}^2\mathrm{\parallel }\widehat{R}{\mathrm{\parallel }}_F^2\mathrm{\parallel }{\widehat{R}}^{-1}{\mathrm{\parallel }}_F^2<1.$ If the test fails, return uncertified-newton-solve. Otherwise Computed reduced QR least squares is backward stable for full column rank applies: the returned direction solves a nearby system with a nonsingular perturbed factor. Compute exactly and report the direction residual both against the rounded stored Newton system and against the original exact rational Newton equations.

For LP/QP mode set $\bar{\alpha }=\min (1,{\min }_{\mathrm{\Delta }{x}_i<0}\frac{-x_i}{\mathrm{\Delta }{x}_i},{\min }_{\mathrm{\Delta }{s}_i<0}\frac{-s_i}{\mathrm{\Delta }{s}_i}),$ where an empty minimum is $+\mathrm{\infty }$, and use $\alpha =\eta \bar{\alpha }$. In symmetric-cone mode compute a positive safe fraction $\beta$ for each primal-slack and dual factor, and use $\alpha =\eta$ times the minimum of those bounds. For an orthant factor use the same coordinate ratios. For a Lorentz factor at $(t,u)\in \mathrm{int}\mathcal{L}$ with direction $(\dot{t},v)$, set $q=t^2-\mathrm{\parallel }u{\mathrm{\parallel }}_2^2$ and $c=\mathrm{\mid }\dot{t}\mathrm{\mid }+{\sum }_i\mathrm{\mid }{v}_i\mathrm{\mid }$; take $\beta =1$ if $c=0$, and otherwise $\beta =\min (1,q\mathrm{/}(2tc))$. For a positive-semidefinite factor at $Y\succ 0$ with symmetric direction $D$, let $L(M)={\sum }_{i,j}\mathrm{\mid }{M}_{ij}\mathrm{\mid }$, compute $Y^{-1}$ by exact rational elimination, and set $C=L(Y^{-1})L(D)$; take $\beta =1$ if $C=0$ and otherwise $\beta =\min (1,1\mathrm{/}C)$. These calculations are finite and use exact rational comparisons; their safety is proved by the fraction-to-boundary theorem on the A page.

Round $\alpha$ downward to a positive normal representable value; if none exists, return lost-range. Update each stored coordinate as $\mathrm{fl}(v+\mathrm{fl}(\alpha \mathrm{\Delta }v))$ under the operationwise model. Any failed operation returns lost-range. Before accepting the candidate update, verify exactly that LP/QP variables $x,s$ remain positive, or that cone variables remain interior and $s=h-Gx$ still holds. If an audit fails, reject the candidate and return boundary-status-unavailable with the previous accepted iterate. Otherwise accept it, increment the update count, and return to the audit loop. Every factorization, inverse, rank, and safe-fraction calculation is bounded by finite matrix dimensions, and the outer loop is bounded by $N$. No failure status is an infeasibility conclusion; Phase I is a separate model.