対象者標準

対象者別の入口 →

対象者の設定はこのブラウザの中にだけ保存されます。

言語 JA

SVD

中級

02.7.2

SVD

最終更新 2020-02-12 読了時間 5 分

まとめ

特異値分解（SVD）は、どんな形の行列 \(A\) でも「回転 \(U\)・各軸方向の伸縮 \(\Sigma\)・回転 \(V^\top\)」の 3 つの積 \(A = U \Sigma V^\top\) に分ける手法です。伸縮量である特異値 \(\sigma_1 \ge \sigma_2 \ge \dots\) は「データの中で強いパターンから順」に並ぶので、上位 \(k\) 個だけを残した低ランク近似 \(A_k\) は、元の行列を少ない数値で表しながら骨格を保ちます。このページでは 163 × 372 ピクセルの文字画像を scipy.linalg.svd で分解し、rank を 1 から 50 まで変えて復元すると、ぼやけた画像から細部まで戻る様子を確かめます。標準化したデータ行列に SVD をかけると \(V\) の列が PCA の主成分と一致するので、PCA の内部計算としても、推薦システムや LSA の基盤としても使われます。

主成分分析（PCA）の概念を先に学ぶと理解がスムーズです

直感 #

グレースケールの画像は、ピクセルの明るさを並べた 163 × 372 の行列です。この行列を「ぼやけた基本パターン」の重ね合わせとして表し直すのが SVD です。1 枚目のパターン \(u_1 v_1^\top\) は文字のだいたいの位置を表す大きな濡れ跡のような絵で、2 枚目、3 枚目と足していくほど線が細く見えてきます。各パターンにどれだけの重みを掛けるかが特異値 \(\sigma_i\) で、重みの大きい順に 10 枚ほど残せば、163 × 372 = 60,636 個の数値を 10 × (163 + 372 + 1) = 5,360 個で近似できます。

行列 A（163×372） SVD 分解 A = UΣVᵀ 特異値 σ を降順に並べる上位 k 個を選ぶ低ランク近似 A_k = Σ σ_i u_i v_iᵀ 圧縮・ノイズ除去

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：点を軸に射影した散らばりが最大の向きが第 1 主成分。情報を最も失わない 1 本の軸。SVD の第 1 特異ベクトル \(v_1\) は、標準化したデータではこの第 1 主成分と同じ向きになる。

詳細な解説 #

1. 直感：行列を「基本パターン」に分解する #

SVD は、データ行列 \(A\) を「基底ベクトル」と「重み」に分解して表す方法です。画像なら「ぼやけた基本的なパターン」を重ね合わせて元画像を作るイメージで、推薦システムのユーザー × アイテム行列なら「好みのパターン」の重ね合わせになります。特異値が大きい要素から順に「情報量の多いパターン」を表すので、上から何枚残すかが圧縮率を決めます。

2. 数式でみる SVD #

任意の行列 \(A \in \mathbb{R}^{m \times n}\) は、3 つの行列の積に分解できます。

$$ A = U \Sigma V^\top $$

\(U \in \mathbb{R}^{m \times m}\) は左特異ベクトルを列に持つ直交行列で、画像なら「縦方向のパターン」です。\(\Sigma \in \mathbb{R}^{m \times n}\) は特異値を対角成分にもつ対角行列で、各パターンの重みです。\(V \in \mathbb{R}^{n \times n}\) は右特異ベクトルを列に持つ直交行列で、画像なら「横方向のパターン」です。特異値 \(\sigma_1 \ge \sigma_2 \ge \dots\) は「データの強いパターン」順に並ぶので、上位の特異値だけ残すと「低ランク近似」ができます。

3. 実験用のデータ（画像） #

「実験」という文字画像をグレースケール化し、163 × 372 の行列として扱います。

1
2
3
4
5
6
7
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
from PIL import Image

img = Image.open("./sample.png").convert("L").resize((163, 372)).rotate(90, expand=True)
img

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

3-4 行目分解は scipy、画像の読み込みは PIL
scipy.linalg.svd は numpy の np.linalg.svd と同じ分解を返しますが、後で使う linalg.diagsvd（特異値ベクトルから \(\Sigma\) 行列を作る関数）が同じモジュールにあるので scipy を選んでいます。PIL.Image は画像ファイルを読み込み、ピクセルの行列に変換するために使います。
6 行目グレースケール 1 枚を行列に見立てる
convert("L") でカラー 3 チャンネルを明るさ 1 チャンネルにまとめ、1 枚の画像を 1 つの行列にします。resize((163, 372)) と rotate(90, expand=True) で最終的に 163 行 × 372 列の横長の行列になり、この行数と列数が次のブロックで表示する \(U\)（163 × 163）と \(V^\top\)（372 × 372）の形を決めます。./sample.png は手元の任意の文字画像に置き換えて試せます。
7 行目ノートブックで画像をそのまま表示
最後の行に変数名 img だけを書くと、Jupyter などのノートブック環境では画像がそのまま表示されます。分解の前に「元の画像がどう見えるか」を確認しておき、後の rank=1〜50 の復元と見比べる基準にします。

