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.

• 200_success: Please fix your indentation. The easiest way to post code is to paste it into the question editor, highlight it, and press Ctrl-K to mark it as a code block.
• Gareth Rees: Did you try scipy.integrate.nquad?
• BlackMath: @GarethRees Not really. Can it be used for an arbitrary number of nested integrals?
• Gareth Rees: See the documentation: "Function signature should be func(x0, x1, ..., xn, t0, t1, ..., tm). Integration is carried out in order. That is, integration over x0 is the innermost integral, and xn is the outermost."
• BlackMath: @GarethRees Doesn't sound like what I am looking for. I don't want to specify the functions and limits each time I change K. I need a dynamic program where all I change is K. That is why I am trying to program the numerical integration from scratch. Otherwise, I would have used the built-in numerical integrations in MATLAB or any other platform. But these functions are not flexible for an arbitrary number of nested integrals.

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.