GeometricMedian.frink

View or download GeometricMedian.frink in plain text format


/** This helps find the geometric median of a set of n-dimensional points.
    The geometric median is also called the spatial median, the
    Fermat-Torricelli-Weber point, etc.  The geometric median is the point
    minimizing the *sum* of the distances to the given points.

    As an example, if you wanted several people to meet by flying to a place
    that minimizes the total miles flown by everyone put together, you want the
    geometric median.

    This also allows the points to be "weighted", meaning that the cost for
    one person to travel can be made higher.

    This follows the algorithm described in:

    Vardi, Y, and C H Zhang. “The multivariate L1-median and associated data
    depth.” Proceedings of the National Academy of Sciences of the United
    States of America vol. 97,4 (2000): 1423-6. doi:10.1073/pnas.97.4.1423

    https://www.ncbi.nlm.nih.gov/pmc/articles/PMC26449/

    which is a minor modification of the Weiszfeld algorithm.

    You generally just want to call GeometricMedian.solve[points]

    TODO:  If you specify weights, this often doesn't converge.  Try to fix this. 
*/


class GeometricMedian
{
   /** Solves for the geometric median of an array of points and their
       optional weights. 
   
       params:
         points:  An array of n-dimensional points.

         weights:  An optional array of weights for each point.  The number of
                  weights should be the same as the number of points.  If
                  weights is undef, all points are weighted equally.

       returns
          a point containing the GeometricMedian
   */

   class solve[points, weights=undef, relativeError = 1e-10] :=
   {
      y = average[points]
      if length[points] == 1
         return points@0

      if length[points] == 2
      {
         [y, r] = nextPointVardi[points, y, weights]
         return y
      }

      if weights != undef and length[points] != length[weights]
      {
         println["GeometricMedian.solve:  points and weights do not have the same length."]
         return undef
      }
         
      // y = points@0
      r = undef
      lastr = undef
      lasty = undef
      do
      {
         lastr = r
         lasty = y
         [y, r] = nextPointVardi[points, y, weights]
//         println["r is $r"]
      } while lastr == undef or r <= lastr

//      println["Bailing at r=$r\nlasty is $lasty\n    y is $y"]
      if (r > lastr)
         return lasty
      
      return y
   }
   
   /** Perform one round of several calculations.  This takes an array of
       n-dimensional points and an array of weights with the same length.
       (If the weights are undef, they will all be weighted equally.)
       This is actually Weizsfeld's algorithm if you iterate it (and if you
       never have a case where the workingPoint is the same one of the points.

       params:
         points:  An array of n-dimensional points.

         workingPoint:  The current guess for the geometric median.  This is y
                  in the paper.

         weights:  An optional array of weights for each point.  The number of
                  weights should be the same as the number of points.  If
                  weights is undef, all points are weighted equally.

       returns: [T, r, etaY]
         T:  The value of T from eq. 2.4 in the paper above.  This is a point
             approximating the next point.
         r:  A dimensionless term approximately representing the remaining error
             This should monotonically decrease to zero if iterated using the
             Vardi algorithm.
         etaY:  The "multiplicity" at y 
   */

   class calcOnce[points, workingPoint, weights=undef] :=
   {
      // We do this to maintain units of measure in the result.
      vecsum = undef
      sum = undef
      rsum = undef
      etaY = 0
      
      for i = 0 to length[points]-1
      {
         point = points@i
         if weights == undef
            weight = 1
         else
            weight = weights@i
         
         distance = distance[point, workingPoint]  // ||y-x_i||
         // println["Distance is $distance"]
         
         // TODO:  Check if distance is zero magnitude.  This may mean that
         // the points are colinear.
         if point != workingPoint
         {
            //println["Vecsum is $vecsum, point is $point, distance is $distance"]
            vecterm = multiplyScalar[point, weight/distance]
            vecsum = addVectors[vecsum, vecterm]

            // println["sum is $sum, weight is $weight, distance is $distance"]

            sumterm = weight / distance
            if sum == undef
               sum = sumterm
            else
               sum = sum + sumterm

            // Eq. 2.7
            rsumterm = multiplyScalar[subtractVectors[point, workingPoint],
                                      weight/distance]
            rsum = addVectors[rsum, rsumterm]
         } else
            etaY = weight
      }
      