4. 特異値分解を実行 #

画像を numpy 配列にして linalg.svd にかけ、3 つの行列の形を表示します。

1
2
3
4
X = np.asarray(img)
U, Sigma, VT = linalg.svd(X, full_matrices=True)

print(f"A: {X.shape}, U: {U.shape}, Σ:{Sigma.shape}, V^T:{VT.shape}")

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1 行目画像を 163 × 372 の数値行列にする
np.asarray(img) で PIL の画像が 0〜255 の整数を並べた 2 次元配列になります。この配列が数式の \(A\) で、行がピクセルの縦位置、列が横位置、値が明るさです。
2 行目full_matrices=True で正方の U と Vᵀ を得る
linalg.svd は \(U\)・特異値の 1 次元配列・\(V^\top\) の 3 つを返します。full_matrices=True だと \(U\) が 163 × 163、\(V^\top\) が 372 × 372 の正方行列になり、数式の定義と同じ形です。False にすると 372 列のうち使われない 209 列が省かれて \(V^\top\) が 163 × 372 になり、メモリを節約できます。
4 行目特異値は min(m, n) = 163 個
出力の Σ:(163,) は、特異値が行数と列数の小さい方の 163 個しかないことを示します。Sigma は対角行列ではなく特異値だけを並べた 1 次元配列なので、行列として掛けるには次のブロックの diagsvd で対角行列に戻す必要があります。

A: (163, 372), U: (163, 163), Σ:(163,), V^T:(372, 372)

5. 低ランク近似で画像を復元 #

特異値の上位 rank 個だけを残して \(U_k \Sigma_k V_k^\top\) を計算し、rank を 1 から 50 まで変えて復元画像を並べます。

1
2
3
4
5
6
7
8
9
for rank in [1, 2, 3, 4, 5, 10, 20, 50]:
    U_i = U[:, :rank]
    Sigma_i = np.matrix(linalg.diagsvd(Sigma[:rank], rank, rank))
    VT_i = VT[:rank, :]
    temp_image = np.asarray(U_i * Sigma_i * VT_i)

    plt.title(f"rank={rank}")
    plt.imshow(temp_image, cmap="gray")
    plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1 行目rank を 1 から 50 まで振る
[1, 2, 3, 4, 5, 10, 20, 50] の 8 段階で復元します。1〜5 を細かく刻んでいるのは、最初の数枚のパターンで画像の骨格がどれだけ決まるかを見せるためで、50 は 163 個の特異値の 3 割に当たります。
2-4 行目上位 k 本の列・行だけを切り出す
U[:, :rank] は左特異ベクトルの先頭 rank 列、VT[:rank, :] は右特異ベクトルの先頭 rank 行です。linalg.diagsvd(Sigma[:rank], rank, rank) は上位 rank 個の特異値を対角に置いた rank × rank 行列を作ります。この 3 つの積が数式の \(A_k = \sum_{i=1}^k \sigma_i u_i v_i^\top\) で、保持する数値は rank × (163 + 372 + 1) 個です。
5 行目np.matrix の * は行列積
Sigma_i を np.matrix にしているので * が行列積になり、U_i * Sigma_i * VT_i で 163 × 372 の近似画像が得られます。np.ndarray のままだと * は要素ごとの積になって形が合わずエラーになるので、@ 演算子で書くのが現在の推奨です。
7-9 行目グレースケールで復元結果を見る
cmap="gray" で明るさとして表示し、タイトルに rank を入れます。rank=1 ではぼやけた帯、rank=5 で文字の形が読め、rank=20 以上で元画像とほぼ区別がつかなくなります。どの rank で「十分」と感じるかが、この画像に必要な情報量の目安です。

ランクが大きいほど細部が復元され、ランクが小さいほどぼやける。

6. 特異ベクトルの解釈 #

各特異値に対応する \(U, V\) の列は「画像のパターン」を表しています。\(u_1\) は「最も大きな構造」、\(u_2\) は「補助的な構造」です。次のコードは、1 番目から 5 番目までの成分 \(\sigma_i u_i v_i^\top\) を 1 枚ずつ描き、その合計も表示します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
total = np.zeros((163, 372))
for rank in [1, 2, 3, 4, 5]:
    U_i = U[:, :rank]
    Sigma_i = np.matrix(linalg.diagsvd(Sigma[:rank], rank, rank))
    VT_i = VT[:rank, :]

    if rank > 1:
        for ri in range(rank - 1):
            Sigma_i[ri, ri] = 0

    temp_image = np.asarray(U_i * Sigma_i * VT_i)
    total += temp_image

    plt.figure(figsize=(5, 5))
    plt.suptitle(f"$u_{rank}$ の寄与")
    plt.subplot(211)
    plt.imshow(temp_image, cmap="gray")
    plt.subplot(212)
    plt.plot(VT[0])
    plt.show()

