Implementing numerical integration in Python

Question

I have this Python code that implements a rectangular numerical integration. It evaluates the (K-1)-dimensional integral for arbitrary integer \$K \geq 1\$

$$\int_{u_K = 0}^{\gamma}\int_{u_{K-1} = 0}^{\gamma-u_K}\cdots\int_{u_2}^{\gamma-u_K-\cdots-u_3}F_U(\gamma-\sum_{k=2}^Ku_k)f_U(u_2)\cdots f_U(u_K)\,du_2\cdots du_K$$

where \$F_U\$ corresponds to the cdf function in the code, and \$f_U\$ to the pdf. I implemented it using recursion as follows:

    #************************** Import necessary libraries******************************************
import numpy as np
import matplotlib.pyplot as plt
import time
#******************************Set the constant scalars and vectors***************************
start_time = time.time()
KU = 3
eta_dB = 3
eta = 10**(eta_dB/10)
ExpanF = 50

tdB = np.arange(-5,11,4)
tVec = 10**(tdB/10)
thVec = (tVec/eta)*(ExpanF-eta*(KU-1))

N = 10000 # For simulation

du = 0.01

#******************************Define functions to be used***************************************
#Define the CDF of U
def CDF(u):
      return 1-1/(u+1)

#Define the PDF of U
def pdf(u):
     return 1/((1+u))**2

def FK(h, k):
    #print(f'h is {h}, and k is {k}')
    if k == 1:
        res = CDF(h)
    else:
        #n_iter = int(h/du)
        res =  0
        u = 0
        while u < h:
            res += FK(h-u, k-1)*pdf(u)*du
            u += du
    return res

#*******************Find the numerical and simulation values of the integral******************
ResultNum = []
ResultSim = []

if (ExpanF-eta*(KU-1)) > 0:
    for t in thVec:
    # Numerical
        ResultNum.append(1-FK(t, KU))
    # Simulation
        count = 0
        for n in range(0, N):
            if np.sum(np.random.exponential(1, (1, KU))/np.random.exponential(1, (1, KU))) <= t:
                count += 1
        ResultSim.append(1-count/N)

The parameter du is the resolution of the numerical integral, and the smaller, the more accurate the numerical approximation is. However, decreasing du comes at the expense of more computational time using the above code. For example I tried to run the above code on my machine Intel i5 @2.5 GH with 6 GB of RAM for du=0.0001, and it took 15 hours and didn't finish, at which case I had to abort it.

Is there anyway the above code can be optimized to run significantly faster? I tried C++, and I got a 30x speed factor using the same code and algorithm, but I was wondering if I can get a similar or better speed up factor in Python.

Answer 1

I have this Python code that implements a rectangular numerical integration.

In what sense "rectangular"? Or is that as opposed to contour integration?

It evaluates the (K-1)-dimensional integral for arbitrary integer \$K \geq 1\$

$$\int_{u_K = 0}^{\gamma}\int_{u_{K-1} = 0}^{\gamma-u_K}\cdots\int_{u_2}^{\gamma-u_K-\cdots-u_3}F_U(\gamma-\sum_{k=2}^Ku_k)f_U(u_2)\cdots f_U(u_K)\,du_2\cdots du_K$$

where \$F_U\$ corresponds to the cdf function in the code, and \$f_U\$ to the pdf.

It would probably help to rewrite it in fully symmetric form: $$\underset{\substack{u_i \ge 0 \\ u_1 + \cdots + u_K \le \gamma}}{\int \cdots \int} f_U(u_1)\cdots f_U(u_K) \textrm{d}u_1 \cdots \textrm{d}u_K$$

Then you can exploit the symmetry to divide the polytope of interest into equivalent chunks and only integrate each chunk once.

---

I implemented it using recursion as follows:

    #************************** Import necessary libraries******************************************
import numpy as np
import matplotlib.pyplot as plt
import time
#******************************Set the constant scalars and vectors***************************
start_time = time.time()

So far I understand.

KU = 3
eta_dB = 3
eta = 10**(eta_dB/10)
ExpanF = 50

tdB = np.arange(-5,11,4)
tVec = 10**(tdB/10)
thVec = (tVec/eta)*(ExpanF-eta*(KU-1))

What do any of those mean? "constant scalars and vectors" is not an adequate description.

N = 10000 # For simulation

du = 0.01

#******************************Define functions to be used***************************************
#Define the CDF of U
def CDF(u):
      return 1-1/(u+1)

#Define the PDF of U
def pdf(u):
     return 1/((1+u))**2

If \$f_U\$ is that simple, it might make more sense to tackle it with symbolic integration rather than quadrature.

def FK(h, k):
    #print(f'h is {h}, and k is {k}')
    if k == 1:
        res = CDF(h)
    else:
        #n_iter = int(h/du)
        res =  0
        u = 0
        while u < h:
            res += FK(h-u, k-1)*pdf(u)*du
            u += du
    return res

This could use some comments. It seems to be first order quadrature. If you use a higher order, maybe a larger step size will give you equally accurate results much faster.

#*******************Find the numerical and simulation values of the integral******************
ResultNum = []
ResultSim = []

if (ExpanF-eta*(KU-1)) > 0:
    for t in thVec:
    # Numerical
        ResultNum.append(1-FK(t, KU))
    # Simulation
        count = 0
        for n in range(0, N):
            if np.sum(np.random.exponential(1, (1, KU))/np.random.exponential(1, (1, KU))) <= t:
                count += 1
        ResultSim.append(1-count/N)

I'm missing a lot of context here.