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.

So, you would have an array `M3` of shape `(100000,100000,100000)` as the input? Isn't that too big to fit into RAM? — Divakar, Jul 16 '16 at 18:10
@UlderiqueDemoitre If any of answers below has solved your question please consider [accepting it](http://meta.stackexchange.com/q/5234/179419) by clicking the check-mark. This indicates to the wider community that you've found a solution and gives some reputation to both the answerer and yourself. There is no obligation to do this. — Alicia Garcia-Raboso, Jul 25 '16 at 23:45

score 4 · Answer 1 · edited May 23 '17 at 12:01

Einstein summation convention? Use np.einsum!

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. 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.

Sorry, I made an edit to my post to avoid the swapaxes, both for performance and being concise. Could you please update the timings test? I am not expecting huge changes, but mostly to relate to the appropriate post. Thanks! — Divakar, Jul 16 '16 at 18:34
@Divakar: Done! There's a difference of a factor of 10 with your previous code. — Alicia Garcia-Raboso, Jul 16 '16 at 21:00

Divakar · Answer 2 · 2016-07-16T18:13:16.257

4

Here's an approach using a series of dot-products -

# Get partial products and thus reach to final output
p1 = np.tensordot(W,M3,axes=(0,0))
p2 = np.tensordot(p1,W,axes=(1,0))
X3out = np.tensordot(p2,W,axes=(1,0))

edited Jul 16 '16 at 18:13

answered Jul 16 '16 at 17:39

Divakar

218,885
19
262
358

Multilinear maps in python using numpy

2 Answers2