朱莉娅与Matlab.一个简单的测试用例

Question

我很久以前就开始学习Julia了,我决定在Julia和Matlab之间进行一个简单的比较,用一个简单的代码来计算一组高维点的欧几里德距离矩阵.

任务很简单,可分为两种情况:

情况1:给定nxd矩阵形式的两个数据集,比如X1和X2,计算X1中每个点与X2中所有点之间的成对欧几里德距离.如果X1的大小为n1×d,并且X2的大小为n2×d,则得到的欧几里德距离矩阵D的大小为n1×n2.在一般设置中,矩阵D不对称,并且对角元素不等于零.

情况2:给定nxd矩阵X形式的一个数据集,计算X中所有n个点之间的成对欧几里德距离.得到的欧几里德距离矩阵D的大小为nxn,对称,主对角线上的元素为零.

我在Matlab和Julia中实现这些函数的方法如下.请注意,没有一个实现依赖于任何类型的循环,而是简单的线性代数运算.另请注意,使用这两种语言的实现非常相似.

在为这些实现运行任何测试之前,我的期望是Julia代码将比Matlab代码快得多,并且显着提高.令我惊讶的是,事实并非如此!

我的实验参数在下面给出了代码.我的机器是MacBook Pro.(15"2015年中),2.8 GHz Intel Core i7(四核)和16 GB 1600 MHz DDR3.

Matlab版本:R2018a

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

结果列于下表(1)中.

表1:计算两个不同数据集(第1列)之间以及一个数据集(第2列)中所有成对点之间的欧几里德距离矩阵的30次试验的平均时间(以标准偏差计).

          Two Datasets  ||  One Dataset     

---

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

朱莉娅V1:5.38(0.17)秒.4.74(0.05)秒

Julia V2:5.2(0.1)秒.     

---

我没想到两种语言之间存在这种显着差异.我希望Julia比Matlab更快,或者至少和Matlab一样快.令人惊讶的是,在这项特定任务中,Matlab的速度几乎是Julia的2.5倍.我不想基于这些结果得出任何早期结论,原因很少.

首先,虽然我认为我的Matlab实现尽可能好,但我想知道我的Julia实现是否是这项任务的最佳实现.我还在学习Julia,我希望有一个更高效的Julia代码可以为这个任务带来更快的计算时间.特别是,朱莉娅在这项任务中的主要瓶颈在哪里？或者,为什么Matlab在这种情况下有优势？

其次,我目前的Julia软件包基于MacOS的通用和标准BLAS和LAPACK软件包.我想知道基于英特尔MKL的基于BLAS和LAPACK的JuliaPro是否会比我正在使用的当前版本更快.这就是为什么我选择从StackOverflow上更多知识渊博的人那里获得一些反馈的原因.

的第三原因是,我想知道为朱编译时间是否包括在示于表1(第二和第三行)的定时,以及是否有更好的方法来评估的功能的执行时间.

我将很感激对前三个问题的任何反馈.

谢谢!

提示:此问题已被识别为StackOverflow上另一个问题的可能重复.但是,这并非完全正确.这个问题有三个方面,如下面的答案所反映.首先,是的,问题的一部分与OpenBLAS与MKL的比较有关.其次,事实证明,如其中一个答案所示,也可以改进实施.最后,使用BenchmarkTools.jl可以改进julia代码本身的基准标记.

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

第一个问题:如果我像这样重写julia距离函数:

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

我的执行时间缩短了40%.

对于单个数据集,您可以节省更多,如下所示:

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

第三个问题:您应该使用BenchmarkTools.jl进行基准测试,并执行这样的基准测试(记住$变量插值):

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

此外,你不应该使用浮动功能,如下所示:X.^2.0.它写得更快,同样正确X.^2.

对于乘法而言,2.0 .* X和之间没有速度差异2 .* X,但您仍然应该更喜欢使用整数,因为它更通用.例如,如果X有Float32元素,乘以2.0将将数组提升为Float64s,而乘以2将保留eltype.

最后,请注意,在新版本的Matlab中,您只需添加Mx1带1xN数组的数组即可获得广播行为.没有必要先通过乘以扩展它们ones(...).