r - R コードで Householder QR 因数分解関数を作成する

Question

Rで行列のQR因数分解を見つけるためのコードに取り組んでいます.

X <- structure(c(0.8147, 0.9058, 0.127, 0.9134, 0.6324, 0.0975, 0.2785, 
0.5469, 0.9575, 0.9649, 0.1576, 0.9706, 0.9572, 0.4854, 0.8003
), .Dim = c(5L, 3L))


myqr <- function(A) {
  n <- nrow(A)
  p <- ncol(A)
  Q <- diag(n)
  Inp <- diag(nrow = n, ncol = p)

  for(k in c(1:ncol(A))) {
    # extract the kth column of the matrix
    col<-A[k:n,k]
    # calculation of the norm of the column in order to create the vector
    norm1<-sqrt(sum(col^2))
    # Define the sign positive if a1 > 0 (-) else a1 < 0(+)  
    sign <- ifelse(col[1] >= 0, -1, +1)  
    # Calculate of the vector a_r
    a_r <- col - sign * Inp[k:n,k] * norm1
    # beta = 2 / ||a-r||^2  
    beta <- 2 / sum(t(a_r) %*% a_r)
    # the next line of code calculates the matrix Q in every step
    Q <- Q - beta *Q %*% c(rep(0,k-1),a_r) %*% t(c(rep(0,k-1),a_r))    
    # calculates the matrix R in each step
    A[k:n,k:p] <- A[k:n,k:p] - beta * a_r %*% t(a_r) %*% A[k:n,k:p]
    }

  list(Q=Q,R=A)
  }

しかし、ここではH、世帯主の反射を表す行列をすべてのステップで計算していません。またA、すべてのステップで行列を計算していません。

としてH = I - 2 v v'、 を掛けるQと得られます

QH = Q - 2 (Qv) v'    // multiplication on the left
HQ = Q - 2 v (Q'v)'    // multiplication on the right

これで、この操作はすべてのステップで機能するはずです。ただし、最初のマトリックスHと2番目のマトリックスを考慮すると、H1これらのマトリックスは最初のマトリックスよりも小さくなります。それを避けるために、次のコード行を使用しました。

 Q <- Q - beta * Q %*% c(rep(0,k-1),a_r) %*% t(c(rep(0,k-1),a_r))

しかし、すべてのステップでゼロa_rの最初のエントリを持つ新しいベクトルを生成すると、コードがうまく機能する理由がわかりません。k