Sunday, August 17, 2014

Quicksort speed, just in time compiling and vectorizing

I was reading the Julia documentation the other day. They do speed comparisons to other languages. Obviously R does not come out very well. The R code for quicksort is here and I noticed it was not vectorized at all. So I wondered if it could be improved. A quick check on wikipedia showed that the algorithm displayed by wikipedia is the one used by the Julia team. By abandoning this structure and using a just in time compiler some extra speed can be achieved. Additional bonus is that the algorithm actually got more simple. It should be noted (wikipedia) 'Quicksort can be implemented with an in-place partitioning algorithm, so the entire sort can be done with only O(log n) additional space'. Such an implementation was used by the Julia team, however profiling seems to show that limited additional space was not achieved with the R implementation.

Algorithm

To quote wikipedia again:

Pick an element, called a pivot, from the array.

Reorder the array so that all elements with values less than the pivot come before the pivot, while all elements with values greater than the pivot come after it (equal values can go either way). After this partitioning, the pivot is in its final position. This is called thepartitionoperation.

Recursivelyapply the above steps to the sub-array of elements with smaller values and separately to the sub-array of elements with greater values.

Variations in implementation

After coding the first version, it was tried to make faster variations. The first three version are intended to have less comparisons than wfqs1(). The last one is intended to build a bit more vectorizing in qsort().

wfqs2 <- function(x) {

if (length(x)<2) return(x)

ipivot <- sample(length(x),1)

pivot <- x[ipivot]

c(wfqs2(x[x<pivot]),pivot,wfqs2(x[-ipivot][x[-ipivot]>=pivot]))

}

wfqs3 <- function(x) {

if (length(x)<2) return(x)

ipivot <- sample(length(x),1)

pivot <- x[ipivot]

split <- x<pivot

split2 <- !split

split2[ipivot] <- FALSE

c(wfqs3(x[split]),pivot,wfqs3(x[split2]))

}

wfqs4 <- function(x) {

if (length(x)<2) return(x)

split <- x<x[1]

c(wfqs4(x[split]),x[1],wfqs4(x[!split][-1]))

}

wfqsx = function(a) {

qsort_kernel = function(lo, hi) {

if(lo>=hi) return()

if(hi-lo==1) {

if(a[lo]>a[hi]) {

t <- a[lo]

a[lo] <<- a[hi]

a[hi] <<- t

}

return()

}

goUp <- a[(lo+1):hi]>a[lo]

up <- which(goUp)

do <- which(!goUp)

pivottarget <- lo+length(do)

up <- up[up<=length(do)]+lo

do <- do[do>length(do)]+lo

t <- a[do]

a[do] <<- a[up]

a[up] <<- t

t <- a[pivottarget]

a[pivottarget] <<- a[lo]

a[lo] <<- t

qsort_kernel(lo,pivottarget-1)

qsort_kernel(pivottarget+1, hi)

}

qsort_kernel(1, length(a))

return(a)

}

Analysis

Speed without JIT

The first figure shows that sorting time is roughly equivalent to the number of items to be sorted. The exception there is base:sort() where around 100 items there is no size effect. Even though not very pronounced on this scale, qsort() is slowest and wfqs4() is fastest after base:sort().

The second figure shows speed relative to base:sort(). This shows that qsort() is takes approximately 6 times as much time as wfqs4().

Profiling

Profiling seems to show that time is spend because it is spend. Most of it is not spend on any function which Rprof() registers. Regarding memory usage, if mem.total is anything to go by, qsort() uses quite more than wfqs4().

##### qsort #############################################

$by.self

self.time self.pct total.time total.pct mem.total

"qsort_kernel" 6.72 86.15 7.80 100.00 1671.4

"<GC>" 0.68 8.72 0.68 8.72 112.2

"+" 0.16 2.05 0.16 2.05 31.5

"<" 0.12 1.54 0.12 1.54 21.5

"-" 0.06 0.77 0.06 0.77 5.4

"floor" 0.04 0.51 0.04 0.51 10.2

"<=" 0.02 0.26 0.02 0.26 5.3

$by.total

total.time total.pct mem.total self.time self.pct

"qsort_kernel" 7.80 100.00 1671.4 6.72 86.15

"<Anonymous>" 7.80 100.00 1671.4 0.00 0.00

"eval" 7.80 100.00 1671.4 0.00 0.00

"qsort" 7.80 100.00 1671.4 0.00 0.00

"<GC>" 0.68 8.72 112.2 0.68 8.72

"+" 0.16 2.05 31.5 0.16 2.05

"<" 0.12 1.54 21.5 0.12 1.54

"-" 0.06 0.77 5.4 0.06 0.77

"floor" 0.04 0.51 10.2 0.04 0.51

"<=" 0.02 0.26 5.3 0.02 0.26

$sample.interval

[1] 0.02

$sampling.time

[1] 7.8

##### wfqs4 ###########################################

$by.self

self.time self.pct total.time total.pct mem.total

"wfqs4" 0.98 75.38 1.30 100.00 258.0

"<" 0.10 7.69 0.12 9.23 26.2

"<GC>" 0.08 6.15 0.08 6.15 11.8

"c" 0.06 4.62 0.08 6.15 12.0

"!" 0.04 3.08 0.04 3.08 9.0

"-" 0.02 1.54 0.02 1.54 4.6

".External" 0.02 1.54 0.02 1.54 0.0

$by.total

total.time total.pct mem.total self.time self.pct

"wfqs4" 1.30 100.00 258.0 0.98 75.38

"<Anonymous>" 1.30 100.00 258.0 0.00 0.00

"eval" 1.30 100.00 258.0 0.00 0.00

"<" 0.12 9.23 26.2 0.10 7.69

"<GC>" 0.08 6.15 11.8 0.08 6.15

"c" 0.08 6.15 12.0 0.06 4.62

"!" 0.04 3.08 9.0 0.04 3.08

"-" 0.02 1.54 4.6 0.02 1.54

".External" 0.02 1.54 0.0 0.02 1.54

"runif" 0.02 1.54 0.0 0.00 0.00

$sample.interval

[1] 0.02

$sampling.time

[1] 1.3

JIT compiling

The JIT is a great equalizer. Everybody profits with the exception of base:sort(). It is kind of obvious base:sort() does not profit from the JIT, the part that does the work should be in machine code. In the R code, it seems the various implementations are much closer to each other. Small refinements are apparently swamped by the JIT. Nevertheless, wfsq4() is still fastest after base:sort(). It takes half the time of qsort().

Discussion

It cannot be stated enough: When programming in R, vectorize. Vectors make compact code. Vectors are efficient; a vectorized calculation which does more actual computations can beat complex loops which avoid computations. The speed effect is reduced when the JIT is used. Starting the JIT then should be the first step when an R algorithm is too slow. The second step is starting the profiler and looking for (vectorized) optimization. When these two fail it is time to roll out the big guns; (inline) C code, faster computer or switching to a faster language altogether.

The algorithms chosen by the Julia team were delliberately chosen to highlight areas where Julia performs better than R.

tl;dr: One side of the argument is that it is disingenuous to pick an algorithm that reflects badly on R, when smallish fixes can improve the performance substantially. The other side of the argument is that in Julia, you don't need to worry as much about how you contruct your algorithm, since it doesn't have the same performance blindspots as R.

Wiekvoet

Wiekvoet is about R, JAGS, STAN, and any data I have interest in. Topics range from sensometrics, statistics, chemometrics and biostatistics. For comments or suggestions please email me at wiekvoet at xs4all dot nl.