Coverage for src/cvx/quadprog/

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 needs no factorisation at all:

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 costs ``O(n^2)``, and the ``n^2`` is all in the first step: ``J`` is dense

21once the first insertion has touched it, so ``x_u = J J^T a`` is two full

22matrix-vector products and does not shrink with the active set. Everything after

23it is ``O(nk + k^2)`` -- one product or gather for ``C_A^T x_u``, two triangular

24solves of order ``k``, one product with ``J_{:,:k}``. Against the ``O(n^2 k)``

25that installing the same set from scratch would cost, the saving is a factor of

26order ``k``.

28Verifying the answer is separate, since the KKT check has to look at every

29constraint and not only the active ones. On a bound-constrained family both that

30check and the ``C_A^T x_u`` above are gathers rather than products -- see

31:func:`~cvx.quadprog._structure._slack_evaluator` and :meth:`Sweep._active_product`

32-- so the whole hit costs ``O(n^2 + m)``. On a dense ``C`` the verification adds

33``O(nm)``, which dominates once ``m`` exceeds ``n``.

35That point is the answer exactly when the KKT conditions hold -- every multiplier

36on an inequality non-negative, no inactive constraint violated -- which is checked

37rather than assumed. When the check fails the active set is *repaired* rather than

38abandoned: multipliers that have gone negative mark constraints that no longer

39belong, dropping them restores the dual feasibility the iteration requires, and it

40resumes from there instead of from the unconstrained minimum. Only if everything

41is dropped does that amount to a cold solve. Either way a :class:`Sweep` never

42returns a different answer from :func:`~cvx.quadprog.solve_qp`; it is only faster.

44Why a class rather than a ``warm_start=`` argument to ``solve_qp``: the cached

45factors are valid only for the ``G`` and ``C`` they were built from, and a function

46cannot check that a caller passed the same ones without an ``O(n^2)`` comparison

47that would cost more than it saves. Owning the data makes the mismatch

48unrepresentable.

49"""

51# G, C, R, J and A are the names from Goldfarb & Idnani (1983), as everywhere else

52# in this package; lowercasing them would obscure the correspondence to the paper.

53# ruff: noqa: N803, N806

55from typing import NamedTuple

57import numpy as np

58from scipy.linalg.blas import dtpsv

60from . import _threads

61from ._base import _EMPTY, VSMALL, Solution, _WarmEntry

62from ._setup import _factorize, _validate

63from ._solve import _solve_with_factors

64from ._steps import _drop_constraint

65from ._structure import _analyse_constraints, _default_constraints, _slack_evaluator

67__all__ = ["Sweep"]

69# How far a multiplier or a slack must be on the wrong side of zero before the

70# cached active set is judged stale. Measured over 480 sweep steps: where the

71# cached set was still optimal the worst multiplier was +3.9e-4 and the worst

72# slack +1.6e-4; where it was not, they reached -6.0e-2 and -3.1e-3. Four orders

73# of magnitude of separation, so this threshold decides nothing delicate.

74#

75# It is also one-sided in the safe direction. Too strict merely falls back to a

76# full solve; only too loose could return a non-optimal point, and "too loose"

77# here would mean crossing four orders of magnitude.

78_STALE_MARGIN = 32.0

81class _Cache(NamedTuple):

82 """The factorisation a previous solve ended on.

84 Held as one object so that a single ``is None`` test narrows all three for the

85 type checker, and so that they cannot get out of step with one another.

87 Attributes:

88 J: Inverse Cholesky factor as the iteration left it.

89 R: Packed triangular factor of the active constraint normals.

90 iact: 1-based active set the factors correspond to.

91 """

93 J: np.ndarray

94 R: np.ndarray

95 iact: np.ndarray

98class Sweep:

99 """Solve a family of QPs sharing ``G``, ``C``, ``b`` and ``meq``.

100

101 Only the linear term changes between calls. The first call solves from

102 scratch; later ones reuse the factorisation when the active set still holds.

103

104 >>> import numpy as np

105 >>> from cvx.quadprog import Sweep, solve_qp

106 >>> G = np.eye(3)

107 >>> C = np.array([[-4.0, 2.0, 0.0], [-3.0, 1.0, -2.0], [0.0, 0.0, 1.0]])

108 >>> b = np.array([-8.0, 2.0, 0.0])

109 >>> sweep = Sweep(G, C, b)

110 >>> a = np.array([0.0, 5.0, 0.0])

111 >>> bool(np.allclose(sweep.solve(a).x, solve_qp(G, a, C, b).x))

112 True

113 >>> bool(np.allclose(sweep.solve(1.01 * a).x, solve_qp(G, 1.01 * a, C, b).x))

114 True

