Optimizing a numerical integration that involves matrices using sympy

The expression I'm integrating is

Where the latent variables are sampled from the a multivariate distribution as following:

And the random variables w_n given a latent variable is sampled from a multivariate distribution as following:

Where D is a covariance matrix describing the diffusion coefficients that we're attempting to estimate, and the v is the known variance described by a positive semi-definite matrix

The code integrates over r_1,r_2,r_3 and then derive the expression with respect to D then solve for D, although the code doesn't raise an error but it takes an extremely long time as in it would take >10 minutes to produce the next step in the integration and evaluation so my question is;

Can anything be optimized in the code? Am I missing some crucial step that makes the problem easier to solve? if so, what is that step? Is my code even correct, knowing the multivariate probability distributions described above?

My code is as following:

import sympy as sp
from sympy.vector import CoordSys3D
from sympy.matrices import Matrix

Coord = CoordSys3D('Coord')

N = 3

r_1 = Matrix([1,1,1])

#sp.pprint(r_n)

V = sp.MatrixSymbol("V",N,N)
V = Matrix([[1,0,0],[1,1,0],[1,0,1]])
V = V.T * V

#sp.pprint(V)

D_1 = sp.symbols("D_1")
D_2 = sp.symbols("D_2")
D_3 = sp.symbols("D_3")
D = Matrix([[D_1,D_1*D_2,D_1*D_3],[D_2*D_3,D_2,D_3*D_3],[D_3*D_1,D_3*D_2,D_3]])
#sp.pprint(D)

w_1 = Matrix([1,1,1])
w_2 = Matrix([2,2,2])
w_3 = Matrix([3,3,3]) 
w_n = [w_1,w_2,w_1]
r_n_list = [sp.symbols(f'r_{i+1}') for i in range(N)]
r_n = Matrix([r_n_list[0],r_n_list[1],r_n_list[2]])

t=0

def expo_power(i=0,v=V,t=t,D=D):
    sp.pprint(w_n[i])
    sp.pprint(r_n)
    q = w_n[i] - r_n
    sp.pprint(type(sympy.Transpose(q)))
    sp.pprint(type(sympy.Inverse(V)))
    sp.pprint(type(q))
    f = (sympy.Transpose(q) * sympy.Inverse(V) * (q)/2)
    """sp.pprint(type(sympy.Transpose(q)))
    sp.pprint(type(sympy.Inverse(V)))
    sp.pprint(type(q))"""
    return f
sp.pprint(expo_power())

def expo_power_sum(lower_bound = 2, upper_bound=N,t=t,D=D,v=V):
    f = expo_power(lower_bound-1,t=t,D=D,v=V)
    for i in range(lower_bound, N):
        f += expo_power(i,t=t,D=D,v=V)  
    return f

def P_w_n_P_r_n(i=0,v=V,D=D,t=t):  
    return (1/((sp.Determinant(V) ** (1/2) * (2 * sp.pi) ** (3/2)) * \
                sp.Determinant(2 * D * t) ** (1/2) * (2 * sp.pi) ** (3/2)) ** (N-1)) \
        * sp.exp(-expo_power_sum(lower_bound=2,upper_bound=N,t=t,D=D,v=V)) 

def P_w_i_r_i(i = 0, v=V):
    return (1/(sp.Determinant(V) ** (1/2) * (2 * sp.pi) ** (3/2))) \
        * sp.exp(-expo_power(i=i,v=V))

def P_r_n_r_n_1(i =0, m = r_n[0], v =V):
    return (1/sp.Determinant(2 * D * t) ** (1/2) * (2 * sp.pi) ** (3/2)) \
        * sp.exp(-expo_power(i=i,v=V))

def integrand(v = V,t=t,lower = 2):
    f = P_w_n_P_r_n(i=lower,v=V,D=D,t=t) * P_w_i_r_i(v=v) * P_r_n_r_n_1(i=0,v=v)
    #f = f * P_w_i_r_i(v=v)
    #f = f * P_r_n_r_n_1(i=r_n[0],v=v)
    sp.pprint(P_w_n_P_r_n(i=lower,v=V,D=D,t=t))
    sp.pprint(P_w_i_r_i(v=v))
    sp.pprint(P_r_n_r_n_1(i=0,v=v))
    return f

def integrate(v=V,t=t, lower_bound = -sp.oo, upper_bound= sp.oo):
    M = integrand(v=v,t=t)

    for i in range(N):  
        inte = sp.Integral(M,(r_n[i],lower_bound,upper_bound))
        #sp.pprint(inte)
        inte = inte.doit()
        inte = inte.evalf()
        M = inte
        sp.pprint(inte)
    
    
    k = inte
    sp.pprint(k)
    sp.pprint(type(M))
    
    d = sp.diff(k,D)
    
    sol = sp.solve(d,D)
    sp.pprint(sol)

integrate(t=1,lower_bound=0,upper_bound=1)```