2010-11-21 24 views
3

私は行列の平方根をとろうとしています。行列Bがありますので、B*B=Aです。私が見つけた方法のどれも実績を上げていません。作業行列平方根

まず私はWikipediaに、この式を見つけました:

設定Y_0 = AZ_0 = Iを、その後の繰り返し:

Y_{k+1} = .5*(Y_k + Z_k^{-1}),

Z_{k+1} = .5*(Z_k + Y_k^{-1}).

その後YBに収束しなければなりません。

しかし、(逆行列のためにnumpyの使用)Pythonでアルゴリズムを実装する

は私にゴミ結果与えた :

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

私は応援行列広場について多くを知らないが、私は上記よりも良好に機能するアルゴリズムが存在する必要があります理解します?

答えて

7

(1)行列[1,2,3,4]の平方根は、その行列の固有値が負であるため、何か複雑なものを与えるべきです。だからあなたの解決策は、まずは正しいとは言えません。

(2)linalg.sqrtmは行列ではなく配列を返します。したがって、それらを掛けるために*を使用することは良い考えではありません。あなたの場合、ソリューションは正しいですが、あなたはそれを見ていません。

編集あなたはそれが正しいのです参照してくださいよ、次のことを試してください。

asmatrix(scipy.linalg.sqrtm(matrix('1,2;3,4')))**2 
+0

ありがとうございます。完璧に動作します。 'sqrtm(matrix( '1,2,3,4')** 2)'がうまくいかない理由はまだ分かりません。 –

+0

これはうまくいきます...解決策を二乗することを試みてください、http://www.mathworks.com/help/techdoc/ref/sqrtmを見てください。htmlの – steabert

+1

右、そこにはたくさんのソリューションがあり、sqrtmは単に私以外のものを選択します。 –

3

あなたのマトリックス[1 2; 3 4]は正ではないので、実際の行列の領域では問題の解決策はありません。

2

あなたがやっている行列の平方根の目的は何ですか?実際のアプリケーションでは、行列が実際に対称正定値(例えば共分散)になる可能性があるので、複素数に遭遇すべきではありません。あなたはここを参照してください、スケールLU分解と同様に、コレスキー分解を計算することができ、その場合には

http://en.wikipedia.org/wiki/Cholesky_decomposition

別の実用的な例としては、あなたの行列が回転している場合は、あなたが最初の行列ログに分解することができるだけで割るです2をログ空間に入れてから、行列指数を使ってローテーションに戻ってください。「ジェネリック行列平方根」を求めるのは奇妙に思えますが、特定のアプリケーションをより深く理解したいと思うかもしれません。