Multilinear maps in python using numpy
Question: In python I have a three dimensional array `T` with size `(n,n,n)`, and a 2-D
array, with size `(n,k)`.

In linear algebra, the multilinear map defined by `T` and applied to `W`, in
code would be:

    
    
    X3 = np.zeros((k,k,k))
        for i in xrange(k):
           for j in xrange(k):
                for t in xrange(k):
                   for l in xrange(n):
                     for m in xrange(n)
                        for h in xrange(n):
                            X3[i, j, t] += M3[l, m, h] * W[l, i] * W[m, j] * W[h, t]
    

See

<https://en.wikipedia.org/wiki/Multilinear_map>

For a reference.

This is very slow. I am wondering if there exist any alternative or any pre
build function in numpy that can speed up the operations.


Answer: Einstein summation convention? Use
[`np.einsum`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html)!

    
    
    X3 = np.einsum('lmh,li,mj,ht->ijt', M3, W, W, W)
    

EDIT:

Out of curiosity, I just ran some benchmarks comparing `np.einsum` to
Divakar's approach. The difference is hugely against `np.einsum`!

    
    
    import numpy as np
    
    def approach1(a, b):
        x = np.einsum('lmh,li,mj,ht->ijt', a, b, b, b)
    
    def approach2(a, b):
        p1 = np.tensordot(b, a, axes=(0,0))
        p2 = np.tensordot(p1, b, axes=(1,0))
        x = np.tensordot(p2, b, axes=(1,0))
    
    n = 100
    k = 10
    a = np.random.random((n, n, n))
    b = np.random.random((n, k))
    %timeit approach1(a, b)      # => 1 loop, best of 3: 26 s per loop
    %timeit approach2(a, b)      # => 100 loops, best of 3: 4.23 ms per loop
    

There's some discussion about it in [this
question](http://stackoverflow.com/questions/20149201/why-is-numpys-einsum-
slower-than-numpys-built-in-functions). It all seems to come down to the
generality that `np.einsum` tries to achieve -- at the cost of being able to
offload computations to low-level linear algebra packages.

