vettoriale - python matrix product
modo numericamente stabile per moltiplicare le matrici di probabilità di log in numpy (2)
Ho bisogno di prendere il prodotto matrice di due matrici NumPy (o altri array 2d) che contengono probabilità di log. Il modo ingenuo np.log(np.dot(np.exp(a), np.exp(b))) non è preferito per ovvi motivi.
utilizzando
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)
funziona ma funziona circa 100 volte più lentamente di np.log(np.dot(np.exp(a), np.exp(b)))
utilizzando
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
o altre combinazioni di tile e risagoma funzionano ma funzionano anche più lentamente del ciclo precedente a causa delle quantità proibitive di memoria richieste per matrici di input di dimensioni realistiche.
Attualmente sto pensando di scrivere un'estensione NumPy in C per calcolarlo, ma ovviamente preferirei evitarlo. C'è un modo stabilito per farlo, o qualcuno sa di un modo meno intensivo di memoria per eseguire questo calcolo?
EDIT: grazie a larsman per questa soluzione (vedi sotto per la derivazione):
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
Un rapido confronto di questo metodo con il metodo pubblicato sopra ( logdot_old ) usando la funzione magic %timeit fornisce quanto segue:
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
Ovviamente il metodo dei larsman cancella il mio!
Stai accedendo alle colonne di res e b , che ha una scarsa localizzazione di riferimento . Una cosa da provare è archiviarli in ordine di colonna principale .
logsumexp funziona valutando il lato destro dell'equazione
log(∑ exp[a]) = max(a) + log(∑ exp[a - max(a)])
Vale a dire, estrae il massimo prima di iniziare a sommare, per evitare l'overflow in exp . Lo stesso può essere applicato prima di fare prodotti a punti vettoriali:
log(exp[a] ⋅ exp[b])
= log(∑ exp[a] × exp[b])
= log(∑ exp[a + b])
= max(a + b) + log(∑ exp[a + b - max(a + b)]) { this is logsumexp(a + b) }
ma prendendo una svolta diversa nella derivazione, otteniamo
log(∑ exp[a] × exp[b])
= max(a) + max(b) + log(∑ exp[a - max(a)] × exp[b - max(b)])
= max(a) + max(b) + log(exp[a - max(a)] ⋅ exp[b - max(b)])
La forma finale ha un prodotto punto vettoriale nelle sue viscere. Si estende anche facilmente alla moltiplicazione della matrice, quindi otteniamo l'algoritmo
def logdotexp(A, B):
max_A = np.max(A)
max_B = np.max(B)
C = np.dot(np.exp(A - max_A), np.exp(B - max_B))
np.log(C, out=C)
C += max_A + max_B
return C
Ciò crea due temporanei A -sized e due B -sized, ma uno di ciascuno può essere eliminato da
exp_A = A - max_A
np.exp(exp_A, out=exp_A)
e allo stesso modo per B (Se le matrici di input possono essere modificate dalla funzione, tutti i provvisori possono essere eliminati).
