所有行对的快速差异与 NumPy

Question

我使用的算法要求每个示例都有一个矩阵，比如Xi，即ai x b，对于每个O(n^2) 示例对，我发现每一行之间的差异Xiu - Xjv ，然后对外积sum_u sum_v np.outer(Xiu - Xjv, Xiu - Xjv)求和。

不幸的是，这个内部双和相当慢，并导致运行时间在大型数据集上失控。现在我只是使用 for 循环来做到这一点。有没有一些pythonic方法来加速这个内部操作？我一直在想一个，但没有运气。

为了澄清，对于每个n 示例，都有一个矩阵Xi，其维度为ai x b，其中ai 对于每个示例都不同。对于每一对(Xi, Xj)，我想遍历两个矩阵之间的所有O(ai * bi)对行并找到Xiu - Xjv，将其外积与自身np.outer(Xiu - Xjv, Xiu - Xjv)相乘，最后将所有这些外积相加。

例如：假设 D 是 [[1,2],[3,4]]，我们只是对两个矩阵都使用它。

然后，即一步可能是np.outer(D[0] - D[1], D[0] - D[1])，这将是矩阵[4,4],[4,4]。

很简单，(0,0) 和 (1,1) 只是 0 个矩阵，而 (0,1) 和 (1,0) 都是 4 个矩阵，所以对所有四个外积的最终总和应该是[[8,8],[8,8]]。

Answer 1

好的，这很有趣。我仍然忍不住想，这一切都可以通过对numpy.tensordot 的一次巧妙调用来完成，但无论如何这似乎已经消除了所有 Python 级循环：

import numpy

def slow( a, b=None ):

    if b is None: b = a
    a = numpy.asmatrix( a )
    b = numpy.asmatrix( b )

    out = 0.0
    for ai in a:
        for bj in b:
            dij = bj - ai
            out += numpy.outer( dij, dij )
    return out

def opsum( a, b=None ):

    if b is None: b = a
    a = numpy.asmatrix( a )
    b = numpy.asmatrix( b )

    RA, CA = a.shape
    RB, CB = b.shape    
    if CA != CB: raise ValueError( "input matrices should have the same number of columns" )

    out = -numpy.outer( a.sum( axis=0 ), b.sum( axis=0 ) );
    out += out.T
    out += RB * ( a.T * a )
    out += RA * ( b.T * b )
    return out

def test( a, b=None ):
    print( "ground truth:" )
    print( slow( a, b ) )
    print( "optimized:" )
    print( opsum( a, b ) )  
    print( "max abs discrepancy:" )
    print( numpy.abs( opsum( a, b ) - slow( a, b ) ).max() )
    print( "" )

# OP example
test( [[1,2], [3,4]] )

# non-symmetric example
a = [ [ 1, 2, 3 ], [-4, 5, 6 ], [7, -8, 9 ], [ 10, 11, -12 ] ]
a = numpy.matrix( a, dtype=float )
b = a[ ::2, ::-1 ] + 15
test( a, b )

# non-integer example
test( numpy.random.randn( *a.shape ), numpy.random.randn( *b.shape ) )

使用该（相当随意的）示例输入，opsum 的计时（使用 IPython 中的 timeit opsum(a,b) 测量）看起来仅比 slow 好 3 到 5 倍。但当然它的扩展性要好得多：将数据点的数量扩大 100 倍，将特征的数量扩大 10 倍，然后我们已经看到了速度提高 10,000 倍。