Coverage for src/cvx/quadprog/_solve.py: 100%
160 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-06 11:28 +0000
1"""The Goldfarb/Idnani dual active-set method for strictly convex QPs.
3A NumPy/SciPy reimplementation of the ``quadprog`` package, which itself
4descends from Berwin Turlach's Fortran translation of the algorithm in [1].
6The method is *dual* feasible throughout: it starts at the unconstrained
7minimum, which satisfies the dual conditions trivially, and drives the primal
8infeasibility to zero one constraint at a time. Because every iterate is dual
9feasible, the objective increases monotonically and the iteration count is
10bounded by the number of constraints -- no phase-1 problem is needed.
12References:
13 [1] D. Goldfarb and A. Idnani (1983). A numerically stable dual method for
14 solving strictly convex quadratic programs. Mathematical Programming,
15 27, 1-33.
16"""
18# G, C, R and J are the names used in Goldfarb & Idnani (1983) and in the
19# reference implementation's public signature `solve_qp(G, a, C, b, meq)`.
20# Lowercasing them would obscure the correspondence to the paper and break
21# drop-in compatibility, so the pep8-naming rules are waived here. TRY003 goes
22# with them: the ValueError messages are reproduced verbatim from the reference
23# so that callers matching on the text keep working.
24#
25# These live in the files rather than in a [lint.per-file-ignores] block because
26# ruff.toml is template-owned -- a local edit to it is reverted by the next
27# `/rhiza:update` sync and flagged as non-template by stage_synced.py.
28# ruff: noqa: N803, N806, TRY003
30from collections.abc import Callable
31from typing import Any, NamedTuple, cast
33import numpy as np
34import scipy.linalg
35import scipy.sparse
37from ._qr import qr_delete, qr_insert
39__all__ = ["Solution", "solve_qp"]
41# scipy resolves its LAPACK wrappers at run time, so they carry no useful static
42# signature: scipy-stubs types get_lapack_funcs as returning a function *or* a
43# list of them, depending on whether one name or a sequence was asked for. We ask
44# for one name, so it is one function -- the cast records that rather than
45# leaving every call site to assert it.
46_LapackFn = Callable[..., Any]
48# Prototype array fixing the precision the wrappers are resolved for.
49_F64 = np.empty(0, dtype=np.float64)
52def _calculate_vsmall() -> float:
53 """Return an upper bound on the relative precision of the arithmetic.
55 Gleaned from Powell's ZQPCVX routine: double the value until it is large
56 enough to perturb 1.0 when scaled by both 0.1 and 0.2. Computed once at
57 import time.
59 Returns:
60 A small positive number, of the order of the machine epsilon.
61 """
62 vsmall = 1e-60
63 while True:
64 vsmall += vsmall
65 if vsmall * 0.1 + 1.0 > 1.0 and vsmall * 0.2 + 1.0 > 1.0:
66 return vsmall
69VSMALL = _calculate_vsmall()
71# Packed triangular solve, resolved once. This runs once per iteration and is
72# the reason R is stored packed: `ap` is an unshaped rank-1 argument, so passing
73# the whole array with n=nact reads the leading triangle in place. The dense
74# equivalent, trtrs on R[:nact, :nact], is handed a strided view and copies it
75# every call -- 77 us against 7.5 us at n = 700. Calling BLAS directly also
76# skips scipy.linalg.solve_triangular's per-call validation, whose check_finite
77# scans the whole array.
78_TPSV = cast("_LapackFn", scipy.linalg.get_blas_funcs("tpsv", (_F64,)))
80# Triangular inverse. Resolved the same way rather than reached as
81# scipy.linalg.lapack.dtrtri: the per-precision wrappers are generated at import
82# time, so no static checker can see them, and get_lapack_funcs is the documented
83# entry point that also picks the precision to match the input.
84_TRTRI = cast("_LapackFn", scipy.linalg.get_lapack_funcs("trtri", (_F64,)))
86# Returned for the dual step direction while the active set is still empty.
87_EMPTY = np.zeros(0)
90class Solution(NamedTuple):
91 """The outcome of a quadratic program.
93 Iterating over an instance yields the same six values, in the same order, as
94 the tuple returned by ``quadprog.solve_qp``, so it is a drop-in replacement.
96 Attributes:
97 x: ``(n,)`` minimiser of the constrained problem.
98 f: Value of the objective at ``x``.
99 xu: ``(n,)`` minimiser of the unconstrained problem, ``G^-1 a``.
100 iterations: ``(2,)`` count of constraints added to the active set (once
101 per outer iteration) and of constraints removed from it.
102 lagrangian: ``(m,)`` Lagrange multipliers, zero for inactive
103 constraints.
104 iact: 1-based indices of the constraints active at the solution.
105 """
107 x: np.ndarray
108 f: float
109 xu: np.ndarray
110 iterations: np.ndarray
111 lagrangian: np.ndarray
112 iact: np.ndarray
115def solve_qp(
116 G: np.ndarray,
117 a: np.ndarray,
118 C: np.ndarray | None = None,
119 b: np.ndarray | None = None,
120 meq: int = 0,
121 factorized: bool = False,
122) -> Solution:
123 r"""Solve a strictly convex quadratic program.
125 .. math::
126 \min_x \tfrac{1}{2} x^T G x - a^T x \quad\text{subject to}\quad C^T x \ge b
128 The first ``meq`` constraints are treated as equalities.
130 Args:
131 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
132 If ``factorized`` is True, pass :math:`R^{-1}` instead, where
133 :math:`G = R^T R` with ``R`` upper triangular.
134 a: ``(n,)`` vector of the linear term.
135 C: ``(n, m)`` constraint matrix, one column per constraint. Defaults to
136 a single inactive constraint, giving the unconstrained problem.
137 b: ``(m,)`` right-hand side of the constraints.
138 meq: Number of leading constraints to treat as equalities.
139 factorized: Whether ``G`` holds :math:`R^{-1}` rather than :math:`G`.
141 Returns:
142 A :class:`Solution` with the minimiser, the objective value, the
143 unconstrained minimiser, the iteration counts, the Lagrange multipliers
144 and the active set.
146 Raises:
147 ValueError: If the shapes are inconsistent, if ``meq`` is out of range,
148 if ``G`` is not positive definite, or if the constraints admit no
149 solution.
150 """
151 G = np.asarray(G, dtype=np.float64)
152 a = np.asarray(a, dtype=np.float64)
154 if C is None and b is None:
155 C, b, meq = np.zeros((len(G), 1)), -np.ones(1), 0
156 elif C is None or b is None:
157 raise ValueError("C and b must be given together")
159 C = np.asarray(C, dtype=np.float64)
160 b = np.asarray(b, dtype=np.float64)
162 n, q = _validate(G, a, C, b, meq)
163 r = min(n, q)
165 # Initialisation. We want xv to hold G^-1 a, the unconstrained minimum, and
166 # J to hold R^-1, so that J J^T = G^-1.
167 J, xv = _factorize(G, a, factorized)
169 # The objective at the unconstrained minimum. Kept as a running total: each
170 # step updates it in closed form rather than re-evaluating the quadratic.
171 obj = -float(a @ xv) / 2.0
172 xu = xv.copy()
174 # The norm of each column of C, used to scale the pivoting rule so that the
175 # choice of constraint is invariant to how each one happens to be scaled.
176 # A zero-norm column reads 0 >= b, which no x can influence; scoring it as
177 # infinitely violated sends the solver to the infeasibility verdict instead
178 # of dividing by zero below.
179 nbv = np.sqrt(np.sum(C * C, axis=0))
180 degenerate = nbv == 0.0
181 nbv_safe = np.where(degenerate, 1.0, nbv)
183 # Sparsity of C, detected once. Bound constraints make most columns a single
184 # scaled unit vector, which turns three of the per-iteration operations from
185 # O(n) or O(n*q) work into scalar indexing.
186 single, srow, sval = _analyse_constraints(C)
187 slack_of = _slack_evaluator(C, single)
189 # Upper triangular, stored as packed columns -- see the note in _qr.
190 R = np.zeros(r * (r + 1) // 2)
191 uv = np.zeros(r) # dual variables of the active constraints
192 iact = np.zeros(q, dtype=np.int64) # 1-based, first nact entries valid
193 lagr = np.zeros(q)
194 nact = 0
195 iter_full, iter_partial = 0, 0
197 while True:
198 iter_full += 1
200 # The slack of every constraint. Slacks of active constraints are forced
201 # to exactly zero as a safeguard against rounding error.
202 sv = slack_of(xv) - b
203 sv[np.abs(sv) < VSMALL] = 0.0
204 sv[iact[:nact] - 1] = 0.0
206 iadd = _choose_constraint(sv, nbv_safe, degenerate, meq)
208 if iadd == 0:
209 # Every constraint is satisfied, so we are at the optimum.
210 lagr[iact[:nact] - 1] = uv[:nact]
211 iterations = np.array([iter_full, iter_partial], dtype=np.int64)
212 return Solution(xv, obj, xu, iterations, lagr, iact[:nact])
214 # An equality constraint may be violated from either side. When its
215 # slack is positive we have to step in the opposite direction.
216 slack = float(sv[iadd - 1])
217 reverse_step = slack > 0.0
218 u = 0.0
220 # A column holding a single scaled unit vector e_row lets the three
221 # products against it below be read off by index instead of computed.
222 unit = bool(single[iadd - 1])
223 if unit:
224 row, val = int(srow[iadd - 1]), float(sval[iadd - 1])
225 else:
226 normal = C[:, iadd - 1]
228 # Inner loop: walk towards the constraint boundary, dropping active
229 # constraints whose multipliers would otherwise turn negative.
230 while True:
231 # dv = J^T n, split as (d_1, d_2) at the size of the active set.
232 # For a unit column this is one scaled row of J, O(n) not O(n^2).
233 dv = val * J[row, :] if unit else J.T @ normal
235 # zv = J_2 d_2 is the step direction of the primal variable, the
236 # component of the constraint normal orthogonal to the active set.
237 zv = J[:, nact:] @ dv[nact:]
239 # rv = R^-1 d_1 is the negated step direction of the dual variable.
240 # Solved on a copy: dv is still needed intact for qr_insert below.
241 rv = _TPSV(nact, R, dv[:nact].copy(), overwrite_x=True) if nact else _EMPTY
243 # The largest step t1 that keeps the dual variables non-negative,
244 # and the constraint idel that would be the first to bind at zero.
245 t1, idel = _dual_step_limit(uv, rv, iact, nact, meq, reverse_step)
246 t1inf = idel == 0
248 # The step t2 that brings the slack of the entering constraint to
249 # zero. ztn is the rate of change of that slack.
250 t2inf = abs(float(zv @ zv)) <= VSMALL
251 if not t2inf:
252 ztn = val * float(zv[row]) if unit else float(zv @ normal)
253 t2 = abs(slack) / ztn
255 if t1inf and t2inf:
256 # We can step infinitely far: the dual is unbounded, so the
257 # primal is infeasible.
258 raise ValueError("constraints are inconsistent, no solution")
260 full_step = not t2inf and (t1inf or t1 >= t2)
261 step_length = t2 if full_step else t1
262 step = -step_length if reverse_step else step_length
264 if not t2inf:
265 xv += step * zv
266 obj += step * ztn * (step / 2.0 + u)
268 uv[:nact] -= step * rv
269 u += step
271 if full_step:
272 break
274 # Only a partial step: drop constraint idel from the active set.
275 qr_delete(nact, idel, J, R)
276 uv[idel - 1 : nact - 1] = uv[idel:nact].copy()
277 iact[idel - 1 : nact - 1] = iact[idel:nact].copy()
278 uv[nact - 1], iact[nact - 1] = 0.0, 0
279 nact -= 1
280 iter_partial += 1
282 if not t2inf:
283 # We moved in primal space, so the slack we are closing has
284 # changed and must be recomputed.
285 reached = val * float(xv[row]) if unit else float(xv @ normal)
286 slack = reached - float(b[iadd - 1])
288 # The entering constraint now holds with equality: add it.
289 nact += 1
290 uv[nact - 1], iact[nact - 1] = u, iadd
291 qr_insert(nact, dv, J, R)
294def _analyse_constraints(C: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
295 """Locate the columns of ``C`` that hold a single scaled unit vector.
297 Bound constraints -- the overwhelmingly common shape, and what
298 ``C = [I, -I]`` is -- make every column one nonzero. Recognising that lets
299 the products against a constraint normal become scalar indexing rather than
300 length-``n`` reductions.
302 Args:
303 C: ``(n, m)`` constraint matrix, one column per constraint.
305 Returns:
306 A boolean mask of the single-nonzero columns, the row index of that
307 nonzero per column, and its value. Entries of the latter two are
308 meaningless where the mask is False.
309 """
310 nonzero = C != 0.0
311 single = nonzero.sum(axis=0) == 1
312 # argmax gives the first nonzero row, which for a single-nonzero column is
313 # the only one. Columns failing the mask still index safely, just uselessly.
314 row = np.argmax(nonzero, axis=0)
315 return single, row, C[row, np.arange(C.shape[1])]
318def _slack_evaluator(C: np.ndarray, single: np.ndarray) -> Callable[[np.ndarray], np.ndarray]:
319 """Return the cheapest available way to evaluate ``C.T @ x``.
321 This runs once per outer iteration over all ``m`` constraints, so on a
322 box-constrained problem -- where ``m = 2n`` -- the dense product is the
323 single largest cost in the solver. Three strategies, in preference order:
325 * every column a single nonzero: one gather, ``O(m)``;
326 * sparse enough to pay for the bookkeeping: a CSR product, ``O(nnz)``;
327 * otherwise the dense product, transposed once here rather than per call.
329 Args:
330 C: ``(n, m)`` constraint matrix.
331 single: Mask of the columns holding exactly one nonzero.
333 Returns:
334 A callable mapping ``x`` to ``C.T @ x``.
335 """
336 n, m = C.shape
338 if m and single.all():
339 nonzero = C != 0.0
340 row = np.argmax(nonzero, axis=0)
341 val = C[row, np.arange(m)]
342 return lambda x: val * x[row]
344 if np.count_nonzero(C) * 4 <= n * m:
345 # csr_matrix multiplication is compiled, so this beats the dense product
346 # well before the matrix is especially sparse.
347 ct = scipy.sparse.csr_matrix(C.T)
348 return lambda x: ct @ x
350 dense = np.ascontiguousarray(C.T)
351 return lambda x: dense @ x
354def _validate(G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray, meq: int) -> tuple[int, int]:
355 """Check that the problem data is dimensionally consistent.
357 Args:
358 G: ``(n, n)`` matrix of the quadratic term.
359 a: ``(n,)`` vector of the linear term.
360 C: ``(n, m)`` constraint matrix.
361 b: ``(m,)`` right-hand side of the constraints.
362 meq: Number of leading constraints treated as equalities.
364 Returns:
365 The number of variables and the number of constraints.
367 Raises:
368 ValueError: If any shape disagrees, or if ``meq`` is out of range.
369 """
370 if G.ndim != 2 or G.shape[0] != G.shape[1]:
371 raise ValueError(f"G must be a square matrix. Received shape={G.shape}")
372 n = G.shape[0]
373 if a.shape != (n,):
374 raise ValueError(f"G and a must have the same dimension. Received G as {G.shape} and a as {a.shape}")
375 if C.ndim != 2 or C.shape[0] != n:
376 raise ValueError(f"G and C must have the same first dimension. Received G as {G.shape} and C as {C.shape}")
377 q = C.shape[1]
378 if b.shape != (q,):
379 raise ValueError(
380 f"The number of columns of C must match the length of b. Received C as {C.shape} and b as {b.shape}"
381 )
382 if not 0 <= meq <= q:
383 raise ValueError(f"meq must satisfy 0 <= meq <= {q}. Received {meq}")
384 return n, q
387def _factorize(G: np.ndarray, a: np.ndarray, factorized: bool) -> tuple[np.ndarray, np.ndarray]:
388 """Return the inverse Cholesky factor of ``G`` and the unconstrained minimum.
390 Args:
391 G: ``(n, n)`` positive definite matrix, or its inverse Cholesky factor
392 :math:`R^{-1}` when ``factorized`` is True.
393 a: ``(n,)`` vector of the linear term.
394 factorized: Whether ``G`` already holds :math:`R^{-1}`.
396 Returns:
397 ``J``, an upper triangular array with ``J J^T = G^-1``, and the
398 unconstrained minimiser ``G^-1 a``. ``J`` is a fresh writable array; the
399 caller updates it in place.
401 Raises:
402 ValueError: If ``G`` is not positive definite.
403 """
404 # Fortran order throughout: the updates in _qr work on column blocks of J,
405 # which are then contiguous and can go straight to BLAS.
406 if factorized:
407 J = np.asfortranarray(np.triu(G))
408 return J, J @ (J.T @ a)
410 # check_finite=False skips a full scan of each array on the way in. A
411 # non-finite G surfaces as the "not positive definite" error below instead of
412 # its own exception, which matches the reference: it does not check either.
413 try:
414 R = scipy.linalg.cholesky(G, lower=False, check_finite=False)
415 except scipy.linalg.LinAlgError as exc:
416 raise ValueError("matrix G is not positive definite") from exc
418 xv = scipy.linalg.cho_solve((R, False), a, check_finite=False)
419 J, info = _TRTRI(R, lower=0)
420 if info != 0: # pragma: no cover
421 # Defensive: trtri fails only on an exactly zero diagonal entry, which
422 # a successful Cholesky has already ruled out.
423 raise ValueError("matrix G is not positive definite")
424 return np.asfortranarray(np.triu(J)), xv
427def _choose_constraint(sv: np.ndarray, nbv_safe: np.ndarray, degenerate: np.ndarray, meq: int) -> int:
428 """Return the 1-based index of the most violated constraint, or 0 if none.
430 Violations are measured relative to the norm of the constraint normal, so
431 the choice does not depend on the scaling of individual constraints. An
432 equality constraint counts as violated in either direction.
434 Scanning for the largest violation is equivalent to taking an ``argmax``,
435 which resolves ties towards the lowest index just as a left-to-right scan
436 with a strict improvement test does.
438 Args:
439 sv: Slack of each constraint.
440 nbv_safe: Norm of each constraint normal, with zeros replaced by one.
441 degenerate: Mask of the constraints whose normal has zero norm.
442 meq: Number of leading constraints treated as equalities.
444 Returns:
445 The 1-based index of the constraint to add, or 0 at the optimum.
446 """
447 # An inequality is violated when its slack is negative; an equality whenever
448 # its slack is nonzero.
449 violation = -sv
450 np.abs(violation[:meq], out=violation[:meq])
452 score = violation / nbv_safe
453 # A zero-norm normal cannot be satisfied by any step, so rank it first.
454 if degenerate.any():
455 score = np.where(degenerate & (violation > 0.0), np.inf, score)
457 iadd = int(np.argmax(score))
458 return iadd + 1 if score[iadd] > 0.0 else 0
461def _dual_step_limit(
462 uv: np.ndarray,
463 rv: np.ndarray,
464 iact: np.ndarray,
465 nact: int,
466 meq: int,
467 reverse_step: bool,
468) -> tuple[float, int]:
469 """Return the largest dual-feasible step and the constraint that limits it.
471 Stepping along ``-rv`` drives the multipliers of the active inequality
472 constraints towards zero. The first one to reach zero caps the step, since a
473 negative multiplier would be dual infeasible. Equality constraints are
474 exempt: their multipliers are unrestricted in sign.
476 Args:
477 uv: Dual variables of the active constraints.
478 rv: Negated step direction of the dual variables.
479 iact: 1-based indices of the active constraints.
480 nact: Size of the active set.
481 meq: Number of leading constraints treated as equalities.
482 reverse_step: Whether the step is taken in the negative direction.
484 Returns:
485 The step limit and the 1-based position in the active set of the
486 constraint that attains it. The position is 0 when no constraint limits
487 the step, in which case the limit is meaningless.
488 """
489 # Working with the signed direction lets one comparison serve both cases and
490 # makes the eligible entries positive, so no separate abs is needed.
491 direction = -rv[:nact] if reverse_step else rv[:nact]
492 eligible = (iact[:nact] > meq) & (direction > 0.0)
493 if not eligible.any():
494 return 0.0, 0
496 # The inner where keeps the division clear of the ineligible entries; the
497 # outer one pushes them above any real ratio so argmin skips them.
498 ratio = np.where(eligible, uv[:nact] / np.where(eligible, direction, 1.0), np.inf)
499 idel = int(np.argmin(ratio))
500 return float(ratio[idel]), idel + 1