ハウスホルダー法による行列三重対角化の可視化

17843 ワード
R 線形代数 R テキストリンク

ヤコビ法による行列対角化の可視化の時と同様に、ハウスホルダー法による三重対角化の様子を可視化する。

処理としては、ヤコビ法と似ていて、行列Aにハウスホルダー行列Pを繰り返し両側からかける。

PをAの両側からかけることで、元行列の特定の行・列ベクトルに鏡像変換をかけたことになる。

行列三重対角化以外に、QR分解のための手法としても用いられる。
cf. https://ja.wikipedia.org/wiki/ハウスホルダー変換

library("fields")
library("stringr")

HouseholderMatrix <- function(A, step){
    x <- A[, step]
    y <- rep(0, length(x))
    y[1:step] <- A[1:step, step]
    s <- sqrt(sum(A[setdiff(1:nrow(A), 1:step), step]^2))
    y[step+1] <- - s
    u <- (x - y) / norm(as.matrix(x - y), "F")
    P <- diag(1, nrow(A))
    P <- P - 2 * u %*% t(u)
    P
}

plotA <- function(A, step){
  image.plot(log10(A[, ncol(A):1]+1),
    main=paste0("k = ", step),
    xaxt="n", yaxt = "n")
  if(step != nrow(A)-1){
      pos <- seq(0, 1, length = nrow(A))
      text(pos[step], 1-pos[(step+2):length(pos)], "×")
      text(pos[(step+2):length(pos)], 1-pos[step], "×")    
  }
}

Householder <- function(A, viz=FALSE, viz.dir=""){
  if(viz){
    if(viz.dir==""){
      plotA(A, 0)
      Sys.sleep(0.2)
    }else{
      png(file=paste0(viz.dir, "/Step0000.png"))
      plotA(A, 0)
      dev.off()
    }
  }
  U <- diag(1, nrow(A))
  step <- 1
  for(k in seq_len(nrow(A)-1)){
    P <- HouseholderMatrix(A, step)
    U <- U %*% P
    A <- P %*% A %*% P # Update
    if(viz){
      if(viz.dir==""){
        plotA(A, step)
        Sys.sleep(0.2)
      }else{
        png(file=paste0(viz.dir, "/Step",
          str_sub(paste0("000", step), start = -4),
          ".png"))
        plotA(A, step)
        dev.off()
      }
    }
    step <- step + 1
  }
  list(U=U, Sigma=A)
}

# 30×30行列で実験
dir.create("HouseHolder", showWarnings=FALSE)
A <- matrix(runif(30*30), nrow=30, ncol=30)
A <- t(A) %*% A # 対称行列にする
out <- Householder(A, viz=TRUE, viz.dir="./Householder")

xは次のステップで0になる予定箇所。
行列Aが徐々に三重対角要素を残して疎な行列になっていく様子がわかる。

ヤコビ法の対角化と比べると、三重対角化は(行列Aの次数)-1分だけ処理をすればいいので、計算がすぐ終わることがわかる。
この三重対角行列にさらに二分法など別の手法を適用することにより、元の行列Aの固有値・固有ベクトルがもとめる。
密行列と比べて要素数も少ないので計算がかなり楽になる。

参考

計算による線形代数
Matrix Computation