Комбинаторная оптимизация метрики расстояния

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, и г может быть около 5-20. к произвольно, но обычно может быть до 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

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

Answer 4

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

, потому что я и J всегда одинаковы в каждом цикле. Здесь нет необходимости использовать izip.

Во-вторых, относительно линии

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

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

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

или просто

A = numpy.sqrt(sum(A)) 

, потому что я думаю, что преобразование А к массиву numpy для использования функции суммы numpy будет медленнее, чем просто использовать встроенную функцию суммы Python, но я могу ошибаться. Кроме того, если это действительно евклидово расстояние, которое вы хотите, тогда вы будете делать меньше sqrt таким образом.

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

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

Answer 5

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

Добавлено 18may, после этого вырезается реализация.
@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]))