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

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.

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)

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}}]

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

IANATopologist, and I don't know how clear I'm making myself in this, but for what it's worth, these are some thoughts on the matter:

Using mass and gravity to calculate this sort of thing might indeed be elegant -- ISTR that there are a number of libraries and efficient algorithms to find the gravity vectors for any number of masses.

If you were using a spherical map, I'd suggest finding within the sphere the actual center of gravity for your N mass points. You then draw a line from the center outward through this inner center of gravity to find the point on the sphere's surface where your mass points wish to congregate.

However, a toroidal map makes this difficult.

My suggestion, then, is to flatten and copy your map to give you a 3 x 3 quilt of maps (using an infinite field of maps will give better results, but might be overkill). I'll assign coordinates (0, 0) to (2, 2) to them, with (1, 1) being your source map. Find the point(s) to which the mass points of your inner map (1, 1) are attracted -- if they all go towards the middle of your map, fine: you found your center of gravity. If not, if one of the points close to the edge is going towards some mass accumulation outside of your inner map, say into map (2, 1), then discard this mass point when calculating your center of gravity. Instead you use the mass point from the opposite map ((0, 1) in this case) that wants to wander over into your middle map.

Adding the acceleration vectors for these mass points gives you the center of gravity on your torus. Done.