python3 - ¿Hay un producto de puntos numpy/scipy, que calcule solo las entradas diagonales del resultado?

python3 numpy

Question

Imagina tener 2 matrices numpy:

> A, A.shape = (n,p)
> B, B.shape = (p,p)

Normalmente, p es un número menor (p <= 200), mientras que n puede ser arbitrariamente grande.

Estoy haciendo lo siguiente:

result = np.diag(A.dot(B).dot(A.T))

Como puede ver, solo guardo las n entradas diagonales, sin embargo, hay una matriz intermedia (nxn) calculada a partir de la cual solo se guardan las entradas diagonales.

Deseo una función como diag_dot (), que solo calcula las entradas diagonales del resultado y no asigna la memoria completa.

Un resultado sería:

> result = diag_dot(A.dot(B), A.T)

¿Existe una funcionalidad prefabricada como esta y se puede hacer de manera eficiente sin la necesidad de asignar la matriz intermedia (nxn)?

Answer 1

Creo que lo conseguí por mi cuenta, pero sin embargo compartiré la solución:

Ya que obtenemos solo las diagonales de una multiplicación matricial.

> Z = N.diag(X.dot(Y))

es equivalente a la suma individual del producto escalar de filas de X y columnas de Y, la declaración anterior es equivalente a:

> Z = (X * Y.T).sum(-1)

Para las variables originales esto significa:

> result = (A.dot(B) * A).sum(-1)

Por favor, corrígeme si me equivoco, pero esto debería ser ...

Answer 2

Puede obtener casi cualquier cosa que haya soñado con numpy.einsum . Hasta que empiezas a entenderlo, básicamente parece un vudú negro ...

>>> a = np.arange(15).reshape(5, 3)
>>> b = np.arange(9).reshape(3, 3)

>>> np.diag(np.dot(np.dot(a, b), a.T))
array([  60,  672, 1932, 3840, 6396])
>>> np.einsum(''ij,ji->i'', np.dot(a, b), a.T)
array([  60,  672, 1932, 3840, 6396])
>>> np.einsum(''ij,ij->i'', np.dot(a, b), a)
array([  60,  672, 1932, 3840, 6396])

EDITAR Usted puede realmente obtener todo en un solo disparo, es ridículo ...

>>> np.einsum(''ij,jk,ki->i'', a, b, a.T)
array([  60,  672, 1932, 3840, 6396])
>>> np.einsum(''ij,jk,ik->i'', a, b, a)
array([  60,  672, 1932, 3840, 6396])

EDITAR Usted no quiere dejar que se calcule demasiado por sí solo aunque ... Agregó la respuesta del OP a su propia pregunta para comparación también.

n, p = 10000, 200
a = np.random.rand(n, p)
b = np.random.rand(p, p)

In [2]: %timeit np.einsum(''ij,jk,ki->i'', a, b, a.T)
1 loops, best of 3: 1.3 s per loop

In [3]: %timeit np.einsum(''ij,ij->i'', np.dot(a, b), a)
10 loops, best of 3: 105 ms per loop

In [4]: %timeit np.diag(np.dot(np.dot(a, b), a.T))
1 loops, best of 3: 5.73 s per loop

In [5]: %timeit (a.dot(b) * a).sum(-1)
10 loops, best of 3: 115 ms per loop

Answer 3

Una respuesta peatonal, que evita la construcción de grandes matrices intermedias es:

result=np.empty( [n.], dtype=A.dtype )
for i in xrange(n):
    result[i]=A[i,:].dot(B).dot(A[i,:])