ねほり.com

何もないから何かみつかる

主成分分析を固有値問題として解いた画像処理例

   

主成分分析の話になると大抵、固有ベクトルか特異値分解のどちらかが出てきます。

行列には固有値や固有ベクトルがあるというのも、統計には分散や共分散があるというのも、理系の大卒なら誰でも知っています。

とはいえ「共分散行列の固有ベクトルが主成分の単位ベクトルになる」といわれると不思議な気がします。

ですが、あるデータのど真ん中を通る直線として、最小2乗法のように2乗誤差を最小にする直線を求める代わりに、その直線方向の分散が最大になる(その直線上に射影した場合の分散が最大になる)直線を求めると、同じ結果になるのです。

主成分分析(principal component analysis:PCA)とは?

教師なし学習の一つです。

データの分散(ばらつき)が大きいところ(主成分)を見つける操作です。

既存の特徴を組み合わせて分散が大きくなる新たな尺度となる特徴を見つけます。

例えば、2つの特徴を組み合わせて1つの特徴でその対象を上手く捉えることができたら、 パラメータを減らせます。

PCAのアルゴリズム

  1. 全データの重心を求める(平均値)
  2. 重心からデータの分散(ばらつき)が最大となる方向を見つける
  3. 新しいデータ表現軸として、2.で求めた方向を基底とする
  4. 上記でとった軸と直交する方向に対して分散が最大となる方向を探す。
  5. 2.~3.をデータの次元分だけ繰り返す

PCAの目的

  • データの特徴抽出(データのバラつきが大きい部分に着目することでよいデータを識別しやすくする)
  • データの次元圧縮(データのバラつきが少ない部分はデータに共通するパターンなのであまり意味をなさない(無視))
  • 多次元特徴量の可視化(多次元データは人間には認識不可、データのバラつきが大きところを見ることでデータの関係性を把握)

[その1]主成分分析がなぜ分散共分散行列を対角化する固有値問題となるか?

まず、最小2乗法で考えてみる。

簡単のために2次元で考える。n個のサンプルデータを

    \[X= \left( \begin{array}{ccccc} x_{11} & ..& x_1_i& ..& x_1n \\ x_{21} & ..& x_2i& ..&x_2n\\ \end{array} \right)\]

とし、第1主成分の単位ベクトルを

    \[P= \left( \begin{array}{c} p_{1} \\ p_{2} \\ \end{array} \right)\]

とすると、Xに対応する主成分軸上の第1主成分Yは

    \[Y= P^{T}X\]

であり、そのYを元の座標系に戻したものX¯は

    \[\bar{X}=PY= PP^{T}X\]

である(高校で習った一次変換)。

このX¯が、Xを第1主成分の軸上に射影したものであり、これとXとの距離が、最小にしたい誤差ということになる。

その誤差Eを、Xを直交座標とした場合の距離の2乗とすると、

であり、p12+p22=1に注意すると、これは

[その2]主成分分析がなぜ分散共分散行列を対角化する固有値問題となるか?

データ

    \[x_{i}=(x_{i1},..,x_{id})^{\mathrm{T}} (i=1,..,N)\]

の分散が最大になる方向を求める(Tは転置行列)。

  • N 個のデータからなるデータ行列:

        \[X=(x_{1},..,x_{N})^{\mathrm{T}}\]

  • N 個のデータの平均ベクトル:

        \[\bar{x}=(\bar{x}_{1},..,\bar{x}_{d})^{\mathrm{T}}\]

  • 平均ベクトルを引き算したデータ行列:

        \[\bar{X}=(x_{1} - \bar{x}..,x_{N} - \bar{x})^{\mathrm{T}}\]

とすると、平均化したデータ行列の分散は

    \[\sum = Var({\bar{X}}) = \frac{1}{N}\bar{X}^{\mathrm{T}}\bar{X}\]

で定義される。

単位ベクトル: a=(a1,⋯,ad)T とすると、N 個のデータ xi−x¯ の単位ベクトル a への正射影は

    \[s_i = (s_{1}, .. ,s_{d})^{\mathrm{T}} = \bar{Xa}\]

となる。 この変換後のデータの分散は、

    \[Var({{s}}) = \frac{1}{N}{s}^Ts = \frac{1}{N}(\bar{Xa})^{\mathrm{T}}(\bar{Xa}) = \frac{1}{N}{a}^T\bar{X}^{\mathrm{T}}\bar{Xa} = a^TVar({{\bar{X}}})a\]

となる。

この分散が最大となる単位ベクトル a は、係数ベクトルa のノルムを1となる制約があること利用して、 ラグランジェの未定乗数法を使って求める。

    \[E(a) = a^TVar({{\bar{X}}})a - \lambda(a^T a - 1)\]

を最大にするaを見つければよい。

λはラグランジェ未定定数である。

aで微分して0としておけば、

    \[\frac{\partial E(a)}{\partial a} = 2Var (\bar{X})a - 2 \lambda a = 0\]

より、

    \[Var ({\bar{X}})a = \lambda a\]

この式は元のデータの共分散行列に関する固有値問題を解くことに等しい。

つまり、分散最大となる単位ベクトル a は固有値問題を解いて求めた固有値・固有ベクトルの中で、 最大固有値に対応する固有ベクトルを aとすればよい。

なお、最大固有値とその固有ベクトルを求める方法としては、ベキ乗法(power iteration)がある。

実験例

