Избегайте вложенных циклов, использующих списки и / или карту

Question

В течение нескольких дней я боролся с тем, как оптимизировать (не только сделать его более приятным) 3 вложенных циклов , содержащих условных и вызов функции внутри. Сейчас у меня есть следующее:

def build_prolongation_operator(p,qs):
    '''
    p: dimension of the coarse basis
    q: dimension of the fine basis

    The prolongation operator describes the relationship between
    the coarse and fine bases:    
    V_coarse = np.dot(V_fine, I)
    '''

    q = sum(qs)

    I = np.zeros([q, p])

    for i in range(0, q):
        for j in range(0, p):
            for k in range(0, qs[j]):
                # if BV i is a child of j, we set I[i, j] = 1
                if i == f_map(j, k, qs):
                    I[i, j] = 1
                    break

    return I

, где f_map:

def f_map(i, j, q):
    '''
    Mapping which returns the index k of the fine basis vector which
    corresponds to the jth child of the ith coarse basis vector.    
    '''

    if j < 0 or j > q[i]:
        print('ERROR in f_map')
        return None

    result = j

    for k in range(0, i):
        result += q[k]

    return result

При профилировании всего моего кода я получаю, что build_prolongation_operator вызывается 45 раз и f_map приблизительно 8,5 миллиона раз !!

Вот картинка:

Я пытался сделать то же самое с пониманием списка и картой, но безуспешно.

Вот пример входных данных, которые build_prolongation_operator ожидает:

p = 10
qs = randint(3, size=p)

Answer 1

Вот оптимизированный цикл, предложенный Раунаком Джайном в их ответом

for j in range(0,p):
    for k in range(0, qs[j]):
        # if BV i is a child of j, we set I[i,j] = 1
            val = f_map(j,k,qs)
            if I[val, j] == 0:
                I[val, j] = 1

и здесь это функция f_map, где я отредактировал имена аргументов, чтобы отразить имена, используемые вызывающей стороной

def f_map(j,k,qs):
    if k < 0 or k > qs[j]:
        print('ERROR in f_map')
        return None
    result = k
    for i in range(0, j):
        result += qs[i]
    return result

Начнем с того, что он всегда равен 0 ≤ k < qs[j] из-за определения цикла в k, поэтому мы можем безопасно удалить проверку работоспособности и написать

def f_map(j,k,qs):
    result = k
    for i in range(0, j):
        result += q[i]
    return result

Теперь, это строка документа встроенного sum 1020 *

Подпись: сумма (повторяемая, начало = 0, /)
Строка документации:
Возвращает сумму начального значения (по умолчанию: 0) плюс итерируемое число

Когда итерация пуста, вернуть начальное значение.
Эта функция предназначена специально для использования с числовыми значениями и может отклонять нечисловые типы.
Тип: встроенная_функция_метод

Очевидно, что мы можем написать

def f_map(j,k,qs):
    return sum(qs[:j], k)

и очевидно также, что мы можем сделать без вызова функции

for j in range(0,p):
    for k in range(0, qs[j]):
        # if BV i is a child of j, we set I[i,j] = 1
            val = sum(qs[:j], k)
            if I[val, j] == 0:
                I[val, j] = 1

Вызов встроенного должен быть более эффективным, чем вызов функции и цикл, не так ли?

---

Обращаясь к замечанию Безумного Физика

Мы можем предварительно вычислить частичные суммы qs, чтобы получить дальнейшее ускорение

sqs = [sum(qs[:i]) for i in range(len(qs))] # there are faster ways...
...
for j in range(0,p):
    for k in range(0, qs[j]):
        # if BV i is a child of j, we set I[i,j] = 1
            val = k+sqs[j]
            if I[val, j] == 0:
                I[val, j] = 1

Answer 2

Попробуйте и проверьте, работает ли это,

for j in range(0,p):
    for k in range(0, qs[j]):
        # if BV i is a child of j, we set I[i,j] = 1
        val = f_map(j,k,qs)
        if I[val, j] == 0:
            I[val, j] = 1

Answer 3

Во-первых, вам действительно не нужен p в качестве параметра вашей функции: len(qs) нужно вызывать только один раз, и это очень дешево.Если ваши входные данные всегда являются пустым массивом (и в данных обстоятельствах нет причин, по которым это не должно быть), qs.size также подойдет.

Давайте начнем с переписывания f_map.В цикле есть только кумулятивная сумма qs (но начинающаяся с нуля), которую вы можете предварительно вычислить один раз (или хотя бы один раз за вызов внешней функции).

def f_map(i, j, cumsum_q):
    return j + cumsum_q[i]

Где cumsum_q будет определен в build_prolongation_operator как

