Question

Я начал изучать Юлию не так давно, и я решил сделать простой Сравнение Юлии и Матлаба по простому коду для вычисления евклидова матрицы расстояний из множества точек высокой размерности.

Задача проста и может быть разделена на два случая:

Случай 1: Для двух наборов данных в виде n x d матриц, скажем, X1 и X2, вычислить евклидово расстояние между каждой точкой в X1 и всеми точками в X2. Если X1 имеет размер n1 x d, а X2 имеет размер n2 x d, то результирующая евклидова матрица расстояний D будет иметь размер n1 x n2. В общем случае матрица D не является симметричной, а диагональные элементы не равны нулю.

Случай 2: Учитывая один набор данных в виде nxd матрицы X, вычислите попарно евклидово расстояние между всеми n точками в X. Результирующая евклидова матрица расстояний D будет иметь размер nxn, симметричный , с нулевыми элементами на главной диагонали.

Моя реализация этих функций в Matlab и в Julia приведена ниже. Обратите внимание, что ни одна из реализаций не опирается на циклы любого рода, а скорее на простые операции линейной алгебры. Также обратите внимание, что реализация с использованием обоих языков очень похожа.

Мои ожидания до запуска каких-либо тестов для этих реализаций состоят в том, что код Julia будет намного быстрее, чем код Matlab, и со значительным отрывом. К моему удивлению, это был не тот случай!

Параметры для моих экспериментов приведены ниже с кодом. Моя машина - MacBook Pro. (15 "середина 2015 г.) с 2,8 ГГц Intel Core i7 (Quad Core) и 16 ГБ 1600 МГц 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)

Результаты приведены в таблице (1) ниже.

Таблица 1: Среднее время в секундах (со стандартным отклонением) за 30 испытаний для вычисления евклидовых матриц расстояний между двумя различными наборами данных (столбец 1), и между всеми попарными точками в одном наборе данных (столбец 2).

          Two Datasets  ||  One Dataset

Matlab: 2,68 (0,12) сек. 1,88 (0,04) сек.

Юлия V1: 5,38 (0,17) сек. 4,74 (0,05) сек.

Юлия V2: 5,2 (0,1) сек.

Я не ожидал такой существенной разницы между двумя языками. Я ожидал, что Джулия будет быстрее, чем Матлаб, или, по крайней мере, так же быстро, как Матлаб. Было действительно удивительно видеть, что Matlab почти в 2,5 раза быстрее, чем Джулия в этой конкретной задаче. Я не хотел делать какие-либо ранние выводы на основе этих результатов по нескольким причинам.

Сначала , хотя я думаю, что моя реализация Matlab настолько хороша, насколько это возможно, мне интересно, является ли моя реализация Julia лучшей для этой задачи. Я все еще изучаю Julia, и я надеюсь, что есть более эффективный код Julia, который может дать более быстрое время вычисления для этой задачи. В частности, где для Юлии главное узкое место в этой задаче? Или почему у Matlab есть преимущество в этом случае?

Второй , мой текущий пакет Julia основан на стандартных и стандартных пакетах BLAS и LAPACK для MacOS. Мне интересно, будет ли JuliaPro с BLAS и LAPACK на базе Intel MKL быстрее, чем текущая версия, которую я использую. Вот почему я решил получить отзывы от более знающих людей о StackOverflow.

третья причина в том, что мне интересно, было ли время компиляции для Джулии включены в сроки, показанные в Таблице 1 (2-я и 3-я строки), и есть ли лучший способ оценить время выполнения функции.

Буду признателен за любые отзывы о моих предыдущих трех вопросах.

Спасибо!

Подсказка: этот вопрос был определен как возможная копия другого вопроса в StackOverflow. Однако это не совсем так. Этот вопрос имеет три аспекта, отраженных в ответах ниже. Во-первых, да, одна часть вопроса связана со сравнением OpenBLAS и MKL. Во-вторых, оказывается, что реализация также может быть улучшена, как показано одним из ответов. И наконец, тестирование самого кода julia может быть улучшено с помощью BenchmarkTools.jl.

* MATLAB 1063 *

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

DNF · Answer 1 · 06 июля 2018

Первый вопрос: Если я переписываю функцию расстояния Юлия следующим образом:

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 увеличит массив до Float64 с, а умножение на 2 сохранит eltype.

И, наконец, обратите внимание, что в новых версиях Matlab вы также можете добиться широковещательного поведения, просто добавив Mx1 массивы с 1xN массивами. Нет необходимости сначала расширять их, умножая на ones(...).

Юлия против Матлаба. Простой тестовый пример

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Юлия против Матлаба. Простой тестовый пример

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Похожие темы