matrix - JULIA 上的 Inv() 与 '\'

Question

因此，我正在编写一个程序，该程序要求我每次在重复大约 1000 次的循环中分配和反转大小为 100 x 100 的大矩阵 M。

我最初使用的是 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() 比 '\' 运算符占用的内存更少，尽管 '\' 运算符最终更快。

这是因为我在 test2() 中使用了额外的单位矩阵 --> Diagonal(ones(100)) 吗？这是否意味着每次运行循环时，都会分配新的内存部分来存储单位矩阵？

我的原始矩阵𝐌是一个三角矩阵。反转具有如此大量零的矩阵是否会花费更多的内存分配？对于这样的矩阵，使用什么更好: inv() 或 '\' 函数还是有其他更好的策略？

P.S:julia 中的求逆矩阵与 C 和 Python 等其他语言相比如何？当我在用 C 编写的旧程序中运行相同的算法时，花费的时间要少得多......所以我想知道 inv() 函数是否是这里的罪魁祸首。

---

编辑:

正如所指出的，我在输入 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() 中所做的那样。它还增强了功能。

Answer 1

您提出的问题很棘手，并且取决于上下文。我可以根据您的情况回答您，但如果您发布真正的问题，答案可能会改变。

因此，对于您的问题，代码并不等效，因为首先您在 inv(B) 中使用了一些矩阵 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

现在，第二件事是您的代码使用对角矩阵，并且 Julia 足够聪明，可以为此类矩阵的 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

然后自己查找代码。

我假设在你的实际练习中你没有对角矩阵。特别是如果你有 block 矩阵，你可以看看 https://github.com/JuliaArrays/BlockArrays.jl包，因为它可能有适合您的用例的优化方法。您也可以看看https://github.com/JuliaMatrices/BandedMatrices.jl .

总而言之 - 您可以期望 Julia 尝试针对特定用例高度优化代码，因此为了获得用例的明确答案，问题的详细规范非常重要。如果您愿意分享，可以给出更具体的答案。