Coverage for src/cvxball/solver.py: 100%
120 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-09-07 19:33 +0000
1"""Convex utilities for computing the minimum enclosing circle/ball.
3:func:`min_circle_active_set` runs an active-set QP method on the *dual* of the
4enclosing-ball problem -- a QP over the unit simplex -- keeping a support set of
5points on the ball's boundary. Each iteration costs one small dense linear
6solve, and it terminates at an exact vertex of the dual feasible set rather than
7at an interior-point tolerance.
9It is the default of the two solvers this package ships, being the faster on
10every cloud measured and the one that returns a checkable dual;
11:mod:`cvxball.fischer_gaertner_kutz` is the other, and reaches the same answer
12from the primal-feasible side. NumPy carries the method; SciPy
13carries the factorisation of its support at large `d`, through the compiled
14Givens updates in :mod:`cvxball._frame` -- see :data:`_MAINTAIN_MIN_DIM` for
15where that starts to pay and why it is not used below it.
16The cone-program route that used to sit beside it -- assembling the
17second-order-cone program by hand and handing it to Clarabel -- now lives in
18``experiments/clarabel_ball.py``, because that is what it had become: the
19reference this method is measured against rather than a second way to get an
20answer. Moving it there is what lets Clarabel be a development dependency, so
21installing this package pulls in NumPy, SciPy, and nothing else.
22"""
24import numpy as np
26from cvxball._frame import _MaintainedFace
29def _validate(points: np.ndarray) -> np.ndarray:
30 """Check that ``points`` is a usable point cloud, and return it as floats.
32 Every rejection here is one a solver would otherwise take: a 1-D array used to
33 fail when ``points.shape`` was unpacked into two names, and an empty or
34 non-finite cloud used to reach the cone solver and come back as
35 ``DualInfeasible`` or ``NumericalError`` -- a status describing the program
36 rather than the input that produced it. Refusing the same cases with the
37 reason keeps the caller's error about the caller's data. The Clarabel route in
38 ``experiments/clarabel_ball.py`` imports this function for that reason.
40 Args:
41 points: The candidate array of shape ``(n, d)``.
43 Returns:
44 ``points`` as a float64 array, ready for the solver.
46 Raises:
47 ValueError: If ``points`` is not two-dimensional, holds no points, has
48 no coordinates, or contains a NaN or an infinity.
49 """
50 array = np.asarray(points, dtype=np.float64)
52 if array.ndim != 2:
53 raise ValueError(f"points must be a 2-D (n, d) array, got {array.ndim}-D of shape {array.shape}") # noqa: TRY003
54 if array.shape[0] == 0:
55 raise ValueError("points is empty: the smallest enclosing ball of no points is undefined") # noqa: TRY003
56 if array.shape[1] == 0:
57 raise ValueError("points has no coordinates: shape (n, 0) describes no ambient space") # noqa: TRY003
58 if not np.isfinite(array).all():
59 bad = int(np.count_nonzero(~np.isfinite(array)))
60 raise ValueError(f"points must be finite: found {bad} NaN or infinite coordinate(s)") # noqa: TRY003
62 return array
65# --- Active-set (support-set) method ------------------------------------------
66#
67# The enclosing-ball problem has the concave dual
68#
69# maximise sum_i u_i ||p_i||^2 - ||sum_i u_i p_i||^2 over u >= 0, sum_i u_i = 1,
70#
71# a quadratic program over the unit simplex whose optimal centre is x = sum_i u_i p_i
72# and whose optimal value is the squared radius. Writing the convex negated dual as
73# g(u) = ||x||^2 - sum_i u_i ||p_i||^2, its KKT conditions read
74#
75# ||p_i - x|| == R for every i with u_i > 0 (support points sit on the boundary)
76# ||p_i - x|| <= R for every i with u_i == 0 (every other point is enclosed),
77#
78# which is precisely the geometric optimality certificate for the smallest enclosing
79# ball. The routines below run a primal active-set QP method on that dual: the free
80# set is the current support, every subproblem is a small dense linear system, and
81# each iteration either adds the farthest violating point or drops a support point
82# whose weight reaches zero.
84# Weights live on the unit simplex, so this is an absolute floor on a support weight.
85_DROP_TOL = 1e-12
86# Slack, relative to the squared radius, before a point counts as outside the ball.
87_FEAS_RTOL = 1e-9
88# Second, scale-aware slack, in multiples of the rounding error of a squared distance.
89# That error tracks the squared extent of the centred cloud, so a constant absolute
90# tolerance would carry units of length^2 and would silently declare optimality on
91# small-magnitude inputs.
92_FEAS_NOISE = 64.0
93# Safety net only: the method is finite, so hitting this means numerical trouble.
94_MAX_ITER_PER_POINT = 50
95# Ambient dimension from which carrying the factorisation across iterations starts
96# to pay. Below it the SciPy update calls cost more than the rebuild they save --
97# measured at 1000 standard normal points on an Apple M4 Pro, best of many runs:
98#
99# d rebuild maintained
100# 20 0.79 ms 1.05 ms 0.75x -- maintaining loses
101# 50 2.18 ms 2.52 ms 0.86x
102# 100 4.78 ms 4.60 ms 1.04x -- parity, near enough
103# 250 20.6 ms 16.2 ms 1.27x -- and it grows from here
104# 8000 11.3 s 3.19 s 3.54x
105#
106# So this is an empirical constant, not a derived one, and it is a threshold on
107# the *dimension* rather than on the support size because the choice has to be
108# made once, before the support exists. `maintain=` overrides it either way.
109_MAINTAIN_MIN_DIM = 100
112def _sq_dist(points: np.ndarray, centre: np.ndarray) -> np.ndarray:
113 """Compute squared distances from every row of ``points`` to ``centre``.
115 Differencing before squaring is what keeps this accurate. Expanding instead to
116 ``||p||^2 - 2 p'x + ||x||^2`` costs a cancellation between terms of size
117 ``||x||^2``, so for a cloud sitting far from the origin the digits that survive
118 are exactly the ones the answer needs — the part of size ``R^2``.
120 Args:
121 points: A ``(n, d)`` array of points.
122 centre: A ``(d,)`` centre.
124 Returns:
125 The ``(n,)`` array of squared distances.
126 """
127 offsets = points - centre
128 squared: np.ndarray = np.einsum("ij,ij->i", offsets, offsets)
129 return squared
132def _affine_null_space(face: np.ndarray) -> np.ndarray:
133 """Find weight directions that reshuffle the support without moving the centre.
135 A weight update ``p`` leaves both the simplex constraint and the centre
136 ``x = face.T @ u`` untouched exactly when ``sum(p) == 0`` and ``face.T @ p == 0``.
137 Eliminating ``p_0 = -(p_1 + ... + p_{m-1})`` collapses that pair of conditions to
138 the single condition ``edges.T @ p[1:] == 0`` on the edge matrix
139 ``edges[j] = q_j - q_0``, so such directions exist exactly when the face is
140 affinely *dependent*. Along one of them the quadratic term of ``g`` is frozen, so
141 ``g`` is *linear* and the equality-constrained subproblem is unbounded — the active
142 set then has to shrink instead of jumping to a minimiser.
144 Testing the rank on the edges rather than on ``[ones; face.T]`` is what keeps this
145 scale-invariant. With the row of ones stacked on, its unit entries dominate the
146 singular values, and every cloud whose extent falls below ``eps`` relative to 1
147 would be misjudged affinely dependent.
149 Args:
150 face: A ``(m, d)`` array of the points currently in the support set.
152 Returns:
153 An ``(m, q)`` array whose columns span those directions. It has ``q == 0``
154 columns exactly when the points of ``face`` are affinely independent — the
155 regular case, where the subproblem instead has the unique solution computed
156 by :func:`_face_weights`.
157 """
158 edges = face[1:] - face[0]
159 if edges.shape[0] == 0:
160 return np.zeros((1, 0))
162 # Only the left factor is read. When ``edges`` is no taller than it is wide the
163 # reduced decomposition already returns all ``m - 1`` of its columns -- a complete
164 # orthonormal basis of the edge row space *and* its null space -- so asking for the
165 # full one there buys nothing and costs the ``d x d`` right factor, built and then
166 # discarded. That discarded factor dominates everything at large ``d``: 2 GB per
167 # call at ``d = 16000``, and about forty times the cost of the whole solve.
168 #
169 # The guard is not decoration. Between a drop and the add that follows, the support
170 # can reach ``d + 2`` points, and on such a face the reduced left factor is
171 # ``(m - 1) x d`` -- one column short of spanning, and the column it is short of is
172 # exactly the null direction this function exists to find. Reporting that face as
173 # affinely independent would send the caller to :func:`_face_weights` with a
174 # singular system.
175 complete = edges.shape[0] > edges.shape[1]
176 left, singular_values, _ = np.linalg.svd(edges, full_matrices=complete)
177 cutoff = max(edges.shape) * float(np.finfo(np.float64).eps) * float(singular_values[0])
178 rank = int(np.count_nonzero(singular_values > cutoff))
180 # Lift each edge-space null vector z back to a weight direction [-sum(z), z].
181 tail = left[:, rank:]
182 return np.vstack([-tail.sum(axis=0, keepdims=True), tail])
185def _face_weights(face: np.ndarray) -> np.ndarray:
186 """Solve the subproblem: put every point of one face on a common sphere.
188 The centre is written as ``x = q_0 + D y``, where ``D`` holds the edges
189 ``q_j - q_0`` and so confines ``x`` to the affine hull of the face. Equating
190 the distances from ``x`` to ``q_0`` and to each ``q_j`` then collapses to the
191 tiny normal-equation system ``(D' D) y = h`` with ``h_j = ||q_j - q_0||^2 / 2``,
192 which is non-singular precisely because the face is affinely independent.
194 Args:
195 face: A ``(m, d)`` array of affinely independent points.
197 Returns:
198 The ``(m,)`` weights, summing to one, that express the circumcentre as an
199 affine combination of ``face``. A negative entry means the circumcentre
200 lies outside the simplex, so the caller has to drop a support point
201 instead of taking the full step.
202 """
203 edges = face[1:] - face[0]
204 if edges.shape[0] == 0:
205 return np.ones(1)
207 gram = edges @ edges.T
208 rhs = 0.5 * np.einsum("ij,ij->i", edges, edges)
209 y = np.linalg.solve(gram, rhs)
210 return np.concatenate(([1.0 - float(y.sum())], y))
213def _shrink(support: "_RebuiltFace | _MaintainedFace", keep: np.ndarray) -> None:
214 """Drop every support position whose weight has fallen to zero.
216 Positions are removed back to front so the earlier ones keep their indices,
217 which matters because a maintained factorisation is repaired per removal.
219 Args:
220 support: The face to shrink, in place.
221 keep: Boolean mask over the current support; ``False`` entries go.
222 """
223 for position in sorted(np.flatnonzero(~keep).tolist(), reverse=True):
224 support.remove(position)
227class _RebuiltFace:
228 """The support set, with its subproblem recomputed from the points each time.
230 The counterpart of :class:`cvxball._frame._MaintainedFace`, and the cheaper of
231 the two whenever the support is small: there is no factorisation to carry, so
232 nothing to repair, and the whole cost is two small dense decompositions that
233 NumPy dispatches straight into LAPACK. Below :data:`_MAINTAIN_MIN_DIM` that
234 beats paying SciPy's per-update overhead.
235 """
237 def __init__(self, points: np.ndarray, seed: int) -> None:
238 """Start from the one-point support ``{seed}``.
240 Args:
241 points: The ``(n, d)`` cloud, already centred and scaled.
242 seed: Index of the point the support starts as.
243 """
244 self._points = points
245 self.support: list[int] = [seed]
246 self.fallbacks = 0
248 @property
249 def face(self) -> np.ndarray:
250 """Return the ``(m, d)`` array of support points."""
251 return self._points[self.support]
253 def insert(self, index: int) -> None:
254 """Take ``points[index]`` into the support.
256 Args:
257 index: Index into the cloud of the entering point.
258 """
259 self.support.append(index)
261 def remove(self, position: int) -> None:
262 """Drop the support point at ``position``.
264 Args:
265 position: Index within the support list.
266 """
267 del self.support[position]
269 def null_space(self) -> np.ndarray:
270 """Return the weight directions that leave the centre fixed."""
271 return _affine_null_space(self.face)
273 def circumcentre_weights(self) -> np.ndarray:
274 """Return the barycentric weights of the face's circumcentre."""
275 return _face_weights(self.face)
278def min_circle_active_set(
279 points: np.ndarray,
280 verbose: bool = False,
281 maintain: bool | None = None,
282) -> tuple[float, np.ndarray]:
283 """Compute the smallest enclosing circle with an active-set method.
285 An active-set QP method on the dual of the enclosing-ball problem, in place of
286 handing a cone program to a conic solver. It
287 maintains a *support set* of points held on the ball's boundary and repeatedly
289 1. centres the ball on that support set by solving one small linear system
290 (:func:`_face_weights`),
291 2. shrinks the support when it cannot hold — either because a weight would turn
292 negative, or because the set has become affinely dependent
293 (:func:`_affine_null_space`) — moving as far as non-negativity allows, and
294 3. adds the farthest point that is still outside the ball.
296 Each subproblem is a ``k x k`` solve with ``k <= d``, so the cost per iteration
297 is driven by the dimension rather than by the number of points, and the method
298 stops at an exact vertex of the dual feasible set instead of at an
299 interior-point tolerance.
301 Args:
302 points: A numpy array of shape ``(n, d)`` where *n* is the number of
303 points and *d* is the ambient dimension.
304 verbose: If ``True``, print the support size and radius per iteration.
305 Defaults to ``False``.
306 maintain: Whether to carry the support's factorisation across iterations
307 and repair it, rather than rebuilding it each time. ``None``,
308 the default, decides on the ambient dimension: repairing costs
309 ``O(d r)`` against ``O(d r^2)`` to rebuild, which is decisive by
310 ``d = 250`` and a net loss below ``d = 100`` where SciPy's
311 per-update overhead exceeds what it saves (see
312 :data:`_MAINTAIN_MIN_DIM`). Both settings compute the same ball;
313 pass one explicitly only to measure the difference.
315 Returns:
316 A tuple ``(radius, center)`` where *radius* is the optimal enclosing
317 radius (float) and *center* is a numpy array of shape ``(d,)``.
319 Raises:
320 ValueError: If ``points`` is not a finite, non-empty ``(n, d)`` array
321 (see :func:`_validate`), or if the iteration limit is reached
322 without the optimality certificate holding.
324 Example:
325 Three points forming a right triangle, whose smallest enclosing circle is
326 the one on its hypotenuse. Both values are pinned to full precision here,
327 where the cone program in ``experiments/clarabel_ball.py`` can only pin its
328 centre to three decimals — this method stops at an exact vertex of the
329 dual feasible set, so on an input whose answer is exactly representable it
330 returns that answer bit-for-bit.
332 >>> import numpy as np
333 >>> from cvxball import min_circle_active_set
334 >>> points = np.array([[0, 0], [1, 0], [0, 1]])
335 >>> radius, center = min_circle_active_set(points)
336 >>> radius == 2**0.5 / 2
337 True
338 >>> center
339 array([0.5, 0.5])
340 """
341 pts = _validate(points)
342 n = pts.shape[0]
344 # Work in coordinates centred on the cloud, undoing the shift on the way out.
345 # Every quantity below -- the subproblem's Gram matrix, the squared distances, the
346 # rounding floor -- is then governed by the *extent* of the cloud rather than by
347 # its distance from an arbitrary origin. Without this, a cloud of extent 1 sitting
348 # a million units out gets a rounding floor of the same order as its own radius,
349 # and the optimality test below then accepts a ball that is visibly too small.
350 shift = pts.mean(axis=0)
351 pts = pts - shift
353 # Recentring fixes the origin but not the magnitude, and this method squares
354 # everything it touches: the Gram matrix, the squared distances, the noise floor.
355 # Squaring halves the usable exponent range, so a cloud of extent 1e-160 -- whose
356 # coordinates are ordinary doubles -- has a Gram matrix of order 1e-320, deep in
357 # the subnormals, and the solve for the circumcentre then overflows to infinity.
358 # The far end fails too, and worse: past 1e+154 the squares saturate and the
359 # method returns a radius of zero rather than raising.
360 #
361 # So normalise the extent to order one. The factor is a *power of two*, which
362 # makes the rescaling exact in binary floating point: not one bit of the answer
363 # moves, so the method keeps returning representable answers bit-for-bit (see the
364 # example above), and the whole iteration runs where doubles are well behaved.
365 # Measuring the extent as a largest coordinate rather than a largest norm is what
366 # keeps the measurement itself in range -- a norm would already have squared.
367 largest = float(np.abs(pts).max())
368 if largest == 0.0:
369 # Every point is the same point, so the ball is that point with radius zero.
370 return 0.0, shift
371 exponent = int(np.frexp(largest)[1])
372 pts = np.ldexp(pts, -exponent)
374 sq_norms: np.ndarray = np.einsum("ij,ij->i", pts, pts)
376 # Warm start with one point and whatever sits farthest from it: that puts both
377 # ends of a near-diameter into the support straight away, which is usually
378 # where the optimum keeps them.
379 seed = int(np.argmax(_sq_dist(pts, pts[0])))
380 # Rounding floor of a squared distance, set by the magnitude of the coordinates
381 # that go into it rather than by any fixed constant.
382 noise_floor = _FEAS_NOISE * float(np.finfo(np.float64).eps) * float(sq_norms.max())
384 if maintain is None:
385 maintain = pts.shape[1] >= _MAINTAIN_MIN_DIM
386 support = _MaintainedFace(pts, seed) if maintain else _RebuiltFace(pts, seed)
387 free = support.support
388 weights = np.ones(1)
390 iteration_limit = _MAX_ITER_PER_POINT * (n + 1)
391 for iteration in range(iteration_limit):
392 face = support.face
393 null_space = support.null_space()
395 if null_space.size:
396 # Affinely dependent support: g is linear on the null space, so follow
397 # steepest descent there until some weight is driven down to zero.
398 centre = face.T @ weights
399 # grad_i = ||x||^2 - ||p_i - x||^2, less the constant ||x||^2: null
400 # directions sum to zero, so dropping it leaves every projection alone
401 # while avoiding the cancellation of ||x||^2 against ||p_i||^2.
402 gradient = -_sq_dist(face, centre)
403 descent = -(null_space @ (null_space.T @ gradient))
404 # Rescale to a unit-max direction. The gradient carries units of length^2,
405 # so comparing it against the dimensionless weight tolerance below would
406 # make the whole method scale-dependent: on a cloud of extent 1e-20 no
407 # component would ever look binding. Only the direction matters here --
408 # `alpha` cancels any positive factor -- so normalising costs nothing.
409 largest = float(np.abs(descent).max())
410 step = descent / largest if largest > 0.0 else descent
411 else:
412 target = support.circumcentre_weights()
413 if target.min() >= -_DROP_TOL:
414 # The subproblem solution is feasible: take it, then test the KKT
415 # certificate and, if it fails, free the most violated point.
416 weights = np.maximum(target, 0.0)
417 weights /= weights.sum()
418 centre = face.T @ weights
419 dist_sq = _sq_dist(pts, centre)
420 radius_sq = float(dist_sq[free].max())
421 radius = float(np.sqrt(max(radius_sq, 0.0)))
422 if verbose:
423 print(f"[{iteration:4d}] support={len(free):3d} radius={radius:.12g}")
425 worst = int(np.argmax(dist_sq))
426 if dist_sq[worst] <= radius_sq * (1.0 + _FEAS_RTOL) + noise_floor:
427 # Undo the exact power-of-two normalisation, then the shift.
428 return float(np.ldexp(radius, exponent)), np.ldexp(centre, exponent) + shift
430 keep = weights > _DROP_TOL
431 _shrink(support, keep)
432 support.insert(worst)
433 weights = np.append(weights[keep], 0.0)
434 continue
436 step = target - weights
438 # Longest step along `step` that keeps every weight non-negative. A step
439 # of zero is impossible: the newly freed point always has step > 0, so the
440 # support strictly shrinks here and `g` strictly decreases.
441 binding = step < -_DROP_TOL
442 ratios = np.where(binding, -weights / np.where(binding, step, -1.0), np.inf)
443 alpha = float(ratios.min())
444 if not np.isfinite(alpha):
445 # Nothing blocks the step, so the support cannot shrink. Unreachable in
446 # exact arithmetic; bail out rather than propagate a non-finite weight.
447 raise ValueError("active-set method stalled: no support point blocks the step") # noqa: TRY003
448 weights = np.maximum(weights + alpha * step, 0.0)
450 keep = weights > _DROP_TOL
451 _shrink(support, keep)
452 weights = weights[keep]
453 weights /= weights.sum()
454 if verbose:
455 print(f"[{iteration:4d}] support={len(free):3d} drop step={alpha:.6g}")
457 raise ValueError(f"active-set method did not converge in {iteration_limit} iterations") # noqa: TRY003