cumsum_q = np.roll(np.cumsum(qs), 1)
cumsum_q[0] = 0

Я уверен, что вы могли бы оценить полезность наличия в f_map того же набора имен переменных, что и в build_prolongation_operator.Чтобы сделать это еще проще, мы можем просто полностью удалить f_map и использовать вместо этого выражение, которое оно представляет, в ваших условиях:

if i == k + cumsum_q[j]:
    I[i, j] = 1

Цикл над k тогда означает «если i равно * 1024»* для any k ", установите элемент в 1. Если мы переписываем условие как i - cumsum_q[j] == k, вы можете видеть, что нам вообще не нужен цикл над k.i - cumsum_q[j] будет равно некоторым k в диапазоне [0, qs[j]), если оно неотрицательно и строго меньше qs[j].Вы можете просто проверить

if i >= cumsum_q[j] and i - cumsum_q[j] < qs[j]:
    I[i, j] = 1

Это сокращает ваш цикл до одной итерации на элемент матрицы.Вы не можете сделать лучше, чем это:

def build_prolongation_operator_optimized(qs):
    '''
    The prolongation operator describes the relationship between
    the coarse and fine bases:    
    V_coarse = np.dot(V_fine, I)
    '''
    qs = np.asanyarray(qs)
    p = qs.size
    cumsum_q = np.roll(np.cumsum(qs), 1)
    q = cumsum_q[0]
    cumsum_q[0] = 0

    I = np.zeros([q, p])

    for i in range(0, q):
        for j in range(0, p):
            # if BV i is a child of j, we set I[i, j] = 1
            if 0 <= i - cumsum_q[j] < qs[j]:
                I[i, j] = 1
    return I

Теперь, когда вы теперь знаете формулу для каждой ячейки, вы можете просто вычислить целую матрицу для вас по существу в одной строке с использованием широковещания:

def build_prolongation_operator_numpy(qs):
    qs = np.asanyarray(qs)
    cumsum_q = np.roll(np.cumsum(qs), 1)
    q = cumsum_q[0]
    cumsum_q[0] = 0
    i_ = np.arange(q).reshape(-1, 1)  # Make this a column vector
    return (i_ >= cumsum_q) & (i_ - cumsum_q < qs)

Я запустил небольшую демонстрацию, чтобы (A) предлагаемые решения получили тот же результат, что и ваш оригинал, и (B) работали быстрее:

In [1]: p = 10
In [2]: q = np.random.randint(3, size=p)

In [3]: ops = (
...     build_prolongation_operator(p, qs),
...     build_prolongation_operator_optimized(qs),
...     build_prolongation_operator_numpy(qs),
...     build_prolongation_operator_RaunaqJain(p, qs),
...     build_prolongation_operator_gboffi(p, qs),
... )

In [4]: np.array([[(op1 == op2).all() for op1 in ops] for op2 in ops])
Out[4]: 
array([[ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True]])

In [5]: %timeit build_prolongation_operator(p, qs)
321 µs ± 890 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [6]: %timeit build_prolongation_operator_optimized(qs)
75.1 µs ± 1.85 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [7]: %timeit build_prolongation_operator_numpy(qs)
24.8 µs ± 77.7 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [8]: %timeit build_prolongation_operator_RaunaqJain(p, qs)
28.5 µs ± 1.55 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [9]: %timeit build_prolongation_operator_gboffi(p, qs)
31.8 µs ± 772 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [10]: %timeit build_prolongation_operator_gboffi2(p, qs)
26.6 µs ± 768 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Как вы можете видеть,самый быстрый вариант - полностью векторизованный, но параметры @ RaunaqJain и @ gboffi * идут очень близко за секунду.

Примечание

Мое векторизованное решение создает логический массив.Если вы не хотите этого, либо используйте I.astype(...) для преобразования в нужный тип dtype, либо просмотрите его как массив np.uint8: I.view(dtype=np.uint8).

Answer 4

Я не знаю о базах и операторах продолжения, но вы должны сосредоточиться на самом алгоритме.Это почти всегда хороший совет, когда речь идет об оптимизации.

Вот, наверное, суть - и если нет, то это то, с чего можно начать: вычисление f_map не зависит от i, но выпересчитать его для каждого значения i.Поскольку i находится в диапазоне от нуля до суммы значений в qs, вы сохраните огромное количество повторных вычислений, кэшируя результаты;Google "Python memoize", и он практически напишет сам.Исправьте это, и вы, вероятно, сделали это, без какой-либо микрооптимизации.

Вам понадобится достаточно места для хранения значений max(p) * max(qs[j]), но из числа вызовов, о которых вы сообщаете, это не должно быть проблемой,