Matlab/Ottava/Numpia differenza numerica

Question

Sto eseguendo il porting di un algoritmo che funziona in Matlab per numpy e ho osservato uno strano comportamento. Il segmento rilevante di codice è

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' 

Questo codice, quando corro con Matlab, fornisce le seguenti matrici:

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 

noti che V1 e V2 sono gli stessi, come previsto.

Quando lo stesso codice viene eseguito in Octave, fornisce:

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 

Nel numpy, il segmento diventa

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

che emette

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

Quindi, sia Octave e numpy fornisce la stessa risposta, ma è molto diversa da quella di Matlab. Il primo punto è V1! = V2, che non sembra giusto. L'altro punto è che, sebbene gli elementi non diagonali siano molti ordini di grandezza più piccoli di quelli diagonali, questo sembra causare alcuni problemi nel mio algoritmo.

Qualcuno conosce il modo insensato e Octave si comporta in questo modo?

Answer 1

Utilizzano i doppi internamente, con una precisione di circa 15 cifre. Le tue operazioni matematiche probabilmente superano questo, il che causa risultati matematicamente sbagliati.

Merita una lettura: http://floating-point-gui.de/

Edit: Dal docs ho capito che non ci sia una maggiore precisione disponibile per Numpy. Sembra che SymPy sia may give you the needed precision - se quella libreria funziona anche per te.

Answer 2

Per quel che vale, io ottengono risultati identici a MATLAB su un sistema a 64 bit:

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

Se si utilizza un sistema a 32-bit (o si dispone di una versione a 32 bit di pitone e numpy installato su un sistema a 64 bit) si incontrano problemi di precisione e si ottiene una risposta diversa, come indicato da @Lucero di seguito. Potresti provare esplicitamente a specificare i float a 64 bit in quel caso (le operazioni saranno più lente, però). Ad esempio, prova a utilizzare np.array(..., dtype=np.float64).

Se pensi di aver bisogno di ulteriori precison, è possibile utilizzare np.longdouble (Idem come su un sistema a 64 bit e np.float96 a 32 bit), ma questo potrebbe non essere supportato su tutte le piattaforme e le molte funzioni di algebra lineare troncherà le cose tornano alla precisione nativa.

Inoltre, quale libreria BLAS stai usando? I risultati di numpy e di ottava sono probabilmente gli stessi perché usano la stessa libreria BLAS.

Infine, è possibile semplificare il codice NumPy verso il basso per:

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 è molto più vicino al 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]] 

I termini diagonale fuori di riga e colonna 2 sono diversi, ma ha gli stessi 0 altrove.

Con il doppio punto, P.dot(A.T) crea errori di arrotondamento quando si aggiungono i prodotti. Quelli si propagano nel prossimo dot. einsum genera tutti i prodotti, con una sola sommatoria. Sospetto che l'interprete MATLAB riconosca questa situazione ed esegua un calcolo speciale progettato per minimizzare gli errori di arrotondamento.

Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix - è una domanda più recente con la stessa spiegazione.