python - una - Numpy dot demasiado inteligente sobre multiplicaciones simétricas

transpuesta de una matriz en python numpy

Question

¿Alguien sabe acerca de la documentación de este comportamiento?

import numpy as np
A = np.random.uniform(0,1,(10,5))
w = np.ones(5)
Aw = A*w
Sym1 = Aw.dot(Aw.T)
Sym2 = (A*w).dot((A*w).T)
diff = Sym1 - Sym2

diff.max () está cerca de la precisión de la máquina, no es cero , por ejemplo, 4.4e-16.

Esto (la discrepancia de 0) suele ser bueno ... en un mundo de precisión finita no debería sorprendernos.

Además, supongo que Numpy está siendo inteligente con los productos simétricos, para ahorrar flops y garantizar una salida simétrica ...

Pero trato con sistemas caóticos, y esta pequeña discrepancia se nota rápidamente cuando se depura . Entonces me gustaría saber exactamente qué está pasando.

Answer 1

Este comportamiento es el resultado de un cambio introducido para NumPy 1.11.0, en la solicitud de extracción # 6932 . De las notas de la versión para 1.11.0:

Anteriormente, las operaciones de gemm BLAS se usaban para todos los productos de matriz. Ahora, si el producto de la matriz se encuentra entre una matriz y su transposición, utilizará las operaciones BLAS de jarabe para aumentar el rendimiento. Esta optimización se ha extendido a @, numpy.dot, numpy.inner y numpy.matmul.

En los cambios para ese PR, uno encuentra este comentario :

/*
* Use syrk if we have a case of a matrix times its transpose.
* Otherwise, use gemm for all other cases.
*/

Entonces NumPy está haciendo una comprobación explícita para el caso de una matriz multiplicada por su transposición, y llamando a una función BLAS subyacente diferente en ese caso. Como @hpaulj notas en un comentario, tal comprobación es barata para NumPy, ya que una matriz 2d transpuesta es simplemente una vista en la matriz original, con forma invertida y pasos, por lo que basta con verificar algunas piezas de metadatos en las matrices ( en lugar de tener que comparar los datos de la matriz real).

Aquí hay un caso un poco más simple que muestra la discrepancia. Tenga en cuenta que usar una .copy en uno de los argumentos para dot es suficiente para vencer a la envoltura especial de NumPy.

import numpy as np
random = np.random.RandomState(12345)
A = random.uniform(size=(10, 5))
Sym1 = A.dot(A.T)
Sym2 = A.dot(A.T.copy())
print(abs(Sym1 - Sym2).max())

Supongo que una ventaja de esta carcasa especial, más allá del potencial obvio de aceleración, es que está garantizado (espero, pero en la práctica dependerá de la implementación BLAS) para obtener un resultado perfectamente simétrico cuando se usa syrk , en lugar de una matriz que es meramente simétrica hasta error numérico. Como una prueba (por cierto no muy buena) para esto, probé:

import numpy as np
random = np.random.RandomState(12345)
A = random.uniform(size=(100, 50))
Sym1 = A.dot(A.T)
Sym2 = A.dot(A.T.copy())
print("Sym1 symmetric: ", (Sym1 == Sym1.T).all())
print("Sym2 symmetric: ", (Sym2 == Sym2.T).all())

Resultados en mi máquina:

Sym1 symmetric: True
Sym2 symmetric: False

Answer 2

Sospecho que esto tiene que ver con la promoción de registros de punto flotante intermedio con una precisión de 80 bits. Un poco lo que confirma esta hipótesis es que si usamos menos flotadores siempre obtenemos 0 en nuestros resultados, ala

A = np.random.uniform(0,1,(4,2))
w = np.ones(2)
Aw = A*w
Sym1 = Aw.dot(Aw.T)
Sym2 = (A*w).dot((A*w).T)
diff = Sym1 - Sym2
# diff is all 0''s (ymmv)