Rで線形代数:固有値と固有ベクトル

# べき乗法により、絶対値最大の固有値とその固有ベクトルを求める関数
power_method <- function(A, num_simulations = 100, tol = 1e-10) {
  # ランダムな初期ベクトルを生成
  b_k <- rnorm(ncol(A))

  for (i in 1:num_simulations) {
    # 1. Ax を計算
    b_k1 <- A %*% b_k

    # 2. 正規化（長さを1にする）
    b_k1_norm <- sqrt(sum(b_k1^2))
    b_k <- b_k1 / b_k1_norm
  }

  # 収束した固有ベクトル b_k を使ってレイリー商から固有値を計算
  # λ = (t(v)Av) / (t(v)v)  (vが単位ベクトルなら t(v)Av)
  lambda <- t(b_k) %*% A %*% b_k

  return(list(eigenvalue = as.numeric(lambda), eigenvector = b_k))
}

# ---- べき乗法のテストと比較 ----
A_test <- matrix(c(2, 1, 1, 2), nrow = 2) # 固有値は3と1
cat("--- べき乗法のテスト (A_test) ---\n")
cat("元の行列 A_test:\n")
print(A_test)

# 実装
result_pm <- power_method(A_test)
cat("\n[実装 (べき乗法)]\n")
cat("絶対値最大の固有値:", result_pm$eigenvalue, "\n")
cat("対応する固有ベクトル:\n")
print(result_pm$eigenvector)

# Rの組み込み関数
eigen_builtin_test <- eigen(A_test)
cat("\n[Rの組み込み関数 eigen()]\n")
cat("全ての固有値:", eigen_builtin_test$values, "\n")
# 最大の固有値と対応する固有ベクトルを抽出
max_eigenvalue_builtin <- eigen_builtin_test$values[1] # eigen()は降順で返す
max_eigenvector_builtin <- eigen_builtin_test$vectors[, 1]
cat("絶対値最大の固有値:", max_eigenvalue_builtin, "\n")
cat("対応する固有ベクトル:\n")
print(max_eigenvector_builtin)

# 比較
cat("\n[比較]\n")
cat("固有値は等しいか:", all.equal(result_pm$eigenvalue, max_eigenvalue_builtin), "\n")
cat("固有ベクトルは（符号の違いは別として）等しいか:", all.equal(abs(result_pm$eigenvector), abs(matrix(max_eigenvector_builtin))), "\n")

# QR法により、すべての固有値と固有ベクトルを求める関数
# 固有値がすべて実数である場合を想定しています。
eigen_qr_method <- function(A, num_simulations = 100) {
  # 入力が正方行列である必要があります。
  if (nrow(A) != ncol(A)) stop("入力は正方行列である必要があります。")

  Ak <- A
  n <- ncol(A)

  # 全てのQ行列の積を保存するための行列Pを初期化
  # Pは最終的に固有ベクトルを列として持つ行列になる
  P <- diag(1, n)

  for (i in 1:num_simulations) {
    # 1. AkのQR分解
    qr_res <- qr(Ak)
    Q <- qr.Q(qr_res)
    R <- qr.R(qr_res)

    # 2. Ak+1 = R %*% Q で更新
    Ak <- R %*% Q

    # 3. Pに行列Qを右から掛けて更新
    P <- P %*% Q
  }

  # 収束した行列の対角成分が固有値
  eigenvalues <- diag(Ak)

  # Pが固有ベクトル行列
  eigenvectors <- P

  return(list(values = eigenvalues, vectors = eigenvectors))
}

# ---- QR法のテストと比較 ----
# 対称行列を作成
A_sym <- matrix(c(
  6, 2, -1,
  2, 5, 1,
  -1, 1, 4
), nrow = 3, byrow = TRUE)

cat("--- QR法による全固有値・固有ベクトルの計算 ---\n")
cat("元の行列 A:\n")
print(A_sym)


# --- 実装による計算 ---
cat("\n\n--- 実装 (QR法) による結果 ---\n")
result_qr <- eigen_qr_method(A_sym)

# QR法の結果はソートされていないので、比較のために降順にソートする
# まずソート順序を取得
order_qr <- order(result_qr$values, decreasing = TRUE)
# その順序に従って固有値と固有ベクトルを並べ替える
sorted_values_qr <- result_qr$values[order_qr]
sorted_vectors_qr <- result_qr$vectors[, order_qr]

cat("固有値 (降順ソート):\n")
print(sorted_values_qr)
cat("\n固有ベクトル (固有値の大きい順に並べ替え):\n")
print(sorted_vectors_qr)


# --- Rの組み込み関数による計算 ---
cat("\n\n--- Rの組み込み関数 eigen() による結果 ---\n")
eigen_builtin <- eigen(A_sym)
cat("固有値:\n")
print(eigen_builtin$values)
cat("\n固有ベクトル:\n")
print(eigen_builtin$vectors)


# --- 最終的な比較 ---
cat("\n\n--- 最終比較 ---\n")
# 固有値の比較
cat(
  "固有値の集合は等しいか:",
  all.equal(sorted_values_qr, eigen_builtin$values), "\n"
)

# 固有ベクトルの比較
# 符号の違いを吸収するために絶対値で比較する
cat(
  "固有ベクトルは（符号の違いは別として）等しいか:",
  all.equal(abs(sorted_vectors_qr), abs(eigen_builtin$vectors)), "\n"
)


# --- 検算: Av = λv が成り立つか ---
cat("\n\n--- 検算 (実装の結果) ---\n")
# 1番目の固有値・固有ベクトルで確認
lambda1 <- sorted_values_qr[1]
v1 <- sorted_vectors_qr[, 1]
Av1 <- A_sym %*% v1
lambda1_v1 <- lambda1 * v1

cat("A %*% v1:\n")
print(Av1)
cat("λ1 * v1:\n")
print(lambda1_v1)
cat("Av = λv は成り立つか:", all.equal(Av1, matrix(lambda1_v1)), "\n")