Benchmarking mir.ndslice + lubeck against numpy and Julia

[p.shkadzko]

Today I decided to write a couple of benchmarks to compare D mir 
with lubeck against Python numpy, then I also added Julia 
snippets. The results appeared to be quite interesting.


Allocation and SVD of [5000 x 10] matrix:

+--------+--------------------+------------+
|  lang  |        libs        | time (sec) |
+--------+--------------------+------------+
| Python | numpy+scipy        |        0.5 |
| Julia  | LinearAlgebra      |     0.0014 |
| D      | mir.ndslice+lubeck |        0.6 |
+--------+--------------------+------------+


Allocation and SVD of [5000 x 100] matrix:

+--------+--------------------+------------+
|  lang  |        libs        | time (sec) |
+--------+--------------------+------------+
| Python | numpy+scipy        |        4.5 |
| Julia  | LinearAlgebra      |      0.024 |
| D      | mir.ndslice+lubeck |        5.2 |
+--------+--------------------+------------+


Allocation and SVD of [5000 x 1000] matrix:

+--------+--------------------+------------+
|  lang  |        libs        | time (sec) |
+--------+--------------------+------------+
| Python | numpy+scipy        |        1.4 |
| Julia  | LinearAlgebra      |          1 |
| D      | mir.ndslice+lubeck |      12.85 |
+--------+--------------------+------------+


Allocation and SVD of [500 x 10000] matrix:

+--------+--------------------+------------+
|  lang  |        libs        | time (sec) |
+--------+--------------------+------------+
| Python | numpy+scipy        |       2.34 |
| Julia  | LinearAlgebra      |        1.1 |
| D      | mir.ndslice+lubeck |         25 |
+--------+--------------------+------------+


Matrices allocation and dot product A [3000 x 3000] * B [3000 x 
3000]

+--------+--------------------+------------+
|  lang  |        libs        | time (sec) |
+--------+--------------------+------------+
| Python | numpy+scipy        |       0.62 |
| Julia  | LinearAlgebra      |      0.215 |
| D      | mir.ndslice+lubeck |        1.5 |
+--------+--------------------+------------+


D lubeck's svd method it quite slow and gets even slower with 
growth of the second dimension while scipy svd becomes 
surprisingly faster. Dot product unfortunately is also 
disappointing. I can only complement Julia on such amazing 
results.

Below is the code I was using.

Allocation and SVD of [A x B] matrix:
Python
------
import numpy as np
from scipy.linalg import svd
import timeit

def svd_fun():
     data = np.random.randint(0, 1024, 5000000).reshape((5000, 
1000))
     u, s, v = svd(data)

print(timeit.timeit(svd_fun, number=1))

Julia
-----
using LinearAlgebra
using BenchmarkTools

function svdFun()
     a = rand([0, 1024], 5000, 1000)
     res = svd(a)
end

@btime svdFun()


D mir.ndslice+lubeck
--------------------
/+dub.sdl:
dependency "mir" version="~>3.2.0"
dependency "lubeck" version="~>1.1.7"
libs "lapack" "openblas"
+/
import std.datetime;
import std.datetime.stopwatch : benchmark;
import std.stdio;
import std.random : Xorshift, unpredictableSeed, uniform;
import std.array: array;
import std.range: generate, take;


import mir.ndslice;
import mir.math.common : optmath;
import lubeck;


void svdFun()
{
     Xorshift rnd;
     rnd.seed(unpredictableSeed);
     auto matrix = generate(() => uniform(0, 1024, rnd))
         .take(5000_000)
         .array
         .sliced(5000, 1000);
     auto svdResult = matrix.svd;
}

void main()
{
     auto svdTime = benchmark!(svdFun)(1);
     writeln(svdTime);

}

Matrices allocation and dot product A [3000 x 3000] * B [3000 x 
3000]
Python
------
def dot_fun():
     a = np.random.random(9000000).reshape((3000, 3000))
     b = np.random.random(9000000).reshape((3000, 3000))
     c = np.dot(a, b)
print(timeit.timeit(dot_fun, number=10)/10)

Julia
-----
function dotFun()
     a = rand([0, 1.0], 3000, 3000)
     b = rand([0, 1.0], 3000, 3000)
     c = dot(a, b)
end
@btime dotFun()

D mir.ndslice+lubeck
--------------------
static @optmath auto rndMatrix(T)(const T maxN, in int dimA, in 
int dimB) {
     Xorshift rnd;
     rnd.seed(unpredictableSeed);
     const amount = dimA * dimB;
     return generate(() => uniform(0, maxN, rnd))
         .take(amount)
         .array
         .sliced(dimA, dimB);
}

static @optmath T fmuladd(T, Z)(const T a, Z z){
     return a + z.a * z.b;
}

void dotFun() {
     auto matrixA = rndMatrix!double(1.0, 3000, 3000);
     auto matrixB = rndMatrix!double(1.0, 3000, 3000);
     auto zipped = zip!true(matrixA, matrixB);
     auto dot = reduce!fmuladd(0.0, zipped);

}

void main()
{
     auto dotTime = benchmark!(dotFun)(1);
     writeln(dotTime);
}