How to change double sum to be computed in manageable time?

Hello.

I am trying to find efficient solution for computation of double sum. To be more specific I have algorithm like following.

i=1:X;
for j=1:Y
A(j) = sum(f1(i).^j);  %f1() is some function
end
Result=sum(A);

The above is the idea that is working. However since my Y=100K or so and X=10K the computation is taking so long. Moreover, I would need to do the computation for diferent Y=1:100K and observe the results. That's way too impossible with this approach.

Can anyone give me suggestions how to do this task in real time?

Answer 1

[Star Strider]

Unless you have a specific reason for storing the A vector, one way to speed it up might be:

i=1:X;
A = 0;
for j=1:Y
    A = sum(f1(i).^j) + A;  %f1() is some function
end
Result = A;

The reason is that for each new value for A(j), MATLAB has to allocate and fill a new memory location. This takes time. You might be able to get around this by preallocating the A vector if you need to save it for use later in your program by putting:

A = zeros(size(Y));

in place of the A = 0; line, and using the code you posted with A(j) defined in every step of the loop.

( NOTE: it’s best not to use i and j as loop indices because MATLAB uses them for its imaginary operators. Using them as loop indices could confuse things. )

[Peter]

Thank you for your suggestion. I indeed don't need A to have stored as a vector. So I have changed it at least to save some operations. However using tic/toc I have noticed only very little improvement.

I was wondering whether there is an option to avoid using loop and use some matlab features instead.

[Star Strider]

I just now noticed that f1 doesn’t change within the loop, so (unless I’m missing something) you can define it outside the loop so you don’t have to keep evaluating it with every iteration of the loop. Your code then becomes:

X = 10;
Y =  7;
i = 1:X;
f1 = @(i) sin(2*pi*i./X) + 1;       % Create ‘f1’
f1vct = f1(i);                      % Create f1 vector (doesn’t change in the loop)
A = 0;
for j = 1:Y
    A = sum(f1vct .^ j) + A;
end

The compiler might pick up on this itself and optimise it, but I’m guessing that it doesn’t. I don’t see how you can avoid the j loop.

Does that go faster?

------

AFTERTHOUGHT — I don’t know what f1 is, but if it goes to zero as j increases, you could test for an incremental value of A to see how it compares to the running sum. If it’s < 1E-15 of the running sum, you can probably exit the j loop without seriously affecting accuracy. The code for that might look like this:

X = 1000;
Y =  100000;
i = 1:X;
f1 = @(x) i/1E+5;            % Create ‘f1’
f1vct = f1(i);               % Create f1 vector (doesn’t change in the loop)
Result = 0;
for j = 1:Y
    A = sum(f1vct .^ j);
    if j > 1 && A/Result < 1E-15
        break
    elseif j >= 1
        Result = Result + A;
    end
end

[Peter]

The truth is that function was just simplification

f1 = a.*(1-(1-a).^j)

a and b are functions. But using the idea I changed it to

f1 = a-a*(1-a).^j
Computation of a and (1-a) prior to loop reduced time rather significantly, approx. by 25-35 %.

Thank you.

[Star Strider]

My pleasure!

Answer 2

[Roger Stafford]

Your addition is in the wrong order, Peter. If you do the summation with respect to j first, you can take advantage of a formula for the sum of a geometric series and avoid actually doing the summation. The formula is:

 a^1 + a^2 + a^3 + a^4 + ... + a^n =
 (a-a^(n+1))/(1-a)

which certainly saves a lot of additions. Hence your code could be:

F = f1(1:X);
Result = sum((F-F.^(Y+1))./(1-F));

You will notice that this form also avoids evaluating f1(i) for the same i repeatedly instead of just once.

There is one problem however. You state that Y can be as large as 100000. To avoid overflow to infinity or underflow to zero in matlab's 'double' numbers, the values produced by f1 need to be extremely close to 1, and this will lead to substantial loss of accuracy in the above formula because you will be subtracting two quantities which are very close to each other in the denominator, 1-F.

I can conceive of a compromise between adding 100000 numbers on the one hand and suffering large errors on the other, by using the above formula for, say, every 100 values and adding the resulting 1000 values, or some such strategy.

[Peter]

Roger, I'm so grateful that you remind me of geometric series. I did not see it there however I have already revised information about it and did some changes to formula. The great thing was that a was:

a = (1-p)^j

I start counting j from 0

j = 0:n-1

so using geometric series I ended up

sum((1-p)^j) = (1-(1-p)^n)/(1-(1-p)) = (1-(1-p)^n)/p)

So now it was really easy to write down in code.

Thank you very much.