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

加藤あいの結婚に衝撃を受けたのですが, みんな関心なくてひとり浮いてました。

あっ, どーも僕です。

会社員と会社役員の差って思ったよりもデカそうです。

追記

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

Rで普通の数値計算その４

やっぱりeigenはすごいなって話

どうやら, 数値計算初心者が固有値＆固有ベクトルと計算するときは次の３つが基本らしいです。

• 二分法　＋　シフト逆反復法 (myEigen.power.final)
• QR分解　＋　シフト逆反復法 (myEigen.HouseQR)
• ヤコビ法 (myEigen.Jacobi)

一応, 手書きで３つとも実装した（Rで普通の数値計算その１～その４のどっかにスクリプトあります。）ので簡単にeigenと比較してみました.

そもそも, みっつともちゃんと実装できているのか？

適当な10次実対称行列を突っ込んでみた結果がこれです.
みっつともeigenと同じ結果になりちゃんと計算できていることがわかります.

N <- 10
X <- matrix (rnorm (N * N), N, N)
print (myEigen.power.final (A = tcrossprod (X), upper = 1000, lower = 0, epsilon = 1e-10)$values)
 [1] 34.0735802 20.2278087 19.1208443 15.2971445
 [5]  9.1134511  6.4262798  1.8800622  0.9762552
 [9]  0.8072248  0.0110320
print (myEigen.HouseQR(tcrossprod (X), epsilon = 1e-10)) # 固有ベクトルは計算していない.
[1] 34.0735802 20.2278087 19.1208443 15.2971445
 [5]  9.1134511  6.4262798  1.8800622  0.9762552
 [9]  0.8072248  0.0110320
print (myEigen.Jacobi (tcrossprod (X), epsilon = 1e-10)$values)
 [1] 34.0735802 20.2278087 19.1208443 15.2971445
 [5]  9.1134511  6.4262798  1.8800622  0.9762552
 [9]  0.8072248  0.0110320
print (eigen(tcrossprod(X))$values)
 [1] 34.0735802 20.2278087 19.1208443 15.2971445
 [5]  9.1134511  6.4262798  1.8800622  0.9762552
 [9]  0.8072248  0.0110320

n = 3, 4, ..., 10のとき

n次の実対称行列を突っ込んだときの計算時間を測ってみました.
計算回数は一回です.

結果のグラフがこちら.
QR分解以外はほとんど時間がかかっていないことがわかります.
QR分解はヤコビ法よりも速いらしいので*1, どうやらわたしの実装があまそうです。

n = 10, 20, ..., 100のとき

前の計算でQR分解の実装があまいとことがわかったので, ここではQR分解は抜きでやってみます。前と同様に計算回数は1回です。

結果のグラフがこちら.(色が前と変わってますが....)
ヤコビ法 ＞ 二分法という計算時間になっています。そして, eigenは0秒のままです. すごいです。この結果だけみると手書きなら二分法が一番ということになりそうです.

小さい固有値があるとき

上の結果では, 二分法が最速となりました. なので, 二分法だけ実装すればいいのかというとどうやらそうでもなさそうです.

特異値分解を実装していたときに, svdのhelp では次のような行列にたいしてexamleが行われていました.

hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
X <- hilbert(9)[, 1:6]

で, これの外積の固有値を大きい方から6個ほしいのですが, 5・6個目がものすごく小さくなり, 二分法ではヤコビ法よりもeigenとの差があります. 大きな差には見えないかもしれませんが, 特異値分解経由で一般逆行列とかを計算してみるとMASSのginvと結果が変わってきてしまうことがわかります.

 print (myEigen.power.final (A = tcrossprod (X), upper = 1000, lower = 0, epsilon = 1e-12)$values)
[1] 2.783669e+00 7.693560e-02 4.944940e-04 1.176559e-06 1.051925e-09
[6] 5.136355e-13 2.568177e-13 1.284089e-13 6.420444e-14

print (myEigen.Jacobi (tcrossprod (X), epsilon = 1e-12)$values)
[1] 2.783669e+00 7.693560e-02 4.944940e-04 1.176559e-06 1.052216e-09
[6] 2.156099e-13 5.626614e-14 2.035486e-15 1.287170e-16

print (eigen(tcrossprod(X))$values)
[1]  2.783669e+00  7.693560e-02  4.944940e-04  1.176559e-06  1.052216e-09
[6]  2.740107e-13  1.791695e-17  1.189029e-17 -2.706279e-17

