Coverage for src/cvx/quadprog/_solve.py: 100%
106 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-15 05:55 +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# Measured, and deliberately left alone. `radon cc` puts this function at B (10) -- the
262# worst block in the package, against an average of A (3.96) over 54 blocks with nothing
263# at C or worse -- and `radon mi` puts this module lowest of the ten, in the mid-20s of
264# band A. Those MI decimals are quoted as a band rather than a figure on purpose: MI folds
265# in Halstead volume, so a comment like this one lowers the very number it reports. The
266# CC is comment-invariant; the MI is not.
267#
268# It stays one function because the thing it transcribes is one thing: the outer loop of
269# Goldfarb & Idnani (1983), whose steps share the working state (`J`, packed `R`, `nact`,
270# the active index vector) and read in the paper's order. Splitting it would move that
271# state into arguments threaded through helpers that are only ever called once, in
272# sequence, from here -- trading a legible transcription for a less legible one and making
273# the correspondence to the paper harder to check, which is the property this file is
274# organised around. The per-iteration work already lives in `_steps.py` and `_qr.py`; what
275# is left is the loop itself.
276#
277# So the B (10) is recorded rather than removed. The signal to revisit is a rank change --
278# this block reaching C, or the module leaving MI band A (below 20) -- not the decimals.
279def _solve_with_factors(
280 G: np.ndarray,
281 a: np.ndarray,
282 C: np.ndarray | None = None,
283 b: np.ndarray | None = None,
284 meq: int = 0,
285 factorized: bool = False,
286 check_finite: bool = False,
287 warm: _WarmEntry | None = None,
288) -> tuple[Solution, np.ndarray, np.ndarray]:
289 r"""Solve a strictly convex quadratic program, returning the factorisation too.
291 .. math::
292 \min_x \tfrac{1}{2} x^T G x - a^T x \quad\text{subject to}\quad C^T x \ge b
294 The first ``meq`` constraints are treated as equalities.
296 Args:
297 G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
298 If ``factorized`` is True, pass :math:`R^{-1}` instead, where
299 :math:`G = R^T R` with ``R`` upper triangular.
300 a: ``(n,)`` vector of the linear term.
301 C: ``(n, m)`` constraint matrix, one column per constraint. Defaults to
302 a single inactive constraint, giving the unconstrained problem.
303 b: ``(m,)`` right-hand side of the constraints.
304 meq: Number of leading constraints to treat as equalities.
305 factorized: Whether ``G`` holds :math:`R^{-1}` rather than :math:`G`.
306 check_finite: Whether to reject NaN and infinity in the inputs. **Off by
307 default**, matching the reference, which does not scan either: the
308 check is :math:`O(n^2)` on ``G`` and callers that already validate
309 their data should not pay for it. Left off, a non-finite ``G`` is not
310 diagnosed and what happens next belongs to the LAPACK build --
311 Accelerate reports a failed factorisation, OpenBLAS propagates NaNs
312 into the result. Neither returns a finite wrong answer, but only one
313 of them raises, so a program that must behave identically on every
314 platform should pass True.
315 warm: A dual-feasible state to resume from, in place of the cold start at
316 the unconstrained minimum. See :class:`_WarmEntry` for the invariant
317 it must satisfy, which is the caller's to establish; from there the
318 iteration cannot tell a resumed state from a cold one.
320 Returns:
321 A :class:`Solution` with the minimiser, the objective value, the
322 unconstrained minimiser, the iteration counts, the Lagrange multipliers
323 and the active set; together with ``J``, the inverse Cholesky factor as
324 the iteration left it, and ``R``, the packed triangular factor of the
325 active constraint normals. Both are live internal buffers, freshly
326 allocated by this call and not aliased to anything the caller passed in.
328 Raises:
329 ValueError: If the shapes are inconsistent, if ``meq`` is out of range,
330 if ``G`` is not positive definite, if the constraints admit no
331 solution, or if ``check_finite`` is set and any input holds a
332 non-finite value.
333 """
334 G = np.asarray(G, dtype=np.float64)
335 a = np.asarray(a, dtype=np.float64)
336 C, b, meq = _default_constraints(G, C, b, meq)
338 n, q = _validate(G, a, C, b, meq, check_finite)
339 r = min(n, q)
341 # The norm of each column of C, used to scale the pivoting rule so that the
342 # choice of constraint is invariant to how each one happens to be scaled.
343 # A zero-norm column reads 0 >= b, which no x can influence; scoring it as
344 # infinitely violated sends the solver to the infeasibility verdict instead
345 # of dividing by zero below.
346 # einsum rather than `np.sum(C * C, axis=0)`, which materialises a whole
347 # n x m temporary to reduce it away again -- 0.32 ms against 0.19 ms at
348 # n = 800, and 3.48 ms against 0.69 ms at n = 1400, where the temporary is
349 # 31 MB and stops fitting anywhere useful (#108).
350 #
351 # The two are not bit-identical on a dense C: the reduction order differs, so
352 # a column norm can land one ulp apart. That cannot change what the solver
353 # computes -- nbv only scales the selection score below, and the arithmetic of
354 # the solve itself never sees it -- and it cannot change which constraint is
355 # chosen except between two whose scores already agree to within an ulp,
356 # where either is a legitimate choice. The two shapes where selection order is
357 # load-bearing are both exactly equal under either form: a single-nonzero
358 # column has one term and no order to differ over, and duplicated or dependent
359 # columns produce identical values to each other under both, so argmax still
360 # resolves the tie towards the lowest index.
361 nbv = np.sqrt(np.einsum("ij,ij->j", C, C))
362 degenerate = nbv == 0.0
363 nbv_safe = np.where(degenerate, 1.0, nbv)
365 # Sparsity of C, detected once. Bound constraints make most columns a single
366 # scaled unit vector, which turns three of the per-iteration operations from
367 # O(n) or O(n*q) work into scalar indexing.
368 single, srow, sval = _analyse_constraints(C)
369 slack_of = _slack_evaluator(C, single, srow, sval)
371 if warm is None:
372 # Cold start. xv holds G^-1 a, the unconstrained minimum, and J holds
373 # R^-1 so that J J^T = G^-1; the active set is empty, which is trivially
374 # dual feasible and is the whole reason this method needs no phase 1. The
375 # objective is kept as a running total, each step updating it in closed
376 # form rather than re-evaluating the quadratic. R is upper triangular
377 # stored as packed columns -- see the note in _qr.
378 J, xv = _factorize(G, a, factorized)
379 obj = -float(a @ xv) / 2.0
380 xu = xv.copy()
381 R = np.zeros(r * (r + 1) // 2)
382 uv = np.zeros(r) # dual variables of the active constraints
383 iact = np.zeros(q, dtype=np.int64) # 1-based, first nact entries valid
384 nact = 0
385 else:
386 # Resuming from a state a caller already holds. It must satisfy the same
387 # invariant the cold start gets for free -- see _WarmEntry -- and from
388 # here the loop cannot tell the two apart.
389 J, R, iact, nact, xv, uv, obj, xu = warm
391 lagr = np.zeros(q)
392 iter_full, iter_partial = 0, 0
394 # Constraints found violated only by rounding at the current xv, which the
395 # iteration can neither enforce nor draw a conclusion from -- see
396 # _is_spurious_violation. Held 0-based, first nign entries valid, and reset
397 # whenever xv moves, so nothing is masked on the strength of a stale iterate.
398 ignored = np.zeros(q, dtype=np.int64)
399 nign = 0
401 while True:
402 iter_full += 1
404 # The slack of every constraint. Slacks of active constraints are forced
405 # to exactly zero as a safeguard against rounding error.
406 sv = slack_of(xv) - b
407 sv[np.abs(sv) < VSMALL] = 0.0
408 sv[iact[:nact] - 1] = 0.0
409 sv[ignored[:nign]] = 0.0
411 iadd = _choose_constraint(sv, nbv_safe, degenerate, meq)
413 if iadd == 0:
414 # Every constraint is satisfied, so we are at the optimum.
415 lagr[iact[:nact] - 1] = uv[:nact]
416 iterations = np.array([iter_full, iter_partial], dtype=np.int64)
417 return Solution(xv, obj, xu, iterations, lagr, iact[:nact]), J, R
419 # An equality constraint may be violated from either side. When its
420 # slack is positive we have to step in the opposite direction.
421 slack = float(sv[iadd - 1])
422 reverse_step = slack > 0.0
423 u = 0.0
425 unit, row, val, normal = _entering(C, single, srow, sval, iadd)
427 # Inner loop: walk towards the constraint boundary, dropping active
428 # constraints whose multipliers would otherwise turn negative.
429 while True:
430 dv, zv, rv, ztn = _step_directions(J, R, nact, unit, val, row, normal)
432 # The largest step t1 that keeps the dual variables non-negative,
433 # and the constraint idel that would be the first to bind at zero.
434 t1, idel = _dual_step_limit(uv, rv, iact, nact, meq, reverse_step)
436 if _is_spurious_violation(ztn, idel, slack, nbv_safe[iadd - 1], xv, b[iadd - 1]):
437 # Satisfied to within the accuracy of xv, but the primal cannot
438 # move and no multiplier can be reduced. Enforcing it would be a
439 # no-op and concluding infeasibility from it would be wrong, so
440 # set it aside and let the outer loop take the next candidate.
441 ignored[nign] = iadd - 1
442 nign += 1
443 break
445 step, full_step = _step_choice(ztn, slack, t1, idel == 0, reverse_step)
447 if ztn is not None:
448 xv += step * zv
449 obj += step * ztn * (step / 2.0 + u)
450 # xv moved, so every slack set aside against the old one is
451 # stale and must be measured again.
452 nign = 0
454 uv[:nact] -= step * rv
455 u += step
457 if full_step:
458 # The entering constraint now holds with equality: add it.
459 nact += 1
460 uv[nact - 1], iact[nact - 1] = u, iadd
461 qr_insert(nact, dv, J, R)
462 break
464 # Only a partial step: drop constraint idel from the active set.
465 nact = _drop_constraint(idel, nact, uv, iact, J, R)
466 iter_partial += 1
468 if ztn is not None:
469 # We moved in primal space, so the slack we are closing has
470 # changed and must be recomputed.
471 reached = val * float(xv[row]) if unit else float(xv @ normal)
472 slack = reached - float(b[iadd - 1])
475def _is_spurious_violation(
476 ztn: float | None, idel: int, slack: float, normal_norm: float, xv: np.ndarray, rhs: float
477) -> bool:
478 """Return whether a stuck iteration reflects rounding rather than infeasibility.
480 The iteration is *stuck* when the primal cannot move (``ztn`` is None, so
481 the entering normal already lies in the span of the active set) and no
482 active multiplier can be reduced (``idel`` is 0). Goldfarb and Idnani's
483 conclusion from that pair is that the dual is unbounded and the primal
484 therefore infeasible -- but the argument assumes the entering constraint is
485 genuinely violated. When its violation is the size of the rounding in
486 ``xv``, it is not, and the conclusion does not follow.
488 That is reachable here rather than being theoretical. ``qr_insert`` reduces
489 with a Householder reflection where the reference chases Givens rotations;
490 the two agree in exact arithmetic (see the README) but round differently, so
491 an iterate that the reference leaves 4.68 * eps inside a constraint can land
492 8 * eps outside it -- either side of the fixed ``VSMALL`` snap applied to the
493 slacks in ``solve_qp``.
495 Args:
496 ztn: Rate at which the entering constraint's slack closes, or None when
497 the primal cannot move.
498 idel: 1-based position of the constraint limiting the dual step, or 0
499 when nothing limits it.
500 slack: Current slack of the entering constraint.
501 normal_norm: Norm of the entering constraint's normal, zeros replaced
502 by one.
503 xv: Current primal iterate, whose magnitude sets the scale of the
504 rounding the slack inherits.
505 rhs: The entering constraint's right-hand side.
507 Returns:
508 True when the violation is indistinguishable from rounding, so the
509 constraint should be set aside rather than treated as proof of
510 infeasibility.
511 """
512 if ztn is not None or idel != 0:
513 return False
515 # Reached only on the stuck path, so the O(n) norm is off the hot loop. The
516 # slack inherits the error in xv rather than merely the error of its own dot
517 # product, so the scale that matters is ||c|| ||x||, not the size of the
518 # terms that formed it -- for a constraint that x sits on, those are already
519 # at the noise floor and say nothing.
520 scale = normal_norm * float(np.max(np.abs(xv))) + abs(rhs)
521 return abs(slack) <= _NOISE_MARGIN * VSMALL * max(scale, 1.0)
524def _entering(
525 C: np.ndarray, single: np.ndarray, srow: np.ndarray, sval: np.ndarray, iadd: int
526) -> tuple[bool, int, float, np.ndarray]:
527 """Return how to read the entering constraint's normal.
529 A column holding a single scaled unit vector ``e_row`` lets the three products
530 against it in :func:`_step_directions` be read off by index instead of
531 computed.
533 Both branches bind all four values, so ``_step_directions`` can take a fixed
534 signature -- but only one branch builds the dense view, and that asymmetry is
535 the point. On a box-constrained problem every column is a unit column, so
536 hoisting ``C[:, iadd - 1]`` out of the branch would construct a strided view
537 per outer iteration that nothing ever reads; measured at ~2% on ``n = 10``.
538 Lifting the branch into this function keeps that property, since the view is
539 still built only where it is used.
541 Args:
542 C: ``(n, m)`` constraint matrix.
543 single: Mask of the columns holding exactly one nonzero.
544 srow: Row index of that nonzero per column.
545 sval: Value of that nonzero per column.
546 iadd: 1-based index of the entering constraint.
548 Returns:
549 Whether the column is a scaled unit vector, the row its nonzero occupies,
550 that nonzero's value, and the dense normal -- the last three meaningful
551 only in the branch that binds them.
552 """
553 if single[iadd - 1]:
554 return True, int(srow[iadd - 1]), float(sval[iadd - 1]), _EMPTY
555 return False, 0, 0.0, C[:, iadd - 1]
558def _choose_constraint(sv: np.ndarray, nbv_safe: np.ndarray, degenerate: np.ndarray, meq: int) -> int:
559 """Return the 1-based index of the most violated constraint, or 0 if none.
561 Violations are measured relative to the norm of the constraint normal, so
562 the choice does not depend on the scaling of individual constraints. An
563 equality constraint counts as violated in either direction.
565 Scanning for the largest violation is equivalent to taking an ``argmax``,
566 which resolves ties towards the lowest index just as a left-to-right scan
567 with a strict improvement test does.
569 Args:
570 sv: Slack of each constraint.
571 nbv_safe: Norm of each constraint normal, with zeros replaced by one.
572 degenerate: Mask of the constraints whose normal has zero norm.
573 meq: Number of leading constraints treated as equalities.
575 Returns:
576 The 1-based index of the constraint to add, or 0 at the optimum.
577 """
578 # An inequality is violated when its slack is negative; an equality whenever
579 # its slack is nonzero.
580 violation = -sv
581 np.abs(violation[:meq], out=violation[:meq])
583 score = violation / nbv_safe
584 # A zero-norm normal cannot be satisfied by any step, so rank it first.
585 if degenerate.any():
586 score = np.where(degenerate & (violation > 0.0), np.inf, score)
588 iadd = int(np.argmax(score))
589 return iadd + 1 if score[iadd] > 0.0 else 0