用NumPy实现所有行对的快速差异

Question

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

不幸的是，这个内部双和相当慢，并导致运行时间螺旋失控的大型数据集。现在我只是用for循环来完成这个任务。有什么节能型的方法来加速这种内部操作吗？我一直在努力想出一个，却没有运气。

为了澄清，对于每个n示例，都有一个带有维度ai x b的矩阵Xi，其中每个示例的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倍。