まえがき
記事の概要
今回は特異値分解を利用した 行列の低ランク近似 について解説していきます。
この記事には
- 数式によるアルゴリズムの解説
- 低ランク近似に関係する数学的な話
- プログラム実装例
が含まれています。それぞれ自分に必要な部分を参照してもらえればいいなと思っています。
前提知識
大学1年生で勉強する「線形代数学」の基本的な事項は既知として書いていきます。
キーワードは (行列の)ランク / 特異値分解 です。
特異値分解について理解してから進みたいという人は、過去に書いた 特異値分解の記事 を参考にしてもらえるといいかなと思います。
導入
行列の低ランク近似って何？
$m \times n$ 行列 $A$ と自然数 $r \ (\le m, n)$ が与えられたとき、 $A$ をランク $r$ 以下の行列 $\tilde{A}$ で近似する問題を 低ランク近似 と言います。 $A$ の最適な近似は特異値分解を利用して求めることができます。
低ランク近似ができると何が嬉しいの？
以降の節で説明するように、低ランク近似した行列はベクトルを使って効率的に表現できるようになります。一般にサイズの大きい行列の解析はとても大変なので、低ランク近似によってより小さなサイズの表現に変換するということが行われます。
ちなみに、「情報が欠落しないようにしつつ、効率的な表現に変換する」という考え方は主成分分析と共通するところがあります。
理論
基礎的な話
行列のランク
与えられた $A$ に対して、次を満たす自然数 $r$ を $A$ の ランク(階数) と言います。
- 行基本変形によって $A$ を階段行列に変換したとき、その段数が $r$ である。
これと同値な定義として、以下のようなものがあります。ここでは同値性の証明はしません。
- $A$ の一次独立な列ベクトルの最大本数が $r$ である。
- $A$ の $r$ 次小行列式が0でなく、(もし存在すれば) $r+1$ 次小行列式がすべて0である。
- $r = \text{dim} \text{Imf} A$ が成立する。
ランクの性質のうち、この記事で利用するものを挙げておきます。なお、各行列は矛盾なく和や積が定義できるようなサイズであるとしてください。
- $\text{rank}(A B) \le \min{\text{rank}(A), \text{rank}(B)}$
- $\text{rank}(C + D) \le {\text{rank}(C) + \text{rank}(D)}$
$n$ 次元列ベクトルが $n \times 1$ 行列、 $n$ 次元行ベクトルが $1 \times n$ 行列であると解釈すれば、定義から明らかに $n$ 次元ベクトルのランクは1になります。行列積のランクの性質から、 $m$ 次元ベクトル ${\bf x}$ と $n$ 次元ベクトル ${\bf y}$ に対して
$$
\text{rank}({\bf x} {\bf y}^T) = 1
$$
となることがわかります。さらに行列和のランクの性質から、ランク $r$ の行列は $r$ 個のランク1行列の和で書き表すことができることがわかります。クリックすると簡単な証明を見ることができます。
ランク $r$ の行列がランク1の行列 $r$ 個の和で表せることの証明
$$
{\bf a_i} = \sum_{j=1}^r c_{ij} {\bf e_j}
$$
として表すことができます。この列ベクトルを束ねて行列表記にすれば
$$
A = \sum_{j=1}^r {\bf e_j} (c_{1j}, \dots, c_{rj})
$$
となります。ベクトル積のランクは1(以下)なので、これで主張が示されました。
ちなみに、ランク $r$ の行列はランク1の行列 $r-1$ 個以下の和として表すことはできません。これは行列積のランクの性質から明らかだと思います。
低ランク近似問題
では、行列の低ランク近似問題を定式化します。
与えられた $m \times n$ 行列 $A$ と自然数 $r$ に対して、 次の制約付き最適化問題を 行列の低ランク近似問題 と言います。
$$
\min_{\tilde{A}} \| A - \tilde{A} \|^2_F \text{ s.t. } \text{rank}(\tilde{A}) \le r
$$
ここで $\| \cdot \|_F$ はフロべニウスノルム
$$
\| A \|_F \equiv \sqrt{\sum_{ij} \| a_{ij} \|^2}
$$
です。 $r = \text{rank}(A)$ の場合には $\tilde{A} = A$ として目的関数を0にできるため、通常 $r < \text{rank}(A)$ の場合を考えます。
ランク $r$ の行列 $\tilde{A}$ は $2r$ 個のベクトル ${\bf b_1}, \dots, {\bf b_r}, {\bf c_1}, \dots, {\bf c_r}$ によって
$$
\tilde{A} = \sum_{i=1}^r {\bf b_i} {\bf c_i}^T
$$
と表せることを前の節で確認しました。各 ${\bf b_i}, {\bf c_i}$ がそれぞれ正規化されているという仮定をおけば、係数 $d_1, \dots, d_r$ を用いて
$$
\tilde{A} = \sum_{i=1}^r d_i {\bf b_i} {\bf c_i}^T
$$
と書くこともできます。
ここで天下り的ですが、 $A$ の特異値分解について考えてみます。 $m \times n$ 行列 $A$ の特異値分解は、 $m$ 次直交行列 $U$, $m \times n$ 対角行列 $\Sigma$, $n$ 次直交行列 $V$ を用いて
$$
A = U \Sigma V^T
$$
と表されます。 $A$ のランクを $R$ とすると、 $\Sigma$ の対角成分には $r$ 個の特異値 $\sigma_1 \ge \dots \ge \sigma_r \ge 0$ が並んでいます。そこで、行列 $U, V$ を $U = ({\bf u_1}, \dots, {\bf u_m}), V = ({\bf v_1}, \dots, {\bf v_n})$ と表すことにすれば
\begin{eqnarray*}
A
&=& ({\bf u_1}, \dots, {\bf u_m}) \Sigma ({\bf v_1}, \dots, {\bf v_n})^T \\
&=& \sum_{i=1}^R \sigma_i {\bf u_i} {\bf v_i}^T
\end{eqnarray*}
となります。これは先ほど見た式 $\tilde{A} = \sum_{i=1}^r d_i {\bf b_i} {\bf c_i}^T$ によく似ていることがわかります。
実は、行列 $A$ をランク $r$ (以下)の行列 $\tilde{A}$ で近似するとき、フロべニウスノルムを最小にするような $\tilde{A}$ は
$$
\tilde{A} = \sum_{i=1}^r \sigma_i {\bf u_i} {\bf v_i}^T
$$
となることが知られています。つまり、 $A$ を特異値分解したときに得られる上位 $r$ 個の特異値の成分の和をとったものが最適な近似になっているということです。このときのフロべニウスノルムを計算してみます。
\begin{eqnarray*}
\| A - \tilde{A} \|^2_F
&=& \| U \Sigma V^T - \tilde{A} \|^2_F \\
&=& \left\| \sum_{i=1}^{\text{rank}(A)} \sigma_i {\bf u_i} {\bf v_i}^T - \sum_{i=1}^r \sigma_i {\bf u_i} {\bf v_i}^T \right\|^2_F \\
&=& \left\| \sum_{i=r+1}^{\text{rank}(A)} \sigma_i {\bf u_i} {\bf v_i}^T \right\|^2_F
\end{eqnarray*}
フロべニウスノルムの性質 $\| A \|^2_F = \text{Tr}(A^T A)$ を利用すると
\begin{eqnarray*}
\| A - \tilde{A} \|^2_F
&=& \text{Tr} \left( \left(\sum_{i=r+1}^{\text{rank}(A)} \sigma_i {\bf u_i} {\bf v_i}^T \right)^T \left( \sum_{j=r+1}^{\text{rank}(A)} \sigma_j {\bf u_j} {\bf v_j}^T \right) \right) \\
&=& \text{Tr} \left( \sum_{i=r+1}^{\text{rank}(A)} \sum_{j=r+1}^{\text{rank}(A)} \sigma_i \sigma_j \left( {\bf u_i} {\bf v_i}^T \right)^T \left( {\bf u_j} {\bf v_j}^T \right) \right) \\
&=& \text{Tr} \left( \sum_{i=r+1}^{\text{rank}(A)} \sum_{j=r+1}^{\text{rank}(A)} \sigma_i \sigma_j ({\bf v_i} {\bf u_i}^T ({\bf u_j} {\bf v_j}^T) \right) \\
&=& \text{Tr} \left( \sum_{i=r+1}^{\text{rank}(A)} \sum_{j=r+1}^{\text{rank}(A)} \delta_{ij} \sigma_i \sigma_j ({\bf v_i} {\bf v_j}^T) \right) \\
&=& \text{Tr} \left( \sum_{i=r+1}^{\text{rank}(A)} \sum_{j=r+1}^{\text{rank}(A)} \delta_{ij} \sigma_i \sigma_j ({\bf v_i}^T {\bf v_j}) \right) \\
&=& \sum_{i=r+1}^{\text{rank}(A)} \sum_{j=r+1}^{\text{rank}(A)} \delta_{ij} \sigma_i \sigma_j \\
&=& \sum_{i=r+1}^{\text{rank}(A)} \sigma_i^2
\end{eqnarray*}
となることがわかります。特異値の2乗が無視できるほど小さくなるような $r$ をとれば、ほとんど誤差を出すことなく低ランク行列で近似できることがわかります。
発展的な話
低ランク近似と特異値分解の関係(命題)
行列 $A$ の $r$ ランク近似問題を制約付き最適化問題として定式化したものを再掲します。
$$
\min_{\tilde{A}} \| A - \tilde{A} \|^2_F \text{ s.t. } \text{rank}(\tilde{A}) \le r
$$
これまで見てきたように、この最適化問題の最適解 $\tilde{A}^{\ast}$ は $A$ の特異値 $\sigma_i$ および特異ベクトル ${\bf u_i}, {\bf v_i}$ によって特徴付けられて
$$
\tilde{A}^{\ast} = \sum_{i=1}^r \sigma_i {\bf u_i} {\bf v_i}
$$
となります。 $\tilde{A}^{\ast}$ のランクは $r$ 以下になるため、許容解になる(すなわち制約条件を満たす)のは明らかです。あとは、 $\tilde{A}^{\ast}$ が他の許容解 $\tilde{A}$ よりも良い解であることを示せば十分です。
しかし、これではあまりにも天下り的です。ということで、以下の項では上の最適化問題を1から解くことで最適解を導出します。そして、その流れで自然に特異値・特異ベクトルが現れることを示したいと思います。
階数分解
証明に進む前に、 階数分解 について少しだけ説明します。ランク $r$ の $m \times n$ 行列 $A$ に対して
$$
A = B C \text{ s.t. } \text{rank}(B) = \text{rank}(C) = r
$$
を満たすような $m \times r$ 行列 $B$, $r \times n$ 行列 $C$ を見つけることを $A$ の階数分解と言います。階数分解は任意の行列に対して実行可能です。
階数分解が実行可能であることは実際にアルゴリズムを構成することで示すことができます。本題から少し逸れた話になるので、証明は折りたたんでおきます。
階数分解が実行可能であることの証明
ランク $r$ の $m \times n$ 行列 $A$ に対して階数分解を実行することを考えます。
まず $m \times r$ 行列 $B$ を構成します。
$A$ のランクは $r$ なので、 $A$ の列ベクトル ${\bf a_1}, \dots, {\bf a_n}$ から一次独立な $r$ 本のベクトルの集合 ${ {\bf a_{i_1}}, \dots, {\bf a_{i_r}} }$ をとってくることができます。これを並べて $B \equiv ({\bf a_{i_1}}, \dots, {\bf a_{i_r}})$ と定義します。構成の仕方から明らかに $B$ のランクは $r$ になります。
次に $r \times n$ 行列 $C$ を構成します。
$B$ を構成する際にとった一次独立な $r$ 本のベクトルの集合 ${ {\bf a_{i_1}}, \dots, {\bf a_{i_r}} }$ は ${ {\bf a_1}, \dots, {\bf a_n} }$ が張る空間の基底(の1つ)になっています。各 ${\bf a_i}$ はこの基底の線形結合によって表すことができます。
$$
{\bf a_i} = \sum_{j=1}^r c_{ij} {\bf a_{i_j}}
$$
この列ベクトルを束ねると
\begin{eqnarray*}
A
&=& \sum_{j=1}^r {\bf a_{i_j}} (c_{1j}, \dots, c_{nj}) \\
&=& \sum_{j=1}^r {\bf b_j} {\bf c_j}
\end{eqnarray*}
となります。ここで ${\bf c_j} \equiv (c_{1j}, \dots, c_{nj})$ としています。これを並べて $C \equiv ({\bf c_1}, \dots, {\bf c_r})^T$ とすれば
\begin{eqnarray*}
A =
\left(
\begin{array}{ccc}
{\bf b_1} & \dots & {\bf b_r}
\end{array}
\right)
\left(
\begin{array}{c}
{\bf c_1} \\
\vdots \\
{\bf c_r}
\end{array}
\right)
= B C
\end{eqnarray*}
とすることができます。今定義した $C$ に関して、ランクの性質から
$$
r = \text{rank}(A) = \text{rank}(B C) \le \min{(\text{rank}(B), \text{rank}(C))} = \min{(r, \text{rank}(C))} \
\text{rank}(C) \le \min{(r, n)} = r
$$
となり、 $\text{rank}(C) = r$ が成立します。
以上より、条件を満たすような $A$ の階数分解 $A = B C$ を得ることができました。
このアルゴリズムでは $B$ を直交行列としてとることができます。というのは、 $A$ から抽出した線形独立なベクトルの集合に対してグラム・シュミットの直交化法を適用することで、正規直交基底を得ることができるからです。(この操作は以降の操作に影響しません。)そのため、低ランク近似問題の最適解を求める際に $B$ を直交行列としてとれることを認めることにします。
低ランク近似と特異値分解の関係(証明)
行列 $A$ の $r$ ランク近似問題
$$
\min_{\tilde{A}} \| A - \tilde{A} \|^2_F \text{ s.t. } \text{rank}(\tilde{A}) \le r
$$
の最適解 $\tilde{A}^{\ast}$ が
$$
\tilde{A}^{\ast} = \sum_{i=1}^r \sigma_i {\bf u_i} {\bf v_i}
$$
で与えられることを証明します。
任意の( $r$ ランク)行列 $A$ に対して( $r$ ランク)行列 $B, C$ が存在して $A = B C$ と階数分解できること、逆に任意の $r$ ランク行列 $B, C$ の積 $A = B C$ が $r$ ランク行列になることを踏まえると、行列 $A$ の低ランク近似問題は次の最適化問題に書き換えることができます。
$$
\min_{\tilde{A}} \| A - B C \|^2_F \text{ s.t. } B^T B = I_m
$$
$B$ の直交性が制約条件になることに注意してください。1
制約付き最適化問題なので、ラグランジュの未定乗数法を用いて最適解を求めます。制約条件を成分ごとに考えると $m^2$ 個の等式制約があるので、ラグランジュ関数 ${\mathcal L}$ は
$$
{\mathcal L} \equiv \| A - B C \|^2_F + \sum_{i=1}^m \sum_{j=1}^m \lambda_{ij} ({\bf b_i}^T {\bf b_j} - \delta_{ij})
$$
となります。ここで $\lambda_{ij}$ はラグランジュ乗数で、等式制約の形から $\lambda_{ij} = \lambda_{ji}$ が成立します。$\lambda_{ij}$ を $(i, j)$ 成分とする行列 $\Lambda$ を定義すると
\begin{eqnarray*}
{\mathcal L}
&=& \| A - B C \|^2_F + \sum_{i=1}^m \sum_{j=1}^m \lambda_{ij} ({\bf b_i}^T {\bf b_j} - \delta_{ij}) \\
&=& \text{Tr}{ (A - B C)^T (A - B C) } + \sum_{i=1}^m \sum_{j=1}^m \lambda_{ij} [B^T B - I_m]_{ji} \\
&=& \text{Tr}{ (A^T A) } - 2 \text{Tr}{ (A^T B C) } + \text{Tr}{ ((B C)^T (B C)) } + \text{Tr}{ (\Lambda (B^T B - I_m)) } \\
&=& \text{Tr}{ (A^T A) } - 2 \text{Tr}{ (A^T B C) } + \text{Tr}{ (C^T B^T B C) } + \text{Tr}{ (\Lambda (B^T B - I_m)) } \\
&=& \text{Tr}{ (A^T A) } - 2 \text{Tr}{ (A^T B C) } + \text{Tr}{ (C^T C) } + \text{Tr}{ (\Lambda (B^T B - I_m)) }
\end{eqnarray*}
と変形できます。 ${\mathcal L}$ を $B$ で微分すると
\begin{eqnarray*}
\frac{\partial}{\partial B} {\mathcal L}({\bf w})
&=& - 2 \frac{\partial}{\partial B} \text{Tr}{ (B C A^T)} + \frac{\partial}{\partial B} \text{Tr}{ (B \Lambda B^T) } \\
&=& - 2 (C A^T)^T + B (\Lambda + \Lambda^T) \\
&=& - 2 (A C^T - B \Lambda)
\end{eqnarray*}
となります。 ${\mathcal L}$ を $C$ で微分すると
\begin{eqnarray*}
\frac{\partial}{\partial C} {\mathcal L}({\bf w})
&=& - 2 \frac{\partial}{\partial C} \text{Tr}{ (C A^T B)} + \frac{\partial}{\partial C} \text{Tr}{ (C^T C) } \\
&=& - 2 (A^T B)^T + 2 C \\
&=& - 2 (B^T A - C)
\end{eqnarray*}
となります。これらを ${\bf 0}$ とおくことで
\begin{eqnarray*}
A C^T &=& B \Lambda \\
B^T A &=& C
\end{eqnarray*}
を得ることができます。第2式を第1式に代入することで
$$
A A^T B = B \Lambda
$$
という式を導出できます。第1式に左から $A^T$ を作用させることで
$$
A^T A C^T = A^T B \Lambda = (B^T A)^T \Lambda = C^T \Lambda
$$
という式を導出できます。さらに第1式に左から $B^T$ を作用させると
$$
\Lambda = (B^T A) C^T = C C^T
$$
となります。以上をまとめると
\begin{eqnarray*}
A A^T B &=& B \Lambda \\
A^T A C^T &=& C^T \Lambda \\
\Lambda &=& C C^T
\end{eqnarray*}
という等式を得ることができました。ここまで頑張って導出した等式を満たす解 $B, C, \Lambda$ は停留点になっています。このときの目的関数値は
\begin{eqnarray*}
\text{Tr}{ (A - B C)^T (A - B C) }
&=& \text{Tr}{ (A^T A) } - 2 \text{Tr}{ (A^T B C) } + \text{Tr}{ ( C^T B^T B C) } \\
&=& \text{Tr}{ (A^T A) } - 2 \text{Tr}{ (C^T C) } + \text{Tr}{ ( C^T C) } \\
&=& \text{Tr}{ (A^T A) } - \text{Tr}{ (C C^T) } \\
&=& \text{Tr}{ (A^T A) } - \text{Tr}{ \Lambda }
\end{eqnarray*}
となります。
最後に $\Lambda$ ついて考察しましょう。これまでは $\Lambda$ を単なる対称行列として扱ってきましたが、実は最適性条件を満たす解 $\Lambda$ として対角行列をとることができます。
最適性条件を満たす解を1つとり、それを $\tilde{B}, \tilde{C}, \tilde{\Lambda}$ と書くことにします。 $\tilde{\Lambda}$ は対称行列なので、$m$ 次直交行列 $Q$ によって対角化できて
$$
\Lambda^{\ast} = Q^T \tilde{\Lambda} Q
$$
という対角行列 $\Lambda^{\ast}$ にすることができます。
この $Q$ を用いて、 $B^{\ast}, C^{\ast}$ を
\begin{eqnarray*}
B^{\ast} &\equiv& B Q^T \\
C^{\ast} &\equiv& Q C
\end{eqnarray*}
と定義します。新しく導入した $B^{\ast}, C^{\ast}, \Lambda^{\ast}$ は最適性条件を満たしており、これも停留点になっていることが確認できます。
したがって、 ラグランジュ乗数(を成分とする行列)が対角行列であるような停留点 $B^{\ast}, C^{\ast}, \Lambda^{\ast}$ が存在することが言えました。最適性条件は
\begin{eqnarray*}
A A^T B^{\ast} &=& B^{\ast} \Lambda^{\ast} \\
A^T A {C^{\ast}}^T &=& {C^{\ast}}^T \Lambda^{\ast}
\end{eqnarray*}
でしたが、特異値・特異ベクトルの定義から $B^{\ast}$ は左特異(列)ベクトルを並べた行列、 $C^{\ast}$ は右特異(行)ベクトルを並べた行列になっており、 $\Lambda^{\ast}$ の対角成分には特異値が並んでいます。
このときの目的関数値が
$$
\text{Tr}{ (A^T A) } - \text{Tr}{ \Lambda^{\ast} }
$$
で与えられることから、これを最小化するには $\Lambda$ の対角成分に $r$ 個の最大特異値を並べれば良いです。よって最小値は
$$
\text{Tr}{ (A^T A) } - \text{Tr}{ \Lambda^{\ast} }
= \sum_{i=1}^n \sigma_i - \sum_{i=1}^r \sigma_i
= \sum_{i=r+1}^n \sigma_i
$$
と計算することができます。(特異ベクトルは特異値に対応する形で選ぶことにします。)2
長くなってしまいましたが、これで行列の低ランク近似と特異値分解の関係を明らかにすることができました。この重要な性質は Eckart-Youngの定理 と呼ばれています。多くの文献では天下り的に特異値・特異ベクトルが登場しますが、今回示したような証明を見てみるとその気持ちが少しわかるかもしれません。
実装
行列の低ランク近似のPythonによる実装はGitHubにあります。
特異値分解のプログラム は以前自作したものを利用しています。低ランク近似のアルゴリズムそのものは単純なので、コードの内容はすぐに理解できると思います。Pythonで実装しているので実行速度はかなり遅いです。
あとがき
証明が意外と重くて大変でしたが、記事を書いたおかげでより理解が深まった気がします。計算が多い証明は眺めていてもあまり頭に入ってこないと思うので、一度自分でやってみるといいです。
参考文献・リンク
行列の低ランク近似
その他
-
「$C$ の各行は線形独立である」という制約条件を外した状態で最適化した結果得られる最適解がその制約を満たしていることから、その解は制約を加えた状態でも最適である、という性質を利用しています。本文中では $C$ の制約を雑に扱っていますが、こうした事情から $C$ については制約がないものとして進めています。 ↩
-
今回の証明では「最適解の1つとして特異値・特異ベクトルを利用した行列 $B^{\ast}, C^{\ast}$ が存在する」ことを示しました。より厳密には「 $\Lambda$ が対角行列とは限らない場合でも、目的関数値が最適値以上になる」ことを言う必要があります。とはいえ、少し考えると $\Lambda$ が対角行列である場合のみを調べれば十分であることが分かります。行列の相似を同値関係とする同値類を考えると、ある同値類に含まれるすべての行列はトレースが等しくなります。よってトレースの最大値/最小値を考える上では代表元に関してのみ考慮すれば良いです。代表元として対角行列をとることにすれば主張が示されたことになります。 ↩