特異値分解と行列の低ランク近似

はじめまして。テクノロジーインテリジェンス部で主にデータ周りの業務を担当しているshmです。今回、並河さんからご指名を頂き、記事を執筆させて頂くことになりました。

さてデータといえば、最近では統計処理ソフトR、Pythonの数値計算ライブラリNumPyや機械学習ライブラリscikit-learnなど高度なデータ分析を行うための環境を簡単に（しかもフリーで！）手に入れることが出来るようになりました。以前であれば大学の研究者やアナリストといった専門家向けのマニアックなツールとして知る人ぞ知る存在だったはずですが、近年のブームの中、書籍やネットで紹介されることで、その存在が広く知れ渡るようになりました。これらのツールは一定の書式で値を入力するだけで簡単に解析結果を返してくれる大変便利なものですが、得られた結果が妥当なものであるかどうかは別途、評価しなければならないのが普通です。その際、適切な評価を行う為にもツールが前提としている理論について把握しておくことが望ましいのですが、実際にはブラックボックスとして使われることも多いと思われます。

そこで本記事では、こうしたツールが提供する様々な機能の中から、広い応用範囲を持つ特異値分解にスポットを当てて解説したいと思います。

本記事の構成は以下の通りです。

ベクトルと行列
固有値分解
特異値分解
行列の低ランク近似
最後に

はじめにベクトルや行列について、本記事の内容を理解するために必要なことに絞って簡単に復習します。次に特異値分解を理解するために必要となる固有値分解を紹介した後*1、特異値分解について証明付きで解説します*2。最後に応用として、行列の低ランク近似について取り上げます。行列の低ランク近似は、データを要約して次元を圧縮するためにしばしば使われますが、この記事では行列が近似されていく様子を画像データを使い視覚的に確認します。

ベクトルと行列

以下、 $\mathbf{K}$ は複素数全体 $\mathbf{C}$ または実数全体 $\mathbf{R}$ を表すものとします。

$\mathbf{K}$ の元 $a_{i}$ を $(a_{1},\ldots,a_{n})$ の様に並べて出来た集合を $\mathbf{K}^{n}$ とします。 $\mathbf{K}^{n}$ に和とスカラー倍 \begin{gather*}(a_{1},\ldots,a_{n})+(b_{1},\ldots,b_{n}) :=(a_{1}+b_{1},\ldots,a_{n}+b_{n})\\\lambda(a_{1},\ldots,a_{n}):=(\lambda a_{1},\ldots,\lambda a_{n})\end{gather*} を定義して線型空間の構造を与えるとき、 $\mathbf{K}^{n}$ の各要素 $(a_{1},\ldots,a_{n})$ をベクトルと呼びます。

次に、 $mn$ 個の $\mathbf{K}$ の要素 $a_{ij}$ $(i=1,2,\ldots,m,\,j=1,2,\ldots,n)$ を縦横２次元に並べたものを $m\times n$ 行列と呼び、 \[\left(\begin{array}{ccc}a_{11} & \cdots & a_{1n}\\\vdots & \ddots & \vdots\\a_{m1} & \cdots & a_{mn}\end{array}\right)\] と書きます。これをより簡潔に $(a_{ij})$ と書くことがあります。さらに \[\left(\begin{array}{cc}A_{p,q} & A_{p,n-q}\\A_{m-p,q} & A_{m-p,n-q}\end{array}\right)\] （ $A_{r,s}$ は $r\times s$ 行列を表す）の様にブロックに分けて書くことがあります。

良く目にするタイプの行列には特に名前が付けられています。

定義

$n\times n$ 行列を $n$ 次の正方行列と呼ぶ。
対角要素以外 0 の行列を対角行列と呼ぶ。対角行列の内、対角要素が全て 1 であるものを単位行列と呼び $I$ と書く。
$m \times n$ 行列 $A=(a_{ij})$ の行と列を入れ替えて出来る $n\times m$ 行列 $({a}_{ji})$ を転置行列と呼び ${}^{t}\! A$ と書く。転置した上で複素共役を取った行列 ${}^{t}\! \bar{A}=(\bar{a}_{ji})$ を随伴行列と呼び $A^{*}$ と書く。
正方行列 $A$ が $A^{*}=A$ を満たすとき、 $A$ をエルミート行列と呼ぶ。
正方行列 $U$ が $U^{*}U=UU^{*}=I$ を満たすとき、 $U$ をユニタリ行列と呼ぶ。
正方行列 $P$ が $P^{2}=P=P^{*}$ を満たすとき $P$ を射影行列と呼ぶ。

