¿Por qué np.dot es impreciso? (matrices n-dim)

Supongamos que tomamos np.dotdos 'float32'matrices 2D:

res = np.dot(a, b)   # see CASE 1
print(list(res[0]))  # list shows more digits

[-0.90448684, -1.1708503, 0.907136, 3.5594249, 1.1374011, -1.3826287]

Números. Excepto que pueden cambiar:

CASO 1 : rebanadaa

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(6, 6).astype('float32')

for i in range(1, len(a)):
    print(list(np.dot(a[:i], b)[0])) # full shape: (i, 6)

[-0.9044868,  -1.1708502, 0.90713596, 3.5594249, 1.1374012, -1.3826287]
[-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
[-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]

Los resultados difieren, a pesar de que el corte impreso se deriva exactamente de los mismos números multiplicados.

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(1, 6).astype('float32')

for i in range(1, len(a)):
    a_flat = np.expand_dims(a[:i].flatten(), -1) # keep 2D
    print(list(np.dot(a_flat, b)[0])) # full shape: (i*6, 6)

[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]

CASO 3 : control más fuerte; establecer todos los entrantes no involucrados a cero : agregar a[1:] = 0al código CASO 1. Resultado: persisten discrepancias.

CASO 4 : verificar índices distintos de [0]; por ejemplo [0], los resultados comienzan a estabilizar un número fijo de ampliaciones de matriz desde su punto de creación. Salida

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(6, 6).astype('float32')

for j in range(len(a) - 2):
    for i in range(1, len(a)):
        res = np.dot(a[:i], b)
        try:    print(list(res[j]))
        except: pass
    print()

Por lo tanto, para el caso 2D * 2D, los resultados difieren, pero son consistentes para 1D * 1D. Según algunas de mis lecturas, esto parece derivarse de 1D-1D usando una suma simple, mientras que 2D-2D usa una adición 'más elegante' que aumenta el rendimiento que puede ser menos precisa (por ejemplo, la adición por pares hace lo contrario). Sin embargo, no puedo entender por qué las discrepancias desaparecen en el caso 1 una vez que ase corta un "umbral" establecido; cuanto más grande ay más btarde, este umbral parece mentir, pero siempre existe.

Todos dijeron: ¿por qué es np.dotimpreciso (e inconsistente) para los arreglos ND-ND? Git relevante

Otros detalles :

• Entorno : sistema operativo Win-10, Python 3.7.4, Spyder 3.3.6 IDE, Anaconda 3.0 2019/10
• CPU : i7-7700HQ 2.8 GHz
• Numpy v1.16.5

Posible biblioteca culpable : Numpy MKL - también bibliotecas BLASS; gracias a Bi Rico por señalar

Código de prueba de esfuerzo : como se señaló, las discrepancias exacerban la frecuencia con matrices más grandes; si el anterior no es reproducible, el siguiente debería serlo (si no, pruebe con atenuaciones más grandes). Mi salida

np.random.seed(1)
a = (0.01*np.random.randn(9, 9999)).astype('float32') # first multiply then type-cast
b = (0.01*np.random.randn(9999, 6)).astype('float32') # *0.01 to bound mults to < 1

for i in range(1, len(a)):
    print(list(np.dot(a[:i], b)[0]))

Gravedad del problema : las discrepancias que se muestran son 'pequeñas', pero ya no lo son cuando se opera en una red neuronal con miles de millones de números multiplicados en unos pocos segundos y billones en todo el tiempo de ejecución; la precisión del modelo reportado difiere en 10s enteros de porcentajes, por este hilo .

A continuación se muestra un gif de matrices resultantes de alimentar a un modelo de lo que es, básicamente a[0], w / len(a)==1vs len(a)==32:

OTROS resultados de PLATAFORMAS , de acuerdo y gracias a las pruebas de Paul :

Caso 1 reproducido (en parte) :

• Google Colab VM - Intel Xeon 2.3 G-Hz - Jupyter - Python 3.6.8
• Escritorio Win-10 Pro Docker - Intel i7-8700K - jupyter / scipy-notebook - Python 3.7.3
• Ubuntu 18.04.2 LTS + Docker - AMD FX-8150 - jupyter / scipy-notebook - Python 3.7.3

Nota : estos producen un error mucho menor que el que se muestra arriba; dos entradas en la primera fila están apagadas en 1 en el dígito menos significativo de las entradas correspondientes en las otras filas.

Caso 1 no reproducido :

• Ubuntu 18.04.3 LTS - Intel i7-8700K - IPython 5.5.0 - Python 2.7.15+ y 3.6.8 (2 pruebas)
• Ubuntu 18.04.3 LTS - Intel i5-3320M - IPython 5.5.0 - Python 2.7.15+
• Ubuntu 18.04.2 LTS - AMD FX-8150 - IPython 5.5.0 - Python 2.7.15rc1

Notas :

• Los entornos de jupyter y cuaderno de notas de Colab vinculados muestran una discrepancia mucho menor (y solo para las dos primeras filas) que la observada en mi sistema. Además, el Caso 2 nunca (todavía) mostró imprecisión.
• Dentro de esta muestra muy limitada, el entorno Jupyter actual (Dockerized) es más susceptible que el entorno IPython.
• np.show_config()demasiado largo para publicar, pero en resumen: las envolturas de IPython están basadas en BLAS / LAPACK; Colab está basado en OpenBLAS. En IPython Linux envs, las bibliotecas BLAS están instaladas en el sistema; en Jupyter y Colab, provienen de / opt / conda / lib

ACTUALIZACIÓN : la respuesta aceptada es precisa, pero amplia e incompleta. La pregunta permanece abierta para cualquiera que pueda explicar el comportamiento a nivel de código , es decir, un algoritmo exacto utilizado por np.dot, y cómo explica las 'inconsistencias consistentes' observadas en los resultados anteriores (también vea los comentarios). Aquí hay algunas implementaciones directas más allá de mi desciframiento: sdot.c - arraytypes.c.src

Esto parece una imprecisión numérica inevitable. Como se explica aquí , NumPy utiliza un método BLAS altamente optimizado y cuidadosamente ajustado para la multiplicación de matrices . Esto significa que probablemente la secuencia de operaciones (suma y productos) seguida para multiplicar 2 matrices, cambia cuando cambia el tamaño de la matriz.

Intentando ser más claros, sabemos que, matemáticamente , cada elemento de la matriz resultante puede calcularse como el producto de puntos de dos vectores (secuencias de números de igual longitud). Pero no es así como NumPy calcula un elemento de la matriz resultante. De hecho, existen algoritmos más eficientes pero complejos, como el algoritmo de Strassen , que obtienen el mismo resultado sin calcular directamente el producto de punto fila-columna.

Cuando se utilizan dichos algoritmos, incluso si el elemento C ij de una matriz resultante C = AB se define matemáticamente como el producto de punto de la fila i-ésima de A con la columna j-ésima de B , si multiplica una matriz A2 que tiene misma fila i-ésima que A con una matriz B2 que tiene la misma columna j-ésima que B , el elemento C2 ij se calculará realmente siguiendo una secuencia diferente de operaciones (que depende de la totalidad de A2 y B2 matrices), posiblemente conduciendo a diferentes errores numéricos.

Es por eso que, incluso si matemáticamente C ij = C2 ij (como en su CASO 1), la diferente secuencia de operaciones seguida por el algoritmo en los cálculos (debido al cambio en el tamaño de la matriz) conduce a diferentes errores numéricos. El error numérico explica también los resultados ligeramente diferentes según el entorno y el hecho de que, en algunos casos, para algunos entornos, el error numérico puede estar ausente.