矢量化 - Sum 和 Bessel 函数

Question

谁能帮助向量化这个 Matlab 代码？具体问题是向量输入的求和和贝塞尔函数。 谢谢！

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];
for ii = 1:length(rho_g)
    for jj = 1:length(phi_g)
        % Coordinates
        rho_o = rho_g(ii);
        phi_o = phi_g(jj);
        % factors
        fc = cos(n.*(phi_o-phi_s));
        fs = sin(n.*(phi_o-phi_s));

        Ez_t(ii,jj) = sum(tau.*besselj(n,k(3)*rho_s).*besselh(n,2,k(3)*rho_o).*fc);
    end
end

Answer 1

您可以尝试对这段代码进行矢量化，这可能通过一些bsxfun 左右来实现，但是很难理解代码，而且它是否会运行得更快是个问题，因为您的代码已经使用了矢量内循环中的数学（即使您的向量只有长度 3）。生成的代码将变得非常难以阅读，因此当您在 2 年后查看它时，您或您的同事将不知道它的作用。

在将时间浪费在矢量化上之前，更重要的是了解loop invariant code motion，它很容易应用于您的代码。一些观察：

• 你不使用fs，所以删除它。
• 术语tau.*besselj(n,k(3)*rho_s) 不依赖于任何循环变量ii 和jj，因此它是常量。在循环之前计算一次。
• 您可能应该预先分配矩阵Ez_t。
• 在循环期间唯一改变的项是fc，它依赖于jj，和besselh(n,2,k(3)*rho_o)，它依赖于ii。我猜后者的计算时间要多得多，所以最好不要在内循环中计算 N*N 次，而在外循环中只计算 N 次。如果基于jj 的计算需要更多时间，您可以在ii 和jj 上交换for 循环，但这里似乎并非如此。

结果代码看起来像这样（未经测试）：

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

% constant part, does not depend on ii and jj, so calculate only once!
temp1 = tau.*besselj(n,k(3)*rho_s); 

Ez_t = nan(length(rho_g), length(phi_g)); % preallocate space
for ii = 1:length(rho_g)
    % calculate stuff that depends on ii only
    rho_o = rho_g(ii);
    temp2 = besselh(n,2,k(3)*rho_o);

    for jj = 1:length(phi_g)
        phi_o = phi_g(jj);
        fc = cos(n.*(phi_o-phi_s));
        Ez_t(ii,jj) = sum(temp1.*temp2.*fc);
    end
end

Answer 2

初始化-

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

嵌套循环形式（从您的代码中复制并在此处显示仅供比较）-

for ii = 1:length(rho_g)
    for jj = 1:length(phi_g)
        % Coordinates
        rho_o = rho_g(ii);
        phi_o = phi_g(jj);
        % factors
        fc = cos(n.*(phi_o-phi_s));
        fs = sin(n.*(phi_o-phi_s));

        Ez_t(ii,jj) = sum(tau.*besselj(n,k(3)*rho_s).*besselh(n,2,k(3)*rho_o).*fc);
    end
end

矢量化解 -

%%// Term - 1
term1 = repmat(tau.*besselj(n,k(3)*rho_s),[N*N 1]);

%%// Term - 2
[n1,rho_g1] = meshgrid(n,rho_g);
term2_intm = besselh(n1,2,k(3)*rho_g1);
term2 = transpose(reshape(repmat(transpose(term2_intm),[N 1]),N,N*N));

%%// Term -3
angle1 = repmat(bsxfun(@times,bsxfun(@minus,phi_g,phi_s')',n),[N 1]);
fc = cos(angle1);

%%// Output
Ez_t = sum(term1.*term2.*fc,2);
Ez_t = transpose(reshape(Ez_t,N,N));

关于这种向量化或代码简化的注意事项 –

1. “fs”不会改变脚本 Ez_t 的输出，因此现在可以将其删除。
2. 输出似乎是“Ez_t”，这需要代码中的三个基本术语： tau.*besselj(n,k(3)*rho_s)、besselh(n,2,k(3)*rho_o) 和 fc >。这些分别作为项 1、2 和 3 分别计算用于矢量化。
3. 所有这三个术语的大小似乎都是 1xN。因此，我们的目标是在没有循环的情况下计算这三个项。现在，这两个循环各运行 N 次，因此我们的总循环数为 NxN。因此，与这些项在嵌套循环中时相比，我们必须将每个此类项中的数据设为 NxN 倍。
4. 这基本上是这里完成的矢量化的本质，因为这三个术语由“term1”、“term2”和“fc”本身表示。

Answer 3

为了给出一个自成一体的答案，我把原来的初始化复制过来

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

并生成一些缺失数据（n维的k(3)和rho_s和phi_s）

rho_s = rand(size(n));
phi_s = rand(size(n));
k(3) = rand(1);

那么您可以使用多维数组计算相同的 Ez_t：

[RHO_G, PHI_G, N] = meshgrid(rho_g, phi_g, n);
[~, ~, TAU] = meshgrid(rho_g, phi_g, tau);
[~, ~, RHO_S] = meshgrid(rho_g, phi_g, rho_s);
[~, ~, PHI_S] = meshgrid(rho_g, phi_g, phi_s);

FC = cos(N.*(PHI_G - PHI_S));
FS = sin(N.*(PHI_G - PHI_S)); % not used

EZ_T = sum(TAU.*besselj(N, k(3)*RHO_S).*besselh(N, 2, k(3)*RHO_G).*FC, 3).';

之后您可以检查两个矩阵是否相同

norm(Ez_t - EZ_T)