Coverage for src/cvx/quadprog/_qr.py: 100%
57 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-06 11:28 +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: at ``n =
18700`` that measured 77 us per solve against 7.5 us for ``tpsv``, and the solve
19runs once per iteration.
21The price is paid in :func:`qr_delete`, which mixes two *rows* across a range of
22columns. Column strides grow with the column index, so that is a gather rather
23than a slice -- see the index arithmetic there.
25Choice of transformation
26------------------------
27The reference reduces an incoming column with a chain of Givens rotations, one
28per trailing component. :func:`qr_insert` instead applies a single Householder
29reflection, which is what the same reduction costs in one pair of BLAS calls
30rather than ``n - r`` Python-level ones.
32This is not merely an equivalent-cost rearrangement, and it is worth being
33precise about why it is legitimate. The two reductions produce *different* ``Q``
34and ``R``: signs along the diagonal of ``R``, and hence the signs of some columns
35of ``Q``, differ. The solver is nonetheless unaffected, because the quantities it
36actually consumes are invariant to that choice. With ``A`` the active normals and
37``G^-1 = J J^T``,
39 R^T R = A^T J J^T A = A^T G^-1 A and R^T d_1 = A^T G^-1 n
41so the dual step direction
43 rv = R^-1 d_1 = (R^T R)^-1 R^T d_1 = (A^T G^-1 A)^-1 A^T G^-1 n
45is a function of ``A``, ``n`` and ``G`` alone. Replacing ``R`` by ``S R`` for any
46sign matrix ``S`` also replaces ``d_1`` by ``S d_1``, and the two cancel exactly
47(a sign flip is exact in IEEE arithmetic). The primal direction
48``zv = J_2 d_2`` is invariant for the same reason: flipping the sign of a column
49of ``J`` flips the matching entry of ``d``, and their product is unchanged.
51:func:`qr_delete` keeps the Givens chase, which is inherently sequential: each
52rotation's parameters depend on the previous one having been applied.
53"""
55# Q and R are the names this factorisation has in every reference on the subject,
56# and the ones the reference implementation uses. Lowercasing them to satisfy
57# pep8-naming would obscure that. The exemption lives here rather than in a
58# [lint.per-file-ignores] block because ruff.toml is template-owned and a local
59# edit to it is reverted by the next `/rhiza:update` sync.
60# ruff: noqa: N803
62import math
63from collections.abc import Callable
64from typing import Any, cast
66import numpy as np
67import scipy.linalg
69__all__ = ["qr_delete", "qr_insert"]
71# Rank-1 update, resolved once. Called directly rather than through np.outer so
72# that the update lands in Q's own buffer instead of an O(n * k) temporary. The
73# cast records that asking for a single name yields a single function, which
74# scipy-stubs cannot express -- see the matching note in _solve.py.
75_GER = cast("Callable[..., Any]", scipy.linalg.get_blas_funcs("ger", (np.empty(0, dtype=np.float64),)))
78def qr_insert(r: int, av: np.ndarray, Q: np.ndarray, R: np.ndarray) -> None:
79 """Append ``av`` to ``R`` as its ``r``-th column, keeping ``R`` triangular.
81 An orthogonal transformation is applied to ``av`` to annihilate the
82 components beyond the ``r``-th, and the same transformation is applied to the
83 columns of ``Q``.
85 ``R`` is upper triangular of order ``r - 1`` on entry and of order ``r`` on
86 exit. All three arrays are modified in place.
88 Args:
89 r: 1-based size of the active set *after* the insertion.
90 av: Length-``n`` vector to append. Overwritten.
91 Q: ``(n, n)`` array whose columns receive the transformation. Should be
92 Fortran-ordered so that the column block is contiguous.
93 R: packed upper triangular array receiving the new column.
94 """
95 # Only columns r-1 .. n-1 take part. Columns of R already in place are
96 # untouched: their entries at those positions are exact zeros, and any
97 # combination of exact zeros is an exact zero.
98 av[r - 1] = _reflect(av[r - 1 :], Q[:, r - 1 :])
100 # Column r-1 holds rows 0 .. r-1, so it is r values at offset (r-1)r/2.
101 start = (r - 1) * r // 2
102 R[start : start + r] = av[:r]
105def _reflect(v: np.ndarray, block: np.ndarray) -> float:
106 """Reduce ``v`` to a multiple of ``e_1``, transforming ``block`` to match.
108 Applies the Householder reflection ``W = I - 2 u u^T / u^T u`` that maps ``v``
109 onto ``alpha e_1``, updating ``block`` in place as ``block @ W``. The sign of
110 ``alpha`` follows ``v[0]``, matching the reference implementation's Givens
111 chain, and its leading entry is formed by a cancellation-free rearrangement
112 of ``v[0] - alpha``.
114 Args:
115 v: Vector to reduce. Not modified.
116 block: ``(n, len(v))`` column block to which the reflection is applied.
117 Modified in place.
119 Returns:
120 ``alpha``, the single surviving component of ``v``, carrying its sign.
121 """
122 head = float(v[0])
123 tail = v[1:]
124 tail_sq = float(tail @ tail)
126 if tail_sq == 0.0:
127 # Already a multiple of e_1, so the reflection is the identity. This also
128 # covers the len(v) == 1 case, where there is nothing to annihilate.
129 return head
131 norm = math.sqrt(head * head + tail_sq)
132 sign = 1.0 if head >= 0.0 else -1.0
133 alpha = sign * norm
135 # u = v - alpha e_1. Writing the leading entry as
136 # head - sign*norm = sign*(head^2 - norm^2)/(|head| + norm)
137 # avoids the cancellation that the direct subtraction suffers when v is
138 # already close to a positive multiple of e_1.
139 u = v.copy()
140 u[0] = -sign * tail_sq / (abs(head) + norm)
141 beta = tail_sq + u[0] * u[0]
143 # block @ (I - 2 u u^T / beta), as a matrix-vector product and a rank-1
144 # update -- two BLAS calls, independent of len(v).
145 w = block @ u
146 if block.flags.f_contiguous:
147 # dger writes through to block's buffer, which is a view into Q.
148 _GER(-2.0 / beta, w, u, a=block, overwrite_a=True)
149 else:
150 # A non-Fortran-ordered block would be copied by dger, losing the
151 # update, so fall back to an explicit (allocating) rank-1 update.
152 block -= np.outer(w, u * (2.0 / beta))
153 return alpha
156def qr_delete(r: int, col: int, Q: np.ndarray, R: np.ndarray) -> None:
157 """Drop the ``col``-th column of ``R``, restoring upper triangular form.
159 Orthogonal transformations are applied to the rows of ``R`` to bring it back
160 to upper triangular form, and the same transformations are applied to the
161 columns of ``Q``.
163 ``R`` is upper triangular of order ``r`` on entry and of order ``r - 1`` on
164 exit. Entries outside the leading ``(r - 1, r - 1)`` block are left stale;
165 the caller shrinks the active set accordingly, so they are never read before
166 being overwritten by :func:`qr_insert`.
168 Args:
169 r: 1-based size of the active set *before* the deletion.
170 col: 1-based index of the column of ``R`` to drop.
171 Q: ``(n, n)`` array whose columns receive the transformations.
172 R: packed upper triangular array to be updated.
173 """
174 for i in range(col, r):
175 # On this iteration, reduce the (i, i) element of R to zero,
176 # then move column i to position i - 1.
177 diagonal = i * (i + 1) // 2 + i
178 if R[diagonal] == 0.0: # pragma: no cover
179 # Defensive, and unreachable in exact arithmetic: a diagonal entry
180 # vanishes only if the active constraint normals lose independence,
181 # which the solver's step rules prevent. Kept because the reference
182 # has it, and cancellation could in principle produce a true zero.
183 continue
185 # The transformation mixes rows i - 1 and i of R over columns i .. r - 1.
186 # Consecutive rows of one column are adjacent, but the column offsets
187 # grow, so addressing a row across columns needs explicit indices.
188 columns = np.arange(i, r)
189 lower = columns * (columns + 1) // 2 + i
190 upper = lower - 1
192 if R[diagonal - 1] == 0.0:
193 # Nothing to combine, so the reflection degenerates to a swap.
194 # Fancy indexing reads copies, so the two assignments cannot alias.
195 R[upper], R[lower] = R[lower], R[upper]
196 _swap(Q[:, i - 1], Q[:, i])
197 else:
198 gc, gs = _reflection_2x2(R[diagonal - 1], R[diagonal])
199 # Rows of R and columns of Q take the same 2x2 reflection: it is
200 # symmetric, so transforming J's columns by it transforms J^T's rows
201 # by it too, which is what keeps J^T A == [[R], [0]] intact.
202 first, second = R[upper], R[lower]
203 R[upper] = gc * first + gs * second
204 R[lower] = gs * first - gc * second
205 _mix(Q[:, i - 1], Q[:, i], gc, gs)
207 # Move column i left into slot i - 1, keeping its rows 0 .. i-1. The
208 # entry just zeroed is dropped with the slot it vacates.
209 R[(i - 1) * i // 2 : (i - 1) * i // 2 + i] = R[diagonal - i : diagonal]
212def _reflection_2x2(x: float, y: float) -> tuple[float, float]:
213 """Return the reflection ``[[c, s], [s, -c]]`` that annihilates ``y``.
215 The result is symmetric and orthogonal, so it is its own inverse. The sign of
216 the hypotenuse follows ``x``, so the surviving component keeps its sign.
218 Args:
219 x: Component to be preserved.
220 y: Component to be annihilated. Must be nonzero.
222 Returns:
223 The cosine and sine of the reflection.
224 """
225 h = math.hypot(x, y)
226 if x < 0.0:
227 h = -h
228 return x / h, y / h
231def _mix(first: np.ndarray, second: np.ndarray, gc: float, gs: float) -> None:
232 """Apply ``[[gc, gs], [gs, -gc]]`` to a pair of vectors in place.
234 Args:
235 first: First vector, overwritten with ``gc * first + gs * second``.
236 second: Second vector, overwritten with ``gs * first - gc * second``.
237 gc: Cosine of the reflection.
238 gs: Sine of the reflection.
239 """
240 combined = gc * first + gs * second
241 second *= -gc
242 second += gs * first
243 first[:] = combined
246def _swap(first: np.ndarray, second: np.ndarray) -> None:
247 """Exchange the contents of two vectors in place.
249 Args:
250 first: First vector.
251 second: Second vector.
252 """
253 tmp = first.copy()
254 first[:] = second
255 second[:] = tmp