行列には以下の様に和と積が定義されます： $(a_{ij})$ 、 $(b_{ij})$ を $m\times n$ 行列、 $(c_{ij})$ 、 $(d_{ij})$ をそれぞれ $p\times r$ 、 $r\times q$ 行列として \begin{gather*}(a_{ij})+(b_{ij}):=(a_{ij}+b_{ij}),\\(c_{ij})(d_{ij}):=\left(\sum_{k=1}^{r}c_{ik}d_{kj}\right).\end{gather*}

行列の内、 $m\times1$ 行列や $1\times n$ 行列はしばしばベクトルと同一視され、それぞれ $m$ 次の列ベクトル、 $n$ 次の行ベクトルと呼ばれます。

なお $a:=(a_{1},\ldots,a_{n})$ 、 $b=(b_{1},\ldots,b_{n})$ をベクトルとするとき、\[(a,b)=\sum_{i=1}^{n}\bar{a}_{i}b_{i}\]によって標準内積と呼ばれる演算が定義されます*3が、 $a$ および $b$ を列ベクトルと同一視するとき\[(a,b)={}^{t}\bar{a}b\]（右辺は行列としての積）が成り立ちます。また $\sqrt{(a,a)}$ を $\|a\|$ と書くことにします。

さて線型空間 $\mathbf{K}^{n}$ のどの要素も $\{e_{i}\}_{i=1}^{n}$ 、 $e_{i}\in\mathbf{K}^{n}$ に関する和とスカラー倍で書き表すことができるとき $\{e_{i}\}_{i=1}^{n}$ を基底と呼びます。特に $(e_{i},e_{j})=\delta_{ij}$ が成り立つとき $\{e_{i}\}_{i=1}^{n}$ を正規直交基底と呼びます。

固有値分解

正方行列 $A$ に対して、\[Ax=\lambda x\] を満たすスカラー $\lambda$ と列ベクトル $x\neq0$ が存在するとき、 $\lambda$ を $A$ の固有値、 $x$ を $A$ の固有ベクトルと呼びます。また $\lambda$ の固有ベクトル全体から作られる線型空間を $\lambda$ に対する固有空間と呼びます。

特異値分解を理解する上で、エルミート行列に対する固有値分解が鍵となるので結果を紹介します。

定理１（固有値分解） $n$ 次エルミート行列 $A$ に対して、適切なユニタリ行列 $U$ を選んで $A$ を対角化することが出来る： \[U^{*}AU=\left(\begin{array}{ccc}\lambda_{1} & & O\\ & \ddots\\O & & \lambda_{n}\end{array}\right).\] また $U$ を $n$ 次元列ベクトル $u_{i}$ を用いて $U=(u_{1} \ldots u_{n})$ と表すと、上式は \[A=\sum_{i=1}^{n}\lambda_{i}u_{i}{}^{t}\bar{u}_{i}\] と書き表すことが出来る。 $U$ がユニタリ行列より $\{u_{i}\}_{i=1}^{n}$ は線型空間 $\mathbf{K}^{n}$ の正規直交基底である。

固有値分解に現れる $n$ 次正方行列 $u_{i}{}^{t}\bar{u}_{i}$ は $A$ の固有ベクトル $u_{i}$ への射影行列になっています。実際、任意の $x\in\mathbf{K}^{n}$ に対して\[(u_{i}{}^{t}\bar{u}_{i})x=(u_{i},x)u_{i}\]であることから直接計算によって確かめることが出来ます。

