Julia vs. MATLAB - Distance Matrix - Run Time Test

Question

I started learning Julia not a long time ago and I decided to do a simple comparison between Julia and Matlab on a simple code for computing Euclidean distance matrices from a set of high dimensional points.

The task is simple and can be divided into two cases:

Case 1: Given two datasets in the form of n x d matrices, say X1 and X2, compute the pair wise Euclidean distance between each point in X1 and all the points in X2. If X1 is of size n1 x d, and X2 is of size n2 x d, then the resulting Euclidean distance matrix D will be of size n1 x n2. In the general setting, matrix D is not symmetric, and diagonal elements are not equal to zero.

Case 2: Given one dataset in the form of n x d matrix X, compute the pair wise Euclidean distance between all the n points in X. The resulting Euclidean distance matrix D will be of size n x n, symmetric, with zero elements on the main diagonal.

My implementation of these functions in Matlab and in Julia is given below. Note that none of the implementations rely on loops of any sort, but rather simple linear algebra operations. Also, note that the implementation using both languages is very similar.

My expectations before running any tests for these implementations is that the Julia code will be much faster than the Matlab code, and by a significant margin. To my surprise, this was not the case!

The parameters for my experiments are given below with the code. My machine is a MacBook Pro. (15" Mid 2015) with 2.8 GHz Intel Core i7 (Quad Core), and 16 GB 1600 MHz DDR3.

Matlab version: R2018a

Julia version: 0.6.3
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
LAPACK: libopenblas64_
LIBM: libopenlibm
LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

The results are given in Table (1) below.

Table 1: Average time in seconds (with standard deviation) over 30 trials for computing Euclidean distance matrices between two different datasets (Col. 1), and between all pairwise points in one dataset (Col. 2).

          Two Datasets  ||  One Dataset     

---

Matlab:      2.68 (0.12) sec.      1.88 (0.04) sec.

Julia V1:      5.38 (0.17) sec.      4.74 (0.05) sec.

Julia V2:      5.2 (0.1) sec.     

---

I was not expecting this significant difference between both languages. I expected Julia to be faster than Matlab, or at least, as fast as Matlab. It was really a surprise to see that Matlab is almost 2.5 times faster than Julia in this particular task. I didn't want to draw any early conclusions based on these results for few reasons.

First, while I think that my Matlab implementation is as good as it can be, I'm wondering whether my Julia implementation is the best one for this task. I'm still learning Julia and I hope there is a more efficient Julia code that can yield faster computation time for this task. In particular, where is the main bottleneck for Julia in this task? Or, why does Matlab have an edge in this case?

Second, my current Julia package is based on the generic and standard BLAS and LAPACK packages for MacOS. I'm wondering whether JuliaPro with BLAS and LAPACK based on Intel MKL will be faster than the current version I'm using. This is why I opted to get some feedback from more knowledgeable people on StackOverflow.

The third reason is that I'm wondering whether the compile time for Julia was included in the timings shown in Table 1 (2nd and 3rd rows), and whether there is a better way to assess the execution time for a function.

I will appreciate any feedback on my previous three questions.

Thank you!

Hint: This question has been identified as a possible duplicate of another question on StackOverflow. However, this is not entirely true. This question has three aspects as reflected by the answers below. First, yes, one part of the question is related to the comparison of OpenBLAS vs. MKL. Second, it turns out that the implementation as well can be improved as shown by one of the answers. And last, bench-marking the julia code itself can be improved by using BenchmarkTools.jl.

MATLAB

num_trials = 30;
dim = 1000;
n1 = 10000;
n2 = 10000;

T = zeros(num_trials,1);

XX1 = randn(n1,dim);
XX2 = rand(n2,dim);

%%% DIFEERENT MATRICES
DD2ds = zeros(n1,n2);

for (i = 1:num_trials)
  tic;
  DD2ds = distmat_euc2ds(XX1,XX2);
  T(i) = toc;
end

mt = mean(T);
st = std(T);

fprintf(1,'\nDifferent Matrices:: dim: %d, n1 x n2: %d x %d -> Avg. Time %f (+- %f) \n',dim,n1,n2,mt,st);

%%% SAME Matrix 
T = zeros(num_trials,1);

DD1ds = zeros(n1,n1);

for (i = 1:num_trials)
  tic;
  DD1ds = distmat_euc1ds(XX1);
  T(i) = toc;
end

mt = mean(T);
st = std(T);

fprintf(1,'\nSame Matrix:: dim: %d, n1 x n1 : %d x %d -> Avg. Time %f (+- %f) \n\n',dim,n1,n1,mt,st);

distmat_euc2ds.m

function [DD] = distmat_euc2ds (XX1,XX2)
    n1 = size(XX1,1);
    n2 = size(XX2,1);
    DD = sqrt(ones(n1,1)*sum(XX2.^2.0,2)' + (ones(n2,1)*sum(XX1.^2.0,2)')' - 2.*XX1*XX2');
