【问题标题】:Combinatorial optimisation of a distance metric距离度量的组合优化
【发布时间】:2011-02-19 06:17:25
【问题描述】:

我有一组轨迹,由沿轨迹的点组成,并具有与每个点关联的坐标。我将这些存储在 3d 数组中(轨迹、点、参数)。我想找到在这些轨迹的可能成对组合之间具有最大累积距离的一组 r 轨迹。我认为正在工作的第一次尝试如下所示:

max_dist = 0
for h in itertools.combinations ( xrange(num_traj), r):
    for (m,l) in itertools.combinations (h, 2):
        accum = 0.
        for ( i, j ) in itertools.izip ( range(k), range(k) ):
            A = [ (my_mat[m, i, z] - my_mat[l, j, z])**2 \
                    for z in xrange(k) ]
            A = numpy.array( numpy.sqrt (A) ).sum()
            accum += A
    if max_dist < accum:
        selected_trajectories = h

这需要很长时间,因为 num_traj 可以在 500-1000 左右,而 r 可以在 5-20 左右。 k 是任意的,但通常可以达到 50。

为了超级聪明,我将所有内容都放入了两个嵌套列表推导式中,大量使用了 itertools:

chunk = [[ numpy.sqrt((my_mat[m, i, :] - my_mat[l, j, :])**2).sum() \
        for ((m,l),i,j) in \
        itertools.product ( itertools.combinations(h,2), range(k), range(k)) ]\
        for h in itertools.combinations(range(num_traj), r) ]

除了不可读(!!!)之外,它还需要很长时间。任何人都可以提出任何改进这一点的方法吗?

