I have the following code which needs to calculate exponential of a 3D matrix. It works fine but I wonder if there is anyway to make it faster. (It currently needs more than 2 seconds. I am trying to make it run under half a second if possible.)
Result = zeros(200,50,300);
for i=1:30
delta = i*randn(200,50,300);
X = exp(1i*2*pi*delta);
Result = Result + X;
end
Any help would be appreciated. Thanks in advance.
First of all, I don't think that this sequence of computations can be made a lot faster. The following might be a tiny bit faster, but certainly not by the amount you're looking for:
dim = [200,50,300]; % given
N = prod(dim); % total number of samples
M = 30; % number of iterations in for-loop
phase = bsxfun(@times, randn(N,M), [1:M]); % scaled phases
Result = reshape(sum(exp(1i*2*pi*phase),2), dim); % sum of exp and reshape
EDIT 1 Start:
As @horchler pointed out in comments below, this method is actually slower than the original method on R2016b. He also suggested using a faster random generation scheme, which I tried and observed significant improvements. A similar amount of improvement can also be had by eliminating some of the temporary variables.
s = RandStream('dsfmt19937','Seed',1);
for i=1:30
Result = Result + exp(1i*2*pi*i*randn(s,200,50,300));
end
My timing results for various stages of optimization, on R2016b, are as follows:
Other random number generation schemes, such as shr3cong can also be tried in an effort to speed it up further.
EDIT 1: End
A different approach: Let me suggest a different approach that would generate the desired output in a different way, but such that it would have the same statistical properties as your output.
Result = sqrt(15) * (randn(200,50,300) + 1i*randn(200,50,300));
Let's try to justify this now. Firstly, we can argue that since the output is a sum of 30 random processes, by central limit theorem (CLT), the output would be Gaussian distributed. CLT only applies in loose sense here, since 30 is not infinite and the processes are not identically distributed. But as we will see shortly, it is still a very good approximation. Moreover, the real and imaginary terms are also independent, due to the sum over 30 independent complex numbers. I am not going to try to prove this here, but we will do some statistical checks.
Once we settle on independent Gaussian distributions, the analysis becomes much simpler. A Gaussian distribution can be defined by just two parameters: mean and the variance. Let's estimate them individually:
Mean: Since the phases are randomly distributed and cover a region much larger than 2*pi, the means for both the real and imaginary terms are 0.
Variance: The variance of sine/cosine of large random phase distribution is 0.5. So the variance for the sum of 30 sines/cosines would be 15. This is the reason for sqrt(15) term in the formula.
Statistical analysis: To verify that all the above assumptions and approximations are reasonable, lets perform some statistical analysis.
Firstly, let's check the distribution:
figure;
xGrid = (-15 : 0.1 : 15);
histogram(real(Result(:)), xGrid, 'Normalization','pdf', 'EdgeColor', 'None');
hold on;
plot(xGrid, normpdf(xGrid, 0, sqrt(15)), 'r', 'LineWidth', 2);
legend({'Simulated histogram', 'Gaussian pdf'});
title('Distribution of the real term');
The histogram of the imaginary term (not shown here) also looks identical. This test verifies the assumptions of Gaussian distribution with zero mean and 15 variance.
Finally, let's check the independence between the real and imaginary terms.
covar = cov(real(Result(:)), imag(Result(:)));
disp(covar);
% 14.9968 0.0036
% 0.0036 14.9936
Two things can be noticed: (1) As stated earlier, the variance of the real and imaginary terms each is about 15. (2) the covariance between the real and imaginary terms is much smaller compared to the individual variance. This supports the thesis of their being independent.
