最初這是mathematica.SE 中出現的一個問題，但由于討論涉及多種編程語言，我認為最好重新表述一下并將其發布在這里。
簡而言之，michalkvasnicka在下面的 MATLAB 示例中發現
s = 15000;
tic
% for-loop version
H = zeros(s,s);
for c = 1:s
for r = 1:s
H(r,c) = 1/(r c-1);
end
end
toc
%Elapsed time is 1.359625 seconds.... For-loop
tic;
% vectorized version
c = 1:s;
r = c';
HH=1./(r c-1);
toc
%Elapsed time is 0.047916 seconds.... Vectorized
isequal(H,HH)
矢量化代碼段比純 for 回圈代碼段快 25 倍以上。雖然我無法訪問 MATLAB，因此無法自己測驗樣本，但時間1.359625似乎表明它是在普通 PC 上測驗的，就像我的一樣。
但是我無法用其他語言（如 fortran 或 julia）重現時間！（我們知道，他們都以數值計算的性能而聞名。好吧，我承認我絕不是 Fortran 或 julia 的專家。）
以下是我用于測驗的樣品。我正在使用配備 i7-8565U CPU、Win 10 的筆記本電腦。
復式
fortran 代碼使用 gfortran 編譯（TDM-GCC-10.3.0-2，帶有編譯選項 -Ofast）。
program tst
use, intrinsic :: iso_fortran_env
implicit none
integer,parameter::s=15000
integer::r,c
real(real64)::hmn(s,s)
do r=1,s
do c=1, s
hmn(r,c)=1._real64/(r c - 1)
end do
end do
print *, hmn(s,s)
end program
編譯時間：0.2057823秒
執行時間：0.7179657秒
朱莉婭
Julia 的版本是 1.6.3。
@time (s=15000; Hmm=[1. /(r c-1) for r=1:s,c=1:s];)
計時：0.7945998秒
問題來了：
MATLAB 的時序可靠嗎？
如果第一個問題的答案是肯定的，那么我們如何
0.05使用 julia、fortran 或任何其他編程語言重現性能（對于 2 GHz CPU，時間應該在幾秒鐘左右）？
uj5u.com熱心網友回復：
tic/toc應該沒問題，但看起來時間被記憶體預分配所扭曲。
我可以重現與您的 MATLAB 示例類似的時間，但是
第一次運行（
clear作業區）- 回圈方法需要 2.08 秒
- 矢量化方法需要 1.04 秒
- 矢量化可節省 50% 的執行時間
第二次運行（作業區未清除）
- 回圈方法需要 2.55 秒
- 矢量化方法需要 0.065 秒
- 矢量化“節省”了 97.5% 的執行時間
我的猜測是，由于回圈方法通過 顯式創建了一個新矩陣zeros，因此每次運行時都會從頭開始重新分配記憶體，并且您看不到后續運行的速度提高。
但是，當HH保留在記憶體中并且該HH=___行輸出相同大小的矩陣時，我懷疑 MATLAB 正在執行一些巧妙的記憶體分配以加快操作速度。
我們可以通過以下測驗來證明這個理論：
Test Num | Workspace cleared | s | Loop (sec) | Vectorised (sec)
1 | Yes | 15000 | 2.10 | 1.41
2 | No | 15000 | 2.73 | 0.07
3 | No | 15000 | 2.50 | 0.07
4 | No | 15001 | 2.74 | 1.73
請參閱測驗 2 和 3 之間的差異，這就是為什么timeit對平均運行時間有幫助的原因（參見腳注）。測驗 3 和測驗 4 之間的輸出大小差異非常小，但對于矢量化方法，執行時間恢復到與測驗 1 中類似的幅度，這表明重新分配以創建HH大部分時間的成本。
腳注：tic/ tocMATLAB 中的計時可以通過使用內置timeit函式來改進，該函式基本上取多次運行的平均值。從timeit雖然的作業中觀察到的一件有趣的事情是，它通過多次呼叫它來明確地“預熱”（參考評論）tic/toc函式。從清晰的作業區（沒有中間代碼）運行tic/toc幾次時，您可以看到第一次呼叫比后續呼叫花費的時間更長，因為必須有一些開銷來初始化計時器。
uj5u.com熱心網友回復：
只是在 Julia 方面添加 - 確保您使用BenchmarkTools基準測驗，將要進行基準測驗的代碼包裝在函式中，以免在全域范圍內進行基準測驗，并插入您傳遞給@btime.
這是我將如何做到的：
julia> s = 15_000;
julia> function f_loop!(H)
for c ∈ 1:size(H, 1)
for r ∈ 1:size(H, 1)
H[r, c] = 1 / (r c - 1)
end
end
end
f_loop! (generic function with 1 method)
julia> function f_vec!(H)
c = 1:size(H, 1)
r = c'
H .= 1 ./ (r . c .- 1)
end
f_vec! (generic function with 1 method)
julia> H = zeros(s, s);
julia> using BenchmarkTools
julia> @btime f_loop!($H);
625.891 ms (0 allocations: 0 bytes)
julia> H = zeros(s, s);
julia> @btime f_vec!($H);
625.248 ms (0 allocations: 0 bytes)
所以這兩個版本同時出現，這是我對這樣一個簡單的操作所期望的，正確的型別推斷代碼應該編譯成大致相同的機器代碼。
uj5u.com熱心網友回復：
對于 Julia，@time包括僅一次運行中涉及的任何編譯時間，并且往往變化很大。BenchmarkTools 多次運行該運算式，并@btime報告最短時間：
julia> @time (s=15000; [1/(r c-1) for r in 1:s, c in 1:s]); # as above
0.268203 seconds (276.67 k allocations: 1.691 GiB, 5.94% gc time, 11.74% compilation time)
julia> using BenchmarkTools
julia> @btime (s=15000; [1/(r c-1) for r in 1:s, c in 1:s]);
214.573 ms (2 allocations: 1.68 GiB)
julia> @time begin s=15000; 1 ./ ((1:s) . (1:s)' .- 1) end; # broadcasting instead of a comprehension
0.274145 seconds (619.10 k allocations: 1.708 GiB, 9.55% gc time, 20.02% compilation time)
julia> @btime begin s=15000; 1 ./ ((1:s) . (1:s)' .- 1) end;
167.883 ms (2 allocations: 1.68 GiB)
如上所述，大部分時間都是關于分配陣列。單獨計時有點棘手，因為有時可能會很懶惰：
julia> mat = @btime zeros(15_000, 15_000); # explicitly filled with zeros
132.559 ms (2 allocations: 1.68 GiB)
julia> @btime $mat .= 1 ./ ((1:15_000) . (1:15_000)' .- 1);
52.753 ms (0 allocations: 0 bytes)
julia> mat2 = @btime Array{Float64}(undef, 15_000, 15_000); # sometimes lazier
3.583 μs (2 allocations: 1.68 GiB)
julia> @btime mat2 .= 1 ./ ((1:15_000) . (1:15_000)' .- 1) setup=(mat2=Array{Float64}(undef, 15_000, 15_000)) evals=1;
169.519 ms (0 allocations: 0 bytes)
julia> @btime mat .= 1 ./ ((1:15_000) . (1:15_000)' .- 1) setup=(mat=zeros(Float64, 15_000, 15_000)) evals=1;
54.374 ms (0 allocations: 0 bytes)
撰寫一組明確的for回圈應該與廣播的速度大致相同。一種方法是：
julia> using Tullio
julia> @btime @tullio $mat[r,c] = 1/(r c-1) threads=false;
52.872 ms (0 allocations: 0 bytes)
這些操作都不是多執行緒的。這很容易添加，但對如此廉價的操作無濟于事，至少在我的計算機上是這樣。
uj5u.com熱心網友回復：
我希望以下修改后的基準可以為這個問題帶來一些新的亮點：
s = 15000;
tic
% for-loop version
H = zeros(s,s);
for i =1:10
for c = 1:s
for r = 1:s
H(r,c) = H(r,c) 1/(r c-1 i);
end
end
end
toc
tic;
% vectorized version
HH = zeros(s,s);
c = 1:s;
r = c';
for i=1:10
HH= HH 1./(r c-1 i);
end
toc
isequal(H,HH)
在這種情況下，通過在每次 for 回圈（通過“i”）迭代時更改矩陣 H (HH) 來避免任何型別的“兌現”。
在這種情況下，我們得到：
Elapsed time is 3.737275 seconds. (for-loop)
Elapsed time is 1.143387 seconds. (vectorized)
因此，由于矢量化，性能仍有提升（~ 3 倍），這可能是通過矢量化 Matlab 命令的隱式多執行緒實作來完成的。
是的，tic/toc vs timeit 并不是嚴格一致的，但是整體的計時功能非常相似。
轉載請註明出處，本文鏈接：https://www.uj5u.com/qiye/399408.html
標籤：performance matlab math julia vectorization
上一篇：在字串中無限匹配比較運算子