Rで線形代数:固有値分解

# -------------------------------------------------------------------
# 1. 固有値・固有ベクトルを求めるQR法
# -------------------------------------------------------------------
eigen_qr_method <- function(A, num_simulations = 100) {
  if (nrow(A) != ncol(A)) stop("入力は正方行列である必要があります。")
  Ak <- A
  n <- ncol(A)
  P <- diag(1, n)
  for (i in 1:num_simulations) {
    qr_res <- qr(Ak)
    Q <- qr.Q(qr_res)
    R <- qr.R(qr_res)
    Ak <- R %*% Q
    P <- P %*% Q
  }
  return(list(values = diag(Ak), vectors = P))
}

# -------------------------------------------------------------------
# 2. 固有値分解を組み立てる関数
# -------------------------------------------------------------------
eigen_decomposition <- function(A) {
  # まず、固有値と固有ベクトルを求める
  eigen_res <- eigen_qr_method(A)

  # 固有値から対角行列Dを作成
  D <- diag(eigen_res$values)

  # 固有ベクトル行列P
  P <- eigen_res$vectors

  # Pの逆行列を計算 (組み込み関数solveを使用)
  # solve()が失敗する場合、行列は正則でない（対角化不可能）
  P_inv <- tryCatch(
    {
      solve(P)
    },
    error = function(e) {
      stop("固有ベクトル行列Pが正則でなく、逆行列を計算できませんでした。行列は対角化不可能です。")
    }
  )

  # 結果をリストで返す
  return(list(P = P, D = D, P_inv = P_inv))
}


# ===================================================================
# テストと比較
# ===================================================================
# 必ず対角化可能な対称行列をサンプルとします
A <- matrix(c(
  6, 2, -1,
  2, 5, 1,
  -1, 1, 4
), nrow = 3, byrow = TRUE)

cat("--- 固有値分解のテスト ---\n")
cat("元の行列 A:\n")
print(A)

# --- 実装による分解 ---
cat("\n\n--- 実装による結果 ---\n")
decomp_s <- eigen_decomposition(A)
cat("固有ベクトル行列 P:\n")
print(decomp_s$P)
cat("\n固有値行列 D:\n")
print(decomp_s$D)
cat("\nPの逆行列 P⁻¹:\n")
print(decomp_s$P_inv)

# 検算
A_reconstructed <- decomp_s$P %*% decomp_s$D %*% decomp_s$P_inv
cat("\n検算: P %*% D %*% P⁻¹ は A と等しいか:", all.equal(A, A_reconstructed), "\n")


# --- Rの組み込み関数による結果との比較 ---
cat("\n\n--- Rの組み込み関数 eigen() を使った分解 ---\n")
eigen_b <- eigen(A)
P_b <- eigen_b$vectors
D_b <- diag(eigen_b$values)
P_inv_b <- solve(P_b)
A_reconstructed_b <- P_b %*% D_b %*% P_inv_b
cat("検算: P %*% D %*% P⁻¹ は A と等しいか:", all.equal(A, A_reconstructed_b), "\n")

# 注意:
# 固有値・固有ベクトルの順序や、固有ベクトルの符号はアルゴリズムによって異なる場合があります。
# 対称行列ですので P⁻¹ は t(P) と等しくなることも確認できます。
cat("\n対称行列の特性: P⁻¹ は t(P) と等しいか:", all.equal(decomp_s$P_inv, t(decomp_s$P)), "\n")