python - tutorial - forma numéricamente estable para multiplicar matrices de probabilidad de registro en numpy

tamaño de un arreglo numpy

Question

Necesito tomar el producto matricial de dos matrices NumPy (u otras matrices de 2d) que contengan probabilidades de registro. La manera ingenua np.log(np.dot(np.exp(a), np.exp(b))) no es preferida por razones obvias.

Utilizando

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) 

funciona pero se ejecuta aproximadamente 100 veces más lento que np.log(np.dot(np.exp(a), np.exp(b)))

Utilizando

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

u otras combinaciones de mosaicos y remodelación también funcionan, pero se ejecutan incluso más lento que el ciclo anterior debido a las cantidades prohibitivamente grandes de memoria requeridas para matrices de entrada de tamaño realista.

Actualmente estoy considerando escribir una extensión NumPy en C para calcular esto, pero por supuesto prefiero evitar eso. ¿Hay una forma establecida de hacerlo, o alguien sabe de una manera menos intensiva en memoria de realizar este cálculo?

EDITAR: Gracias a larsmans por esta solución (ver a continuación la derivación):

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

Una comparación rápida de este método con el método publicado anteriormente ( logdot_old ) utilizando la función %timeit mágica de %timeit produce lo siguiente:

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

¡Obviamente el método de Larsman borra el mío!

Answer 1

Está accediendo a las columnas de res y b , que tiene poca localidad de referencia . Una cosa para intentar es almacenar estos en orden de columna mayor .

Answer 2

logsumexp funciona evaluando el lado derecho de la ecuación

log(∑ exp[a]) = max(a) + log(∑ exp[a - max(a)])

Es decir, saca el máximo antes de comenzar a sumar, para evitar el desbordamiento en exp . Lo mismo se puede aplicar antes de hacer productos vector dot:

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

pero al tomar un giro diferente en la derivación, obtenemos

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 final tiene un producto vector dot en sus entrañas. También se extiende fácilmente a la multiplicación de matrices, por lo que obtenemos el 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

Esto crea dos temporarios A y dos B -dimensionados, pero uno de cada uno puede ser eliminado por

exp_A = A - max_A
np.exp(exp_A, out=exp_A)

y de manera similar para B (Si las matrices de entrada pueden ser modificadas por la función, se pueden eliminar todos los temporales).