title: efficiently computing parafac / CP product in numpy
body: <p>This question focuses on numpy.</p>
  
  <p>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</p>
  
  <p>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)...</p>
  
  <p>So if I have two factors, I end up with simple matrix product.</p>
  
  <p>Of course I can compute this "by hand" through outer products, something like:</p>
  
  <pre><code>    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
  </code></pre>
  
  <p>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 ?</p>
  
  <p>Cheers, thanks</p>