¿Cómo reducir el tiempo de ejecución en mi código de python?

Actualmente estoy trabajando en un proyecto que requiere que ejecute una base de código completo de Python. Para fines de investigación, necesito ejecutar el código lo más rápido posible. Sin embargo, soy bastante nuevo en programación y no tengo idea de cómo reducir el tiempo de ejecución. Así que espero que alguien pueda ayudarme en eso. Cualquier consejo sería apreciado. Aquí está parte de mi base de código, que utiliza muchos bucles “for” anidados, por lo que podría aumentar significativamente el tiempo de ejecución.

 def a_j(r, a, A): # el factor de Claussius-Mossotti, determinado por una matriz simétrica (3 x 3) tal que (A)_i^T = A_i
        alph = np.array([[0,0,0],[0,0,0],[0,0,0]],complex)
        for i in range(3):
            for j in range(3):
                alph[i,j] = (r * a * A[i,j])
        return alph

def W_ext(x, k, rho, alpha, A): # término de interacción partícula-partícula
    n = x.shape[0] # el número de vectores x
    result = np.zeros([3*n,3*n],complex)
    u = np.zeros((n, 3)) # u = x - x'

    for i in range(n):
        for j in range(n):
            if i != j:
                u[i] = x[i] - x[j]
                block_result = a_j(rho[i], alpha, A) * G((u[i]), k) * a_j(rho[j], alpha, A) 
                for m in range(3):
                    for l in range(3):
                        result[3*i + m, 3*j + l] = block_result[m,l]
    return result.imag

def A_ext(rho, a, A): # término de una sola partícula
    n = rho.shape[0]
    result = np.zeros([3*n,3*n],complex)

    for i in range(n):
        for j in range(n):
            if i == j:
                block_result = a_j(rho[i], a, A).imag
                for m in range(3):
                    for l in range(3):
                        result[3*i + m, 3*j + l] = block_result[m,l]
    return result # matriz (3 x 3)

def P_ext(e, A, W, omega):
    eT = np.matrix.getH(e)
    mm1 = np.matmul(A, e)
    mm2 = np.matmul(W, e)
    extinction = (np.dot(eT, mm1) + np.dot(eT, mm2)) * (omega/2.0)
    return extinction

#ABSORPTION

def W_abs(x, k, rho, alpha, A, chi): # término de interacción partícula-partícula
    n = x.shape[0] 
    result = np.zeros([3*n,3*n],complex)
    u = np.zeros((n, 3)) 

    for i in range(n):
        for j in range(n):
            if i != j:
                u[i] = x[i] - x[j]
                block_result = np.matrix.getH(a_j(rho[i], alpha, A)) * (1.0 / np.conjugate(chi)).imag * a_j(rho[i], alpha, A) * G((u[i]), k) * a_j(rho[j], alpha, A) 
                for m in range(3):
                    for l in range(3):
                        result[3*i + m, 3*j + l] = block_result[m,l]        

    return 2.0 * result.real  # matriz (3 x 3)

def A_abs(rho, a, A, chi): # término de una sola partícula
    n = rho.shape[0]
    result = np.zeros([3*n,3*n],complex)

    for i in range(n):
        for j in range(n):
            if i == j:
                block_result = np.matrix.getH(a_j(rho[i], a, A)) * (1.0 / np.conjugate(chi)).imag * a_j(rho[i], a, A)
                for m in range(3):
                    for l in range(3):
                        result[3*i + m, 3*j + l] = block_result[m,l]
    return result # matriz (3 x 3)