Rで線形代数:極分解

Rで 線形代数:極分解 を試みます。

1. 行列の極分解の説明

極分解とは、任意の正方行列 A を、正規直交行列 U と対称半正定値行列 P の積に分解することです。

A = UP

これは、複素数の極形式 z = re^(iθ) と類比できます。

  • r（大きさ、非負の実数）が、対称半正定値行列 P（伸縮）に対応します。
  • e^(iθ)（回転、大きさ1）が、正規直交行列 U（回転）に対応します。

行列の各要素

  • 行列 A: 分解したい n x n の正方行列。逆行列を持つ（正則である）必要はありません。
  • 行列 U (正規直交行列):

    • U は n x n の正規直交行列（直交行列とも呼ばれます）です。
    • t(U)U = U t(U) = I（単位行列）を満たします。
    • 幾何学的には、回転または鏡映 を表します。det(U)=+1 なら純粋な回転、det(U)=-1 なら回転と鏡映の組み合わせです。
  • 行列 P (対称半正定値行列):

    • P は n x n の対称 (P = t(P)) かつ半正定値な行列です。
    • 半正定値とは、すべての固有値が非負（0以上）であることを意味します。
    • 幾何学的には、空間の伸縮（ストレッチ） を表します。P は、どの方向にどれだけ伸び縮みするかを決定します。

極分解の考え方：SVDとの関係

極分解は、特異値分解（SVD） を使って導出できます。

A のSVDを A = U_svd Σ t(V_svd) とします。

この式を変形します。

A = U_svd Σ t(V_svd)

= U_svd Σ (V_svd t(V_svd)) t(V_svd) （V_svdは正規直交行列なので V_svd t(V_svd) = I を挿入）

= (U_svd t(V_svd)) (V_svd Σ t(V_svd))

ここで、2つの部分に分けて見てみます。

  1. U = U_svd t(V_svd):

    • U_svd と V_svd はどちらも正規直交行列です。正規直交行列の積（や転置）もまた正規直交行列になります。
    • したがって、この部分が極分解の U（回転行列）になります。
  2. P = V_svd Σ t(V_svd):

    • この P は対称行列です (t(P) = t(t(V_svd) Σ t(V_svd)) = V_svd t(Σ) t(t(V_svd)) = V_svd Σ t(V_svd) = P)。
    • また、P の固有値は Σ の対角成分（つまり特異値）となり、すべて非負です。したがって、P は対称半正定値行列です。
    • この部分が極分解の P（伸縮行列）になります。

このように、SVDを計算すれば、その構成要素を組み合わせることで極分解が求められます。

右極分解と左極分解

分解の順序によって2つの種類があります。

  • 右極分解: A = UP (先に伸縮、次に回転)
  • 左極分解: A = P'U (先に回転、次に伸縮)
  • P' は P とは異なる伸縮行列 (P' = U P t(U)) ですが、U は同じ回転行列です。

2. R言語による実装

# -------------------------------------------------------------------
# 行列の極分解を行う関数
# (SVDの計算にはRの組み込み関数 svd() を利用)
# -------------------------------------------------------------------
polar_decomposition <- function(A) {
  # 入力は正方行列である必要があります
  if (nrow(A) != ncol(A)) {
    stop("入力は正方行列である必要があります。")
  }

  # 1. まず、Aの特異値分解(SVD)を行う
  svd_res <- svd(A)
  U_svd <- svd_res$u
  V_svd <- svd_res$v

  # 特異値から対角行列Σを作成
  # svd()が返すdはベクトルですので、対角行列に変換する
  Sigma <- diag(svd_res$d)

  # 2. SVDの結果から極分解のUとPを計算

  U <- U_svd %*% t(V_svd)

  P <- V_svd %*% Sigma %*% t(V_svd)

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

# ===================================================================
# テストと比較
# ===================================================================
# テスト用の行列
A <- matrix(c(
  1, 2, 3,
  4, 5, 6,
  7, 8, 0
), nrow = 3, byrow = TRUE)

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

# --- 実装による分解 ---
cat("\n\n--- 実装による結果 ---\n")
decomp_p <- polar_decomposition(A)
cat("回転行列 U:\n")
print(decomp_p$U)
cat("\n伸縮行列 P:\n")
print(decomp_p$P)

# --- 検算 ---
# 1. UP = A が成り立つか
A_reconstructed <- decomp_p$U %*% decomp_p$P
cat("\n\n--- 検算 ---\n")
cat("U %*% P は A と等しいか:", all.equal(A, A_reconstructed), "\n")

# 2. Uが正規直交行列か
cat("Uは正規直交行列か (t(U)U = I):", all.equal(t(decomp_p$U) %*% decomp_p$U, diag(3)), "\n")

# 3. Pが対称行列か
cat("Pは対称行列か (P = t(P)):", isSymmetric(decomp_p$P), "\n")

# 4. Pが半正定値か (すべての固有値が非負か)
p_eigenvalues <- eigen(decomp_p$P, only.values = TRUE)$values
cat("Pは半正定値か (全固有値 >= 0):", all(p_eigenvalues >= -1e-10), "\n")


cat("\n\n--- (参考)polar{epca}による結果 ---\n")
cat("「The function returns the u matrix」(ヘルプより)\n")
epca::polar(A)
--- 極分解のテスト ---
元の行列 A:
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    0


--- 実装による結果 ---
回転行列 U:
           [,1]       [,2]       [,3]
[1,] -0.6339607 0.67427574  0.3787427
[2,]  0.5243259 0.01474809  0.8513900
[3,]  0.5684859 0.73833239 -0.3628899

伸縮行列 P:
         [,1]     [,2]     [,3]
[1,] 5.442744 5.901595 1.244073
[2,] 5.901595 7.328951 2.111316
[3,] 1.244073 2.111316 6.244568


--- 検算 ---
U %*% P は A と等しいか: TRUE 
Uは正規直交行列か (t(U)U = I): TRUE 
Pは対称行列か (P = t(P)): TRUE 
Pは半正定値か (全固有値 >= 0): TRUE 


--- (参考)polar{epca}による結果 ---
「The function returns the u matrix」(ヘルプより)
           [,1]       [,2]       [,3]
[1,] -0.6339607 0.67427574  0.3787427
[2,]  0.5243259 0.01474809  0.8513900
[3,]  0.5684859 0.73833239 -0.3628899

以上です。