私は行列の平方根をとろうとしています。行列B
がありますので、B*B=A
です。私が見つけた方法のどれも実績を上げていません。作業行列平方根
まず私はWikipediaに、この式を見つけました:
設定Y_0 = A
とZ_0 = I
を、その後の繰り返し:
Y_{k+1} = .5*(Y_k + Z_k^{-1}),
Z_{k+1} = .5*(Z_k + Y_k^{-1}).
その後Y
がB
に収束しなければなりません。
は私にゴミ結果与えた :
>>> def denbev(Y,Z,n):
if n == 0: return Y,Z
return denbev(.5*(Y+Z**-1), .5*(Z+Y**-1), n-1)
>>> denbev(matrix('1,2;3,4'), matrix('1,0;0,1'), 3)[0]**2
matrix([[ 1.31969074, 1.85986159],
[ 2.78979239, 4.10948313]])
>>> denbev(matrix('1,2;3,4'), matrix('1,0;0,1'), 100)[0]**2
matrix([[ 1.44409972, 1.79685675],
[ 2.69528512, 4.13938485]])
あなたが見ることができるように、100回の反復は、3回の反復よりも悪化し結果を与え、いずれの結果も40%の誤差マージンに達しません。
それから私はscipyのダウンロードsqrtm方法を試してみましたが、それがさらに悪化した:
>>> scipy.linalg.sqrtm(matrix('1,2;3,4'))**2
array([[ 0.09090909+0.51425948j, 0.60606061-0.34283965j],
[ 1.36363636-0.77138922j, 3.09090909+0.51425948j]])
>>> scipy.linalg.sqrtm(matrix('1,2;3,4')**2)
array([[ 1.56669890+0.j, 1.74077656+0.j],
[ 2.61116484+0.j, 4.17786374+0.j]])
私は応援行列広場について多くを知らないが、私は上記よりも良好に機能するアルゴリズムが存在する必要があります理解します?
ありがとうございます。完璧に動作します。 'sqrtm(matrix( '1,2,3,4')** 2)'がうまくいかない理由はまだ分かりません。 –
これはうまくいきます...解決策を二乗することを試みてください、http://www.mathworks.com/help/techdoc/ref/sqrtmを見てください。htmlの – steabert
右、そこにはたくさんのソリューションがあり、sqrtmは単に私以外のものを選択します。 –