Coverage for src/cvx/quadprog/_solve.py: 100%
115 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-10 13:50 +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.
12What is here is the driver and the outer loop: choosing the constraint to enter,
13and deciding when the problem is infeasible. One pass of the inner loop lives in
14:mod:`._steps`, the work before the first iteration in :mod:`._setup`, the
15constraint-shape detection in :mod:`._structure`, and the records they are all
16written in terms of in :mod:`._base`.
18References:
19 [1] D. Goldfarb and A. Idnani (1983). A numerically stable dual method for
20 solving strictly convex quadratic programs. Mathematical Programming,
21 27, 1-33.
22"""
24# G, C, R and J are the names used in Goldfarb & Idnani (1983) and in the
25# reference implementation's public signature `solve_qp(G, a, C, b, meq)`.
26# Lowercasing them would obscure the correspondence to the paper and break
27# drop-in compatibility, so the pep8-naming rules are waived here. TRY003 goes
28# with them: the ValueError messages are reproduced verbatim from the reference
29# so that callers matching on the text keep working.
30#
31# These live in the files rather than in a [lint.per-file-ignores] block because
32# ruff.toml is template-owned -- a local edit to it is reverted by the next
33# `/rhiza:update` sync and flagged as non-template by stage_synced.py.
34# ruff: noqa: N803, N806
36import numpy as np
38from . import _pdas, _threads
39from ._base import _EMPTY, VSMALL, Solution, _WarmEntry
40from ._qr import qr_insert
41from ._setup import _factorize, _validate
42from ._steps import _drop_constraint, _dual_step_limit, _step_choice, _step_directions
43from ._structure import _analyse_constraints, _default_constraints, _slack_evaluator
45__all__ = ["Solution", "solve_qp"]
48# How far above the arithmetic's noise floor a violation must sit before the
49# solver is willing to call the problem infeasible. Used only by
50# _is_spurious_violation, on the one path that would otherwise raise; constraint
51# selection keeps the reference's VSMALL snap untouched.
52#
53# The constant is loose on purpose, and can afford to be. Selection has to
54# separate "rounding" from "violated but tiny", which admits no safe margin --
55# a real violation can be arbitrarily small. This test separates "rounding" from
56# "provably infeasible", and infeasibility is macroscopic: the violation is set
57# by the geometry of the constraints, not by the arithmetic. Any threshold
58# between the two works, so there is nothing here to tune. 32 leaves two decades
59# over the worst residual observed (8 * eps against VSMALL's 6.43 * eps, #36)
60# and stays some thirteen orders below a genuine infeasibility.
61_NOISE_MARGIN = 32.0
64# What is left here is the dual method's add/drop state machine: an outer loop
65# choosing the most violated constraint, an inner loop walking to its boundary
66# while dropping constraints whose multipliers would turn negative. The work
67# inside each pass is delegated -- _step_directions, _dual_step_limit,
68# _step_choice and _drop_constraint -- so this function reads as the algorithm
69# Goldfarb and Idnani specify rather than as the arithmetic implementing it.
70#
71# The inner loop's helpers are called once per *inner* iteration, which is the
72# hot path, so the split was benchmarked rather than assumed. Twelve interleaved
73# A/B rounds on box-constrained problems put the cost at about 1% for n <= 25
74# and nothing measurable from n = 50 up; the paired per-round differences
75# straddle zero (-1.5% to +3.0%), so 1% is the right order but the sign is only
76# just resolvable. Two extra Python calls per iteration are small against the
77# ~18 NumPy dispatches that already dominate at those sizes.
78#
79# That is the price of every block in the package rating B or better. If a
80# future change makes small-n dispatch matter more than it does today, inlining
81# _step_choice back into the loop is the first thing to undo.
82#
83# One thing the split must not do is hoist `C[:, iadd - 1]` out of the `unit`
84# branch to give the helpers a uniform signature: on a box-constrained problem
85# every column is a unit column, so that view would be built once per outer
86# iteration and never read. That regressed n = 10 by ~2% when tried.
87def solve_qp(
88 G: np.ndarray,
89 a: np.ndarray,
90 C: np.ndarray | None = None,
91 b: np.ndarray | None = None,
92 meq: int = 0,
93 factorized: bool = False,
94 check_finite: bool = False,
95 fast: bool = False,
96 blas_threads: int | None = None,
97) -> Solution:
98 r"""Solve a strictly convex quadratic program.
100 Minimises :math:`\tfrac{1}{2} x^T G x - a^T x` subject to
101 :math:`C^T x \ge b`, with the first ``meq`` constraints held as equalities.
103 The example below is chosen so its answer can be written down rather than
104 discovered: with ``G`` the identity the objective separates, and each
105 coordinate reduces to minimising :math:`\tfrac{1}{2} x_i^2 - a_i x_i` over
106 :math:`x_i \ge 0`. That is the unconstrained minimiser with its negative
107 entries clipped to zero -- the projection of ``a`` onto the non-negative
108 orthant.
110 >>> import numpy as np
111 >>> from cvx.quadprog import solve_qp
112 >>> G = np.eye(3)
113 >>> a = np.array([1.0, -2.0, 3.0])
114 >>> C = np.eye(3)
115 >>> b = np.zeros(3)
116 >>> solution = solve_qp(G, a, C, b)
117 >>> bool(np.allclose(solution.x, [1.0, 0.0, 3.0]))
118 True
120 The other fields describe the same solve. ``xu`` is the unconstrained
121 minimiser :math:`G^{-1} a`, which here is ``a`` itself; ``iact`` names the
122 constraints that ended up binding, 1-based -- only the second, since it is
123 the only one ``xu`` violates; and ``f`` is the objective at ``x``.
125 >>> bool(np.allclose(solution.xu, a))
126 True
127 >>> solution.iact.tolist()
128 [2]
129 >>> round(solution.f, 12)
130 -5.0
132 Args:
133 G: See :func:`_solve_with_factors`.
134 a: See :func:`_solve_with_factors`.
135 C: See :func:`_solve_with_factors`.
136 b: See :func:`_solve_with_factors`.
137 meq: See :func:`_solve_with_factors`.
138 factorized: See :func:`_solve_with_factors`.
139 check_finite: See :func:`_solve_with_factors`.
140 fast: Offer the problem to the primal-dual active-set path in
141 :mod:`._pdas` before walking it exactly. That path guesses the whole
142 active set at once and is checked against the KKT conditions, so it
143 returns **the same minimiser or nothing at all** -- when it declines,
144 the exact walk runs and the result is bit-for-bit what it would have
145 been. Measured 1.0x to 5.0x faster, growing with ``n``, because the
146 exact walk's iteration count grows with the active set where this
147 stays at two to four repairs.
149 It is not *uniformly* faster, which is the other reason it is opt-in.
150 Where the exact walk happens to converge in one or two iterations --
151 a box-constrained problem whose unconstrained minimum is nearly
152 feasible, say -- there is nothing to save, and the factorisation and
153 certificate this path pays for anyway make it up to 20% slower. Those
154 are also the cheapest solves there are, so the loss is a handful of
155 microseconds against the hundreds this saves elsewhere.
157 Two reported fields differ when the fast path answers, which is why
158 this is off by default. ``iterations`` counts working-set additions
159 and removals of a *different algorithm*, so it no longer matches the
160 reference implementation's, and ``iact`` is ordered by constraint
161 index rather than by insertion. ``x``, ``f``, ``xu`` and
162 ``lagrangian`` are unaffected. The path is skipped entirely when
163 ``factorized`` is set, since the certificate needs ``G`` itself.
164 blas_threads: Cap the BLAS thread count for the duration of this call, via
165 a scoped `threadpoolctl <https://github.com/joblib/threadpoolctl>`_
166 context that restores the previous limits on exit. Requires
167 ``threadpoolctl``, an optional dependency; a no-op in effect on
168 Accelerate, which exposes no thread knob to set.
170 **Left unset, threading is touched only where it has been measured to
171 be catastrophic**: on Linux, against an OpenBLAS build, with more
172 threads configured than there are physical cores, and only once ``n``
173 is large enough for the collapse to be reachable -- at which point the
174 count is capped to the physical core count. Everywhere else nothing is
175 changed on the caller's behalf, because there is no default worth
176 having: the best count differs by BLAS in opposite directions -- the
177 fast path wants 4 threads on OpenBLAS, where 16 reads 0.05x, and 16 on
178 MKL, where it is still improving -- and by path, since every
179 contributed Windows exact-path sweep is best at 1. See
180 :func:`~cvx.quadprog._threads.auto_cap_threads` for the gate, and #66
181 for the measurements behind it.
183 Set it explicitly to override that, in either direction: an explicit
184 count is used as given and the automatic gate is not consulted. Worth
185 doing on MKL, where more threads than cores is not the trap it is on
186 OpenBLAS, or to pin a solve to 1. Not worth doing around a small solve:
187 ``threadpoolctl`` costs ~100 microseconds against a 0.2 ms solve at
188 ``n = 10``, and for a batch of solves one context around the batch is
189 cheaper than one per call.
191 Returns:
192 The solution.
194 This is a thin wrapper over :func:`_solve_with_factors`, which additionally
195 returns the factorisation it ends on. Nothing about the solve differs; the
196 factors are simply discarded here, because for a single problem they are dead
197 state and ``J`` alone is ``n^2`` doubles -- 15.7 MB at ``n = 1400``, against
198 the 33 KB of the :class:`Solution` itself. :class:`~cvx.quadprog.Sweep` keeps
199 them instead, which is the whole reason the split exists.
201 Setting ``fast`` additionally offers the problem to :mod:`._pdas` first. See
202 the argument's own documentation for what that changes and what it does not.
204 Raises:
205 ValueError: As :func:`_solve_with_factors`, and if ``blas_threads`` is not
206 at least 1.
207 ImportError: If ``blas_threads`` is given and ``threadpoolctl`` is not
208 installed.
209 """
210 if blas_threads is not None:
211 with _threads.limit(blas_threads):
212 return _dispatch(G, a, C, b, meq, factorized, check_finite, fast)
214 n = np.shape(G)[0] if len(np.shape(G)) > 0 else 0
215 auto_threads = _threads.auto_cap_threads(n, fast=fast)
216 if auto_threads is not None:
217 with _threads.limit(auto_threads):
218 return _dispatch(G, a, C, b, meq, factorized, check_finite, fast)
220 return _dispatch(G, a, C, b, meq, factorized, check_finite, fast)
223def _dispatch(
224 G: np.ndarray,
225 a: np.ndarray,
226 C: np.ndarray | None,
227 b: np.ndarray | None,
228 meq: int,
229 factorized: bool,
230 check_finite: bool,
231 fast: bool,
232) -> Solution:
233 """Offer the problem to the fast path if asked, then walk it exactly.
235 Split out of :func:`solve_qp` only so that ``blas_threads`` can wrap both
236 paths in one context manager without the body being written twice. Every
237 argument means what it does there.
239 Args:
240 G: See :func:`_solve_with_factors`.
241 a: See :func:`_solve_with_factors`.
242 C: See :func:`_solve_with_factors`.
243 b: See :func:`_solve_with_factors`.
244 meq: See :func:`_solve_with_factors`.
245 factorized: See :func:`_solve_with_factors`.
246 check_finite: See :func:`_solve_with_factors`.
247 fast: See :func:`solve_qp`.
249 Returns:
250 The solution.
251 """
252 if fast and not factorized and C is not None and b is not None:
253 solution = _pdas._fast_solution(G, a, C, b, meq, check_finite)
254 if solution is not None:
255 return solution
257 solution, _J, _R = _solve_with_factors(G, a, C, b, meq, factorized, check_finite)
258 return solution
261# The outer loop, with the inner loop in `_walk` and the cold start in `_start`. As one
262# function it rated B (10) under `radon cc`, the worst block in the package; split, the
263# loop is A (3) and `_walk` A (5) (#139). The working state stays in the arrays the
264# two share and mutates in place, so what crosses the call is the scalars -- `nact`,
265# `obj`, `nign` -- and the entering constraint's description.
266#
267# The split costs one call per *outer* iteration, never per inner one, and no record:
268# `_entering` already returned a tuple, which `_walk` takes as it is. #98 tried the same
269# split around an allocated record and measured 3-4% for n <= 100. This version was
270# measured against the single function in 15 interleaved A/B rounds per size on
271# box-constrained problems (arm64, Accelerate): paired medians from 0.997 to 1.013 at
272# n = 10 to 200, both run orders, with every per-round range straddling 1. That is
273# nothing resolvable, and `benchmarks/ref_probe.py` reproduces the README's exact-path
274# rows and the n ~ 135 crossover unchanged. Results, iteration counts and infeasibility
275# verdicts are bit-for-bit those of the single function, on 3000 random problems and
276# two `Sweep` families.
277def _solve_with_factors(
278 G: np.ndarray,
279 a: np.ndarray,
280 C: np.ndarray | None = None,
281 b: np.ndarray | None = None,
282 meq: int = 0,
283 factorized: bool = False,
284 check_finite: bool = False,
285 warm: _WarmEntry | None = None,
286) -> tuple[Solution, np.ndarray, np.ndarray]:
287 r"""Solve a strictly convex quadratic program, returning the factorisation too.
289 .. math::
290 \min_x \tfrac{1}{2} x^T G x - a^T x \quad\text{subject to}\quad C^T x \ge b
292 The first ``meq`` constraints are treated as equalities.
294 Args:
295 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
296 If ``factorized`` is True, pass :math:`R^{-1}` instead, where
297 :math:`G = R^T R` with ``R`` upper triangular.
298 a: ``(n,)`` vector of the linear term.
299 C: ``(n, m)`` constraint matrix, one column per constraint. Defaults to
300 a single inactive constraint, giving the unconstrained problem.
301 b: ``(m,)`` right-hand side of the constraints.
302 meq: Number of leading constraints to treat as equalities.
303 factorized: Whether ``G`` holds :math:`R^{-1}` rather than :math:`G`.
304 check_finite: Whether to reject NaN and infinity in the inputs. **Off by
305 default**, matching the reference, which does not scan either: the
306 check is :math:`O(n^2)` on ``G`` and callers that already validate
307 their data should not pay for it. Left off, a non-finite ``G`` is not
308 diagnosed and what happens next belongs to the LAPACK build --
309 Accelerate reports a failed factorisation, OpenBLAS propagates NaNs
310 into the result. Neither returns a finite wrong answer, but only one
311 of them raises, so a program that must behave identically on every
312 platform should pass True.
313 warm: A dual-feasible state to resume from, in place of the cold start at
314 the unconstrained minimum. See :class:`_WarmEntry` for the invariant
315 it must satisfy, which is the caller's to establish; from there the
316 iteration cannot tell a resumed state from a cold one.
318 Returns:
319 A :class:`Solution` with the minimiser, the objective value, the
320 unconstrained minimiser, the iteration counts, the Lagrange multipliers
321 and the active set; together with ``J``, the inverse Cholesky factor as
322 the iteration left it, and ``R``, the packed triangular factor of the
323 active constraint normals. Both are live internal buffers, freshly
324 allocated by this call and not aliased to anything the caller passed in.
326 Raises:
327 ValueError: If the shapes are inconsistent, if ``meq`` is out of range,
328 if ``G`` is not positive definite, if the constraints admit no
329 solution, or if ``check_finite`` is set and any input holds a
330 non-finite value.
331 """
332 G = np.asarray(G, dtype=np.float64)
333 a = np.asarray(a, dtype=np.float64)
334 C, b, meq = _default_constraints(G, C, b, meq)
336 n, q = _validate(G, a, C, b, meq, check_finite)
337 r = min(n, q)
339 # The norm of each column of C, used to scale the pivoting rule so that the
340 # choice of constraint is invariant to how each one happens to be scaled.
341 # A zero-norm column reads 0 >= b, which no x can influence; scoring it as
342 # infinitely violated sends the solver to the infeasibility verdict instead
343 # of dividing by zero below.
344 # einsum rather than `np.sum(C * C, axis=0)`, which materialises a whole
345 # n x m temporary to reduce it away again -- 0.32 ms against 0.19 ms at
346 # n = 800, and 3.48 ms against 0.69 ms at n = 1400, where the temporary is
347 # 31 MB and stops fitting anywhere useful (#108).
348 #
349 # The two are not bit-identical on a dense C: the reduction order differs, so
350 # a column norm can land one ulp apart. That cannot change what the solver
351 # computes -- nbv only scales the selection score below, and the arithmetic of
352 # the solve itself never sees it -- and it cannot change which constraint is
353 # chosen except between two whose scores already agree to within an ulp,
354 # where either is a legitimate choice. The two shapes where selection order is
355 # load-bearing are both exactly equal under either form: a single-nonzero
356 # column has one term and no order to differ over, and duplicated or dependent
357 # columns produce identical values to each other under both, so argmax still
358 # resolves the tie towards the lowest index.
359 nbv = np.sqrt(np.einsum("ij,ij->j", C, C))
360 degenerate = nbv == 0.0
361 nbv_safe = np.where(degenerate, 1.0, nbv)
363 # Sparsity of C, detected once. Bound constraints make most columns a single
364 # scaled unit vector, which turns three of the per-iteration operations from
365 # O(n) or O(n*q) work into scalar indexing.
366 single, srow, sval = _analyse_constraints(C)
367 slack_of = _slack_evaluator(C, single, srow, sval)
369 J, R, iact, nact, xv, uv, obj, xu = _start(G, a, factorized, warm, r, q)
371 lagr = np.zeros(q)
372 iter_full, iter_partial = 0, 0
374 # Constraints found violated only by rounding at the current xv, which the
375 # iteration can neither enforce nor draw a conclusion from -- see
376 # _is_spurious_violation. Held 0-based, first nign entries valid, and reset
377 # whenever xv moves, so nothing is masked on the strength of a stale iterate.
378 ignored = np.zeros(q, dtype=np.int64)
379 nign = 0
381 while True:
382 iter_full += 1
384 # The slack of every constraint. Slacks of active constraints are forced
385 # to exactly zero as a safeguard against rounding error.
386 sv = slack_of(xv) - b
387 sv[np.abs(sv) < VSMALL] = 0.0
388 sv[iact[:nact] - 1] = 0.0
389 sv[ignored[:nign]] = 0.0
391 iadd = _choose_constraint(sv, nbv_safe, degenerate, meq)
393 if iadd == 0:
394 # Every constraint is satisfied, so we are at the optimum.
395 lagr[iact[:nact] - 1] = uv[:nact]
396 iterations = np.array([iter_full, iter_partial], dtype=np.int64)
397 return Solution(xv, obj, xu, iterations, lagr, iact[:nact]), J, R
399 entering = _entering(C, single, srow, sval, iadd)
400 slack, normal_norm, rhs = float(sv[iadd - 1]), float(nbv_safe[iadd - 1]), float(b[iadd - 1])
401 nact, obj, nign, partial = _walk(
402 J, R, uv, iact, nact, xv, obj, ignored, nign, meq, iadd, slack, entering, normal_norm, rhs
403 )
404 iter_partial += partial
407def _start(G: np.ndarray, a: np.ndarray, factorized: bool, warm: _WarmEntry | None, r: int, q: int) -> _WarmEntry:
408 """Return the dual-feasible state the outer loop begins from.
410 Args:
411 G: See :func:`_solve_with_factors`.
412 a: See :func:`_solve_with_factors`.
413 factorized: See :func:`_solve_with_factors`.
414 warm: See :func:`_solve_with_factors`. Returned as it is when given.
415 r: Largest possible active set, ``min(n, q)``.
416 q: Number of constraints.
418 Returns:
419 The state to iterate from, in the fields of :class:`_WarmEntry`.
420 """
421 if warm is not None:
422 # Resuming from a state a caller already holds. It must satisfy the same
423 # invariant the cold start gets for free -- see _WarmEntry -- and from
424 # here the loop cannot tell the two apart.
425 return warm
427 # Cold start. xv holds G^-1 a, the unconstrained minimum, and J holds R^-1 so
428 # that J J^T = G^-1; the active set is empty, which is trivially dual feasible
429 # and is the whole reason this method needs no phase 1. The objective is kept
430 # as a running total, each step updating it in closed form rather than
431 # re-evaluating the quadratic. R is upper triangular stored as packed columns
432 # -- see the note in _qr.
433 J, xv = _factorize(G, a, factorized)
434 obj = -float(a @ xv) / 2.0
435 R = np.zeros(r * (r + 1) // 2)
436 uv = np.zeros(r) # dual variables of the active constraints
437 iact = np.zeros(q, dtype=np.int64) # 1-based, first nact entries valid
438 return _WarmEntry(J, R, iact, 0, xv, uv, obj, xv.copy())
441def _walk(
442 J: np.ndarray,
443 R: np.ndarray,
444 uv: np.ndarray,
445 iact: np.ndarray,
446 nact: int,
447 xv: np.ndarray,
448 obj: float,
449 ignored: np.ndarray,
450 nign: int,
451 meq: int,
452 iadd: int,
453 slack: float,
454 entering: tuple[bool, int, float, np.ndarray],
455 normal_norm: float,
456 rhs: float,
457) -> tuple[int, float, int, int]:
458 """Walk towards the entering constraint's boundary: the inner loop.
460 Each pass takes the largest step that keeps every multiplier non-negative.
461 If that reaches the boundary the constraint is added and the walk ends; if a
462 multiplier hits zero first, its constraint is dropped and the walk goes on.
464 The arrays -- ``J``, ``R``, ``uv``, ``iact``, ``xv`` and ``ignored`` -- are
465 the outer loop's working state and are updated in place. The scalars cannot
466 be, so the ones the walk changes come back as the return value.
468 Args:
469 J: Inverse Cholesky factor, updated for the active set.
470 R: Packed triangular factor of the active constraint normals.
471 uv: Multipliers of the active constraints.
472 iact: 1-based active set, first ``nact`` entries valid.
473 nact: Size of the active set.
474 xv: The iterate.
475 obj: Objective value at ``xv``.
476 ignored: 0-based constraints set aside as violated only by rounding,
477 first ``nign`` entries valid -- see :func:`_is_spurious_violation`.
478 nign: Number of valid entries in ``ignored``.
479 meq: Number of leading constraints treated as equalities.
480 iadd: 1-based index of the entering constraint.
481 slack: Its slack at ``xv``.
482 entering: How to read its normal, as returned by :func:`_entering`.
483 normal_norm: Norm of its normal, zeros replaced by one.
484 rhs: Its right-hand side.
486 Returns:
487 ``nact``, ``obj`` and ``nign`` as the walk leaves them, and the number of
488 partial steps it took.
489 """
490 unit, row, val, normal = entering
491 # An equality constraint may be violated from either side. When its slack is
492 # positive we have to step in the opposite direction.
493 reverse_step = slack > 0.0
494 u = 0.0
495 partial = 0
497 while True:
498 dv, zv, rv, ztn = _step_directions(J, R, nact, unit, val, row, normal)
500 # The largest step t1 that keeps the dual variables non-negative, and the
501 # constraint idel that would be the first to bind at zero.
502 t1, idel = _dual_step_limit(uv, rv, iact, nact, meq, reverse_step)
504 if _is_spurious_violation(ztn, idel, slack, normal_norm, xv, rhs):
505 # Satisfied to within the accuracy of xv, but the primal cannot move
506 # and no multiplier can be reduced. Enforcing it would be a no-op and
507 # concluding infeasibility from it would be wrong, so set it aside and
508 # let the outer loop take the next candidate.
509 ignored[nign] = iadd - 1
510 return nact, obj, nign + 1, partial
512 step, full_step = _step_choice(ztn, slack, t1, idel == 0, reverse_step)
514 if ztn is not None:
515 xv += step * zv
516 obj += step * ztn * (step / 2.0 + u)
517 # xv moved, so every slack set aside against the old one is stale and
518 # must be measured again.
519 nign = 0
521 uv[:nact] -= step * rv
522 u += step
524 if full_step:
525 # The entering constraint now holds with equality: add it.
526 nact += 1
527 uv[nact - 1], iact[nact - 1] = u, iadd
528 qr_insert(nact, dv, J, R)
529 return nact, obj, nign, partial
531 # Only a partial step: drop constraint idel from the active set.
532 nact = _drop_constraint(idel, nact, uv, iact, J, R)
533 partial += 1
534 slack = _slack_after_drop(slack, ztn, entering, xv, rhs)
537def _slack_after_drop(
538 slack: float, ztn: float | None, entering: tuple[bool, int, float, np.ndarray], xv: np.ndarray, rhs: float
539) -> float:
540 """Return the entering constraint's slack after a partial step.
542 Args:
543 slack: Its slack before the step.
544 ztn: Rate at which the step closed it, or None when the primal did not move.
545 entering: How to read its normal, as returned by :func:`_entering`.
546 xv: The iterate after the step.
547 rhs: Its right-hand side.
549 Returns:
550 ``slack`` itself when the primal did not move, and otherwise the slack
551 measured again at the new ``xv``.
552 """
553 if ztn is None:
554 return slack
555 unit, row, val, normal = entering
556 reached = val * float(xv[row]) if unit else float(xv @ normal)
557 return reached - rhs
560def _is_spurious_violation(
561 ztn: float | None, idel: int, slack: float, normal_norm: float, xv: np.ndarray, rhs: float
562) -> bool:
563 """Return whether a stuck iteration reflects rounding rather than infeasibility.
565 The iteration is *stuck* when the primal cannot move (``ztn`` is None, so
566 the entering normal already lies in the span of the active set) and no
567 active multiplier can be reduced (``idel`` is 0). Goldfarb and Idnani's
568 conclusion from that pair is that the dual is unbounded and the primal
569 therefore infeasible -- but the argument assumes the entering constraint is
570 genuinely violated. When its violation is the size of the rounding in
571 ``xv``, it is not, and the conclusion does not follow.
573 That is reachable here rather than being theoretical. ``qr_insert`` reduces
574 with a Householder reflection where the reference chases Givens rotations;
575 the two agree in exact arithmetic (see the README) but round differently, so
576 an iterate that the reference leaves 4.68 * eps inside a constraint can land
577 8 * eps outside it -- either side of the fixed ``VSMALL`` snap applied to the
578 slacks in ``solve_qp``.
580 Args:
581 ztn: Rate at which the entering constraint's slack closes, or None when
582 the primal cannot move.
583 idel: 1-based position of the constraint limiting the dual step, or 0
584 when nothing limits it.
585 slack: Current slack of the entering constraint.
586 normal_norm: Norm of the entering constraint's normal, zeros replaced
587 by one.
588 xv: Current primal iterate, whose magnitude sets the scale of the
589 rounding the slack inherits.
590 rhs: The entering constraint's right-hand side.
592 Returns:
593 True when the violation is indistinguishable from rounding, so the
594 constraint should be set aside rather than treated as proof of
595 infeasibility.
596 """
597 if ztn is not None or idel != 0:
598 return False
600 # Reached only on the stuck path, so the O(n) norm is off the hot loop. The
601 # slack inherits the error in xv rather than merely the error of its own dot
602 # product, so the scale that matters is ||c|| ||x||, not the size of the
603 # terms that formed it -- for a constraint that x sits on, those are already
604 # at the noise floor and say nothing.
605 scale = normal_norm * float(np.max(np.abs(xv))) + abs(rhs)
606 return abs(slack) <= _NOISE_MARGIN * VSMALL * max(scale, 1.0)
609def _entering(
610 C: np.ndarray, single: np.ndarray, srow: np.ndarray, sval: np.ndarray, iadd: int
611) -> tuple[bool, int, float, np.ndarray]:
612 """Return how to read the entering constraint's normal.
614 A column holding a single scaled unit vector ``e_row`` lets the three products
615 against it in :func:`_step_directions` be read off by index instead of
616 computed.
618 Both branches bind all four values, so ``_step_directions`` can take a fixed
619 signature -- but only one branch builds the dense view, and that asymmetry is
620 the point. On a box-constrained problem every column is a unit column, so
621 hoisting ``C[:, iadd - 1]`` out of the branch would construct a strided view
622 per outer iteration that nothing ever reads; measured at ~2% on ``n = 10``.
623 Lifting the branch into this function keeps that property, since the view is
624 still built only where it is used.
626 Args:
627 C: ``(n, m)`` constraint matrix.
628 single: Mask of the columns holding exactly one nonzero.
629 srow: Row index of that nonzero per column.
630 sval: Value of that nonzero per column.
631 iadd: 1-based index of the entering constraint.
633 Returns:
634 Whether the column is a scaled unit vector, the row its nonzero occupies,
635 that nonzero's value, and the dense normal -- the last three meaningful
636 only in the branch that binds them.
637 """
638 if single[iadd - 1]:
639 return True, int(srow[iadd - 1]), float(sval[iadd - 1]), _EMPTY
640 return False, 0, 0.0, C[:, iadd - 1]
643def _choose_constraint(sv: np.ndarray, nbv_safe: np.ndarray, degenerate: np.ndarray, meq: int) -> int:
644 """Return the 1-based index of the most violated constraint, or 0 if none.
646 Violations are measured relative to the norm of the constraint normal, so
647 the choice does not depend on the scaling of individual constraints. An
648 equality constraint counts as violated in either direction.
650 Scanning for the largest violation is equivalent to taking an ``argmax``,
651 which resolves ties towards the lowest index just as a left-to-right scan
652 with a strict improvement test does.
654 Args:
655 sv: Slack of each constraint.
656 nbv_safe: Norm of each constraint normal, with zeros replaced by one.
657 degenerate: Mask of the constraints whose normal has zero norm.
658 meq: Number of leading constraints treated as equalities.
660 Returns:
661 The 1-based index of the constraint to add, or 0 at the optimum.
662 """
663 # An inequality is violated when its slack is negative; an equality whenever
664 # its slack is nonzero.
665 violation = -sv
666 np.abs(violation[:meq], out=violation[:meq])
668 score = violation / nbv_safe
669 # A zero-norm normal cannot be satisfied by any step, so rank it first.
670 if degenerate.any():
671 score = np.where(degenerate & (violation > 0.0), np.inf, score)
673 iadd = int(np.argmax(score))
674 return iadd + 1 if score[iadd] > 0.0 else 0