# Integrate the product of four infinite series functions

I wrote the following code to integrate the product of four functions.

Two of the functions are hypergeometric functions and the other two are incomplete gamma functions. I have written each function as an infinite series, one per nested loop. The first loop (variable k) is for time.

y1 = 0;
for k = 1:L-1
for i = 0:N
total1 = i*log(x) + 2*(n-m-i+k)*log(sig);
for j = 0:N
total2 = j*log(y) + 2*(n-m-j+k)*log(sig);
for p = 0:100
total3 = 0;
for q = 0:100
temp = (n-m+i+k+p)*log(a);
temp = temp - 2*(2*n - 2*m + i + j + 2*k + p + q)*log(sig);
temp = temp - func1(k, j+1, q+1) + func2(i+1, j+1, p+1, q+1, k);
temp = temp + (i + j + 2*k + p + q)*log(sig^2/(1+a));
temp = temp + log(gammainc(b*(1 + x)/sig^2,i + j + p + q));
total3 = total3 + exp(temp);
end
total4 = func3(k, i+1) + (n - m + i + k + p)*log(a/sig^2);
total4 = total4 + func5(k, i+1, p+1);
total4 = total4 + (n - m + i + k + p + 2)*log(sig^2/a);
total4 = total4 + log(gammainc(a*b/sig^2, i + k + p)));
total4 = exp(total4) - total3;
y1 = y1 + exp(total1 + total2 - func5(k, i+1, p+1) + log(total4));
end
end
end
end


Thereafter, I removed the innermost loop the following way to improve speed:

y1 = 0;
for k = 1:L-1
for i = 0:N
total1 = i*log(x) + 2*(n-m-i+k)*log(sig);
for j = 0:N
total2 = j*log(y) + 2*(n-m-j+k)*log(sig);
for p = 0:100
q = 0:100;
temp = (n-m+i+k+p)*log(a);
temp = temp - 2*(2*n - 2*m + i + j + 2*k + p + q)*log(sig);
temp = temp - reshape(func1(k, j+1, q+1),1,101)
temp = temp + reshape(func2(i+1, j+1, p+1, q+1, k), 1, 101) ;
temp = temp + (i + j + 2*k + p + q)*log(sig^2/(1+a));
temp = temp + log(gammainc(b*(1 + x)/sig^2,i + j + p + q));
total3 = total3 + exp(temp);
total4 = func3(k, i+1) + (n - m + i + k + p)*log(a/sig^2);
total4 = total4 + reshape(func5(k, i+1, p+1), 1, 101);
total4 = total4 + (n - m + i + k + p + 2)*log(sig^2/a);
total4 = total4 + log(gammainc(a*b/sig^2,i + k + p)));
total4 = exp(total4) - total3;
y1 = y1 + exp(total1 + total2 - reshape(func5(k, i+1, p+1),1,101) + log(total4));
end
end
end
end


Can you suggest further changes to improve speed?

Basically this is the function I am trying to integrate:

\begin{eqnarray} F &=& \sum\limits_{k=1}^{L-1}\sum\limits_{i=1}^{\infty}\sum\limits_{j=1}^{\infty}\sum\limits_{p=1}^{\infty}\left(\frac{x^i\sigma^{2\left(n-m-i+k\right)}}{\left(n-m+1\right)_ii!}\frac{y^j\sigma^{2\left(n-m-j+k\right)}}{\left(n-m+1\right)_jj!}\frac{1}{s\left(s+1\right)\cdots\left(s+p\right)}\right)\nonumber\\ &&\left[\left(\Gamma\left(n-m+j+k\right)\left(\frac{a}{\sigma^2}\right)^{n-m+i+k+p}\int_0^{b*}e^{-\frac{ab}{\sigma^2}}b^{n-m+i+k+p+1}db\right)\right]\nonumber\\ &-&\left(\sum\limits_{q=0}^{\infty}\frac{a^{n-m+i+k+p}}{\sigma^{2\left(2n-2m+i+j+2k+p+q\right)}}\frac{1}{s(s+1)\cdots(s+q)}\int_0^{b*}e^{-\frac{ab+b}{\sigma^2}b^{2n-2m+i+j+p+q+2k+2}}db\right) \end{eqnarray}

• I edited the question. The innermost loop of the code calculates the innermost sum i.e. $\sum\limits_{q=0}^{100}$. I took the number of terms to be 100 for each sum. Commented Jan 30, 2015 at 0:01
• I know it is a bit difficult to grasp the expression. It is product of four series. And it is taking a lot of time for MATLAB to calculate it. In fact, F is an element of a 4 X 4 matrix. So, I have to find F for 16 times just to get one evaluation point. And I have to get at least 100 evaluation points. :( Commented Jan 30, 2015 at 0:06

Caveat: I am not a Matlab expert, nor am I competent at an advanced undergraduate level in mathematics.

## Efficiency

Removing data structure impedance and redundant or unnecessary calculations is likely to improve the running time of the code under review.

1. func5 [whatever the func that is] is called twice with exactly the same arguments.

2. Each time it is wrapped in reshape which adds to the cost.

3. Indeed func1, func2, and func5 are defined each returning a data structure that requires reshaping. This can be avoided by either working with their data structure directly in the code under review or redefining the functions to return the data structure used by the code.

4. (n - m + i + k + p)*log(a/sig^2) is calculated twice. Sure in one location it needs to be doubled.

5. Each time through the inner loop (n - m + i + k + p) increases by one. There's no need to recalculate it in full. It can be passed as a parameter from the enclosing function.

6. a [defined elsewhere?] is a constant within the code under review. Therefore log(a/sig^2) is also a constant. Yet is being calculated twice in the inner loop.

7. (log sig) [sig defined elsewhere?] is calculated for each iteration at the each level of the loop. It's a constant.

8. There are probably other places where redundant computations can be avoided.

## Code Clarity

The code is hard to read.

The more that the names of items in the code reflect the way we would label the pieces of the problem domain [mathematics here], the easier it becomes to reason about the way the code implements the problem.

sigma_k better expresses the mathematics than total1. And though I know temp = temp + ... does something, but it's not obvious what. function1_result = start_value + ... is better. partialgamma1_result = start_value + ... is better.

## Final thoughts

Giving things meaningful names helps a lot. logsig = log(sig) seems like more work because it is more typing [by 8 more keystrokes]. But it's a false economy because debugging and refactoring code for speed is orders of magnitude more difficult than banging on the keyboard.