余談ですが、行列の固有値は $n$ 次の代数方程式\[\det(\lambda I-A)=0\]（固有方程式または特性方程式と呼ばれています）の解であることから、固有値を求める為に固有方程式を解くことを考えるかもしれません*4。しかしそれが実行可能なのは $n=4$ までで、 $n$ が 5 以上になると方程式の係数に対する四則演算や根号の組合せによって解を求めることが（一般には）出来なくなります*5。当然、計算機を使って固有値を求める場合にもこの制約は付いて回るため、一般に固有値を求めるアルゴリズムは反復解法に限られることになります。

特異値分解

固有値分解は正方行列に対する分解ですが、全ての正方行列に対して固有値分解が可能であるわけではありません *6。一方、特異値分解はどのような行列に対しても実行可能です。

定理２（特異値分解） $A$ を $m$ 行 $n$ 列の行列とする。このとき $m$ 次のユニタリ行列 $U$ と $n$ 次のユニタリ行列 $V$ が存在して \[U^{*}AV=D:=\left(\begin{array}{cccc}\sigma_{1} & & O\\ & \ddots & & O_{r,n-r}\\O & & \sigma_{r}\\ & O_{m-r,r} & & O_{m-r,n-r}\end{array}\right)\]と書くことができる。ただし \[\sigma_{1}\ge\sigma_{2}\ge\ldots\ge\sigma_{r}>\sigma_{r+1}=\ldots=\sigma_{n}=0.\] $\sigma_{i}$ を $A$ の特異値と呼ぶ。また $U$ 、 $V$ を $m$ 次元列ベクトル $u_{i}$ 、 $n$ 次元列ベクトル $v_{i}$ を用いて $U=(u_{1} \dots u_{m})$ 、 $V=(v_{1} \dots v_{n})$ と表すと、上式は\[A=\sum_{i=1}^{\min\{m,n\}}\sigma_{i}u_{i}{}^{t}\bar{v}_{i}\]と書き表すことが出来る。ここで $U$ 、 $V$ がユニタリ行列より $\{u_{i}\}_{i=1}^{m}$ 、 $\{v_{i}\}_{i=1}^{n}$ は線型空間 $\mathbf{K}^{m}$ 、 $\mathbf{K}^{n}$ の正規直交基底である。

証明 $A^{*}A$ はエルミート行列であるから対角化可能である。 $A^{*}A$ の任意の固有値を $\lambda$ 、それに対する長さ１の固有ベクトルを $x$ とすると \[\lambda=\lambda\|x\|^{2}=\left(x,A^{*}Ax\right)=\|Ax\|^{2}\ge0\]より $A^{*}A$ の固有値は非負であるから、 $A^{*}A$ の固有値分解は \[A^{*}A=\sum_{i=1}^{n}\sigma_{i}^{2}v_{i}{}^{t}\bar{v}_{i}\]と書ける。ただし $\sigma_{i}\ge0$ は $\sigma_{i}\ge\sigma_{i+1}$ となるように並べるものとする。
さて、 $u_{i}$ を \[u_{i}:=\frac{1}{\sigma_{i}}Av_{i},\quad i=1,2,\ldots,r\]で定義すると、 $\{u_{i}\}_{i=1}^{r}$ は\[(u_{i},u_{j})=\frac{1}{\sigma_{i}\sigma_{j}}\left(v_{i},A^{*}Av_{j}\right)=\delta_{ij}\] を満たす、つまり長さが１で互いに直交しているベクトルの集合である。さらに $\{u_{i}\}_{i=r+1}^{m}$ を適当に定め、 $\{u_{i}\}_{i=1}^{m}$ が $\mathbf{K}^{m}$ の正規直交基底になるようにする。
$u_{i}$ ( $i=1,2,\ldots,r$ ) の定義、および $\sigma_{i}=0$ ( $i=r+1,\ldots,\min\{m,n\}$ )から \[\sum_{i=1}^{\min\{m,n\}}\sigma_{i}u_{i}{}^{t}\bar{v}_{i} =\sum_{i=1}^{r}\sigma_{i}u_{i}{}^{t}\bar{v}_{i}=A\left(\sum_{i=1}^{r}v_{i}{}^{t}\bar{v}_{i}\right)=AP\]が成り立つ。ここで $P:=\sum_{i=1}^{r}v_{i}{}^{t}\bar{v}_{i}$ とした。
一方、任意の $x\in\mathbf{K}^{n}$ に対して \[\begin{split}A^{*}A(1-P)x & =\left(\sum_{i=1}^{n}\sigma_{i}^{2}v_{i}{}^{t}\bar{v}_{i}\right)\left(\sum_{j=r+1}^{n}v_{j}{}^{t}\bar{v}_{j}\right)x\\ & =\sum_{i=r+1}^{n}\sigma_{i}^{2}\left(v_{i},x\right)v_{i}\\ & =0\end{split}\] より\[\|A(I-P)x\|^{2}=\left((I-P)x,A^{*}A(1-P)x\right)=0,\]つまり $A(I-P)=0$ が成り立つ。
以上より\[A=AP+A(I-P)=\sum_{i=1}^{\min\{m,n\}}\sigma_{i} u_{i}{}^{t}\bar{v}_{i}\]となる。（証明終）

