cvx-quadprog¶

A pure NumPy/SciPy implementation of the Goldfarb/Idnani dual active-set method for strictly convex quadratic programs. It is a reimplementation of quadprog, which wraps C code descended from Berwin Turlach's Fortran translation of the original algorithm.

No compiler, no Cython, no build step — just NumPy and SciPy.

The problem¶

Minimise

$$\tfrac{1}{2} x^T G x - a^T x \quad \text{subject to} \quad C^T x \ge b$$

with G symmetric positive definite. The first meq constraints are treated as equalities.

Note the two conventions inherited from the original: the linear term is subtracted, and constraints are given column-wise (C is n × m, one column per constraint) as >=.

Usage¶

import numpy as np
from cvx.quadprog import solve_qp

G = np.eye(3)
a = np.array([0.0, 5.0, 0.0])
C = np.array([[-4.0, 2.0, 0.0], [-3.0, 1.0, -2.0], [0.0, 0.0, 1.0]])
b = np.array([-8.0, 2.0, 0.0])

solution = solve_qp(G, a, C, b)

solution.x  # array([0.47619048, 1.04761905, 2.09523810])
solution.f  # -2.380952380952381
solution.xu  # array([0., 5., 0.])  the unconstrained minimiser
solution.iterations  # array([3, 0])  constraints added, constraints dropped
solution.lagrangian  # array([0., 0.23809524, 2.09523810])
solution.iact  # array([3, 2])  1-based indices of the active set

Solution is a NamedTuple yielding those six values in the order returned by quadprog.solve_qp, so existing tuple-unpacking code keeps working:

x, f, xu, iterations, lagrangian, iact = solve_qp(G, a, C, b)

If C and b are omitted the unconstrained problem is solved. Passing factorized=True means G holds $R^{-1}$ rather than $G$, where $G = R^T R$ with R upper triangular — useful when a cheaper factorisation is available, for instance when G is banded.

Infeasible constraints, a non-positive-definite G, and inconsistent shapes all raise ValueError.

Why the dual method¶

The algorithm starts at the unconstrained minimum $G^{-1} a$, which is dual feasible by construction, and adds the most violated constraint one at a time. Every iterate stays dual feasible, so the objective increases monotonically and no phase-1 feasibility problem is required. Constraints whose multipliers would turn negative are dropped along the way.

The factorisation of the active constraint normals is carried between iterations and updated orthogonally rather than recomputed, which is what makes each iteration $O(n^2)$ and the method numerically stable. Insertions use a Householder reflection and deletions a Givens chase — see Performance.

Agreement with the C implementation¶

tests/test_against_c.py runs both implementations on the same problems and compares every return value. Across a wider sweep of 4000 random problems (2 ≤ n ≤ 11, up to 14 constraints, mixed equalities):

Matching the iteration counts exactly means the two follow the same active-set path, adding and dropping the same constraints in the same order — a much stronger statement than agreeing on the final answer.

Deliberate deviations¶

• Cholesky and triangular inversion use LAPACK (via SciPy) instead of the hand-rolled routines in linear-algebra.c. A matrix that is positive definite only marginally may therefore be accepted by one and rejected by the other. Input arrays are not scanned for NaN/inf, matching the reference; a non-finite G surfaces as the "not positive definite" error.
• Constraint insertion uses a Householder reflection rather than a chain of Givens rotations, so Q and R differ by column and row signs. See Performance for why the solver is indifferent to this.
• Inputs are never destroyed. The C routine overwrites G and a.
• R uses the reference's packed-column layout, for the reason given under Performance — not merely to halve the memory.
• Summation order differs wherever a loop became a NumPy dot product, so results agree to floating-point tolerance rather than bit for bit. The objective is accumulated incrementally by both, as in the original. Measuring each against a direct re-evaluation at its own minimiser over 2164 problems, the worst-case drift is somewhat smaller here — 1.5e-8 absolute (7.4e-15 relative) against 3.7e-8 (1.8e-14) — but neither dominates problem by problem: the reference is the closer of the two on 801 problems, this implementation on 782, with 581 ties.
• Extra validation: meq is range-checked, and passing C without b is an error rather than a crash.

Where the two may legitimately differ¶

Duplicated or linearly dependent constraints make the dual solution non-unique: the multiplier can sit on either copy. Both implementations return a valid KKT point, but not necessarily the same one, and lagrangian/iact differ accordingly. x and f are unaffected. tests/test_against_c.py covers this case by verifying the KKT conditions rather than demanding an identical dual.

Performance¶

Box-constrained problems (n variables, 2n constraints), per solve. Timings are the best of five batches, after a warm-up call, on an arm64 machine with Python 3.12 / NumPy 2.5.1 against quadprog 0.1.13:

The crossover sits at n ≈ 160 — measured by sweeping the interval, where the ratio passes 1.0 between n = 150 (1.03×) and n = 160 (0.94×). Below it, cost is dominated by per-call NumPy dispatch: about 15 µs per iteration spread over roughly 18 array operations, against ~6 µs for C to do an entire n = 10 solve. That is a floor set by the interpreter, not by the algorithm.

