为什么在 Matlab 脚本中使用 gpuArray 时 a*b*a 比 (a'*(a*b)')' 花费的时间更长？
why does a*b*a take longer than (a'*(a*b)')' when using gpuArray in Matlab scripts?
下面的代码以两种不同的方式对 gpuArrays a 和 b 执行相同的操作。第一部分计算 (a'*(a*b)')' ，而第二部分计算 a*b*a。然后验证结果是一样的。
%function test
clear
rng('default');rng(1);
a=sprand(3000,3000,0.1);
b=rand(3000,3000);
a=gpuArray(a);
b=gpuArray(b);
tic;
c1=gather(transpose(transpose(a)*transpose(a*b)));
disp(['time for (a''*(a*b)'')'': ' , num2str(toc),'s'])
clearvars -except c1
rng('default');
rng(1)
a=sprand(3000,3000,0.1);
b=rand(3000,3000);
a=gpuArray(a);
b=gpuArray(b);
tic;
c2=gather(a*b*a);
disp(['time for a*b*a: ' , num2str(toc),'s'])
disp(['error = ',num2str(max(max(abs(c1-c2))))])
%end
然而，计算 (a'*(a*b)')' 大约比计算 a*b*a 快 4 倍。这是 Nvidia K20 上 R2018a 中上述脚本的输出（我尝试过具有类似行为的不同版本和不同 GPU）。
>> test
time for (a'*(a*b)')': 0.43234s
time for a*b*a: 1.7175s
error = 2.0009e-11
更奇怪的是，如果上面脚本的第一行和最后一行没有注释（将其变成一个函数），那么两者都会花费更长的时间（~1.7s 而不是~0.4s）。下面是这种情况的输出：
>> test
time for (a'*(a*b)')': 1.717s
time for a*b*a: 1.7153s
error = 1.0914e-11
我想知道是什么导致了这种行为，以及如何在更短的时间内执行 a*b*a 或 (a'*(a*b)')' 或两者（即 ~0.4s 而不是 ~1.7s ) 在 matlab 函数中而不是在脚本中。
编辑 2 我可能是对的，请参阅
编辑：他们使用MAGMA，这是专栏专业。我的回答不成立，但我会把它留在这里一段时间，以防它能帮助破解这种奇怪的行为。
下面的答案是错误的
这是我的猜测，如果不知道 MATLAB 的代码，我不能 100% 告诉你。
假设： MATLAB 的并行计算代码使用 CUDA 库，而不是他们自己的。
重要信息
- MATLAB 是列专业，CUDA 是行专业。
- 没有二维矩阵，只有带有 2 个索引的一维矩阵
为什么这很重要？好吧，因为 CUDA 是高度优化的代码，它使用内存结构来最大化每个内核的缓存命中率（GPU 上最慢的操作是读取内存）。这意味着标准的 CUDA 矩阵乘法代码将利用内存读取的顺序来确保它们相邻。但是，行优先的相邻内存不是列优先的。
因此，作为编写软件的人，有 2 个解决方案
- 在 CUDA 中编写您自己的列主代数库
- 从 MATLAB 中取出每个 input/output 并转置它（即从列优先转换为行优先）
他们已经完成了第2点，假设有一个用于MATLAB并行处理工具箱的智能JIT编译器（合理假设），对于第二种情况，需要a和b，转置他们，做数学运算，并在 gather.
时转置输出
然而，在第一种情况下，您已经不需要转置输出，因为它在内部已经转置并且 JIT 捕获了这一点，因此它不会调用 gather(transpose( XX )) 它只是跳过输出转置 is side .与 transpose(a*b) 相同。请注意 transpose(a*b)=transpose(b)*transpose(a)，因此突然不需要转置（它们都在内部被跳过）。转置是一项代价高昂的操作。
确实这里有一件奇怪的事情：把代码变成一个函数突然变慢了。我最好的猜测是，由于 JIT 在不同情况下的行为不同，它不会捕获所有这些 transpose 内部的东西，只是执行所有操作，从而失去速度。
有趣的观察：在我的 PC 中 CPU 与 GPU 花费的时间相同 a*b*a。
GPU 上两个稀疏矩阵的乘法似乎存在问题。全矩阵稀疏的时间比稀疏稀疏快 1000 倍以上。一个简单的例子：
str={'sparse*sparse','sparse*full'};
for ii=1:2
rng(1);
a=sprand(3000,3000,0.1);
b=sprand(3000,3000,0.1);
if ii==2
b=full(b);
end
a=gpuArray(a);
b=gpuArray(b);
tic
c=a*b;
disp(['time for ',str{ii},': ' , num2str(toc),'s'])
end
在您的上下文中，它是最后一个乘法。为了演示我将 a 替换为重复的 c，并乘以它两次，一次作为稀疏矩阵，一次作为完整矩阵。
str={'a*b*a','a*b*full(a)'};
for ii=1:2
%rng('default');
rng(1)
a=sprand(3000,3000,0.1);
b=rand(3000,3000);
rng(1)
c=sprand(3000,3000,0.1);
if ii==2
c=full(c);
end
a=gpuArray(a);
b=gpuArray(b);
c=gpuArray(c);
tic;
c1{ii}=a*b*c;
disp(['time for ',str{ii},': ' , num2str(toc),'s'])
end
disp(['error = ',num2str(max(max(abs(c1{1}-c1{2}))))])
我可能错了，但我的结论是 a * b * a涉及两个稀疏矩阵的乘法（a 和 a 再次）并且没有得到很好的处理，而使用 transpose() 方法将过程分为两个阶段乘法，在none中有两个稀疏矩阵。
我联系了 Mathworks 技术支持，Rylan 终于阐明了这个问题。 （感谢 Rylan！）他的完整回复如下。函数与脚本问题似乎与某些优化有关，matlab 自动应用于未按预期工作的函数（但不是脚本）。
瑞兰的回复：
感谢您在此问题上的耐心等待。我咨询了 MATLAB GPU 计算开发人员以更好地理解这一点。
这个问题是由于MATLAB在遇到矩阵乘法和转置等特定操作时进行的内部优化引起的。在执行 MATLAB 函数（或匿名函数）而不是脚本时，可能会专门启用其中一些优化。
当您的初始代码从脚本执行时，不会执行特定的矩阵转置优化，这导致 'res2' 表达式比 'res1' 表达式更快:
n = 2000;
a=gpuArray(sprand(n,n,0.01));
b=gpuArray(rand(n));
tic;res1=a*b*a;wait(gpuDevice);toc % Elapsed time is 0.884099 seconds.
tic;res2=transpose(transpose(a)*transpose(a*b));wait(gpuDevice);toc % Elapsed time is 0.068855 seconds.
但是，当将上述代码放在 MATLAB 函数文件中时，会进行额外的矩阵转置时间优化，这会导致 'res2' 表达式通过不同的代码路径（并且不同的 CUDA 库函数调用）与从脚本调用的同一行相比。因此，当从函数文件调用时，此优化会为 'res2' 行生成较慢的结果。
为避免函数文件中出现此问题，需要以阻止 MATLAB 应用此优化的方式拆分转置和乘法运算。分隔 'res2' 语句中的每个子句似乎就足够了：
tic;i1=transpose(a);i2=transpose(a*b);res3=transpose(i1*i2);wait(gpuDevice);toc % Elapsed time is 0.066446 seconds.
在上面的行中，'res3' 是从两个中间矩阵生成的：'i1' 和 'i2'。从脚本执行时，性能（在我的系统上）似乎与 'res2' 表达式相当；此外，'res3' 表达式在从 MATLAB 函数文件执行时也显示出类似的性能。但是请注意，可以使用额外的内存来存储初始数组的转置副本。如果您在系统上看到不同的性能行为，请告诉我，我可以进一步调查。
此外，'res3' 操作在使用 'gputimeit' 函数测量时也显示出更快的性能。有关这方面的更多信息，请参阅随附的 'testscript2.m' 文件。我还附上了 'test_v2.m'，这是对 Stack Overflow post.
中 'test.m' 函数的修改
感谢您向我报告此问题。对于由此问题造成的任何不便，我深表歉意。我已经创建了一个内部错误报告来通知 MATLAB 开发人员有关此行为的信息。他们可能会在 MATLAB 的未来版本中为此提供修复。
由于您还有一个关于比较使用 'gputimeit' 与使用 'tic' 和 'toc' 的 GPU 代码性能的其他问题，我只想提供一个MATLAB GPU 计算开发人员之前提到的建议。通常最好在 'tic' 语句之前调用 'wait(gpuDevice)' 以确保前几行的 GPU 操作不会在下一行的测量中重叠。例如，在以下行中：
b=gpuArray(rand(n));
tic; res1=a*b*a; wait(gpuDevice); toc
如果在'tic'之前没有调用'wait(gpuDevice)'，从上一行构造'b'数组所花费的一些时间可能会重叠并得到计入执行 'res1' 表达式所花费的时间。这将是首选：
b=gpuArray(rand(n));
wait(gpuDevice); tic; res1=a*b*a; wait(gpuDevice); toc
除此之外，我没有发现您使用 'tic' 和 'toc' 函数的方式有任何具体问题。但是请注意，通常建议使用 'gputimeit' 而不是直接使用 'tic' 和 'toc' 进行 GPU 相关的分析。
我会继续并暂时关闭此案例，但如果您对此有任何疑问，请告诉我。
%testscript2.m
n = 2000;
a = gpuArray(sprand(n, n, 0.01));
b = gpuArray(rand(n));
gputimeit(@()transpose_mult_fun(a, b))
gputimeit(@()transpose_mult_fun_2(a, b))
function out = transpose_mult_fun(in1, in2)
i1 = transpose(in1);
i2 = transpose(in1*in2);
out = transpose(i1*i2);
end
function out = transpose_mult_fun_2(in1, in2)
out = transpose(transpose(in1)*transpose(in1*in2));
end
.
function test_v2
clear
%% transposed expression
n = 2000;
rng('default');rng(1);
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);
tic;
c1 = gather(transpose( transpose(a) * transpose(a * b) ));
disp(['time for (a''*(a*b)'')'': ' , num2str(toc),'s'])
clearvars -except c1
%% non-transposed expression
rng('default');
rng(1)
n = 2000;
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);
tic;
c2 = gather(a * b * a);
disp(['time for a*b*a: ' , num2str(toc),'s'])
disp(['error = ',num2str(max(max(abs(c1-c2))))])
%% sliced equivalent
rng('default');
rng(1)
n = 2000;
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);
tic;
intermediate1 = transpose(a);
intermediate2 = transpose(a * b);
c3 = gather(transpose( intermediate1 * intermediate2 ));
disp(['time for split equivalent: ' , num2str(toc),'s'])
disp(['error = ',num2str(max(max(abs(c1-c3))))])
end
下面的代码以两种不同的方式对 gpuArrays a 和 b 执行相同的操作。第一部分计算 (a'*(a*b)')' ，而第二部分计算 a*b*a。然后验证结果是一样的。
%function test
clear
rng('default');rng(1);
a=sprand(3000,3000,0.1);
b=rand(3000,3000);
a=gpuArray(a);
b=gpuArray(b);
tic;
c1=gather(transpose(transpose(a)*transpose(a*b)));
disp(['time for (a''*(a*b)'')'': ' , num2str(toc),'s'])
clearvars -except c1
rng('default');
rng(1)
a=sprand(3000,3000,0.1);
b=rand(3000,3000);
a=gpuArray(a);
b=gpuArray(b);
tic;
c2=gather(a*b*a);
disp(['time for a*b*a: ' , num2str(toc),'s'])
disp(['error = ',num2str(max(max(abs(c1-c2))))])
%end
然而，计算 (a'*(a*b)')' 大约比计算 a*b*a 快 4 倍。这是 Nvidia K20 上 R2018a 中上述脚本的输出（我尝试过具有类似行为的不同版本和不同 GPU）。
>> test
time for (a'*(a*b)')': 0.43234s
time for a*b*a: 1.7175s
error = 2.0009e-11
更奇怪的是，如果上面脚本的第一行和最后一行没有注释（将其变成一个函数），那么两者都会花费更长的时间（~1.7s 而不是~0.4s）。下面是这种情况的输出：
>> test
time for (a'*(a*b)')': 1.717s
time for a*b*a: 1.7153s
error = 1.0914e-11
我想知道是什么导致了这种行为，以及如何在更短的时间内执行 a*b*a 或 (a'*(a*b)')' 或两者（即 ~0.4s 而不是 ~1.7s ) 在 matlab 函数中而不是在脚本中。
编辑 2 我可能是对的，请参阅
编辑：他们使用MAGMA，这是专栏专业。我的回答不成立，但我会把它留在这里一段时间，以防它能帮助破解这种奇怪的行为。
下面的答案是错误的
这是我的猜测，如果不知道 MATLAB 的代码，我不能 100% 告诉你。
假设： MATLAB 的并行计算代码使用 CUDA 库，而不是他们自己的。
重要信息
- MATLAB 是列专业，CUDA 是行专业。
- 没有二维矩阵，只有带有 2 个索引的一维矩阵
为什么这很重要？好吧，因为 CUDA 是高度优化的代码，它使用内存结构来最大化每个内核的缓存命中率（GPU 上最慢的操作是读取内存）。这意味着标准的 CUDA 矩阵乘法代码将利用内存读取的顺序来确保它们相邻。但是，行优先的相邻内存不是列优先的。
因此，作为编写软件的人，有 2 个解决方案
- 在 CUDA 中编写您自己的列主代数库
- 从 MATLAB 中取出每个 input/output 并转置它（即从列优先转换为行优先）
他们已经完成了第2点，假设有一个用于MATLAB并行处理工具箱的智能JIT编译器（合理假设），对于第二种情况，需要a和b，转置他们，做数学运算，并在 gather.
然而，在第一种情况下，您已经不需要转置输出，因为它在内部已经转置并且 JIT 捕获了这一点，因此它不会调用 gather(transpose( XX )) 它只是跳过输出转置 is side .与 transpose(a*b) 相同。请注意 transpose(a*b)=transpose(b)*transpose(a)，因此突然不需要转置（它们都在内部被跳过）。转置是一项代价高昂的操作。
确实这里有一件奇怪的事情：把代码变成一个函数突然变慢了。我最好的猜测是，由于 JIT 在不同情况下的行为不同，它不会捕获所有这些 transpose 内部的东西，只是执行所有操作，从而失去速度。
有趣的观察：在我的 PC 中 CPU 与 GPU 花费的时间相同 a*b*a。
GPU 上两个稀疏矩阵的乘法似乎存在问题。全矩阵稀疏的时间比稀疏稀疏快 1000 倍以上。一个简单的例子：
str={'sparse*sparse','sparse*full'};
for ii=1:2
rng(1);
a=sprand(3000,3000,0.1);
b=sprand(3000,3000,0.1);
if ii==2
b=full(b);
end
a=gpuArray(a);
b=gpuArray(b);
tic
c=a*b;
disp(['time for ',str{ii},': ' , num2str(toc),'s'])
end
在您的上下文中，它是最后一个乘法。为了演示我将 a 替换为重复的 c，并乘以它两次，一次作为稀疏矩阵，一次作为完整矩阵。
str={'a*b*a','a*b*full(a)'};
for ii=1:2
%rng('default');
rng(1)
a=sprand(3000,3000,0.1);
b=rand(3000,3000);
rng(1)
c=sprand(3000,3000,0.1);
if ii==2
c=full(c);
end
a=gpuArray(a);
b=gpuArray(b);
c=gpuArray(c);
tic;
c1{ii}=a*b*c;
disp(['time for ',str{ii},': ' , num2str(toc),'s'])
end
disp(['error = ',num2str(max(max(abs(c1{1}-c1{2}))))])
我可能错了，但我的结论是 a * b * a涉及两个稀疏矩阵的乘法（a 和 a 再次）并且没有得到很好的处理，而使用 transpose() 方法将过程分为两个阶段乘法，在none中有两个稀疏矩阵。
我联系了 Mathworks 技术支持，Rylan 终于阐明了这个问题。 （感谢 Rylan！）他的完整回复如下。函数与脚本问题似乎与某些优化有关，matlab 自动应用于未按预期工作的函数（但不是脚本）。
瑞兰的回复：
感谢您在此问题上的耐心等待。我咨询了 MATLAB GPU 计算开发人员以更好地理解这一点。
这个问题是由于MATLAB在遇到矩阵乘法和转置等特定操作时进行的内部优化引起的。在执行 MATLAB 函数（或匿名函数）而不是脚本时，可能会专门启用其中一些优化。
当您的初始代码从脚本执行时，不会执行特定的矩阵转置优化，这导致 'res2' 表达式比 'res1' 表达式更快:
n = 2000;
a=gpuArray(sprand(n,n,0.01));
b=gpuArray(rand(n));
tic;res1=a*b*a;wait(gpuDevice);toc % Elapsed time is 0.884099 seconds.
tic;res2=transpose(transpose(a)*transpose(a*b));wait(gpuDevice);toc % Elapsed time is 0.068855 seconds.
但是，当将上述代码放在 MATLAB 函数文件中时，会进行额外的矩阵转置时间优化，这会导致 'res2' 表达式通过不同的代码路径（并且不同的 CUDA 库函数调用）与从脚本调用的同一行相比。因此，当从函数文件调用时，此优化会为 'res2' 行生成较慢的结果。
为避免函数文件中出现此问题，需要以阻止 MATLAB 应用此优化的方式拆分转置和乘法运算。分隔 'res2' 语句中的每个子句似乎就足够了：
tic;i1=transpose(a);i2=transpose(a*b);res3=transpose(i1*i2);wait(gpuDevice);toc % Elapsed time is 0.066446 seconds.
在上面的行中，'res3' 是从两个中间矩阵生成的：'i1' 和 'i2'。从脚本执行时，性能（在我的系统上）似乎与 'res2' 表达式相当；此外，'res3' 表达式在从 MATLAB 函数文件执行时也显示出类似的性能。但是请注意，可以使用额外的内存来存储初始数组的转置副本。如果您在系统上看到不同的性能行为，请告诉我，我可以进一步调查。
此外，'res3' 操作在使用 'gputimeit' 函数测量时也显示出更快的性能。有关这方面的更多信息，请参阅随附的 'testscript2.m' 文件。我还附上了 'test_v2.m'，这是对 Stack Overflow post.
中 'test.m' 函数的修改感谢您向我报告此问题。对于由此问题造成的任何不便，我深表歉意。我已经创建了一个内部错误报告来通知 MATLAB 开发人员有关此行为的信息。他们可能会在 MATLAB 的未来版本中为此提供修复。
由于您还有一个关于比较使用 'gputimeit' 与使用 'tic' 和 'toc' 的 GPU 代码性能的其他问题，我只想提供一个MATLAB GPU 计算开发人员之前提到的建议。通常最好在 'tic' 语句之前调用 'wait(gpuDevice)' 以确保前几行的 GPU 操作不会在下一行的测量中重叠。例如，在以下行中：
b=gpuArray(rand(n));
tic; res1=a*b*a; wait(gpuDevice); toc
如果在'tic'之前没有调用'wait(gpuDevice)'，从上一行构造'b'数组所花费的一些时间可能会重叠并得到计入执行 'res1' 表达式所花费的时间。这将是首选：
b=gpuArray(rand(n));
wait(gpuDevice); tic; res1=a*b*a; wait(gpuDevice); toc
除此之外，我没有发现您使用 'tic' 和 'toc' 函数的方式有任何具体问题。但是请注意，通常建议使用 'gputimeit' 而不是直接使用 'tic' 和 'toc' 进行 GPU 相关的分析。
我会继续并暂时关闭此案例，但如果您对此有任何疑问，请告诉我。
%testscript2.m
n = 2000;
a = gpuArray(sprand(n, n, 0.01));
b = gpuArray(rand(n));
gputimeit(@()transpose_mult_fun(a, b))
gputimeit(@()transpose_mult_fun_2(a, b))
function out = transpose_mult_fun(in1, in2)
i1 = transpose(in1);
i2 = transpose(in1*in2);
out = transpose(i1*i2);
end
function out = transpose_mult_fun_2(in1, in2)
out = transpose(transpose(in1)*transpose(in1*in2));
end
.
function test_v2
clear
%% transposed expression
n = 2000;
rng('default');rng(1);
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);
tic;
c1 = gather(transpose( transpose(a) * transpose(a * b) ));
disp(['time for (a''*(a*b)'')'': ' , num2str(toc),'s'])
clearvars -except c1
%% non-transposed expression
rng('default');
rng(1)
n = 2000;
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);
tic;
c2 = gather(a * b * a);
disp(['time for a*b*a: ' , num2str(toc),'s'])
disp(['error = ',num2str(max(max(abs(c1-c2))))])
%% sliced equivalent
rng('default');
rng(1)
n = 2000;
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);
tic;
intermediate1 = transpose(a);
intermediate2 = transpose(a * b);
c3 = gather(transpose( intermediate1 * intermediate2 ));
disp(['time for split equivalent: ' , num2str(toc),'s'])
disp(['error = ',num2str(max(max(abs(c1-c3))))])
end