特異値分解を詳しく解説

はじめに

特異値分解は機械学習ではよく使われるテクニックだが、大学の教養過程で使うような線形代数の教科書には載っていなかったりする。
書店で色々探した結果、こちらの本に分かりやすい解説があり、理解が深まったのでまとめておく。
線形代数セミナー: 射影,特異値分解,一般逆行列　金谷 健一 (著)

今回は特異値分解の解説のみで、具体的な活用法はまた別の機会に書く。

定義

${\rm rank}(\boldsymbol A)=r$ の $m \times n$ 行列 $\boldsymbol A$ を考える。 ただし、$r < m < n$とする。
$\rm rank$って何？という方はあまり気にしなくても大丈夫ですが、気になる方はこちらの記事を参照。
このとき、$\boldsymbol A$の特異値分解は下記のように表される。

$$\boldsymbol A = \boldsymbol U \boldsymbol \Gamma \boldsymbol V^T \tag{a}$$

ここで、$\boldsymbol U, \boldsymbol V$は直交行列であり、$\boldsymbol \Gamma$は対角行列である。$\boldsymbol \Gamma$の対角成分$\sigma_1,..., \sigma_r$が行列$\boldsymbol A$の特異値と呼ばれる。特異値は0より大きい値である（0も特異値とする場合もある）。

ここで、$\boldsymbol U = (\boldsymbol u_1 \cdots \boldsymbol u_r), \boldsymbol V = (\boldsymbol v_1 \cdots \boldsymbol v_r)$とすると、$\boldsymbol u_i, \boldsymbol v_i$ （$\boldsymbol U, \boldsymbol V$の列ベクトル）はそれぞれ左特異ベクトル、右特異ベクトルと呼ぶ。$\boldsymbol U, \boldsymbol V$がそれぞれ$\boldsymbol \Gamma$の左、右に来るためだと覚えておけばよい。

特異値・特異ベクトルの定義は下記のようにも書ける。固有値・固有ベクトルの定義（$\boldsymbol A \boldsymbol v = \lambda \boldsymbol v$）に似ているが、左辺と右辺に含まれるベクトルがそれぞれ $\boldsymbol u, \boldsymbol v$ と異なるベクトルになっていることに注意。

\begin{align}
\boldsymbol A \boldsymbol v &= \sigma \boldsymbol u \\
\boldsymbol A^T \boldsymbol u &= \sigma \boldsymbol v \tag{b} \\
(\sigma > 0, \boldsymbol u &\neq \boldsymbol 0, \boldsymbol v \neq \boldsymbol 0)
\end{align}

ここで、特異値と固有値の関係を知るために、少し式変形してみる。
(b)式にそれぞれ、$\boldsymbol A^T, \boldsymbol A$を左からかけると、

\boldsymbol A^T \boldsymbol A \boldsymbol v = \sigma \boldsymbol A^T \boldsymbol u = \sigma^2 \boldsymbol v \\
\boldsymbol A \boldsymbol A^T \boldsymbol u = \sigma \boldsymbol A \boldsymbol v = \sigma^2 \boldsymbol u

上式の最左辺と最右辺を比べると、前述した固有ベクトルの定義（$\boldsymbol A \boldsymbol v = \lambda \boldsymbol v$）と同じ形になっている。
よって、$\boldsymbol u, \boldsymbol v$はそれぞれ、$\boldsymbol A \boldsymbol A^T, \boldsymbol A^T \boldsymbol A$ の固有ベクトルであり、固有値はいずれも$\sigma^2$となる。

つまり、

ということが言える。

特異値分解のイメージとしては、対称行列（行と列の数が等しい行列）で定義される固有値を、一般の行列に拡張したものだと考えると良さそう。実際、対称行列の場合は固有値と特異値は一致する。

（参考）

ここで、前述の(a)式と(b)式が等価であることを示しておく。まず、(a)式は、以下のように変形できる。

\begin{align}
\boldsymbol A &= \boldsymbol U \boldsymbol \Gamma \boldsymbol V^T \\
&= \left(\boldsymbol u_1 \cdots \boldsymbol u_r \right) 
\left(
    \begin{array}{cccc}
      \sigma_1 & & \\
      & \ddots & \\
      & & \sigma_r
    \end{array}
\right)
\left(
    \begin{array}{cccc}
      \boldsymbol v_1^T \\
      \vdots \\
      \boldsymbol v_r^T
    \end{array}
\right) \\
&= \sigma_1 \boldsymbol u_1 \boldsymbol v_1^T + \cdots + \sigma_r \boldsymbol u_r \boldsymbol v_r^T \tag{c}
\end{align}

