Каков наилучший способ вычисления трассировки матричного продукта в numpy?

Если у меня есть массивы numpy A и B, тогда я могу вычислить след их матричного продукта с помощью:

tr = numpy.linalg.trace(A.dot(B))

Однако матричное умножение A.dot(B) излишне вычисляет все недиагональные элементы в матричном произведении, когда в трассировке используются только диагональные элементы. Вместо этого я мог бы сделать что-то вроде:

tr = 0.0
for i in range(n):
    tr += A[i, :].dot(B[:, i])

но это выполняет цикл в коде Python и не так очевидно, как numpy.linalg.trace.

Есть ли лучший способ вычислить трассировку матричного произведения массивов numpy? Какой самый быстрый или самый идиоматический способ сделать это?

Ответы

Ответ 1

Вы можете улучшить решение @Bill, уменьшив промежуточное хранилище только на диагональные элементы:

from numpy.core.umath_tests import inner1d

m, n = 1000, 500

a = np.random.rand(m, n)
b = np.random.rand(n, m)

# They all should give the same result
print np.trace(a.dot(b))
print np.sum(a*b.T)
print np.sum(inner1d(a, b.T))

%timeit np.trace(a.dot(b))
10 loops, best of 3: 34.7 ms per loop

%timeit np.sum(a*b.T)
100 loops, best of 3: 4.85 ms per loop

%timeit np.sum(inner1d(a, b.T))
1000 loops, best of 3: 1.83 ms per loop

Другой вариант - использовать np.einsum и не иметь явного промежуточного хранилища:

# Will print the same as the others:
print np.einsum('ij,ji->', a, b)

В моей системе он работает немного медленнее, чем при использовании inner1d, но может не выполняться для всех систем, см. этот вопрос:

%timeit np.einsum('ij,ji->', a, b)
100 loops, best of 3: 1.91 ms per loop

Ответ 2

Из Википедии вы можете рассчитать трассу, используя продукт hadamard (поэлементное умножение):

# Tr(A.B)
tr = (A*B.T).sum()

Я думаю, что это требует меньше вычислений, чем numpy.trace(A.dot(B)).

Изменить:

Пробежал несколько таймеров. Этот способ намного быстрее, чем при использовании numpy.trace.

In [37]: timeit("np.trace(A.dot(B))", setup="""import numpy as np;
A, B = np.random.rand(1000,1000), np.random.rand(1000,1000)""", number=100)
Out[38]: 8.6434469223022461

In [39]: timeit("(A*B.T).sum()", setup="""import numpy as np;
A, B = np.random.rand(1000,1000), np.random.rand(1000,1000)""", number=100)
Out[40]: 0.5516049861907959

Ответ 3

Обратите внимание, что один небольшой вариант - взять точечный продукт vec торизированных матриц. В python векторизация выполняется с помощью .flatten('F'). Это немного медленнее, чем брать сумму продукта Адамара, на моем компьютере, поэтому это худшее решение, чем wflynny, но я думаю, что это интересно, поскольку в некоторых ситуациях, на мой взгляд, это может быть более интуитивно понятным. Например, лично я нахожу, что для нормального распределения матрицы для меня легче понять векторизованное решение.

Сравнение скорости в моей системе:

import numpy as np
import time

N = 1000

np.random.seed(123)
A = np.random.randn(N, N)
B = np.random.randn(N, N)

tart = time.time()
for i in range(10):
    C = np.trace(A.dot(B))
print(time.time() - start, C)

start = time.time()
for i in range(10):
    C = A.flatten('F').dot(B.T.flatten('F'))
print(time.time() - start, C)

start = time.time()
for i in range(10):
    C = (A.T * B).sum()
print(time.time() - start, C)

start = time.time()
for i in range(10):
    C = (A * B.T).sum()
print(time.time() - start, C)

Результат:

6.246593236923218 -629.370798672
0.06539678573608398 -629.370798672
0.057890892028808594 -629.370798672
0.05709719657897949 -629.370798672