python - float - Resultado inesperado con+= en matrices NumPy

numpy array type float

Question

Estoy creando matrices / matrices simétricas en Python con NumPy, usando un método estándar:

x = rand(500,500)
x = (x+x.T)
all(x==x.T)
> True

Ahora seamos inteligentes

x = rand(500,500)
x += x.T
all(x==x.T)
> False

¿Esperar lo?

x==x.T
> array([[ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       ..., 
       [False, False, False, ...,  True,  True,  True],
       [False, False, False, ...,  True,  True,  True],
       [False, False, False, ...,  True,  True,  True]], dtype=bool)

Los segmentos superior izquierdo e inferior derecho son simétricos. ¿Qué pasa si elegí una matriz más pequeña?

x = rand(50,50)
x += x.T
all(x==x.T)
> True

DE ACUERDO....

x = rand(90,90)
x += x.T
all(x==x.T)
> True

x = rand(91,91)
x += x.T
all(x==x.T)
> False

Y solo para estar seguro ...

x = rand(91,91)
x = (x+x.T)
all(x==x.T)
> True

¿Se trata de un error o estoy a punto de aprender algo alocado sobre las matrices += y NumPy?

Answer 1

El detalle de implementación mencionado por otros se llama buffering . Puede leer más sobre esto en los documentos en iteración de matriz .

Si observas tu ejemplo fallido con un poco más de detalle:

>>> a = np.random.rand(91, 91)
>>> a += a.T
>>> a[:5, -1]
array([ 0.83818399,  1.06489316,  1.23675312,  0.00379798,  1.08967428])
>>> a[-1, :5]
array([ 0.83818399,  1.06489316,  1.75091827,  0.00416305,  1.76315071])

Entonces, el primer valor incorrecto es 90*91 + 2 = 8192 , lo cual, como era de esperar, es lo que obtenemos de:

>>> np.getbufsize()
8192

Y también podemos ponerlo más alto, y luego:

>>> np.setbufsize(16384)  # Must be a multiple of 16
8192  # returns the previous buffer size
>>> a = np.random.rand(91, 91)
>>> a += a.T
>>> np.all(a == a.T)
True

Aunque ahora:

>>> a = np.random.rand(129, 129)
>>> a += a.T
>>> np.all(a == a.T)
False

Esto es, por supuesto, una cosa peligrosa en la que confiar para la corrección, ya que de hecho es un detalle de implementación sujeto a cambios.

Answer 2

El problema es que la adición no se hace "a la vez"; xT es una vista de x para que comience a agregar a cada elemento de x , se muta xT . Esto ensucia las adiciones posteriores.

El hecho de que funcione para tamaños por debajo (91, 91) es casi definitivamente un detalle de implementación. Utilizando

x = numpy.random.rand(1000, 1000)
x += x.T.copy()
numpy.all(x==x.T)
#>>> True

arregla eso, pero entonces realmente no tienes ningún beneficio de espacio.

Answer 3

La operación de transpose devuelve una vista de la matriz, lo que significa que no se asigna ninguna nueva matriz. Lo que, a su vez, significa que está leyendo y modificando la matriz al mismo tiempo. Es difícil decir por qué funcionan algunos tamaños o algunas áreas del resultado, pero lo más probable es que tenga que ver con la forma en que Nympy se ocupa de la adición de matrices (tal vez hace copias de submatrices) y / o vistas de matrices (tal vez para tamaños pequeños crea) una nueva matriz).

La operación x = x + xT funciona porque allí estás creando una nueva matriz y luego asignando a x , por supuesto.