是否有可能使这个向量化的 MATLAB 代码更快？

Question

在我的程序中，我需要计算总和：

.

我用 C 和 z 的新值计算了这个总和 2500 次。

参数z 可能是一个向量。我写了简单的for循环和矢量化版本代码如下：

K = 200;
n_z = 40000;
C = ones(K,1); % an example, in real life the arey some coefficients (at next call will be new)
k = 0:K-1;
z = linspace(0, 2*pi, n_z); % at next call will be new

tic;
    my_sum_for = zeros(1, K);
    for i=1:n_z
       my_sum_for(i) = C' * tan(k' * z(i));
    end
toc; % Elapsed time is 1.820485 seconds.

tic;
     my_sum = C' * tan(k' * z);
toc; % Elapsed time is 0.160924 seconds.

矢量化版本更快，但还不够。 是否可以改进矢量化版本？

在 Dominique Jacquel 的回答之后，我有了这个矢量化版本，它更快：

    K = 200;
    n_z = 40000;
    C = ones(K,1)';  % an example, in real life they are some coefficients (at next call will be new)
    k = (0:K-1)';
    z = linspace(0, 2*pi, n_z);  % at next call will be new

    tic;
        my_sum_for = zeros(1, K);
        for i=1:n_z
           my_sum_for(i) = C * tan(k * z(i));
        end
    toc; % Elapsed time is 1.521587 seconds.

    tic;
         my_sum = C * tan(k * z);
    toc; % Elapsed time is 0.125468 seconds.

是否有可能进一步改进矢量化版本（bsxfun、arrayfun 之类的）？ 250 秒的时间对我来说仍然很慢（占所有计算的 75%）。

Answer 1

我认为您已经非常接近这里的硬件限制了。 Matlab 中的矩阵乘法是使用BLAS libraries 完成的，在性能方面已被证明难以击败。

AFAIK，正切函数有实际的dedicated hardware 来计算它的值。此外，Matlab 会自动将大型矩阵的三角函数分布到多个内核上，因此基本上没有什么需要改进的地方。

另外，如果我错了，请纠正我，但考虑到可能的数据开销和内存问题，我认为这种计算在 GPU 上实际上会更慢。

现在如果你比较这些：

tic;
for ii = 1:10
    my_sum = C * tan(k * z);
end
toc 

tic;
for ii = 1:10
    my_sum_notan = C * k * z;
end
toc

你会看到所有的痛苦都来自于正切函数，所以你最好专注于这一点。正如您所读到的here，基本上只有在牺牲一些准确性的情况下才能加快触发函数的速度。

你应该问自己这些问题：

1. 您是否愿意放弃全双精度，或者说，6 位“足够接近”？

2. 您不能重新制定问题以便之后计算切线吗？还是以前？或者无论如何，元素数量要少得多？在上述问题中，这显然是不可能的，但我不知道完整的代码——可能有一些很好的三角标识可以适用于您的问题。

3. 鉴于以上所有因素，进一步优化所需的工作量真的超过了更长的运行时间吗？与编写实现残缺但快速的触发函数的自定义、难以移植的 MEX 函数相比，250 秒听起来还不错。

Answer 2

尽可能多地预先进行矩阵运算（在这种情况下为转置）以节省循环时间

K = 200;
n_z = 40000;
C = ones(K,1)';
k = (0:K-1)';
z = linspace(0, 2*pi, n_z);

tic;
    my_sum_for = zeros(1, K);
    for i=1:n_z
       my_sum_for(i) = C * tan(k * z(i));
    end
toc

tic;
     my_sum = C * tan(k * z);
toc;

我之前的执行次数

Elapsed time is 1.266158 seconds.
Elapsed time is 0.531173 seconds.

之后

Elapsed time is 0.496803 seconds.
Elapsed time is 0.185396 seconds.

Answer 3

这是一个在我的电脑上运行速度稍快的版本：

k = repmat((0:K-1)', 1, n_z);
z = repmat(linspace(0, 2*pi, n_z), K, 1);
C = ones(1, K);
tic
my_sum = C*tan(k.*z);
toc

本质上，我直接对矩阵进行运算，而不是 k 和 z 的外积。

第一版

Elapsed time is 0.652923 seconds.
Elapsed time is 0.240300 seconds.

在 Dominique Jacquel 的回答之后

Elapsed time is 0.376218 seconds.
Elapsed time is 0.214047 seconds.

我的版本

Elapsed time is 0.168535 seconds.

您可能需要添加 repmats 的成本，但也许您只能这样做一次，我不知道其余的代码。

我完全同意 Rody Oldenhuis 的观点。大部分工作在于切线函数。我可以告诉你更多。 k.*z 计算非常有效，不能有太大改进。如果你计算内存带宽，它在我的电脑上大约为 10GB/s。我能得到的峰值大约是 16GB/s，所以很接近。那里没有多少可能性。与 C*T 相同。那是一个简单的 BLAS2 矩阵向量乘法，它是有内存限制的。对于您所显示的系统大小，MATLAB 开销并不太大。

编辑：正如 Rody 提到的，新版本的 MATLAB 已经对 tan() 进行了并行化。所以这里也不多。

您只能希望改进 tan() - 可能通过并行运行它。毕竟，这是一项非常容易并行化的任务……考虑将其导出到 MEX 文件，该文件将使用 OpenMP。非常简单的工作，如果您有几个备用内核，则可以大大加快速度。