さらに、前述の通り、$\{\boldsymbol u_1 \cdots \boldsymbol u_r\}, \{\boldsymbol v_1 \cdots \boldsymbol v_r\}$は対象行列の固有ベクトルなので、正規直交基底として選ぶことができる（参考リンク1、参考リンク2）。すなわち、$\boldsymbol u_i^T \boldsymbol u_j = \delta_{ij}, \boldsymbol v_i^T \boldsymbol v_j = \delta_{ij}$ （ただし$\delta_{ij}$はクロネッカーのデルタと呼ばれ、$i=j$の場合に1, $i \neq j$の場合に0となる）である。

この性質を踏まえ、(c)式の両辺に右から例えば$\boldsymbol v_1$をかけると、

\begin{align}
\boldsymbol A \boldsymbol v_1 &= \sigma_1 \boldsymbol u_1 \boldsymbol v_1^T \boldsymbol v_1 + \cdots + \sigma_r \boldsymbol u_r \boldsymbol v_r^T \boldsymbol v_1 \\
&= \sigma_1 \boldsymbol u_1
\end{align}

となり、(b)の1番目の式と一致する。$\boldsymbol v_2, \dots, \boldsymbol v_r$についても同様。
また、(c)式を転置した、
$$\boldsymbol A^T = \sigma_1 \boldsymbol v_1 \boldsymbol u_1^T + \cdots + \sigma_r \boldsymbol v_r \boldsymbol u_r^T$$

に対して右から$\{\boldsymbol u_1 \cdots \boldsymbol u_r\}$をかけることで、(b)の2番目の式と一致することも示せる。

表記法

特異値分解の表記方法はいくつかあって、自分はちょっと混乱した。

まず、０を特異値に含めて、$\boldsymbol \Gamma$を $m \times m$ の行列とする場合もある（下図(2)）。ただし、 $\mathbf u, \mathbf v$ は特異値０に対応する左特異ベクトル、右特異ベクトルを並べた行列である。表記(1)から拡張されたグレーの部分は、結局かけ合わせて０になるので、実質的な影響はない。

また、$\boldsymbol \Gamma$を対角行列ではなく、$m \times n$の行列とした表記もある（下図(3)）。これは、$\boldsymbol \Gamma$を 零行列でさらに拡張して$m \times n$に、$V$も$n \times n$に拡張したものだと考えれば良い。やはり、グレーの部分は掛け合わせて０になるので、実質的な影響はない。

Pythonコード

pythonでは、numpyのlinalgモジュールにあるsvdという関数で特異値分解を実行することができる。
この関数が返す $\boldsymbol U, \boldsymbol V$ はデフォルトでは表記(3)に対応する対称行列である。ただし、full_matrices=Falseとすると、表記(2)に対応する $\boldsymbol U, \boldsymbol V$ が得られる。
また、$\boldsymbol \Gamma$については、行列ではなく（０を含む）特異値の配列が返される。行列$\boldsymbol A$のランクは２なので、０でない特異値の数は２つとなっている。

特異値分解の$\boldsymbol \Gamma$にするためには、np.diagで対角行列にする必要がある。
確認のため、$\boldsymbol U, \boldsymbol \Gamma, \boldsymbol V$の積をとってみると、確かに元の行列$\boldsymbol A$と一致することが分かる。

import numpy as np
from numpy.linalg import svd, matrix_rank

# rank=2の　３ｘ4行列を作成
A = np.array([[2, 4, 1, 3], [1, 5, 3, 2], [5, 7, 0, 7]])

print('matrix A\n', A)
print('rank: ', matrix_rank(A))

# singular value decomposition
u, s, vh = svd(A)
print('\nSVD result')
print('shape of u, s, vh:', u.shape, s.shape, vh.shape)
print('singular values:', s.round(2))

# full_matrices=Falseの場合
u, s, vh = svd(A, full_matrices=False)
print('\nSVD result　(full_matrices: False)')
print('shape of u, s, vh:', u.shape, s.shape, vh.shape)

# 復元
# A_re = (u @ np.diag(s, -1)[1:] @ vh).round(2)
A_re = (u @ np.diag(s) @ vh).round(2)
print('\nreconstructed A:\n', A_re)

matrix A
 [[2 4 1 3]
 [1 5 3 2]
 [5 7 0 7]]
rank:  2

SVD result
shape of u, s, vh: (3, 3) (3,) (4, 4)
singular values: [13.39  3.58  0.  ]

SVD result　(full_matrices: False)
shape of u, s, vh: (3, 3) (3,) (3, 4)

reconstructed A:
 [[ 2.  4.  1.  3.]
 [ 1.  5.  3.  2.]
 [ 5.  7. -0.  7.]]