How to compute double sum in Matlab where upper bound of second sum is lower bound of first sum?

Question

I am having big difficulties implementing a function in Matlab that calls other functions that I have programmed in different .m files. The part where I am stuck is the part where a sum over different values inputted in another function is taken, where there is also a sum inside the other function. The lower bound of the first sum is the upper bound of the second sum.

The functions are:

I have the Hh(n,x) function working properly for n inputted as a vector and x inputted as a scalar. Because of the vectorized input of n, the sum over Hi inside the In function can be quickly calculated by calling sum(Hh(0:n,x)).

I want to do the same for the In function, but because n inside In now ranges from 0 to k-1 and in the outer function k ranges from 1 to n, I do not know how to evaluate this double sum, where the inner sum has the lower bound of the outer sum as an upper bound. I want to evaluate this double sum as efficiently as possible, since later I want to do many simulations with these formulas. Now I am evaluating the function In n times, storing each value in a vector and then taking the sum, which is very computationally intense...

My Matlab code for the In function is:

function in = In(n,c,alphaa,betaa, delta)
ie = 0:n;
in = -(exp(alphaa*c)/alphaa)...
     .*sum((betaa/alphaa).^(n-ie).*Hh(ie,betaa*c-delta))...
     -(betaa/alphaa).^(n+1)
end

The Matlab code for the outer function, lets call it function f for now is:

function f = F(n,a,mu,sigma,eta1,T)
for k = 1: n
    vector(k) =  In(k-1,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T)));
end
f = sum(vector);
end

How can I make the input of n inside In vectorized, so that I do not have to store all inputted n values separately and then calculate the sum, but calculate the sum directly for inputted vector n.

Any help is appreciated since I am seriously stuck at the moment! Thank you in advance!

Answer 1

So to answer specifically your question (input of n vectorized), I would do the following:

% 1. create a vector ivect going from 0 to k-1 for every k, and a vector kvect corresponding to the iterate ivect:

ivect=[];
kvect=[];
for k=1:n
    ivect = [ivect 0:k-1];
    kvect = [kvect (k-1)*ones(1,k)];
end

% 2. create the function that is inside the two sums, depending on your iterates i and k:

    function in = In2(k,ie,c,alphaa,betaa, delta)
in = -(exp(alphaa*c)/alphaa)*((betaa/alphaa).^(k-ie).*func(ie,betaa*c-delta))-(betaa/alphaa).^(k+1)./(k+1);
end

% 3. then sum:

f=sum(In2(kvect,ivect,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T))))

At first I thought this method would reduce the computational time in the computation of the sum (of course point 1. is costly, but can be done one time and be used for several simulations without the need for reconstructing these vectors), but even comparing just the lines:

tic
for k = 1: n
    vector(k) =  In(k-1,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T)));
end
f = sum(vector)
toc

with

tic
f=sum(In2(kvect,ivect,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T))))
toc

I obtain 1.029899 seconds for the first method, and 1.684432 seconds for the second one, with n = 5000. But fortunately, I found why! This is due to the fact that this method actually engenders more artihmetic manipulations, as in the function In2 I oblige my code to compute -(betaa/alphaa).^(k+1)./(k+1) for every k and i, while this term just depends on k. So now, lets define functions In3 and In4 such that:

function in = In3(k,ie,c,alphaa,betaa, delta)
in = -(exp(alphaa*c)/alphaa)*((betaa/alphaa).^(k-ie).*Hh(ie,betaa*c-delta));    end

function in = In4(k,alphaa,betaa)
    in=-(betaa/alphaa).^(k+1);
end

and call

tic
f=sum(In3(kvect,ivect,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T)))) +sum(In4([1:n],-eta1,-1/(sigma*sqrt(T))))
toc

And voila, I obtain the same result in 0.929530 seconds for N=5000. So to sum up, you gain a bit of time with this last method, but to be completely fair I would say that you have two drawbacks:

1. The time computation of ivect and kvect at the beginning - but I insist these vectors can be constructed once for all, stored, and used for several simulations-

2. The memory cost: ivect and kvect are of size n(n+1)/2, which is a lot for large n

Hope it helped you, and anyway thanks for the interesting question which will probably serve me in my own research!