Coverage for src/cvx/quadprog/_qr.py: 100%
57 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-15 05:55 +0000
1"""Orthogonal updates of a QR factorisation held as an explicit pair (Q, R).
3The Goldfarb/Idnani algorithm maintains a factorisation of the matrix of active
4constraint normals. Every iteration either appends a column (a constraint enters
5the active set) or drops one (a constraint leaves), so the factorisation is
6updated rather than recomputed from scratch.
8``R`` is stored as packed columns, as in the reference C implementation: the
9entry in column ``j`` and row ``i`` (both 0-based, ``i <= j``) lives at flat
10offset ``j * (j + 1) // 2 + i``, so column ``j`` occupies one contiguous run of
11``j + 1`` values.
13That layout is not just about the factor of two in memory. The solver's hot
14operation is a triangular solve against the *leading* ``nact`` columns, and in
15this form that submatrix is the leading ``nact * (nact + 1) // 2`` entries --
16contiguous, so BLAS ``tpsv`` reads it in place. The same leading block of a dense
17``(r, r)`` array is strided, which forces a full copy on every call. The solve
18runs once per iteration, so the difference is paid every time.
20The gap is two effects rather than one, and `benchmarks/layout_probe.py`
21separates them. At ``nact = 800`` in an array twice as wide: 590 us strided,
2299 us once the array is Fortran-ordered so no copy is needed, and 31 us packed.
23Avoiding the copy is worth 5.9x, and the packed routine is worth a further 3.2x
24on top of that -- ``tpsv`` is a level-2 BLAS call reading its argument in place,
25where ``trtrs`` is a general LAPACK routine that no dense layout can talk out of
26being. 18.8x altogether, none of it arithmetic.
28The control that separates the two has a trap in it worth naming, because we fell
29in. It must be **Fortran**-ordered. ``np.ascontiguousarray`` gives C order, which
30LAPACK copies exactly as it copies a strided view, so a control built that way
31measures one copy against another, shows no difference, and appears to prove that
32the copy is not the cost.
34The price is paid in :func:`qr_delete`, which mixes two *rows* across a range of
35columns. Column strides grow with the column index, so that is a gather rather
36than a slice -- see the index arithmetic there.
38Choice of transformation
39------------------------
40The reference reduces an incoming column with a chain of Givens rotations, one
41per trailing component. :func:`qr_insert` instead applies a single Householder
42reflection, which is what the same reduction costs in one pair of BLAS calls
43rather than ``n - r`` Python-level ones.
45This is not merely an equivalent-cost rearrangement, and it is worth being
46precise about why it is legitimate. The two reductions produce *different* ``Q``
47and ``R``: signs along the diagonal of ``R``, and hence the signs of some columns
48of ``Q``, differ. The solver is nonetheless unaffected, because the quantities it
49actually consumes are invariant to that choice. With ``A`` the active normals and
50``G^-1 = J J^T``,
52 R^T R = A^T J J^T A = A^T G^-1 A and R^T d_1 = A^T G^-1 n
54so the dual step direction
56 rv = R^-1 d_1 = (R^T R)^-1 R^T d_1 = (A^T G^-1 A)^-1 A^T G^-1 n
58is a function of ``A``, ``n`` and ``G`` alone. Replacing ``R`` by ``S R`` for any
59sign matrix ``S`` also replaces ``d_1`` by ``S d_1``, and the two cancel exactly
60(a sign flip is exact in IEEE arithmetic). The primal direction
61``zv = J_2 d_2`` is invariant for the same reason: flipping the sign of a column
62of ``J`` flips the matching entry of ``d``, and their product is unchanged.
64:func:`qr_delete` keeps the Givens chase, which is inherently sequential: each
65rotation's parameters depend on the previous one having been applied.
66"""
68# Q and R are the names this factorisation has in every reference on the subject,
69# and the ones the reference implementation uses. Lowercasing them to satisfy
70# pep8-naming would obscure that. The exemption lives here rather than in a
71# [lint.per-file-ignores] block because ruff.toml is template-owned and a local
72# edit to it is reverted by the next `/rhiza:update` sync.
73# ruff: noqa: N803
75import math
77import numpy as np
78from scipy.linalg.blas import dger, drot
80__all__ = ["qr_delete", "qr_insert"]
82# `dger` and `drot` are named directly rather than resolved through
83# get_blas_funcs, for the reason given in _solve.py: every array reaching them is
84# already float64, so there is no precision left to select, and the wrappers are
85# the identical objects either way. Doing so also drops the cast that was
86# flattening their signatures to Callable[..., Any].
87#
88# dger is the rank-1 update, called rather than np.outer so that the update lands
89# in Q's own buffer instead of an O(n * k) temporary. drot is the plane rotation
90# _mix uses to apply the delete step's 2x2 to a pair of Q's columns without
91# allocating; see the note there on why a rotation suffices for a reflection.
94def qr_insert(r: int, av: np.ndarray, Q: np.ndarray, R: np.ndarray) -> None:
95 """Append ``av`` to ``R`` as its ``r``-th column, keeping ``R`` triangular.
97 An orthogonal transformation is applied to ``av`` to annihilate the
98 components beyond the ``r``-th, and the same transformation is applied to the
99 columns of ``Q``.
101 ``R`` is upper triangular of order ``r - 1`` on entry and of order ``r`` on
102 exit. All three arrays are modified in place.
104 Args:
105 r: 1-based size of the active set *after* the insertion.
106 av: Length-``n`` vector to append. Overwritten.
107 Q: ``(n, n)`` array whose columns receive the transformation. Should be
108 Fortran-ordered so that the column block is contiguous.
109 R: packed upper triangular array receiving the new column.
110 """
111 # Only columns r-1 .. n-1 take part. Columns of R already in place are
112 # untouched: their entries at those positions are exact zeros, and any
113 # combination of exact zeros is an exact zero.
114 av[r - 1] = _reflect(av[r - 1 :], Q[:, r - 1 :])
116 # Column r-1 holds rows 0 .. r-1, so it is r values at offset (r-1)r/2.
117 start = (r - 1) * r // 2
118 R[start : start + r] = av[:r]
121def _reflect(v: np.ndarray, block: np.ndarray) -> float:
122 """Reduce ``v`` to a multiple of ``e_1``, transforming ``block`` to match.
124 Applies the Householder reflection ``W = I - 2 u u^T / u^T u`` that maps ``v``
125 onto ``alpha e_1``, updating ``block`` in place as ``block @ W``. The sign of
126 ``alpha`` follows ``v[0]``, matching the reference implementation's Givens
127 chain, and its leading entry is formed by a cancellation-free rearrangement
128 of ``v[0] - alpha``.
130 Args:
131 v: Vector to reduce. Not modified.
132 block: ``(n, len(v))`` column block to which the reflection is applied.
133 Modified in place.
135 Returns:
136 ``alpha``, the single surviving component of ``v``, carrying its sign.
137 """
138 head = float(v[0])
139 tail = v[1:]
140 tail_sq = float(tail @ tail)
142 if tail_sq == 0.0:
143 # Already a multiple of e_1, so the reflection is the identity. This also
144 # covers the len(v) == 1 case, where there is nothing to annihilate.
145 return head
147 norm = math.sqrt(head * head + tail_sq)
148 sign = 1.0 if head >= 0.0 else -1.0
149 alpha = sign * norm
151 # u = v - alpha e_1. Writing the leading entry as
152 # head - sign*norm = sign*(head^2 - norm^2)/(|head| + norm)
153 # avoids the cancellation that the direct subtraction suffers when v is
154 # already close to a positive multiple of e_1.
155 u = v.copy()
156 u[0] = -sign * tail_sq / (abs(head) + norm)
157 beta = tail_sq + u[0] * u[0]
159 # block @ (I - 2 u u^T / beta), as a matrix-vector product and a rank-1
160 # update -- two BLAS calls, independent of len(v).
161 w = block @ u
162 if block.flags.f_contiguous:
163 # dger writes through to block's buffer, which is a view into Q.
164 dger(-2.0 / beta, w, u, a=block, overwrite_a=True)
165 else:
166 # A non-Fortran-ordered block would be copied by dger, losing the
167 # update, so fall back to an explicit (allocating) rank-1 update.
168 block -= np.outer(w, u * (2.0 / beta))
169 return alpha
172def qr_delete(r: int, col: int, Q: np.ndarray, R: np.ndarray) -> None:
173 """Drop the ``col``-th column of ``R``, restoring upper triangular form.
175 Orthogonal transformations are applied to the rows of ``R`` to bring it back
176 to upper triangular form, and the same transformations are applied to the
177 columns of ``Q``.
179 ``R`` is upper triangular of order ``r`` on entry and of order ``r - 1`` on
180 exit. Entries outside the leading ``(r - 1, r - 1)`` block are left stale;
181 the caller shrinks the active set accordingly, so they are never read before
182 being overwritten by :func:`qr_insert`.
184 Args:
185 r: 1-based size of the active set *before* the deletion.
186 col: 1-based index of the column of ``R`` to drop.
187 Q: ``(n, n)`` array whose columns receive the transformations.
188 R: packed upper triangular array to be updated.
189 """
190 for i in range(col, r):
191 # On this iteration, reduce the (i, i) element of R to zero,
192 # then move column i to position i - 1.
193 diagonal = i * (i + 1) // 2 + i
194 if R[diagonal] == 0.0: # pragma: no cover
195 # Defensive, and unreachable in exact arithmetic: a diagonal entry
196 # vanishes only if the active constraint normals lose independence,
197 # which the solver's step rules prevent. Kept because the reference
198 # has it, and cancellation could in principle produce a true zero.
199 continue
201 # The transformation mixes rows i - 1 and i of R over columns i .. r - 1.
202 # Consecutive rows of one column are adjacent, but the column offsets
203 # grow, so addressing a row across columns needs explicit indices.
204 columns = np.arange(i, r)
205 lower = columns * (columns + 1) // 2 + i
206 upper = lower - 1
208 if R[diagonal - 1] == 0.0:
209 # Nothing to combine, so the reflection degenerates to a swap.
210 # Fancy indexing reads copies, so the two assignments cannot alias.
211 R[upper], R[lower] = R[lower], R[upper]
212 _swap(Q[:, i - 1], Q[:, i])
213 else:
214 gc, gs = _reflection_2x2(R[diagonal - 1], R[diagonal])
215 # Rows of R and columns of Q take the same 2x2 reflection: it is
216 # symmetric, so transforming J's columns by it transforms J^T's rows
217 # by it too, which is what keeps J^T A == [[R], [0]] intact.
218 first, second = R[upper], R[lower]
219 R[upper] = gc * first + gs * second
220 R[lower] = gs * first - gc * second
221 _mix(Q[:, i - 1], Q[:, i], gc, gs)
223 # Move column i left into slot i - 1, keeping its rows 0 .. i-1. The
224 # entry just zeroed is dropped with the slot it vacates.
225 R[(i - 1) * i // 2 : (i - 1) * i // 2 + i] = R[diagonal - i : diagonal]
228def _reflection_2x2(x: float, y: float) -> tuple[float, float]:
229 """Return the reflection ``[[c, s], [s, -c]]`` that annihilates ``y``.
231 The result is symmetric and orthogonal, so it is its own inverse. The sign of
232 the hypotenuse follows ``x``, so the surviving component keeps its sign.
234 Args:
235 x: Component to be preserved.
236 y: Component to be annihilated. Must be nonzero.
238 Returns:
239 The cosine and sine of the reflection.
240 """
241 h = math.hypot(x, y)
242 if x < 0.0:
243 h = -h
244 return x / h, y / h
247def _mix(first: np.ndarray, second: np.ndarray, gc: float, gs: float) -> None:
248 """Apply ``[[gc, gs], [gs, -gc]]`` to a pair of vectors in place.
250 BLAS offers no 2x2 reflection, only ``rot``'s rotation
251 ``[[gc, -gs], [gs, gc]]``. The two agree on the first output, and the
252 rotation's second output is the exact negation of the reflection's, so one
253 sign flip -- exact in IEEE -- recovers the reflection. That is worth doing
254 because spelling the arithmetic in NumPy allocates two ``n``-vectors per
255 call, and this runs once per step of the chase in :func:`qr_delete`.
257 Args:
258 first: First vector, overwritten with ``gc * first + gs * second``.
259 second: Second vector, overwritten with ``gs * first - gc * second``.
260 gc: Cosine of the reflection.
261 gs: Sine of the reflection.
262 """
263 if first.dtype == np.float64 and first.flags.contiguous and second.flags.contiguous:
264 # drot writes through to the columns of Q, which are contiguous when Q is
265 # Fortran-ordered as the solver builds it. On a strided view f2py copies
266 # instead and silently drops the overwrite, hence the guard -- the same
267 # hazard _reflect handles for dger.
268 drot(first, second, gc, gs, overwrite_x=True, overwrite_y=True)
269 second *= -1.0
270 else:
271 combined = gc * first + gs * second
272 second *= -gc
273 second += gs * first
274 first[:] = combined
277def _swap(first: np.ndarray, second: np.ndarray) -> None:
278 """Exchange the contents of two vectors in place.
280 Args:
281 first: First vector.
282 second: Second vector.
283 """
284 tmp = first.copy()
285 first[:] = second
286 second[:] = tmp