10

R : で非正方線形システムを解く方法はA X = B?

(系に解がないか、解が無限にある場合)

例 :

A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
B=matrix(c(-17,28,11),3,1,T)

A
     [,1] [,2] [,3] [,4]
[1,]    0    1   -2    3
[2,]    5   -3    1   -2
[3,]    5   -2   -1    1


B
     [,1]
[1,]  -17
[2,]   28
[3,]   11
4

2 に答える 2

11

行列 A の行数が列数よりも多い場合、最小二乗法を使用する必要があります。

行列 A の行数が列数より少ない場合、特異値分解を実行する必要があります。各アルゴリズムは、仮定を使用してソリューションを提供するために最善を尽くします。

SVD をソルバーとして使用する方法を示すリンクを次に示します。

http://www.ecse.rpi.edu/~qji/CV/svd_review.pdf

それをあなたの問題に適用して、うまくいくか見てみましょう:

入力行列Aと既知の RHS ベクトルB:

> A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
> B=matrix(c(-17,28,11),3,1,T)
> A
     [,1] [,2] [,3] [,4]
[1,]    0    1   -2    3
[2,]    5   -3    1   -2
[3,]    5   -2   -1    1
> B
     [,1]
[1,]  -17
[2,]   28
[3,]   11

A行列を分解しましょう:

> asvd = svd(A)
> asvd
$d
[1] 8.007081e+00 4.459446e+00 4.022656e-16

$u
           [,1]       [,2]       [,3]
[1,] -0.1295469 -0.8061540  0.5773503
[2,]  0.7629233  0.2908861  0.5773503
[3,]  0.6333764 -0.5152679 -0.5773503

$v
            [,1]       [,2]       [,3]
[1,]  0.87191556 -0.2515803 -0.1764323
[2,] -0.46022634 -0.1453716 -0.4694190
[3,]  0.04853711  0.5423235  0.6394484
[4,] -0.15999723 -0.7883272  0.5827720

> adiag = diag(1/asvd$d)
> adiag
          [,1]      [,2]        [,3]
[1,] 0.1248895 0.0000000 0.00000e+00
[2,] 0.0000000 0.2242431 0.00000e+00
[3,] 0.0000000 0.0000000 2.48592e+15

鍵は次のとおりです。 の 3 番目の固有値dは非常に小さいです。逆に、 の対角要素adiagは非常に大きいです。解く前に、ゼロに等しく設定します。

> adiag[3,3] = 0
> adiag
          [,1]      [,2] [,3]
[1,] 0.1248895 0.0000000    0
[2,] 0.0000000 0.2242431    0
[3,] 0.0000000 0.0000000    0

それでは、解を計算してみましょう (上記のリンクのスライド 16 を参照してください)。

> solution = asvd$v %*% adiag %*% t(asvd$u) %*% B
> solution
          [,1]
[1,]  2.411765
[2,] -2.282353
[3,]  2.152941
[4,] -3.470588

解決策が得られたので、それを元に戻して、同じ結果が得られるかどうかを確認しましょうB

> check = A %*% solution
> check
     [,1]
[1,]  -17
[2,]   28
[3,]   11

それはBあなたが始めた側なので、私たちは良いと思います。

これは、AMS からのもう 1 つの素晴らしい SVD ディスカッションです。

http://www.ams.org/samplings/feature-column/fcarc-svd

于 2013-11-04T12:13:14.680 に答える