Rで普通の数値解析その３
あっ, どーも僕です.
そういえば, 2019年にハンドボール女子世界選手権が日本開催に決定したのでここで宣伝しときますね！！
日本代表女子の試合はyoutubeに動画が上がっています！！
ハンドボールをもっと楽しみたい人は本場のハンドボールを見てください！！
full match のところにあるパリサンジェルマン対バルセロナの試合がオススメです！
EHFTV | LAOLA1.tv
Rで普通の数値解析その３
追記
スクリプトの最新版はこちらから。
Rで普通の数値解析
それはおいといて, 今回の修行本はこちら.
- 作者: 名取亮
- 出版社/メーカー: 朝倉書店
- 発売日: 1993/05
- メディア: 単行本
- この商品を含むブログを見る
こちらから、次の３つを実装してみまいした.
これの中で, 特にハウスホルダー変換＋二分法からのシフト逆反復による固有値と固有ベクトルが, 数値解析の小技集的な存在なので, 線形計算をはじめるときはこれを目標にいろいろ実装してみると楽しいかと思います.
また, いままでの普通のスクリプトをまとめてgoogle driveにアップしました. まとめて落とせるのでよかったら見てみてください.（ファイル名は適当ですが・・・）
MTA法とかも入っています。まだ, バグがあるかもなんで（今日もハウスホルダー変換を直しました)もしありましたら報告してもらえると助かります.
まとめて公開と思ったのですが, どうもプロパティが削除できなかったのでまた今度にしますね。
スクリプト
まず, ハウスホルダー変換＋二分法からのシフト逆反復による固有値と固有ベクトルを載せておきますね.
これを計算するのに, ハウスホルダー変換, 二分法, LU分解, 前進代入, 後退代入, ピポットの部分選択, 逆反復法といままで作ってきたスクリプトがちょこちょこ生かされています.
irisの分散共分散行列で動作確認しています. とりあえずRのeigenと同じ結果が得られることがわかります.
ハウスホルダー変換＋二分法からのシフト逆反復で固有値と固有ベクトル(修正あり)
# ハウスホルダー法 + 二分法 で固有値を求めたのちに # 元の行列を固有値でシフトさせて, 逆反復法で最小固有ベクトルをともめる # シフト付き逆反復法という名前らしい.... # Aは実対称行列 myEigen.power.final <- function (A, lower = 0, upper = NULL, epsilon = 1e-10) { if (!exists ("myEigen.Dichotomy")) { source ("myEigen_Dichotomy.r") } if (!exists ("myPartPivotSelect")) { source ("myPartPivotSelect.r") } if (!exists("myLUdecomposition")) { source ("myLUdecomposition.r") } # 固有値 eigenValues <- myEigen.Dichotomy (myEigen.Householder (A), lower = lower, upper = upper, epsilon = epsilon)[[1]] num <- length (eigenValues) eigenVectors <- matrix (0, num, num) D <- diag (num) # 固有値ごとに固有ベクトルを逆反復法で計算する for (i in 1:num) { # 固有値でシフトして行置換 tmp1 <- myPartPivotSelect (A - eigenValues [i] * D) A.shifted <- tmp1$mat P <- tmp1$P # LU分解 tmp2 <- myLUdecomposition (A.shifted, LDU = FALSE) U <- tmp2$U L <- tmp2$L #loop rtn <- rnorm (num) for (j in 1:5) { # ハウスホルダー変換＋二分法の精度がいいので # ここは3回くらい計算すればいいらしい y <- myForwardSolve (cbind (L, rtn)) v <- myBackSolve (cbind (U, y)) u <- v / max (abs (v)) rtn <- crossprod (t (P), u) } # 固有ベクトルはユークリッドノルムで正規化 eigenVectors[, i] <- u / sqrt (c (crossprod (u))) } return (list (values = eigenValues, vectors = eigenVectors)) } # 動作確認 # 値そのものがあっているが, 符号が安定しない. # 符号が違うときがある理由がよくわからん....... A <- cov (iris[, -5]) print (myEigen.power.final (A)$values) $values [1] 4.22824171 0.24267075 0.07820950 0.02383509 $vectors [,1] [,2] [,3] [,4] [1,] 0.36138659 -0.65658877 -0.58202985 -0.3154872 [2,] -0.08452251 -0.73016143 0.59791083 0.3197231 [3,] 0.85667061 0.17337266 0.07623608 0.4798390 [4,] 0.35828920 0.07548102 0.54583143 -0.7536574 print (eigen (A)) $values [1] 4.22824171 0.24267075 0.07820950 0.02383509 $vectors [,1] [,2] [,3] [,4] [1,] 0.36138659 -0.65658877 -0.58202985 0.3154872 [2,] -0.08452251 -0.73016143 0.59791083 -0.3197231 [3,] 0.85667061 0.17337266 0.07623608 -0.4798390 [4,] 0.35828920 0.07548102 0.54583143 0.7536574
LU分解による逆反復法
# LU分解を利用した逆反復法 # つまり, 逆行列を使わないで, 最小固有値を求める. myEigen.power3 <- function (A, epsilon = 1e-10) { if (!exists ("myLUdecomposition")) { source ("myLUdecomposition.r") } if (!exists ("myBackSolve")) { source ("myBackSolve") } if (!exists ("myForwardSolve")) { source ("myForwardSolve") } # Au_k-1 = u_kでなくて, Au_k = u_k-1 でu_kを求める # AをLU分解しておくことで, 計算量を減らせる tmp <- myLUdecomposition (A, LDU = FALSE) U <- tmp$U L <- tmp$L #loop gain <- TRUE nr <- nrow (A) rtn <- rnorm (nr); rtn <- rtn / max (rtn) while (gain) { y <- myForwardSolve (cbind (L, rtn)) v <- myBackSolve (cbind (U, y)) u <- v / max (abs (v)) gain <- max (abs (u - rtn)) > epsilon rtn <- u } rtn <- rtn / sqrt (crossprod (rtn)) val <- crossprod (rtn, crossprod (A, rtn)) # v の最大値がAの最小固有値 # rtnがAの最小固有値に属する固有ベクトル return (list (val = val, vector = rtn)) } # 動作確認 A <- matrix (c(3, .01, .1, .01, 2, .01, .1, .01, 1), 3, 3, byrow = TRUE) print (myEigen.power3 (A)) print (eigen (A))
不完全コレスキー分解
# 不完全コレスキー分解 # http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%CF%A2%CE%A91%BC%A1%CA%FD%C4%F8%BC%B0%A1%A7%BB%B0%B3%D1%CA%AC%B2%F2 # http://emperor.yz.yamagata-u.ac.jp/pdf_dir/cholesky.pdf # A は実対称行列 # A が大規模疎行列のときに特に有効らしい. 例題がないので検証はどうする.... myICDecomposition <- function (A) { nr <- nrow (A) nc <- ncol (A) if (nr != nc) { warning ("またのお越しをお待ちしております. ") } # 一列目の処理 L <- matrix (0, nr, nc) L[ ,1] <- A[ , 1] d <- numeric (nr) d[1] <- 1 / L[1, 1] # 二列目以降 for (i in 2:(nr-1)) { for (j in i:nr) { if (abs (A[j, i]) > 1e-10) { tmp <- crossprod (L[i, 1:(i-1), drop = TRUE], d[1:(i-1)] * L[j, 1:(i-1), drop = TRUE]) lld <- A[j, i] - tmp L[j, i] <- lld if (i == j) { d[i] <- 1 / L[i, i] } } else { next } } } # 最後列 L[nr, nr] <- A[nr, nr] - crossprod (L[nr, 1:(nr-1)], d[1:(nr-1)] * L[nr, 1:(nr-1)]) d[nr] <- 1 / L[nr, nr] D <- diag (d) return (list (L = L, D = D)) } # 動作確認 # すうがくぶっくす 線形計算のやつとは同じものができてる A <- matrix (c (4, -1, -1, -0, -1, 4, 0, -1, -1, 0, 4, -1, 0, -1, -1, 4), 4, 4, byrow = TRUE) print (myICDecomposition (A))
