Matlab/Octave/Numpy différence numérique

Question

Je suis portage d'un algorithme qui fonctionne dans Matlab à numpy et observé un comportement étrange. Le segment concerné du code est

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' 

Ce code, quand je lance avec Matlab, fournit les matrices suivantes:

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 

Notez que V1 et V2 sont les mêmes, comme prévu.

Lorsque le même code fonctionne dans Octave, il fournit:

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, le segment devient

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

qui sort

[[ 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]] 

Ainsi, les deux Octave et numpy fournit la même réponse, mais c'est très différent de celui de Matlab. Le premier point est que V1! = V2, ce qui ne semble pas correct. L'autre point est que, bien que les éléments non-diagonaux soient beaucoup plus petits que les diagonaux, cela semble causer un problème dans mon algorithme.

Est-ce que quelqu'un sait très peu et Octave se comporte de cette façon?

Answer 1

Ils utilisent le double en interne, qui n'ont que 15 chiffres de précision. Vos opérations mathématiques dépassent probablement ceci, ce qui provoque des résultats mathématiquement faux.

Vaut lire: http://floating-point-gui.de/

Edit: De l'docs je suppose qu'il n'y a pas une plus grande précision disponible pour Numpy. Il semble que SymPy si may give you the needed precision - si cette bibliothèque fonctionne pour vous aussi.

Answer 2

Pour ce que ça vaut, j'obtenir des résultats identiques à Matlab sur un système 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 vous utilisez un système 32 bits (ou si vous avez une version 32 bits de python et numpy installé sur un système 64 bits) vous rencontrerez des problèmes de précision, et obtiendrez une réponse différente, comme indiqué par @Lucero ci-dessous. Vous pourriez essayer de spécifier explicitement les flottants de 64 bits dans ce cas (les opérations seront cependant plus lentes). Par exemple, essayez d'utiliser np.array(..., dtype=np.float64).

Si vous pensez que vous avez besoin precison supplémentaire, vous pouvez utiliser np.longdouble (Comme sur un système 64 bits et np.float96 sur 32 bits), mais cela ne peut être pris en charge sur toutes les plates-formes et de nombreuses fonctions d'algèbre linéaire tronquent les choses à la précision native.

En outre, quelle bibliothèque BLAS utilisez-vous? Les résultats numpy et octave sont probablement les mêmes car ils utilisent la même bibliothèque BLAS.

Enfin, vous pouvez simplifier votre code numpy jusqu'à:

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) 

Answer 3

np.einsum est beaucoup plus proche de 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]] 

Les termes hors diagonale en ligne et la colonne 2 sont différents, mais il a les mêmes 0 ailleurs.

Avec le double point, P.dot(A.T) crée des erreurs d'arrondi lors de l'ajout des produits. Ceux se propagent dans le prochain dot. einsum génère tous les produits, avec une seule sommation. Je soupçonne que l'interpréteur MATLAB reconnaît cette situation, et effectue un calcul spécial conçu pour minimiser les erreurs d'arrondi.

Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix - est une question plus récente avec la même explication.