115 """

116

117 def __init__(

118 self,

119 G: np.ndarray,

120 C: np.ndarray | None = None,

121 b: np.ndarray | None = None,

122 meq: int = 0,

123 check_finite: bool = False,

124 blas_threads: int | None = None,

125 ) -> None:

126 """Fix the part of the problem that does not vary.

127

128 Args:

129 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.

130 C: ``(n, m)`` constraint matrix, one column per constraint. Defaults to

131 the unconstrained problem.

132 b: ``(m,)`` right-hand side of the constraints.

133 meq: Number of leading constraints to treat as equalities.

134 check_finite: Whether to reject NaN and infinity in ``G``, ``C`` and

135 ``b``, and in each ``a`` passed to :meth:`solve`. Off by default,

136 matching :func:`~cvx.quadprog.solve_qp`.

137 blas_threads: Cap the BLAS thread count for the expensive parts of this

138 sweep, as :func:`~cvx.quadprog.solve_qp`'s argument of the same name

139 does for one solve: the factorisation below, and every

140 :meth:`solve` that misses the cache. A hit is deliberately left

141 outside the context, which costs ~100 microseconds against an

142 ``O(n^2)`` recovery -- and against the handful of microseconds a

143 hit actually takes at the small ``n`` this class is most worth

144 using at, where that arithmetic is still below the dispatch

145 overhead.

146

147 Decided once here rather than per call, because ``n`` is fixed for

148 this object's lifetime and so the automatic gate's answer is too.

149 Left unset, that gate is consulted exactly as it is for

150 ``solve_qp`` -- see there for what it does and does not change, and

151 :func:`~cvx.quadprog._threads.auto_cap_threads` for the conditions.

152

153 Raises:

154 ValueError: If the shapes are inconsistent, if ``meq`` is out of range,

155 if ``G`` is not positive definite, or if ``blas_threads`` is not at

156 least 1.

157 ImportError: If ``blas_threads`` is given and ``threadpoolctl`` is not

158 installed.

159 """

160 G = np.asarray(G, dtype=np.float64)

161 self.C, self.b, self.meq = _default_constraints(G, C, b, meq)

162 self.n, self._q = _validate(G, np.zeros(len(G)), self.C, self.b, self.meq, check_finite)

163 self._check_finite = check_finite

164 self.G = G

165

166 # C, b and meq are fixed for this object's lifetime, so the shape analysis

167 # runs once here and is amortised over every call -- where solve_qp has to

168 # pay it per solve. Before this, the hit path re-derived the slacks with a

169 # dense `C.T @ x` and never reached the bound-constraint gather at all,

170 # which on a box family is 13% of a hit at n = 800 (#109).

171 self._single, self._srow, self._sval = _analyse_constraints(self.C)

172 self._slack_of = _slack_evaluator(self.C, self._single, self._srow, self._sval)

173

174 # An explicit count is used as given; otherwise the automatic gate decides,

175 # and it is asked once because `n` cannot change under it. None means "leave

176 # the process alone", which is what `scoped_limit` turns into a no-op.

177 #

178 # Sweep is the API most exposed to the OpenBLAS collapse -- large problems,

179 # solved repeatedly -- and until #107 it was the one path with no guard,

180 # because it calls `_solve_with_factors` below the level solve_qp installs

181 # the context at.

182 self._blas_threads = blas_threads if blas_threads is not None else _threads.auto_cap_threads(self.n, fast=False)

183

184 # The Cholesky is a property of G alone, so it is done once here and every

185 # later cold solve is handed the factor instead, via `factorized=True`.

186 # That is the same reuse the reference package offers, and it is the part

187 # of the saving that applies even when the active set does change.

188 #

189 # It is also the single largest BLAS call this object ever makes, at

190 # O(n^3), so it is inside the cap. A bad `blas_threads` therefore raises

191 # here, at construction, rather than at the first solve.

192 with _threads.scoped_limit(self._blas_threads):

193 self._Rinv, _xu = _factorize(G, np.zeros(self.n), False)

194 self._cache: _Cache | None = None

195 self.hits = 0

196 self.misses = 0

197

198 def solve(self, a: np.ndarray) -> Solution:

199 """Solve for a new linear term.

200

201 Args:

202 a: ``(n,)`` vector of the linear term.

203

204 Returns:

205 The same :class:`~cvx.quadprog.Solution` that

206 :func:`~cvx.quadprog.solve_qp` would return for this problem, except

207 that ``iterations`` is ``(0, 0)`` when the cached factorisation was

208 reused outright -- no active-set iteration was performed.

209

210 Raises:

211 ValueError: If ``a`` has the wrong shape, if the constraints admit no

212 solution, or if ``check_finite`` is set and ``a`` holds a

213 non-finite value.

