Rで線形代数:LU分解

# LU分解を実装する関数 (ピボット選択付き PA=LU)
fun_lu <- function(A) {
  n <- nrow(A)

  # L, U, P を初期化
  L <- diag(0, n) # まずはゼロ行列として開始
  U <- A # UはAのコピーとして開始
  P <- diag(1, n) # Pは単位行列として開始

  # ガウスの消去法を実行
  # kはピボットの列を指す
  for (k in 1:(n - 1)) {
    # 1. ピボット選択 (Partial Pivoting)
    # k列目のk行目以降で絶対値が最大の要素を探す
    pivot_val <- max(abs(U[k:n, k]))

    # 最大値がある行のインデックスを取得 (kからの相対位置なのでk-1を足す)
    max_row_index <- which.max(abs(U[k:n, k])) + k - 1

    # 2. 行の交換 (U, P, そしてLの計算済み部分)
    if (max_row_index != k) {
      # Uの行を交換
      temp_row_U <- U[k, ]
      U[k, ] <- U[max_row_index, ]
      U[max_row_index, ] <- temp_row_U

      # Pの行を交換 (行の入れ替え操作を記録)
      temp_row_P <- P[k, ]
      P[k, ] <- P[max_row_index, ]
      P[max_row_index, ] <- temp_row_P

      # Lの交換は少し注意が必要。これまでに計算済みの部分(k列より左)を交換する
      if (k > 1) {
        temp_row_L <- L[k, 1:(k - 1)]
        L[k, 1:(k - 1)] <- L[max_row_index, 1:(k - 1)]
        L[max_row_index, 1:(k - 1)] <- temp_row_L
      }
    }

    # 3. 前進消去
    # ピボット行より下の行の、k列目の要素を0にしていく
    for (i in (k + 1):n) {
      # 消去係数を計算し、Lに格納する
      L[i, k] <- U[i, k] / U[k, k]

      # Uの行から、(ピボット行 * 係数) を引く
      # これにより U[i, k] がゼロになる
      U[i, k:n] <- U[i, k:n] - L[i, k] * U[k, k:n]
    }
  }

  # Lの対角成分を1にする (ドゥーリトル法)
  diag(L) <- 1

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

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

# テスト用の行列を作成
A <- matrix(c(
  1, 1, 1,
  2, 2, 5,
  4, 6, 8
), nrow = 3, byrow = TRUE)

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

# 実装でLU分解を実行
lu_result <- fun_lu(A)

cat("\n--- 実装による結果 ---\n")
cat("P (置換行列):\n")
print(lu_result$P)
cat("\nL (下三角行列):\n")
print(lu_result$L)
cat("\nU (上三角行列):\n")
print(lu_result$U)

# 検算: PA = LU が成り立つか確認
PA <- lu_result$P %*% A
LU <- lu_result$L %*% lu_result$U

cat("\n--- 検算 ---\n")
cat("P %*% A の結果:\n")
print(PA)
cat("\nL %*% U の結果:\n")
print(LU)

# all.equalで浮動小数点数の誤差を考慮して比較
cat("\nPA と LU は等しいか:", all.equal(PA, LU), "\n")

lu_builtin_result <- Matrix::lu(A)
expanded_lu <- Matrix::expand(lu_builtin_result)
expanded_lu

# as.matrix()で標準の行列に変換

P_builtin <- as.matrix(expanded_lu$P)
L_builtin <- as.matrix(expanded_lu$L)
U_builtin <- as.matrix(expanded_lu$U)

cat("P (組み込み関数):\n")
print(P_builtin)
cat("\nL (組み込み関数):\n")
print(L_builtin)
cat("\nU (組み込み関数):\n")
print(U_builtin)