Coverage for src/cvx/quadprog/_pdas.py: 100%
151 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 06:10 +0000
1r"""A primal-dual active-set fast path, checked against the KKT conditions.
3The dual method in :mod:`._solve` reaches the active set one constraint at a
4time, so its iteration count grows with the size of that set: a
5budget-plus-bounds problem at ``n = 100`` takes 74 outer iterations, and at these
6sizes a solve costs very nearly its count of interpreter-level operations rather
7than its flops.
9This module guesses the whole active set instead, solves one dense KKT system for
10it, and repairs the guess from the signs that come back -- a multiplier wanting to
11go negative does not belong in the set, a violated constraint does. That converges
12in two to four repairs across every family measured, independent of ``n``, so it
13trades flops, which are free at these sizes, for dispatches, which are not.
15The trade is only sound because the answer is *checked*. Primal-dual active set is
16not globally convergent: it can cycle, it can stabilise on a set whose KKT system
17is singular, and -- measured, not hypothesised -- it can stabilise on a point that
18is simply not the minimiser. Every candidate therefore goes through
19:func:`_certified` before it is returned, and anything that fails is discarded so
20:func:`~cvx.quadprog.solve_qp` can fall back to the exact walk. For a strictly
21convex program the KKT conditions are sufficient, so a point that passes is the
22unique minimiser and needs no second opinion.
23"""
25# G and C are the names used in Goldfarb & Idnani (1983) and in the reference
26# implementation's public signature, and CA, Y are the matrices of the KKT system
27# as it is conventionally written. Lowercasing them here would break the
28# correspondence with _solve.py, which waives the same rules for the same reason.
29# ruff: noqa: N803, N806
31from typing import NamedTuple
33import numpy as np
34import scipy.linalg as sla
36from ._base import Solution
38# Below this many variables the fast path is not attempted. Its guess is likeliest
39# to be linearly dependent when there are few variables to spread the active set
40# over, and that is also where it wins least: measured over 300 instances per
41# size, the certified fraction is 100% from twelve variables up on every family
42# tried, against 23% at n = 3 and 48% at n = 5 on equality-constrained problems,
43# where the exact walk costs under 100 us anyway.
44_MIN_VARIABLES = 12
46# Cap on set repairs before the attempt is abandoned. Two to four is the observed
47# range; anything near this bound is not converging and is better handed over.
48_MAX_REPAIRS = 30
50# Relative tolerance for deciding set membership. This one is not delicate: it
51# steers which set is tried next, and a bad choice costs a repair or a fallback,
52# never a wrong answer -- that is what the certificate is for.
53_SET_TOL = 1e-10
55# Relative tolerance of the certificate itself, which *is* delicate, since it is
56# the only thing standing between a non-optimal point and the caller. Measured
57# over 1164 converged attempts, points that satisfied the conditions did so with a
58# residual of at most 4e-15, while the two that did not missed by 1e-1 and worse.
59# Anything between those extremes separates them; this sits six orders above the
60# worst good residual and eight below the best bad one.
61_CERTIFY_TOL = 1e-9
64# Work, as n * k, above which the working-set system is formed from one half of
65# the Cholesky factorisation rather than by applying both. The split form does
66# half the flops of the two-sided one but makes five scipy calls where it makes
67# three, so which wins is a question of size. Timed on the system alone at
68# k = n/3, the split form loses by 22% at n = 25 (n * k = 200), by 11% at n = 50
69# and by 17% at n = 75 (1875), then wins by 8% at n = 100 (3300), 28% at n = 200
70# and 40% at n = 400.
71#
72# The value is not fitted to that crossover, because it does not have to be. The
73# n * k a real solve presents is far from continuous: over box, budget-plus-bounds
74# and dense-C families it is at most 5600 at n = 100 and at least 9800 at n = 200,
75# so every gate in that gap sends the same solves down the same paths. This sits
76# in the middle of the gap, where a family whose active set is a little larger or
77# smaller than those measured does not change which side it falls on. End to end
78# the difference at n <= 100 is in any case below what this machine can resolve;
79# the gate is there so that the flop-free case cannot pay for dispatches it does
80# not need, not because the small sizes measured a loss.
81_SPLIT_MIN_WORK = 6000
83# Ceiling on the column cache below, in entries. The cache has the same shape as
84# C, so it can at most double what the solver holds for the constraints, and this
85# bounds that trade at 128 MB rather than letting it scale with the problem. Above
86# the ceiling the cache is simply not built and every repair re-solves, which is
87# the behaviour this constant replaced.
88_CACHE_MAX_ENTRIES = 16_000_000
91class _Reuse(NamedTuple):
92 """Per-solve state that a repair can read instead of recomputing.
94 Two things survive from one repair to the next. ``C.T @ xu`` is fixed for the
95 whole attempt, and every repair needs the entries of it its working set names;
96 rebuilding those from ``C_A`` copied an ``n`` by ``k`` block of ``C`` per
97 repair, 0.11 ms at ``n = 800``, to arrive at numbers already computed while
98 seeding.
100 The columns of ``U^-T C`` are the larger saving. Each repair solves for its
101 own working set, and consecutive working sets overlap heavily -- measured over
102 five instances per cell, 58% to 69% of the columns a solve asks for on box and
103 dense-C families were solved on an earlier repair, and 29% to 37% on
104 budget-plus-bounds. Those solves are the largest single cost in this module,
105 ``n^2 k`` flops against ``n k^2`` for the dual Hessian, so the repeated ones
106 are worth keeping.
108 Reuse is exact, not approximate: a column of ``U^-T C`` does not depend on the
109 working set it was solved for, so a cached column is the column the repair
110 would have computed. ``Z`` and ``have`` are None together on a problem where
111 holding them would cost more memory than :data:`_CACHE_MAX_ENTRIES` allows,
112 and every column is then re-solved, which is what this class replaced.
113 """
115 ctxu: np.ndarray
116 Z: np.ndarray | None
117 have: np.ndarray | None
120def _reuse_state(n: int, m: int, ctxu: np.ndarray) -> _Reuse:
121 """Return the reusable state for one attempt, with or without a column cache.
123 Args:
124 n: Number of variables.
125 m: Number of constraints.
126 ctxu: ``C.T @ xu``, computed while seeding.
128 Returns:
129 A :class:`_Reuse` holding an empty ``(n, m)`` cache, or one whose cache is
130 None when that shape exceeds :data:`_CACHE_MAX_ENTRIES`.
131 """
132 if n * m > _CACHE_MAX_ENTRIES:
133 return _Reuse(ctxu, None, None)
134 return _Reuse(ctxu, np.empty((n, m)), np.zeros(m, dtype=bool))
137class Attempt(NamedTuple):
138 """A candidate solution that has already passed the KKT certificate.
140 Attributes:
141 x: ``(n,)`` minimiser.
142 xu: ``(n,)`` unconstrained minimiser, ``G^-1 a``.
143 lagrangian: ``(m,)`` multipliers, zero off the active set.
144 active: ``(m,)`` boolean mask of the active constraints.
145 added: Constraints added to the working set, summed over all repairs.
146 dropped: Constraints removed from it, summed over all repairs.
147 """
149 x: np.ndarray
150 xu: np.ndarray
151 lagrangian: np.ndarray
152 active: np.ndarray
153 added: int
154 dropped: int
157def attempt(G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray, meq: int) -> Attempt | None:
158 """Try to solve by primal-dual active set, returning None if anything is off.
160 Every rejection path -- too small, singular, cycling, not converging, or
161 failing the certificate -- returns None rather than raising, because the
162 caller's response to all of them is the same: solve it the exact way.
164 Args:
165 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
166 a: ``(n,)`` vector of the linear term.
167 C: ``(n, m)`` constraint matrix, one column per constraint.
168 b: ``(m,)`` right-hand side of ``C.T @ x >= b``.
169 meq: Number of leading constraints held as equalities.
171 Returns:
172 A certified :class:`Attempt`, or None if the fast path did not produce
173 one.
174 """
175 seeded = _seed(G, a, C, b, meq)
176 if seeded is None:
177 return None
178 cho, xu, active, scale, ctxu = seeded
179 m = C.shape[1]
180 reuse = _reuse_state(C.shape[0], m, ctxu)
182 seen: set[bytes] = set()
183 added, dropped = int(active.sum()), 0
184 least_index = False
185 steps, limit = 0, _MAX_REPAIRS
186 while steps < limit:
187 steps += 1
188 step = _working_set_solve(cho, xu, C, b, active, m, reuse)
189 if step is None:
190 return None
191 x, lagrangian = step
193 slack = C.T @ x - b
194 following = _repair(active, lagrangian, slack, meq, _SET_TOL * scale, least_index)
196 if np.array_equal(following, active):
197 if not _certified(G, a, C, b, meq, x, lagrangian):
198 return None
199 return Attempt(x, xu, lagrangian, active, added, dropped)
201 key = following.tobytes()
202 if key in seen:
203 if least_index:
204 return None
205 # The block exchange is going round in circles. Drop to one index at
206 # a time, lowest first, and give it room to walk there.
207 least_index = True
208 limit = steps + m
209 following = _repair(active, lagrangian, slack, meq, _SET_TOL * scale, True)
210 key = following.tobytes()
211 # No second guard here: flipping one index always changes the set, and
212 # if that set has been seen before the check at the top catches it on
213 # the next pass, by which time `least_index` is set and it gives up.
215 seen.add(key)
216 added += int((following & ~active).sum())
217 dropped += int((~following & active).sum())
218 active = following
220 return None
223def _seed(
224 G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray, meq: int
225) -> tuple[tuple[np.ndarray, bool], np.ndarray, np.ndarray, float, np.ndarray] | None:
226 """Factorise ``G`` and pick the working set to start from, or decline.
228 The guess is the equalities plus whatever the unconstrained minimiser
229 violates, which is already the answer when it violates nothing.
231 Args:
232 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
233 a: ``(n,)`` vector of the linear term.
234 C: ``(n, m)`` constraint matrix.
235 b: ``(m,)`` right-hand side.
236 meq: Number of leading constraints held as equalities.
238 Returns:
239 The Cholesky factorisation, the unconstrained minimiser, the starting
240 working set, the scale the tolerances are measured against and
241 ``C.T @ xu``; or None if the problem is one the fast path does not take.
243 The last of those is returned rather than discarded because every repair
244 needs ``C_A^T xu`` for its right-hand side, and that is this vector indexed by
245 the working set. Rebuilding it from ``C_A`` cost a fancy-index copy of an
246 ``n`` by ``k`` block per repair -- 0.11 ms at ``n = 800`` -- to recompute
247 numbers already in hand.
248 """
249 n, m = C.shape
250 if n < _MIN_VARIABLES or m == 0 or meq > n:
251 return None
253 try:
254 cho = sla.cho_factor(G, check_finite=False)
255 xu = sla.cho_solve(cho, a, check_finite=False)
256 except (np.linalg.LinAlgError, ValueError):
257 return None
259 scale = max(1.0, float(np.abs(b).max(initial=0.0)))
260 ctxu = C.T @ xu
261 active = np.zeros(m, dtype=bool)
262 active[:meq] = True
263 active |= ctxu < b - _SET_TOL * scale
264 return cho, xu, active, scale, ctxu
267def _repair(
268 active: np.ndarray,
269 lagrangian: np.ndarray,
270 slack: np.ndarray,
271 meq: int,
272 tol: float,
273 least_index: bool,
274) -> np.ndarray:
275 """Return the working set to try next.
277 The block rule exchanges every index that violates its sign condition at
278 once, which is what converges in two to four repairs when it converges at
279 all. Exchanging a batch can also over-shoot -- a drop can remove support that
280 a later add restores, returning to a set already visited -- and that is what
281 the least-index rule is for: flip only the lowest-indexed offender, the
282 anti-cycling device of Bland and of Murty's least-index rule for
283 complementarity problems. It is slower per step and it is not a termination
284 proof here, since the general constraints leave no P-matrix to appeal to (see
285 :func:`_working_set_solve`), but it makes progress where the block rule
286 merely oscillates.
288 Args:
289 active: Current working set.
290 lagrangian: Multipliers at the current point.
291 slack: ``C.T @ x - b`` at the current point.
292 meq: Number of leading constraints held as equalities.
293 tol: Absolute tolerance for a sign being meant.
294 least_index: Whether to exchange one index rather than all of them.
296 Returns:
297 The next working set.
298 """
299 following = active.copy()
300 following[meq:] = ((lagrangian[meq:] > -tol) & active[meq:]) | (slack[meq:] < -tol)
301 if not least_index:
302 return following
304 offenders = np.flatnonzero(following != active)
305 following = active.copy()
306 if offenders.size:
307 first = int(offenders[0])
308 following[first] = not following[first]
309 return following
312def _working_set_solve(
313 cho: tuple[np.ndarray, bool],
314 xu: np.ndarray,
315 C: np.ndarray,
316 b: np.ndarray,
317 active: np.ndarray,
318 m: int,
319 reuse: _Reuse | None = None,
320) -> tuple[np.ndarray, np.ndarray] | None:
321 """Minimise with the working set held as equalities.
323 Stationarity gives ``x = xu + G^-1 C_A nu``, and substituting it into
324 ``C_A^T x = b_A`` leaves ``(C_A^T G^-1 C_A) nu = b_A - C_A^T xu``. That matrix
325 is positive definite exactly when ``C_A`` has full column rank, so its
326 Cholesky doubles as the rank test: a guess that made the working set linearly
327 dependent fails here instead of returning nonsense.
329 Args:
330 cho: Cholesky factorisation of ``G``, from ``scipy.linalg.cho_factor``.
331 xu: Unconstrained minimiser.
332 C: Constraint matrix.
333 b: Right-hand side.
334 active: Boolean mask of the working set.
335 m: Total number of constraints.
336 reuse: State carried across the repairs of one attempt, see
337 :class:`_Reuse`. None computes everything from ``C_A``, which is what
338 a caller outside :func:`attempt` gets.
340 Returns:
341 The minimiser and the full multiplier vector, or None if the working set
342 was rank deficient.
344 Which of two algebraically identical forms computes that is decided by
345 :data:`_SPLIT_MIN_WORK`; the larger one is the subject of :func:`_half_solve`.
347 Every scipy call here passes ``check_finite=False``, as
348 :func:`~cvx.quadprog._setup._factorize` does on the exact path and for the
349 reason given there: the reference implementation does not check either, and
350 scanning an ``n`` by ``k`` array on the way into every repair cost 7% of this
351 path at ``n = 800``. A non-finite entry that reaches here is not diagnosed,
352 and cannot produce a finite wrong answer, since the certificate is what
353 decides whether the result is returned at all.
354 """
355 lagrangian = np.zeros(m)
356 idx = np.flatnonzero(active)
357 if idx.size == 0:
358 return xu, lagrangian
360 if C.shape[0] * idx.size < _SPLIT_MIN_WORK:
361 # Small enough that the dispatches cost more than the flops they save.
362 CA = C[:, idx]
363 Y = sla.cho_solve(cho, CA, check_finite=False)
364 nu = _multipliers(CA.T @ Y, b[idx] - CA.T @ xu)
365 if nu is None:
366 return None
367 lagrangian[idx] = nu
368 return xu + Y @ nu, lagrangian
370 if reuse is None:
371 Z, rhs = _half_solve(cho, C[:, idx]), b[idx] - C[:, idx].T @ xu
372 else:
373 Z, rhs = _columns(cho, C, idx, reuse), b[idx] - reuse.ctxu[idx]
375 nu = _multipliers(_gram(Z), rhs)
376 if nu is None:
377 return None
378 lagrangian[idx] = nu
379 return xu + _half_solve_back(cho, Z @ nu), lagrangian
382def _columns(cho: tuple[np.ndarray, bool], C: np.ndarray, idx: np.ndarray, reuse: _Reuse) -> np.ndarray:
383 """Return ``U^-T C_A``, solving only for the columns not already held.
385 Args:
386 cho: Cholesky factorisation of ``G``.
387 C: ``(n, m)`` constraint matrix.
388 idx: 0-based indices of the working set.
389 reuse: State carried across the repairs of one attempt.
391 Returns:
392 The ``(n, k)`` block, gathered from the cache where it is available and
393 solved for and recorded where it is not. A ``reuse`` whose cache was
394 declined for its size re-solves the whole block, which is the same
395 arithmetic by a slower route.
396 """
397 if reuse.Z is None or reuse.have is None:
398 return _half_solve(cho, C[:, idx])
399 fresh = idx[~reuse.have[idx]]
400 if fresh.size:
401 reuse.Z[:, fresh] = _half_solve(cho, C[:, fresh])
402 reuse.have[fresh] = True
403 gathered: np.ndarray = reuse.Z[:, idx]
404 return gathered
407def _half_solve(cho: tuple[np.ndarray, bool], B: np.ndarray) -> np.ndarray:
408 """Return ``U^-T B``, one triangular solve rather than the two ``G^-1 B`` needs.
410 With ``G = U^T U`` the dual Hessian factors as
411 ``C_A^T G^-1 C_A = (U^-T C_A)^T (U^-T C_A)``, so the working-set system can be
412 formed from one half of the Cholesky factorisation instead of applying both
413 halves and then multiplying by ``C_A^T``. That halves the flops of the largest
414 term in a repair, and it is the same identity Section 3 of the accompanying
415 paper uses to relate ``R`` to the working set.
417 Args:
418 cho: Cholesky factorisation of ``G``, from ``scipy.linalg.cho_factor``.
419 B: ``(n, k)`` array, or an ``(n,)`` vector.
421 Returns:
422 ``U^-T B``, or ``L^-1 B`` when scipy handed back a lower factor.
423 """
424 factor, lower = cho
425 return sla.solve_triangular(factor, B, lower=lower, trans=0 if lower else 1, check_finite=False)
428def _half_solve_back(cho: tuple[np.ndarray, bool], v: np.ndarray) -> np.ndarray:
429 """Return ``U^-1 v``, the other half of the same factorisation.
431 Applied once per repair to a single vector, which recovers
432 ``G^-1 C_A nu = U^-1 (U^-T C_A nu)`` from the ``Z`` that
433 :func:`_half_solve` already produced, so the iterate costs one triangular
434 solve on a vector rather than an ``n`` by ``k`` product.
436 Args:
437 cho: Cholesky factorisation of ``G``.
438 v: ``(n,)`` vector.
440 Returns:
441 ``U^-1 v``, or ``L^-T v`` when scipy handed back a lower factor.
442 """
443 factor, lower = cho
444 return sla.solve_triangular(factor, v, lower=lower, trans=1 if lower else 0, check_finite=False)
447def _gram(Z: np.ndarray) -> np.ndarray:
448 """Return the upper triangle of ``Z^T Z``, at half the flops of the product.
450 The dual Hessian is symmetric, so a general product computes every off-diagonal
451 entry twice. ``syrk`` computes one triangle, and ``cho_factor(..., lower=False)``
452 reads only that triangle, so the other one is never needed.
454 Args:
455 Z: ``(n, k)`` array.
457 Returns:
458 ``(k, k)`` array whose upper triangle holds ``Z^T Z``.
459 """
460 result: np.ndarray = sla.blas.dsyrk(1.0, Z, trans=1, lower=0)
461 return result
464def _multipliers(H: np.ndarray, rhs: np.ndarray) -> np.ndarray | None:
465 """Solve the dual Hessian system, or decline a working set that is dependent.
467 ``H`` is positive definite exactly when the working set has full column rank,
468 so its Cholesky doubles as the rank test -- a dependent guess fails here
469 rather than returning something plausible. Only the upper triangle is read,
470 which is what lets :func:`_gram` fill one triangle and leave the other alone.
472 Args:
473 H: ``(k, k)`` dual Hessian, upper triangle significant.
474 rhs: ``(k,)`` right-hand side.
476 Returns:
477 The multipliers of the working-set constraints, or None if ``H`` was not
478 positive definite.
479 """
480 try:
481 return sla.cho_solve(sla.cho_factor(H, lower=False, check_finite=False), rhs, check_finite=False)
482 except (np.linalg.LinAlgError, ValueError):
483 return None
486def _certified(
487 G: np.ndarray,
488 a: np.ndarray,
489 C: np.ndarray,
490 b: np.ndarray,
491 meq: int,
492 x: np.ndarray,
493 lagrangian: np.ndarray,
494) -> bool:
495 """Return whether the KKT conditions hold, which for this problem is proof.
497 The program is strictly convex, so these conditions are sufficient and not
498 merely necessary: a point satisfying them is *the* minimiser. Stationarity is
499 checked against ``G`` directly rather than trusted from the construction,
500 since the construction is exactly what an ill-conditioned working set
501 corrupts.
503 Args:
504 G: Matrix of the quadratic term.
505 a: Vector of the linear term.
506 C: Constraint matrix.
507 b: Right-hand side.
508 meq: Number of leading constraints held as equalities.
509 x: Candidate minimiser.
510 lagrangian: Candidate multipliers.
512 Returns:
513 True when every condition holds to :data:`_CERTIFY_TOL`.
514 """
515 scale = max(
516 1.0,
517 float(np.abs(a).max(initial=0.0)),
518 float(np.abs(b).max(initial=0.0)),
519 float(np.abs(lagrangian).max(initial=0.0)),
520 )
521 tol = _CERTIFY_TOL * scale
522 slack = C.T @ x - b
523 return bool(
524 np.all(np.abs(G @ x - a - C @ lagrangian) <= tol)
525 and np.all(np.abs(slack[:meq]) <= tol)
526 and np.all(slack[meq:] >= -tol)
527 and np.all(lagrangian[meq:] >= -tol)
528 and np.all(np.abs(lagrangian[meq:] * slack[meq:]) <= tol)
529 )
532def _fast_solution(
533 G: np.ndarray,
534 a: np.ndarray,
535 C: np.ndarray,
536 b: np.ndarray,
537 meq: int,
538 check_finite: bool,
539) -> Solution | None:
540 """Assemble a :class:`Solution` from a certified fast-path attempt.
542 Anything malformed returns None rather than raising, so that the message the
543 caller sees for a bad problem is the one the exact path raises, unchanged.
545 Args:
546 G: ``(n, n)`` matrix of the quadratic term.
547 a: ``(n,)`` vector of the linear term.
548 C: ``(n, m)`` constraint matrix.
549 b: ``(m,)`` right-hand side.
550 meq: Number of leading constraints held as equalities.
551 check_finite: Whether to reject non-finite input.
553 Returns:
554 The solution, or None if the fast path declined the problem.
555 """
556 G = np.asarray(G, dtype=np.float64)
557 a = np.asarray(a, dtype=np.float64)
558 C = np.asarray(C, dtype=np.float64)
559 b = np.asarray(b, dtype=np.float64)
561 if meq < 0 or not _shapes_agree(G, a, C, b):
562 return None
563 if check_finite and not all(bool(np.isfinite(array).all()) for array in (G, a, C, b)):
564 return None
566 found = attempt(G, a, C, b, meq)
567 if found is None:
568 return None
570 return Solution(
571 x=found.x,
572 f=float(found.x @ G @ found.x) / 2.0 - float(a @ found.x),
573 xu=found.xu,
574 iterations=np.array([found.added, found.dropped], dtype=np.int64),
575 lagrangian=found.lagrangian,
576 iact=np.flatnonzero(found.active).astype(np.int64) + 1,
577 )
580def _shapes_agree(G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray) -> bool:
581 """Return whether the four arrays describe a well-formed program.
583 This is deliberately not the full validation :func:`_validate` performs. It
584 only has to be strict enough that the fast path never works on nonsense; a
585 problem it turns away is then rejected, with the proper message, by the exact
586 path that follows.
588 Args:
589 G: Matrix of the quadratic term.
590 a: Vector of the linear term.
591 C: Constraint matrix.
592 b: Right-hand side.
594 Returns:
595 True when the shapes are mutually consistent.
596 """
597 return (
598 G.ndim == 2
599 and G.shape[0] == G.shape[1]
600 and a.shape == (G.shape[0],)
601 and C.ndim == 2
602 and C.shape[0] == G.shape[0]
603 and b.shape == (C.shape[1],)
604 )