JuliaのSVD：データ分析、機械学習、画像処理への応用

特異値分解とは

与えられた行列 A (m x n) を、以下の3つの行列の積に分解することを指します。

A = U * Σ * V'

ここで、

• V'
  V の複素共役転置 (または転置)。
• V
  n x n のユニタリー行列 (または直交行列)。右特異ベクトルを列ベクトルとして持ちます。
• Σ
  m x n の対角行列。特異値を対角成分に持ち、他の要素は0です。特異値は通常、降順に並んでいます。
• U
  m x m のユニタリー行列 (または直交行列)。左特異ベクトルを列ベクトルとして持ちます。

LinearAlgebra.SVD の使い方

Julia で LinearAlgebra.SVD を使うには、まず LinearAlgebra モジュールを読み込む必要があります。

using LinearAlgebra

そして、分解したい行列を svd 関数に渡します。

A = [1.0 2.0; 3.0 4.0; 5.0 6.0]  # 例の行列
F = svd(A)

svd 関数は、SVD 型のオブジェクトを返します。このオブジェクトには、分解された行列 U, Σ, V が格納されています。

U = F.U  # 左特異ベクトル
S = F.S  # 特異値 (ベクトル)
V = F.V  # 右特異ベクトル
Σ = diagm(S) # 対角行列 Σ を作る場合

F.S は特異値を格納したベクトルを返します。対角行列 Σ を明示的に作りたい場合は、diagm(S) を使います。

LinearAlgebra.SVD のオプション

svd 関数にはいくつかのオプションがあります。

• compute_uv::Bool: デフォルトは true で、UとVを計算します。false にすると、特異値のみを計算します。
• full::Bool: デフォルトは true で、全ての特異ベクトルを計算します。false にすると、計算量を削減できます (特に大きな行列の場合)。

例：

F = svd(A; full=false) # 特異値のみを計算
F = svd(A; compute_uv=false) # 特異値のみを計算

特異値分解の応用

特異値分解は、以下のような様々な応用があります。

• 推薦システム
  協調フィルタリングなどに使われます。
• 画像処理
  画像の圧縮やノイズ除去などに使われます。
• 線形方程式の求解
  特異値分解を使って、線形方程式の最小二乗解を求められます。
• 行列の近似
  特異値の小さいものを無視することで、元の行列を低ランクの行列で近似できます。
• 次元削減
  特異値の大きいものからいくつかを選び、対応する特異ベクトルを使って低次元の近似空間を構成します (主成分分析など)。

---

一般的なエラーと原因

1. ArgumentError: A must be a real or complex matrix: 入力行列 A が実数または複素数行列でない場合に発生します。行列の要素がシンボリック型や他の非数値型になっていないか確認してください。

2. DimensionMismatch: 特異値分解の結果を利用する際に、行列の次元が一致しない場合に発生します。例えば、U、Σ、V の次元が期待するものと異なる場合によく起こります。特に、full=false オプションを使用した場合に、U や V のサイズが小さくなることに注意が必要です。

3. SingularException: 行列が特異行列 (正方行列で逆行列を持たない) または特異に近い行列の場合に発生することがあります。このような行列では、特異値分解の結果が数値的に不安定になる可能性があります。

4. DomainError: 入力行列の要素が不正な値 (NaN, Inf など) を含んでいる場合に発生することがあります。

トラブルシューティング

1. 入力行列の確認
   まず、入力行列 A の型、次元、要素の値を確認します。typeof(A)、size(A)、display(A) などを使って、行列が意図した通りに作成されているか確認してください。特に、数値型であることを確認することが重要です。

2. 特異値の確認
   svd(A) の結果として得られる特異値 F.S を確認します。特異値が極端に小さい値を含んでいる場合、行列が特異に近い可能性があります。minimum(F.S) や maximum(F.S) などを使って、特異値の分布を確認すると良いでしょう。

3. full オプションの確認
   full=false オプションを使用している場合、U や V のサイズが小さくなることに注意してください。必要な部分行列だけが計算されるため、次元が合わないというエラーが発生しやすくなります。必要に応じて full=true に変更するか、U や V の適切な部分を取り出すようにしてください。

