Diferencia numérica Matlab/Octave/Numpy

Estoy portando un algoritmo que funciona en Matlab para nublarse y observé un comportamiento extraño. El segmento relevante de código esDiferencia numérica Matlab/Octave/Numpy

P = eye(4)*1e20; 
A = [1 -0.015 -0.025 -0.035; 0.015 1 0.035 -0.025; 0.025 -0.035 1 0.015; 0.035 0.025 -0.015 1]; 
V1 = A*(P*A') 
V2 = (A*P)*A'

Este código, cuando funciono con Matlab, proporciona las siguientes matrices:

V1 = 1.0021e+20   0 -8.0000e+00   0 
       0 1.0021e+20   0   0 
    -8.0000e+00   0 1.0021e+20   0 
       0   0   0 1.0021e+20 


V2 = 1.0021e+20   0 -8.0000e+00   0 
       0 1.0021e+20   0   0 
    -8.0000e+00   0 1.0021e+20   0 
       0   0   0 1.0021e+20

Nota que V1 y V2 son los mismos, como se esperaba.

Cuando el mismo código se ejecuta en Octave, que dispone lo siguiente:

V1 = 1.0021e+20 4.6172e+01 -1.3800e+02 1.8250e+02 
    -4.6172e+01 1.0021e+20 -1.8258e+02 -1.3800e+02 
    1.3801e+02 1.8239e+02 1.0021e+20 -4.6125e+01 
    -1.8250e+02 1.3800e+02 4.6125e+01 1.0021e+20 

V2 = 1.0021e+20 -4.6172e+01 1.3801e+02 -1.8250e+02 
    4.6172e+01 1.0021e+20 1.8239e+02 1.3800e+02 
    -1.3800e+02 -1.8258e+02 1.0021e+20 4.6125e+01 
    1.8250e+02 -1.3800e+02 -4.6125e+01 1.0021e+20

En numpy, el segmento se convierte en

from numpy import array, dot, eye 
A = numpy.array([[1, -0.015, -0.025, -0.035],[0.015, 1, 0.035, -0.025],[0.025, -0.035, 1, 0.015],[0.035, 0.025, -0.015, 1]]) 
P = numpy.eye(4)*1e20 
print numpy.dot(A,numpy.dot(P,A.transpose())) 
print numpy.dot(numpy.dot(A,P),A.transpose())

que da salida a

[[ 1.00207500e+20 4.61718750e+01 -1.37996094e+02 1.82500000e+02] 
[ -4.61718750e+01 1.00207500e+20 -1.82582031e+02 -1.38000000e+02] 
[ 1.38011719e+02 1.82386719e+02 1.00207500e+20 -4.61250000e+01] 
[ -1.82500000e+02 1.38000000e+02 4.61250000e+01 1.00207500e+20]] 
[[ 1.00207500e+20 -4.61718750e+01 1.38011719e+02 -1.82500000e+02] 
[ 4.61718750e+01 1.00207500e+20 1.82386719e+02 1.38000000e+02] 
[ -1.37996094e+02 -1.82582031e+02 1.00207500e+20 4.61250000e+01] 
[ 1.82500000e+02 -1.38000000e+02 -4.61250000e+01 1.00207500e+20]]

Así, tanto Octave y numpy proporciona la misma respuesta, pero es muy diferente de la de Matlab. El primer punto es que V1! = V2, que no parece correcto. El otro punto es que, aunque los elementos no diagonales son muchos órdenes de magnitud más pequeños que los diagonales, esto parece estar causando algún problema en mi algoritmo.

¿Alguien sabe cómo numpy y Octave se comporta de esta manera?

Fuente

2012-09-07 user1655812

Usan dobles internamente, que solo tienen unos 15 dígitos de precisión. Sus operaciones matemáticas probablemente exceden esto, lo que causa los resultados matemáticamente incorrectos.

vale la pena leer: http://floating-point-gui.de/

Editar: Desde el docs deduzco que no hay una precisión más alta disponible para Numpy. Parece que SymPy aunque may give you the needed precision - si esa biblioteca también funciona para usted.

Fuente

2012-09-07 22:16:26 Lucero

Eso no es del todo cierto. Hay un tipo de datos float128, sin embargo, creo que su precisión no siempre está bien definida. – seberg

@Sebastian, no encontré referencia alguna para un tipo float128 - only complex128 (porque esos son dos float64 vistos como un número con la parte real e imaginaria). http://docs.scipy.org/doc/numpy/user/basics.types.html – Lucero

Sí ... eso es porque float128 solo está disponible dependiendo de la computadora en la que se esté ejecutando. Pero en la PC habitual es. – seberg

Por si sirve de algo, yo obtener resultados idénticos a Matlab en un sistema de 64 bits:

[[ 1.00207500e+20 0.00000000e+00 -8.00000000e+00 0.00000000e+00] 
[ 0.00000000e+00 1.00207500e+20 0.00000000e+00 0.00000000e+00] 
[ -8.00000000e+00 0.00000000e+00 1.00207500e+20 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00207500e+20]] 
[[ 1.00207500e+20 0.00000000e+00 -8.00000000e+00 0.00000000e+00] 
[ 0.00000000e+00 1.00207500e+20 0.00000000e+00 0.00000000e+00] 
[ -8.00000000e+00 0.00000000e+00 1.00207500e+20 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00207500e+20]] 
[[ 1.00207500e+20 0.00000000e+00 -8.00000000e+00 0.00000000e+00]

Si está utilizando un sistema de 32 bits (o si tiene una versión de 32 bits de pitón y numpy instalado en un sistema de 64 bits) se encontrará con problemas de precisión y obtendrá una respuesta diferente, como lo señala @Lucero a continuación. Puede intentar especificar explícitamente flotantes de 64 bits en ese caso (aunque las operaciones serán más lentas). Por ejemplo, intente usar np.array(..., dtype=np.float64).

Si usted piensa que necesita precison adicional, puede utilizar np.longdouble (Igual que en un sistema de 64 bits y np.float96 de 32 bits), pero esto no puede ser apoyado en todas las plataformas y muchas funciones de álgebra lineal truncará cosas de vuelta a la precisión nativa.

Además, ¿qué biblioteca BLAS estás utilizando? Los resultados numpy y octava son probablemente los mismos porque están usando la misma biblioteca BLAS.

Por último, se puede simplificar el código numpy a:

import numpy as np 
A = np.array([[1,  -0.015, -0.025, -0.035], 
       [0.015, 1,  0.035, -0.025], 
       [0.025, -0.035, 1,  0.015], 
       [0.035, 0.025, -0.015, 1]]) 
P = np.eye(4)*1e20 
print A.dot(P.dot(A.T)) 
print A.dot(P).dot(A.T)

Fuente

2012-09-08 00:08:19

Realmente no veo el punto de 32 frente a 64 bits (todos los matlab + numpy usan doble precisión por defecto). Sin embargo, la biblioteca de Blas es el punto probable, cuando con ATLAS obtuve el mismo resultado que @ user1655812 sin obtener el exacto. Para arrojar algo más, 'np.einsum' podría hacer el mismo resultado mientras evita ATLAS aquí tal vez. – seberg

La diferencia entre 32 bits y 64 bits es suficiente para aparecer en este caso. En cualquier caso, recibo una respuesta significativamente diferente, pero no es suficiente para explicar por completo los resultados del PO. Estoy de acuerdo, es probable que se deba a la biblioteca BLAS, pero no pensé en probarlo. (¡Me alegra que lo haya hecho!) Mis resultados no están usando ATLAS. (¡Y buen punto sobre einsum!) –

Lo sentimos, pero siempre son puntos flotantes de 64 bits, el sistema no importa. – seberg

np.einsum es mucho más cerca de los MATLAB

In [1299]: print(np.einsum('ij,jk,lk',A,P,A)) 
[[ 1.00207500e+20 0.00000000e+00 -5.07812500e-02 0.00000000e+00] 
[ 0.00000000e+00 1.00207500e+20 5.46875000e-02 0.00000000e+00] 
[ -5.46875000e-02 5.46875000e-02 1.00207500e+20 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00207500e+20]]

Los términos diagonales fuera de la fila y la columna 2 son diferentes, pero tiene los mismos 0 en otros lugares.

Con el punto doble, P.dot(A.T) crea errores de redondeo cuando agrega los productos. Esos se propagan en el próximo dot. einsum genera todos los productos, con solo un resumen. Sospecho que el intérprete de MATLAB reconoce esta situación y realiza un cálculo especial diseñado para minimizar los errores de redondeo.

Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix - es una pregunta más reciente con la misma explicación.

Fuente

2015-08-07 15:47:50 hpaulj

Diferencia numérica Matlab/Octave/Numpy

Respuesta

Cuestiones relacionadas