"Center of Mass" between a set of points on a Toroidally-Wrapped Map that minimizes average distance to all points

Question

edit As someone has pointed out, what I'm looking for is actually the point minimizing total geodesic distance between all other points

---

My map is topographically similar to the ones in Pac Man and Asteroids. Going past the top will warp you to the bottom, and going past the left will warp you to the right.

Say I have two points (of the same mass) on the map and I wanted to find their center of mass. I could use the classical definition, which basically is the midpoint.

However, let's say the two points are on opposite ends of the mass. There is another center of mass, so to speak, formed by wrapping "around". Basically, it is the point equidistant to both other points, but linked by "wrapping around" the edge.

Example

b . O . . a . . O .

Two points O. Their "classical" midpoint/center of mass is the point marked a. However, another midpoint is also at b (b is equidistant to both points, by wrapping around).

In my situation, I want to pick the one that has lower average distance between the two points. In this case, a has an average distance between the two points of three steps. b has an average distance of two steps. So I would pick b.

One way to solve for the two-point situation is to simply test both the classical midpoint and the shortest wrapped-around midpoint, and use the one that has a shorter average distance.

However! This does not easily generalize to 3 points, or 4, or 5, or n points.

Is there a formula or algorithm that I could use to find this?

(Assume that all points will always be of equal mass. I only use "center of mass" because it is the only term I knew to loosely describe what I was trying to do)

If my explanation is unclear, I will try to explain it better.

Answer 1

The notion of center of mass is a notion relevant on affine spaces. The n-dimensional torus has no affine structure.

What you want is a point which minimizes (geodesic) distance to all the other points.

I suggest the following: let x_1...x_n be a collection of points on the d-dimensional torus (or any other metric space for that purpose).

Your problem:

find a point mu such that sum(dist(mu, x_k)^2) is minimal.

In the affine-euclidian case, you get the usual notion of center of mass back.

This is a problem you will be able to solve (for instance, there are probably better options) with the conjugate gradient algorithm, which performs well in this case. Beware that you need moderate n (say n < 10^3) since the algorithm needs n^2 in space and n^3 in time.

Perhaps better suited is the Levenberg-Marquardt algorithm, which is tailored for minimization of sum of squares.

Note that if you have a good initial guess (eg. the usual center of mass of the points seen as points in R^d instead of the torus) the method will converge faster.

Edit: If (x1...xd) and (y1...yd) are points on the torus, the distance is given by dist(x, y)^2 = alpha1^2 + ... + alphad^2

where alphai = min((xi - yi) mod 1, (yi - xi) mod 1)

Answer 2

I made a little program to check the goodness of the involved functions and found that you should be very carefull with the minimization process.

Below you can see two sets of plots showing the points distribution, the function to minimize in the euclidean case, and the one corresponding to the "toric metric".

As you may see, the euclidean distance is very well-behaved, while the toric present several local minima that difficult the finding of the global minima. Also, the global minimum in the toric case is not unique.

Just in case, the program in Mathematica is:

Clear["Global`*"];

(*Define  non wrapping distance for dimension n*)
nwd[p1_, p2_, n_] := (p1[[n]] - p2[[n]])^2;

(*Define wrapping distance for dimension n *)
wd[p1_, p2_, max_,n_] := (max[[n]] - Max[p1[[n]], p2[[n]]] + Min[p1[[n]], p2[[n]]])^2;

(*Define minimal distance*)
dist[p1_, p2_, max_] :=
  Min[nwd[p1, p2, 1], wd[p1, p2, max, 1]] +
  Min[nwd[p1, p2, 2], wd[p1, p2, max, 2]];

(*Define Euclidean distance*)
euclDist[p1_, p2_, max_] := nwd[p1, p2, 1] + nwd[p1, p2, 2];

(*Set torus dimensions *)
MaxX = 20; 
MaxY = 15;

(*Examples of Points sets *)
lCircle = 
  Table[{10 Cos[fi] + 10, 5 Sin[fi] + 10}, {fi, 0, 2 Pi - .0001, Pi/20}];

lRect = Join[
   Table[{3, y}, {y, MaxY - 1}],
   Table[{MaxX - 1, y}, {y, MaxY - 1}],
   Table[{x, MaxY/2}, {x, MaxY - 1}],
   Table[{x, MaxY - 1}, {x, MaxX - 1}],
   Table[{x, 1}, {x, MaxX - 1}]];

(*Find Euclidean Center of mass *)
feucl = FindMinimum[{Total[
    euclDist[#, {a, b}, {MaxX, MaxY}] & /@ lRect], 0 <= a <= MaxX, 
             0 <= b <= MaxY}, {{a, 10}, {b, 10}}]

(*Find Toric Center of mass *)
ftoric = FindMinimum[{Total[dist[#, {a, b}, {MaxX, MaxY}] & /@ lRect],
         0 <= a <= MaxX, 0 <= b <= MaxY}, {{a, 10}, {b, 10}}]

Answer 3

In the 1 dimensional case, your problem would be analagous to finding a mean angle. The mean of angles a and b can be computed by

mean = remainder( a + remainder( b-a, C)/2.0, C) where C is the measure of a whole circle (ie 2*PI if you're using radians).

If you have n angles a[], the mean can be computed by

mean = a[0]; for i=1..n mean=remainder( mean + remainder( a[i]-mean, C)/(i+1), C)

So I reckon

meanX = X[0]; meanY = Y[0]

for i=1..n

 meanX = remainder( meanX + remainder( X[i]-meanX, W)/(i+1), W)
 meanY = remainder( meanY + remainder( Y[i]-meanY, H)/(i+1), H)

might do the job.

But note that this will result in -W/2<=meanX