在 numpy 中有效地计算 parafac / CP 产品答案

【问题标题】：efficiently computing parafac / CP product in numpy在 numpy 中有效地计算 parafac / CP 产品
【发布时间】：2012-12-07 10:23:55
【问题描述】：

本题关注numpy。

我有一组矩阵，它们共享相同的列数和不同的行数。我们称它们为 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 函数，但我找不到相关函数，有人可以帮助我吗？

干杯，谢谢

【问题讨论】：

标签： python numpy

【解决方案1】：

非常感谢大家的回答，我花了一天的时间终于找到了解决办法，所以我把它贴在这里记录下来

此解决方案需要 numpy 1.6 并使用 einsum，即强大的巫术魔法

基本上，如果您有 factor=[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)

【讨论】：

它确实是强大的巫术，甚至运行速度大约是我制作的两倍
不错。你有比较过这个版本的速度和你原来的速度吗？我尝试使用形状为 (10,3)、(24,3)、(15,3) 和 (75,3) 的四个数组。您的原始版本大约需要 2ms，使用 einsum 的版本大约需要 7.5ms。
看起来 einsum 确实受益于多核架构，而我原来的东西却没有。此外，我通过实验注意到它确实可以更好地扩展（真正感兴趣的案例是针对数千行和类似 50 列的矩阵）。我会试试这个
@antoine 你知道如何将其更改为非负约束吗？
很好的答案。谢谢！

【解决方案2】：

请记住，外部产品是伪装的克罗内克产品，您的问题应该通过这个简单的功能来解决：

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 评论后编辑。

【讨论】：

不起作用...如果将其应用于 2 个大小为 (5,8) 和 (4,8) 的向量，您将得到一个新的 (20, 64) 向量，然后尝试重塑为（5,4）...充其量，您在重塑之前错过了求和步骤，尽管我不太确定结果是否会是所要求的

【解决方案3】：

好的，所以下面的工作。首先是一个正在发生的事情的例子......

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)

【讨论】：