4. compute_uv オプションの確認
   特異値だけが必要で U と V が不要な場合は、compute_uv=false を指定することで計算量を削減できます。しかし、U や V を後で使おうとするとエラーが発生するので、注意が必要です。

5. 数値的安定性
   行列が特異に近い場合、特異値分解の結果が数値的に不安定になることがあります。このような場合は、特異値を小さい値でクリップしたり、他の分解方法を検討したりすることを推奨します。

6. エラーメッセージの確認
   Julia が表示するエラーメッセージをよく読んでください。エラーメッセージには、問題の原因や場所に関するヒントが含まれていることが多いです。

7. デバッグツールの利用
   Julia のデバッガー (@enter マクロなど) を利用して、コードの実行過程をステップバイステップで確認することができます。これにより、エラーが発生した箇所や変数の値を特定しやすくなります。

8. 最小限のコードで再現
   問題を報告する場合、できるだけ最小限のコードで問題を再現できるようにしてください。これにより、他の人が問題を理解しやすくなり、解決策を見つけやすくなります。

例: DimensionMismatch の場合

using LinearAlgebra

A = [1.0 2.0; 3.0 4.0; 5.0 6.0]
F = svd(A; full=false)  # U, V が小さいサイズになる

# U * Σ * V' を計算しようとするとエラーが発生する
# なぜなら、U, Σ, V' のサイズが合わないから
# A ≈ F.U * diagm(F.S) * F.V'  # これはエラーになる

# 正しい計算例
A ≈ F.U * diagm(F.S) * transpose(F.V) # transposeを使う必要がある

F = svd(A; full=true) # full=trueの場合
A ≈ F.U * diagm(F.S) * F.V' # これは正しく計算できる

---

例1: 基本的なSVD分解と再構成

using LinearAlgebra

# 例の行列
A = [1.0 2.0; 3.0 4.0; 5.0 6.0]

# SVD分解
F = svd(A)

# 分解された行列を取得
U = F.U
S = F.S  # 特異値
V = F.V

# 対角行列 Σ を作成
Σ = diagm(S)

# 元の行列を再構成 (近似)
A_reconstructed = U * Σ * V'

# 再構成誤差を確認
println("再構成誤差: ", norm(A - A_reconstructed))

# 結果の確認
println("U:\n", U)
println("Σ:\n", Σ)
println("V:\n", V)
println("A_reconstructed:\n", A_reconstructed)

この例では、まず行列 A を定義し、svd(A) で特異値分解を行っています。分解された U、S (特異値)、V を使って、元の行列を再構成しています。norm(A - A_reconstructed) で再構成誤差を計算し、元の行列と再構成された行列がどれくらい近いかを確認できます。

例2: 特異値を使った次元削減 (主成分分析の例)

using LinearAlgebra

# データ行列 (各行がサンプル、各列が特徴量)
X = randn(100, 5)

# 中心化 (各特徴量の平均を0にする)
X_centered = X .- mean(X, dims=1)

# SVD分解
F = svd(X_centered)

# 上位2つの特異値に対応する特異ベクトルを取得
k = 2  # 削減後の次元数
U_reduced = F.U[:, 1:k]

# 低次元空間への射影
X_reduced = X_centered * U_reduced

# 結果の確認
println("元のデータ行列のサイズ: ", size(X))
println("削減後のデータ行列のサイズ: ", size(X_reduced))

# X_reducedを使って、例えば可視化などを行う

この例では、ランダムなデータ行列 X を作成し、主成分分析 (PCA) を行って次元削減を行っています。svd(X_centered) で中心化されたデータ行列を特異値分解し、上位 k 個の特異ベクトルを使って低次元空間へ射影しています。

例3: 特異値分解を用いた線形方程式の最小二乗解

using LinearAlgebra

# 行列 A とベクトル b
A = [1.0 2.0; 3.0 4.0; 5.0 6.0]
b = [3.0; 7.0; 11.0]

# SVD分解
F = svd(A)

# 最小二乗解を計算
x = F.V * diagm(1 ./ F.S) * F.U' * b

# 結果の確認
println("最小二乗解: ", x)

