Самый быстрый способ суммировать по верхним треугольным элементам с наименьшим количеством памяти

Question

Мне нужно выполнить суммирование вида i<j на симметричных матрицах.Это эквивалентно сумме по верхним треугольным элементам матрицы, исключая диагональ.

Учитывая A симметричный массив N x N, самое простое решение - np.triu(A,1).sum(), однако мне было интересно, существуют ли более быстрые методы, которыетребует меньше памяти.Кажется, что (A.sum() - np.diag(A).sum())/2 быстрее на большом массиве, но как избежать создания даже массива N x 1 из np.diag?Вдвойне вложенный цикл for не требует дополнительной памяти, но в Python это явно не тот путь.

Answer 1

Самый быстрый метод с наименьшим объемом памяти, в чистом виде - это сложение всей суммы и вычитание диагонали.

Это может показаться расточительным с точки зрения FLOPS, но обратите внимание, что теоретическая экономия относительнодля этой реализации есть только фактор 2. Если это что-то значит для вас, вам, вероятно, не следует использовать numpy в первую очередь.

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

Но даже этот «оптимизированный» цикл C может на практике получить задницу от A.sum ().Если A смежна, эта сумма потенциально может отправлять максимально оптимизированный кэш и SIMD-оптимизированный кодовый путь.Скорее всего, любой пользователь vanilly-C, который вы напишете, будет полностью уничтожен A.sum () в тесте.

Answer 2

Вы можете использовать обозначение Эйнштейна для суммирования по диагонали: np.einsum('ii', a) эквивалентно np.diag(a).sum().Для целей бенчмаркинга:

import numpy as np
a = np.arange(25).reshape(5, 5)
%timeit np.einsum('ii', a)
1.72 µs ± 88.6 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit np.diag(a).sum()
3.93 µs ± 29.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Answer 3

Добавляя мои 2 цента к идеям в других ответах и ​​комментариях, вы можете быть заинтересованы в следующем быстродействии для симметричной матрицы 1000x1000.Как видите, метод sum_diag в этом случае выигрывает незначительно.

import numpy as np

N = 1000
a = np.random.randint(-2000,2000,size=(N,N))
A = (a + a.T)/2

def sum_triu(A):
    return np.triu(A,1).sum()

def sum_diag(A):
    return (A.sum() - np.diag(A).sum())/2

def sum_trace(A):
    return (A.sum() - np.trace(A))/2

%timeit sum_triu(A)
# 3.65 ms ± 406 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit sum_diag(A)
# 663 µs ± 88.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit sum_trace(A)
# 732 µs ± 120 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Answer 4

Вы можете заменить np.diag(A).sum() на np.trace(A);это не создаст временный Nx1 массив