我想計算一個二次形式：x' Q y在 Julia 中。
對于這種情況，最有效的計算方法是什么：
- 沒有假設。
Q是對稱的。x和y相同 (x = y)。- 兩者
Q都是對稱的并且x = y。
我知道朱莉婭有dot()。但我想知道它是否比 BLAS 呼叫更快。
uj5u.com熱心網友回復：
現有的答案都很好。但是，還有幾點：
雖然默認情況下 Julia 確實會在這里簡單地呼叫 BLAS，但正如 Oscar 所指出的那樣，一些 BLAS 比其他的更快。特別是，在現代 x86 硬體上，MKL 通常會比 OpenBLAS 快一些。幸運的是，在 Julia 中，手動選擇 BLAS 后端非常容易。只需
using MKL在 Julia 1.7 或更高版本上的 REPL 中鍵入即可通過MKL.jl包將您切換到 MKL 后端。雖然默認情況下不使用它，但現在確實存在一些純 Julia 線性代數包，它們可以匹配甚至擊敗傳統的基于 Fortran 的 BLAS。特別是Octavian.jl：
uj5u.com熱心網友回復：
如果您的矩陣是對稱的，請使用Symmetric包裝器來提高性能（然后呼叫不同的方法）：
julia> a = rand(10000); b = rand(10000);
julia> x = rand(10000, 10000); x = (x x') / 2;
julia> y = Symmetric(x);
julia> @btime dot($a, $x, $b);
47.000 ms (0 allocations: 0 bytes)
julia> @btime dot($a, $y, $b);
27.392 ms (0 allocations: 0 bytes)
如果x與有關選項的討論y參見https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606相同（但總的來說它似乎dot仍然那么快）。
uj5u.com熱心網友回復：
Julia 的 LinearAlgebra 標準庫具有 3-argument 的本機實作dot，以及專門用于對稱/厄米矩陣的版本。您可以在此處和此處查看源代碼。
您可以確認它們沒有使用BenchmarkTools.@btimeor進行分配BenchmarkTools.@ballocated（請記住使用 插入變數$）。矩陣的對稱性被利用了，但是查看源代碼，我看不出如何x == y能夠實作任何嚴重的加速，除了可能節省一些陣列查找。
編輯：要比較 BLAS 版本和本機版本的執行速度，您可以這樣做
1.7.0> using BenchmarkTools, LinearAlgebra
1.7.0> X = rand(100,100); y=rand(100);
1.7.0> @btime $y' * $X * $y
42.800 μs (1 allocation: 896 bytes)
1213.5489200642382
1.7.0> @btime dot($y, $X, $y)
1.540 μs (0 allocations: 0 bytes)
1213.548920064238
這是原生版本的一大勝利。但是，對于更大的矩陣，情況會發生變化：
1.7.0> X = rand(10000,10000); y=rand(10000);
1.7.0> @btime $y' * $X * $y
33.790 ms (2 allocations: 78.17 KiB)
1.2507105095988091e7
1.7.0> @btime dot($y, $X, $y)
44.849 ms (0 allocations: 0 bytes)
1.2507105095988117e7
可能是因為 BLAS 使用執行緒，而dot不是多執行緒。還有一些浮點差異。
uj5u.com熱心網友回復：
您可以使用Tullio.jl撰寫一個優化的回圈，一次完成所有這些。但我認為它不會顯著擊敗 BLAS：
鏈式乘法也很慢，因為它不知道有更好的演算法。
julia> # a, b, x are the same as in Bogumi?'s answer
julia> @btime dot($a, $x, $b);
82.305 ms (0 allocations: 0 bytes)
julia> f(a, x, b) = @tullio r := a[i] * x[i,j] * b[j]
f (generic function with 1 method)
julia> @btime f($a, $x, $b);
80.430 ms (1 allocation: 16 bytes)
添加 LoopVectorization.jl 可能是值得的：
julia> using LoopVectorization
julia> f3(a, x, b) = @tullio r := a[i] * x[i,j] * b[j]
f3 (generic function with 1 method)
julia> @btime f3($a, $x, $b);
73.239 ms (1 allocation: 16 bytes)
但我不知道如何處理對稱情況。
julia> @btime dot($a, $(Symmetric(x)), $b);
42.896 ms (0 allocations: 0 bytes)
雖然可能有線性代數技巧可以用 Tullio.jl 智能地減少它。
在這類問題中，基準權衡就是一切。
轉載請註明出處，本文鏈接：https://www.uj5u.com/yidong/407604.html
標籤：
上一篇：在svg圓圈內顯示影像-（關閉）