Download or view 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]
}
}
Download or view 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 19760 days, 21 hours, 58 minutes ago.