# 解の検証
println("Ax - b: ", A * x - b) # 残差が小さければOK

この例では、特異値分解を使って線形方程式 Ax = b の最小二乗解を求めています。特異値分解の結果を使って、最小二乗解を直接計算することができます。

using LinearAlgebra

A = randn(100, 10) # 大きな行列

# full=false で計算
F_reduced = svd(A; full=false)

# UとVのサイズが小さくなっていることを確認
println("Uのサイズ: ", size(F_reduced.U))
println("Vのサイズ: ", size(F_reduced.V))
println("Sのサイズ: ", size(F_reduced.S))

# full=true で計算
F_full = svd(A; full=true)

# UとVのサイズが大きくなっていることを確認
println("Uのサイズ: ", size(F_full.U))
println("Vのサイズ: ", size(F_full.V))
println("Sのサイズ: ", size(F_full.S))

# 時間計測
@time svd(A; full=false)
@time svd(A; full=true)

---

特異値分解 (SVD) の計算方法の選択

Julia の svd 関数は、様々なアルゴリズムを内部で使用して SVD を計算します。LinearAlgebra パッケージには、いくつかの関連関数も存在します。

• bidiagonalize(A)
  A を二重対角行列に変換します。SVD の計算の前処理として使われることがあります。
• svdfact(A)
  特異値分解の因子を計算します。svd と似ていますが、因子の型が異なります。必要に応じて、U, S, V を取り出すことができます。
• svdvals(A)
  特異値のみを計算する場合に使用します。svd(A; compute_uv=false) と同じ結果になりますが、より直接的な関数です。U や V を計算する必要がない場合に高速です。
• svd(A)
  通常の SVD。多くの場合はこれで十分です。

特異値分解の近似計算

大規模な行列の場合、完全な SVD 計算は計算コストが高くなることがあります。近似的な SVD を計算する方法も存在します。

• ランダム化 SVD
  特異値分解を高速に近似計算する方法です。Arpack.jlなどのパッケージで利用できます。大規模なデータに対して有効です。

using Arpack

# 例: ランダム化 SVD (k は近似ランク)
k = 10
F = svds(A; nsv=k) # nsvは取得したい特異値の数
U = F.U
S = F.S
V = F.V

# または
F = svds(A, k) # nsv=kの省略形

• Interpolative Decomposition (ID)
  行列の低ランク近似を効率的に計算する方法です。これも大規模な行列に対して有効です。

固有値分解 (EVD) の利用

行列が正方行列で、かつ対称行列（またはエルミート行列）の場合、特異値分解の代わりに固有値分解 (EVD) を使うことができます。EVD は SVD よりも計算コストが低い場合があります。

• eigvecs(A)
  正方行列 A の固有ベクトルのみを計算します。
• eigvals(A)
  正方行列 A の固有値のみを計算します。
• eigen(A)
  正方行列 A の固有値と固有ベクトルを計算します。

対称行列 A に対しては、svd(A) の代わりに eigen(A) を使うと、U と V が同じになり（符号の違いはありえます）、特異値は固有値の絶対値になります。

using LinearAlgebra

A = [1.0 2.0; 2.0 1.0] # 対称行列

# EVD
F_eig = eigen(A)
λ = F_eig.values # 固有値
v = F_eig.vectors # 固有ベクトル

# SVD
F_svd = svd(A)
σ = F_svd.S # 特異値
U = F_svd.U
V = F_svd.V

println("EVD eigenvalues: ", λ)
println("SVD singular values: ", σ)
println("EVD eigenvectors: ", v)
println("SVD left singular vectors: ", U)
println("SVD right singular vectors: ", V)

# 対称行列では、固有値の絶対値は特異値に等しい
@assert abs.(λ) ≈ σ
# UとVは（符号を除いて）等しい
@assert abs.(U) ≈ abs.(V)

特異値分解の必要性の再検討

そもそも、特異値分解が本当に必要かどうかを再検討することも重要です。問題によっては、他の方法でより効率的に解決できる場合があります。例えば、線形方程式を解く場合、svd ではなく、LU分解やQR分解の方が適している場合があります。