Комбинаторная оптимизация метрики расстояния - PullRequest
1 голос
/ 13 мая 2010

У меня есть набор траекторий, состоящий из точек вдоль траектории и с координатами, связанными с каждой точкой. Я храню их в трехмерном массиве (траектория, точка, парам). Я хочу найти множество 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) ]

Помимо того, что это довольно нечитабельно (!!!), это также занимает много времени. Кто-нибудь может предложить какие-либо способы улучшить это?

Ответы [ 5 ]

3 голосов
/ 13 мая 2010

Вместо того, чтобы пересчитывать расстояние между каждой парой траекторий по запросу, вы можете начать с вычисления расстояния между всеми парами траекторий. Вы можете хранить их в словаре и искать их по мере необходимости.

Таким образом, ваш внутренний цикл for (i,j) ... будет заменен поиском с постоянным временем.

2 голосов
/ 14 мая 2010

Вот несколько интересных мест и предложений в дополнение к тому, что все остальные упоминали. (Между прочим, предложение 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 будет медленнее, чем просто использование встроенной функции суммы Python, но я могу ошибаться. Кроме того, если вам действительно нужно евклидово расстояние, то вы будете делать меньше sqrt таким образом.

В-третьих, вы понимаете, сколько потенциальных комбинаций вы пытаетесь повторить? Для худшего случая с num_traj = 1000 и r = 20, по моей оценке, это примерно 6,79E42 комбинаций. Это довольно сложно с вашим текущим методом. Даже для наилучшего случая num_traj = 500 и r = 5 это комбинации 1.28E12, что довольно мало, но не невозможно. Это реальная проблема, с которой вы столкнулись, потому что, следуя советам математиков, первые два момента, которые я упомянул, не очень важны.

Что вы можете сделать тогда? Ну, тебе нужно быть немного умнее. Мне пока не ясно, что было бы отличным методом для этого. Я предполагаю, что вам нужно каким-то образом сделать алгоритм эвристическим. Одна мысль у меня была в том, чтобы попробовать метод динамического программирования с эвристическим подходом. Для каждой траектории вы можете найти сумму или среднее значение расстояний для каждой ее пары с другой траекторией и использовать ее в качестве меры пригодности. Некоторые из траекторий с наименьшими показателями пригодности могут быть отброшены перед переходом к трио. Затем вы можете сделать то же самое с трио: найти сумму или среднее значение накопленных расстояний для всех трио (среди оставшихся возможных траекторий), с которыми связана каждая траектория, и использовать это в качестве меры пригодности, чтобы решить, какие из них отбрасывать, прежде чем двигаться дальше до четверки. Это не гарантирует оптимального решения, но оно должно быть достаточно хорошим и значительно снизит временную сложность решения, на мой взгляд.

2 голосов
/ 13 мая 2010

Вы можете отказаться от вычисления квадратного корня при расчете расстояния ... максимальная сумма также будет иметь максимальную квадратную сумму, хотя это только приводит к постоянному ускорению.

1 голос
/ 15 мая 2010

Это звучит как проблема "взвешенной клики": например, найдите. r = 5 человек в сети с максимальной совместимостью / максимальной суммой весов C (5,2) пар.
Алгоритм "взвешенной клики" Google - "перколяция кликов" & rarr; 3k хитов.
НО я бы пошел с методом Джастина Пила потому что это понятно и контролируемо
(возьмите n2 лучшие пары, из них лучшие n3 тройки ... отрегулируйте n2 n3 ... для легкого компромисса времени выполнения / качества результатов.)

Добавлено 18 мая, после реализации следует сокращение.
@ Хосе, было бы интересно посмотреть, какая последовательность 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]))
1 голос
/ 14 мая 2010

Скорее всего, это будет длиться вечно, так как ваш алгоритм занимает около ~ O( C( N, r ) * r^2 ), где C( N, r ) - это N, выберите r. Для меньших r (или N) это может быть хорошо, но если вам абсолютно необходимо найти максимум, а не использовать эвристику приближения, вы должны попробовать переходить и связываться с различными стратегиями. Это может сработать для меньших r и сэкономит вам немного на ненужных пересчетах.

...