numpy有效计算parafac/CP产品

我有一组矩阵，它们共享相同数量的列并具有不同的行数。让我们称他们为A，B，C，d等，并让它们的尺寸是IaxK IbxK，IcxK等

我要的是有效地计算IaxIbxIc ......张量P定义如下： P（IA， ib，ic，id，ie，...）= \ sum_k A（ia，k）B（ib，k）C（ic，k）...

所以如果我有两个因素，用简单的矩阵产品。

当然，我可以在这个“手”，通过外积计算，是这样的：

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

不过，我想的东西太多计算效率更高，比如依靠内置numpy的功能，但我无法找到相关功能，有人可以帮助我吗？

干杯，感谢

来源

2012-12-07 antoine

谢谢大家多为你的答案，我已经花了一天就这个问题和我最终找到了解决办法，所以我在这里发布备案

这种解决方案需要numpy的1.6和利用einsum，这是强大的巫术

基本上，如果有系数= [A，B，C，d]其中A，B，C和d与相同数量的列的矩阵，那么你将计算使用PARAFAC模型：

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

所以，一条线！

在一般情况下，我结束了这一点：

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)

来源

2012-12-07 23:50:19 antoine

它确实是强大的伏都教，甚至运行的速度是我制作的 – Jaime

的两倍。你有没有比较过这个版本的速度与原来的？我尝试使用四个阵列形状（10,3），（24,3），（15,3）和（75,3）。您的原始版本大约需要2ms，使用'einsum'的版本大约需要7.5ms。 –

看起来爱因斯确实从多核架构中受益，而我的原始东西没有。此外，我已经通过实验注意到它的扩展性更好（真正感兴趣的情况相当于数千行的矩阵和50列的矩阵）。我会试试这个 – antoine

记住有了这样的外部产品克罗内克产品变相您的问题应该通过这个简单的功能来解决：

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需要真正的数和矩阵列表的列表。

由Jaime评论后编辑。

来源

2012-12-07 21:52:21 Gawron

不工作...如果你把它应用到尺寸的2个向量（5,8）和（4,8），你得到一个新的（20,64），然后尝试重塑成（5,4）......最好的情况是，在重塑之前你错过了总和步骤，尽管我不确定结果会是什么问 – Jaime

好吧，下面的工作。首先这是怎么回事的制定例子...

a = np.random.rand(5, 8) 
b = np.random.rand(4, 8) 
c = np.random.rand(3, 8) 
ret = np.ones(5,4,3,8) 
ret *= a.reshape(5,1,1,8) 
ret *= b.reshape(1,4,1,8) 
ret *= c.reshape(1,1,3,8) 
ret = ret.sum(axis=-1)

和拥有全部功能

def tensor(elems) : 
    cols = elems[0].shape[-1] 
    n_elems = len(elems) 
    ret = np.ones(tuple([j.shape[0] for j in elems] + [cols])) 
    for j,el in enumerate(elems) : 
     ret *= el.reshape((1,) * j + (el.shape[0],) + 
          (1,) * (len(elems) - j - 1) + (cols,)) 
    return ret.sum(axis=-1)

来源

2012-12-07 21:54:49 Jaime

numpy有效计算parafac/CP产品

回答

相关问题