行列の低ランク近似

特異値分解の応用として、行列の低ランク近似を取り上げます。

行列 $A$ の特異値分解\[A=\sum_{i=1}^{\min\{m,n\}}\sigma_{i}u_{i}{}^{t}\bar{v}_{i}\]において、途中の $k$ 番目まで和をとって出来る行列 \[A_{k}=\sum_{i=1}^{k}\sigma_{i}u_{i}{}^{t}\bar{v}_{i},\quad k=1,2,\ldots,\min\{m,n\}\] を考えてみましょう。実は、こうして作られる $A_{k}$ が行列 $A$ の持つ情報を要約したもの（低ランク近似）になっているのですが、このことを画像データを使って視覚的に確認してみたいと思います。

まず次の画像（750 $\times$ 1125ピクセル、8ビットグレイスケール）の輝度を画像処理ソフトを使ってテキストデータ（000.csv）として出力します。

f:id:astamuse:20170613111922j:plain:w400

画像の輝度が作る行列を $A$ とします。 $A$ は 750 行 1125 列の行列であり、 $A$ の各成分は 0 から 255 までの整数となっています。この行列 $A$ を特異値分解して行列 $A_{k}$ を計算します。 $A_{k}$ を算出するために用いた R のコードは以下の通りです。

library(stringr)
file_path <- "C:\\work\\"

# 特異値分解
A <- read.table(str_c(file_path,"000.csv"))
svd_A <- svd(A)

# 近似行列の計算
for (k in c(1,5,10,15,20,25,30,35,40,45,50,75,100,750)){
  
  d <- diag(svd_A$d)[1:k, 1:k, drop=F]
  u <- svd_A$u[,1:k, drop=F]
  v <- svd_A$v[,1:k, drop=F]

  write.table(round(u %*% d %*% t(v)),
              str_c(str_c(file_path,formatC(k,width=3,flag="0")),".csv"),
              sep="\t", row.names=FALSE, col.names=FALSE)
}

出力されたcsvファイルをグレイスケール画像に変換したものが以下の図になります。

f:id:astamuse:20170613115432j:plain:w400

$k=1$ の場合

f:id:astamuse:20170613115555j:plain:w400

$k=5$ の場合

f:id:astamuse:20170613115626j:plain:w400

$k=10$ の場合

f:id:astamuse:20170613115643j:plain:w400

$k=20$ の場合

f:id:astamuse:20170613121003j:plain:w400

$k=30$ の場合

f:id:astamuse:20170613115659j:plain:w400

$k=50$ の場合

実験の結果からは、 $k$ を大きくするにつれて徐々に元の画像が再現されていく様子が確認出来ます。また $k=750$ に比べてそれほど大きくない $k$ で元の画像の情報が再現されていることも分かります。

さて、ここで $A$ が画像データではない場合（例えば統計調査によって得られたデータ等）を考えましょう。このとき $A_{k}$ が $A$ に徐々に近づいていくという現象を画像データの時と同様に確認することが出来るのでしょうか？画像データの時とは異なり、視覚的に確認するという手は使えません。

実は、近似の度合いを行列に対するノルム*7を使って定量的に評価することが出来ます。

