Rで線形代数:行列のランク

# ガウスの消去法により行列のランクを計算する関数
# A: 入力行列
# tol: 浮動小数点数の計算誤差を考慮するための許容誤差
fun_rank <- function(A, tol = 1e-10) {
  # 行列をコピーして、元の行列を変更しないようにする
  M <- A

  # 行数と列数を取得
  m <- nrow(M)
  n <- ncol(M)

  # ランク（非ゼロ行の数）を初期化
  rank <- 0
  # 現在注目しているピボット行
  pivot_row <- 1

  # 各列について処理を進める
  for (j in 1:n) {
    # 全ての行がピボットとして確定したら終了
    if (pivot_row > m) {
      break
    }

    # --- 部分ピボット選択 (Partial Pivoting) ---
    # 数値的安定性のために、pivot_row以降で絶対値が最大の要素をピボットとして探す
    i <- pivot_row
    max_row <- i
    max_val <- abs(M[i, j])

    # pivot_rowより下の行を探索
    if (i < m) {
      for (k in (i + 1):m) {
        if (abs(M[k, j]) > max_val) {
          max_val <- abs(M[k, j])
          max_row <- k
        }
      }
    }

    # ピボット（列の最大値）が許容誤差より小さい場合、
    # この列は線形従属とみなし、次の列に移る
    if (max_val < tol) {
      next
    }

    # --- 行の交換 ---
    # 見つけた最大値を持つ行を現在のpivot_rowと交換する
    if (max_row != pivot_row) {
      temp_row <- M[pivot_row, ]
      M[pivot_row, ] <- M[max_row, ]
      M[max_row, ] <- temp_row
    }

    # --- 前進消去 (Forward Elimination) ---
    # pivot_rowより下のすべての行について、j列目の要素を0にする
    # （ただし、pivot_rowが最終行でない場合）
    if (pivot_row < m) {
      for (i in (pivot_row + 1):m) {
        # 消去するための係数を計算
        factor <- M[i, j] / M[pivot_row, j]
        # 行の基本変形（下の行から、ピボット行の倍数を引く）
        M[i, ] <- M[i, ] - factor * M[pivot_row, ]
      }
    }

    # ピボットが見つかったので、ランクを1増やし、次のピボット行に進む
    rank <- rank + 1
    pivot_row <- pivot_row + 1
  }

  return(rank)
}