Rで線形代数:シュミットの直交化法

# --------------------------------------------------
# シュミットの直交化法のシミュレーション
# --------------------------------------------------

# --- 準備 ---

# 1. 元となる線形独立なベクトルの組を3つ定義
v1 <- c(1, 1, 0)
v2 <- c(1, 0, 1)
v3 <- c(0, 1, 1)

# 元のベクトルをまとめた行列 V
V <- cbind(v1, v2, v3)
cat("元のベクトル (行列Vの列):\n")
print(V)

# ベクトルの内積を計算するヘルパー関数
dot_product <- function(a, b) {
  sum(a * b)
}

# 射影ベクトルを計算するヘルパー関数
# v の u 上への射影を計算
projection <- function(v, u) {
  scalar_projection <- dot_product(v, u) / dot_product(u, u)
  return(scalar_projection * u)
}


# --- Step by Step で直交化を実行 ---

cat("\n--- 直交化のプロセス ---\n")

# u1 の計算
u1 <- v1
cat("u1 (v1をそのまま使用):\n")
print(u1)

# u2 の計算
proj_u1_v2 <- projection(v2, u1)
u2 <- v2 - proj_u1_v2
cat("\nu2 (v2 - proj_u1(v2)):\n")
print(u2)

# u3 の計算
proj_u1_v3 <- projection(v3, u1)
proj_u2_v3 <- projection(v3, u2)
u3 <- v3 - proj_u1_v3 - proj_u2_v3
cat("\nu3 (v3 - proj_u1(v3) - proj_u2(v3)):\n")
print(u3)

# --- 確認1: u1, u2, u3 が互いに直交しているか？ ---
cat("\n--- 確認1: u1, u2, u3の直交性の確認 (内積が0か？) ---\n")
cat("u1・u2 =", round(dot_product(u1, u2), 10), "\n")
cat("u1・u3 =", round(dot_product(u1, u3), 10), "\n")
cat("u2・u3 =", round(dot_product(u2, u3), 10), "\n")
cat("→ すべての内積が0になり、互いに直交していることが確認できました。\n")


# --- 最終ステップ: 正規化 (長さを1にする) ---
cat("\n--- 最終ステップ: 正規化 ---\n")

# ノルム（ベクトルの長さ）を計算する関数
norm_vec <- function(v) {
  sqrt(sum(v^2))
}

e1 <- u1 / norm_vec(u1)
e2 <- u2 / norm_vec(u2)
e3 <- u3 / norm_vec(u3)

cat("正規直交ベクトル e1:\n")
print(e1)
cat("正規直交ベクトル e2:\n")
print(e2)
cat("正規直交ベクトル e3:\n")
print(e3)

# これらの正規直交ベクトルを列とする行列Qを作成
Q <- cbind(e1, e2, e3)
cat("\n作成された行列 Q (列がe1, e2, e3):\n")
print(Q)

# --- 確認2: 行列Qが直交行列になっているか？ ---
cat("\n--- 確認2: Qが直交行列か？ (Q^T * Q = I) ---\n")
QtQ <- t(Q) %*% Q
print(round(QtQ, 10))
cat("→ 単位行列となり、Qが直交行列であることが確認できました。\n")

# --- QR分解との比較 ---
cat("\n--- QR分解の結果との比較 ---\n")
Q_qr <- qr.Q(qr(V))
cat("Rのqr()関数で得られた直交行列:\n")
print(Q_qr)
cat("→ 自作したQと（1列目の符号は異なっていますが）同じ基底を張っていることがわかります。\n")
# ※QR分解のアルゴリズムによっては、ベクトルの符号が逆になることがありますが、
#   張られる空間や直交性は同じです。