Мне нужно взять матричное произведение двух матриц NumPy (или других 2d-массивов), содержащих логарифмические вероятности. Наивный способ np.log(np.dot(np.exp(a), np.exp(b)))
не является предпочтительным по очевидным причинам.
Использование
from scipy.misc import logsumexp
res = np.zeros((a.shape[0], b.shape[1]))
for n in range(b.shape[1]):
# broadcast b[:,n] over rows of a, sum columns
res[:, n] = logsumexp(a + b[:, n].T, axis=1)
работает, но работает примерно в 100 раз медленнее, чем np.log(np.dot(np.exp(a), np.exp(b)))
Использование
logsumexp((tile(a, (b.shape[1],1)) + repeat(b.T, a.shape[0], axis=0)).reshape(b.shape[1],a.shape[0],a.shape[1]), 2).T
или другие комбинации плитки и перестройки также работают, но работают еще медленнее, чем цикл выше из-за непомерно больших объемов памяти, необходимых для матриц входного сигнала с реалистичным размером.
В настоящее время я рассматриваю возможность написания расширения NumPy в C, чтобы вычислить это, но, конечно, я бы предпочел избежать этого. Есть ли установленный способ сделать это, или кто-нибудь знает об менее интенсивном в памяти способе выполнения этого вычисления?
EDIT: Благодаря larsmans для этого решения (см. Ниже для вывода):
def logdot(a, b):
max_a, max_b = np.max(a), np.max(b)
exp_a, exp_b = a - max_a, b - max_b
np.exp(exp_a, out=exp_a)
np.exp(exp_b, out=exp_b)
c = np.dot(exp_a, exp_b)
np.log(c, out=c)
c += max_a + max_b
return c
Быстрое сравнение этого метода с вышеописанным методом (logdot_old
) с использованием функции iPython magic %timeit
дает следующее:
In [1] a = np.log(np.random.rand(1000,2000))
In [2] b = np.log(np.random.rand(2000,1500))
In [3] x = logdot(a, b)
In [4] y = logdot_old(a, b) # this takes a while
In [5] np.any(np.abs(x-y) > 1e-14)
Out [5] False
In [6] %timeit logdot_old(a, b)
1 loops, best of 3: 1min 18s per loop
In [6] %timeit logdot(a, b)
1 loops, best of 3: 264 ms per loop
Очевидно, метод larsmans уничтожает мой!