定理３ $\| \cdot \|$ を $m\times n$ 行列のユニタリ不変ノルムとする。つまり任意の $m$ 次、 $n$ 次ユニタリ行列 $U$ 、 $V$ に対して\[\|UAV\|=\|A\|\]が成り立つとする。このとき $A_{k}$ による $A$ の近似は、ユニタリ不変ノルムに関して、階数が $k$ 以下の行列による近似の中で最良の近似になっている： \[\min_{\substack{B\in\mathrm{Mat}(m,n);\\\mathrm{rank}(B)\le k}}\|A-B\|=\|A-A_{k}\|.\] 特に作用素ノルム\[\|A\|_{op}:=\sup_{x\in\mathbf{K}^{n}\backslash\{0\}}\frac{\|Ax\|}{\|x\|}\]に関して \[\min_{\substack{B\in\mathrm{Mat}(m,n);\\\mathrm{rank}(B)\le k}}\|A-B\|_{op}=\|A-A_{k}\|_{op}=\sigma_{k+1},\] またフロベニウスノルム（ヒルベルト-シュミットノルム） \[\|A\|_{F}:=\sqrt{\sum_{i=1}^{m}\sum_{j=1}^{n}|A_{ij}|^{2}}=\sqrt{\sum_{i=1}^{\min\{m,n\}}\sigma_{i}^{2}}\]に関して \[\min_{\substack{B\in\mathrm{Mat}(m,n);\\\mathrm{rank}(B)\le k}}\|A-B\|_{F}=\|A-A_{k}\|_{F} =\sqrt{\sum_{i=k+1}^{\min\{m,n\}}\sigma_{i}^{2}}\]が成り立つ。

特異値が降順に並んでいたことを思い出すと、上記の定理により $k$ を大きくするに従って $A_{k}$ は $A$ に徐々に近づいていくことが分かります。

最後に

本記事を執筆するにあたり参考にした文献を紹介します。

新井仁之, 線形代数基礎と応用, 日本評論社, 2006

画像データを使った近似行列に関する議論は、この書籍の16.4節を参考にさせて頂きました。500ページ超とかなりボリュームのある書籍ですが、記述自体はとても丁寧であり、読み進めるのに苦労することは（多分）無いと思います。

さて、特異値分解について駆け足で見てきましたが如何だったでしょうか。記事の中でベクトルや行列について簡単に復習はしましたが、実際はこれらの概念に馴染んでいないと読み通すのは難しかったのではないでしょうか。

ベクトルや行列に馴染みの無い方は、これを機会に勉強を始めてみては如何でしょうか。馴染むまでには多少時間は掛かるかもしれませんが、機械学習や統計学の書籍を読むのに必須の知識であり、投資しただけのリターンは確実にあると思います。

*1:固有値分解についての証明は行いません。大学初年時の講義で使われるような線型代数学の教科書であれば、大抵は載っていますのでそちらを参照してください。

*2:特異値分解は内容的には線型代数学の領域に属しますが、線型代数学の書籍の多くが大学初年時に行われる授業の教科書として編纂されているということもあり、応用色が強い特異値分解について、証明付きで解説している書籍はそう多くはありません。

*3:数学の書籍では $a$ 側ではなく $b$ 側で複素共役をとる流儀が多いと思われます。ちなみに $a$ 側で複素共役をとる流儀は物理学の書籍で多く見られます。

*4:実際、線型代数学の教科書には固有方程式を解いて固有値を求める演習問題があります。

*5:いわゆるガロア理論による帰結です。

*6:行列 $A$ が対角化可能であるためには、 $A$ が正規行列、つまり $A^{*}A=AA^{*}$ を満たすことが必要十分であることが知られています。当然ですが、エルミート行列は正規行列です。

*7:行列 $A$ 、 $B$ に対して、 $\|A\|\ge 0$ ( $A\neq 0$ )、 $\|\alpha A \|=|\alpha|\|A\|$ ( $\alpha\in\mathbf{K}$ )、 $\|A+B\|\le\|A\|+\|B\|$ を満たす $\|\cdot\|$ のことをノルムと言います。