plt.imshow(total)

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1-2 行目成分を 1 枚ずつ足し込む土台を用意する
total は 163 × 372 のゼロ行列で、各ランクの成分画像を足し込んでいきます。rank を 1 から 5 まで回すので、最終的に total は rank=5 の低ランク近似と同じ画像になります。
7-9 行目k 番目以外の特異値を 0 にして 1 成分だけ残す
Sigma_i[ri, ri] = 0 で、rank 番目より前の対角成分を 0 にしています。こうすると U_i * Sigma_i * VT_i は \(\sigma_k u_k v_k^\top\) の 1 成分だけになり、「k 枚目のパターンが画像のどこを担当しているか」が見えます。この行を省くと前のブロックと同じ累積の復元画像になってしまいます。
14-20 行目パターンと右特異ベクトルを上下に並べる
上段 subplot(211) が k 番目の成分画像、下段 subplot(212) が VT[0]（第 1 右特異ベクトル）の折れ線です。VT[0] は横方向 372 ピクセルの「明るさの基本パターン」で、文字がある列で値が動きます。本来は VT[rank - 1] を描く方が各成分と対応しますが、ここでは第 1 成分の横パターンを基準として毎回表示しています。
22 行目5 成分の合計を確認する
plt.imshow(total) は 5 成分を足した画像で、前のブロックの rank=5 と一致するはずです。成分を 1 枚ずつ見たあとに合計を見ることで、「低ランク近似 = パターンの重ね合わせ」が目で確かめられます。cmap を指定していないので既定のカラーマップで表示されます。

7. 実務での応用 #

SVD の応用は 4 つの分野で定番です。画像圧縮では、上の例のように rank=20 程度の少数の特異値だけで近似して容量を削減します。ノイズ除去では、小さい特異値はノイズに対応することが多いので、それを 0 にすると画像や信号のざらつきが減ります。推薦システムでは、ユーザー × アイテムの評価行列を分解し、\(U\) と \(V\) の行を「潜在的な好み」のベクトルとして使います。自然言語処理では LSA（潜在意味解析）が単語 × 文書行列を分解します。同じ話題の単語が LSA の座標上で近くに集まるのが利点です。

特異値 \(\sigma_k\) と \(\sigma_{k+1}\) が似た大きさの場合、どこで打ち切るかによって結果が大きく変わります。累積寄与率（上位k個の特異値の二乗和 / 全特異値の二乗和）を確認し、80〜90%を目安にランクを決めてください。また、欠損値（NaN）がある行列にはそのまま適用できないため、事前に補完が必要です。

大規模データにはsklearn.decomposition.TruncatedSVDを使うとメモリ効率よく上位k成分だけを計算できます。テキストデータの潜在意味解析（LSA）では、TF-IDF行列にTruncatedSVDを適用するのが定番です。

まとめ #

SVD は任意の行列を直交行列 \(U\)・特異値の対角行列 \(\Sigma\)・直交行列 \(V^\top\) の積に分解し、特異値の大きい順にデータの強いパターンを取り出します。163 × 372 の文字画像では rank=5 で文字が読め、rank=20 で元とほぼ同じになりました。上位の特異値だけを残す低ランク近似で、圧縮やノイズ除去ができます。標準化したデータ行列に SVD をかけると \(V\) の列が PCA の主成分と一致するので、scikit-learn の PCA も内部では SVD を使っています。打ち切るランクは特異値の 2 乗和の累積寄与率 80〜90% を目安に決め、大規模な疎行列には TruncatedSVD を使います。

発展：SVDと他手法との関係 #

低ランク近似の最適性として、ランク \(k\) の近似 \(A_k\) は次のように表せ、同じランクのどんな行列よりも Frobenius ノルムで元の \(A\) に近い（Eckart–Young の定理）ことが知られています。

$$ A_k = \sum_{i=1}^k \sigma_i u_i v_i^\top $$

PCA との関係では、データ行列を標準化して SVD を行えば \(V\) の列が PCA の主成分と一致し、特異値の 2 乗を \(n\) で割ったものが固有値（各主成分の分散）になります。計算の工夫として、数万 × 数万の大規模データでは全成分を求めるのではなく、ランダム化 SVD や TruncatedSVD で上位 \(k\) 成分だけを効率よく計算します。

主成分分析（PCA） — SVD を活用した次元削減
Kernel-PCA — カーネル行列の固有値分解で非線形構造へ広げる拡張
主成分回帰 — SVD/PCAの回帰応用

ランクと SVD #

スライダーで打ち切りランク n_components を 2, 5, 10, 20, 40 と変えると、手書き数字（10 クラス）の 2D 射影と保持される分散の割合がどう変わるかを確認できます。2 成分で全分散の 24%、5 で 47%、10 で 67%、20 で 85%、40 で 98% を保持し、20 付近で 90% の目安に届きます。