The most efficient way to do multiple summation in Matlab?

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.

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!

[Jan]

It depends on what f is. There is no way to implement this general formula in efficient code.

What are the values x and k?

[Shuangfeng Jiang]

Thanks for your reply!

the f is:

f=@(i1, i2, i3, x, k) (x.^3).*(1-x).^(i1+i2+i3+k);

where k is around 3000 and x is around 3e-3.

[Jan]

Is x a scalar or vector?

[Shuangfeng Jiang]

x is a scalar

Answer 1

[Jan]

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.

[Shuangfeng Jiang]

Thanks Jan, your answers are impressive!

btw, it seems that your codes also work when x and k are vectors rather than scalars.

Answer 2

[Jan]

% 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

[Shuangfeng Jiang]

Thanks for your reply!

Yeah I believe this loop should be feasible, but can we vectorise this summation calculation somehow to shorten the executing time?

Cheers

[Jan]

I've inserted the function in the loop and ran some speed comparisons. Avoiding repeated calculations of the same power operation is much more important than a vectorization.

I post a 2nd answer to reduce the clutter.