14
$\begingroup$

I am playing with julialang and loving it, but at the same time noticed that they have a benchmark comparison with Mathematica. I have already submited a better version of recursion_fibonacci and now am looking into recursion_quicksort. I would like to find an implementation in Mathematica that is comparable to C-lang.

Currently, the benchmark uses the following code

(* numeric vector sort *) ClearAll[qsort]; (* qsort[ain_, loin_, hiin_] := Module[ {a = ain, i = loin, j = hiin, lo = loin, hi = hiin, pivot}, While[ i < hi, pivot = a[[BitShiftRight[lo + hi] ]]; While[ i <= j, While[a[[i]] < pivot, i++]; While[a[[j]] > pivot, j--]; If[ i <= j, a[[{i,j}]] = a[[{j, i}]]; i++; j--; ]; ]; If[ lo < j, a = qsort[a, lo, j] ]; {lo, j} = {i, hi}; ]; a ]; *) qsort = Compile[ {{ain, _Real, 1}, {loin, _Integer}, {hiin, _Integer}}, Module[ {a = ain, i = loin, j = hiin, lo = loin, hi = hiin, pivot}, While[ i < hi, pivot = a[[ Floor[(lo + hi)/2] ]]; While[ i <= j, While[a[[i]] < pivot, i++]; While[a[[j]] > pivot, j--]; If[ i <= j, a[[{i,j}]] = a[[{j, i}]]; i++; j--; ]; ]; If[ lo < j, a[[lo;;j]] = qsort[ a[[lo;;j]], 1, j - lo + 1] ]; {lo, j} = {i, hi}; ]; a ] ]; ClearAll[sortperf]; sortperf[n_] := Module[{vec = RandomReal[1, n]}, qsort[vec, 1, n]]; test[OrderedQ[sortperf[5000]] ]; timeit[sortperf[5000], "recursion_quicksort"];

where there is compiled and uncompiled versions of quicksort algorithm. On my laptop the compiled version takes 10.3ms, while uncompiled version takes 137.8ms. I think there is space for improvement since the inbuilt method Sort[] takes only 0.379ms.

How do we speed-up the quicksort algorithm? Bonus points if we don't use Compile[]

Helper functions to run the code above

Needs["CCompilerDriver`"]; If[ Length[CCompilers[]] > 0, $CompilationTarget = "C" ]; ClearAll[timeit]; SetAttributes[timeit, HoldFirst]; timeit[ex_, name_String] := Module[ {t}, t = Infinity; Do[ t = Min[t, N[First[AbsoluteTiming[ex]]]]; , {i, 1, 5} ]; If[$printOutput, Print["mathematica,", name, ",", t*1000]; ]; ]; ClearAll[test]; SetAttributes[test, HoldFirst]; test[ex_] := Assert[ex]; On[Assert];
Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ I would hardly expect any top-level interpreted code to compete with a highly optimized low-level implementation (i.e. the built in), especially if you don't want to compile it. I understand your interest in doing that for recreational / educational purposes though. What strategies have you tried / would you like to implement to speed that up? $\endgroup$ Commented Jun 26, 2018 at 16:13
  • 3
    $\begingroup$ For the recursive_fibonacci, I assume you use memoization and this is not what the Julia-team wanted to measure with this test. This test is supposed to measure the timing of recursive calls which are horribly slow in Mathematica. I've been there myself when I saw the testing code and asked, who on earth would implement fib like that. As it turns out, for their purpose it is the right way to do it. $\endgroup$ Commented Jun 26, 2018 at 16:18
  • $\begingroup$ @halirutan I put the function into Module[], which should not allow any memorisation. Commit is here. Any thoughts? $\endgroup$ Commented Jun 26, 2018 at 16:26
  • 1
    $\begingroup$ Also the matmul example in the Julia benchmarks is complete bogos. Total[Unitize[A.ConjugateTranspose[A] - 200], 2] == 0 performs the task in a tenth of time. Benchmarking matrix multiplication is ridiculous anyways as any reasonable language would delegate that to BLAS routines. The mandel implementation is not Listable. Using Compile with options CompilationTarget -> "C", RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed" would make it 10(!) times faster, bringing us much closer to the C performance. $\endgroup$ Commented Jun 27, 2018 at 9:03
  • 1
    $\begingroup$ @Karolis "Without this trick Mathematica is painfully slow" -- that is "a trick" that comes from the inherent properties of Mathematica's language. It can be seen as a manifestation of its core design principles. (Rule-based programming.) Other languages do not allow so readily dynamic increments of function definitions. $\endgroup$ Commented Jan 12, 2020 at 15:13

