Rで線形代数:直交行列

# --------------------------------------------------
# 直交行列のシミュレーション
# --------------------------------------------------

# --- シナリオ1: 直交行列の作成と性質の確認 ---
cat("--- シナリオ1: 直交行列の作成と性質の確認 ---\n")

# 1. 元となる適当な行列を作成
seed <- 20250624
set.seed(seed)
A <- matrix(rnorm(9), nrow = 3)

# 2. QR分解を用いて直交行列Qを抽出
#    qr()関数でQR分解を行い、qr.Q()でQを取り出す
Q <- qr.Q(qr(A))

cat("生成された直交行列 Q:\n")
print(Q)

# 3. 性質の確認
# (a) Q^T * Q = I となるか？
cat("\n(a) Q^T * Q の計算結果:\n")
QtQ <- t(Q) %*% Q
print(round(QtQ, 10)) # 計算誤差を考慮して丸める
cat("→ 単位行列になりました。\n")

# (b) Qの転置行列が、Qの逆行列と一致するか？
cat("\n(b) 転置行列と逆行列の比較:\n")
Q_transpose <- t(Q)
Q_inverse <- solve(Q)

cat("Qの転置行列:\n")
print(round(Q_transpose, 5))
cat("Qの逆行列:\n")
print(round(Q_inverse, 5))
cat("→ 両者が一致していることがわかります。\n")

# (c) 行列式が 1 または -1 になるか？
cat("\n(c) 行列式の計算結果:\n")
det_Q <- det(Q)
cat("det(Q) =", det_Q, "\n")
cat("→ 行列式が 1 になりました (これは鏡映を含む変換を意味します)。\n\n")


# --- シナリオ2: 幾何学的な意味の可視化 ---
cat("--- シナリオ2: 幾何学的な意味の可視化 ---\n")

# 2次元の例で可視化する
# (a) 回転を表す直交行列 (det=1)
theta <- pi / 6 # 30度
R_matrix <- matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), nrow = 2)
cat("回転行列 R (det=", round(det(R_matrix), 2), ") Figure 1:\n")
print(R_matrix)

# (b) 鏡映を表す直交行列 (det=-1)
#    y軸に対する鏡映
F_matrix <- matrix(c(-1, 0, 0, 1), nrow = 2)
cat("\n鏡映行列 F (det=", round(det(F_matrix), 2), ") Figure 2:\n")
print(F_matrix)

# 描画用の関数
plot_rigid_transform <- function(M, shape) {
  transformed_shape <- t(M %*% t(shape))

  all_points <- rbind(shape, transformed_shape)
  plot_range <- range(all_points) * 1.2
  plot(NULL,
    xlim = plot_range, ylim = plot_range, asp = 1,
    main = paste("Transformation by Orthogonal Matrix (det=", round(det(M), 2), ")"), xlab = "", ylab = ""
  )

  polygon(shape, col = "#AAAAFF80", border = "blue")
  polygon(transformed_shape, col = "#FFAAAA80", border = "red")

  abline(h = 0, v = 0, lty = 2, col = "grey")
  legend("topleft",
    legend = c("Original", "Transformed"),
    fill = c("#AAAAFF80", "#FFAAAA80"), bty = "n"
  )
}

# 変換対象の図形 (非対称なL字型)
l_shape <- rbind(c(0, 0), c(1, 0), c(1, 0.5), c(0.5, 0.5), c(0.5, 2), c(0, 2), c(0, 0))

# 回転のプロット
plot_rigid_transform(R_matrix, l_shape)

# 鏡映のプロット
plot_rigid_transform(F_matrix, l_shape)