2

乱数の行列を用意し、その逆数を計算し、それを元の行列で行列乗算します。これにより、理論的には単位行列が得られます。どうすればnumpyそれを私に任せることができますか?

import numpy

A = numpy.zeros((100,100))
E = numpy.zeros((100,100))

size = 100

for i in range(size):
    for j in range(size):
        A[i][j]+=numpy.random.randint(10)
        if i == j:
            E[i][j]+=1

A_inv = numpy.linalg.linalg.inv(A)
print numpy.dot(A, A_inv)

コードを実行すると生成されます

[me]machine @ numeric $ python rand_diag.py 
[[  1.00000000e+00  -7.99360578e-15  -1.14491749e-16 ...,   3.81639165e-17
   -4.42701431e-15   1.17961196e-15]
[ -5.55111512e-16   1.00000000e+00  -2.22044605e-16 ...,  -3.88578059e-16
    1.33226763e-15  -8.32667268e-16]

結果が単位行列であることは明らかですが、正確ではないため、print numpy.dot(A, A_inv) == E明らかにFalse. 私は線形代数を練習し、マシンが限界に達する行列のサイズを見つけようとしてこれを行っています。を取得するTrueことは、教訓的に魅力的です。

編集：

設定size=10000、メモリ不足

[me]machine @ numeric $ Python(794) malloc:
***mmap(size=800002048) failed (error code=12)
*** error: can\'t allocate region
*** set a breakpoint in malloc_error_break to debug
Traceback (most recent call last):
File "rand_diag.py", line 14, in <module>     A_inv = numpy.linalg.linalg.inv(A)
File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 445, in inv
return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 323, in solve
a, b = _fastCopyAndTranspose(t, a, b)
File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 143, in _fastCopyAndTranspose
cast_arrays = cast_arrays + (_fastCT(a),)
MemoryError

[1]+  Exit 1                  python rand_diag.py

より多くのメモリを割り当てるにはどうすればよいですか? また、これを並行して実行するにはどうすればよいですか (私は 4 つのコアを持っています)。

4

4 に答える 4

8

得Trueは教訓的に魅力的ですが、浮動小数点計算の現実から切り離されることにもなります。

浮動小数点を扱うときは、不正確な結果だけでなく、発生するあらゆる種類の数値の問題に備える必要があります。

What Every Computer Scientist Should Know About Floating-Point Arithmeticを読むことを強くお勧めします。

あなたの特定のケースでは、それA * inv(A)が単位行列に十分近いことを確認するために、行列ノルムを計算し、numpy.dot(A, A_inv) - Eそれが十分に小さいことを確認できます。

補足として、ループを使用してAandを設定する必要はありませんE。代わりに、単に使用できます

A = numpy.random.randint(0, 10, (size,size))
E = numpy.eye(size)
于 2013-03-29T16:28:02.480 に答える
4

すでに指摘されている点のほとんどに同意します。ただし、個々の非対角要素を見るのではなく、それらの実効値の合計を取ることをお勧めします。これはある意味で、不完全な計算の結果として非対角項に漏れた「エネルギー」を反映しています。次に、この RMS 数を対角項の合計で割ると、逆関数がどれだけうまく機能したかを示すメトリックが得られます。たとえば、次のコード:

import numpy
import matplotlib.pyplot as plt
from numpy import mean, sqrt
N = 1000
R = numpy.zeros(N)

for size in range(50,N,50):

  A = numpy.zeros((size, size))
  E = numpy.zeros((size, size))

  for i in range(size):
      for j in range(size):
          A[i][j]+=numpy.random.randint(10)
          if i == j:
              E[i][j]=1

  A_inv = numpy.linalg.linalg.inv(A)
  D = numpy.dot(A, A_inv) - E
  S = sqrt(mean(D**2))
  R[size] = S/size
  print "size: ", size, "; rms is ",  S/size

plt.plot(range(50,N,50), R[range(50, N, 50)])
plt.ylabel('RMS fraction')
plt.show()

rms エラーは、アレイのサイズが 950x950 のサイズまでかなり安定していることを示しています (少し遅くなります...)。ただし、それは決して「正確」ではなく、いくつかの外れ値があります (おそらく行列がほぼ特異である場合 - これはランダム行列で発生する可能性があります)。

プロットの例 (実行するたびに、少し異なって見えます):

ここに画像の説明を入力

于 2013-03-29T16:55:47.827 に答える
2

最後に、答えを丸めることができます

m = np.round(m, decimals=10)

または、それらが大きく異なるかどうかを確認します。

np.abs(A*A.I - i).mean() < 1e-10

小さな数字を消したい場合。

これをnumpy.matrixクラスで実装します。

import numpy

size = 100
A = numpy.matrix(numpy.random.randint(0,10,(size,)*2))
E = numpy.eye(size)

print A * A.I
print np.abs(A * A.I - E).mean() < 1e-10
于 2013-03-29T16:27:12.470 に答える
0

あなたの問題は、一般的なフロート比較の問題に減らすことができます。このような配列を比較する正しい方法は次のとおりです。

EPS = 1e-8  # for example
(np.abs(numpy.dot(A, A_inv) - E) < EPS).all()
于 2013-03-29T16:29:02.600 に答える