まとめ

eigenスゲーーーーー

次の修行本

ほんとは, 森先生とかの数値計算本に手を出すべきなのですが, fortranのコード読むのが大変なので, 線形計算はこの本くらいまでにしときます。
どちらも版落ちなのでアマゾンで合わせて3000円で買えました！！

行列の固有値―最新の解法と応用

Matrix Computations (Johns Hopkins Studies in the Mathematical Sciences) (Johns Hopkins Series in the Mathematical Sciences)

スクリプト

なにを参考にしたのか忘れましたが新たに３つのスクリプト

• QR分解経由で固有値
• 特異値分解
• 特異値分解経由で一般逆行列

QR分解経由で固有値計算

# QR分解経由で固有値を計算する
# 対称じゃなくても行けるはず
myEigen.HouseQR <- function (A, epsilon = 1e-12) {
    if (!exists("myQRdecompositionHouseholder")) {
        source ("myQRdecompotion_householder.r")
    }
    if (!exists("myEigen.Householder")) {
        source ("myEigen_Householder.r")
    }
    
    # ハウスホルダー変換による三重対角化
    # Aが非対称のときにはヘッセンベルグ行列になっている
    A.prime <- myEigen.Householder (A)
    
    # loops
    gain   <- TRUE
    values <- diag (A.prime) 
    while (gain) {
        # もっと高速化したいときは, QR分解でQだけ計算して, new = t(Q) old Qで計算
        tmp     <- myQRdecompositionHouseholder (A.prime)
        A.prime <- crossprod (t (tmp$R), tmp$Q)
        d       <- diag (A.prime)
        gain    <- max (abs (d - values)) > epsilon
        values  <- d
    }
    
    return (list (values = d))
}

# 動作確認
A <- cov (iris[, -5])
print (myEigen.HouseQR (A))
print (eigen (A)$values)

特異値分解

# http://www.kde.cs.tut.ac.jp/~atsushi/?p=468
# 特異値分解
# Aは実長方行列　
# num.sigvalは特異値をいくつ計算するか
mySVDdecomposition <- function (A, 
                                using.eigen = FALSE,
                                num.sigval = NULL,
                                epsilon = 1e-15) {
    nr <- nrow (A)
    nc <- ncol (A)
    if (is.null (num.sigval)) {
        num.sigval <- nc
    }
    if (!exists ("myEigen.Jacobi")) {
        source ("myEigen_Jacobi.r")
    }
    
    # 長い方で正方行列を作成して, それを固有値分解する
    # 固有値が0に近いと, ちゃんと計算できない(ハウスホルダー＋二分法での固有値計算が)....
    # ランク落ちの基準ってなに
    # 小さい固有値まで必要なときは, eigenを使う
    B <- tcrossprod (A)
    if (using.eigen) {
        B.eigen <- eigen (B)
    } else {
        # 自作だとヤコビが一番精度よく(小さい固有値が)計算できるのでヤコビで固有値を計算
        B.eigen <- myEigen.Jacobi (B, epsilon = epsilon)
    }
    d  <- sqrt (B.eigen$values[1:num.sigval])
    U  <- B.eigen$vectors[, 1:num.sigval]
    Vt <- sweep (crossprod (U, A), 1, d, "/")
    return (list (d = d, U = U, V = t (Vt)))
}

# 動作確認
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
X <- hilbert(9)[, 1:6]
print(s <- svd(X))
print (mySVDdecomposition (X, using.eigen = TRUE))

特異値分解経由で一般逆行列

# 特異値分解を利用した一般逆行列
# ペンローズの収束計算よりも精度はいいもよう..
myInverseGMatrix <- function (A, 
                              epsilon = 1e-12) {
    if (!exists ("mySVDdecomposition")) {
        source ("mySVDdecomposition.r")
    }
    tmp <- mySVDdecomposition (A, epsilon)
    U   <- tmp$U
    V   <- tmp$V
    W.t <- diag ((1 / tmp$d))
    return (crossprod (t (V), tcrossprod (W.t, U)))
    
}
# 動作確認
#N <- 5
#B <- matrix (rnorm (N * N), N, N)
# 有効数字は3桁くらいかな?
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
X <- hilbert(9)[, 1:6]
print (myInverseGMatrix (X))
print (ginv (X))

*1:山本哲郎：数値解析[増訂版], サイエンス社, p117