6

Matlab で動作するアルゴリズムを 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'

このコードを Matlab で実行すると、次の行列が提供されます。

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

予想どおり、V1 と V2 は同じであることに注意してください。

同じコードを Octave で実行すると、以下が提供されます。

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

numpy では、セグメントは次のようになります

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

出力する

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

したがって、Octave と numpy の両方で同じ答えが得られますが、Matlab のものとは大きく異なります。最初のポイントは、V1 != V2 であり、正しくないようです。もう 1 つのポイントは、非対角要素は対角要素よりも桁違いに小さいにもかかわらず、これが私のアルゴリズムで何らかの問題を引き起こしているように思われることです。

numpy と Octave がこのように動作する方法を知っている人はいますか?

4

3 に答える 3

6

内部的にはdoubleを使用しており、精度は約15桁です。数学演算がこれを超える可能性があり、数学的に間違った結果が発生します。

読む価値があります：http：//floating-point-gui.de/

編集：ドキュメントから、Numpyで利用できるより高い精度はないことがわかりました。ただし、SymPyを使用すると、必要な精度が得られる可能性があります。そのライブラリが機能する場合も同様です。

于 2012-09-07T22:16:26.593 に答える
4

どんなに価値があっても、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 ビット バージョンの python と numpy が 64 ビット システムにインストールされている場合)、精度の問題が発生し、@Lucero で指摘されているように、別の答えが得られます。下。その場合は、64 ビット浮動小数点を明示的に指定してみてください (ただし、操作は遅くなります)。たとえば、 を使用してみてくださいnp.array(..., dtype=np.float64)。

追加の精度が必要だと思われる場合は、 ( 64 ビット システムおよび32 ビット システムとnp.longdouble同じ) を使用できますが、これはすべてのプラットフォームでサポートされているわけではなく、多くの線形代数関数はネイティブの精度に切り捨てます。np.float128np.float96

また、どのBLASライブラリを使用していますか? numpy と octave の結果は、同じ BLAS ライブラリを使用しているため、おそらく同じです。

最後に、numpy コードを次のように単純化できます。

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)
于 2012-09-08T00:08:19.977 に答える
0

np.einsumMATLAB に非常に近い

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

行 2 と列 2 の非対角項は異なりますが、他の場所では同じ 0 です。

二重ドットを使用P.dot(A.T)すると、製品を追加するときに丸め誤差が生じます。それらは次の に伝播しdotます。 einsum1 回の合計で、すべての積を生成します。MATLAB インタープリターがこの状況を認識し、丸め誤差を最小限に抑えるように設計された特別な計算を実行すると思われます。

Numpy Matrix Multiplication U*B*UT Results in Nonsymmetric Matrix - 同じ説明の最近の質問です。

于 2015-08-07T15:47:50.267 に答える