야코비법에 의한 행렬 대각화의 시각화 의 경우와 마찬가지로, 하우스 홀더법에 의한 삼중 대각화의 모습을 가시화한다.

처리로서는, 야코비법과 유사하고, 행렬 A에 하우스 홀더 행렬 P를 반복하여 양측으로부터 걸린다.

P를 A의 양측으로부터 곱하는 것으로, 원행렬의 특정의 행·열 벡터에 거울상 변환을 가한 것이 된다.

행렬 삼중 대각화 외에도 QR 분해를위한 방법으로도 사용됩니다.
cf. https://ko.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