Calculate a tricky sum in Matlab

问题:

With the following variables:

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>;

how do I compute the sum over m and n of S_mn*exp(1i*(m*theta + n*phi)), i.e.

I've thought of things like

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

but that requires theta and phi to have the same number of elements as m and n, and it gives me just one element in return - I want a matrix the the size of meshgrid(theta,phi), regardless of the sizes of theta and phi (i.e. I want to be able to evaluate the sum as a function of theta and phi).

How do I do this calculation in matlab?

回答1:

Since I don't know what S is...

S = randn(4,32);

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

Works fine for me.

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

If you now wish to do this for a large set of values phi and theta, a pair of loops now are the trivial solution. Or, you can do it all in one computation, although the arrays will get larger. Still not hard. WTP?

You do realize that both meshgrid and ndgrid take more than just two arguments? So it is time to learn how to use bsxfun, and then squeeze.

[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));

Or do this, which will be a bit faster. The previous computation took my machine .07 seconds. This last one took .05, so some savings by using bsxfun heavily.

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));

If you need to do the above 2000 times, so it should take 100 seconds to do. WTP? Get some coffee and relax.

回答2:

First save the size of each variable:

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

Use ngrid function like this:

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

This will give you an array with 4 dimensions (one for each of your variables: theta, phi, m, n). Now you can calculate this:

m.*theta + n.*phi

Now you need to make S_mn have 4 dimensions with sizes size_theta, size_phi, size_m, size_n like this:

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

Now you can calculate your sum like this:

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

Finally you can sum along the last 2 dimensions (m and n) to get an array with 2 dimensions with size size_theta by size_phi:

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

Note: I don't have access to Matlab right now, so I can't test if this actually works.

回答3:

There are several ways you could go about this.

One way is to create a function(-handle) that returns the sum as a function of theta and phi, and then use arrayfun to do the sums. Another is to fully vectorize the computation, though that will use more memory.

The arrayfun version:

[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);

The vectorized version:

[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);