R言語には逆行列を解く関数がある
連立方程式のときも使ったsolve関数を使うと、逆行列を求めることができる
例えば \[ \begin{equation} A= \begin{pmatrix} 2 &4 \\ 1 &3 \end{pmatrix} \Longrightarrow A^{-1}= \begin{pmatrix} 1.5 &-2 \\ -0.5 &1 \end{pmatrix} \end{equation} \]
すげー楽
例えば \[ \begin{equation} A= \begin{pmatrix} 2 &4 \\ 1 &3 \end{pmatrix} \Longrightarrow A^{-1}= \begin{pmatrix} 1.5 &-2 \\ -0.5 &1 \end{pmatrix} \end{equation} \]
> # 行列を生成
> A <- matrix(c(2, 1, 4, 3),2,2)
> A # 確認
[,1] [,2]
[1,] 2 4
[2,] 1 3
>
> #逆行列を求める組み込み関数
> solve(A)
[,1] [,2]
[1,] 1.5 -2
[2,] -0.5 1
すげー楽
solve関数で一撃計算できない行列
当然、正則でない行列は逆行列が存在しないため、計算できないが
以下のパターンもエラーになってしまう
例えば \[ \begin{equation} A= \begin{pmatrix} 10^{100} &4 \\ 1 &3 \end{pmatrix} \end{equation} \] を考える。
行列式は$3 \times 10^{100}-4$となり、正則行列となるので、逆行列をもつ行列だが、solve関数で求めるとエラーとなる
原因はおそらく、「.Machine＄double.ep」の値
組み込み関数であるため、どういうロジックで計算しているのか不明だが、ヘルプを見る限り、「.Machine＄double.ep」の値より、小さい数字に関しては計算しないようになっているっぽい
そこで、第二引数として、toleranceを設定することで、計算下限値を調整すればなんとかなる
これで対応できた。
しかし、この対応法に気が付かず、自作で逆行列を求める関数を作ってしまった...
逆行列を解く関数を自作する
逆行列と言えば、
掃き出し法で$(A|E)$の形から、基本変形を駆使して、$(E|B)$に持っていくと$B$は$A$の逆行列になっている！
ってのを使うね
他の言語なら、間違いなく掃き出し法だが、
R言語だと、余因子行列が一瞬で生成できるから、余因子行列から生成する方法をとる
問題なさそう
でもはやり、組み込み関数に比べたら処理速度落ちるだろうから、toleranceを設定する方法とるなぁ
以下のパターンもエラーになってしまう
例えば \[ \begin{equation} A= \begin{pmatrix} 10^{100} &4 \\ 1 &3 \end{pmatrix} \end{equation} \] を考える。
行列式は$3 \times 10^{100}-4$となり、正則行列となるので、逆行列をもつ行列だが、solve関数で求めるとエラーとなる
> # 係数行列を生成
> A <- matrix(c(10^100,1,4,3),2,2)
> A # 確認
[,1] [,2]
[1,] 1e+100 4
[2,] 1e+00 3
>
> #逆行列を求める組み込み関数
> solve(A)
Error in solve.default(A) :
system is computationally singular: reciprocal condition number = 3e-100
原因はおそらく、「.Machine＄double.ep」の値
組み込み関数であるため、どういうロジックで計算しているのか不明だが、ヘルプを見る限り、「.Machine＄double.ep」の値より、小さい数字に関しては計算しないようになっているっぽい
そこで、第二引数として、toleranceを設定することで、計算下限値を調整すればなんとかなる
> #逆行列を求める組み込み関数にtoleranceを設定
> solve(A, tol = 10^(-100))
[,1] [,2]
[1,] 1.000000e-100 -1.333333e-100
[2,] -3.333333e-101 3.333333e-01
これで対応できた。
しかし、この対応法に気が付かず、自作で逆行列を求める関数を作ってしまった...
逆行列を解く関数を自作する
逆行列と言えば、
掃き出し法で$(A|E)$の形から、基本変形を駆使して、$(E|B)$に持っていくと$B$は$A$の逆行列になっている！
ってのを使うね
他の言語なら、間違いなく掃き出し法だが、
R言語だと、余因子行列が一瞬で生成できるから、余因子行列から生成する方法をとる
my_solve <- function(A) {
# 入力データチェック
len <- length(A[,1]) # 行列の大きさ
if(length(A[1,]) != len) return (print("正方行列じゃないよ"))
det_A <- det(A) # 行列式
if(det_A == 0) return (print("正則行列じゃないよ"))
# 2×2行列の場合、余因子行列の行列式が求められないため、ベタ打ち
if(len == 2) return (matrix(1/det_A * c(A[2, 2], -A[1, 2], -A[2, 1], A[1, 1]), 2, 2, T))
# 出力用逆行列生成
A_ <- matrix(0, len, len)
# 余因子行列作成
for(i in 1 : len) {
for(j in 1 : len) {
A_[j, i] <- ((-1)^(i + j)) * det(A[-i, -j])
}
}
# 全体を行列式で割れば完成！
return (A_ / det_A)
}
検算してみる
> my_solve(A)
[,1] [,2]
[1,] 1.000000e-100 -1.333333e-100
[2,] -3.333333e-101 3.333333e-01
# いい感じの結果
# 他にも検算してみる
# 逆行列は元の行列と掛けて単位行列になるか見ればいいので
> for(n in 2:10) {
+ X <- matrix(runif(n * n), n, n)
+ print(round(X %*% my_solve(X),10))
+ }
[,1] [,2]
[1,] 1 0
[2,] 0 1
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
[2,] 0 1 0 0 0
[3,] 0 0 1 0 0
[4,] 0 0 0 1 0
[5,] 0 0 0 0 1
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 0 0 1 0 0
[5,] 0 0 0 0 1 0
[6,] 0 0 0 0 0 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 1 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0
[3,] 0 0 1 0 0 0 0
[4,] 0 0 0 1 0 0 0
[5,] 0 0 0 0 1 0 0
[6,] 0 0 0 0 0 1 0
[7,] 0 0 0 0 0 0 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0
[5,] 0 0 0 0 1 0 0 0
[6,] 0 0 0 0 0 1 0 0
[7,] 0 0 0 0 0 0 1 0
[8,] 0 0 0 0 0 0 0 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0
[7,] 0 0 0 0 0 0 1 0 0
[8,] 0 0 0 0 0 0 0 1 0
[9,] 0 0 0 0 0 0 0 0 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0
[9,] 0 0 0 0 0 0 0 0 1 0
[10,] 0 0 0 0 0 0 0 0 0 1
問題なさそう
でもはやり、組み込み関数に比べたら処理速度落ちるだろうから、toleranceを設定する方法とるなぁ
0 件のコメント:
コメントを投稿