距离度量的组合优化

Question

我有一组轨迹，由轨迹上的点以及与每个点关联的坐标组成。我将它们存储在 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) ]

除了非常难以阅读（！！！）之外，它还采用了很长一段时间。谁能提出任何改进方法？

Answer 1

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

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

Answer 2

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

Answer 3

除了其他人提到的内容之外，这里还有一些兴趣点和建议。 （顺便说一句，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 函数，但我可能是错的。另外，如果它确实是您想要的欧几里德距离，那么您将以这种方式执行更少的开方。

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

那你能做什么呢？好吧，你需要更聪明一点。我还不清楚什么是一个很好的方法来实现这一点。我猜你需要以某种方式使算法启发式。我的一个想法是尝试一种带有启发式的动态编程方法。对于每个轨迹，您可以找到它与另一个轨迹的每个配对的距离总和或平均值，并将其用作适应度度量。在进入三重奏之前，可以放弃一些适应度最低的轨迹。然后，您可以对三重奏做同样的事情：找到每个轨迹所涉及的所有三重奏（在剩余的可能轨迹中）的累积距离的总和或平均值，并将其用作适应度度量来决定在继续之前丢弃哪些轨迹到四人一组。它不能保证最佳解决方案，但它应该相当好，并且我相信它会大大降低解决方案的时间复杂度。

Answer 4

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

Answer 5

这听起来像是一个“加权派系”问题：找到例如
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]))