Pythonで特異値分解(SVD)を理解する

数学 Python

以前の投稿( 論文メモ: Item2Vec: Neural Item Embedding for Collaborative Filtering - け日記 )で比較対象になっていた特異値分解(SVD)についてまとめ、Pythonで実装してみます。

SVDとは

特異値分解(singular value decomposition: SVD)は、数学的にはM×N(M行N列)の行列を分解する方法の一つです。
コンピュータサイエンスでは、行列で表現される特徴(情報検索の分野では、例えば文書毎の単語の出現頻度の行列など)を可能な限り損ねること無く次元を圧縮するために利用され、多次元の特徴を扱う画像処理、自然言語処理、検索や推薦に応用されています。
他の次元圧縮手法としては、主成分分析(PCA)や非負値行列因子分解(NMF)などがあります。

SVDで何ができるのか

任意のM×Nの行列Cは、以下の形に分解することができます。

C=U \Sigma V^T

  • 左特異行列Uは、M×Mの正方行列で、各列は {CC^T}と直交する固有ベクトル
  • 右特異行列Vは、N×Nの正方行列で、各列は {C^TC}と直交する固有ベクトル
  • ΣはN×M行列で、対角行列成分が {\sigma_1,\dots,\sigma_r} (ただし、 {\sigma_1\geq\dots\geq\sigma_r>0} )で、それ以外が0
    • {\sigma_i}は、行列Cの特異値とよばれる
    • rは行列Cのランク(一次独立な列または行ベクトルの数)

このとき、特異値σは {CC^T}の固有値の平方根となります。

  • {CC^T=U \Sigma V^T V \Sigma^T U^T=U \Sigma\Sigma^T U^T}となる
  • {CC^T}は正方対称行列で、対称対角化定理(正方対称行列Sは {Q \Lambda Q^T}に分解でき、このときのΛはSの固有値になる)により、 {\Sigma\Sigma^T}の対角要素は {CC^T}の固有値 {\lambda_1 \dots \lambda_r}にあたる
  • Σは対角行列のため {\Sigma\Sigma^T=\Sigma^2}であるため、 {\sigma_i=\sqrt{\lambda_i}}となる

具体的な例を挙げます。

以下の3つの文章から単語(ここでは名詞と動詞に限定します)を抽出して、出現頻度行列Cを作ります。

I drive my car.
I read this document.
I read this paper.

行列Cは、行は単語(上から"I", "drive", "car", "read", "document", "paper")、列は文書で、値は単語の出現頻度にそれぞれ対応しています。

$$ C=\left( \begin{matrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \right) $$

これを特異値分解( {C=U \Sigma V^T} )するとこうなります。

  • UのN列目以降、ΣがN行目以降は削除してます(ΣのN行目以降は全て0のため)

$$ C=\left( \begin{matrix} 0.717 & -0.158 & 0 \\ 0.192 & -0.59 & 0 \\ 0.192 & -0.59 & 0 \\ 0.525 & 0.432 & 0 \\ 0.262 & 0.216 & -0.707 \\ 0.262 & 0.216 & 0.707 \end{matrix} \right) \left( \begin{matrix} 2.394 & 0 & 0 \\ 0 & 1.506 & 0 \\ 0 & 0 & 1 \end{matrix} \right) \left( \begin{matrix} 0.46 & 0.628 & 0.628 \\ -0.888 & 0.325 & 0.325 \\ 0 & -0.707 & 0.707 \end{matrix} \right) $$

実践では、データ量が減らせることもあって、 {U\Sigma}と {V^T}の2つの行列で持つことが多いとのことです。

$$ C=\left( \begin{matrix} 1.716 & -0.238 & 0 \\ 0.46 & -0.888 & 0 \\ 0.46 & -0.888 & 0 \\ 1.256 & 0.65 & 0 \\ 0.628 & 0.325 & -0.707 \\ 0.628 & 0.325 & 0.707 \end{matrix} \right) \left( \begin{matrix} 0.46 & 0.628 & 0.628 \\ -0.888 & 0.325 & 0.325 \\ 0 & -0.707 & 0.707 \end{matrix} \right) $$

numpyでSVDを実装するとこんな感じになりますが、numpy.linarg.svdがありますので、実践ではそちらを使いましょう。

import numpy as np
import pandas as pd


np.set_printoptions(precision=3, suppress=True)

# 分解する行列
c = np.array([[1,1,1],[1,0,0,],[1,0,0],[0,1,1],[0,1,0],[0,0,1]])

# 特異値行列Σ(sigma) → 右特異行列V(v) → 左特異行列U(u) の順に求める

# C^TCの固有値と固有ベクトルの計算
ctc = np.dot(c.T, c)
eigen_values, eigen_vectors = np.linalg.eig(ctc)

# 特異値の計算
singular_values = np.sqrt(eigen_values)
singular_index = np.argsort(singular_values)[::-1]

# 特異値行列の計算
sigma = np.diag(singular_values[singular_index])

# 右特異行列の計算
v = eigen_vectors[:,singular_index]

# 左特異行列の計算
u = np.array([np.dot(c, v[:,i]) / sigma.diagonal()[i] for i in range(len(sigma.diagonal()))]).T

SVDを使うとどうして次元圧縮できるのか

{C=U \Sigma V^T}においてΣは対角行列のため、対角要素σが大きければ、Cに与える影響も大きくなります(逆に言いますと、σが小さくなれば影響も小さくなります)。
そこで、上位k(<M)個のσを残し、それ以外は0とした {C_k=U_k \Sigma_k V_k^T}で近似することで、次元を削減します。 低ランク近似と呼ばれており、小さな特異値から削減することで、 {C-C_k}の各要素の2乗誤差が最小になることが知られています。

上の例でk=2として {C_k=U_k \Sigma_k V_k^T}を求めると、以下の形になり、概ね行列Cと {C_k}の値が近いことがわかります。
またこの近似によって、 {C_k}の5行目("document")と6行目("paper")の値はいずれも(0, 0.5, 0.5)となっています。 同一の文書に現れるわけではない"document"と"paper"が、SVDによって類似した単語と判定できるようになった、と言うことができます。
SVDを使ってレコメンドなど欠測値が多いデータの補完に用いられているケースもあります。

$$ C_k=\left( \begin{matrix} 0.717 & -0.158 \\ 0.192 & -0.59 \\ 0.192 & -0.59 \\ 0.525 & 0.432 \\ 0.262 & 0.216 \\ 0.262 & 0.216 \end{matrix} \right) \left( \begin{matrix} 2.394 & 0 \\ 0 & 1.506 \end{matrix} \right) \left( \begin{matrix} 0.46 & 0.628 & 0.628 \\ -0.888 & 0.325 & 0.325 \end{matrix} \right)=\left( \begin{matrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 0.5 & 0.5 \\ 0 & 0.5 & 0.5 \end{matrix} \right) $$

Pythonの実装は以下のとおりです。

# U_kの計算(k列目以降を削除)
u_k = np.delete(u, range(k, u.shape[1]), axis=1)

# Σ_kの計算(k×k行列にする)
sigma_k = np.delete(sigma, range(k, sigma.shape[0]), axis=0)
sigma_k = np.delete(sigma_k, range(k, sigma.shape[1]), axis=1)

# V_kの計算(k列目以降を削除
v_k = np.delete(v, range(k, v.shape[1]), axis=1)

c_k = u_k.dot(sigma_k.dot(v_k.T))

参考文献

理解・実装にあたり、こちら2冊が参考になりました。

岩波データサイエンス Vol.5

岩波データサイエンス Vol.5

情報検索の基礎

情報検索の基礎