end

distmat_euc1ds.m

function [DD] = distmat_euc1ds (XX)
    n1 = size(XX,1);
    GG = XX*XX';
    DD = sqrt(ones(n1,1)*diag(GG)' + diag(GG)*ones(1,n1) - 2.*GG);
end

JULIA

include("distmat_euc.jl")

num_trials = 30;

dim = 1000;
n1 = 10000;
n2 = 10000;

T = zeros(num_trials);

XX1 = randn(n1,dim)
XX2 = rand(n2,dim)

DD = zeros(n1,n2)

# Euclidean Distance Matrix: Two Different Matrices V1
# ====================================================
for i = 1:num_trials
    tic()
    DD = distmat_eucv1(XX1,XX2)
    T[i] = toq();
end

mt = mean(T)
st = std(T)

println("Different Matrices V1:: dim:$dim, n1 x n2: $n1 x $n2 -> Avg. Time $mt (+- $st)")


# Euclidean Distance Matrix: Two Different Matrices V2
# ====================================================
for i = 1:num_trials
    tic()
    DD = distmat_eucv2(XX1,XX2)
    T[i] = toq();
end

mt = mean(T)
st = std(T)

println("Different Matrices V2:: dim:$dim, n1 x n2: $n1 x $n2 -> Avg. Time $mt (+- $st)")


# Euclidean Distance Matrix: Same Matrix V1
# =========================================
for i = 1:num_trials
    tic()
    DD = distmat_eucv1(XX1)
    T[i] = toq();
end

mt = mean(T)
st = std(T)

println("Same Matrix V1:: dim:$dim, n1 x n2: $n1 x $n2 -> Avg. Time $mt (+- $st)")

distmat_euc.jl

function distmat_eucv1(XX1::Array{Float64,2},XX2::Array{Float64,2})
    (num1,dim1) = size(XX1)
    (num2,dim2) = size(XX2)

    if (dim1 != dim2)
        error("Matrices' 2nd dimensions must agree!")
    end

    DD = sqrt.((ones(num1)*sum(XX2.^2.0,2)') +
         (ones(num2)*sum(XX1.^2.0,2)')' - 2.0.*XX1*XX2');
end


function distmat_eucv2(XX1::Array{Float64,2},XX2::Array{Float64,2})
    (num1,dim1) = size(XX1)
    (num2,dim2) = size(XX2)

    if (dim1 != dim2)
        error("Matrices' 2nd dimensions must agree!")
    end

    DD = (ones(num1)*sum(Base.FastMath.pow_fast.(XX2,2.0),2)') +
         (ones(num2)*sum(Base.FastMath.pow_fast.(XX1,2.0),2)')' -
         Base.LinAlg.BLAS.gemm('N','T',2.0,XX1,XX2);

    DD = Base.FastMath.sqrt_fast.(DD)
end


function distmat_eucv1(XX::Array{Float64,2})
    n = size(XX,1)
    GG = XX*XX';
    DD = sqrt.(ones(n)*diag(GG)' + diag(GG)*ones(1,n) - 2.0.*GG)
end

Answer 1

First question: If I re-write the julia distance function like so:

function dist2(X1::Matrix, X2::Matrix)
    size(X1, 2) != size(X2, 2) && error("Matrices' 2nd dimensions must agree!")
    return sqrt.(sum(abs2, X1, 2) .+ sum(abs2, X2, 2)' .- 2 .* (X1 * X2'))
end

I shave >40% off the execution time.

For a single dataset you can save a bit more, like this:

function dist2(X::Matrix)
    G = X * X'
    dG = diag(G)
    return sqrt.(dG .+ dG' .- 2 .* G)
end

Third question: You should do your benchmarking with BenchmarkTools.jl, and perform the benchmarking like this (remember $ for variable interpolation):

julia> using BenchmarkTools
julia> @btime dist2($XX1, $XX2);

Additionally, you should not do powers using floats, like this: X.^2.0. It is faster, and equally correct to write X.^2.

For multiplication there is no speed difference between 2.0 .* X and 2 .* X, but you should still prefer using an integer, because it is more generic. As an example, if X has Float32 elements, multiplying with 2.0 will promote the array to Float64s, while multiplying with 2 will preserve the eltype.

And finally, note that in new versions of Matlab, too, you can get broadcasting behaviour by simply adding Mx1 arrays with 1xN arrays. There is no need to first expand them by multiplying with ones(...).