---
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