Coverage for src/cvx/quadprog/_pdas.py: 100%
108 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-09 18:50 +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
64class Attempt(NamedTuple):
65 """A candidate solution that has already passed the KKT certificate.
67 Attributes:
68 x: ``(n,)`` minimiser.
69 xu: ``(n,)`` unconstrained minimiser, ``G^-1 a``.
70 lagrangian: ``(m,)`` multipliers, zero off the active set.
71 active: ``(m,)`` boolean mask of the active constraints.
72 added: Constraints added to the working set, summed over all repairs.
73 dropped: Constraints removed from it, summed over all repairs.
74 """
76 x: np.ndarray
77 xu: np.ndarray
78 lagrangian: np.ndarray
79 active: np.ndarray
80 added: int
81 dropped: int
84def attempt(G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray, meq: int) -> Attempt | None:
85 """Try to solve by primal-dual active set, returning None if anything is off.
87 Every rejection path -- too small, singular, cycling, not converging, or
88 failing the certificate -- returns None rather than raising, because the
89 caller's response to all of them is the same: solve it the exact way.
91 Args:
92 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
93 a: ``(n,)`` vector of the linear term.
94 C: ``(n, m)`` constraint matrix, one column per constraint.
95 b: ``(m,)`` right-hand side of ``C.T @ x >= b``.
96 meq: Number of leading constraints held as equalities.
98 Returns:
99 A certified :class:`Attempt`, or None if the fast path did not produce
100 one.
101 """
102 seeded = _seed(G, a, C, b, meq)
103 if seeded is None:
104 return None
105 cho, xu, active, scale = seeded
106 m = C.shape[1]
108 seen: set[bytes] = set()
109 added, dropped = int(active.sum()), 0
110 least_index = False
111 steps, limit = 0, _MAX_REPAIRS
112 while steps < limit:
113 steps += 1
114 step = _working_set_solve(cho, xu, C, b, active, m)
115 if step is None:
116 return None
117 x, lagrangian = step
119 slack = C.T @ x - b
120 following = _repair(active, lagrangian, slack, meq, _SET_TOL * scale, least_index)
122 if np.array_equal(following, active):
123 if not _certified(G, a, C, b, meq, x, lagrangian):
124 return None
125 return Attempt(x, xu, lagrangian, active, added, dropped)
127 key = following.tobytes()
128 if key in seen:
129 if least_index:
130 return None
131 # The block exchange is going round in circles. Drop to one index at
132 # a time, lowest first, and give it room to walk there.
133 least_index = True
134 limit = steps + m
135 following = _repair(active, lagrangian, slack, meq, _SET_TOL * scale, True)
136 key = following.tobytes()
137 # No second guard here: flipping one index always changes the set, and
138 # if that set has been seen before the check at the top catches it on
139 # the next pass, by which time `least_index` is set and it gives up.
141 seen.add(key)
142 added += int((following & ~active).sum())
143 dropped += int((~following & active).sum())
144 active = following
146 return None
149def _seed(
150 G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray, meq: int
151) -> tuple[tuple[np.ndarray, bool], np.ndarray, np.ndarray, float] | None:
152 """Factorise ``G`` and pick the working set to start from, or decline.
154 The guess is the equalities plus whatever the unconstrained minimiser
155 violates, which is already the answer when it violates nothing.
157 Args:
158 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
159 a: ``(n,)`` vector of the linear term.
160 C: ``(n, m)`` constraint matrix.
161 b: ``(m,)`` right-hand side.
162 meq: Number of leading constraints held as equalities.
164 Returns:
165 The Cholesky factorisation, the unconstrained minimiser, the starting
166 working set and the scale the tolerances are measured against; or None if
167 the problem is one the fast path does not take.
168 """
169 n, m = C.shape
170 if n < _MIN_VARIABLES or m == 0 or meq > n:
171 return None
173 try:
174 cho = sla.cho_factor(G)
175 xu = sla.cho_solve(cho, a)
176 except (np.linalg.LinAlgError, ValueError):
177 return None
179 scale = max(1.0, float(np.abs(b).max(initial=0.0)))
180 active = np.zeros(m, dtype=bool)
181 active[:meq] = True
182 active |= C.T @ xu < b - _SET_TOL * scale
183 return cho, xu, active, scale
186def _repair(
187 active: np.ndarray,
188 lagrangian: np.ndarray,
189 slack: np.ndarray,
190 meq: int,
191 tol: float,
192 least_index: bool,
193) -> np.ndarray:
194 """Return the working set to try next.
196 The block rule exchanges every index that violates its sign condition at
197 once, which is what converges in two to four repairs when it converges at
198 all. Exchanging a batch can also over-shoot -- a drop can remove support that
199 a later add restores, returning to a set already visited -- and that is what
200 the least-index rule is for: flip only the lowest-indexed offender, the
201 anti-cycling device of Bland and of Murty's least-index rule for
202 complementarity problems. It is slower per step and it is not a termination
203 proof here, since the general constraints leave no P-matrix to appeal to (see
204 :func:`_working_set_solve`), but it makes progress where the block rule
205 merely oscillates.
207 Args:
208 active: Current working set.
209 lagrangian: Multipliers at the current point.
210 slack: ``C.T @ x - b`` at the current point.
211 meq: Number of leading constraints held as equalities.
212 tol: Absolute tolerance for a sign being meant.
213 least_index: Whether to exchange one index rather than all of them.
215 Returns:
216 The next working set.
217 """
218 following = active.copy()
219 following[meq:] = ((lagrangian[meq:] > -tol) & active[meq:]) | (slack[meq:] < -tol)
220 if not least_index:
221 return following
223 offenders = np.flatnonzero(following != active)
224 following = active.copy()
225 if offenders.size:
226 first = int(offenders[0])
227 following[first] = not following[first]
228 return following
231def _working_set_solve(
232 cho: tuple[np.ndarray, bool],
233 xu: np.ndarray,
234 C: np.ndarray,
235 b: np.ndarray,
236 active: np.ndarray,
237 m: int,
238) -> tuple[np.ndarray, np.ndarray] | None:
239 """Minimise with the working set held as equalities.
241 Stationarity gives ``x = xu + G^-1 C_A nu``, and substituting it into
242 ``C_A^T x = b_A`` leaves ``(C_A^T G^-1 C_A) nu = b_A - C_A^T xu``. That matrix
243 is positive definite exactly when ``C_A`` has full column rank, so its
244 Cholesky doubles as the rank test: a guess that made the working set linearly
245 dependent fails here instead of returning nonsense.
247 Args:
248 cho: Cholesky factorisation of ``G``, from ``scipy.linalg.cho_factor``.
249 xu: Unconstrained minimiser.
250 C: Constraint matrix.
251 b: Right-hand side.
252 active: Boolean mask of the working set.
253 m: Total number of constraints.
255 Returns:
256 The minimiser and the full multiplier vector, or None if the working set
257 was rank deficient.
258 """
259 lagrangian = np.zeros(m)
260 idx = np.flatnonzero(active)
261 if idx.size == 0:
262 return xu, lagrangian
264 CA = C[:, idx]
265 Y = sla.cho_solve(cho, CA)
266 try:
267 nu = sla.cho_solve(sla.cho_factor(CA.T @ Y), b[idx] - CA.T @ xu)
268 except (np.linalg.LinAlgError, ValueError):
269 return None
271 lagrangian[idx] = nu
272 return xu + Y @ nu, lagrangian
275def _certified(
276 G: np.ndarray,
277 a: np.ndarray,
278 C: np.ndarray,
279 b: np.ndarray,
280 meq: int,
281 x: np.ndarray,
282 lagrangian: np.ndarray,
283) -> bool:
284 """Return whether the KKT conditions hold, which for this problem is proof.
286 The program is strictly convex, so these conditions are sufficient and not
287 merely necessary: a point satisfying them is *the* minimiser. Stationarity is
288 checked against ``G`` directly rather than trusted from the construction,
289 since the construction is exactly what an ill-conditioned working set
290 corrupts.
292 Args:
293 G: Matrix of the quadratic term.
294 a: Vector of the linear term.
295 C: Constraint matrix.
296 b: Right-hand side.
297 meq: Number of leading constraints held as equalities.
298 x: Candidate minimiser.
299 lagrangian: Candidate multipliers.
301 Returns:
302 True when every condition holds to :data:`_CERTIFY_TOL`.
303 """
304 scale = max(
305 1.0,
306 float(np.abs(a).max(initial=0.0)),
307 float(np.abs(b).max(initial=0.0)),
308 float(np.abs(lagrangian).max(initial=0.0)),
309 )
310 tol = _CERTIFY_TOL * scale
311 slack = C.T @ x - b
312 return bool(
313 np.all(np.abs(G @ x - a - C @ lagrangian) <= tol)
314 and np.all(np.abs(slack[:meq]) <= tol)
315 and np.all(slack[meq:] >= -tol)
316 and np.all(lagrangian[meq:] >= -tol)
317 and np.all(np.abs(lagrangian[meq:] * slack[meq:]) <= tol)
318 )
321def _fast_solution(
322 G: np.ndarray,
323 a: np.ndarray,
324 C: np.ndarray,
325 b: np.ndarray,
326 meq: int,
327 check_finite: bool,
328) -> Solution | None:
329 """Assemble a :class:`Solution` from a certified fast-path attempt.
331 Anything malformed returns None rather than raising, so that the message the
332 caller sees for a bad problem is the one the exact path raises, unchanged.
334 Args:
335 G: ``(n, n)`` matrix of the quadratic term.
336 a: ``(n,)`` vector of the linear term.
337 C: ``(n, m)`` constraint matrix.
338 b: ``(m,)`` right-hand side.
339 meq: Number of leading constraints held as equalities.
340 check_finite: Whether to reject non-finite input.
342 Returns:
343 The solution, or None if the fast path declined the problem.
344 """
345 G = np.asarray(G, dtype=np.float64)
346 a = np.asarray(a, dtype=np.float64)
347 C = np.asarray(C, dtype=np.float64)
348 b = np.asarray(b, dtype=np.float64)
350 if meq < 0 or not _shapes_agree(G, a, C, b):
351 return None
352 if check_finite and not all(bool(np.isfinite(array).all()) for array in (G, a, C, b)):
353 return None
355 found = attempt(G, a, C, b, meq)
356 if found is None:
357 return None
359 return Solution(
360 x=found.x,
361 f=float(found.x @ G @ found.x) / 2.0 - float(a @ found.x),
362 xu=found.xu,
363 iterations=np.array([found.added, found.dropped], dtype=np.int64),
364 lagrangian=found.lagrangian,
365 iact=np.flatnonzero(found.active).astype(np.int64) + 1,
366 )
369def _shapes_agree(G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray) -> bool:
370 """Return whether the four arrays describe a well-formed program.
372 This is deliberately not the full validation :func:`_validate` performs. It
373 only has to be strict enough that the fast path never works on nonsense; a
374 problem it turns away is then rejected, with the proper message, by the exact
375 path that follows.
377 Args:
378 G: Matrix of the quadratic term.
379 a: Vector of the linear term.
380 C: Constraint matrix.
381 b: Right-hand side.
383 Returns:
384 True when the shapes are mutually consistent.
385 """
386 return (
387 G.ndim == 2
388 and G.shape[0] == G.shape[1]
389 and a.shape == (G.shape[0],)
390 and C.ndim == 2
391 and C.shape[0] == G.shape[0]
392 and b.shape == (C.shape[1],)
393 )