【问题讨论】:

    标签: python algorithm numpy combinatorics itertools


    【解决方案1】:

    这听起来像是一个“加权集团”问题:找到例如 网络中 r=5 人,具有最大兼容性/C(5,2) 对权重的最大总和。
    谷歌“加权集团”算法 -“集团渗透”→ 3k 点击。
    但我会选择贾斯汀皮尔的方法 因为它是可以理解和可控的
    (取n2个最好的对,从中选出最好的n3个三元组...... 调整 n2 n3 ... 以轻松权衡运行时间/结果质量。)

    添加于 5 月 18 日,随后是一个实施的削减。
    @Jose,看看什么 nbest[] 序列适合你会很有趣。

    #!/usr/bin/env python
    """ cliq.py: grow high-weight 2 3 4 5-cliques, taking nbest at each stage
        weight ab = dist[a,b] -- a symmetric numpy array, diag << 0
        weight abc, abcd ... = sum weight all pairs
        C[2] = [ (dist[j,k], (j,k)) ... ]  nbest[2] pairs
        C[3] = [ (cliqwt(j,k,l), (j,k,l)) ... ]  nbest[3] triples
        ...
        run time ~ N * (N + nbest[2] + nbest[3] ...)
    
    keywords: weighted-clique heuristic python
    """
    # cf "graph clustering algorithm"
    
    from __future__ import division
    import numpy as np
    
    __version__ = "denis 18may 2010"
    me = __file__.split('/') [-1]
    
    def cliqdistances( cliq, dist ):
        return sorted( [dist[j,k] for j in cliq  for k in cliq if j < k], reverse=True )
    
    def maxarray2( a, n ):
        """ -> max n [ (a[j,k], (j,k)) ...]  j <= k, a symmetric """
        jkflat = np.argsort( a, axis=None )[:-2*n:-1]
        jks = [np.unravel_index( jk, a.shape ) for jk in jkflat]
        return [(a[j,k], (j,k)) for j,k in jks if j <= k] [:n]
    
    def _str( iter, fmt="%.2g" ):
        return " ".join( fmt % x  for x in iter )
    
    #...............................................................................
    
    def maxweightcliques( dist, nbest, r, verbose=10 ):
    
        def cliqwt( cliq, p ):
            return sum( dist[c,p] for c in cliq )  # << 0 if p in c
    
        def growcliqs( cliqs, nbest ):
            """ [(cliqweight, n-cliq) ...] -> nbest [(cliqweight, n+1 cliq) ...] """
                # heapq the nbest ? here just gen all N * |cliqs|, sort
            all = []
            dups = set()
            for w, c in cliqs:
                for p in xrange(N):
                        # fast gen [sorted c+p ...] with small sorted c ?
                    cp = c + [p]
                    cp.sort()
                    tup = tuple(cp)
                    if tup in dups:  continue
                    dups.add( tup )
                    all.append( (w + cliqwt(c, p), cp ))
            all.sort( reverse=True )
            if verbose:
                print "growcliqs: %s" % _str( w for w,c in all[:verbose] ) ,
                print " best: %s" % _str( cliqdistances( all[0][1], dist )[:10])
            return all[:nbest]
    
        np.fill_diagonal( dist, -1e10 )  # so cliqwt( c, p in c ) << 0
        C = (r+1) * [(0, None)]  # [(cliqweight, cliq-tuple) ...]
            # C[1] = [(0, (p,)) for p in xrange(N)]
        C[2] = [(w, list(pair)) for w, pair in maxarray2( dist, nbest[2] )]
        for j in range( 3, r+1 ):
            C[j] = growcliqs( C[j-1], nbest[j] )
        return C
    
    #...............................................................................
    if __name__ == "__main__":
        import sys
    
        N = 100
        r = 5  # max clique size
        nbest = 10
        verbose = 0
        seed = 1
        exec "\n".join( sys.argv[1:] )  # N= ...
        np.random.seed(seed)
        nbest = [0, 0, N//2] + (r - 2) * [nbest]  # ?
    
        print "%s  N=%d  r=%d  nbest=%s"  % (me, N, r, nbest)
    
            # random graphs w cluster parameters ?
        dist = np.random.exponential( 1, (N,N) )
        dist = (dist + dist.T) / 2
        for j in range( 0, N, r ):
            dist[j:j+r, j:j+r] += 2  # see if we get r in a row
        # dist = np.ones( (N,N) )
    
        cliqs = maxweightcliques( dist, nbest, r, verbose )[-1]  # [ (wt, cliq) ... ]
    
        print "Clique weight,  clique,  distances within clique"
        print 50 * "-"
        for w,c in cliqs:
            print "%5.3g  %s  %s" % (
                w, _str( c, fmt="%d" ), _str( cliqdistances( c, dist )[:10]))
    

    【讨论】:

      【解决方案2】:

      除了其他人提到的内容之外,这里还有一些兴趣点和建议。 (顺便说一句,mathmike 建议生成所有所有对距离的查找列表是您应该立即实施的建议。它消除了算法复杂度中的 O(r^2)。)

      首先是线条

      for ( i, j ) in itertools.izip ( range(k), range(k) ):
          A = [ (my_mat[m, i, z] - my_mat[l, j, z])**2 \
              for z in xrange(k) ]
      

      可以替换为

      for i in xrange(k):
          A = [ (my_mat[m, i, z] - my_mat[l, i, z])**2 \
              for z in xrange(k) ]
      

      因为 i 和 j 在每个循环中总是相同的。这里根本不需要使用 izip。

      第二,关于线路

      A = numpy.array( numpy.sqrt (A) ).sum()
      

      您确定这是您要计算的方式吗?可能是这样,但这让我觉得很奇怪,因为如果这更像是向量之间的欧几里德距离,那么这条线将是:

      A = numpy.sqrt (numpy.array( A ).sum())
      

      或者只是

      A = numpy.sqrt(sum(A))
      

      因为我认为将 A 转换为 numpy 数组以使用 numpy 的 sum 函数会比仅使用 Python 内置的 sum 函数要慢,但我可能错了。此外,如果它确实是您想要的欧几里德距离,那么您将通过这种方式减少 sqrt。

      第三,您是否意识到您可能会尝试迭代多少种潜在组合?对于 num_traj = 1000 和 r = 20 的最坏情况,据我估计,这大约是 6.79E42 个组合。这对于您当前的方法非常棘手。即使对于 num_traj = 500 和 r = 5 的最佳情况,这也是 1.28E12 组合,这是相当多的,但并非不可能。这是您在这里遇到的真正问题,因为根据 mathmike 的建议,我提到的前两点并不是很重要。

      那你能做什么?好吧,你需要更聪明一点。我还不清楚什么是一个很好的方法。我猜你需要以某种方式使算法启发式。我的一个想法是尝试一种带有启发式的动态编程方法。对于每个轨迹,您可以找到它与另一个轨迹的每一对的距离的总和或平均值,并将其用作适应度度量。在转向三重奏之前,可以放弃一些具有最低健身措施的轨迹。然后,您可以对三重奏做同样的事情:找到每个轨迹所涉及的所有三重奏(在剩余的可能轨迹中)的累积距离的总和或平均值,并将其用作适应度测量来决定在继续之前放弃哪些到四人一组。不能保证最优解,但应该是相当不错的,相信会大大降低解的时间复杂度。

      【讨论】:

        【解决方案3】:

        无论如何,它可能需要永远,因为您的算法大约需要 ~ O( C( N, r ) * r^2 ),其中 C( N, r ) 是 N 选择 r。对于较小的 r(或 N),这可能没问题,但如果您绝对需要找到最大值,而不是使用近似启发式,您应该尝试使用不同的策略进行分支和绑定。这可能适用于较小的 r,并且可以为您节省很多不必要的重新计算。

        【讨论】:

          【解决方案4】:

          您可以在距离计算中放弃平方根计算...最大和也将具有最大平方和,尽管这只会产生恒定的加速。

          【讨论】:

            【解决方案5】:

            您可以从计算所有轨迹对之间的距离开始,而不是按需重新计算每对轨迹之间的距离。您可以将它们存储在字典中并根据需要进行查找。

            这样,您的内循环 for (i,j) ... 将被替换为恒定时间查找。

            【讨论】:

            • 这是加速的主要来源!谢谢!
            猜你喜欢
            • 1970-01-01
            • 2021-09-08
            • 1970-01-01
            • 2015-03-13
            • 1970-01-01
            • 2015-08-25
            • 1970-01-01
            • 1970-01-01
            • 1970-01-01
            相关资源
            最近更新 更多