I am looking for the optimal way of solving a discrete arbitrarily-dimensional integral with numpy. Thing is, this integral operates numbers from matrices of different shape configurations, so indexing is a little tricky. What I currently have is the following:
'''
For whatever integer N >= 2: Hxy.ndim=N, Hx.ndim=N-1, Hy.ndim=1.
"nbins_xy" has the upper limit of all variables, lower limit is always 0.
'''
# All possible values my variables can take at any time
# (eg. for a triple integral[(1,2,3), ..., (2,4,1)]).
possible_values = [range(i) for i in nbins_xy]
# The actual integral, list comprehension + sum of its elements.
integral_naive = sum([Hxy[indices] * np.log(Hxy[indices] / Hx[indices[:-1]] / Hy[indices[-1:]])
for indices in itertools.product(*possible_values)])
I'm interested in boosting speed but also not having to store the list comprehension in memory would be good. Also I believe there must be a way of changing that for loop for a proper index array, but I haven't been able to figure it out. Solutions with torch tensors or similar are also welcome.
Could you please point me in the right direction? Thanks in advance!