Rで線形代数:コレスキー分解

# コレスキー分解を実装する関数
fun_cholesky <- function(A) {
  n <- nrow(A)

  # --- 前提条件のチェック ---
  # 1. 正方行列か？
  if (ncol(A) != n) {
    stop("入力は正方行列である必要があります。")
  }
  # 2. 対称行列か？ (isSymmetricは誤差を考慮してくれる)
  if (!isSymmetric(A, tol = 1e-10)) {
    stop("入力された行列は対称行列ではありません。")
  }

  # Lを初期化 (ゼロ行列で)
  L <- matrix(0, n, n)

  # --- コレスキー分解の計算 ---
  # Lの要素を左上の要素から順に計算していく
  for (i in 1:n) {
    for (j in 1:i) {
      # Lの対角要素 (j == i) の場合
      if (j == i) {
        # A[i,i] から、それまでに計算したLのi行目の要素の二乗和を引く
        sum_k <- 0
        if (i > 1) {
          sum_k <- sum(L[i, 1:(i - 1)]^2)
        }

        # 平方根の中身を計算
        val_under_sqrt <- A[i, i] - sum_k

        # 3. 正定値性のチェック (平方根の中身が負にならないか)
        if (val_under_sqrt < 0) {
          stop("入力された行列は正定値ではありません。")
        }

        L[i, i] <- sqrt(val_under_sqrt)
      } else { # Lの非対角要素 (j < i) の場合

        # A[i,j] から、Lのi行目とj行目の内積を引く
        sum_k <- 0
        if (j > 1) {
          sum_k <- sum(L[i, 1:(j - 1)] * L[j, 1:(j - 1)])
        }

        L[i, j] <- (A[i, j] - sum_k) / L[j, j]
      }
    }
  }

  return(L)
}

# ---- テスト ----

# テスト用の対称正定値行列を作成
A <- matrix(c(
  4, 12, -16,
  12, 37, -43,
  -16, -43, 98
), nrow = 3, byrow = TRUE)

cat("--- コレスキー分解のテスト ---\n")
cat("元の行列 A:\n")
print(A)

# 実装でコレスキー分解を実行
cat("\n--- 実装による結果 ---\n")
L <- fun_cholesky(A)
cat("L (実装):\n")
print(L)


# 検算: L %*% t(L) が A に戻るか確認
A_reconstructed <- L %*% t(L)
cat("\n検算: L %*% t(L) の結果（Aに復元されるか）:\n")
print(A_reconstructed)
cat("AとL %*% t(L)は等しいか:", all.equal(A, A_reconstructed), "\n")


# 参考: Rの組み込み関数によるコレスキー分解
cat("\n\n--- Rの組み込み関数 chol() による結果 ---\n")
# chol() はデフォルトで「上三角行列」を返すので、転置してLにする
L_builtin <- t(chol(A))
cat("L (組み込み関数):\n")
print(L_builtin)


# --- 最終的な比較 ---
cat("\n\n--- 最終比較 ---\n")
cat("実装Lと組み込み関数のLは等しいか:", all.equal(L, L_builtin), "\n")