有一个更快速的通过空矩阵乘法初始化数组的方法吗？（Matlab）
有一个更快速的通过空矩阵乘法初始化数组的方法吗？（Matlab）
我偶然发现 Matlab 在处理空矩阵时的奇怪方式（在我看来）。例如，如果两个空矩阵相乘，结果是：
zeros(3,0)*zeros(0,3) ans = 0 0 0 0 0 0 0 0 0
这使我非常惊讶，但是快速搜索让我找到了上述链接，对此有了解释。
然而，没有什么准备能让我接受下面这个观察结果。我问自己，这种乘法相比于只使用 zeros(n) 函数来说，效率有多高，比如说用于初始化？我使用 timeit 得出了答案：
f=@() zeros(1000)
timeit(f)
ans =
0.0033
与：
g=@() zeros(1000,0)*zeros(0,1000)
timeit(g)
ans =
9.2048e-06
都得出了相同的结果，即一个由 double 类型的零组成的 1000x1000 矩阵，但是空矩阵乘法的速度快了约 350 倍！（使用 tic、toc 和一个循环也会得到类似的结果）
这怎么可能？是 timeit 或者 tic,toc 有所欺骗，还是我发现了一种更快的初始化矩阵的方法？（这是在 Matlab 2012a 上、在 win7-64 机器上、使用 intel-i5 650 3.2Ghz 进行的）
编辑：
在阅读了你们的反馈后，我更仔细地研究了这个奇怪的现象，在两台不同的计算机上进行了测试（虽然 Matlab 版本相同，都是 2012a）。测试代码考察了运行时间和矩阵 n 的大小之间的关系，结果如下：
![\"enter](\"https://i.stack.imgur.com/DSZn7.png\")
生成这个代码的方法与之前一样，使用timeit。但是，使用tic和toc的循环看起来是相同的。因此，对于小的尺寸，zeros(n)是可比的。然而，在n=400左右，空矩阵乘法的性能有一个跳跃。我用来生成这个图的代码是：
n=unique(round(logspace(0,4,200)));
for k=1:length(n)
f=@() zeros(n(k));
t1(k)=timeit(f);
g=@() zeros(n(k),0)*zeros(0,n(k));
t2(k)=timeit(g);
end
loglog(n,t1,'b',n,t2,'r');
legend('zeros(n)','zeros(n,0)*zeros(0,n)',2);
xlabel('matrix size (n)'); ylabel('time [sec]');
您有没有经历过这种情况？
编辑#2：
顺便说一句，不需要使用空矩阵乘法即可获得此效果。可以简单地执行以下操作：
z(n,n)=0;
其中n>前面图表中看到的某个阈值矩阵大小，并使用timeit获得与空矩阵乘法完全相同的效率曲线。
![\"enter](\"https://i.stack.imgur.com/EZVOI.png\")
这是一个改善代码效率的例子：
n = 1e4;
clear z1
tic
z1 = zeros( n );
for cc = 1 : n
z1(:,cc)=cc;
end
toc % Elapsed time is 0.445780 seconds.
%%
clear z0
tic
z0 = zeros(n,0)*zeros(0,n);
for cc = 1 : n
z0(:,cc)=cc;
end
toc % Elapsed time is 0.297953 seconds.
但是，使用z(n,n)=0;可以得到与zeros(n)相似的结果。
结果可能有点误导性。当你将两个空矩阵相乘时，生成的矩阵并不会立即“分配”和“初始化”，而是在你第一次使用它时延迟分配（有点像延迟评估）。
同样地，当你在索引超出范围以扩展变量时（像增长），如果是数值数组就会把任何缺失的条目填充为零（我会在后面讨论非数值情况）。当然，以这种方式增大矩阵并不会覆盖现有元素。
因此，尽管它可能看起来更快，但你只是推迟了分配时间直到你实际上第一次使用矩阵。最终，你将拥有类似于从一开始就进行分配的时间测量。
以下示例展示了这种行为，与一些其他替代方案进行比较：
N = 1000; clear z tic, z = zeros(N,N); toc tic, z = z + 1; toc assert(isequal(z,ones(N))) clear z tic, z = zeros(N,0)*zeros(0,N); toc tic, z = z + 1; toc assert(isequal(z,ones(N))) clear z tic, z(N,N) = 0; toc tic, z = z + 1; toc assert(isequal(z,ones(N))) clear z tic, z = full(spalloc(N,N,0)); toc tic, z = z + 1; toc assert(isequal(z,ones(N))) clear z tic, z(1:N,1:N) = 0; toc tic, z = z + 1; toc assert(isequal(z,ones(N))) clear z val = 0; tic, z = val(ones(N)); toc tic, z = z + 1; toc assert(isequal(z,ones(N))) clear z tic, z = repmat(0, [N N]); toc tic, z = z + 1; toc assert(isequal(z,ones(N)))
结果表明，如果你将每种情况下两个指令的经过时间相加，最终得到的总时间测量相似：
// zeros(N,N) Elapsed time is 0.004525 seconds. Elapsed time is 0.000792 seconds. // zeros(N,0)*zeros(0,N) Elapsed time is 0.000052 seconds. Elapsed time is 0.004365 seconds. // z(N,N) = 0 Elapsed time is 0.000053 seconds. Elapsed time is 0.004119 seconds.
其他的时间测量是：
// full(spalloc(N,N,0)) Elapsed time is 0.001463 seconds. Elapsed time is 0.003751 seconds. // z(1:N,1:N) = 0 Elapsed time is 0.006820 seconds. Elapsed time is 0.000647 seconds. // val(ones(N)) Elapsed time is 0.034880 seconds. Elapsed time is 0.000911 seconds. // repmat(0, [N N]) Elapsed time is 0.001320 seconds. Elapsed time is 0.003749 seconds.
这些测量的时间非常短，可能不太准确，所以你可能需要循环运行这些命令几千次，并取平均数。有时，运行已保存的M函数比在脚本或命令行提示符中运行要快，因为只有这种方法才会发生某些优化……
无论何种情况，分配通常只需进行一次，因此多花30ms无关紧要 🙂
对于单元格数组或结构数组，可以看到类似的行为。考虑以下例子：
N = 1000;
tic, a = cell(N,N); toc
tic, b = repmat({[]}, [N,N]); toc
tic, c{N,N} = []; toc
它输出的结果是：
Elapsed time is 0.001245 seconds. Elapsed time is 0.040698 seconds. Elapsed time is 0.004846 seconds.
请注意，即使它们全部相等，它们也会占用不同数量的内存：
>> assert(isequal(a,b,c)) >> whos a b c Name Size Bytes Class Attributes a 1000x1000 8000000 cell b 1000x1000 112000000 cell c 1000x1000 8000104 cell
实际上，此处情况要稍微复杂一些，因为MATLAB可能是在所有单元格中共享同一个空矩阵，而不是创建多个副本。
单元格数组a实际上是未初始化单元格的数组（一个空指针数组），而b是一个单元格数组，其中每个单元格都是一个空数组[]（通过数据共享，只有第一个单元格b{1}指向[]，所有其他单元格都引用第一个单元格）。最终的数组c类似于a（未初始化单元格），但最后一个单元格包含一个空的数值矩阵[]。
我查看了libmx.dll导出的C函数列表（使用Dependency Walker工具），发现了一些有趣的事情。
- 有未记录的函数可以创建未初始化的数组：mxCreateUninitDoubleMatrix、mxCreateUninitNumericArray和mxCreateUninitNumericMatrix。实际上，File Exchange上有一个提交利用这些函数来提供一个更快的替代zeros函数的方法。
- 存在一个未记录的函数称为mxFastZeros。在网上搜索时，我发现你也在MATLAB Answers上发布了这个问题，并获得了一些很好的答案。James Tursa（上面提到的UNINIT的同一作者）给出了如何使用这个未记录的函数的示例。
- libmx.dll与tbbmalloc.dll共享库链接。这是Intel TBB可扩展内存分配器。该库提供优化并行应用程序的等效内存分配函数（malloc、calloc、free）。请记住，许多MATLAB函数是自动多线程的，因此当矩阵大小足够大时（这里是Loren Shure最近的评论）zeros（..）可能是多线程的，并且正在使用Intel的内存分配器。
关于内存分配器的最后一点，你可以编写类似于@PavanYalamanchili所做的C/C++类似的基准测试，并比较可用的各种分配器。例如像这个。请记住，MEX文件具有稍微更高的内存管理开销，因为MATLAB会自动使用mxCalloc、mxMalloc或mxRealloc函数释放MEX文件中分配的任何内存。值得一提的是，在旧版本中还可以更改内部内存管理器。
编辑：
下面是一个更全面的基准测试，用于比较所讨论的替代方案。它特别表明，一旦你强调使用整个分配的矩阵，所有三种方法都处于平等的地位，差异微不足道。
function compare_zeros_init()
iter = 100;
for N = 512.*(1:8)
% ZEROS(N,N)
t = zeros(iter,3);
for i=1:iter
clear z
tic, z = zeros(N,N); t(i,1) = toc;
tic, z(:) = 9; t(i,2) = toc;
tic, z = z + 1; t(i,3) = toc;
end
fprintf('N = %4d, ZEROS = %.9f\n', N, mean(sum(t,2)))
% z(N,N)=0
t = zeros(iter,3);
for i=1:iter
clear z
tic, z(N,N) = 0; t(i,1) = toc;
tic, z(:) = 9; t(i,2) = toc;
tic, z = z + 1; t(i,3) = toc;
end
fprintf('N = %4d, GROW = %.9f\n', N, mean(sum(t,2)))
% ZEROS(N,0)*ZEROS(0,N)
t = zeros(iter,3);
for i=1:iter
clear z
tic, z = zeros(N,0)*zeros(0,N); t(i,1) = toc;
tic, z(:) = 9; t(i,2) = toc;
tic, z = z + 1; t(i,3) = toc;
end
fprintf('N = %4d, MULT = %.9f\n\n', N, mean(sum(t,2)))
end
end
下面是按矩阵大小递增的时间平均值。我在R2013a中进行了测试。
>> compare_zeros_init N = 512, ZEROS = 0.001560168 N = 512, GROW = 0.001479991 N = 512, MULT = 0.001457031 N = 1024, ZEROS = 0.005744873 N = 1024, GROW = 0.005352638 N = 1024, MULT = 0.005359236 N = 1536, ZEROS = 0.011950846 N = 1536, GROW = 0.009051589 N = 1536, MULT = 0.008418878 N = 2048, ZEROS = 0.012154002 N = 2048, GROW = 0.010996315 N = 2048, MULT = 0.011002169 N = 2560, ZEROS = 0.017940950 N = 2560, GROW = 0.017641046 N = 2560, MULT = 0.017640323 N = 3072, ZEROS = 0.025657999 N = 3072, GROW = 0.025836506 N = 3072, MULT = 0.051533432 N = 3584, ZEROS = 0.074739924 N = 3584, GROW = 0.070486857 N = 3584, MULT = 0.072822335 N = 4096, ZEROS = 0.098791732 N = 4096, GROW = 0.095849788 N = 4096, MULT = 0.102148452
很奇怪，我看到f函数比g函数更快，但是对于我来说它们都是相同的。也许是MATLAB的不同版本？
>> g = @() zeros(1000, 0) * zeros(0, 1000);
>> f = @() zeros(1000)
f =
@()zeros(1000)
>> timeit(f)
ans =
8.5019e-04
>> timeit(f)
ans =
8.4627e-04
>> timeit(g)
ans =
8.4627e-04
编辑，你可以再f和g函数结尾处加上+1，看看你得到的时间吗。
编辑Jan 6，2013 7:42 EST
我是远程使用计算机的，所以很抱歉图片质量不好（不得不盲目生成它们）。
计算机配置：
i7 920. 2.653 GHz. Linux. 12 GB RAM. 8MB 缓存。
看起来，即使我可以访问的机器也显示出这种行为，只是大小更大一些（在1979和2073之间）。我现在想不出一个空矩阵乘法在更大的大小时为什么会更快。
我会进行更多的调查，然后回来。
编辑Jan 11，2013
在@EitanT的帖子之后，我想要进行更多的挖掘。我写了一些C代码，以查看MATLAB可能如何创建零矩阵。以下是我使用的C++代码。
int main(int argc, char **argv)
{
for (int i = 1975; i <= 2100; i+=25) {
timer::start();
double *foo = (double *)malloc(i * i * sizeof(double));
for (int k = 0; k < i * i; k++) foo[k] = 0;
double mftime = timer::stop();
free(foo);
timer::start();
double *bar = (double *)malloc(i * i * sizeof(double));
memset(bar, 0, i * i * sizeof(double));
double mmtime = timer::stop();
free(bar);
timer::start();
double *baz = (double *)calloc(i * i, sizeof(double));
double catime = timer::stop();
free(baz);
printf("%d, %lf, %lf, %lf\n", i, mftime, mmtime, catime);
}
}
以下是结果。
$ ./test 1975, 0.013812, 0.013578, 0.003321 2000, 0.014144, 0.013879, 0.003408 2025, 0.014396, 0.014219, 0.003490 2050, 0.014732, 0.013784, 0.000043 2075, 0.015022, 0.014122, 0.000045 2100, 0.014606, 0.014480, 0.000045
正如你所看到的，calloc（第4列）似乎是最快的方法。它在2025和2050之间也变得更快了（我假设它会在2048左右？）。
现在我回到MATLAB去检查同样的事情。以下是结果。
>> test 1975, 0.003296, 0.003297 2000, 0.003377, 0.003385 2025, 0.003465, 0.003464 2050, 0.015987, 0.000019 2075, 0.016373, 0.000019 2100, 0.016762, 0.000020
看起来f()和g()都在较小的大小时使用calloc（<2048？）。但是在较大的大小时，f()（zeros(m，n)）开始使用malloc+ memset，而g()（zeros（m，0）* zeros（0，n））仍然使用calloc。
因此，以下是解释分歧的原因：
zeros(..)在更大的尺寸下开始使用不同（更慢？）的方案。calloc的行为也有些出乎意料，并带来了性能的改善。
这是在Linux上的行为。有没有人可以在另一台机器上（也许是不同的操作系统）进行同样的实验，看看实验是否成立？