如何在matlab中对双重求和进行向量化

Question

我正在尝试为这个等式生成矢量化代码(以及同一形式的多个其他代码):

这将被评估,ode45并且唯一的变化部分是Rmn(t)我预先计算贝塞尔和sin函数.目前我的代码看起来像这样:

for iR = 1:nR
    p(:, iR) = sum(Rmn.*J0m(:, :, iR),1)*sinanz;
end

M,N是我总和中的术语数R,Z是我正在使用的数量r和z坐标.Rmn是一个M*N矩阵.J0m是一个M*N*R数组.这是一个M*R重复的矩阵N.sinanz是一个N*Z矩阵.J0m并sinanz预先计算,不要改变.

这工作但很慢,所以我试图优化它.我认为第一步是减少J0m所以它只是m*R但我无法弄清楚如何.我正在寻找有关如何执行此操作的任何建议以及有关如何优化此操作的任何建议.

Answer 1

您可能已经知道,您应该p在循环之前预先分配:

p = zeros(Z, nR);

这可以防止数组p在每次迭代中增长,从而极大地加速循环.

你可以通过以下方式对整个事物进行矢量化bsxfun:

% C ? M×N×R array of all products Rmn·J0m
C = bsxfun(@times, Rmn, J0m);

% D ? M×N×R×Z array of all products C·sinanz
D = bsxfun(@times, C, permute(sinanz, [3 1 4 2]));

% Sum in M and N directions, and re-order
p = permute(sum(sum(D,1),2), [4 3 1 2]);

但我怀疑它会更快; MATLAB(读取:BLAS)与2D矩阵相当快,但对于更多D阵列通常不是很好.

我建议你读一下bsxfun; 它也是以你描述的方式将J0m数组减少到M× R的方法.

当然,您可以permute通过正确定义变量来摆脱s,所以让我们在循环代码和矢量化代码的"理想"版本中进行一个小测试:

%% Initialize some variables

% Number of tests
TESTS = 1e4;

% Your dimensions
M = 5;   nR = 4;
N = 2;   Z  = 3;

% Some dummy data
Rmn    = rand(M,N);
sinanz = rand(N,Z);
J0m    = rand(M,nR); % NOTE: J0m doesn't have to be 3D anymore by using bsxfun

%% Test 1: your own version, optimized

% Start test
start = tic;

% Repeat the test a couple of times to get a good average
for ii = 1:TESTS
    p1 = zeros(Z,nR);
    for iR = 1:nR
        p1(:, iR) = sum( bsxfun(@times,Rmn,J0m(:, iR)), 1 )*sinanz;
    end    
end
stop = toc(start);

% Average execution time
fprintf(1, 'Average time of looped version: %f seconds.\n', stop/TESTS);

%% Vectorized version, using 4D arrays: 

% Don't make the permutes part of the test
J0m    = permute(J0m, [1 3 2]);
sinanz = permute(sinanz, [3 1 4 2]);

% Start test
start = tic;

% Repeat the test a couple of times to get a good average
for ii = 1:TESTS
    C  = bsxfun(@times, Rmn, J0m);
    D  = bsxfun(@times, C, sinanz);
    p2 = sum(sum(D,1),2);

end
stop = toc(start);

% Average execution time
fprintf(1, 'Average time of vectorized version: %f seconds.\n', stop/TESTS);

%% Check for equality

maxDifference = max(max(p1 - squeeze(p2)'))

结果:

Average time of looped version    : 0.000054 seconds.
Average time of vectorized version: 0.000023 seconds.
maxDifference =
    4.440892098500626e-016

这看起来还不错!但是,使用

M = 50;   nR = 40;
N = 20;   Z  = 30;

你得到

Average time of looped version    : 0.000708 seconds.
Average time of vectorized version: 0.009835 seconds.
maxDifference =
    8.526512829121202e-014

因此矢量化版本比循环变体慢一个数量级.

当然,您的里程可能会有所不同,但带回家的信息是,您应该预计这种差异会随着尺寸的增加而变得越来越糟.

所以,最后:

• 如果尺寸很大,请坚持使用.如果他们是小,也使用了循环-它只是SO对眼睛比矢量版本更容易:)
• 确保预先分配p变量
• 使用bsxfun来减少内存占用(但不会提高速度远)