214 """

215 a = np.asarray(a, dtype=np.float64)

216 warm = None

217 cache = self._cache

218 if cache is not None and self._usable(a):

219 hit = self._reuse(a, cache)

220 if hit is not None:

221 self.hits += 1

222 return hit

223 # The cached set is stale, but it is still a far better place to start

224 # than the unconstrained minimum: repairing it into a dual-feasible

225 # state costs a few drops, where a cold solve re-walks the whole set.

226 warm = self._repair(a, cache)

227

228 self.misses += 1

229 # Only the miss is wrapped. A hit is an O(n^2) recovery plus a KKT check,

230 # and entering a threadpoolctl context costs ~100 microseconds, which would

231 # be a tax on exactly the path this class exists to make cheap.

232 with _threads.scoped_limit(self._blas_threads):

233 solution, J, R = _solve_with_factors(

234 self._Rinv, a, self.C, self.b, self.meq, True, self._check_finite, warm

235 )

236 self._cache = _Cache(J, R, solution.iact)

237 return solution

238

239 def _usable(self, a: np.ndarray) -> bool:

240 """Return whether the cache may be consulted at all for this ``a``.

241

242 Args:

243 a: ``(n,)`` vector of the linear term.

244

245 Returns:

246 False when ``a`` is the wrong length, or when ``check_finite`` is set

247 and it is not finite -- every KKT comparison against NaN is False, so

248 without this the fast path would *accept* a non-finite point instead

249 of rejecting it. Falling back lets the full solve raise, which is what

250 the caller asked for.

251 """

252 if len(a) != self.n:

253 return False

254 return not (self._check_finite and not np.isfinite(a).all())

255

256 def _recover(

257 self, a: np.ndarray, J: np.ndarray, R: np.ndarray, iact: np.ndarray, nact: int

258 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:

259 """Return the minimiser over a given active set, and its multipliers.

260

261 Costs ``O(n^2)``, all of it in ``x_u = J J^T a`` below: the factors already

262 encode everything about ``G`` and the active constraints, so only the

263 right-hand side has changed, but ``J`` is a dense ``(n, n)`` and those two

264 products do not shrink with the active set. The rest is ``O(nk + k^2)``.

265

266 Args:

267 a: ``(n,)`` linear term.

268 J: Inverse Cholesky factor for this active set.

269 R: Packed triangular factor for this active set.

270 iact: 1-based active set, first ``nact`` entries valid.

271 nact: Size of the active set.

272

273 Returns:

274 ``(x, lam, xu)`` -- the minimiser subject to the active set held as

275 equalities, its multipliers, and the unconstrained minimiser.

276 """

277 xu = J @ (J.T @ a)

278 if nact == 0:

279 # Distinct arrays even though the values coincide: a resumed solve

280 # updates the iterate in place, and would otherwise corrupt ``xu``

281 # along with it. The cold path copies here for the same reason.

282 return xu.copy(), _EMPTY, xu

283 active = iact[:nact] - 1

284 y = dtpsv(nact, R, self.b[active] - self._active_product(active, xu), lower=0, trans=1, overwrite_x=True)

285 x = xu + J[:, :nact] @ y

286 lam = dtpsv(nact, R, y.copy(), lower=0, trans=0, overwrite_x=True)

287 return x, lam, xu

288

289 def _active_product(self, active: np.ndarray, xu: np.ndarray) -> np.ndarray:

290 """Return ``C_A^T xu`` for the active columns, by gather where it can.

291

292 Where every active column holds a single nonzero the product is ``k``