Above the crossover this implementation wins, because the reference's linear-algebra.c uses hand-rolled scalar loops for its dot products and axpys, while the work here is expressed as BLAS calls that reach tuned, vectorised kernels.

How the inner loop is organised¶

The reference reduces each incoming constraint normal with a chain of Givens rotations — one per trailing component, each touching every row of Q. In Python that is O(n) interpreter round-trips per insertion, and it dominated everything else (85% of runtime at n = 150).

qr_insert instead applies a single Householder reflection, which performs the same reduction in one matrix-vector product plus one rank-1 update. The rank-1 update goes through BLAS dger directly into Q's buffer, so no O(n·k) temporary is allocated.

This is safe despite producing a different Q and R than the reference (some diagonal signs differ), because the quantities the solver consumes are invariant to the choice of reduction:

$$rv = R^{-1} d_1 = (A^T G^{-1} A)^{-1} A^T G^{-1} n$$

depends only on A, n and G. Replacing R by SR for a sign matrix S also replaces d₁ by Sd₁, and the two cancel exactly. zv = J_2 d_2 is invariant for the same reason. The measured iteration counts confirm it: they still match the reference exactly on every problem tested, including n up to 220 in the test suite.

qr_delete keeps the Givens chase, which is inherently sequential — each rotation's parameters depend on the previous one having been applied.

The triangular solve¶

Each iteration solves R rv = d₁ for the dual step direction. With R held as a dense (r, r) array, the active block R[:nact, :nact] is a strided view, so handing it to LAPACK forces a full copy — about 1 MB per iteration at n = 700. The copy, not the arithmetic, was the cost:

So R uses the reference's packed-column layout instead: column j is j + 1 contiguous values at offset j(j+1)/2, which makes the leading nact triangle the leading nact(nact+1)/2 entries — contiguous by construction, and readable in place by BLAS tpsv with no copy at any active-set size.

Measured over the whole solve at n = 700, that operation went from 15.5 ms (21% of runtime) to 2.0 ms (3.5%). The cost is borne by qr_delete, which mixes two rows across a range of columns: column offsets grow, so that becomes a gather rather than a slice.

Constraint structure¶

A bound constraint is one nonzero in its column of C, and a box-constrained problem is nothing but bounds. Three of the per-iteration products then stop being reductions and become indexing, so solve_qp detects the structure once, per column:

Detection is per column rather than all-or-nothing because the useful case is mixed: mean-variance carries a dense budget column (Σx = 1) beside 2n bounds. An all-or-nothing test would see that one dense column and send the whole problem down the slow path. The slack product has its own three-way choice — all unit, sparse (a compiled CSR product), or dense.

This is a fast path around arithmetic the dense path would do anyway, so it cannot change the answer, and the differential tests against the C implementation cover box, mixed budget-plus-bounds, and fully dense C.

Where the remaining time goes at n = 700, after both optimisations:

Keeping Q implicit: measured, and rejected¶

qr_insert is now the whole game, and it updates J explicitly on every insertion. The obvious next move is the one LAPACK's geqrf/ormqr make: store the Householder vectors and never form Q. Insertion then costs nothing at all, because dv already is the reduced column — the reflection is read off it and appended.

It was prototyped, checked against this implementation on 300 problems (identical iteration counts, worst |Δx| 4.7e-12), and measured. It is 2.4–2.6× slower, even with zero deletions:

The flop count says why. Per iteration at active size k:

Removing the insertion does not remove its work, it relocates it. Applying an implicitly stored Q costs O(nk), and the solver applies it twice per iteration — which is exactly what the explicit update pays once. Forming J amortises the accumulated Q into a single dense matrix, so every later application is one gemv regardless of k. That is the whole reason to form it. Implicit storage wins when Q is applied rarely relative to the number of reflections; here it is applied twice per reflection, which is the worst case.

Deletion is the second, independent objection. A Givens chase cannot be absorbed into a stored Householder chain, so a deletion becomes a refactorisation — 82% of runtime on a problem with 200 of them, and 2.48× slower overall. Deletions are rare in practice (0% of steps on box and budget-plus-bounds problems, 2.2% on random dense C), so a hybrid would have been viable had the insertion side won — but it does not.

Accuracy is unaffected. Over 3000 random problems the worst relative KKT stationarity residual is 7.3e-13 here against 8.8e-13 for the reference, and this implementation is strictly the more accurate of the two on 1035 problems to the reference's 869.

Layout¶

src/cvx/quadprog/_solve.py   the dual active-set iteration
src/cvx/quadprog/_qr.py      Givens QR insert/delete
tests/test_reference.py      the upstream test suite, ported
tests/test_qr.py             QR update invariants, in isolation
tests/test_against_c.py      differential test vs. the C implementation

947 tests, 100% line and branch coverage of src/.

Reference¶

D. Goldfarb and A. Idnani (1983). A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Programming, 27, 1–33.