Coverage for src/cvx/quadprog/_sweep.py: 100%
87 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-09 18:50 +0000
1r"""Reuse one factorisation across a family of QPs that differ only in ``a``.
3A frontier sweep, a rolling rebalance and a scenario grid all solve the same
4problem repeatedly with a slightly different linear term. Solved independently,
5each one rediscovers an active set it almost always already had: a 1% relative
6perturbation of ``a`` moves 2.4 of 167 active constraints on a box problem and
7none at all on a budget-plus-bounds problem, whose long-only optimum is a vertex
8with under 1% of the variables interior.
10The saving is not in passing the active set back in. Installing a set of size
11``k`` costs ``k`` Householder insertions, ``O(n^2 k)``, which is what the cold
12walk already pays for its own insertions. It is that ``J`` depends only on ``G``
13and ``R`` only on ``G`` and the active set, so across such a family both are
14reusable verbatim and recovering the solution costs ``O(nk)``:
16.. math::
17 x_u = J J^T a, \\quad r = b_A - C_A^T x_u, \\quad R^T y = r, \\quad
18 x = x_u + J_{:,:k}\\, y, \\quad R \\lambda = y
20That point is the answer exactly when the KKT conditions hold -- every multiplier
21on an inequality non-negative, no inactive constraint violated -- which is checked
22rather than assumed. When the check fails the active set is *repaired* rather than
23abandoned: multipliers that have gone negative mark constraints that no longer
24belong, dropping them restores the dual feasibility the iteration requires, and it
25resumes from there instead of from the unconstrained minimum. Only if everything
26is dropped does that amount to a cold solve. Either way a :class:`Sweep` never
27returns a different answer from :func:`~cvx.quadprog.solve_qp`; it is only faster.
29Why a class rather than a ``warm_start=`` argument to ``solve_qp``: the cached
30factors are valid only for the ``G`` and ``C`` they were built from, and a function
31cannot check that a caller passed the same ones without an ``O(n^2)`` comparison
32that would cost more than it saves. Owning the data makes the mismatch
33unrepresentable.
34"""
36# G, C, R, J and A are the names from Goldfarb & Idnani (1983), as everywhere else
37# in this package; lowercasing them would obscure the correspondence to the paper.
38# ruff: noqa: N803, N806
40from typing import NamedTuple
42import numpy as np
43from scipy.linalg.blas import dtpsv
45from ._base import _EMPTY, VSMALL, Solution, _WarmEntry
46from ._setup import _factorize, _validate
47from ._solve import _solve_with_factors
48from ._steps import _drop_constraint
49from ._structure import _default_constraints
51__all__ = ["Sweep"]
53# How far a multiplier or a slack must be on the wrong side of zero before the
54# cached active set is judged stale. Measured over 480 sweep steps: where the
55# cached set was still optimal the worst multiplier was +3.9e-4 and the worst
56# slack +1.6e-4; where it was not, they reached -6.0e-2 and -3.1e-3. Four orders
57# of magnitude of separation, so this threshold decides nothing delicate.
58#
59# It is also one-sided in the safe direction. Too strict merely falls back to a
60# full solve; only too loose could return a non-optimal point, and "too loose"
61# here would mean crossing four orders of magnitude.
62_STALE_MARGIN = 32.0
65class _Cache(NamedTuple):
66 """The factorisation a previous solve ended on.
68 Held as one object so that a single ``is None`` test narrows all three for the
69 type checker, and so that they cannot get out of step with one another.
71 Attributes:
72 J: Inverse Cholesky factor as the iteration left it.
73 R: Packed triangular factor of the active constraint normals.
74 iact: 1-based active set the factors correspond to.
75 """
77 J: np.ndarray
78 R: np.ndarray
79 iact: np.ndarray
82class Sweep:
83 """Solve a family of QPs sharing ``G``, ``C``, ``b`` and ``meq``.
85 Only the linear term changes between calls. The first call solves from
86 scratch; later ones reuse the factorisation when the active set still holds.
88 >>> import numpy as np
89 >>> from cvx.quadprog import Sweep, solve_qp
90 >>> G = np.eye(3)
91 >>> C = np.array([[-4.0, 2.0, 0.0], [-3.0, 1.0, -2.0], [0.0, 0.0, 1.0]])
92 >>> b = np.array([-8.0, 2.0, 0.0])
93 >>> sweep = Sweep(G, C, b)
94 >>> a = np.array([0.0, 5.0, 0.0])
95 >>> bool(np.allclose(sweep.solve(a).x, solve_qp(G, a, C, b).x))
96 True
97 >>> bool(np.allclose(sweep.solve(1.01 * a).x, solve_qp(G, 1.01 * a, C, b).x))
98 True
99 """
101 def __init__(
102 self,
103 G: np.ndarray,
104 C: np.ndarray | None = None,
105 b: np.ndarray | None = None,
106 meq: int = 0,
107 check_finite: bool = False,
108 ) -> None:
109 """Fix the part of the problem that does not vary.
111 Args:
112 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
113 C: ``(n, m)`` constraint matrix, one column per constraint. Defaults to
114 the unconstrained problem.
115 b: ``(m,)`` right-hand side of the constraints.
116 meq: Number of leading constraints to treat as equalities.
117 check_finite: Whether to reject NaN and infinity in ``G``, ``C`` and
118 ``b``, and in each ``a`` passed to :meth:`solve`. Off by default,
119 matching :func:`~cvx.quadprog.solve_qp`.
121 Raises:
122 ValueError: If the shapes are inconsistent, if ``meq`` is out of range,
123 or if ``G`` is not positive definite.
124 """
125 G = np.asarray(G, dtype=np.float64)
126 self.C, self.b, self.meq = _default_constraints(G, C, b, meq)
127 self.n, self._q = _validate(G, np.zeros(len(G)), self.C, self.b, self.meq, check_finite)
128 self._check_finite = check_finite
129 self.G = G
131 # The Cholesky is a property of G alone, so it is done once here and every
132 # later cold solve is handed the factor instead, via `factorized=True`.
133 # That is the same reuse the reference package offers, and it is the part
134 # of the saving that applies even when the active set does change.
135 self._Rinv, _xu = _factorize(G, np.zeros(self.n), False)
136 self._cache: _Cache | None = None
137 self.hits = 0
138 self.misses = 0
140 def solve(self, a: np.ndarray) -> Solution:
141 """Solve for a new linear term.
143 Args:
144 a: ``(n,)`` vector of the linear term.
146 Returns:
147 The same :class:`~cvx.quadprog.Solution` that
148 :func:`~cvx.quadprog.solve_qp` would return for this problem, except
149 that ``iterations`` is ``(0, 0)`` when the cached factorisation was
150 reused outright -- no active-set iteration was performed.
152 Raises:
153 ValueError: If ``a`` has the wrong shape, if the constraints admit no
154 solution, or if ``check_finite`` is set and ``a`` holds a
155 non-finite value.
156 """
157 a = np.asarray(a, dtype=np.float64)
158 warm = None
159 cache = self._cache
160 if cache is not None and self._usable(a):
161 hit = self._reuse(a, cache)
162 if hit is not None:
163 self.hits += 1
164 return hit
165 # The cached set is stale, but it is still a far better place to start
166 # than the unconstrained minimum: repairing it into a dual-feasible
167 # state costs a few drops, where a cold solve re-walks the whole set.
168 warm = self._repair(a, cache)
170 self.misses += 1
171 solution, J, R = _solve_with_factors(self._Rinv, a, self.C, self.b, self.meq, True, self._check_finite, warm)
172 self._cache = _Cache(J, R, solution.iact)
173 return solution
175 def _usable(self, a: np.ndarray) -> bool:
176 """Return whether the cache may be consulted at all for this ``a``.
178 Args:
179 a: ``(n,)`` vector of the linear term.
181 Returns:
182 False when ``a`` is the wrong length, or when ``check_finite`` is set
183 and it is not finite -- every KKT comparison against NaN is False, so
184 without this the fast path would *accept* a non-finite point instead
185 of rejecting it. Falling back lets the full solve raise, which is what
186 the caller asked for.
187 """
188 if len(a) != self.n:
189 return False
190 return not (self._check_finite and not np.isfinite(a).all())
192 def _recover(
193 self, a: np.ndarray, J: np.ndarray, R: np.ndarray, iact: np.ndarray, nact: int
194 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
195 """Return the minimiser over a given active set, and its multipliers.
197 Costs ``O(nk)``: the factors already encode everything about ``G`` and the
198 active constraints, so only the right-hand side has changed.
200 Args:
201 a: ``(n,)`` linear term.
202 J: Inverse Cholesky factor for this active set.
203 R: Packed triangular factor for this active set.
204 iact: 1-based active set, first ``nact`` entries valid.
205 nact: Size of the active set.
207 Returns:
208 ``(x, lam, xu)`` -- the minimiser subject to the active set held as
209 equalities, its multipliers, and the unconstrained minimiser.
210 """
211 xu = J @ (J.T @ a)
212 if nact == 0:
213 # Distinct arrays even though the values coincide: a resumed solve
214 # updates the iterate in place, and would otherwise corrupt ``xu``
215 # along with it. The cold path copies here for the same reason.
216 return xu.copy(), _EMPTY, xu
217 active = iact[:nact] - 1
218 y = dtpsv(nact, R, self.b[active] - self.C[:, active].T @ xu, lower=0, trans=1, overwrite_x=True)
219 x = xu + J[:, :nact] @ y
220 lam = dtpsv(nact, R, y.copy(), lower=0, trans=0, overwrite_x=True)
221 return x, lam, xu
223 def _reuse(self, a: np.ndarray, cache: "_Cache") -> Solution | None:
224 """Return the solution from the cached factorisation, or None if it is stale.
226 Args:
227 a: ``(n,)`` vector of the linear term.
228 cache: The factorisation a previous solve ended on.
230 Returns:
231 A :class:`~cvx.quadprog.Solution` when the cached active set still
232 satisfies the KKT conditions for this ``a``, otherwise None.
233 """
234 J, R, iact = cache
235 x, lam, xu = self._recover(a, J, R, iact, len(iact))
236 lagr = np.zeros(self._q)
237 lagr[iact - 1] = lam
238 return self._verified(a, x, xu, lagr, lam, iact)
240 def _repair(self, a: np.ndarray, cache: "_Cache") -> _WarmEntry:
241 """Turn a stale active set into a dual-feasible state to resume from.
243 A multiplier that has gone negative marks a constraint that no longer
244 belongs in the active set. Dropping it and recomputing is exactly the
245 step the solver's own inner loop takes, and repeating until none is
246 negative restores the invariant the iteration requires. Whatever is left
247 may still be primally infeasible -- constraints outside the set may be
248 violated -- and driving that to zero is what the resumed loop is for.
250 Terminates because each pass either stops or shrinks the active set; in
251 the worst case everything is dropped and the resumed loop starts from the
252 unconstrained minimum, which is the cold start.
254 Args:
255 a: ``(n,)`` linear term.
256 cache: The stale factorisation.
258 Returns:
259 A :class:`~cvx.quadprog._solve._WarmEntry` satisfying that invariant.
260 """
261 # Copied because a repair mutates them, and the cache must survive intact
262 # if the resumed solve then fails.
263 J, R = cache.J.copy(), cache.R.copy()
264 nact = len(cache.iact)
265 iact = np.zeros(self._q, dtype=np.int64)
266 iact[:nact] = cache.iact
267 uv = np.zeros(min(self.n, self._q))
269 while True:
270 x, lam, xu = self._recover(a, J, R, iact, nact)
271 uv[:nact] = lam
272 if nact == 0:
273 break
274 # Equalities carry unrestricted multipliers, so only inequalities can
275 # mark themselves as no longer belonging.
276 candidates = np.where(iact[:nact] > self.meq, lam, np.inf)
277 worst = int(np.argmin(candidates))
278 if candidates[worst] >= 0.0:
279 break
280 nact = _drop_constraint(worst + 1, nact, uv, iact, J, R)
282 obj = 0.5 * float(x @ (self.G @ x)) - float(a @ x)
283 return _WarmEntry(J, R, iact, nact, x, uv, obj, xu)
285 def _verified(
286 self,
287 a: np.ndarray,
288 x: np.ndarray,
289 xu: np.ndarray,
290 lagr: np.ndarray,
291 lam: np.ndarray,
292 iact: np.ndarray,
293 ) -> Solution | None:
294 """Return a Solution if ``x`` satisfies the KKT conditions, else None.
296 For a strictly convex QP the KKT conditions are sufficient, so this is a
297 proof rather than a heuristic: dual feasibility on the inequalities, and
298 primal feasibility of everything not held active.
300 Args:
301 a: ``(n,)`` linear term.
302 x: Candidate minimiser.
303 xu: Unconstrained minimiser.
304 lagr: Full-length multiplier vector.
305 lam: Multipliers of the active constraints only.
306 iact: 1-based active set.
308 Returns:
309 The :class:`~cvx.quadprog.Solution`, or None if the cache is stale.
310 """
311 scale = _STALE_MARGIN * VSMALL * max(1.0, float(np.max(np.abs(x))))
313 if lam.size and np.any(lam[iact > self.meq] < -scale):
314 return None
316 sv = self.C.T @ x - self.b
317 if iact.size:
318 sv[iact - 1] = 0.0
319 if np.any(sv[self.meq :] < -scale) or np.any(np.abs(sv[: self.meq]) > scale):
320 return None
322 obj = 0.5 * float(x @ (self.G @ x)) - float(a @ x)
323 return Solution(x, obj, xu, np.zeros(2, dtype=np.int64), lagr, iact)