Rで普通の数値解析その3

あっ, どーも僕です.


そういえば, 2019年にハンドボール女子世界選手権が日本開催に決定したのでここで宣伝しときますね!!


日本代表女子の試合はyoutubeに動画が上がっています!!




ハンドボールをもっと楽しみたい人は本場のハンドボールを見てください!!
full match のところにあるパリサンジェルマンバルセロナの試合がオススメです!
EHFTV | LAOLA1.tv


Rで普通の数値解析その3

追記

スクリプトの最新版はこちらから。
Rで普通の数値解析





それはおいといて, 今回の修行本はこちら.

線形計算 (すうがくぶっくす)

線形計算 (すうがくぶっくす)

こちらから、次の3つを実装してみまいした.

これの中で, 特にハウスホルダー変換+二分法からのシフト逆反復による固有値固有ベクトルが, 数値解析の小技集的な存在なので, 線形計算をはじめるときはこれを目標にいろいろ実装してみると楽しいかと思います.


また, いままでの普通のスクリプトをまとめて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))