7
$\begingroup$

I have a list of $d$-dimensional data points $x_i$, and I want to find a $d$-ball which minimizes the sum of the squares of distances between the surface of the ball and the points. In other words, I’m trying to find a center point $c$ and a radius $r$ that minimizes

$$\sum_{i} (r-\text{dist}(c,x_i))^2$$

I thought about computing the partial derivatives of the above function and solving them simultaneously but I feel like there must be a more direct method in Mathematica, using FindMinimum or a similar function. What’s the best way to accomplish this?

$\endgroup$
  • $\begingroup$ It sounds like you're need general methodology advice more-so than anything hyper-specific to programming in Mathematica. $\endgroup$ – user6014 Dec 31 '19 at 20:04
  • 1
    $\begingroup$ @user6014 Perhaps, but I already know how to solve this problem - in any other language, I would compute the derivatives by hand and just have the program numerically solve them. I’m wondering whether the optimization built-ins of Mathematica offer an advantage in this particular case $\endgroup$ – Nico A Dec 31 '19 at 20:07
  • $\begingroup$ Optimization methods will outperform solving for partial derivatives of a sum of squares set to zero. If the points tend to be in any sense randomly arranged, good initial values for center and radius might be the barycenter (mean) and the mean of distances to that point. $\endgroup$ – Daniel Lichtblau Jan 1 at 19:28
6
$\begingroup$

kglr has provided an answer if you want the points to be enclosed within the sphere. However, that was not specified in your question, so I propose a direct approach that minimizes the distance between the points and the surface of a sphere, as you mentioned.

Note that this could result in a large sphere whose surface is locally "almost flat" where your points are, thereby minimizing the distance between that locally-almost-planar-surface and the points.

SeedRandom[20191231]
pts = RandomReal[{3, 5}, {20, 3}];

