在 Matlab 中计算一个棘手的总和

Question

具有以下变量：

m = 1:4; n = 1:32;
phi = linspace(0, 2*pi, 100);
theta = linspace(-pi, pi, 50);

S_mn = <a 4x32 coefficient matrix, corresponding to m and n>;

我如何计算 S_mn*exp(1i*(m*theta + n*phi)) 的 m 和 n 的总和，即

我想过像这样的事情

[m, n] = meshgrid(m,n);
[theta, phi] = meshgrid(theta,phi);
r_mn = S_mn.*exp(1i*(m.*theta + n.*phi));
thesum = sum(r_mn(:));

但这需要theta 和phi 具有与m 和n 相同数量的元素，它只给了我一个元素作为回报——我想要一个大小为meshgrid(theta,phi) 的矩阵，无论theta 和phi 的大小如何（即，我希望能够将总和评估为theta 和phi 的函数）。

如何在 matlab 中进行计算？

Answer 1

因为我不知道 S 是什么...

S = randn(4,32);

[m,n] = ndgrid(1:4,1:32);
fun = @(theta,phi) sum(sum(S.*exp(sqrt(-1)*(m*theta + n*phi))));

对我来说很好。

fun(pi,3*pi/2)
ans =
          -15.8643373238676 -      1.45785698818839i

如果您现在希望对大量 phi 和 theta 值执行此操作，那么现在一对循环是简单的解决方案。或者，您可以在一次计算中完成所有操作，尽管数组会变得更大。还是不难。 WTP？

您确实意识到 meshgrid 和 ndgrid 都需要两个以上的参数？所以是时候学习如何使用 bsxfun，然后再挤压。

[m,n,theta,phi] = ndgrid(1:4,1:32,linspace(-pi, pi, 50),linspace(0, 2*pi, 100));
res = bsxfun(@times,S,exp(sqrt(-1)*(m.*theta + n.*phi)));
res = squeeze(sum(sum(res,1),2));

或者这样做，会快一点。之前的计算花费了我的机器 0.07 秒。最后一个花费了 0.05，因此通过大量使用 bsxfun 可以节省一些费用。

m = (1:4)';
n = 1:32;
[theta,phi] = ndgrid(linspace(-pi, pi, 50),linspace(0, 2*pi, 100));
theta = reshape(theta,[1,1,size(theta)]);
phi = reshape(phi,size(theta));
res = bsxfun(@plus,bsxfun(@times,m,theta*sqrt(-1)),bsxfun(@times,n,phi*sqrt(-1)));
res = bsxfun(@times,S,exp(res));
res = squeeze(sum(sum(res,1),2));

如果你需要做2000次以上，那么应该需要100秒。支付意愿？喝杯咖啡放松一下。

Answer 2

先保存每个变量的大小：

size_m = size(m);
size_n = size(n);
size_theta = size(theta);
size_phi = size(phi);

像这样使用ngrid函数：

[theta, phi, m, n] = ngrid(theta, phi, m, n)

这将为您提供一个 4 维数组（每个变量一个：theta、phi、m、n）。现在你可以计算了：

m.*theta + n.*phi

现在你需要让 S_mn 有 4 个维度，大小为 size_theta、size_phi、size_m、size_n，如下所示：

S_tpmn = repmat(S_mn, [size_theta size_phi size_m size_n]);

现在您可以这样计算总和：

aux_sum = S_tpmn.*exp(1i*(m.*theta + n.*phi));

最后，您可以将最后 2 个维度（m 和 n）相加，得到一个大小为 size_theta by size_phi 的 2 个维度的数组：

final_sum = sum(sum(aux_sum, 4), 3);

注意：我现在无法访问 Matlab，所以我无法测试这是否真的有效。

Answer 3

有几种方法可以解决这个问题。

一种方法是创建一个函数（-handle），将总和作为 theta 和 phi 的函数返回，然后使用 arrayfun 进行总和。另一个是完全向量化计算，虽然这会使用更多的内存。

arrayfun 版本：

[m, n] = meshgrid(m,n);

sumHandle = @(theta,phi)sum(reshape(...
    S_mn.*exp(1i(m*theta + n*phi)),...
    [],1)) 

[theta, phi] = meshgrid(theta,phi);

sumAsFunOfThetaPhi = arrayfun(sumHandle,theta,phi);

矢量化版本：

[m, n] = meshgrid(m,n);
m = permute(m(:),[2 4 1 3]); %# vector along dim 3
n = permute(n(:),[2 3 4 1]); %# vector along dim 4

S_mn = repmat( permute(S_mn,[3 4 1 2]), length(theta),length(phi));

theta = theta(:); %# vector along dim 1 (phi is along dim 2 b/c of linspace)

fullSum = S_mn.* exp( 1i*(...
    bsxfun(@plus,...
       bsxfun(@times, m, theta),...
       bsxfun(@times, n, phi),...
    )));

sumAsFunOfThetaPhi = sum(sum( fullSum, 3),4);