Singular matrix error for points in non-generic position

[OnDraganov]

The code assumes that points are in generic position in the following sense:
If the points are d-dimensional, then every p-sphere in R^d passes through at most p+1 points

If this is not the case, for example if we have four points in 3D lying on a common circle, then the difference vectors used in get_circumsphere are linearly dependent, and we get LinAlgError from numpy.linalg.solve.

There is a good use-case for non-generic point clouds: you can restrict the center of the minimized bounding ball to an affine space by reflecting the points through that affine space. For this reason, it would be useful if the code allowed for such point clouds.

Suggested solution:
Replace numpy.linalg.solve by numpy.linalg.lstsq.

The linear equation solved is of the form $U^T U x = b$, where we then care about $c = U x$. Therefore, if there is more solutions, it does not matter which one we take: if $x' = x + k$ is a different solution with $k \in \mathrm{Ker}\ U$, then $U x' = U (x + k) = U x = c$.

To treat the potential breaking case where the equation had no solutions and an approximation was returned, we can just check np.isclose(((a @ x - b) ** 2).sum(), 0.), where a, x, b are the matrix, the returned solution, and the right hand side.

Altogether I would replace line 48 with the following:

...
A = numpy.inner(U, U)
x, *_ = numpy.linalg.lstsq(A, B, rcond=None)
if not numpy.isclose(((A @ x - B) ** 2).sum(), 0):
    raise numpy.linalg.LinAlgError('Linear equation has no solution.')
C = numpy.dot(x, U)
...

[OnDraganov]

I found a better solution -- the code never runs into the singular matrix error if the circle containment is done right: instead of checking just for inequality, we can also check for isclose:

def circle_contains(D, p):
    c, r2 = D
    distance = numpy.square(p - c).sum()
    return distance <= r2 or np.isclose(distance, r2)

[marmakoide]

Can you provide an example of such a configuration of points ? I tried this, which is a degenerate configuration but does not trigger any issues

# Generate points on the unit circle, projected in 3d
n = 4
theta = numpy.linspace(0, 2 * numpy.pi, n, endpoint = False)
S = numpy.stack([
    numpy.cos(theta),
    numpy.sin(theta),
    numpy.zeros(n)
], axis = 1)

# Get the bounding sphere
C, r2 = miniball.get_bounding_ball(S)

# Check that the bounding sphere and the support sphere are correct
assert numpy.allclose(C, numpy.zeros(S.shape[1]))
assert numpy.allclose(r2, 1.)

[OnDraganov]

Sorry for a very long time it took me to react. I don't quite understand what the problem cases are, but below I give a few examples. Note that the error only raises ~10% of the times you run miniball.get_bounding_ball -- you need to be unlucky with the randomizer to put undesirable points together.

np.array([
[0.025728409071614222, 0.372152791836762, 0.3313301933484486], 
[0.2531163257212625, 0.9836838066044385, 0.2053900271644188], 
[0.01222939518218047, 0.8318221612473629, 0.6333422846064242], 
[0.15230829793039352, 0.45195218214902, 0.10645228460244449]
])

np.array([
[0.019703769768801283, 0.22959782389636796, 0.9347328105384431], 
[0.025728409071614222, 0.372152791836762, 0.3313301933484486], 
[0.07256394493759732, 0.18814654709738765, 0.6659776158852926], 
[-0.28448686909485404, 0.9845638879736157, 0.4729161039578733], 
[-0.29051150839766654, 0.8420089200332219, 1.0763187211478678], 
[-0.3373470442636498, 1.0260151647725961, 0.7416712986110239]
])

np.array([
[0.6230334535523369, 0.10461020463586657, 0.6246333540949156],
[0.4368180091809194, 0.47867605481000886, 0.4940056593092361],
[0.19956709139681328, 0.3156469527895308, 0.46603709980425945],
[0.19258810832074413, 0.00570894472836192, 0.77173157151973],
[0.07256394493759732, 0.18814654709738765, 0.6659776158852926],
[0.14079426383352756, 0.46309890200523474, 0.8713935666594159],
[0.3270097082049448, 0.08903305183109245, 1.0020212614450952],
[0.22415372406335649, -0.014249183710294266, 0.7522170325705094],
[0.2311327071394254, 0.29568882435087457, 0.4465225608550386],
[0.2954212758833918, 0.06960900045054658, 0.506756030291579]
])

[akchan]

I found a better solution -- the code never runs into the singular matrix error if the circle containment is done right: instead of checking just for inequality, we can also check for isclose:

def circle_contains(D, p):
    c, r2 = D
    distance = numpy.square(p - c).sum()
    return distance <= r2 or np.isclose(distance, r2)

I have same error LinAlgError: Singular Matrix at the line of numpy.linalg.solve() in the case of 3D points rarely.
The solution proposed by OnDraganov works effectively on my data. The error doesn't appear now.

I would be grateful if the revision request by OnDraganov is approved.

[akchan]

I got another error case. Here is the sample code.

import numpy
import miniball

S = numpy.array([[165.9375, 90.9375, 39.19999933],
 [169.6875, 87.1875, 39.19999933],
 [165.9375, 87.1875, 39.19999933],
 [169.6875, 90.9375, 39.19999933]])

miniball.get_circumsphere(S)