Coverage for src/cvxball/fischer_gaertner_kutz.py: 100%
197 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-09-07 19:33 +0000
1"""The Fischer-Gärtner-Kutz pivoting method, the second solver this package ships.
3Kaspar Fischer, Bernd Gärtner and Martin Kutz, *Fast Smallest-Enclosing-Ball
4Computation in High Dimensions*, ESA 2003, LNCS 2832, 630-641.
6:func:`min_circle_fgk` answers the same question as
7:func:`cvxball.min_circle_active_set`, to the same exactness, and reaches it from
8the opposite side. The active-set method works on the *dual* -- a QP over the unit
9simplex, whose iterates are weights -- and its ball encloses the cloud only once
10the KKT test passes. This method is *primal-feasible throughout*: every iterate is
11an honest enclosing ball, and the algorithm deflates it. Hopp and Reeve proposed
12the same picture as a heuristic; the contribution of Fischer et al. is the proof
13of termination, and the pivot rule that makes it hold under degeneracy.
15**Which one to call.** Either: they agree on the ball and on the support set, and
16on Gaussian clouds from ``d = 1000`` to ``d = 16000`` they are within a factor of
171.1 to 1.6 in time (``experiments/bench_seb.py``). The active-set
18method is the default because it is the faster of the two on every row measured
19and returns the dual weights, which are a certificate the caller can check in one
20pass. Call this one when a *feasible* ball matters before convergence -- its
21iterates enclose the cloud and its radius falls monotonically, where the
22active-set radius rises to the answer from below and its ball encloses nothing
23until the last iteration -- or when the support set and pivot counts of
24:func:`ball_with_counts` are what you are after.
26The state is a pair ``(c, T)`` carrying the invariant (Fig. 2 of the paper)::
28 B(c, T) contains S, T is on the boundary of B(c, T), T affinely independent
30where ``B(c, T)`` is the ball about ``c`` through the farthest point of ``T``.
31By Lemma 1 (Seidel) the ball is optimal exactly when ``c = cc(T)``, the
32circumcentre of ``T``, *and* ``c`` lies in the convex hull of ``T`` -- so the
33algorithm loops until ``c`` is in ``conv(T)``, each iteration being a *dropping*
34phase (only when ``c`` has landed in ``aff(T)``) followed by a *walking* phase:
36- **Dropping.** ``c`` in ``aff(T)`` but not in ``conv(T)`` means some affine
37 coefficient of ``c`` with respect to ``T`` is negative. Drop such a point.
38 Lemmata 4 and 5 together say the dropped point cannot immediately stop the walk
39 that follows, which is what stops the two phases undoing each other.
40- **Walking.** Move ``c`` along the straight line towards ``cc(T)``. Lemma 3 says
41 that segment is orthogonal to ``aff(T)``, that ``T`` stays on the boundary all
42 the way along it, and that the radius strictly decreases -- the deflation. Stop
43 early at the first point of ``S`` to hit the shrinking boundary, and take it
44 into ``T``; otherwise arrive at ``cc(T)``.
46Two things follow from Lemma 3 that this module leans on, both verified
47numerically before they were relied upon:
49- Because ``[c, cc(T)]`` is orthogonal to ``aff(T)``, the circumcentre *is* the
50 orthogonal projection of ``c`` onto ``aff(T)`` (paper, section 4). So no
51 circumsphere is ever solved for -- one projection does the job, which is the
52 reformulation the paper's own implementation uses.
53- The stopping fraction has a closed form. Walking ``c(a) = c + a e`` with
54 ``e = cc(T) - c``, a point ``p`` joins the boundary when ``||p - c(a)||`` equals
55 the common support distance; the ``a^2`` terms cancel and, using
56 ``<q - c, e> = ||e||^2`` for every ``q`` in ``T``, that leaves
58 a_p = (R^2 - ||p - c||^2) / (2 (||e||^2 - <p - c, e>)).
60 Lemma 4's criterion (1) -- ``p`` is "behind" ``aff(T)`` and so uncritical -- is
61 exactly the statement that this denominator is non-positive, so one sign test
62 separates the candidates and the numerator is then non-negative by the
63 invariant. The walk takes ``a = min(1, min_p a_p)``.
65**Scope.** Both halves of the paper are here: the pivoting algorithm of Fig. 2,
66and section 4's dynamic QR-decomposition that makes it fast. :class:`_Frame`
67carries ``Q`` and ``R`` for the edge matrix ``A = [q_1 - q_0, ...]`` across pivots
68and repairs them in ``O(d r)`` as points enter and leave, rather than
69refactorising in ``O(d r^2)``; ``dynamic_qr=False`` selects the rebuild instead,
70which gives the same answers and is the baseline that says what the data
71structure is worth.
73Implementing it at all is only worthwhile because ``scipy.linalg`` exposes the
74updates -- ``qr_insert``, ``qr_delete``, ``qr_update`` -- as compiled
75LAPACK-backed routines. Written as a Python loop the Givens sweeps would lose to a
76single vectorised ``numpy.linalg.qr``, and a measurement of this method against
77the shipped one would then be a measurement of the interpreter. Those three
78routines are already why SciPy is a dependency of this package, for
79:mod:`cvxball._frame`, so this module adds no import that ``import cvxball`` did
80not already pull in.
82Both of the paper's pivot rules are here, selected by ``pivot_rule``:
84- ``"bland"`` -- Bland's rule adapted to this setting: fix an order on ``S``, drop
85 the negative-coefficient point of smallest rank, and admit the smallest-rank
86 point when several stop the walk at once. Theorem 1 proves termination with it,
87 degeneracies included.
88- ``"heuristic"`` (the default, and what the paper's own code runs) -- drop the
89 point of *minimal* coefficient, and among points that stop the walk at
90 effectively the same place admit the one farthest from ``aff(T)``. The paper
91 reports Bland's rule as correct but slow, and this as both faster and more
92 robust to roundoff.
94Like the shipped solver, and unlike the paper's fixed epsilon, this is written to
95be scale- and origin-invariant: the cloud is recentred on its mean before the
96first iteration, and the one tolerance carrying units -- the margin a point must
97clear from ``aff(T)`` before it is allowed into the support, the paper's stability
98threshold -- is sized off the cloud's extent. The coefficient tolerance needs no
99such treatment, being dimensionless and bounded by Lemma 2 (every coefficient of a
100centre in ``conv(T)`` is at most 1/2).
101"""
103from typing import Literal, NamedTuple
105import numpy as np
107# The three update routines are re-exported from a Cython extension, so the type
108# stubs do not declare them on `scipy.linalg` even though they are there at run
109# time -- the same suppression `cvxball._frame` needs, for the same three names.
110from scipy.linalg import (
111 qr_delete, # ty: ignore[unresolved-import]
112 qr_insert, # ty: ignore[unresolved-import]
113 qr_update, # ty: ignore[unresolved-import]
114 solve_triangular,
115)
117from cvxball.solver import _validate
119# How far a point must sit from aff(T), relative to the cloud's extent, before it
120# may enter the support. This is the paper's stability threshold (section 4): it
121# is what keeps T affinely independent in floating point, since a point *on*
122# aff(T) would make the edge matrix rank-deficient and the next projection
123# meaningless. Points behind aff(T) are discarded anyway, by Lemma 4.
124_AFF_RTOL = 1e-13
125# Affine coefficients are dimensionless and, by Lemma 2, at most 1/2 in the
126# optimal configuration -- so this needs no scaling, unlike the margin above.
127_COEFF_TOL = 1e-13
128# The smallest angle, as a reciprocal condition number, at which a point may join
129# the support. `qr_insert` refuses anything below roughly machine epsilon here;
130# this sits far enough above that the update never has to be second-guessed.
131_RCOND_FLOOR = 1e-12
132# Two stopping points count as tied, and so as competing under the pivot rule,
133# when their stopping fractions agree to this much.
134_TIE_RTOL = 1e-10
135# Safety net only. Termination is a theorem under Bland's rule, so reaching this
136# means numerical trouble rather than a slow instance.
137_MAX_ITER_PER_POINT = 50
139PivotRule = Literal["heuristic", "bland"]
142class Ball(NamedTuple):
143 """A ball, plus the work that went into finding it.
145 Attributes:
146 radius: The radius of the enclosing ball.
147 centre: The centre, of shape ``(d,)``.
148 support: Indices into the input of the final support set ``T``, whose
149 circumsphere is the ball. At most ``d + 1`` of them.
150 iterations: How many turns of the main loop ran.
151 drops: How many points left the support set.
152 insertions: How many points entered it.
153 """
155 radius: float
156 centre: np.ndarray
157 support: np.ndarray
158 iterations: int
159 drops: int
160 insertions: int
163class _Frame:
164 """The support set, and the economic ``QR`` of its edge matrix, kept in step.
166 This is section 4 of the paper: rather than refactorise the edge matrix
167 ``A = [q_1 - q_0, ..., q_r - q_0]`` from scratch at every pivot, carry ``Q``
168 and ``R`` and repair them as one point enters or leaves. Rebuilding costs
169 ``O(d r^2)``; each repair costs ``O(d r)``, and with ``r`` in the hundreds at
170 the dimensions this method is built for, that is the difference the paper's
171 engineering is about.
173 Three updates cover every move the algorithm makes, and only the third needs
174 thought:
176 - **A point joins ``T``.** One column is appended: ``qr_insert``.
177 - **A point other than the origin leaves.** One column is dropped:
178 ``qr_delete``.
179 - **The origin leaves.** No single column corresponds to ``q_0``, since every
180 column is measured *from* it, so dropping it changes them all. Promote
181 ``q_1`` to origin and the new edges are ``a_j - a_1``: delete column ``a_1``,
182 then subtract it from what remains, which is a rank-one update. Two compiled
183 calls, both ``O(d r)``, which is the "appropriate rank-1-update" the paper
184 mentions without spelling out.
186 All three come from ``scipy.linalg``, so the Givens sweeps run compiled. That
187 is what makes reproducing section 4 worthwhile here: written as a Python loop
188 they would lose to a single vectorised ``numpy.linalg.qr``, and the comparison
189 would measure the interpreter rather than the data structure. Pass
190 ``dynamic=False`` to get exactly that rebuild-every-pivot behaviour, which is
191 what the two are measured against each other with.
192 """
194 def __init__(self, points: np.ndarray, origin: int, dynamic: bool) -> None:
195 """Start from the one-point support ``{origin}``, whose edge matrix is empty.
197 Args:
198 points: The ``(n, d)`` cloud, in centred coordinates.
199 origin: Index of the single point the support starts as.
200 dynamic: Maintain ``Q`` and ``R`` across changes, rather than
201 refactorising after each one.
202 """
203 self._points = points
204 self._dynamic = dynamic
205 self.support = [origin]
206 self.rebuilds = 0
207 # Updates that could not be applied and were refactorised instead. The
208 # stability threshold makes a dependent insert rare but does not make it
209 # impossible: it bounds the entering point's distance from aff(T) from
210 # below, which is a weaker statement than scipy's condition-number test.
211 self.fallbacks = 0
212 self._q: np.ndarray = np.zeros((points.shape[1], 0))
213 self._r: np.ndarray = np.zeros((0, 0))
215 @property
216 def origin(self) -> np.ndarray:
217 """Return ``q_0``, the point every edge is measured from."""
218 point: np.ndarray = self._points[self.support[0]]
219 return point
221 @property
222 def basis(self) -> np.ndarray:
223 """Return ``Q``: orthonormal columns spanning the direction space of ``aff(T)``."""
224 return self._q
226 def _refactorise(self) -> None:
227 """Rebuild ``Q`` and ``R`` from the current support, from scratch."""
228 self.rebuilds += 1
229 face = self._points[self.support]
230 edges = (face[1:] - face[0]).T
231 if edges.shape[1] == 0:
232 self._q = np.zeros((edges.shape[0], 0))
233 self._r = np.zeros((0, 0))
234 else:
235 self._q, self._r = np.linalg.qr(edges)
237 def _economise(self) -> None:
238 """Trim ``Q`` and ``R`` back to the economic shape after an update.
240 scipy's updates leave ``Q`` as wide as it was, so deleting a column from a
241 square factorisation returns ``Q`` of shape ``(d, d)`` beside an ``R`` of
242 shape ``(d, r-1)`` -- no longer the economic pair. Both consumers break on
243 that, and one of them breaks *silently*: with ``Q`` square, ``Q Q'`` is the
244 identity, so the projection onto ``aff(T)`` would come back as the point
245 itself and the walk direction would collapse to zero, reading as "already
246 on the hull" at every subsequent pivot.
248 ``R`` is upper triangular, so its surplus rows are zero and the surplus
249 columns of ``Q`` multiply nothing: dropping them leaves ``QR = A`` exactly.
250 """
251 columns = len(self.support) - 1
252 if self._q.shape[1] != columns or self._r.shape != (columns, columns):
253 self._q = self._q[:, :columns]
254 self._r = self._r[:columns, :columns]
256 def insert(self, index: int) -> None:
257 """Take ``points[index]`` into the support, appending one edge column.
259 Args:
260 index: Index into the cloud of the entering point.
261 """
262 column = self._points[index] - self.origin
263 self.support.append(index)
264 if not self._dynamic:
265 self._refactorise()
266 elif self._r.shape[0] == 0:
267 # The first edge: its economic QR is just the normalised column, and
268 # scipy has no zero-column factorisation to insert into.
269 length = float(np.linalg.norm(column))
270 self._q = (column / length)[:, None]
271 self._r = np.array([[length]])
272 else:
273 try:
274 self._q, self._r = qr_insert(self._q, self._r, column, self._r.shape[1], which="col")
275 except np.linalg.LinAlgError:
276 # The column is already in the span of Q, so no economic
277 # factorisation can hold it. Rebuilding always can: numpy's QR
278 # accepts a rank-deficient matrix and simply returns a singular R.
279 self.fallbacks += 1
280 self._refactorise()
281 else:
282 self._economise()
284 def remove(self, position: int) -> None:
285 """Drop the support point at ``position``, repairing the factorisation.
287 Args:
288 position: Index *within the support list*, not into the cloud.
289 Position 0 is the origin, and takes the re-origining path.
290 """
291 if not self._dynamic:
292 del self.support[position]
293 self._refactorise()
294 return
296 if position > 0:
297 # An ordinary column, sitting at position - 1 of the edge matrix.
298 del self.support[position]
299 try:
300 self._q, self._r = qr_delete(self._q, self._r, position - 1, which="col")
301 except np.linalg.LinAlgError:
302 self.fallbacks += 1
303 self._refactorise()
304 else:
305 self._economise()
306 return
308 # The origin leaves, so q_1 is promoted and every edge is re-measured from
309 # it. Deleting a_1's column leaves [a_2, ..., a_r]; the rank-one update
310 # then turns those into [a_2 - a_1, ..., a_r - a_1].
311 first_edge = self._points[self.support[1]] - self.origin
312 remaining = self._r.shape[0] - 1
313 if remaining == 0:
314 self._q = np.zeros((self._points.shape[1], 0))
315 self._r = np.zeros((0, 0))
316 del self.support[0]
317 else:
318 del self.support[0]
319 try:
320 self._q, self._r = qr_delete(self._q, self._r, 0, which="col")
321 self._q = self._q[:, :remaining]
322 self._r = self._r[:remaining, :remaining]
323 self._q, self._r = qr_update(self._q, self._r, -first_edge, np.ones(remaining))
324 except (np.linalg.LinAlgError, ValueError):
325 self.fallbacks += 1
326 self._refactorise()
327 else:
328 self._economise()
330 def admits(self, index: int) -> bool:
331 """Report whether ``points[index]`` can join the support safely.
333 The paper's stability threshold is absolute -- it asks that the entering
334 point sit some distance from ``aff(T)`` measured against the cloud's
335 extent. A factorisation cares about something else: the *angle*, i.e. that
336 distance relative to the entering column's own length, which is exactly
337 the reciprocal condition number ``qr_insert`` tests. The two disagree on a
338 cloud carrying a far outlier, where a point can clear the absolute margin
339 and still be numerically inside the span.
341 Testing what the factorisation tests is what keeps the two in step. It
342 matters more here than it would for the shipped solver, because this
343 algorithm has no answer to a dependent support -- Fig. 2's invariant
344 requires ``T`` affinely independent, and there is no null-space descent
345 step to fall back on.
347 Args:
348 index: Index into the cloud of the candidate point.
350 Returns:
351 ``True`` when the candidate is safely off ``aff(T)``.
352 """
353 offset = self._points[index] - self.origin
354 length = float(np.linalg.norm(offset))
355 if length == 0.0:
356 return False
357 residual = offset - self._q @ (self._q.T @ offset)
358 return bool(float(np.linalg.norm(residual)) / length > _RCOND_FLOOR)
360 def direction_to_circumcentre(self, centre: np.ndarray) -> np.ndarray:
361 """Return ``cc(T) - centre``, the walking direction.
363 By Lemma 3(i) that segment is orthogonal to ``aff(T)``, so the
364 circumcentre is the orthogonal projection of ``centre`` onto the hull and
365 ``Q Q'`` is all that is needed to find it.
367 Args:
368 centre: The current centre.
370 Returns:
371 The ``(d,)`` step from ``centre`` to the circumcentre of the support.
372 """
373 offset = centre - self.origin
374 step: np.ndarray = self._q @ (self._q.T @ offset) - offset
375 return step
377 def coefficients(self, centre: np.ndarray) -> np.ndarray:
378 """Express ``centre`` as an affine combination of the support.
380 The paper's route, verbatim: solve ``R x = Q' (centre - q_0)`` by back
381 substitution. The entries of ``x`` are the coefficients of ``q_1, ..., q_r``
382 and the missing one for ``q_0`` follows from their summing to one.
384 Meaningful only when ``centre`` lies in ``aff(T)``, which the caller
385 guarantees by asking solely after a walk has run to completion.
387 Args:
388 centre: The current centre.
390 Returns:
391 The ``(m,)`` coefficients. A negative entry certifies that ``centre``
392 lies outside ``conv(T)``, which is what drives the drop.
393 """
394 if self._r.shape[0] == 0:
395 return np.ones(1)
396 rhs = self._q.T @ (centre - self.origin)
397 try:
398 tail = solve_triangular(self._r, rhs, lower=False)
399 except np.linalg.LinAlgError:
400 # A support can drift into near-dependence despite the stability
401 # threshold, and a *rebuilt* factorisation then puts an exact zero on
402 # R's diagonal where the incrementally updated one keeps it merely
403 # small -- so this fires on `dynamic=False` and not on the maintained
404 # factorisation. The least-squares solution agrees with back
405 # substitution wherever the system is solvable at all, so falling back
406 # costs nothing and keeps the rebuild baseline usable as a comparison.
407 tail = np.linalg.lstsq(self._r, rhs, rcond=None)[0]
408 return np.concatenate(([1.0 - float(tail.sum())], tail))
411def _leaving_point(coefficients: np.ndarray, support: list[int], pivot_rule: PivotRule) -> int | None:
412 """Choose which support point to drop, or report that the ball is optimal.
414 Args:
415 coefficients: The affine coefficients of the centre with respect to the
416 support, as returned by :meth:`_Frame.coefficients`.
417 support: The indices currently in the support, in insertion order.
418 pivot_rule: ``"bland"`` takes the negative-coefficient point of smallest
419 rank in the fixed order on ``S`` -- here the index into the input,
420 which is the arbitrary order Theorem 1 asks us to fix.
421 ``"heuristic"`` takes the most negative coefficient.
423 Returns:
424 The *position within* ``support`` of the point to drop, or ``None`` when
425 every coefficient is non-negative -- the centre is then in ``conv(T)``,
426 which by Lemma 1 is the optimality certificate.
427 """
428 negative = np.flatnonzero(coefficients < -_COEFF_TOL)
429 if negative.size == 0:
430 return None
431 if pivot_rule == "bland":
432 return int(negative[np.argmin(np.asarray(support)[negative])])
433 return int(negative[np.argmin(coefficients[negative])])
436def _entering_point(
437 fractions: np.ndarray,
438 shortest: float,
439 points: np.ndarray,
440 frame: _Frame,
441 pivot_rule: PivotRule,
442) -> int:
443 """Choose which point stopping the walk should enter the support.
445 Args:
446 fractions: The ``(n,)`` stopping fractions, ``inf`` where a point cannot
447 stop the walk at all.
448 shortest: The smallest of them, the distance actually walked.
449 points: The ``(n, d)`` cloud, in centred coordinates.
450 frame: The current support and its factorisation.
451 pivot_rule: ``"bland"`` takes the smallest index among the points that
452 stop the walk at the same place, as Theorem 1 requires.
453 ``"heuristic"`` takes the one farthest from ``aff(T)``, the paper's
454 roundoff-motivated choice: the farther the new point is from the
455 hull it is joining, the better conditioned the enlarged support.
457 Returns:
458 The index into ``points`` of the entering point.
459 """
460 tied = np.flatnonzero(fractions <= shortest + _TIE_RTOL * max(abs(shortest), 1.0))
461 if tied.size == 1 or pivot_rule == "bland":
462 return int(tied[0])
464 basis = frame.basis
465 offsets = points[tied] - frame.origin
466 residuals = offsets - (offsets @ basis) @ basis.T
467 return int(tied[np.argmax(np.einsum("ij,ij->i", residuals, residuals))])
470def ball_with_counts(
471 points: np.ndarray,
472 pivot_rule: PivotRule = "heuristic",
473 dynamic_qr: bool = True,
474 verbose: bool = False,
475) -> Ball:
476 """Solve the smallest enclosing ball by the pivoting method, reporting the work done.
478 Args:
479 points: A ``(n, d)`` array with ``n >= 1``.
480 pivot_rule: ``"heuristic"`` (the paper's own code) or ``"bland"`` (the
481 rule Theorem 1 proves terminating).
482 dynamic_qr: Carry the factorisation across pivots, repairing it in
483 ``O(d r)`` as section 4 does. ``False`` refactorises from scratch at
484 every pivot, in ``O(d r^2)`` -- the same answers, and the baseline the
485 data structure is worth measuring against.
486 verbose: If ``True``, print the phase, support size and radius per turn.
488 Returns:
489 The :class:`Ball`, including the support set and the pivot counts.
491 Raises:
492 ValueError: If ``points`` is not a finite, non-empty ``(n, d)`` array (see
493 :func:`cvxball.solver._validate`), if ``pivot_rule`` is not one of the
494 two rules, or if the iteration limit is reached -- which, termination
495 being a theorem, means numerical trouble rather than a hard instance.
496 """
497 if pivot_rule not in ("heuristic", "bland"):
498 raise ValueError(f"pivot_rule must be 'heuristic' or 'bland', got {pivot_rule!r}") # noqa: TRY003
500 pts = _validate(points)
501 n = pts.shape[0]
503 # Recentre, so that every tolerance below is governed by the extent of the
504 # cloud rather than by its distance from an arbitrary origin, and difference
505 # before squaring for the same reason the shipped solver does.
506 shift = pts.mean(axis=0)
507 pts = pts - shift
508 extent = float(np.sqrt(np.einsum("ij,ij->i", pts, pts).max()))
509 if extent == 0.0:
510 # Every point is the same point; the ball is that point, radius zero.
511 return Ball(0.0, shift.copy(), np.zeros(1, dtype=np.intp), 0, 0, 0)
512 margin = _AFF_RTOL * extent
514 # Initialisation, from Fig. 2: c is any point of S, and T the single point of
515 # S farthest from it -- which is what makes B(c, T) enclose S to begin with.
516 centre = pts[0].copy()
517 offsets = pts - centre
518 frame = _Frame(pts, int(np.argmax(np.einsum("ij,ij->i", offsets, offsets))), dynamic_qr)
519 support = frame.support
520 # Whether c is known to lie in aff(T). It is not, initially: c is one point of
521 # the cloud and aff(T) is a different one. It becomes true exactly when a walk
522 # runs to completion, since the centre is then the projection onto that hull.
523 on_hull = False
525 drops = insertions = 0
526 limit = _MAX_ITER_PER_POINT * (n + 1)
527 for iteration in range(limit):
528 if on_hull:
529 leaving = _leaving_point(frame.coefficients(centre), support, pivot_rule)
530 if leaving is None:
531 # c is in conv(T): by Lemma 1 this ball is SEB(S).
532 radius = float(np.linalg.norm(pts[support] - centre, axis=1).max())
533 if verbose:
534 print(f"[{iteration:4d}] optimal support={len(support):3d} radius={radius:.12g}")
535 return Ball(radius, centre + shift, np.array(support, dtype=np.intp), iteration, drops, insertions)
536 if verbose:
537 print(f"[{iteration:4d}] drop support={len(support):3d} point={support[leaving]}")
538 frame.remove(leaving)
539 drops += 1
540 on_hull = False
541 continue
543 # --- Walking phase ---------------------------------------------------
544 # cc(T) is the orthogonal projection of c onto aff(T), by Lemma 3(i).
545 direction = frame.direction_to_circumcentre(centre)
546 direction_sq = float(direction @ direction)
547 if np.sqrt(direction_sq) <= margin:
548 # Already on the hull -- reached whenever the support has grown to
549 # d + 1 points, whose affine hull is the whole space.
550 on_hull = True
551 continue
553 offsets = pts - centre
554 squared = np.einsum("ij,ij->i", offsets, offsets)
555 radius_sq = float(squared[support].max())
557 # Lemma 4: p can stop the walk only if <p - c, e> < <e, e>, i.e. only if
558 # this denominator is positive. Requiring it to clear `margin * ||e||`
559 # rather than zero is the stability threshold, and it is exactly the
560 # right test: e is orthogonal to aff(T), so denom / ||e|| is p's distance
561 # from aff(T) measured along e, and thus a lower bound on its full
562 # distance from aff(T). Clearing the margin here therefore certifies that
563 # T stays affinely independent when p joins it -- and the points it
564 # excludes have a near-zero denominator, hence an enormous stopping
565 # fraction, so they were never going to be the minimiser anyway.
566 denominator = direction_sq - offsets @ direction
567 stoppers = denominator > margin * np.sqrt(direction_sq)
568 stoppers[support] = False
570 fractions = np.full(n, np.inf)
571 np.divide(
572 np.maximum(radius_sq - squared, 0.0),
573 2.0 * denominator,
574 out=fractions,
575 where=stoppers,
576 )
577 # Walk to the nearest stopper the support can actually take. A candidate
578 # that fails `admits` is numerically inside aff(T): taking it would break
579 # the affine independence Fig. 2's invariant rests on, so it is passed
580 # over and the next-nearest considered. Usually the first one is fine.
581 while True:
582 shortest = float(fractions.min())
583 if shortest >= 1.0:
584 # Nothing stops the walk: the centre reaches cc(T) and lands on the hull.
585 centre = centre + direction
586 on_hull = True
587 if verbose:
588 radius = float(np.linalg.norm(pts[support] - centre, axis=1).max())
589 print(f"[{iteration:4d}] arrive support={len(support):3d} radius={radius:.12g}")
590 break
592 entering = _entering_point(fractions, shortest, pts, frame, pivot_rule)
593 if not frame.admits(entering):
594 fractions[entering] = np.inf
595 continue
597 centre = centre + shortest * direction
598 frame.insert(entering)
599 insertions += 1
600 if verbose:
601 radius = float(np.linalg.norm(pts[support] - centre, axis=1).max())
602 print(f"[{iteration:4d}] insert support={len(support):3d} radius={radius:.12g} step={shortest:.6g}")
603 break
605 raise ValueError(f"pivoting method did not converge in {limit} iterations") # noqa: TRY003
608def min_circle_fgk(points: np.ndarray, verbose: bool = False) -> tuple[float, np.ndarray]:
609 """Compute the smallest enclosing ball, in this package's solver signature.
611 Args:
612 points: A numpy array of shape ``(n, d)`` where *n* is the number of
613 points and *d* is the ambient dimension.
614 verbose: If ``True``, print one line per pivot step. Defaults to ``False``.
616 Returns:
617 A tuple ``(radius, center)``, matching :func:`cvxball.min_circle_active_set`.
619 Raises:
620 ValueError: If ``points`` is not a finite, non-empty ``(n, d)`` array.
622 Example:
623 The right triangle whose smallest enclosing circle is the one on its
624 hypotenuse.
626 The values are rounded here, as they are for the cone program, and the
627 reason is worth stating because it is easy to assume otherwise: this
628 method terminates at an exact *combinatorial* configuration -- the
629 support set ``{(1, 0), (0, 1)}`` -- but the centre it reports is not the
630 exact circumcentre of that set. It is the running sum of the walks that
631 got there, so it lands a few ulp out (four, on this input). The shipped
632 active-set method solves afresh for the centre of its final support and
633 so returns ``sqrt(2) / 2`` bit-for-bit; the difference is one of
634 arithmetic, not of which ball the two methods identify.
636 >>> import numpy as np
637 >>> from cvxball import min_circle_fgk
638 >>> radius, center = min_circle_fgk(np.array([[0, 0], [1, 0], [0, 1]]))
639 >>> round(radius, 12)
640 0.707106781187
641 >>> np.round(center, 12)
642 array([0.5, 0.5])
643 """
644 ball = ball_with_counts(points, verbose=verbose)
645 if verbose:
646 print(f"fgk: iterations={ball.iterations} drops={ball.drops} insertions={ball.insertions}")
647 return ball.radius, ball.centre