0

-1/2 の累乗で行列Aを計算する必要があります。これは基本的に、初期行列の逆行列の平方根を意味します。

A が特異値の場合、Moore-Penrose の一般化逆行列は MASS パッケージのginv関数を使用して計算されます。それ以外の場合、通常の逆行列はsolve関数を使用して計算されます。

行列 A は次のように定義されます。

A <- structure(c(604135780529.807, 0, 58508487574887.2, 67671936726183.9, 
            0, 0, 0, 1, 0, 0, 0, 0, 58508487574887.2, 0, 10663900590720128, 
            10874631465443760, 0, 0, 67671936726183.9, 0, 10874631465443760, 
            11315986615387788, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), .Dim = c(6L, 
                                                                                   6L))

ランクと次元の比較で特異点をチェックします。

rankMatrix(A) == nrow(A)

上記のコードは FALSE を返すため、逆を取得するにはginvを使用する必要があります。A の逆関数は次のとおりです。

A_inv <- ginv(A)

逆行列の平方根は、expm パッケージの sqrtm 関数で計算されます。

library(expm)
sqrtm(A_inv)

関数は次のエラーを返します。

solve.default(X[ii, ii] + X[ij, ij], S[ii, ij] - sumU) のエラー:
Lapack ルーチン zgesv: システムは厳密に特異です

では、この場合の平方根はどのように計算すればよいでしょうか? 行列 A は常に特異であるとは限らないため、問題の一般的な解を提供する必要があることに注意してください。

4

1 に答える 1

6

あなたの質問は、2 つの異なる問題に関連しています。

  1. 逆行列
  2. 行列の平方根

逆

逆行列は特異行列には存在しません。一部のアプリケーションでは、Moore-Penrose またはその他の一般化された逆関数が、逆関数の適切な代替として使用される場合があります。ただし、ほとんどの場合、コンピュータ数値には丸め誤差が発生することに注意してください。これらのエラーにより、特異な行列がコンピューターに規則的に表示される場合や、その逆の場合があります。

あなたが与えたマトリックスのブロック構造を常に示す場合A、その非対角ブロックのみを考慮することをお勧めします

A3 = A[ c( 1, 3, 4 ), c( 1, 3, 4 ) ]

A3
             [,1]         [,2]         [,3]
[1,] 6.041358e+11 5.850849e+13 6.767194e+13
[2,] 5.850849e+13 1.066390e+16 1.087463e+16
[3,] 6.767194e+13 1.087463e+16 1.131599e+16

すべての代わりに、A効率を高め、丸めの問題を減らします。残りの 1 対角のエントリは、平方根の逆数で 1 のままであるため、それらで計算を混乱させる必要はありません。この単純化の影響を理解するには、R が計算できることに注意してください。

A3inv = solve(A3)

計算できなかったのに

Ainv = solve(A)

しかし、以下で明らかになるように、A3inverse は必要ありません。

平方根

原則として、行列の平方根は、行列が対角ヨルダン正規形 ( https://en.wikipedia.org/wiki/Jordan_normal_formA )を持つ場合にのみ存在します。したがって、必要な問題の真に一般的な解決策はありません。

幸いなことに、「ほとんどの」(実数または複素数) 行列が可逆であるように、「ほとんどの」(実数または複素数) 行列は対角複素数ジョルダン正規形を持ちます。この場合、ジョーダン標準形

A3 = T・J・T⁻¹</p>

R で次のように計算できます。

X = eigen(A3)
T = X$vectors
J = Diagonal( x=X$values )

このレシピをテストするには、比較してください

Tinv = solve(T)
T %*% J %*% Tinv

とA3。A3対角ジョーダン正規形の場合、それらは (丸め誤差まで) 一致する必要があります。

は対角なのでJ、その平方根は単純に平方根の対角行列です。

Jsqrt = Diagonal( x=sqrt( X$values ) )

そのようにJsqrt·Jsqrt = J。さらに、これは

(T・Jsqrt・T⁻¹)² = T・Jsqrt・T⁻¹・T・Jsqrt・T⁻¹ = T・Jsqrt・Jsqrt・T⁻¹ = T・J・T⁻¹ = A3

実際に得られるように

√A3 = T·Jsqrt·T⁻¹</p>

またはRコードで

A3sqrt = T %*% Jsqrt %*% Tinv

これをテストするには、計算します。

A3sqrt %*% A3sqrt

と比較してくださいA3。

逆数の平方根

逆数の平方根 (または、同様に、平方根の逆数) は、対角ジョルダン正規形が計算されると簡単に計算できます。J使用する代わりに

Jinvsqrt = Diagonal( x=1/sqrt( X$values ) )

上記と同様に、

A3invsqrt = T %*% Jinvsqrt %*% Tinv

そして観察する

A3·A3invsqrt² = … = T·(J/√J/√J)·T⁻¹ = 1

A3invsqrt が望ましい結果になるように、単位行列を計算します。

A3 が可逆でない場合、一般化された逆 (必ずしもムーア-ペンローズのものではない) は、未定義のすべてのエントリをJinvsqrt0 で置き換えることによって計算できますが、上で述べたように、これは全体的な観点から適切な注意を払って行う必要があります。アプリケーションと丸め誤差に対するその安定性。

A3 が対角ジョーダン標準形を持たない場合、平方根がないため、上記の式は別の結果をもたらします。不運なときにこのケースに遭遇しないようにするには、次のテストを実装するのが最善です

A3invsqrt %*% A3 %*% A3invsqrt

は、1 行列と見なすものに十分近いです (これは、最初に A3 が可逆であった場合にのみ適用されます)。

追伸: の対角線のエントリごとに記号 ± を前に付けることができることに注意してJinvsqrtください。

于 2015-08-10T11:44:59.640 に答える