---
author: Stéphane Laurent
date: '2020-10-22'
highlighter: 'pandoc-solarized'
output:
  html_document:
    highlight: tango
    theme: cosmo
  md_document:
    preserve_yaml: True
    variant: markdown
rbloggers: yes
tags: 'R, haskell, julia, Rcpp, C, special-functions, maths'
title: Haskell is fast
---

*Updated title:* **Haskell is fast, but Julia is faster** (see updates
at the end).

My R package 'HypergeoMat' provides a Rcpp implementation of Koev &
Edelman's algorithm for the evaluation of the hypergeometric function of
a matrix argument.

I also implemented this algorithm in
[Julia](https://gist.github.com/stla/e85e2de1ad9aeeebc01583f1d0b8e1d0#file-hypergeompq9-jl)
and in [Haskell](https://github.com/stla/hypergeomPFQ).

So let us benchmark now.

Here is the hypergeometric function of a matrix argument:

$$
{}_pF_q^{(\alpha)}
\left(\begin{matrix} a_1, \ldots, a_p \\ b_1, \ldots, b_q\end{matrix}; X\right) = 
\sum_{k=0}^\infty\sum_{\kappa \vdash k} 
\frac{{(a_1)}_\kappa^{(\alpha)} \cdots {(a_p)}_\kappa^{(\alpha)}}
{{(b_1)}_\kappa^{(\alpha)} \cdots {(b_q)}_\kappa^{(\alpha)}} 
\frac{C_\kappa^{(\alpha)}(X)}{k!}.
$$

Well, I will not explain this expression. But observe that this is a sum
from $k=0$ to $\infty$. The algorithm evaluates the partial sums of this
series, that is, the sum from $k=0$ to an integer $m$.

My Haskell library generates a shared library (a DLL) which can be
called from R. And one can call Julia from R with the help of the
'XRJulia' package. So we will benchmark the three implementations from
R.

Firstly, let's check that they return the same value:

``` {.r}
library(HypergeoMat)
library(XRJulia)
# source the Julia code
juliaSource("HypergeomPQ09.jl")
# load the Haskell DLL
dll <- "libHypergeom.so"
dyn.load(dll)
.C("HsStart")

a <- c(8, 7, 3)
b <- c(9, 16)
x <- c(0.1, 0.2, 0.3)
alpha <- 2
m <- 5L # `m` is the truncation order

hypergeomPFQ(m, a, b, x, alpha)
# 2.116251
juliaEval("hypergeom(5, [8.0, 7.0, 3.0], [9.0, 16.0], [0.1, 0.2, 0.3], 2.0)")
# 2.116251
.Call("hypergeomR", m, a, b, x, alpha)
# 2.116251
```

Well, the same results. Now, let's run a first series of benchmarks, for
$m=5$.

``` {.r}
library(microbenchmark)
microbenchmark(
  Rcpp =
    hypergeomPFQ(m, a, b, x, alpha),
  Julia =
    juliaEval("hypergeom(5, [8.0, 7.0, 3.0], [9.0, 16.0], [0.1, 0.2, 0.3], 2.0)"),
  Haskell =
    .Call("hypergeomR", m, a, b, x, alpha),
  times = 10
)
```

    Unit: microseconds
        expr      min        lq       mean    median        uq       max neval cld
        Rcpp  356.682   623.807   837.7237   827.402  1084.191  1382.500    10  a 
       Julia 4052.000 47767.565 44725.3895 48845.156 50597.779 51308.089    10   b
     Haskell  610.852  1136.963  1343.7442  1289.435  1504.323  2650.976    10  a 

Should we conclude that Rcpp is the winner, and that Julia is slow?
That's not sure. Observe that the unit of these durations is the
microsecond. Perhaps the call to Julia via `juliaEval` is
time-consuming, as well as the call to the Haskell DLL via `.Call`.

So let us try with $m=40$ now.

``` {.r}
m <- 40L
microbenchmark(
  Rcpp =
    hypergeomPFQ(m, a, b, x, alpha),
  Julia =
    juliaEval("hypergeom(40, [8.0, 7.0, 3.0], [9.0, 16.0], [0.1, 0.2, 0.3], 2.0)"),
  Haskell =
    .Call("hypergeomR", m, a, b, x, alpha),
  times = 10
)
```

    Unit: seconds
        expr       min        lq      mean    median        uq      max neval cld
        Rcpp 25.547556 25.924749 26.130888 26.185776 26.354177 26.47846    10   c
       Julia 18.959032 19.088749 19.191394 19.173662 19.291175 19.62415    10  b 
     Haskell  6.642601  6.653627  6.736082  6.735448  6.760926  6.94283    10 a 

This time, the unit is the second. Haskell is clearly the winner,
followed by Julia.

I'm using Julia 1.2.0, and I have been told that there is a great
improvement of performance in Julia 1.5.0, the latest version. I'll try
with Julia 1.5.0 and then I will update this post to show whether there
is a gain of speed.

One should not conclude from this experiment that Haskell *always* beats
C++. That depends on the algorithm we benchmark. This one intensively
uses recursion, and perhaps Haskell is strong when dealing with
recursion.

Don't forget:

``` {.r}
dyn.unload(dll)
```

Update: Julia 1.5 is amazing
============================

Now I upgraded Julia to the latest version, 1.5.2. The results are
amazing:

    Unit: seconds
        expr       min        lq      mean    median        uq       max neval cld
        Rcpp 23.464676 24.392115 24.860484 24.823062 25.013047 27.437176    10   c
       Julia  2.806364  2.852674  3.101521  2.973963  3.363618  3.897855    10 a
     Haskell  6.912441  7.459939  7.648012  7.674404  7.798719  8.322777    10  b

19 seconds for Julia 1.2.0 and 3 seconds for Julia 1.5.2! It beats
Haskell.

Update: even better
===================

Thanks to some advice I got on
[discourse.julialang.org](https://discourse.julialang.org/), I improved
my [Julia
code](https://gist.github.com/stla/e85e2de1ad9aeeebc01583f1d0b8e1d0#file-hypergeompq10-jl),
and it is faster now:

    Unit: seconds
      expr      min       lq     mean   median       uq      max neval
     Julia 1.499753 1.549549 1.750907 1.658282 1.915167 2.428611    10