python - functions - ¿Hay un método de numpy/scipy dot "mejorado"?
numpy python 3.6 windows (3)
Problema
Me gustaría calcular lo siguiente usando numpy o scipy:
Y = A**T * Q * A
donde A es una matriz mxn , A**T es la transposición de A y Q es una matriz diagonal mxm .
Como Q es una matriz diagonal, almaceno solo sus elementos diagonales como un vector.
Formas de resolver para Y
Actualmente puedo pensar en dos formas de calcular Y :
-
Y = np.dot(np.dot(AT, np.diag(Q)), A)y -
Y = np.dot(AT * Q, A).
Claramente, la opción 2 es mejor que la opción 1, ya que no se debe crear una matriz real con diag(Q) (si esto es lo que numpy realmente hace ...)
Sin embargo, ambos métodos adolecen del defecto de tener que asignar más memoria de la que realmente es necesaria ya que AT * Q y np.dot(AT, np.diag(Q)) deben almacenarse junto con A para calcular Y
Pregunta
¿Hay algún método en numpy / scipy que elimine la asignación innecesaria de memoria extra donde solo pasaría dos matrices A y B (en mi caso B es AT ) y un vector de ponderación Q junto con ella?
Solo quería poner eso en SO, pero esta solicitud de extracción debería ser útil y eliminar la necesidad de una función separada para numpy.dot https://github.com/numpy/numpy/pull/2730 Esto debería estar disponible en numpy 1.7
Mientras tanto, utilicé el ejemplo anterior para escribir una función que puede reemplazar al punto numpy, cualquiera que sea el orden de tus matrices, y hacer la llamada correcta a fblas.dgemm. http://pastebin.com/M8TfbURi
Espero que esto ayude,
(w / r / t la última frase del OP: no conozco un método tan numpy / scipy pero no tengo la pregunta en el título OP (es decir, mejora el rendimiento del punto NumPy) lo que está debajo debería ser de algunos ayuda. En otras palabras, mi respuesta está dirigida a mejorar el rendimiento de la mayoría de los pasos que comprenden su función para Y).
En primer lugar, esto debería darle un impulso notable sobre el método vainilla NumPy dot :
>>> from scipy.linalg import blas as FB
>>> vx = FB.dgemm(alpha=1., a=v1, b=v2, trans_b=True)
Tenga en cuenta que las dos matrices, v1, v2 están ambas en orden C_FORTRAN
Puede acceder al orden de bytes de una matriz NumPy a través del atributo de banderas de una matriz, como sigue:
>>> c = NP.ones((4, 3))
>>> c.flags
C_CONTIGUOUS : True # refers to C-contiguous order
F_CONTIGUOUS : False # fortran-contiguous
OWNDATA : True
MASKNA : False
OWNMASKNA : False
WRITEABLE : True
ALIGNED : True
UPDATEIFCOPY : False
para cambiar el orden de una de las matrices para que ambas estén alineadas, solo llame al constructor de la matriz NumPy, pase la matriz y establezca la bandera de orden apropiada en True
>>> c = NP.array(c, order="F")
>>> c.flags
C_CONTIGUOUS : False
F_CONTIGUOUS : True
OWNDATA : True
MASKNA : False
OWNMASKNA : False
WRITEABLE : True
ALIGNED : True
UPDATEIFCOPY : False
Puede optimizar aún más explotando la alineación de orden de matriz para reducir el consumo de exceso de memoria causado por la copia de las matrices originales .
Pero, ¿por qué se copian las matrices antes de pasarlas a punto ?
El producto punto se basa en las operaciones BLAS. Estas operaciones requieren matrices almacenadas en orden C contiguo; esta restricción hace que las matrices se copien.
Por otro lado, la transposición no afecta a una copia, aunque desafortunadamente devuelve el resultado en orden de Fortran :
Por lo tanto, para eliminar el cuello de botella de rendimiento, debe eliminar el paso de copiar matriz de predicado ; hacer eso solo requiere pasar ambas matrices al punto en C-orden contiguo *.
Entonces, para calcular el punto (AT, A) sin hacer una copia adicional:
>>> import scipy.linalg.blas as FB
>>> vx = FB.dgemm(alpha=1.0, a=A.T, b=A.T, trans_b=True)
En resumen, la expresión anterior (junto con la instrucción de importación de predicados) puede sustituir al punto, para proporcionar la misma funcionalidad pero un mejor rendimiento
puedes vincular esa expresión a una función como esta:
>>> super_dot = lambda v, w: FB.dgemm(alpha=1., a=v.T, b=w.T, trans_b=True)
numpy.einsum es lo que estás buscando:
numpy.einsum(''ij, i, ik -> jk'', A, Q, A)
Esto no necesitará ninguna memoria adicional (aunque normalmente el einsum funciona más lento que las operaciones BLAS)