In the past few months I have written posts about implementing the bubble sort algorithm in different languages. In the mean while I have gotten some feedback and suggestions regarding improvements to the implementation I made, see the end of the post for the new source code of the algorithms. These often had a quite profound effect on performance.

One of the best tips was to use the R compiler, i.e. the `compiler`

package which is now part of the standard R distribution. This works by simply calling `cmpfun`

:

1 |
function_compiled = cmpfun(function) |

The following table presents the timings in microseconds of the different algorithms:

1 2 3 4 5 6 7 8 |
min median max bubble_sort 1506985.326 1571561.4795 1667417.009 bubble_sort_compiled 308573.445 322742.5950 350686.185 bubble_sort_recursive 1524428.777 1579060.6400 1692169.993 bubble_sort_recursive_compiled 1536280.034 1589214.5010 1660944.670 bubble_sort_cpp 791.354 821.7165 1170.818 sort_julia 958.000 1009.000 1092.000 sort 105.771 170.8530 454.844 |

The following things I find striking:

- Compiling the for/while loop based R solution benefits massively form the compiler package, increasing the speed almost 5 fold.
- Compiling the recursive R based solution does not yield any improvements.
- The C++ solution is obviously much faster than any R based solution, increasing the speed between 1900 and 400 times.
- The Julia based solution is almost as fast as the C++ solution, which is very impressive for a high-level programming language.
- The native R
`sort`

is almost 8 times faster than the fastest bubble sort in C++ and Julia, but`sort`

probably uses a faster algorithm that scales O(n log n) in stead of the 0(n^2) of the bubble sort algorithm.

R (recursive implementation). No improvements.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
larger = function(pair) { if(pair[1] > pair[2]) return(TRUE) else return(FALSE) } swap_if_larger = function(pair) { if(larger(pair)) { return(rev(pair)) } else { return(pair) } } swap_pass = function(vec) { for(i in seq(1, length(vec)-1)) { vec[i:(i+1)] = swap_if_larger(vec[i:(i+1)]) } return(vec) } bubble_sort_recursive = function(vec) { new_vec = swap_pass(vec) if(isTRUE(all.equal(vec, new_vec))) { return(new_vec) } else { return(bubble_sort(new_vec)) } } |

R (for/while loop implementation). This implementation was previously not present.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
bubble_sort = function(vec) { no_passes = 0 while(1) { no_swaps = 0 for (j in 1 : (length(vec) - 1 - no_passes)) { if (vec[j] > vec[j + 1]) { s = vec[j] vec[j] = vec[j+1] vec[j+1] = s no_swaps = no_swaps + 1 } } no_passes = no_passes + 1 if(no_swaps == 0) break } vec } |

C++ (linked into R using Rcpp/inline). Precomputing `vec.size()`

outside the loop improved performance by a factor of 2.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
require(inline) ## for cxxfunction() src = 'Rcpp::NumericVector vec = Rcpp::NumericVector(vec_in); double tmp = 0; int n = vec.size(); int no_swaps; int passes; passes = 0; while(true) { no_swaps = 0; for (int i = 0; i < n - 1 - passes; ++i) { if(vec[i] > vec[i+1]) { no_swaps++; tmp = vec[i]; vec[i] = vec[i+1]; vec[i+1] = tmp; }; }; if(no_swaps == 0) break; passes++; }; return(vec);' bubble_sort_cpp = cxxfunction(signature(vec_in = "numeric"), body=src, plugin="Rcpp") |

Julia. Subtly changing the definition of the for loop (`1:(length(vec_in) - 1 - passes)`

vs `[1:(length(vec_in) - 1 - passes)]`

) improved performance two fold.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
function swap_vector(vec_in, passes) no_swaps = 0 for vec_index in 1:(length(vec_in) - 1 - passes) if vec_in[vec_index] > vec_in[vec_index + 1] no_swaps += 1 tmp = vec_in[vec_index] vec_in[vec_index] = vec_in[vec_index + 1] vec_in[vec_index + 1] = tmp end end return no_swaps end function sort_julia(vec_in) passes = 0 while(true) no_swaps = swap_vector(vec_in, passes) if no_swaps == 0 break end passes += 1 end return vec_in end |

Julia is not interpreted — it’s just-in-time compiled.

Bubble sort is slow, requiring n^2 operations rather than n * log n, like any library function would use.

Thanks for the feedback, I will edit in some more details into my post. I chose to try and implement the bubble sort algorithm quite arbitrarily a few months back.

R’s sort uses Shellsort or Quicksort, which you can specify with the ‘method=’ keyword. It appears it may default to Shellsort.

Note also that if the R structure being sorted has names, order will be preserved in the case of ties, which takes additional work but is probably what you want to have happen. I doubt that Julia does this, just as Julia doesn’t handle NA’s by default.

Right, but quick sort isn’t stable. So using base,

`shell sort`

would be the only go to.Just for completeness: shellsort is not stable either, but has been made stable in R at additional cost. We should avoid comparing different algos when talking about speed of languages. We also should avoid toy-algos like bubble sort. With a serious N*log(N) algo like mergesort or quicksort, benchmarking gains more credibility, especially because it allows to seriously compare at different N, where beyond pure CPU the timings of memory allocation and cache-efficiency strongly inluence real-world performance. Finally: those serious divide-and-conquer algos are often for good reasons implemented recursively. If one language can compile a recursive algo bether than another language, that is an important result. BTW: the above function bubble_sort_recursive can by definition not benefit from compilation because it recalls bubble_sort, not itself. However recalling bubble_sort_recursive() or using the R convention Recall() doesn’t seem to speed-up recursive call either. However, this is likely not a problem of recursion itself but a consequence of creating multiple copies of the array to be sorted (R passes parameters by value). Compared to all that copying, the recursive overhead is minor. I pointed out already in 1998 that pass-by-value kills recursion performance in S+ (see http://math.yorku.ca/Who/Faculty/Monette/S-news/0190.html) and created package ref as a workaround (now defunct). Therefore Julia wins, R looses. Anyhow, something happens when compiling recursive R functions: