How to speed up computation of medoids?

Question

I need to compute medoids (based on SquaredEuclideanDistance for large lists of k-vectors (k = 14). My implementations so far are impractically slow.

Here is a simple implementation, where the aim is the clarity of the code rather than its performance¹. The function medoidSimple takes one argument, cluster (a matrix of reals), and returns a pair consisting of cluster's medoid, and its index in cluster:

medoidSimple[cluster_List/;MatrixQ[cluster, (NumberQ[#] && Im[#] == 0)&]] :=
Module[
  {
     totalDistance
   , distances
   , indexOfMin
  }
  , totalDistance =
        Function[element,
                 Total[SquaredEuclideanDistance[element, #]& /@ cluster]]
  ; distances = totalDistance /@ cluster  (* BOTTLENECK *)
  ; indexOfMin = First[Ordering[distances, 1]]
  ; {cluster[[indexOfMin]], indexOfMin}
];

The bottleneck is the $O(N^2)$ computation of all pairwise distances.

I tried to optimize the code above by aborting the computation of the pairwise sums for any element once the cumulative total for that element exceeded the best minimum value found so far. To my surprise, this "optimization" made my code much slower. I give the code for this version (medoid) below, but here are the timing results:

SeedRandom[0];
cluster = RandomReal[{0, 1}, {5000, 14}];

Timing[medoidSimple[cluster];]
(* {34.428, Null} *)

Timing[medoid[cluster];]
(* {69.628, Null} *)

Next I tried compilation, but all my attempts failed for one reason or another. For example:

compiledTotalDistance = Compile[
    {{element, _Real, 1}, {cluster, _Real, 2}}
  , Total[SquaredEuclideanDistance[element, #]& /@ cluster]
];

compiledTotalDistance[First[cluster], cluster];

Mathematica graphics

I've run out of ideas. Any suggestions would be appreciated.

FWIW, here's the code for my failed optimization.

The main function (medoid) uses the helper function totalDistances. The latter takes as input a point (base), a list of points (peers), and an upper bound (limit). The function returns Infinity if the sum of all the pairwise SquaredEuclideanDistance's between base and the points in peers exceeds limit; otherwise it returns this sum. (The computation of the sum is aborted as soon as soon as its value exceeds limit.)

totalDistance[base_, peers_, limit_] := Module[{accumulator}
  , accumulator[partialSum_, element_] := Module[{nextSum}
      , nextSum = partialSum + SquaredEuclideanDistance[base, element]
      ; If[nextSum > limit, Throw[Infinity], nextSum]
    ]
  ; Catch[Fold[accumulator, 0, peers]]
];

medoid[cluster_] := Module[
  {
     accumulator
   , best
   , medoid
   , indexOfMedoid
  }
  , accumulator[bestSoFar_, index_] := Module[
      {
         limit = First[bestSoFar]
       , element = cluster[[index]]
       , candidate
      }
      , candidate = totalDistance[element, cluster, limit]
      ; If[candidate < limit, {candidate, index}, bestSoFar]
    ]
  ; best = Fold[accumulator, {Infinity, Null}, Range[Length[cluster]]]
  ; indexOfMedoid = best[[2]]
  ; medoid = cluster[[indexOfMedoid]]
  ; {medoid, indexOfMedoid}
];

UPDATE

Summary of (absolute) timings:

medoidSimple->27.8218
medoidMeyer1->3.66165
medoidMeyer2->1.99197
medoidCiao->0.67927
medoidBlochwave1->0.697672
medoidBlochwave2->3.94929
medoidBlochwave3->0.609076

The definitions corresponding to the names above follow.

medoidSimple[
   cluster_List /; MatrixQ[cluster, (NumberQ[#] && Im[#] == 0) &]] := 
  Module[{totalDistance, distances, indexOfMedoid}, 
   totalDistance = 
    Function[element, 
     Total[SquaredEuclideanDistance[element, #] & /@ cluster]]; 
   distances = totalDistance /@ cluster; 
   indexOfMedoid = 
    First[Ordering[distances, 1]]; {cluster[[indexOfMedoid]], 
    indexOfMedoid}];

medoidMeyer1[cluster_] := 
  Block[{distances, indexOfMin}, 
   distances = Total@distMatCompiled[cluster];
   indexOfMin = First[Ordering[distances, 1]];
   {cluster[[indexOfMin]], indexOfMin}];

medoidMeyer2 = Compile[{{cluster, _Real, 2}},
  Block[{len = Length[cluster], bestInd = 1, bestVal, partsum, j},
   bestVal = 
    Total[Table[
      Function[x, x.x][cluster[[1]] - cluster[[i]]], {i, len}]];
   Do[
    partsum = 0.;
    j = 1;
    While[j <= len && partsum < bestVal,
     partsum += Function[x, x.x][cluster[[i]] - cluster[[j]]];
     j++
     ];
    If[j == len + 1 && partsum < bestVal, bestVal = partsum; 
     bestInd = i]
    , {i, 2, len}
    ];
   bestInd
   ]
  , CompilationTarget -> "C"
  ];

medoidCiao = 
  With[{m = #, 
     d = Tr /@ 
       DistanceMatrix[#, 
        DistanceFunction -> SquaredEuclideanDistance]}, {First@m[[#]],
        First@#} &@Pick[Range@Length@d, d, Min@d]] &;

medoidBlochwave1[cluster_] := 
  Module[{distances, indexOfMin}, 
   distances = 
    Tr /@ DistanceMatrix[cluster, 
      DistanceFunction -> SquaredEuclideanDistance];
   indexOfMin = First[Ordering[distances, 1]];
   {cluster[[indexOfMin]], indexOfMin}];

medoidBlochwave2[cluster_] := 
  Module[{distances, indexOfMin}, 
   distances = 
    Total@With[{tr = Transpose[cluster]}, 
      Function[point, Total@Abs[(point - tr)^2]] /@ cluster];
   indexOfMin = First[Ordering[distances, 1]];
   {cluster[[indexOfMin]], indexOfMin}];


myDistMatrix = Compile[{{point, _Real, 1}, {tr, _Real, 2}}, 
                 Total[(point - tr)^2], 
                 RuntimeOptions -> "Speed", 
                 RuntimeAttributes -> {Listable}, 
                 Parallelization -> True];
medoidBlochwave3[cluster_] := 
  Module[{distances, indexOfMin}, 
   distances = Tr /@ myDistMatrix[cluster, Transpose[cluster]];
   indexOfMin = First[Ordering[distances, 1]];
   {cluster[[indexOfMin]], indexOfMin}];

TableForm[(# -> 
     First[AbsoluteTiming[ToExpression[#][cluster];]]) & /@ {"medoidSimple", 
   "medoidMeyer1", "medoidMeyer2", "medoidCiao", "medoidBlochwave1", 
   "medoidBlochwave2", "medoidBlochwave3"}]

^{¹ Well, not exactly. The implementation of medoidSimple does do one thing for the sake of performance: it uses the sum of all pairwise distance, rather than their mean value, as the statistic to optimize. In this case, both strategies are equivalent, since all sets of pairwise distances have the same number of elements.}

I've improved it further! Do you have a C compiler on your system? Can you try my latest version with CompilationTarget->"C"? — dr.blochwave, May 24 '16 at 20:39
@kjo Make sure to use AbsoluteTiming[] and not Timing[] - that's why the results seem to have changed! — dr.blochwave, May 24 '16 at 20:53
It's probably the number of cores available. It would be interesting to see how it scales as the number of points increases from 5000. — dr.blochwave, May 24 '16 at 21:03

ciao · Accepted Answer · 2016-05-24T19:57:05.030

5

If you're on 10.3+, this should be faster (it's two orders of magnitude faster than your first example on my loungebook):

medioid=With[{m = #, d = Tr /@ DistanceMatrix[#, DistanceFunction ->SquaredEuclideanDistance]},
             {First@m[[#]], First@#} &@Pick[Range@Length@d, d, Min@d]]&;

edited May 24 '16 at 19:57

answered May 24 '16 at 19:44

ciao

25,774
2
58
139

Thanks. That's insanely fast. I'll post timings when I have a chance. – kjo May 24 '16 at 20:00

score 5 · Answer 2 · edited Apr 13 '17 at 12:56

Update: You can go even faster with Compile, and exploit the Listable and Parallelization attributes to great effect, if you have a multi-core machine:

SeedRandom[0];
cluster = RandomReal[{0, 1}, {5000, 14}];

myDistMatrix = Compile[{{point, _Real, 1}, {tr, _Real, 2}}, 
                 Total[(point - tr)^2], 
                 RuntimeOptions -> "Speed", 
                 RuntimeAttributes -> {Listable}, 
                 Parallelization -> True,
                 CompilationTarget -> "C"];
medioidCompile[cluster_] := 
  Module[{distances, indexOfMin}, 
   distances = Tr /@ myDistMatrix[cluster, Transpose[cluster]];
   indexOfMin = First[Ordering[distances]];
   {cluster[[indexOfMin]], indexOfMin}];

AbsoluteTiming[medioidCompile[cluster];]
(* 0.35 seconds *)
(* Original: 29.05 seconds *)

If you have Mathematica 10.3, you can use the experimental function DistanceMatrix for a speed-up of 50x:

medioidFast[cluster_] := 
  Module[{distances, indexOfMin}, 
   distances = Tr/@DistanceMatrix[cluster, DistanceFunction -> SquaredEuclideanDistance];
   indexOfMin = First[Ordering[distances]];
   {cluster[[indexOfMin]], indexOfMin}];

AbsoluteTiming[medioidFast[cluster];]
(* 0.521 seconds *)
(* Original: 29.05 seconds *)

(* Also, for fun: *)
medioidFast[cluster]; // RepeatedTiming (* 0.522 seconds *)
medioidCiao[cluster]; // RepeatedTiming (* 0.522 seconds *)

For versions prior to 10.3, or if you don't want to use an experimental function, then the following solution isn't bad. Taking inspiration from Leonid's answer here, you can achieve a speed-up of 12x:

medioid[cluster_] := Module[{distances, indexOfMin}, 
         distances = Tr/@With[{tr = Transpose[cluster]}, 
               Function[point, Total[(point - tr)^2]] /@ cluster];
         indexOfMin = First[Ordering[distances, 1]]; 
         {cluster[[indexOfMin]], indexOfMin}];

AbsoluteTiming[medioid[cluster];]
(* 2.44 seconds *)
(* Original: 29.05 seconds *)

score 4 · Answer 3 · edited Apr 13 '17 at 12:55

This is much faster:

distMatCompiled = Compile[{{cluster, _Real, 2}},
  Outer[Function[diff, diff.diff][#1 - #2] &, cluster, cluster, 1, 1]
  , CompilationTarget -> "C"
  ]

medoidCompiled[cluster_] := Block[{distances, indexOfMin},
  distances = Total@distMatCompiled[cluster];
  indexOfMin = First[Ordering[distances, 1]];
  {cluster[[indexOfMin]], indexOfMin}
  ]

On my machine:

SeedRandom[0];
cluster = RandomReal[{0, 1}, {5000, 14}];

AbsoluteTiming[m1 = medoidSimple[cluster]]
(* 48.116032 *)

AbsoluteTiming[m2 = medoidCompiled[cluster]]
(* 6.779268 *)

m1 == m2
(* True *)

The reason you couldn't compile efficiently in your example is that SquaredEuclideanDistance is not in the list of compilable functions.

UPDATE:

Just for fun I wanted to implement your optimization idea:

medoidCompiled2 = Compile[{{cluster, _Real, 2}},
  Block[{len = Length[cluster], bestInd = 1, bestVal, partsum, j},
   bestVal = 
    Total[Table[
      Function[x, x.x][cluster[[1]] - cluster[[i]]], {i, len}]];
   Do[
    partsum = 0.;
    j = 1;
    While[j <= len && partsum < bestVal,
     partsum += Function[x, x.x][cluster[[i]] - cluster[[j]]];
     j++
     ];
    If[j == len + 1 && partsum < bestVal, bestVal = partsum; 
     bestInd = i]
    , {i, 2, len}
    ];
   bestInd
   ]
  , CompilationTarget -> "C"
  ]

It performs a bit better than my first attempt:

AbsoluteTiming[m3 = medoidCompiled2[cluster]]
(* 4.286030 *)

Last[m1] == m3
(* True *)

Of course this cannot beat the insane timings ciao and blochwave have come up with ;)

Mind timing (if you have 10.x) my answer in your environment? — ciao, May 24 '16 at 19:52
Sorry @ciao, I have 10.0.2.0 and DistanceMatrix doesn't seem to be in that one... @blochwave's answer is comparable to mine, beats mine by a little bit. — Marius Ladegård Meyer, May 24 '16 at 19:55
@ciao - added in 10.3 indeed - my update shows timings for an almost equivalent to your answer using Tr. — dr.blochwave, May 24 '16 at 20:16

How to speed up computation of medoids?

UPDATE

3 Answers3