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];
$\endgroup$
  • 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$ – MarcoB Jun 26 '18 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$ – halirutan Jun 26 '18 at 16:18
  • $\begingroup$ @halirutan I put the function into Module[], which should not allow any memorisation. Commit is here. Any thoughts? $\endgroup$ – Karolis Jun 26 '18 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$ – Henrik Schumacher Jun 27 '18 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$ – Anton Antonov Jan 12 at 15:13
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$
  • $\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$ – b3m2a1 Jun 26 '18 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$ – Henrik Schumacher Jun 26 '18 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$ – b3m2a1 Jun 26 '18 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$ – b3m2a1 Jun 26 '18 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$ – Theo Tiger Jun 27 '18 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.