MATLAB的时序可靠吗？如果是，我们可以用 julia、fortran 等重现性能吗？
Is the timing of MATLAB reliable? If yes, can we reproduce the performance with julia, fortran, etc.?
最初这是一个问题in mathematica.SE，但由于讨论涉及多种编程语言，我认为最好改写一下，post 在这里。
简而言之，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
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的时序靠谱吗？
如果第一个问题的答案是肯定的，那么我们如何用 julia 重现性能（对于 2 GHz CPU，时间应该在 0.05 秒左右） 、fortran 或任何其他编程语言？
tic/toc 应该没问题，但看起来时间因内存预分配而出现偏差。
我可以重现与您的 MATLAB 示例类似的时序，但是
首先 运行（clear 工作区）
- 循环方法需要 2.08 秒
- 矢量化方法需要 1.04 秒
- 矢量化节省了 50% 的执行时间
第二个 运行（工作区未清除）
- 循环方法需要 2.55 秒
- 矢量化方法需要 0.065 秒
- 矢量化“节省”了 97.5% 的执行时间
我的猜测是，由于循环方法通过 zeros 显式地创建了一个新矩阵，内存在每个 运行 上都从头开始重新分配，并且您看不到随后的速度改进运行s.
然而，当 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 之间的输出大小差异非常小，但执行时间 returns 与向量化方法的测试 1 中的执行时间相似，表明重新分配以创建 HH花费大部分时间。
脚注：tic/toc MATLAB 中的时序可以通过使用内置的 timeit 函数来改进，该函数本质上是取几个 运行 的平均值.从 timeit 的工作中观察到的一件有趣的事情是，它通过多次调用 tic/toc 函数来明确地“预热”（引用评论）函数。您可以从清晰的工作区（没有中间代码）中多次 运行ning tic/toc 时看到第一次调用比后续调用花费的时间更长，因为必须有一些开销用于初始化计时器。
希望以下修改后的基准能够给问题带来一些新的启示：
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 并不严格一致，但整体计时功能非常相似。
只是在 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)
所以两个版本同时出现，这正是我所期望的这样一个简单的操作，在这种操作中，正确的类型推断代码应该编译成大致相同的机器代码。
此外，这里有一个简单的 python 脚本，它使用 numpy 进行矢量化操作：
from timeit import default_timer
import numpy as np
s = 15000
start = default_timer()
# for-loop
H = np.zeros([s, s])
for c in range(1, s):
for r in range(1, s):
H[r, c] = 1 / (r + c - 1)
end = default_timer()
print(end - start)
start = default_timer()
# vectorized
c = np.arange(1, s).reshape([1, -1])
r = c.T
HH = 1 / (c + r - 1)
end = default_timer()
print(end - start)
for 循环：32.94566780002788 秒
矢量化：0.494859800033737 秒
虽然 for 循环版本非常慢，但矢量化版本比发布的 fortran/julia 倍快。 Numpy 在内部尝试使用特殊的 SIMD 硬件指令来加速向量运算，这可能会产生重大影响。 fortran/julia 编译器可能无法从提供的代码生成这些指令，但 numpy/matlab 可以。然而，Matlab 仍然比 numpy 代码快 10 倍，我认为这不能通过更好地使用 SIMD 指令来解释。相反，他们也可能使用多个线程来并行计算，因为矩阵相当大。
最终，我认为 matlab 数字是合理的，但我不确定它们是如何获得加速的。
最初这是一个问题in mathematica.SE，但由于讨论涉及多种编程语言，我认为最好改写一下，post 在这里。
简而言之，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
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的时序靠谱吗？
如果第一个问题的答案是肯定的，那么我们如何用 julia 重现性能（对于 2 GHz CPU，时间应该在
0.05秒左右） 、fortran 或任何其他编程语言？
tic/toc 应该没问题，但看起来时间因内存预分配而出现偏差。
我可以重现与您的 MATLAB 示例类似的时序，但是
首先 运行（
clear工作区）- 循环方法需要 2.08 秒
- 矢量化方法需要 1.04 秒
- 矢量化节省了 50% 的执行时间
第二个 运行（工作区未清除）
- 循环方法需要 2.55 秒
- 矢量化方法需要 0.065 秒
- 矢量化“节省”了 97.5% 的执行时间
我的猜测是，由于循环方法通过 zeros 显式地创建了一个新矩阵，内存在每个 运行 上都从头开始重新分配，并且您看不到随后的速度改进运行s.
然而，当 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 之间的输出大小差异非常小，但执行时间 returns 与向量化方法的测试 1 中的执行时间相似，表明重新分配以创建 HH花费大部分时间。
脚注：tic/toc MATLAB 中的时序可以通过使用内置的 timeit 函数来改进，该函数本质上是取几个 运行 的平均值.从 timeit 的工作中观察到的一件有趣的事情是，它通过多次调用 tic/toc 函数来明确地“预热”（引用评论）函数。您可以从清晰的工作区（没有中间代码）中多次 运行ning tic/toc 时看到第一次调用比后续调用花费的时间更长，因为必须有一些开销用于初始化计时器。
希望以下修改后的基准能够给问题带来一些新的启示：
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 并不严格一致，但整体计时功能非常相似。
只是在 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)
所以两个版本同时出现，这正是我所期望的这样一个简单的操作，在这种操作中，正确的类型推断代码应该编译成大致相同的机器代码。
此外，这里有一个简单的 python 脚本，它使用 numpy 进行矢量化操作：
from timeit import default_timer
import numpy as np
s = 15000
start = default_timer()
# for-loop
H = np.zeros([s, s])
for c in range(1, s):
for r in range(1, s):
H[r, c] = 1 / (r + c - 1)
end = default_timer()
print(end - start)
start = default_timer()
# vectorized
c = np.arange(1, s).reshape([1, -1])
r = c.T
HH = 1 / (c + r - 1)
end = default_timer()
print(end - start)
for 循环：32.94566780002788 秒
矢量化：0.494859800033737 秒
虽然 for 循环版本非常慢，但矢量化版本比发布的 fortran/julia 倍快。 Numpy 在内部尝试使用特殊的 SIMD 硬件指令来加速向量运算，这可能会产生重大影响。 fortran/julia 编译器可能无法从提供的代码生成这些指令，但 numpy/matlab 可以。然而，Matlab 仍然比 numpy 代码快 10 倍，我认为这不能通过更好地使用 SIMD 指令来解释。相反，他们也可能使用多个线程来并行计算，因为矩阵相当大。
最终，我认为 matlab 数字是合理的，但我不确定它们是如何获得加速的。