python - Khatri product of matrices using np.tensordot -

i'm trying decompose tensor (m, n, o) matrices a(m, r), b (n, r) , c (k, r). known parafac decomposition. tensorly kind of decomposition.

an important step multiply a, b, , c tensor of shape (m, n, o).

tensorly follows:

def kt_to_tensor(a, b, c):     factors = [a, b, c]     r in range(factors[0].shape[1]):         vecs = np.ix_(*[u[:, r] u in factors])         if r:             res += reduce(np.multiply, vecs)         else:             res = reduce(np.multiply, vecs)     return res 

however, package i'm using (autograd) not support np.ix_ operations. wrote simpler definition follows:

def new_kt_to_tensor(a, b, c):     m, n, o = a.shape[0], b.shape[0], c.shape[0]     out = np.zeros((m, n, o))     k_max = a.shape[1]     alpha in range(0, m):         beta in range(0, n):             delta in range(0, o):                 k in range(0, k_max):                     out[alpha, beta, delta]=out[alpha, beta, delta]+ a[alpha, k]*b[beta, k]*c[delta, k]     return out 

however, turns out implementation has aspects autograd not support. however, autograd support np.tensordot.

i wondering how use np.tensordot obtain multiplication. think tensorflow's tf.tensordot have similar functionality.

intended solution should like:

def tensordot_multplication(a, b, c):     """     use np.tensordot     """ 

don't think np.tensordot here, needs spread-out axes don't participate in sum-reductions, have alignment requirement of keeping last axis aligned between 3 inputs while performing multiplication. thus, tensordot, need processing , have more memory requirements there.

i suggest 2 methods - 1 broadcasting , np.einsum.

approach #1 : broadcasting -

(a[:,none,none,:]*b[:,none,:]*c).sum(-1) 

explanation :

• extend a 4d, introducing new axes @ axis=(1,2) none/np.newaxis.

• similarly extend b 3d, introducing new axis @ axis=(1).

• keep c , perform elementwise multiplications resulting in 4d array.

• finally, sum-reduction comes in along last axis of 4d array.

schematically put -

a        : m        r b        :    n     r c        :       k  r  => a*b*c : m  n  k  r => out   : m  n  k    # (sum-reduction along last axis) 

approach #2 : np.einsum -

np.einsum('il,jl,kl->ijk',a,b,c) 

the idea same here previous broadcasting one, string notations helping out in conveying axes info in more concise manner.

broadcasting surely available on tensorflow has tools expand dimensions, whereas np.einsum not.