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!



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).







logarithm