Rで線形代数:QR分解

Rで 線形代数:QR分解 を試みます。

QR分解の説明

QR分解とは、ある行列 A を、正規直交行列 Q と 上三角行列 R の積に分解することです。数式で表すと以下のようになります。

A = QR

ここで、それぞれの行列には次のような性質があります。

行列 `A`

分解したい元の行列です。ここでは、m行 n列の行列 (m >= n) で、列ベクトルが線形独立であるとします。

行列 `Q` (正規直交行列)

m行 n列の行列で、その列ベクトルが互いに正規直交であるという性質を持ちます。
正規直交とは、次の2つの条件を満たすことです。
1. 直交 (Orthogonal): Q のどの2つの異なる列ベクトルを選んでも、その内積（ドット積）は0になります。幾何学的には、ベクトルが互いに直角に交わっていることを意味します。
2. 正規 (Normal): Q のすべての列ベクトルの長さ（ノルム）が1です。
この性質から、Qの転置行列 t(Q) と Q の積は、単位行列 I になります (t(Q) %*% Q = I)。

行列 `R` (上三角行列)

n行 n列の正方行列で、上三角行列です。
上三角行列とは、行列の対角成分より左下の要素がすべて0である行列のことです。

[r_11, r_12, r_13]
[   0, r_22, r_23]
[   0,    0, r_33]

QR分解の考え方：グラム・シュミットの正規直交化法

ここでは、QR分解を求めるためのアルゴリズムに「グラム・シュミットの正規直交化法」を利用します。これは、元の行列 A の列ベクトル（線形独立なベクトルの集まり）を、互いに直交するベクトルの集まりに変換し、最後に各ベクトルの長さを1に正規化する、というプロセスです。

A の1番目の列ベクトル a_1 を取り出し、長さを1にして q_1 を作ります。
A の2番目の列ベクトル a_2 を取り出します。a_2 から、q_1 と同じ方向の成分（a_2 の q_1 への射影）を引き算します。これにより q_1 と直交するベクトルができます。これを正規化して q_2 を作ります。
A の3番目の列ベクトル a_3 から、q_1 方向の成分と q_2 方向の成分を引き算し、残ったベクトルを正規化して q_3 を作ります。
この操作を A のすべての列ベクトルに対して繰り返します。

こうして得られた正規直交ベクトル q_1, q_2, ..., q_n を並べたものが Q 行列です。 R 行列の各要素 r_ij は、この変換プロセスの「係数」を記録したものです。具体的には、R = t(Q) %*% A という関係式から求めることができ、この計算方法により R は上三角行列になります。

R言語による実装

以下に、グラム・シュミットの正規直交化法を用いてQR分解を実装したRコードを示します。

# QR分解を実装する関数
fun_qr <- function(A) {
  # 入力行列Aの次元を取得
  m <- nrow(A)
  n <- ncol(A)

  # QとRを初期化（ゼロ行列で）
  Q <- matrix(0, m, n)
  R <- matrix(0, n, n)

  # グラム・シュミットの正規直交化法
  for (j in 1:n) {
    # 1. 元のベクトルvを取得
    v <- A[, j]

    # 2. 既に計算済みのQの列ベクトルとの直交化
    #    j=1のときはこのループは実行されない
    if (j > 1) {
      for (i in 1:(j - 1)) {
        # R[i, j]は、A[,j]のQ[,i]への射影の係数
        R[i, j] <- t(Q[, i]) %*% A[, j] # 内積を計算

        # A[,j]からQ[,i]方向の成分を引く
        v <- v - R[i, j] * Q[, i]
      }
    }

    # 3. 正規化
    #    vのノルム（長さ）を計算し、Rの対角成分に格納
    R[j, j] <- sqrt(sum(v^2))

    #    vをノルムで割り、Qのj列目のベクトル（正規直交基底）とする
    Q[, j] <- v / R[j, j]
  }

  # 結果をリストで返す
  return(list(Q = Q, R = R))
}

コードのテスト

なお、テストする(QR分解する)行列は列ベクトルが線形独立の行列とします。

3×3の行列

A <- matrix(c(
  1, -1, 4,
  1, 4, -2,
  1, 4, 2,
  1, -1, 0
), nrow = 4, byrow = TRUE)

# QR分解を実行
qr_result <- fun_qr(A)

cat("--- 結果 ---\n")
cat("Q行列:\n")
print(qr_result$Q)
cat("\nR行列:\n")
print(qr_result$R)

--- 結果 ---
Q行列:
     [,1] [,2] [,3]
[1,]  0.5 -0.5  0.5
[2,]  0.5  0.5 -0.5
[3,]  0.5  0.5  0.5
[4,]  0.5 -0.5 -0.5

R行列:
     [,1] [,2] [,3]
[1,]    2    3    2
[2,]    0    5   -2
[3,]    0    0    4

検算1: Q %*% R が A に戻るか確認

A_reconstructed <- qr_result$Q %*% qr_result$R
cat("\n検算1: Q %*% R の結果（Aに復元されるか）:\n")
print(A_reconstructed)
# all.equalで浮動小数点数の誤差を考慮して比較
cat("\nAとQ%*%Rは等しいか:", all.equal(A, A_reconstructed), "\n")


検算1: Q %*% R の結果（Aに復元されるか）:
     [,1] [,2] [,3]
[1,]    1   -1    4
[2,]    1    4   -2
[3,]    1    4    2
[4,]    1   -1    0

AとQ%*%Rは等しいか: TRUE

検算2: Qが正規直交行列か確認 (t(Q) %*% Q が単位行列になるか)

I <- t(qr_result$Q) %*% qr_result$Q
cat("\n検算2: t(Q) %*% Q の結果（単位行列になるか）:\n")
print(I)
cat("\nt(Q)%*%Qは単位行列と等しいか:", all.equal(I, diag(ncol(A))), "\n")


検算2: t(Q) %*% Q の結果（単位行列になるか）:
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

t(Q)%*%Qは単位行列と等しいか: TRUE

参考: Rの組み込み関数によるQR分解

注意: QR分解の結果は一意ではありません(Rの対角成分をすべて正である(r_ii > 0)と制約しない限り)。Qの列ベクトルの符号を反転させ、それに応じてRの行の符号も反転させたものも、有効なQR分解です。そのため、実装と組み込み関数の結果で符号が異なる場合があります。

cat("--- Rの組み込み関数による結果 ---\n")
qr_builtin <- qr(A)
Q_builtin <- qr.Q(qr_builtin)
R_builtin <- qr.R(qr_builtin)

cat("Q行列 (組み込み関数):\n")
print(Q_builtin)
cat("\nR行列 (組み込み関数):\n")
print(R_builtin)

--- Rの組み込み関数による結果 ---
Q行列 (組み込み関数):
     [,1] [,2] [,3]
[1,] -0.5  0.5 -0.5
[2,] -0.5 -0.5  0.5
[3,] -0.5 -0.5 -0.5
[4,] -0.5  0.5  0.5

R行列 (組み込み関数):
     [,1] [,2] [,3]
[1,]   -2   -3   -2
[2,]    0   -5    2
[3,]    0    0   -4

以上です。