efficiently computing parafac / CP product in numpy

Question

This question focuses on numpy.

I have a set of matrices which all share the same number of columns and have different number of rows. Let's call them A, B, C, D, etc and let their dimensions be IaxK IbxK, IcxK, etc

What I want is to efficiently compute the IaxIbxIc... tensor P defined as follow: P(ia,ib,ic,id,ie,...)=\sum_k A(ia,k)B(ib,k)C(ic,k)...

So if I have two factors, I end up with simple matrix product.

Of course I can compute this "by hand" through outer products, something like:

    def parafac(factors,components=None):
        ndims = len(factors)
        ncomponents = factors[0].shape[1]
        total_result=array([])
        if components is None:
            components=range(ncomponents)

        for k in components:
            #for each component (to save memory)
            result = array([])
            for dim in range(ndims-1,-1,-1):
                #Augments model with next dimension
                current_dim_slice=[slice(None,None,None)]
                current_dim_slice.extend([None]*(ndims-dim-1))
                current_dim_slice.append(k)
                if result.size:
                    result = factors[dim].__getitem__(tuple(current_dim_slice))*result[None,...]
                else:
                    result = factors[dim].__getitem__(tuple(current_dim_slice))
            if total_result.size:
                total_result+=result
            else:
                total_result=result
        return total_result

Still, I would like something much more computationally efficient, like relying on builtin numpy functions, but I cannot find relevant functions, can someone help me ?

Cheers, thanks

Answer 1

Thank you all very much for your answers, I've spent the day on this and I eventually found the solution, so I post it here for the record

This solution requires numpy 1.6 and makes use of einsum, which is powerful voodoo magic

basically, if you had factor=[A,B,C,D] with A,B,C and D matrices with the same number of columns, then you would compute the parafac model using:

import numpy
P=numpy.einsum('az,bz,cz,dz->abcd',A,B,C,D)

so, one line!

In the general case, I end up with this:

def parafac(factors):
    ndims = len(factors)
    request=''
    for temp_dim in range(ndims):
        request+=string.lowercase[temp_dim]+'z,'
    request=request[:-1]+'->'+string.lowercase[:ndims]
    return einsum(request,*factors)

Answer 2

Having in mind that outer product is Kronecker product in disguise your problem should be solved by this simple functions:

def outer(vectors):
    shape=[v.shape[0] for v in vectors]
    return reduce(np.kron, vectors).reshape(shape)
def cp2Tensor(l,A):
    terms=[]    
    for r in xrange(A[0].shape[1]):
        term=l[r]*outer([A[n][:,r] for n in xrange(len(A))])
        terms.append(term)
    return sum(terms)

cp2Tensor takes list of real numbers and list of matrices.

Edited after comment by Jaime.