Inv () против '\' на Юлии - PullRequest
       23

Inv () против '\' на Юлии

0 голосов
/ 26 февраля 2019

Итак, я работаю над кодированием программы, которая требует от меня назначать и инвертировать большую матрицу M размером 100 x 100 каждый раз в цикле, который повторяется примерно 1000 раз.

Первоначально я использовал функцию inv (), но, поскольку это занимает много времени, я хочу оптимизировать свою программу, чтобы она работала быстрее.Поэтому я написал некоторый фиктивный код для проверки того, что может замедлить процесс:

function test1()
  for i in (1:100)
    ?=Diagonal(rand(100))  
    inverse_?=inv(B)
 end
end

using BenchmarkTools
@benchmark test1()
---------------------------------
BenchmarkTools.Trial: 
memory estimate:  178.13 KiB
allocs estimate:  400
--------------
minimum time:     67.991 μs (0.00% GC)
median time:      71.032 μs (0.00% GC)
mean time:        89.125 μs (19.43% GC)
maximum time:     2.490 ms (96.64% GC)
--------------
samples:          10000
evals/sample:     1

Когда я использую оператор '\' для оценки обратного:

function test2()
for i in (1:100)
    ?=Diagonal(rand(100))  
   inverse_?=?\Diagonal(ones(100))
 end
end
using BenchmarkTools
@benchmark test2()
-----------------------------------
BenchmarkTools.Trial: 
memory estimate:  267.19 KiB
allocs estimate:  600
--------------
minimum time:     53.728 μs (0.00% GC)
median time:      56.955 μs (0.00% GC)
mean time:        84.430 μs (30.96% GC)
maximum time:     2.474 ms (96.95% GC)
--------------
samples:          10000
evals/sample:     1

Я могуобратите внимание, что inv () занимает меньше памяти, чем оператор '\', хотя оператор '\' в конце быстрее.

Это потому, что я использую дополнительную матрицу идентичности -> Диагональ (единицы)(100)) в test2 ()?Означает ли это, что каждый раз, когда я запускаю свой цикл, для хранения матрицы идентификаторов выделяется новая часть памяти?

Моя исходная матрица ? является трехугольной матрицей.Стоит ли переворачивать матрицу с таким большим количеством нулей больше памяти?Для такой матрицы, что лучше использовать: функцию inv () или '\', или есть какая-то другая лучшая стратегия?

PS: Как инвертирующие матрицы в julia сравниваются с другими языками, такими как C и Python?Когда я запустил тот же алгоритм в моей более старой программе, написанной на C, это заняло значительно меньше времени ... поэтому мне было интересно, была ли здесь функция inv ().


EDIT:

Итак, как было указано, я сделал опечатку при вводе функции test1 ().Это на самом деле

function test1()
for i in (1:100)
   ?=Diagonal(rand(100))  
   inverse_?=inv(?)
  end
end

Однако моя проблема осталась прежней, функция test1 () выделяет меньше памяти, но занимает больше времени:

using BenchmarkTools
@benchmark test1()

>BenchmarkTools.Trial: 
memory estimate:  178.13 KiB
allocs estimate:  400
--------------
minimum time:     68.640 μs (0.00% GC)
median time:      71.240 μs (0.00% GC)
mean time:        90.468 μs (20.23% GC)
maximum time:     3.455 ms (97.41% GC)

samples:          10000
evals/sample:     1

using BenchmarkTools

@benchmark test2()

BenchmarkTools.Trial: 
memory estimate:  267.19 KiB
allocs estimate:  600
--------------
minimum time:     54.368 μs (0.00% GC)
median time:      57.162 μs (0.00% GC)
mean time:        86.380 μs (31.68% GC)
maximum time:     3.021 ms (97.52% GC)
--------------
samples:          10000
evals/sample:     1

Я также проверил некоторые другиеварианты функции test2 ():

function test3()
for i in (1:100)
   ?=Diagonal(rand(100)) 
    ?=Diagonal(ones(100))
   inverse_?=?\?
 end
end


function test4(?)
for i in (1:100)
     ?=Diagonal(rand(100))  
   inverse_?=?\?
end
end

using BenchmarkTools

@benchmark test3()

>BenchmarkTools.Trial: 
 memory estimate:  267.19 KiB
allocs estimate:  600
 --------------
