Explination of numpy's einsum

Question

I am currently doing some studies on computing a 4th order tensor in numpy with the einsum function.

The tensor I am computing is written in Einstein notation and the function einsun does the work perfectly! But I would like to know what it is doing in the following case:

import numpy as np

a=np.array([[2,0,3],[0,1,0],[0, 0, 4]])
b= np.eye(3)

r1=np.einsum("ij,kl->ijkl", a, b)
r2=np.einsum("ik,jl->ijkl", a, b)

in r1 I am basically doing the standard tensor product (equivalent to np.tensordot(a,b,axes=0)).

What about in r2?

I know I can get the value by doing a[:,None,:,None]*b[None,:,None,:] but I do not know what the indexing is doing. Does this operation have a name?

Sorry if this is too basic!

I think of both r1 and r2 as an outer product. The 2nd just reorders the axes of the 4d result. Matrix product can do the same sort of multiplication, but it sums one or more axes (the shared index). In your case there's no summation of products. r1 can also be written with broadcasting - just reorder the None. — hpaulj
np.tensordot reshapes and transposes the inputs till it can do a standard 2d dot product, and then transforms the result back. My memory is that the single number axes parameter is translated into an expanded tuple of axes. I find that the einsum indexing is clearer, and more powerful. — hpaulj
As to your question on what the indexing is doing, it's not really clear what you mean but it sounds like you are asking what None is doing. And it's creating an extra dimension of size 1 so that the arrays become compatible for multiplication. It's called broadcasting. — Ananda

Manuel Oliveira Manuel Oliveira · Accepted Answer · 2020-11-19T11:54:29

I tried to use the transpose definition to change multiple axes.

It works for 'ij,kl -> ijkl' , 'ik,jl->ijkl' ,'kl,ij->ijkl' but fails for 'il,jk->ijkl', 'jl,ik->ijkl'and 'jk,il->ijkl'

import numpy as np

a=np.eye(3)
a[0][0]=2
a[0][-1]=3
a[-1][-1]=4

b=np.eye(3)

def permutation(str_,Arr):
    
    Arr=np.reshape(Arr,[3,3,3,3])
   
    def splitString(str_):
        tmp1=str_.split(',')
        tmp2=tmp1[1].split('->')
        str_idx1=tmp1[0]
        str_idx2=tmp2[0]
        str_idx_out=tmp2[1]
        return str_idx1,str_idx2, str_idx_out
    
    
    idx_a, idx_b, idx_out=splitString(str_)
    
    dict_={'i':0,'j':1,'k':2,'l':3}
    
    def split(word): 
        return [char for char in word]
                
    
    a,b=split(idx_a)
    c,d=split(idx_b)
    
    Arr=np.transpose(Arr,(dict_[a],dict_[b],dict_[c],dict_[d]))

    return Arr


str_='jk,il->ijkl'

d=np.outer(a,b)
f=np.einsum(str_, a,b)

check=permutation(str_,d)

if (np.count_nonzero(f-check)==0):
    print ('Code is working!')

else:
    print("Something is wrong...")

Appreciate your suggestions!

Explination of numpy's einsum

2 Answers