図1に示すclass1とclass2の画像を同一の手法で対象領域と背景領域に分割し、2値画像(黒白画像)を生成してみる。

図1 class1の画像(上)とclass2の画像(下)

多機能性をもつ手法

濃淡ヒストグラムをx軸,y軸,z軸…とおき,多次元空間中で分割面を決定し2値画像を生成する。

3次元空間中にプロットした例を図2,3に示す.

ここで,黄緑色が物体領域,青色が背景領域である(左のクラスが物体領域,右のクラスが背景領域).

図2 平均,コントラスト,分散を軸とした3次元グラフ(class1)

図3 分散,エネルギー,エントロピーを軸とした3次元グラフ(class2)

3次元空間の場合,点群の共分散行列を次のように定義する.

(1)   \begin{eqnarray*} A=\sum \left( \begin{array}{ccc} x'^2_i & x'_iy'_i& x'_iz'_i \\ x'_iy'_i& y'^2_i& y'_iz'_i\\ x'_iz'_i& y'_iz'_i& z'^2_i\\ \end{array} \right) \end{eqnarray*}

ただし,

x'_i=x_i-\bar x
y'_i=y_i-\bar x
z'_i=z_i-\bar x

である.

この最大固有ベクトルを,\be_1 とし,この軸への投影を,

(2)   \begin{eqnarray*} u_i = (x_i~~ y_i~~ z_i) \be_1 \end{eqnarray*}

とする.
u_i にたいして 大津の閾値 の手法で閾値(\theta)を決める.

すなわち,次の判別式で2値化する事になる.

(3)   \begin{eqnarray*} (e_{11}~~e_{12}~~e_{13}) \left( \begin{array}{ccc} x_i\\ y_i\\ z_i\\ \end{array} \right)> \theta  \end{eqnarray*}

アルゴリズム設計

設計したアルゴリズムを図4に示す.

対象と背景を分割するアルゴリズム

まず,あらかじめ平均(mean),コントラスト(contrast),分散(variance) ,エネルギー(energy),エントロピー(entropy)を計算しておく.

次に5×5行列に拡張した式(1)を解き,最大固有ベクトルを求め閾値を求め,2値化処理を行う.

実験結果

class1の結果のヒストグラムを図5に,class2の結果のヒストグラムを図6に示す.

図5 class1のヒストグラム

図6 class2のヒストグラム

本手法で得られたclass1とclass2の2値化結果を図7,図8に示す.

図7 class1の画像結果

図8 class2の画像結果

精度評価

得られた2つの画像の精度の評価を行う.今回は正解画像が提供されているため,正解画像と出力画像との誤り画素数を計測する.

メディアンフィルタを用いる前の精度結果を表1に示す.

表1 本手法で得られた画像の精度

画像 正解画像との画素誤り(画素) 正解画像との画素誤り
class1 6046 2.306
class2 6125 2.337

表1より本手法を適用した結果,class1,class2共に約2.3\%の精度が得られた.

次にこの結果にメディアンフィルタを用いた結果を図9,図10に示す.

また,精度結果は表2に示す.

図9 class1の画像結果

図10 class2の画像結果

表2 メディアンフィルタを用いた結果

画像正解画像との画素誤り(画素)正解画像との画素誤り
class1 4046 1.543
class2 5688 2.170

表2より,メディアンフィルタが最後の処理として有効であることが分かる.

しかしながら,class1画像に平均処理を加えた画像に対し,大津の手法とメディアンフィルタを加えた結果(表3のclass1)とオイルペインティングを用いて最大量の濃度数(手動)で2値化した結果(表3のclass2)比較すると,精度は落ちている.

表3 各画像に対して最適な手法を用いた結果

画像正解画像との画素誤り(画素)正解画像との画素誤り
class1 2671 1.019
class2 2588 0.987

ソース

以下に作成したソース(必要であろうと思われる箇所のみ)。

 - 2019年(社会人15年), 学業, 調査結果

  関連記事

大学卒業して初めて大学研究の論文が論文誌に採録された

2006年04月06日(木) 研究者には戻れない やっと私の大学の研究が論文誌に …

母校に戻りスーツを着て教育実習(数学)(1/4)

教育実習のオリエンテーションが開始された。    教育実習で …

名字(苗字)の由来から先祖調査

庶民が名字をいつから名乗ったか? それは明治8年から。とよく書かれています。 そ …

人生で食べたラーメンのランキング

January 19, 2005 広島ラーメンランキング 最近またラーメンやらお …

エンジニアによる品質保証(品質説明力編)

「バグ出しのテスト」をすれば品質が保証できると思いがちですが、そうではありません …

情報処理学会「CVIM研究会」で「卒論セッション」発表

さて、2月頃には提出が決まっていた情報処理学会「コンピュータビジョンとイメージメ …

「探偵ファイルのオフ会 in 広島 」に参加してみた

調子にのって、また やっちゃいました・・・ 探偵ファイルオフ in 広島 に参加 …

「かごめ歌」を南光坊天海が作ったと言われているが・・・

やった~!10 年ぶりの「大吉」です! 賭け事も「吉」です。では早速祖母の家で賭 …

テレビを捨てた放送局、テレビを捨てた家電メーカー

「テレビの視聴者離れが進み、ネットコンテンツが充実し、テレビはオワコン」 と、よ …

毎年「当たり年」なボージョレ・ヌーボーをランキング

2011年11月02日(水) ボージョレ・ヌーボー・ランキング ボージョレ・ヌー …