The most efficient way to do multiple summation in Matlab?
5 views (last 30 days)
Show older comments
Hi
Recetly I've got a problem about multiple summation calculation in Matlab. As the expression below shows, I would like to numerically evaluate y w.r.t. x and k.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/732879/image.png)
Previously I use the function 'symsum' in a nested way to get an analytical expression of y without the summation signs, then evaluate the expression numerically, but the evaluation would cost much time (up to a couple of days).
Therefore I wonder if there are any other more efficent ways to do so, like vectorise the multiple summation (have got no clue about how to vectorise it)or using orther functions like bsxfun (it is said that this function is effcient in terms of loop?) or cumsum?
Thanks in advance for your help!
4 Comments
Accepted Answer
Jan
on 8 Sep 2021
Edited: Jan
on 8 Sep 2021
After the code runs in half a minute now instead of a couple of days, it is time to examine the numerical stability:
for i3 = 0:i1 + i2
y = y + x3 * xpow(c + i3);
end
Due to the accumulation, y and be larger than the 2nd part of the sum. This reduces the accuracy. It is better to accumulate the inner loop separately. In addition all terms of the sum are multiplied by x^3, so this can be done after the loop once:
s = 0;
for i3 = 0:i1 + i2
s = s + xpow(c + i3);
end
y = y + s;
...
y = y * x^3;
Of course this changes the output slightly:
% For x=3e-3, k=300
% 0.0365429439932915 original
% 0.0365429439930864 separate accumulation - should be more accurate
Now it is time to vectorize the inner loop:
function y = YourSum_vec(x, k)
x_1 = 1 - x;
xpow = zeros(1, 5 * k);
for ii = 1:5 * k
xpow(ii) = x_1 ^ ii;
end
y = 0;
for i1 = 0:k
for i2 = 0:i1
c = i1 + i2 + k;
y = y + sum(xpow(c:c + i1 + i2));
end
end
y = y * x^3;
end
This reduces the runtime from 32.1 seconds to 20.8 seconds.
Compared to several days of the symbolic solution, this is both much faster. The vectorization was not the main part.
More Answers (1)
Jan
on 8 Sep 2021
Edited: Jan
on 8 Sep 2021
% Timings on: Matlab R2018b, Win10, i7:
tic % A small k for bearable run times at first:
YourSum_orig(3e-3, 300) % 0.0365429439932915
toc % 4.236337 seconds
tic
YourSum_inline(3e-3, 300) % 0.0365429439932915
toc % 1.715529 seconds
tic
YourSum_pre(3e-3, 300) % 0.0365429439932915
toc % 0.066605 seconds, 70 times faster!
tic % Now with the full data:
YourSum_pre(3e-3, 3000) % 4.56878584437021e-05
toc % 35.517909 seconds
And the different implementations:
function y = YourSum_orig(x, k) % Original loop
f = @(i1, i2, i3, x, k) (x.^3).*(1-x).^(i1+i2+i3+k);
y = 0;
for i1 = 0:k
for i2 = 0:i1
for i3 = 0:i1 + i2
y = y + f(i1, i2, i3, x, k);
end
end
end
end
Now insert the function in the code and move repated calculations out out the loop:
function y = YourSum_inline(x, k)
y = 0;
x3 = x^3;
for i1 = 0:k
for i2 = 0:i1
for i3 = 0:i1 + i2
y = y + x3 * (1-x) ^ (i1 + i2 + i3 + k);
end
end
end
end
The power operation is very expensive, but repated frequently for the same input. So pre-calclutate it once:
function y = YourSum_pre(x, k)
y = 0;
x3 = x^3;
x_1 = 1 - x;
xpow = zeros(1, 5 * k);
for ii = 1:5 * k
xpow(ii) = x_1 ^ ii;
end
for i1 = 0:k
for i2 = 0:i1
c = i1 + i2 + k;
for i3 = 0:i1 + i2
y = y + x3 * xpow(c + i3);
end
end
end
end
2 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!