python - functions - ¿Hay un método de numpy/scipy dot "mejorado"?

numpy python 3.6 windows

Question

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 :

1. Y = np.dot(np.dot(AT, np.diag(Q)), A) y
2. 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?

Answer 1

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,

Answer 2

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

Answer 3

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)