minimum time:     54.248 μs (0.00% GC)
median time:      57.120 μs (0.00% GC)
mean time:        86.628 μs (32.01% GC)
maximum time:     3.151 ms (97.23% GC)
 --------------
samples:          10000
evals/sample:     1

using BenchmarkTools

@benchmark test4(Diagonal(ones(100)))

>BenchmarkTools.Trial: 
 memory estimate:  179.02 KiB
 allocs estimate:  402
 --------------
 minimum time:     48.556 μs (0.00% GC)
 median time:      52.731 μs (0.00% GC)
 mean time:        72.193 μs (25.48% GC)
 maximum time:     3.015 ms (97.32% GC)
 --------------
 samples:          10000
 evals/sample:     1

test2 () и test3 () эквивалентны.Я понял, что могу выполнить дополнительное выделение памяти в test2 (), передав вместо этого матрицу идентификации как переменную, как это было сделано в test4 ().Это также ускоряет функцию.

1 Ответ

0 голосов
/ 26 февраля 2019

Вопрос, который вы задаете, сложен и зависит от контекста.Я могу ответить вам в вашем контексте, но если вы опубликуете реальную проблему, ответ может измениться.

Так что для вашего вопроса коды не эквивалентны, потому что в первом вы используете матрицу B в inv(B), который не определен (и, вероятно, является глобальным, нестабильным, переменным типа). Если вы измените B на ?, на самом деле первый код будет немного быстрее:

julia> function test1()
         for i in (1:100)
           ?=Diagonal(rand(100))
           inverse_?=inv(?)
        end
       end
test1 (generic function with 1 method)

julia> function test2()
       for i in (1:100)
           ?=Diagonal(rand(100))
          inverse_?=?\Diagonal(ones(100))
        end
       end
test2 (generic function with 1 method)

julia> using BenchmarkTools

julia> @benchmark test1()
BenchmarkTools.Trial:
  memory estimate:  178.13 KiB
  allocs estimate:  400
  --------------
  minimum time:     28.273 μs (0.00% GC)
  median time:      32.900 μs (0.00% GC)
  mean time:        43.447 μs (14.28% GC)
  maximum time:     34.779 ms (99.70% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark test2()
BenchmarkTools.Trial:
  memory estimate:  267.19 KiB
  allocs estimate:  600
  --------------
  minimum time:     28.273 μs (0.00% GC)
  median time:      33.928 μs (0.00% GC)
  mean time:        45.907 μs (15.25% GC)
  maximum time:     34.718 ms (99.74% GC)
  --------------
  samples:          10000
  evals/sample:     1

ТеперьВо-вторых, ваш код использует диагональные матрицы, а Джулия достаточно умна, чтобы иметь специальные методы для inv и \ для матриц этого типа.Их определения следующие:

(\)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .\ Db.diag)

function inv(D::Diagonal{T}) where T
    Di = similar(D.diag, typeof(inv(zero(T))))
    for i = 1:length(D.diag)
        if D.diag[i] == zero(T)
            throw(SingularException(i))
        end
        Di[i] = inv(D.diag[i])
    end
    Diagonal(Di)
end

И вы можете видеть, что такой пример не является полностью репрезентативным для общего случая (если бы матрицы не были диагональными, использовались бы другие методы).Вы можете проверить, какие методы используются следующим образом:

julia> @which �\Diagonal(ones(100))
\(Da::Diagonal, Db::Diagonal) in LinearAlgebra at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\diagonal.jl:493

julia> @which inv(�)
inv(D::Diagonal{T,V} where V<:AbstractArray{T,1}) where T in LinearAlgebra at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\diagonal.jl:496

и посмотреть код самостоятельно.

Я предполагаю, что в вашем реальном упражнении у вас нет диагональных матриц.В частности, если у вас есть блочные матрицы, вы можете взглянуть на пакет https://github.com/JuliaArrays/BlockArrays.jl, так как он мог бы оптимизировать методы для вашего случая использования.Вы также можете взглянуть на https://github.com/JuliaMatrices/BandedMatrices.jl.

Вкратце - вы можете ожидать, что Джулия попытается оптимизировать код для конкретного варианта использования, поэтому для получения подробного ответа для вашего варианта использования подробноспецификация вашей проблемы будет иметь значение.Если вы поделитесь им, можно дать более конкретный ответ.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...