Coverage for src/cvx/quadprog/_steps.py: 100%
52 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-09 18:50 +0000
1"""One pass of the dual method's inner loop.
3Given a direction to move the iterate in, these decide how far it may go
4before a multiplier would turn negative, whether the step is full or partial,
5and which constraint leaves the active set when it is partial.
6"""
8# G, C, R and J are the names used in Goldfarb & Idnani (1983) and in the
9# reference implementation's public signature `solve_qp(G, a, C, b, meq)`.
10# Lowercasing them would obscure the correspondence to the paper, so the
11# pep8-naming rules are waived here, as they are in _solve.py.
12# ruff: noqa: N803, TRY003
14import numpy as np
15from scipy.linalg.blas import dtpsv
17from ._base import _EMPTY, VSMALL
18from ._qr import qr_delete
20# `dtpsv` and `dtrtri` are imported by name rather than resolved through
21# get_blas_funcs/get_lapack_funcs. Those helpers pick a precision from prototype
22# arrays, which is what you want when the caller's dtype is open; here every
23# array is float64 by the time it reaches them -- solve_qp coerces on the way in
24# -- so the choice is already made and resolving it costs an indirection, a
25# module-level prototype array and a cast. Naming the double-precision wrappers
26# yields the identical objects (`get_blas_funcs("tpsv", (f64,)) is dtpsv`) and a
27# sharper static type: mypy reads dtpsv as returning ndarray[float64] where the
28# cast to Callable[..., Any] erased it. Supporting float32 would mean going back.
29#
30# The packed triangular solve runs once per iteration and is the reason R is
31# stored packed: `ap` is an unshaped rank-1 argument, so passing the whole array
32# with n=nact reads the leading triangle in place. The dense equivalent, trtrs on
33# R[:nact, :nact], is handed a strided view and copies it every call -- 77 us
34# against 7.5 us at n = 700. Calling BLAS directly also skips
35# scipy.linalg.solve_triangular's per-call validation, whose check_finite scans
36# the whole array.
38# Active-set size below which the ratio test in _dual_step_limit runs on Python
39# lists rather than arrays. Its work is one division and one argmin over `nact`
40# entries, which for a small active set costs far less than asking NumPy to do
41# it: measured per call at the solver's own call site, the array form takes
42# 2.63 us against 0.75 us for the loop at nact = 4, and 3.29 us against 0.90 us
43# at nact = 9. The two cross where the per-element interpreter cost overtakes
44# NumPy's per-call overhead, a little below sixty.
45_SCALAR_CUTOFF = 50
48def _step_directions(
49 J: np.ndarray, R: np.ndarray, nact: int, unit: bool, val: float, row: int, normal: np.ndarray
50) -> tuple[np.ndarray, np.ndarray, np.ndarray, float | None]:
51 """Return the primal and dual step directions for the entering constraint.
53 Recomputed on every pass of the inner loop because dropping a constraint
54 changes ``J`` and ``R``.
56 Args:
57 J: ``(n, n)`` inverse Cholesky factor, with ``J J^T = G^-1``.
58 R: Packed upper triangular factor of the active constraint normals.
59 nact: Size of the active set.
60 unit: Whether the constraint normal is a single scaled unit vector.
61 val: Its nonzero value, meaningful only when ``unit``.
62 row: The row that nonzero occupies, meaningful only when ``unit``.
63 normal: The constraint normal in dense form.
65 Returns:
66 ``dv``, ``J^T n`` split as ``(d_1, d_2)`` at the size of the active
67 set; ``zv``, the primal step direction; ``rv``, the negated dual step
68 direction; and ``ztn``, the rate at which the entering constraint's
69 slack closes -- ``None`` when the primal cannot move, which is the
70 caller's signal that the step is limited by the dual alone.
71 """
72 # For a unit column this is one scaled row of J, O(n) rather than O(n^2).
73 dv = val * J[row, :] if unit else J.T @ normal
75 # zv = J_2 d_2, the component of the constraint normal orthogonal to the
76 # active set.
77 zv = J[:, nact:] @ dv[nact:]
79 # rv = R^-1 d_1. Solved on a copy: dv is still needed intact for qr_insert.
80 rv = dtpsv(nact, R, dv[:nact].copy(), overwrite_x=True) if nact else _EMPTY
82 if abs(float(zv @ zv)) <= VSMALL:
83 # The primal cannot move, so the entering constraint's slack does not
84 # close at any rate and t2 is infinite.
85 return dv, zv, rv, None
87 ztn = val * float(zv[row]) if unit else float(zv @ normal)
88 return dv, zv, rv, ztn
91def _step_choice(ztn: float | None, slack: float, t1: float, t1inf: bool, reverse_step: bool) -> tuple[float, bool]:
92 """Return the step to take and whether it reaches the entering constraint.
94 Two limits compete: ``t1``, past which an active multiplier would turn
95 negative, and ``t2``, at which the entering constraint's slack closes. The
96 smaller one wins. Reaching ``t2`` ends the inner loop; stopping at ``t1``
97 means dropping a constraint and going round again.
99 Args:
100 ztn: Rate at which the entering constraint's slack closes, or None
101 when the primal cannot move and ``t2`` is therefore infinite.
102 slack: Current slack of the entering constraint.
103 t1: Largest dual-feasible step.
104 t1inf: Whether ``t1`` is unbounded, in which case its value is
105 meaningless.
106 reverse_step: Whether to step in the negative direction, which is the
107 case for an equality constraint violated from above.
109 Returns:
110 The signed step, and whether it is a full step to the constraint.
112 Raises:
113 ValueError: If neither limit is finite, which means the dual is
114 unbounded and so the primal is infeasible.
115 """
116 # Spelled as three outcomes rather than as a boolean built from both
117 # limits, so that the branch establishing t2 is finite is also the branch
118 # that steps to it -- otherwise nothing in the types rules out stepping to
119 # an infinite t2, and a reader has to reconstruct the argument.
120 if ztn is None:
121 if t1inf:
122 # Neither limit is finite: we can step infinitely far, so the dual
123 # is unbounded and the primal is infeasible.
124 raise ValueError("constraints are inconsistent, no solution")
125 step_length, full_step = t1, False
126 else:
127 t2 = abs(slack) / ztn
128 if t1inf or t1 >= t2:
129 step_length, full_step = t2, True
130 else:
131 step_length, full_step = t1, False
133 return (-step_length if reverse_step else step_length), full_step
136def _drop_constraint(idel: int, nact: int, uv: np.ndarray, iact: np.ndarray, J: np.ndarray, R: np.ndarray) -> int:
137 """Remove the ``idel``-th active constraint, closing the gap it leaves.
139 ``uv``, ``iact``, ``J`` and ``R`` are all modified in place.
141 Args:
142 idel: 1-based position in the active set of the constraint to drop.
143 nact: Size of the active set before the drop.
144 uv: Dual variables of the active constraints.
145 iact: 1-based indices of the active constraints.
146 J: ``(n, n)`` inverse Cholesky factor, updated by the QR downdate.
147 R: Packed upper triangular factor, updated by the QR downdate.
149 Returns:
150 The size of the active set after the drop.
151 """
152 qr_delete(nact, idel, J, R)
153 uv[idel - 1 : nact - 1] = uv[idel:nact].copy()
154 iact[idel - 1 : nact - 1] = iact[idel:nact].copy()
155 uv[nact - 1], iact[nact - 1] = 0.0, 0
156 return nact - 1
159def _dual_step_limit(
160 uv: np.ndarray,
161 rv: np.ndarray,
162 iact: np.ndarray,
163 nact: int,
164 meq: int,
165 reverse_step: bool,
166) -> tuple[float, int]:
167 """Return the largest dual-feasible step and the constraint that limits it.
169 Stepping along ``-rv`` drives the multipliers of the active inequality
170 constraints towards zero. The first one to reach zero caps the step, since a
171 negative multiplier would be dual infeasible. Equality constraints are
172 exempt: their multipliers are unrestricted in sign.
174 Below :data:`_SCALAR_CUTOFF` active constraints the array form spends nearly
175 all of its time in NumPy's per-call overhead rather than on the handful of
176 divisions it performs, so the work is handed to :func:`_dual_step_scalar`,
177 which computes the same answer -- ties included -- on Python lists.
179 Args:
180 uv: Dual variables of the active constraints.
181 rv: Negated step direction of the dual variables.
182 iact: 1-based indices of the active constraints.
183 nact: Size of the active set.
184 meq: Number of leading constraints treated as equalities.
185 reverse_step: Whether the step is taken in the negative direction.
187 Returns:
188 The step limit and the 1-based position in the active set of the
189 constraint that attains it. The position is 0 when no constraint limits
190 the step, in which case the limit is meaningless.
191 """
192 if nact <= _SCALAR_CUTOFF:
193 return _dual_step_scalar(uv, rv, iact, nact, meq, reverse_step)
195 # Working with the signed direction lets one comparison serve both cases and
196 # makes the eligible entries positive, so no separate abs is needed.
197 direction = -rv[:nact] if reverse_step else rv[:nact]
198 eligible = (iact[:nact] > meq) & (direction > 0.0)
200 # `where` leaves the ineligible entries at the infinity they were filled
201 # with, so argmin skips them and the division never sees them.
202 ratio = np.full(nact, np.inf)
203 np.divide(uv[:nact], direction, out=ratio, where=eligible)
204 idel = int(np.argmin(ratio))
205 limit = float(ratio[idel])
206 return (0.0, 0) if limit == np.inf else (limit, idel + 1)
209def _dual_step_scalar(
210 uv: np.ndarray,
211 rv: np.ndarray,
212 iact: np.ndarray,
213 nact: int,
214 meq: int,
215 reverse_step: bool,
216) -> tuple[float, int]:
217 """Run the same ratio test on Python lists, for a small active set.
219 ``tolist()`` costs one call per array and then every comparison is an
220 interpreter operation rather than an array one, which wins outright while
221 the active set is smaller than :data:`_SCALAR_CUTOFF`. The list of active
222 indices is built only when there are equalities to exclude: with ``meq == 0``
223 every constraint is an inequality and that test is vacuous, which is the case
224 for box-constrained and dense-``C`` problems.
226 Args:
227 uv: Dual variables of the active constraints.
228 rv: Negated step direction of the dual variables.
229 iact: 1-based indices of the active constraints.
230 nact: Size of the active set.
231 meq: Number of leading constraints treated as equalities.
232 reverse_step: Whether the step is taken in the negative direction.
234 Returns:
235 Exactly what :func:`_dual_step_limit` returns, ties included.
236 """
237 sign = -1.0 if reverse_step else 1.0
238 directions = rv[:nact].tolist()
239 duals = uv[:nact].tolist()
240 active = iact[:nact].tolist() if meq else None
242 limit, idel = np.inf, 0
243 for i in range(nact):
244 direction = sign * directions[i]
245 if direction > 0.0 and (active is None or active[i] > meq):
246 ratio = duals[i] / direction
247 if ratio < limit:
248 limit, idel = ratio, i + 1
249 return (0.0, 0) if idel == 0 else (limit, idel)