ClearAll[distfun]
distfun[centerCoords_?(VectorQ[#, NumericQ] &), radius_?NumericQ, points_] := 
 Module[{rdf},
   rdf = RegionDistance@Sphere[centerCoords, radius];
   Total[rdf@points]
 ]

min = NMinimize[{distfun[Array[x, 3], r, pts], r > 0}, Flatten@{Array[x, 3], r}]

Graphics3D[{
    Red, PointSize[0.02], Point[pts],
    White, Opacity[0.5], Sphere[Array[x, 3], r] /. Last@min
  }, Lighting -> "Neutral"
]

best fit sphere and points in 3D

$\endgroup$
  • $\begingroup$ This is exactly what I was looking for - but is there a solution which works in more than 3 dimensions, without manually writing in (in my case) 784 variables? $\endgroup$ – Nico A Dec 31 '19 at 20:30
  • $\begingroup$ This works in any dimensions. Specify the variables using Array instead of naming them specifically (e.g. Array[x, 3] will give {x[1], x[2], x[3]}). That many variables will probably be slow though. $\endgroup$ – MarcoB Dec 31 '19 at 20:38
  • 1
    $\begingroup$ Great! Didn’t know you could use Array like that. Thank you! $\endgroup$ – Nico A Dec 31 '19 at 20:39
  • 1
    $\begingroup$ +1 Isn't Norm[rdf@points] rather than Total[rdf@points] the correct way to minimize "the sum of the squares of distances between the surface of the ball and the points" $\endgroup$ – Bob Hanlon Dec 31 '19 at 21:37
  • 1
    $\begingroup$ @mikado I think it can be written as a sum of squares of norms. Those last are convex by triangle inequality, and the compositions of squaring and summing maintain convexity. $\endgroup$ – Daniel Lichtblau Jan 2 at 15:04
7
$\begingroup$

You can use BoundingRegion with "MinBall" as the form spec:

SeedRandom[1];
pts = RandomReal[10, {100, 10}];

BoundingRegion[pts, "MinBall"]

Ball[{5.03195, 4.49308, 5.1499, 4.36746, 5.62561, 5.27036, 5.95754, 4.78651, 4.46375, 5.04961}, 10.8419]

$\endgroup$
5
$\begingroup$

Here, as a supplement, my simple straightforward approach, which doesn't need the Region- functionality

SeedRandom[20191231]
pts = RandomReal[{3, 5}, {20, 3}];

opt = NMinimize[{Total@Map[ (Norm[(# - {x1, x2, x3})] - r )^2 &, pts], r > 0}, {r,x1,x2, x3}]
(*{25.4341, {r -> 1.91772, x1 -> 3.9304, x2 -> 3.94589, x3 -> 3.97011}}*)

Graphics3D[{Opacity[.1] , Sphere[{x1, x2, x3}, r], Opacity[1], Red,Point[{x1, x2, x3}], Blue, Point[pts]} /. opt[[2]], Lighting -> "Neutral"]

enter image description here

$\endgroup$
  • $\begingroup$ @DanielLichtblau Sorry, you are right I modify my answer $\endgroup$ – Ulrich Neumann Jan 1 at 21:05
2
$\begingroup$

Ulrich beat me to it, but I'm showing this because I use it on a regular basis to identify True Airspeed at altitude from the sequence of GPS groundspeed/track readings at a constant indicated airspeed, where you have the deal with error from winds at altitude. Once you've converted groundspeed/track to rectilinear coordinates, the radius is the true airspeed.

Penalty function

penalty[c_, r_, ptList_] := Total[(Map[(r - Norm[# - c])^2 &, ptList])];

Exercising it, just need to make sure you get the number of dimensions right.

dim = 2;
center = Array[c,dim];
rectList = {{70., 0.}, {0., 80.}, {-90., 0.}}
Minimize[penalty[center, vtas, rectList], Flatten@{center, vtas}]

(* {3.55269*10^-20, {c[1] -> -10., c[2] -> 0.625, vtas -> 80.0024}} *)

Example with a bunch of points:

enter image description here

$\endgroup$
2
$\begingroup$

Stated differently, the best-fitting hypersphere for the point set is sought. This can be recast as a linear least-squares problem. I show one way to go about that in this 2003 MathGroup post. Here I'll just give the code.

hypersphereFit[data_List] := 
 Module[{mat, rhs, qq, rr, params, cen}, 
  mat = Map[Append[#, 1] &, data];
  rhs = Map[#.# &, data];
  {qq, rr} = QRDecomposition[mat];
  params = LinearSolve[rr, qq.rhs];
  cen = Drop[params, -1]/2;
  {Sqrt[Last[params] + cen.cen], cen}]

Here is an example.

pts = RandomReal[{0, 10}, {20, 3}];
hypersphereFit[pts]

(* Out[773]= {5.01619, {5.6742, 5.47101, 4.33734}} *)

We check this by redoing using nonlinear optimization. First find plausible starting points.

c = Mean[pts];
dists = Sqrt[(c - #).(c - #)] & /@ pts;
r = Mean[dists];

Now form a sum of squares of discrepancies and minimize.

newdists = Sqrt[({cx, cy, cz} - #).({cx, cy, cz} - #)] & /@ pts;
Timing[fm = 
  FindMinimum[
   Total[(newdists^2 - rad^2)^2], {cx, c[[1]]}, {cy, c[[2]]}, {cz, 
    c[[3]]}, {rad, r}]]

(* Out[776]= {0.015625, {1147.63,
{cx -> 5.6742, cy -> 5.47101, 
   cz -> 4.33734, rad -> 5.01619}}} *)

We recovered the same values. Note however that the problem solved here is not quite the same as the one originally stated. Here we square distances and take discrepancies between those and the squared radius. So this minimum is with respect to a slightly different metric.

$\endgroup$
  • $\begingroup$ I really like the linear algebra approach! (+1) $\endgroup$ – MarcoB Jan 2 at 5:40
1
$\begingroup$

Unless I have misunderstood, the problem is not guaranteed to have a unique minimum. I give an example for which this occurs. Take the following set of 4 points in 2D

example = {{3, 4}, {0, 10}, {-3, 4}, {0, -5}}
(* {{3, 4}, {0, 10}, {-3, 4}, {0, -5}} *)

Define the objective function specified by the OP

obj[r_, c_, x_] := Total[(r - Norm[# - c])^2 & /@ x]

Specify the minimum as a function of r

minval[r_?NumericQ] := First[NMinimize[Evaluate[obj[r, {cx, cy}, example]], {cx, cy}]]

We can now see that this objective function has local minima for 2 different values of r. (One is probably at r -> ∞ in this case).

Plot[minval[r], {r, 3, 20}, PlotRange -> All]

enter image description here

Clearly, the objective function specified in the OP is not convex.

Starting the minimisation at different values of r gives different solutions.

FindMinimum[minval[r], {r, 0}]

(* During evaluation of FindMinimum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances.*)
(* {14.8752, {r -> 5.754}} *)

FindMinimum[minval[r], {r, 20}]

(* During evaluation of FindMinimum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances.*)
(* {18.0015, {r -> 664.526}} *)
$\endgroup$
  • $\begingroup$ I see your point about failing convexity. But for practical use, that might be more than anything else an argument for using the different formulation. $\endgroup$ – Daniel Lichtblau Jan 2 at 21:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.