Перспективный подход и код решения
Для повышения эффективности памяти и повышения производительности мы могли бы использовать тензорное умножение матриц поэтапно.
Чтобы проиллюстрировать соответствующие шаги, давайте используем простейшее из решений с np.einsum от @ pv. -
np.einsum('ai,bj,ck,dl,abcd->ijkl', g, g, g, g, T)
Как видно, мы теряем первое измерение из g с тензорным умножением между его четырьмя вариантами и T.
Давайте сделаем это суммирование для умножения тензорной матрицы пошагово. Давайте начнем с первого варианта g и T:
p1 = np.einsum('abcd, ai->bcdi', T, g)
Таким образом, мы получаем тензор измерений в виде строковой записи: bcdi. Следующие шаги будут включать суммирование этого тензора по сравнению с остальными тремя g вариантами, используемыми в первоначальной реализации einsum. Следовательно, следующее сокращение будет -
p2 = np.einsum('bcdi, bj->cdij', p1, g)
Как видно, мы потеряли первые два измерения со строковыми обозначениями: a, b. Мы продолжаем это еще на два шага, чтобы избавиться от c и d тоже и останемся с ijkl в качестве конечного результата, вот так -
p3 = np.einsum('cdij, ck->dijk', p2, g)
p4 = np.einsum('dijk, dl->ijkl', p3, g)
Теперь мы могли бы использовать np.tensordot для этих сумм-сокращений, что было бы гораздо более эффективным.
Окончательная реализация
Таким образом, при переходе на np.tensordot мы получим итоговую реализацию примерно так -
p1 = np.tensordot(T,g,axes=((0),(0)))
p2 = np.tensordot(p1,g,axes=((0),(0)))
p3 = np.tensordot(p2,g,axes=((0),(0)))
out = np.tensordot(p3,g,axes=((0),(0)))
Испытание во время выполнения
Давайте протестируем все подходы на основе NumPy, опубликованные в других публикациях, чтобы решить проблему с производительностью.
Подходит как функции -
def rotT_Philipp(T, g): # @Philipp's soln
gg = np.outer(g, g)
gggg = np.outer(gg, gg).reshape(4 * g.shape)
axes = ((0, 2, 4, 6), (0, 1, 2, 3))
return np.tensordot(gggg, T, axes)
def rotT_Sven(T, g): # @Sven Marnach's soln
Tprime = T
for i in range(4):
slices = [None] * 4
slices[i] = slice(None)
slices *= 2
Tprime = g[slices].T * Tprime
return Tprime.sum(-1).sum(-1).sum(-1).sum(-1)
def rotT_pv(T, g): # @pv.'s soln
return np.einsum('ai,bj,ck,dl,abcd->ijkl', g, g, g, g, T)
def rotT_Divakar(T,g): # Posted in this post
p1 = np.tensordot(T,g,axes=((0),(0)))
p2 = np.tensordot(p1,g,axes=((0),(0)))
p3 = np.tensordot(p2,g,axes=((0),(0)))
p4 = np.tensordot(p3,g,axes=((0),(0)))
return p4
Сроки с исходными размерами набора данных -
In [304]: # Setup inputs
...: T = np.random.rand(3,3,3,3)
...: g = np.random.rand(3,3)
...:
In [305]: %timeit rotT(T, g)
...: %timeit rotT_pv(T, g)
...: %timeit rotT_Sven(T, g)
...: %timeit rotT_Philipp(T, g)
...: %timeit rotT_Divakar(T, g)
...:
100 loops, best of 3: 6.51 ms per loop
1000 loops, best of 3: 247 µs per loop
10000 loops, best of 3: 137 µs per loop
10000 loops, best of 3: 41.6 µs per loop
10000 loops, best of 3: 28.3 µs per loop
In [306]: 6510.0/28.3 # Speedup with the proposed soln over original code
Out[306]: 230.03533568904592
Как обсуждалось в начале этого поста, мы пытаемся добиться эффективности памяти и, следовательно, повышения производительности с ее помощью. Давайте проверим это по мере увеличения размеров набора данных -
In [307]: # Setup inputs
...: T = np.random.rand(5,5,5,5)
...: g = np.random.rand(5,5)
...:
In [308]: %timeit rotT(T, g)
...: %timeit rotT_pv(T, g)
...: %timeit rotT_Sven(T, g)
...: %timeit rotT_Philipp(T, g)
...: %timeit rotT_Divakar(T, g)
...:
100 loops, best of 3: 6.54 ms per loop
100 loops, best of 3: 7.17 ms per loop
100 loops, best of 3: 2.7 ms per loop
1000 loops, best of 3: 1.47 ms per loop
10000 loops, best of 3: 39.9 µs per loop