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][1]. **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 ;) [1]: https://mathematica.stackexchange.com/questions/1096/list-of-compilable-functions