Как оптимизировать математические операции над матрицей в python

Я пытаюсь сократить время функции, которая выполняет серию вычислений с двумя матрицами. В поисках этого я слышал о NumPy, но я действительно не знаю, как применить это к моей проблеме. Кроме того, я думаю, что одна из вещей, замедляющая мою функцию - это наличие множества точек (я слышал об этом в этомэта страница ).

Математика соответствует факторизации для задачи квадратичного присваивания:

Мой код:

    delta = 0
    for k in xrange(self._tam):
        if k != r and k != s:
            delta +=
                self._data.stream_matrix[r][k] \
                * (self._data.distance_matrix[sol[s]][sol[k]] - self._data.distance_matrix[sol[r]][sol[k]]) + \
                self._data.stream_matrix[s][k] \
                * (self._data.distance_matrix[sol[r]][sol[k]] - self._data.distance_matrix[sol[s]][sol[k]]) + \
                self._data.stream_matrix[k][r] \
                * (self._data.distance_matrix[sol[k]][sol[s]] - self._data.distance_matrix[sol[k]][sol[r]]) + \
                self._data.stream_matrix[k][s] \
                * (self._data.distance_matrix[sol[k]][sol[r]] - self._data.distance_matrix[sol[k]][sol[s]])
    return delta

Выполнение этого для задачи размером 20 (матрица 20x20) занимает около 20 сегментов, узкое место в этой функции

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
303878   15.712    0.000   15.712    0.000 Heuristic.py:66(deltaC)

Я пытался подать заявкуmap в цикл for, но поскольку тело цикла не является вызовом функции, это невозможно.

Как я мог сократить время?

EDIT1

Чтобы ответить на комментарий Айкенберга:

sol является перестановкой, например [1,2,3,4]. функция вызывается, когда я генерирую соседние решения, поэтому сосед [1,2,3,4] равен [2,1,3,4]. Я меняю только две позиции в исходной перестановке и затем вызываюdeltaC, который вычисляет факторизацию решения с замененными позициями r, s (в приведенном выше примере r, s = 0,1). Эта перестановка сделана, чтобы избежать расчета полной стоимости соседнего решения. Я полагаю, я могу хранить значенияsol[k,r,s] в локальной переменной, чтобы избежать поиска ее значения в каждой итерации.Я не знаю, если это то, что вы спрашивали в своем комментарии.

EDIT2

Минимальный рабочий пример:

import random


distance_matrix = [[0, 12, 6, 4], [12, 0, 6, 8], [6, 6, 0, 7], [4, 8, 7, 0]]
stream_matrix = [[0, 3, 8, 3], [3, 0, 2, 4], [8, 2, 0, 5], [3, 4, 5, 0]]

def deltaC(r, s, S=None):
    '''
    Difference between C with values i and j swapped
    '''

    S = [0,1,2,3]

    if S is not None:
        sol = S
    else:
        sol = S

    delta = 0

    sol_r, sol_s = sol[r], sol[s]

    for k in xrange(4):
        if k != r and k != s:
            delta += (stream_matrix[r][k] \
                * (distance_matrix[sol_s][sol[k]] - distance_matrix[sol_r][sol[k]]) + \
                stream_matrix[s][k] \
                * (distance_matrix[sol_r][sol[k]] - distance_matrix[sol_s][sol[k]]) + \
                stream_matrix[k][r] \
                * (distance_matrix[sol[k]][sol_s] - distance_matrix[sol[k]][sol_r]) + \
                stream_matrix[k][s] \
                * (distance_matrix[sol[k]][sol_r] - distance_matrix[sol[k]][sol_s]))
    return delta


for _ in xrange(303878):
    d = deltaC(random.randint(0,3), random.randint(0,3))
print d

Теперь я думаю, что лучше использовать NumPy. Я пытался с Matrix (), но не улучшил производительность.

Лучшее решение найдено

Ну, наконец, я смог немного сократить время, комбинируя решение @ TooTone и сохраняя индексы в наборе, чтобы избежать if. Время сократилось с 18 до 8 секунд. Вот код:

def deltaC(self, r, s, sol=None):
    delta = 0
    sol = self.S if sol is None else self.S
    sol_r, sol_s = sol[r], sol[s]

    stream_matrix = self._data.stream_matrix
    distance_matrix = self._data.distance_matrix

    indexes = set(xrange(self._tam)) - set([r, s])

    for k in indexes:
        sol_k = sol[k]
        delta += \
            (stream_matrix[r][k] - stream_matrix[s][k]) \
            * (distance_matrix[sol_s][sol_k] - distance_matrix[sol_r][sol_k]) \
            + \
            (stream_matrix[k][r] - stream_matrix[k][s]) \
            * (distance_matrix[sol_k][sol_s] - distance_matrix[sol_k][sol_r])
    return delta

Я думаю, что для того, чтобы еще больше сократить время, лучше всего написать модуль.

Ответы на вопрос(1)

Ваш ответ на вопрос