matlab - Julia vs.MATLAB-距离矩阵-运行时测试

Question

不久前我开始学习茱莉亚，我决定做一个简单的
Julia与Matlab计算欧几里德数的比较
一组高维点的距离矩阵。
任务简单，可分为两种情况：
例1：给定两个n x d矩阵形式的数据集，比如X1和X2，计算X1中每个点与X2中所有点之间的成对欧氏距离如果X1的大小为n1×d，X2的大小为n2×d，则得到的欧氏距离矩阵d的大小为n1×n2在一般情况下，矩阵D是非对称的，对角线元素不等于零。
情形2：给定一个n x d矩阵x形式的数据集，计算x中所有n个点之间的成对欧氏距离。得到的欧氏距离矩阵d的大小为nxn，对称，主对角线上有零个元素。
下面给出了我在Matlab和Julia中实现这些函数的方法注意，没有一个实现依赖于任何类型的循环，而是简单的线性代数运算另外，请注意，使用这两种语言的实现非常相似。
在对这些实现运行任何测试之前，我的期望是Julia代码将比Matlab代码快得多，并且有很大的优势令我惊讶的是，事实并非如此！
下面给出了我的实验参数和代码我的机器是MacBook Pro（2015年年中15英寸），2.8 GHz Intel Core i7（四核）和16 GB 1600 MHz DDR3。
Matlab版本：R2018a
Julia版本：0.6.3
BLAS:libopenblas（使用64位动态关系Haswell）
拉帕克：libopenblas64_
LIBM:libopenlibm
LLVM:libLLVM-3.9.1（ORCJIT，哈斯韦尔）
结果见下表（1）。
表1：计算两个不同数据集之间欧几里德距离矩阵的30次试验的平均时间（以秒为单位）（标准差）（第1列），
在一个数据集中的所有成对点之间（第2列）。

          Two Datasets  ||  One Dataset     

Matlab:2.68（0.12）秒1.88（0.04）秒。
朱莉娅1:5.38（0.17）秒4.74（0.05）秒。
茱莉亚2:5.2（0.1）秒    
我没想到两种语言之间会有这么大的差异我希望茱莉亚比Matlab快，或者至少和Matlab一样快在这个特殊的任务中，Matlab的速度几乎是Julia的2.5倍，这真是令人惊讶我不想基于这些结果得出任何早期结论，原因不多。
首先，虽然我认为我的Matlab实现是最好的，但我想知道我的Julia实现是否是完成此任务的最佳实现我仍然在学习Julia，我希望有一个更有效的Julia代码，可以产生更快的计算时间来完成这项任务特别是，朱莉娅在这项任务中的主要瓶颈在哪里或者，在这种情况下，为什么Matlab有优势？
其次，我目前的Julia包是基于MacOS的通用和标准BLAS和LAPACK包的我想知道基于Intel MKL的带有BLAS和LAPACK的JuliaPro是否会比我现在使用的版本更快这就是为什么我选择从StackOverflow上的更有知识的人那里获得一些反馈。
第三个原因是我想知道Julia的编译时间是否
包括在表1所示的计时中（第2行和第3行），以及是否有更好的方法来评估函数的执行时间。
我将非常感谢您对我前三个问题的反馈。
谢谢您！
提示：此问题已被标识为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);

蒸馏水

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

距离

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

朱莉娅

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(...)相乘来扩展它们。