1

Imagine that you have a tensor $T_{ijkl}$ (0<=i,j,k,l< N) that has the following property

T_{ijkl} is nonzero if i+j=k+l

which acts as a "selection rule".

You want to contract the tensor T with another tensor u to make a third tensor v as follows

$$v_{ij} = sum_{kl=0}^N T_{ijkl} u_{kl}$$

is there a way to do this construction such that:

a. Uses the tensor contraction machinery of numpy and/or tensorflow? Things like einsum and the like.

b. Uses the explicit symmetry of the tensor, i.e the above operation should require $O(N^3)$ operations, not $O(N^4)$.

One could use einsum directly getting $O(N^4)$ operations:

v = np.einsum('ijkl,ij->kl',T,u)

Or one could do a custom triple loop as

v = np.zeros_like(u)
for i in range(N):
    for j in range(N):
        for k in range(max(1+i+j-N,0),min(i+j,N-1)+1):
            v[i,j]+=T[i,j,k,i+j-k]*u[k,i+j-k]

But neither of this satisfy both a. and b. above

Nicolás Quesada
  • 139
  • 2
  • `v = np.einsum('ijkl,kl->ij', T, u) – hpaulj Jan 03 '20 at 03:58
  • This will imply doing N^4 operations which is suboptimal given the structure of the tensor T. One can do it in ~ 0.66 N^3 operations using 3 loops as shown above. – Nicolás Quesada Jan 03 '20 at 04:02
  • @NicolásQuesada Note that `O(N^4)$ operations` with einsum which is implemented in C is an entirely different ball game than `O(N^3)` with your loop based which is at Python level. If you are not stuck up on computational-complexity, actually timing could clarify/verify those. Also, it seems `np.tensordot` could be used and might be better than `einsum`. – Divakar Jan 03 '20 at 05:31
  • `(T.reshape(i*j, k*l) @ u.reshape(k*l)).reshape(i,j)` might be faster than the `einsum` (that in effect is what `tensordot` does). – hpaulj Jan 03 '20 at 05:40
  • @Divakar Agreed. But If I just take the python sum, turn into a function add a numba decorator suddenly this sum way of doing contractions is faster. – Nicolás Quesada Jan 03 '20 at 16:50
  • BTW: You don't need T at all if you implement it in Numba or Cython. For larger N eg.100 the construction of T is by far the most costly part. Also if you already have a tensor T, calculating just `v[i,j]+=u[k,i+j-k]` in the inner loop will be faster than your example. – max9111 Jan 08 '20 at 12:42
  • @max9111 thanks for your comment! I am confused what do you mean by not needing T? – Nicolás Quesada Jan 09 '20 at 13:06
  • In your loop implementation the accessed elements of T are always one/True (at least I assume this because you wrote that T is a selection rule). Therefore you don't need T at all to calculate the tensor `v` depending on the tensor `u`. `v[i,j]+=u[k,i+j-k]` is enough in the inner loop. – max9111 Jan 09 '20 at 13:12
  • Hi --- I guess is missused the word selection rule. I just meant that the elements that satisfy the rule are nonzero, not that they are all equal. – Nicolás Quesada Jan 10 '20 at 14:08

0 Answers0