Rで 線形代数:特異値分解 を試みます。
1. 特異値分解 (SVD) の説明
特異値分解 (Singular Value Decomposition, SVD) とは、任意の m x n 行列 A を、以下の3つの行列の積に分解することです。行列 A の形（正方行列かどうか）や性質（対称性など）に制約はありません。
A = U Σ t(V)
それぞれの行列の性質は以下の通りです。
行列 U (左特異ベクトル)
-
m x mの正規直交行列です。 -
Uの各列ベクトルu_iは左特異ベクトルと呼ばれます。 - これらは
Aの「出力空間（列空間）」の正規直交基底を形成します。
行列 Σ (シグマ, 特異値)
-
m x nの対角行列です（Aと同じサイズ）。 - 対角成分
σ_i（Σ_ii）は特異値と呼ばれ、非負の実数です (σ₁ ≥ σ₂ ≥ σ₃ ... ≥ 0)。 - 特異値は、行列
Aがベクトルをどれだけ「引き伸ばす」か、その大きさを表します。ゼロでない特異値の数が、その行列のランクと一致します。
行列 V (右特異ベクトル)
-
n x nの正規直交行列です。 -
Vの各列ベクトルv_iは右特異ベクトルと呼ばれます。 - これらは
Aの「入力空間（行空間）」の正規直交基底を形成します。 - 分解式では転置
t(V)が使われることに注意が必要です。
特異値分解の考え方と求め方
特異値分解は、A とその転置 t(A) から作られる2つの対称行列 t(A)A と A t(A) の固有値分解と関連しています。
-
VとΣを見つける:-
n x nの対称行列t(A)Aを考えます。この行列の固有値分解を行います。 -
t(A)Aの固有ベクトルが、Vの列ベクトル（右特異ベクトル）になります。 -
t(A)Aの固有値λ_iの正の平方根が、特異値σ_i = sqrt(λ_i)になります。
-
-
Uを見つける:- 同様に、
m x mの対称行列A t(A)の固有ベクトルが、Uの列ベクトル（左特異ベクトル）になります。 - または、
A = U Σ t(V)の関係式からAV = UΣと変形し、u_i = (1/σ_i) * A * v_iという関係式を使って、計算済みのVとΣからUを求めることもできます。
- 同様に、
2. R言語による実装
fun_svd <- function(A) {
# 入力チェック
if (!is.matrix(A) || !is.numeric(A)) {
stop("引数 A は数値型の行列でなければなりません。")
}
m <- nrow(A)
n <- ncol(A)
# svd()のデフォルト出力に合わせ、Uはm x p, Vはn x pのpを計算
p <- min(m, n)
# ランクを計算
rank <- qr(A)$rank
# ランク0の場合は特別処理
if (rank == 0) {
return(list(
d = rep(0, p),
u = diag(1, m, p),
v = diag(1, n, p)
))
}
# --- U と特異値 d の計算 ---
AAt <- A %*% t(A)
eigen_AAt <- eigen(AAt, symmetric = TRUE)
# 特異値はp個取得する（ランクを超えた分は0になる）
singular_values_sq <- pmax(eigen_AAt$values, 0)
d <- sqrt(singular_values_sq[1:p])
# U は m x p のサイズで取得
U <- eigen_AAt$vectors[, 1:p, drop = FALSE]
# --- V の計算 ---
# Vは A^T A の固有ベクトルから計算する
AtA <- t(A) %*% A
eigen_AtA <- eigen(AtA, symmetric = TRUE)
# V も n x p のサイズで取得
V <- eigen_AtA$vectors[, 1:p, drop = FALSE]
# --- とVの符号（向き）を合わせる -
# A*v_i = d_i*u_i の関係が成り立つようにVの符号を調整する
# ランクがゼロになる部分（d[i]が小さい）は調整不要
for (i in 1:rank) {
# A %*% V[, i] と d[i] * U[, i] の向きが逆の場合にV側を反転させる
if (sum((A %*% V[, i] - d[i] * U[, i])^2) > sum((A %*% V[, i] + d[i] * U[, i])^2)) {
V[, i] <- -V[, i]
}
}
return(list(d = d, u = U, v = V))
}3.実装コードによる特異値分解の実行例
ランク落ちしていない行列A1
cat("--- 行列A1 ---\n")
(A1 <- matrix(c(1, 2, 3, 4, 5, 7), nrow = 2, ncol = 3))
cat("\n--- 行列A1のランク ---\n")
print(qr(A1)$rank)
cat("\n--- 実装による特異値分解の結果 ---\n")
(result1 <- fun_svd(A1))
cat("\n--- U D V^T を計算して元の行列に戻るか確認 ---\n")
(reconstructed_A1 <- result1$u %*% diag(result1$d) %*% t(result1$v))
cat("\nA1とU D V^Tは等しいか:", all.equal(A1, reconstructed_A1), "\n")
cat("\n--- (参考)Rの組み込み関数による特異値分解の結果 ---\n")
svd(A1)--- 行列A1 ---
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 7
--- 行列A1のランク ---
[1] 2
--- 実装による特異値分解の結果 ---
$d
[1] 10.1914283 0.3671377
$u
[,1] [,2]
[1,] 0.5797531 -0.8147922
[2,] 0.8147922 0.5797531
$v
[,1] [,2]
[1,] 0.2167839 0.93892288
[2,] 0.4904541 -0.34146390
[3,] 0.8440732 -0.04273442
--- U D V^T を計算して元の行列に戻るか確認 ---
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 7
A1とU D V^Tは等しいか: TRUE
--- (参考)Rの組み込み関数による特異値分解の結果 ---
$d
[1] 10.1914283 0.3671377
$u
[,1] [,2]
[1,] -0.5797531 -0.8147922
[2,] -0.8147922 0.5797531
$v
[,1] [,2]
[1,] -0.2167839 0.93892288
[2,] -0.4904541 -0.34146390
[3,] -0.8440732 -0.04273442ランク落ちした行列A2
cat("--- 行列A2 ---\n")
cat("\n--- ランク落ちした行列A2 ---\n")
(A2 <- matrix(c(2, 3, 4, 4, 6, 8), nrow = 3, ncol = 2))
cat("\n--- 行列A2のランク ---\n")
print(qr(A2)$rank)
cat("\n--- 実装による特異値分解の結果 ---\n")
(result2 <- fun_svd(A2))
cat("\n--- U D V^T を計算して元の行列に戻るか確認 ---\n")
(reconstructed_A2 <- result2$u %*% diag(result2$d) %*% t(result2$v))
cat("\nA2とU D V^Tは等しいか:", all.equal(A2, reconstructed_A2), "\n")
cat("\n--- (参考)Rの組み込み関数による特異値分解の結果 ---\n")
svd(A2)--- 行列A2 ---
--- ランク落ちした行列A2 ---
[,1] [,2]
[1,] 2 4
[2,] 3 6
[3,] 4 8
--- 行列A2のランク ---
[1] 1
--- 実装による特異値分解の結果 ---
$d
[1] 1.204159e+01 1.685874e-07
$u
[,1] [,2]
[1,] -0.3713907 0.9284767
[2,] -0.5570860 -0.2228344
[3,] -0.7427814 -0.2971125
$v
[,1] [,2]
[1,] -0.4472136 -0.8944272
[2,] -0.8944272 0.4472136
--- U D V^T を計算して元の行列に戻るか確認 ---
[,1] [,2]
[1,] 2 4
[2,] 3 6
[3,] 4 8
A2とU D V^Tは等しいか: TRUE
--- (参考)Rの組み込み関数による特異値分解の結果 ---
$d
[1] 1.204159e+01 5.617334e-16
$u
[,1] [,2]
[1,] -0.3713907 0.9191450
[2,] -0.5570860 -0.3337319
[3,] -0.7427814 -0.2092736
$v
[,1] [,2]
[1,] -0.4472136 0.8944272
[2,] -0.8944272 -0.4472136縦長の行列A3
set.seed(20250622)
cat("--- 行列A3 ---\n")
(A3 <- matrix(sample(x = 0:100, size = 20 * 3), nrow = 20, ncol = 3))
cat("\n--- 行列A3のランク ---\n")
print(qr(A3)$rank)
cat("\n--- 実装による特異値分解の結果 ---\n")
(result3 <- fun_svd(A3))
cat("\n--- U D V^T を計算して元の行列に戻るか確認 ---\n")
(reconstructed_A3 <- result3$u %*% diag(result3$d) %*% t(result3$v))
cat("\nA3とU D V^Tは等しいか:", all.equal(A3, reconstructed_A3), "\n")
cat("\n--- (参考)Rの組み込み関数による特異値分解の結果 ---\n")
(svd_result3 <- svd(A3))--- 行列A3 ---
[,1] [,2] [,3]
[1,] 48 21 34
[2,] 42 70 45
[3,] 12 79 93
[4,] 26 95 40
[5,] 80 91 3
[6,] 39 65 82
[7,] 83 19 64
[8,] 43 96 86
[9,] 78 27 38
[10,] 28 6 62
[11,] 76 7 17
[12,] 75 81 100
[13,] 29 32 5
[14,] 14 30 22
[15,] 54 67 18
[16,] 4 85 87
[17,] 88 55 77
[18,] 25 60 97
[19,] 33 15 10
[20,] 63 92 36
--- 行列A3のランク ---
[1] 3
--- 実装による特異値分解の結果 ---
$d
[1] 419.8827 139.8720 112.8332
$u
[,1] [,2] [,3]
[1,] -0.13619739 0.146760322 0.128208462
[2,] -0.21889362 -0.005374740 -0.137599653
[3,] -0.26556730 -0.363222079 0.038687859
[4,] -0.23087709 -0.126840188 -0.353243503
[5,] -0.23546480 0.333206890 -0.464173149
[6,] -0.26036600 -0.138869125 0.105295026
[7,] -0.21672520 0.264275024 0.365283259
[8,] -0.31760288 -0.177561550 -0.068858959
[9,] -0.18605027 0.307370981 0.155036246
[10,] -0.12987774 -0.044716867 0.360109709
[11,] -0.12364559 0.397207215 0.161473591
[12,] -0.35224146 -0.004232501 0.156346729
[13,] -0.08941882 0.109350356 -0.139118928
[14,] -0.09298611 -0.035816460 -0.048966729
[15,] -0.19005479 0.162766461 -0.257250628
[16,] -0.25675978 -0.401849879 -0.046841004
[17,] -0.29547355 0.193698931 0.211969256
[18,] -0.25772667 -0.266434033 0.204938846
[19,] -0.07549268 0.144520152 0.006942488
[20,] -0.26394989 0.117711246 -0.303805689
$v
[,1] [,2] [,3]
[1,] -0.4912580 0.8559049 0.1615314
[2,] -0.6349248 -0.2249355 -0.7391039
[3,] -0.5962685 -0.4656510 0.6539366
--- U D V^T を計算して元の行列に戻るか確認 ---
[,1] [,2] [,3]
[1,] 48 21 34
[2,] 42 70 45
[3,] 12 79 93
[4,] 26 95 40
[5,] 80 91 3
[6,] 39 65 82
[7,] 83 19 64
[8,] 43 96 86
[9,] 78 27 38
[10,] 28 6 62
[11,] 76 7 17
[12,] 75 81 100
[13,] 29 32 5
[14,] 14 30 22
[15,] 54 67 18
[16,] 4 85 87
[17,] 88 55 77
[18,] 25 60 97
[19,] 33 15 10
[20,] 63 92 36
A3とU D V^Tは等しいか: TRUE
--- (参考)Rの組み込み関数による特異値分解の結果 ---
$d
[1] 419.8827 139.8720 112.8332
$u
[,1] [,2] [,3]
[1,] -0.13619739 -0.146760322 -0.128208462
[2,] -0.21889362 0.005374740 0.137599653
[3,] -0.26556730 0.363222079 -0.038687859
[4,] -0.23087709 0.126840188 0.353243503
[5,] -0.23546480 -0.333206890 0.464173149
[6,] -0.26036600 0.138869125 -0.105295026
[7,] -0.21672520 -0.264275024 -0.365283259
[8,] -0.31760288 0.177561550 0.068858959
[9,] -0.18605027 -0.307370981 -0.155036246
[10,] -0.12987774 0.044716867 -0.360109709
[11,] -0.12364559 -0.397207215 -0.161473591
[12,] -0.35224146 0.004232501 -0.156346729
[13,] -0.08941882 -0.109350356 0.139118928
[14,] -0.09298611 0.035816460 0.048966729
[15,] -0.19005479 -0.162766461 0.257250628
[16,] -0.25675978 0.401849879 0.046841004
[17,] -0.29547355 -0.193698931 -0.211969256
[18,] -0.25772667 0.266434033 -0.204938846
[19,] -0.07549268 -0.144520152 -0.006942488
[20,] -0.26394989 -0.117711246 0.303805689
$v
[,1] [,2] [,3]
[1,] -0.4912580 -0.8559049 -0.1615314
[2,] -0.6349248 0.2249355 0.7391039
[3,] -0.5962685 0.4656510 -0.6539366以上です。
