Matlab/Octave/Numpy 숫자 차이

Question

Matlab에서 작동하는 알고리즘을 numpy로 포팅하고 이상한 동작을 관찰했습니다. 코드의 관련 부분은 내가 matlab에 함께 실행할 때,

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' 

이 코드는 다음과 같은 행렬을 제공 : 예상대로 V1 및 V2가 동일한 지

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 

참고. 동일한 코드 옥타브에서 실행될 때

, 그것이 제공 NumPy와에서는

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 

는 세그먼트 그래서

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

출력

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

되고, 옥타브 NumPy와 두 동일한 답변을 제공하지만 Matlab과 매우 다릅니다. 첫 번째 요점은 V1! = V2이며, 이는 옳지 않은 것입니다. 다른 점은 대각선이 아닌 요소가 대각선의 요소보다 훨씬 더 큰 차수이지만 이것이 알고리즘에서 문제를 일으키는 것 같습니다.

누구나 numpy와 Octave가이 방식으로 작동하는 것을 알고 있습니까?

Answer 1

두 배는 내부적으로 사용되며, 정확도는 약 15 자리입니다. 수학 연산이이를 초과 할 가능성이있어 수학적으로 잘못된 결과가 발생합니다. 내가 NumPy와에 사용할 수있는 높은 정밀도가 없음을 수집 docs에서 : http://floating-point-gui.de/

편집 : 읽기

가치. SymPy (may give you the needed precision) - 해당 라이브러리가 잘 작동하는 것 같습니다. 그것은 가치가 무엇이든 들어

는

Answer 2

, 나는 64 비트 시스템에 대한 MATLAB과 동일한 결과를 얻을 수 :

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

당신은 32 비트 시스템을 사용하고 (또는 당신은 파이썬의 32 비트 버전이있는 경우 numpy가 64 비트 시스템에 설치되어있는 경우) 아래의 @Lucero가 지적한 바와 같이 정밀 문제가 발생하여 다른 대답을 얻게됩니다. 이 경우 64 비트 부동 소수점을 명시 적으로 지정할 수도 있습니다 (그러나 작업은 느려질 것입니다). 예를 들어 np.array(..., dtype=np.float64)을 사용해보세요. 당신이 추가 precison 필요하다고 생각하는 경우

, 당신은 (32 비트에서 64 비트 시스템 및 np.float96에 과 동일) np.longdouble를 사용할 수 있지만이 모든 플랫폼 및 자릅니다 많은 선형 대수 함수에 지원되지 않을 수 있습니다 네이티브 정밀도로 되돌아갑니다.

또한 BLAS 라이브러리는 무엇을 사용하고 있습니까? numpy와 octave 결과는 동일한 BLAS 라이브러리를 사용하기 때문에 동일 할 수 있습니다.

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은 행과 열이가 다른의 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]] 

오프 대각선 조건에 많은 가까이 있지만 :

마지막으로, 당신은 아래로 NumPy와 코드를 단순화 할 수 있습니다 다른 곳에서는 같은 0을가집니다.

두 점으로 표시하면 P.dot(A.T)은 제품을 추가 할 때 반올림 오류를 생성합니다. 그것들은 다음 dot으로 전파됩니다. einsum은 합계를 하나만 사용하여 모든 제품을 생성합니다. 필자는 MATLAB 인터프리터가 이러한 상황을 인식하고 반올림 오류를 최소화하도록 특수 계산을 수행합니다.

Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix - 동일한 설명이있는 최근 질문입니다.