293 multiplications (#109). The test is on the *active* columns rather than on

294 all of ``C``, so a mixed matrix -- a budget row plus bounds -- still takes

295 that path whenever the set happens to be all bounds. It is ``O(k)``

296 against what it guards.

297

298 What it guards is no longer a block of ``C``. Fancy-indexing an ``(n, k)``

299 block out and multiplying against it costs ``O(nk)`` in flops but a copy

300 of the block in bandwidth, and the copy is what dominates: at

301 ``n = 800``, ``m = 400``, ``k = 50`` it measured 0.024 ms against 0.005 ms

302 for evaluating all ``m`` products and keeping ``k`` of them, a product the

303 evaluator of :mod:`._structure` has already chosen the cheapest form for.

304 Doing the arithmetic for constraints whose answers are then discarded is

305 the faster route by a factor of five, and on a reused solve of a

306 dense-``C`` family it was 54% of the whole cost.

307

308 Args:

309 active: 0-based indices of the active constraints.

310 xu: ``(n,)`` unconstrained minimiser.

311

312 Returns:

313 The length-``k`` vector of active constraint values at ``xu``.

314 """

315 # Annotated on the way out for the reason given in _threads.limit: indexing

316 # an ndarray by an ndarray is typed as Any, so returning either expression

317 # directly is an untyped escape under --strict.

318 if self._single[active].all():

319 gathered: np.ndarray = self._sval[active] * xu[self._srow[active]]

320 return gathered

321 product: np.ndarray = self._slack_of(xu)[active]

322 return product

323

324 def _reuse(self, a: np.ndarray, cache: "_Cache") -> Solution | None:

325 """Return the solution from the cached factorisation, or None if it is stale.

326

327 Args:

328 a: ``(n,)`` vector of the linear term.

329 cache: The factorisation a previous solve ended on.

330

331 Returns:

332 A :class:`~cvx.quadprog.Solution` when the cached active set still

333 satisfies the KKT conditions for this ``a``, otherwise None.

334 """

335 J, R, iact = cache

336 x, lam, xu = self._recover(a, J, R, iact, len(iact))

337 lagr = np.zeros(self._q)

338 lagr[iact - 1] = lam

339 return self._verified(a, x, xu, lagr, lam, iact)

340

341 def _repair(self, a: np.ndarray, cache: "_Cache") -> _WarmEntry:

342 """Turn a stale active set into a dual-feasible state to resume from.

343

344 A multiplier that has gone negative marks a constraint that no longer

345 belongs in the active set. Dropping it and recomputing is exactly the

346 step the solver's own inner loop takes, and repeating until none is

347 negative restores the invariant the iteration requires. Whatever is left

348 may still be primally infeasible -- constraints outside the set may be

349 violated -- and driving that to zero is what the resumed loop is for.

350

351 Terminates because each pass either stops or shrinks the active set; in

352 the worst case everything is dropped and the resumed loop starts from the

353 unconstrained minimum, which is the cold start.

354

355 Args:

356 a: ``(n,)`` linear term.

357 cache: The stale factorisation.

358

359 Returns:

360 A :class:`~cvx.quadprog._solve._WarmEntry` satisfying that invariant.

361 """

362 # Copied because a repair mutates them, and the cache must survive intact

363 # if the resumed solve then fails.

364 J, R = cache.J.copy(), cache.R.copy()

365 nact = len(cache.iact)

366 iact = np.zeros(self._q, dtype=np.int64)

367 iact[:nact] = cache.iact

368 uv = np.zeros(min(self.n, self._q))

369

370 while True:

371 x, lam, xu = self._recover(a, J, R, iact, nact)

372 uv[:nact] = lam

373 if nact == 0:

374 break

375 # Equalities carry unrestricted multipliers, so only inequalities can

376 # mark themselves as no longer belonging.

377 candidates = np.where(iact[:nact] > self.meq, lam, np.inf)

378 worst = int(np.argmin(candidates))

379 if candidates[worst] >= 0.0:

380 break

381 nact = _drop_constraint(worst + 1, nact, uv, iact, J, R)

382

383 obj = 0.5 * float(x @ (self.G @ x)) - float(a @ x)

384 return _WarmEntry(J, R, iact, nact, x, uv, obj, xu)

385

386 def _verified(

387 self,

388 a: np.ndarray,

389 x: np.ndarray,

390 xu: np.ndarray,

391 lagr: np.ndarray,

392 lam: np.ndarray,

393 iact: np.ndarray,

394 ) -> Solution | None:

395 """Return a Solution if ``x`` satisfies the KKT conditions, else None.

396

397 For a strictly convex QP the KKT conditions are sufficient, so this is a

398 proof rather than a heuristic: dual feasibility on the inequalities, and

399 primal feasibility of everything not held active.

400

401 Args:

402 a: ``(n,)`` linear term.

403 x: Candidate minimiser.

404 xu: Unconstrained minimiser.

405 lagr: Full-length multiplier vector.

406 lam: Multipliers of the active constraints only.

407 iact: 1-based active set.

408

409 Returns:

410 The :class:`~cvx.quadprog.Solution`, or None if the cache is stale.

411 """

412 scale = _STALE_MARGIN * VSMALL * max(1.0, float(np.max(np.abs(x))))

413

414 if lam.size and np.any(lam[iact > self.meq] < -scale):

415 return None

416

417 # Fresh array from every branch of the evaluator, which matters because the

418 # active entries are forced to zero in place on the next line.

419 sv = self._slack_of(x) - self.b

420 if iact.size:

421 sv[iact - 1] = 0.0

422 if np.any(sv[self.meq :] < -scale) or np.any(np.abs(sv[: self.meq]) > scale):

423 return None

424

425 obj = 0.5 * float(x @ (self.G @ x)) - float(a @ x)

426 return Solution(x, obj, xu, np.zeros(2, dtype=np.int64), lagr, iact)

Coverage for src/cvx/quadprog/_sweep.py: 100%

99 statements