1 Answer 1

Reset to default
18
$\begingroup$

I also got angry about those randomly picked and ill-implemented benchmarks by the Julia team. I appreciate their efforts (jit compilers are useful), but the Fibonacci example was straight away ridiculous.

Here is a compiled quick sort implementation that employs a stack in order to avoid recursive calls. The problem with recursion in CompiledFunctions is that it needs to call the main evaluator - and such calls are notoriously slow. I synthesized this implementation from the benchmark implementation and from this guide on geeksforgeeks.org.

Admittedly, one may argue whether this is still a "pure" Mathematica implementation or not (imho, it is rather C), but the Julia team also employed Compile. Hence, once started, one has to do that correctly.

quicksortiterative = Compile[{{ain, _Real, 1}}, Module[{iter, l, h, a, stack, top, x, i, j, ai, aj}, a = ain; top = 0; stack = Table[0, {Length[a] + 1}]; stack[[++top]] = 1; stack[[++top]] = Length[a]; While[top >= 1, h = Compile`GetElement[stack, top--]; l = Compile`GetElement[stack, top--]; x = Compile`GetElement[a, l + Quotient[h - l, 2]]; i = l; j = h; While[i <= j, ai = Compile`GetElement[a, i]; While[ai < x, ai = Compile`GetElement[a, ++i]]; aj = Compile`GetElement[a, j]; While[aj > x, aj = Compile`GetElement[a, --j]]; If[i <= j, a[[i++]] = aj; a[[j--]] = ai; ]; ]; If[j > l, stack[[++top]] = l; stack[[++top]] = j]; If[i < h, stack[[++top]] = i; stack[[++top]] = h]; ]; a ], CompilationTarget -> "C", RuntimeOptions -> "Speed" ];

To my own surpise, this actually beats Mathematica's built-in Sort function:

n = 1000000; a = RandomReal[{-1, 1}, n]; r1 = quicksortiterative[a]; // RepeatedTiming // First r2 = Sort[a]; // RepeatedTiming // First r1 == r2

0.105

0.127

True

Improve this answer
$\endgroup$
7
  • $\begingroup$ I wonder what sort of generalizations or inefficiencies built-in Sort uses that allows this to be faster... I would have assumed that'd be implemented as efficiently as possible. $\endgroup$ Commented Jun 26, 2018 at 20:16
  • 4
    $\begingroup$ So thought I. But Sort can sort everything. This can only sort real numbers. 20% is not that much of a performance loss if you take Sort's flexibility into account, is it? Moreover, the pre-allocated memory for stack is definately oversized; it should need only $O(\log_2(n))$ of memory. But you would need some routines to handle stack overflow... I am sure that Sort is implemented somewhat cleverer in that regard. $\endgroup$ Commented Jun 26, 2018 at 20:21
  • $\begingroup$ I guess I assumed that Sort would dispatch to the fastest routine if, say, it had a PackedArray of numbers. I suppose it could make sense with respect to a memory/speed tradeoff. $\endgroup$ Commented Jun 26, 2018 at 20:24
  • 3
    $\begingroup$ To my great amusement, on my machine your implementation is exactly 20% faster and uses exactly 20% more memory. (I got that from this: rm1 = MaxMemoryUsed[{rt1, r1} = RepeatedTiming@quicksortiterative[a];]; rm2 = MaxMemoryUsed[{rt2, r2} = RepeatedTiming@Sort[a];]; (rt1 - rt2)/rt1 N@(rm1 - rm2)/rm1 $\endgroup$ Commented Jun 26, 2018 at 20:28
  • 4
    $\begingroup$ Thank you very much for your answer. I think the important message from your post is that the recursive-compiled-function approach used for the benchmark is not "the Mathematica way to do it". Note how this is part of many Q&As here on SE, as is code optimization for performance in general. I feel that developers have to acquire quite a big amount of knowledge about the Mathematica language until they are really proficient in writing efficient code. The benchmark results illustrate this fact, rather than qualifying Mathematica as being slower than other languages. $\endgroup$ Commented Jun 27, 2018 at 13:29

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.