Coverage for src/cvxball/_frame.py: 100%
85 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-09-07 19:33 +0000
1"""A support set whose factorisation is carried across iterations, not rebuilt.
3:func:`cvxball.min_circle_active_set` changes its support by one point per
4iteration, and every iteration needs two things from it: whether the support is
5affinely independent, and where the circumcentre of its face lies. Recomputing
6both from the raw points costs ``O(d r^2)`` for a support of ``r + 1`` points in
7``d`` dimensions. Repairing a factorisation that already exists costs
8``O(d r)``.
10At the dimensions this module exists for, that is the whole running time. It is
11also the entirety of why SciPy is a dependency of this package: ``qr_insert``,
12``qr_delete`` and ``qr_update`` are the compiled Givens updates that make the
13repair possible. Written as a Python loop they would lose to one vectorised
14``numpy.linalg.qr``, and the exercise would be pointless.
16**The algebra.** With ``D`` the ``r x d`` matrix of edges ``q_j - q_0`` and
17``D' = QR`` its economic factorisation, three identities take every per-iteration
18quantity off ``d`` and onto the small ``r x r`` block:
20- ``D D' = R'R``, so ``R`` is already a Cholesky factor of the Gram matrix the
21 circumcentre subproblem solves against: two triangular solves, no ``O(r^2 d)``
22 product.
23- ``null(D') = null(R)``, so the affine-dependence test reads ``R``. Note this is
24 the *right* null space; the left factor, which the edge matrix itself would
25 hand you, is the wrong one.
26- ``||q_j - q_0||^2 = ||R[:,j]||^2``, so even the right-hand side comes from ``R``.
28**Affine dependence is the awkward part**, and it is where this method differs
29from the pivoting method of Fischer, Gärtner and Kutz. Theirs forbids a dependent
30support outright; this one *permits* it and answers it with a descent direction
31in the null space, which means the factorisation has to survive situations an
32economic ``QR`` cannot represent: ``qr_insert`` refuses a column already in the
33span, and a dependent support can outgrow ``d + 1`` and leave ``R`` rectangular.
34Both are handled by refactorising, which ``numpy.linalg.qr`` accepts happily,
35returning the singular ``R`` that :meth:`_MaintainedFace.null_space` then reads.
36Such fallbacks are counted rather than hidden: if they stopped being rare the
37data structure would have stopped paying.
38"""
40import numpy as np
42# The three update routines are re-exported from a Cython extension, so the type
43# stubs do not declare them on `scipy.linalg` even though they are there at run
44# time. Importing them from the private module instead would be worse.
45from scipy.linalg import (
46 qr_delete, # ty: ignore[unresolved-import]
47 qr_insert, # ty: ignore[unresolved-import]
48 qr_update, # ty: ignore[unresolved-import]
49 solve_triangular,
50)
52_EPS = float(np.finfo(np.float64).eps)
55class _MaintainedFace:
56 """The support set, and the economic ``QR`` of its edge matrix, kept in step.
58 Three updates cover every move the solver makes. A point joining the support
59 appends a column, a point other than the origin leaving deletes one, and the
60 origin leaving has no column of its own -- every column being measured *from*
61 it -- so ``q_1`` is promoted and the remaining edges ``a_j - a_1`` come from a
62 deletion plus a rank-one update.
63 """
65 def __init__(self, points: np.ndarray, seed: int) -> None:
66 """Start from the one-point support ``{seed}``, whose edge matrix is empty.
68 Args:
69 points: The ``(n, d)`` cloud, already centred and scaled.
70 seed: Index of the point the support starts as.
71 """
72 self._points = points
73 self.support: list[int] = [seed]
74 self.fallbacks = 0
75 self._q: np.ndarray = np.zeros((points.shape[1], 0))
76 self._r: np.ndarray = np.zeros((0, 0))
78 @property
79 def face(self) -> np.ndarray:
80 """Return the ``(m, d)`` array of support points."""
81 return self._points[self.support]
83 def _refactorise(self) -> None:
84 """Rebuild ``Q`` and ``R`` from the current support."""
85 face = self.face
86 edges = (face[1:] - face[0]).T
87 if edges.shape[1] == 0:
88 self._q = np.zeros((edges.shape[0], 0))
89 self._r = np.zeros((0, 0))
90 else:
91 self._q, self._r = np.linalg.qr(edges)
93 def _economise(self) -> None:
94 """Trim ``Q`` and ``R`` back to the economic shape after an update.
96 SciPy's updates leave ``Q`` as wide as it was, so deleting a column from a
97 square factorisation returns a ``(d, d)`` ``Q`` beside a rectangular
98 ``R``. Left alone that breaks the triangular solve loudly and the
99 projection silently -- with ``Q`` square, ``Q Q'`` is the identity.
100 ``R`` is upper triangular, so the surplus is zero and dropping it leaves
101 ``QR = D'`` exactly.
102 """
103 columns = len(self.support) - 1
104 rows = min(self._points.shape[1], columns)
105 if self._q.shape[1] != rows or self._r.shape != (rows, columns):
106 self._q = self._q[:, :rows]
107 self._r = self._r[:rows, :columns]
109 def insert(self, index: int) -> None:
110 """Take ``points[index]`` into the support.
112 Args:
113 index: Index into the cloud of the entering point.
114 """
115 column = self._points[index] - self.face[0]
116 self.support.append(index)
117 if self._r.shape[0] == 0:
118 length = float(np.linalg.norm(column))
119 if length == 0.0:
120 self._refactorise()
121 return
122 self._q = (column / length)[:, None]
123 self._r = np.array([[length]])
124 return
125 try:
126 self._q, self._r = qr_insert(self._q, self._r, column, self._r.shape[1], which="col")
127 except np.linalg.LinAlgError:
128 # The entering point lies in the affine hull of the support, so the
129 # new column is already in the span of Q and no economic
130 # factorisation can hold it. Rebuilding always can, and the singular
131 # R that comes back is exactly what `null_space` needs to see.
132 self.fallbacks += 1
133 self._refactorise()
134 else:
135 self._economise()
137 def remove(self, position: int) -> None:
138 """Drop the support point at ``position``.
140 Args:
141 position: Index within the support list. Position 0 is the origin,
142 which takes the re-origining path.
143 """
144 if position > 0:
145 del self.support[position]
146 try:
147 self._q, self._r = qr_delete(self._q, self._r, position - 1, which="col")
148 except np.linalg.LinAlgError:
149 self.fallbacks += 1
150 self._refactorise()
151 else:
152 self._economise()
153 return
155 first_edge = self._points[self.support[1]] - self.face[0]
156 remaining = len(self.support) - 2
157 del self.support[0]
158 if remaining <= 0:
159 self._q = np.zeros((self._points.shape[1], 0))
160 self._r = np.zeros((0, 0))
161 return
162 try:
163 self._q, self._r = qr_delete(self._q, self._r, 0, which="col")
164 rows = min(self._points.shape[1], remaining)
165 self._q = self._q[:, :rows]
166 self._r = self._r[:rows, :remaining]
167 self._q, self._r = qr_update(self._q, self._r, -first_edge, np.ones(remaining))
168 except (np.linalg.LinAlgError, ValueError):
169 self.fallbacks += 1
170 self._refactorise()
171 else:
172 self._economise()
174 def null_space(self) -> np.ndarray:
175 """Find weight directions that reshuffle the support without moving the centre.
177 The same object :func:`cvxball.solver._affine_null_space` returns, read off
178 the ``r x r`` block instead of the ``r x d`` edge matrix: ``null(D')`` is
179 ``null(R)``, which is the *right* null space, so the factor to take is
180 ``Vh`` rather than ``U``.
182 Returns:
183 An ``(m, q)`` array of weight directions, empty exactly when the
184 support is affinely independent.
185 """
186 if self._r.shape[0] == 0:
187 return np.zeros((1, 0))
189 _, singular_values, right = np.linalg.svd(self._r, full_matrices=True)
190 cutoff = max(self._r.shape[1], self._points.shape[1]) * _EPS * float(singular_values[0])
191 rank = int(np.count_nonzero(singular_values > cutoff))
192 tail = right[rank:, :].T
193 return np.vstack([-tail.sum(axis=0, keepdims=True), tail])
195 def circumcentre_weights(self) -> np.ndarray:
196 """Solve the subproblem: put every support point on a common sphere.
198 ``D D' = R'R`` is already factored, so this is two triangular solves and
199 never touches ``d`` -- the right-hand side ``||q_j - q_0||^2 / 2`` is the
200 squared column norms of ``R``.
202 Returns:
203 The ``(m,)`` weights, summing to one, that express the circumcentre as
204 an affine combination of the support. A negative entry means the
205 circumcentre lies outside the simplex.
206 """
207 if self._r.shape[0] == 0:
208 return np.ones(1)
210 rhs = 0.5 * np.einsum("ij,ij->i", self._r.T, self._r.T)
211 forward = solve_triangular(self._r, rhs, lower=False, trans="T")
212 tail = solve_triangular(self._r, forward, lower=False)
213 return np.concatenate(([1.0 - float(tail.sum())], tail))