HighwayJohn
HighwayJohn

Reputation: 921

Mean tensor product

I have another question which is related to my last problem( Python tensor product). There I found a mistake in my calculation. With np.tensordot I am calculating the following equation: enter image description here <..> should display the average. In python code it does look like this (ewp is a vector and re a tensor):

q1 = numpy.tensordot(re, ewp, axes=(1, 0))
q2 = numpy.tensordot(q1, ewp, axes=(1, 0))
serc = q2 ** 2

or

serc = numpy.einsum('im, m -> i', numpy.einsum('ilm, l -> im',
numpy.einsum('iklm, k -> ilm', numpy.einsum('ijklm, j -> iklm',
numpy.einsum('ijk, ilm -> ijklm', re, re), ewp), ewp), ewp), ewp)

Now in both python codes I neglect, that all possibilities are multiplied. But of course w_j and w_k are not independent for j=k. In the case, that only j and k are the same we get < w_j*w_j*w_l*w_m> = <w_j>*<w_l>*<w_m>. For j=k=l we get: < w_j*w_j*w_j*w_m> = <w_j>*<w_m>. For j=k=l=m: < w_j*w_j*w_j*w_j> = <w_j>. Only if all variable are different, independence is true and we get: < w_i*w_j*w_l*w_m> = <w_i>*<w_j>*<w_l>*<w_m>. Now this is what the code does for all possibilities. I hope this makes my problem understandable. Now my question is how can I represent this in my code?

Edit: An idea that I have is to first create a 4dim. tensor which represents <w_j w_k w_l w_m>:

wtensor = numpy.einsum('jkl, m -> jklm', numpy.einsum('jk, l -> jkl',
numpy.einsum('j, k -> jk', ewp, ewp), ewp), ewp)

Then I need to change the values which are not idependent. I assume they should be on a diagonal? But I really don't know so much about tensor calculus, so at this point I am struggling. After manipulating the w tensor I would get the result by performing:

serc = numpy.einsum('ijklm, jklm -> i', numpy.einsum('ijk, ilm ->
ijklm', re, re), wtensor)

Edit2: In another post I asked precisely how I can manipulate the 4dim so that it does fit here. Divakar had a really nice solution which can be seen here: Fill a multidimensional array efficiently that have many if else statements

from itertools import product

n_dims = 4 # Number of dims
# Create 2D array of all possible combinations of X's as rows
idx = np.sort(np.array(list(product(np.arange(gn),
repeat=n_dims))),axis=1)
# Get all X's indexed values from ewp array
vals = ewp[idx]
# Set the duplicates along each row as 1s. With the np.prod coming up
next,
#these 1s would not affect the result, which is the expected pattern
here.
vals[:,1:][idx[:,1:] == idx[:,:-1]] = 1
# Perform product along each row and reshape into multi-dim array
out = vals.prod(1).reshape([gn]*n_dims)

The array which I am getting here is wtensor which I now can use in the code above:

serc = numpy.einsum('ijklm, jklm -> i', numpy.einsum('ijk, ilm ->
ijklm', re, re), wtensor)

This gives me finally the result which I wanted and basically answers the question. Although there is one problem. The lenght of ewp which then also defines the size of the tensors shouldn't be larger then 6. Otherwise the code will use a lot of memory. My intention was to use it until a size of 8, so that is unfortunately now my next problem.

Upvotes: 1

Views: 260

Answers (1)

Divakar
Divakar

Reputation: 221664

Well you can do that efficiently with a combination of np.tensordot and np.einsum, like so -

serc = np.einsum('ilm,ilm->i',re,np.tensordot(re,wtensor,axes=[(1,2),(0,1)]))

Upvotes: 2

Related Questions