Оптимизация численного интегрирования, включающего матрицы, с использованием Sympy ⇐ Python

[Anonymous]

Я интегрирую следующее выражение:


Где скрытые переменные выбираются из многомерного распределения следующим образом:


И случайные переменные w_n, заданные скрытой переменной, выбираются из многомерного распределения как следующее:


Где D — это ковариационная матрица, описывающая коэффициенты диффузии, которые мы пытаемся оценить, а v — это известная дисперсия, описываемая положительной полуопределенной матрицей
Код интегрируется по r_1,r_2,r_3, затем выводит выражение относительно D, затем находит D, хотя код не вызывает ошибка, но это занимает очень много времени, как в для выполнения следующего шага в
интеграции и оценке потребуется >10 минут, поэтому у меня вопрос:
Можно ли что-нибудь оптимизировать в коде? Не пропустил ли я какой-то важный шаг, который облегчит решение проблемы? если да, то что это за шаг? Правилен ли мой код, зная описанные выше многомерные распределения вероятностей?
Мой код выглядит следующим образом:
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)
sp.pprint(r_n)
q = w_n - 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,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)```


Подробнее здесь: https://stackoverflow.com/questions/793 ... sing-sympy