Julia 与 MATLAB - 距离矩阵 - 运行 时间测试
Julia vs. MATLAB - Distance Matrix - Run Time Test
不久前我开始学习 Julia，我决定做一个简单的
Julia 和 Matlab 在计算欧几里德的简单代码上的比较
来自一组高维点的距离矩阵。
任务简单，分为两种情况：
案例 1：给定两个 n x d 矩阵形式的数据集，比如 X1 和 X2，计算 X1 中每个点与 X2 中所有点之间的成对欧几里德距离.如果 X1 的大小为 n1 x d，X2 的大小为 n2 x d，则生成的欧几里得距离矩阵 D 的大小将为 n1 x n2。一般情况下，矩阵D不对称，对角线元素不为零
情况 2：给定一个 n x d 矩阵 X 形式的数据集，计算 X 中所有 n 个点之间的成对欧几里得距离。生成的欧几里得距离矩阵 D 将大小为 n x n，对称，主对角线上有零个元素。
我在 Matlab 和 Julia 中实现的这些函数如下所示。请注意，none 的实现依赖于任何类型的循环，而是简单的线性代数运算。另外请注意，使用这两种语言的实现非常相似。
在 运行 对这些实现进行任何测试之前，我的期望是 Julia 代码将比 Matlab 代码快得多，而且差距很大。令我惊讶的是，事实并非如此！
下面的代码给出了我实验的参数。我的机器是 MacBook Pro。 （15" 2015 年中）配备 2.8 GHz 英特尔酷睿 i7（四核）和 16 GB 1600 MHz DDR3。
Matlab版本：R2018a
朱莉娅版本：0.6.3
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
LAPACK：libopenblas64_
LIBM：libopenlibm
LLVM：libLLVM-3.9.1（ORCJIT，haswell）
结果在下面的Table(1)中给出。
Table 1：计算两个不同数据集之间的欧氏距离矩阵的 30 次试验的平均时间（以秒为单位）（带标准差）（第 1 栏），
以及一个数据集中的所有成对点之间（第 2 列）。
Two Datasets || One Dataset
Matlab： 2.68 (0.12) 秒。 1.88 (0.04) 秒
Julia 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 代码可以为这项任务产生更快的计算时间。特别是，Julia 在这个任务中的主要瓶颈在哪里？或者，为什么 Matlab 在这种情况下有优势？
其次，我目前的 Julia 包基于适用于 MacOS 的通用标准 BLAS 和 LAPACK 包。我想知道基于英特尔 MKL 的带有 BLAS 和 LAPACK 的 JuliaPro 是否会比我使用的当前版本更快。这就是为什么我选择从 Whosebug 上更有知识的人那里获得一些反馈。
第三个原因是我想知道 Julia 的编译时间是否是
包括在 Table 1（第 2 和第 3 行）中显示的时间，以及是否有更好的方法来评估函数的执行时间。
对于我之前三个问题的任何反馈，我将不胜感激。
谢谢！
提示：该问题已被确定为可能与 Whosebug 上的另一个问题重复。然而，这并不完全正确。这个问题有以下三个方面的答案。首先，是的，问题的一部分与 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
茱莉亚
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
第一个问题：如果我像这样重写 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(...).
来扩展它们
不久前我开始学习 Julia，我决定做一个简单的 Julia 和 Matlab 在计算欧几里德的简单代码上的比较 来自一组高维点的距离矩阵。
任务简单，分为两种情况：
案例 1：给定两个 n x d 矩阵形式的数据集，比如 X1 和 X2，计算 X1 中每个点与 X2 中所有点之间的成对欧几里德距离.如果 X1 的大小为 n1 x d，X2 的大小为 n2 x d，则生成的欧几里得距离矩阵 D 的大小将为 n1 x n2。一般情况下，矩阵D不对称，对角线元素不为零
情况 2：给定一个 n x d 矩阵 X 形式的数据集，计算 X 中所有 n 个点之间的成对欧几里得距离。生成的欧几里得距离矩阵 D 将大小为 n x n，对称，主对角线上有零个元素。
我在 Matlab 和 Julia 中实现的这些函数如下所示。请注意，none 的实现依赖于任何类型的循环，而是简单的线性代数运算。另外请注意，使用这两种语言的实现非常相似。
在 运行 对这些实现进行任何测试之前，我的期望是 Julia 代码将比 Matlab 代码快得多，而且差距很大。令我惊讶的是，事实并非如此！
下面的代码给出了我实验的参数。我的机器是 MacBook Pro。 （15" 2015 年中）配备 2.8 GHz 英特尔酷睿 i7（四核）和 16 GB 1600 MHz DDR3。
Matlab版本：R2018a
朱莉娅版本：0.6.3
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
LAPACK：libopenblas64_
LIBM：libopenlibm
LLVM：libLLVM-3.9.1（ORCJIT，haswell）
结果在下面的Table(1)中给出。
Table 1：计算两个不同数据集之间的欧氏距离矩阵的 30 次试验的平均时间（以秒为单位）（带标准差）（第 1 栏）， 以及一个数据集中的所有成对点之间（第 2 列）。
Two Datasets || One Dataset
Matlab： 2.68 (0.12) 秒。 1.88 (0.04) 秒
Julia 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 代码可以为这项任务产生更快的计算时间。特别是，Julia 在这个任务中的主要瓶颈在哪里？或者，为什么 Matlab 在这种情况下有优势？
其次，我目前的 Julia 包基于适用于 MacOS 的通用标准 BLAS 和 LAPACK 包。我想知道基于英特尔 MKL 的带有 BLAS 和 LAPACK 的 JuliaPro 是否会比我使用的当前版本更快。这就是为什么我选择从 Whosebug 上更有知识的人那里获得一些反馈。
第三个原因是我想知道 Julia 的编译时间是否是 包括在 Table 1（第 2 和第 3 行）中显示的时间，以及是否有更好的方法来评估函数的执行时间。
对于我之前三个问题的任何反馈，我将不胜感激。
谢谢！
提示：该问题已被确定为可能与 Whosebug 上的另一个问题重复。然而，这并不完全正确。这个问题有以下三个方面的答案。首先，是的，问题的一部分与 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
茱莉亚
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
第一个问题：如果我像这样重写 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(...).