10
$\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
  • 2
    $\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
  • $\begingroup$ Karolis, The memoization that @halirutan mentioned is the recFib[n_] := recFib[n] = ... bit in the linked code. You are saving the results of previous calls to recFib. Wrapping the code in a Module does not change that. Try and calculate recFib[30] with and without the ... recFib[n] = ... bit in the definition, and you will see what a tremendous difference that makes. $\endgroup$ – MarcoB Jun 26 '18 at 19:35
  • 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
12
$\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 wether 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

$\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
  • 2
    $\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
  • 2
    $\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
  • 1
    $\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.