为什么a*b*a在Matlab脚本中使用gpuArray时需要比('*'(a*b)')'更长的时间？

Question

下面的代码执行操作上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)')'速度大约是计算速度的4倍a*b*a.以下是R2018a上Nvidia K20上面脚本的输出(我尝试过不同的版本和不同的GPU,具有相似的行为).

>> test
time for (a'*(a*b)')': 0.43234s
time for a*b*a: 1.7175s
error = 2.0009e-11

甚至更奇怪的是,如果上述脚本的第一行和最后一行是未注释的(它变成一个函数),则两个取较长的时间量(〜1.7S代替〜0.4秒).以下是此案例的输出:

>> test
time for (a'*(a*b)')': 1.717s
time for a*b*a: 1.7153s
error = 1.0914e-11

我想知道是什么导致了这种行为,以及如何在matlab函数内而不是在脚本内部的较短时间内(即~0.4s而不是~1.7s)执行a*b*a或(a'*(a*b)')'两者.

Answer 1

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

在你的上下文中,它是最后一个乘法.演示我用副本c替换a,并将其乘以两次,一次为稀疏,一次为全矩阵.

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()方法将过程分为两阶段乘法,没有其中有两个稀疏矩阵.

Answer 2

我联系了 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”表达式的性能相当；此外，当从 MATLAB 函数文件执行时，“res3”表达式也显示出类似的性能。但请注意，可以使用额外的内存来存储初始数组的转置副本。如果您在系统上发现不同的性能行为，请告诉我，我可以进一步调查此问题。

此外，当使用“gputimeit”函数测量时，“res3”操作也显示出更快的性能。有关详细信息，请参阅随附的“testscript2.m”文件。我还附上了“test_v2.m”，它是您的 Stack Overflow 帖子中“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”函数的方式存在任何具体问题。但请注意，对于 GPU 相关的分析，通常建议使用“gputimeit”而不是直接使用“tic”和“toc”。

我将暂时关闭此案例，但如果您对此还有任何其他疑问，请告诉我。

%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