      T = divideScalar[vecsum, sum]
      r = vectorLength[rsum]
      return [T, r, etaY]
   }

   /** This returns the next point in the iteration using the Vardi-Zhang
       algorithm.  That is, eq. 2.6 in the paper.

       It returns a point that is the next point in the iteration and its error.
   
       returns [T, r]
          T: The next point
          r: The relative error.  This should monotonically decrease toward 0.
   */

   class nextPointVardi[points, workingPoint, weights=undef] :=
   {
      [T, r, etaY] = calcOnce[points, workingPoint, weights]
      ratio = ( r==0 || etaY == 0 ? 0 : etaY/r )
      // println["r is $r, etaY is $etaY, ratio is $ratio"]
      T = addVectors[multiplyScalar[T,(1-ratio)],
                     multiplyScalar[workingPoint, min[1,ratio]]]

      return [T, r]
   }

   /** This returns the next point in the iteration using the Weiszfeld
       algorithm.  This will get stuck if the working point is exactly on
       one of the points.

       returns [T, r]
          T: The next point
          r: The relative error.
   */

   class nextPointWeiszfeld[points, workingPoint, weights=undef] :=
   {
      [T, r, etaY] = calcOnce[points, workingPoint, weights]
      return [T, r]
   }

   /** Return the distance between two n-dimensional points, represented as
       arrays. */

   class distance[x, y] :=
   {
      d2 = 0 (x@0)^2   // Preserve dimensions
      for i = 0 to length[x]-1
         d2 = d2 + (x@i - y@i)^2
      
      return sqrt[d2]
   }

   /** Return the distances between an array of points and the working point
       as an array of distances.
   */

   class distances[points, workingPoint] :=
   {
      result = new array
      for p = points
         result.push[distance[p, workingPoint]]

      return result
   }

   /** Return the length of a single n-dimensional vector. */
   class vectorLength[x] :=
   {
      d2 = 0 (x@0)^2   // Preserve dimensions
      for i = 0 to length[x]-1
         d2 = d2 + (x@i)^2
      
      return sqrt[d2]
   }
   
   /** Multiply a vector (point) by a scalar.  It will be nice when Frink does
       vector multiplication automatically.
   */

   class multiplyScalar[vec, scalar] :=
   {
      result = new array[length[vec]]
      for i = rangeOf[vec]
         result@i = vec@i * scalar

       return result
   }

   /** Divide a vector (point) by a scalar.  It will be nice when Frink does
       vector multiplication automatically.
   */

   class divideScalar[vec, scalar] :=
   {
      result = new array[length[vec]]
      for i = rangeOf[vec]
         result@i = vec@i / scalar

      return result
   }
   
   /** Add vectors. */
   class addVectors[v1, v2] :=
   {
      if v1 == undef
         return v2
      
      result = new array[length[v1]]
      for i = rangeOf[v1]
         result@i = v1@i + v2@i

      return result
   }
   
   /** Subtract vectors. */
   class subtractVectors[v1, v2] :=
   {
      result = new array[length[v1]]
      for i = rangeOf[v1]
         result@i = v1@i - v2@i

      return result
   }

   /** Find the average of a bunch of points.  This gives a good first guess
       to the solver.  Input is an array of points. */

   class average[points] :=
   {
      items = length[points]
      result = undef
      for p = points
         result = addVectors[result, p]
      
      return divideScalar[result, items]
   }
}


View or download GeometricMedian.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 18492 days, 21 hours, 42 minutes ago.