パターン認識

(2) K-平均法

「凝集型クラスタリング ( Agglomerative Clustering )」は「階層的手法 ( Hierarchical Method )」の一種であり、最終的にはクラスタの階層 ( 樹形図 ; Dendrogram ) が生成されます。これに対して、クラスタとしてふさわしい状態を関数で表し、その関数の値が最適になるようにクラスタリングする手法として「非階層手法 ( Non-hierarchical Method )」があります。今回はその代表である「K-平均法 ( K-means )」から始め、さらに「混合ガウス分布 ( Mixtures of Gaussians )」「EM アルゴリズム ( Expectaion-Miximization Algorithm )」について紹介したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 分割最適化クラスタリング ( Partitional Optimization Clustering )

前章で紹介した「凝集型クラスタリング ( Agglomerative Clustering )」は、ノルムが最小となるクラスタどうしを順番に結合することでクラスタの階層を生成する方法でした。階層は「樹形図 ( Dendrogram )」と呼ばれ、任意の高さで断面を見ることでそれに応じた数のクラスタに分割された結果を得ることができるのでした。これに対し、はじめからクラスタの数を決めておき、クラスタとして最も最適な場合を得るための関数を定義して、関数値が最も最適になるような形にクラスタを分割する方法を「分割最適化クラスタリング ( Partitional Optimization Clustering )」または「非階層的クラスタリング ( Non-hierarchical Clustering )」といいます。

仮に、N 個の特徴ベクトルを k 個のクラスタに分割する場合、各場合に対して関数値を計算し、その中から最適解を得ようとすると計算回数は「第 2 種のスターリング数 ( Stirling Numbers of The Second Kind )」と呼ばれる次の値に従います。

NSk = ( 1 / k! )Σm{1→k}( ( -1 )k-mkCmmN )

この値は N に対して指数的に増加することになるため一般には最適解を計算することはできません。その代わりに、局所最適値を求めることで近似解を得る手法が用いられることになります。


2) K-平均法 ( K-means )

分割最適化クラスタリングの中で最も代表的な手法に「K-平均法 ( K-means Clustering )」があります。まずはアルゴリズムから説明します。

  1. k 個の代表点 μi ( i = 1, 2, ... k ) をランダムに選ぶ。
  2. 各特徴ベクトルに対し、最もノルムが最小となる μi を選び、その番号 i のクラスタに割り当てる。
  3. 各クラスタに割り当てられた特徴ベクトルから代表点 μi を再計算する。
  4. クラスタへの割り当て方が変化しなくなるまで 2 から 3 までを繰り返す。

一次元データ { 0, 1, 2.5, 4.9 } を例に、K-平均法を実際に使ってみます。クラスタは 2 つに分けることにしましょう。まず、代表点をランダムに決めます。通常は、実データの中から無作為に抽出することが多いのでその例にならい、1 番目の 0 と 3 番目の 2.5 を代表点とします。代表点に近い側に分類すると、データは

{ 0, 1 } { 2.5, 4.9 }

にそれぞれクラスタリングされます。代表点は 0.5, 3.7 となり、この代表点を使って再分類しても結果は上記と同じとなるので、これで処理は完了したことになります。

この処理が必ず収束するか、すなわちクラスタへの割り当て方が変化しなくなるのかが疑問として残りますが、結論から言うと必ず収束します。これは、後で示したいと思います。


それでは早速、サンプル・プログラムを示したいと思います。

/*
  InitMeans : K-平均法において代表点を初期化する

  mean : 初期化する代表点を保持する配列へのポインタ
  count : 各クラスタに含まれるデータ数(ゼロで初期化)
  p : データの次元数
*/
void InitMeans( std::vector< std::vector< double > >* mean, std::vector< unsigned int >* count,
                std::vector< std::vector< double > >::size_type p )
{
  for ( std::vector< std::vector< double > >::size_type i = 0 ; i < mean->size() ; ++i )
    (*mean)[i].assign( p, 0 );
  count->assign( mean->size(), 0 );
}

/*
  GetMinNorm : 最小のノルムになるクラスタ番号を返す

  mean : 初期化する代表点を保持する配列
  count : 対象のデータ
*/
unsigned int GetMinNorm( const std::vector< std::vector< double > >& mean, const std::vector< double >& data )
{
  L2Distance l2Distance;

  double minNorm = l2Distance( mean[0].begin(), mean[0].end(), data.begin() ); // 最小のノルム
  unsigned int index = 0; // 最小のノルムになるクラスタ番号
  for ( unsigned int i = 1 ; i < mean.size() ; ++i ) {
    double norm = l2Distance( mean[i].begin(), mean[i].end(), data.begin() );
    if ( norm < minNorm ) {
      minNorm = norm;
      index = i;
    }
  }

  return( index );
}

/*
  K-平均法によるクラスタリングを行う

  data : 対象のデータ
  cluster : 各データに対応するクラスタ番号を返す配列へのポインタ
  num : クラスタ数
*/
void Kmeans( const std::vector< std::vector< double > >& data, std::vector< unsigned int >* cluster, unsigned int num )
{
  typedef std::vector< std::vector< double > >::size_type size_type;

  size_type n = data.size(); // データ数
  if ( n == 0 ) return;
  size_type p = data[0].size(); // データの次元数
  for ( size_type i = 1 ; i < num ; ++i )
    assert( data[i].size() == p );

  std::vector< std::vector< double > > mean[2]; // プロトタイプの値(平均値)
  std::vector< unsigned int > count; // 各プロトタイプに割り当てられたデータ数
  mean[0].resize( num );
  mean[1].resize( num );
  count.resize( num );

  // mean[0] の値をデータからランダムに抽出
  for ( unsigned int i = 0 ; i < num ; ++i )
    mean[0][i] = data[i * n / num];

  cluster->assign( n, 0 );
  bool isMatched;    // データの割り当てが固定されたか
  size_t curIdx = 0; // 平均値 mean として利用する側の index
  do {
    // もう片側の mean を初期化
    InitMeans( &mean[( curIdx + 1 ) % 2], &count, p );
    isMatched = true;
    for ( size_type i = 0 ; i < n ; ++i ) {
      unsigned int index = GetMinNorm( mean[curIdx], data[i] ); // 最も近い mean の添字
      // 添字が変化した場合の処理
      if ( (*cluster)[i] != index ) {
        isMatched = false;
        (*cluster)[i] = index;
      }
      // もう片側の mean にデータを加算する
      std::transform( mean[( curIdx + 1 ) % 2][index].begin(), mean[( curIdx + 1 ) % 2][index].end(),
                      data[i].begin(),
                      mean[( curIdx + 1 ) % 2][index].begin(),
                      std::plus< double >() );
      ++( count[index] );
    }
    // 新たな mean の算出
    curIdx = ( curIdx + 1 ) % 2;
    for ( size_type i = 0 ; i < num ; ++i )
      std::transform( mean[curIdx][i].begin(), mean[curIdx][i].end(),
                      mean[curIdx][i].begin(),
                      std::bind2nd( std::divides< double >(), count[i] ) );
  } while ( ! isMatched );
}

メイン・ルーチンは Kmeans になります。処理の流れは前に説明した通りであり難しいところはないと思います。代表値を保持するための変数 mean は二つ分用意されていて、次に使う平均値をもう片側で計算するようにしてあります。また内部で、平均値を初期化するための InitMeans と、最小のノルムになる代表値を求めるための GetMinNorm を呼び出しています。関数 GetMinNorm では前章のサンプル・プログラムで紹介した L2 ノルム計算用の関数オブジェクト L2Distance を使用しています。


下図は、三つの異なる平均を持つ正規分布に従うランダムな点を K-平均法でクラスタリングした結果を示したものです。図から、各分布の中心 ( 点の密集している箇所 ) からある範囲の円内にある点を同じグループとすればうまく分類できるように見え、実際に分割数を 3 とした場合は予想した結果が得られています。しかし、分割数が 2, 4 の場合は思ったような結果になっていません。K-平均法の弱点は、あらかじめ分割数を決めておかなければならないという点です。

図 2-1. K-平均法の処理結果
分割数 3分割数 2分割数 4
3分割 2分割 4分割

K-平均法のもう一つの特徴として、代表値の初期値によって結果が変化するという点です。つまり、得られる結果は極小解 ( ある範囲内での最小解 ) であって、真の最小解であるとは限らないということになります。サンプル・プログラムでは、特徴ベクトル列の先頭から順に均等に代表点を取得しているため同じデータ列に対しては必ず解が一致しますが、代表点を完全にランダムにした場合は下図のように得られる結果はバラバラとなります。

図 2-2. 代表点による K-平均法の処理結果の差異 ( X は代表点 )
123
結果1 結果2 結果3

特徴ベクトル xi を k 番目のクラスタに割り当てた時 rik = 1, rij = 0 ( j ≠ k ) となるような指示変数 ri = { rij } を定義します。このような手法は「1-of-K 符号化法 ( 1-of-K Coding Scheme )」と呼ばれます。K 個の代表点 μj に対し、割り当てられたクラスタとの L2-ノルムの二乗和 J( ri, μj ) は以下のような式で表されます。

J( ri, μj ) = Σi( Σj( rij|| xi - μj ||2 ) )

μj を固定したとき J( ri, μj ) を最小にするには、|| xi - μj ||2 が最小になる最も近い代表点に特徴ベクトルを割り当てればいいことになります。これは、K-平均法の 2 番目のステップにあたります。このとき、ri は最も近い代表点の番号を添え字とする成分だけが 1 で、それ以外は 0 になります。今度は ri を固定して J( ri, μj ) が最小となるような μk を求めます。μk の m 番目の成分 μkm について J の偏微分 ∂J / ∂μkm を求めると、

∂J / ∂μkm=( ∂ / ∂μkmi( Σj( rij|| xi - μj ||2 ) )
=( ∂ / ∂μkmi( Σj( rijΣn( ( xin - μjn )2 ) ) )
=Σi( -2rik( xim - μkm ) )
=-2[ Σi( rikxim ) - μkmΣi( rik ) ]

となるので、J( ri, μj ) の極小値をとる μkm は正規方程式 ∂J / ∂μkm = 0 から求められて

-2[ Σi( rikxim ) - μkmΣi( rik ) ] = 0 より

μkm = Σi( rikxim ) / Σi( rik )

という結果が得られます。すなわち

μk = Σi( rikxi ) / Σi( rik )

であり、分子は k 番目のクラスタに属する特徴ベクトルの総和、分母はその数を意味するので、μk は k 番目のクラスタに属する特徴ベクトルの平均になります。これは、K-平均法の 3 番目のステップを計算したことと同じです。

2, 3 番目のステップで J( ri, μj ) は必ず前の値以下になります。ri の状態は有限個なので、それに応じて μj も有限個の値にしかならず、J( ri, μj ) は必ずある値で変化しなくなります。従って、K-平均法は必ず収束することが保証されています。しかし、前述のとおり初期値によって分類の仕方が変化することからもわかるように、得られた結果は最小解であるとは限りません。


3) 混合ガウス分布 ( Mixtures of Gaussians )

K 通りの事象がそれぞれ πj の確率で発生すると仮定します。但し、

Σj{1→K}( πj ) = 1

を満たし、K 通りの事象の和が全事象を表しているものとします。j 番目の事象に対して確率変数 x が発生する確率は、パラメータ θj を持ったある確率密度関数 pj( x | θj ) で表されるとした時、事象 j でかつ x が発生する同時確率 p( j, x ) は条件付確率の定義式から

p( j, x ) = πjpj( x | θj )

となります。従って、周辺確率 p( x ) は

p( x ) = Σj{1→K}( πjpj( x | θj ) )

となり、これ自身が明らかに確率密度関数としての条件を満たしています。このような確率分布を「混合分布 ( Mixture Distribution )」といいます。その中で

pj( x | μj, Vj ) = N( x | μj, Vj ) ≡ [ 1 / ( 2π )D/2|Vj|1/2 ] exp( -( x - μj, Vj-1( x - μj ) ) / 2 )

としたもの、すなわち正規分布の混合分布を「混合ガウス分布 ( Mixtures of Gaussians )」といいます。

K-平均法では、特徴ベクトルを必ず一つのクラスタに割り当てていました。例えば二つのクラスタの中間あたりにあるような場合でも、一方に存在するものとしてしまいます。各クラスタの発生確率があらかじめ既知であれば、先ほどの例の場合、例えばそれぞれのクラスタに対して等確率で発生するのなら、半分ずつに割り当てをするようなことが可能になります。しかし、発生確率は未知の状態なので、何らかの形で推定しなければなりません。そこで、各特徴ベクトル xi ( i = 1, 2, ... N ) が互いに独立に発生したと仮定して、対数尤度関数

Σi{1→N}( ln p( xi ) ) = Σi{1→N}( ln Σj{1→K}( πjpj( xi | θj ) ) )

が最大になる πk を求めます。但し、ここで制約条件 Σj{1→K}( πj ) = 1 を満たす必要があるので、「ラグランジュの未定乗数法 ( Lagrange Multiplier )」用いて次の式を πk で偏微分します。

Σi{1→N}( ln p( xi ) ) = Σi{1→N}( ln Σj{1→K}( πjpj( xi | θj ) ) ) - λ[ Σj{1→K}( πj ) - 1 ]

但し λ は「ラグランジュ乗数 ( Lagrange Multiplier )」を表します。

( ∂ / ∂πk ){ Σi{1→N}( ln Σj{1→K}( πjpj( xi | θj ) ) ) - λ[ Σj{1→K}( πj ) - 1 ] }
=Σi{1→N}( pk( xi | θk ) / Σj{1→K}( πjpj( xi | θj ) ) ) - λ = 0

より

Σi{1→N}( pk( xi | θk ) / Σj{1→K}( πjpj( xi | θj ) ) ) = λ

という結果が得られますが、両辺に対して πk を掛けて k について和を取ると左辺は N となり、右辺は λ となることから λ = N なので、

Σi{1→N}( pk( xi | θk ) / Σj{1→K}( πjpj( xi | θj ) ) ) = N

となります。両辺に再び πk を掛けて整理すると、

πk = Σi{1→N}( πkpk( xi | θk ) / Σj{1→K}( πjpj( xi | θj ) ) ) / N ≡ Σi{1→N}( γik ) / N

となり、右辺に現れた変数

γik ≡ πkpk( xi | θk ) / Σj{1→K}( πjpj( xi | θj ) )

を「負担率 (Responsibility)」といいます。この値は「ベイズの定理 ( Bayes' Theorem )」から、確率変数が xi であった場合にそれが k 番目のクラスタに属している確率とみなすことができます。従って、その和 Σi{1→N}( γik ) はその確率分だけ各クラスタに値を分配して得られた数と解釈できます。これを Nk とおいて、

πk = Nk / N

とすることで、非常にシンプルでわかりやすい式が得られます。

次に、πk を固定して pk( xi | θk ) のパラメータ θk の最尤解を求めます。pk( xi | θk ) は混合ガウス分布であるとして、まずは μk の要素 μkp ( p = 1, 2, ... D ) で偏微分すると、

( ∂ / ∂μkpi{1→N}( ln Σj{1→K}( πjN( xi | μj, Vj ) ) )
=Σi{1→N}( πk( ∂ / ∂μkp )N( xi | μk, Vk ) / Σj{1→K}( πjN( xi | μj, Vj ) ) )
 
( ∂ / ∂μkp )N( xi | μk, Vk )
=( ∂ / ∂μkp )[ 1 / ( 2π )D/2|Vk|1/2 ] exp( -( xi - μk, Vk-1( xi - μk ) ) / 2 )
=[ 1 / ( 2π )D/2|Vk|1/2 ] exp( -( xi - μk, Vk-1( xi - μk ) ) / 2 )・( ∂ / ∂μkp )[ -( xi - μk, Vk-1( xi - μk ) ) / 2 ]
=N( xi | μk, Vk )・( ∂ / ∂μkp )[ -( xi - μk, Vk-1( xi - μk ) ) / 2 ]
 
( ∂ / ∂μkp )[ -( xi - μk, Vk-1( xi - μk ) ) / 2 ]
=(-1/2)( ∂ / ∂μkpl{1→D}( ( xil - μklm{1→D}( vklm( xim - μkm ) ) )
=(-1/2)[ -Σm{1→D;m≠p}( vkpm( xim - μkm ) ) - Σl{1→D;l≠p}( ( xil - μkl )vklp ) - 2vkpp( xip - μkp ) ]
=m{1→D}( vkpm( xim - μkm ) ) - Σl{1→D}( vklp( xil - μkl ) )
=-v'kp( xi - μk ) - vkpT( xi - μk )

上式において、vklmVk-1 の l 行 m 列の要素を表し、v'kpVk-1 の p 番目の行ベクトル、vkpVk-1 の p 番目の列ベクトルをそれぞれ表します。ここで V-1 が対称行列であることから v'kpT = vkp なので、最終的に

( ∂ / ∂μkp )[ -( xi - μk, Vk-1( xi - μk ) ) / 2 ] = -2vkpT( xi - μk )

となり、この結果から

( ∂ / ∂μkpi{1→N}( ln Σj{1→K}( πjN( xi | μj, Vj ) ) ) = Σi{1→N}( γikvkpT( xi - μk ) )

なので、μk による偏微分は

( ∂ / ∂μki{1→N}( ln Σj{1→K}( πjN( xi | μj, Vj ) ) ) = Σi{1→N}( γikVk-1( xi - μk ) )

となります。これが 0 になるとき、両辺に Vk を掛けて整理すると

Σi{1→N}( γik( xi - μk ) ) = 0 より

μk = Σi{1→N}( γikxi ) / Nk

という結果が得られます。Vk については、μkp と同様の操作により

( ∂ / ∂Vki{1→N}( ln Σj{1→K}( πjN( xi | μj, Vj ) ) )
=Σi{1→N}( πk( ∂ / ∂Vk )N( xi | μk, Vk ) / Σj{1→K}( πjN( xi | μj, Vj ) ) )
=Σi{1→N}( πk( ∂ / ∂Vk )exp( ( -D / 2 )ln 2π - ( 1 / 2 )ln |Vk| - ( xi - μk, Vk( xi - μk ) ) / 2 ) / Σj{1→K}( πjN( xi | μj, Vj ) ) )
=Σi{1→N}( πkN( xi | μk, Vk )( ∂ / ∂Vk )[ -( 1 / 2 )ln |Vk-1| - ( xi - μk, Vk-1( xi - μk ) ) / 2 ] / Σj{1→K}( πjN( xi | μj, Vj ) ) )

と変形することができます。ここで、単一の多変量正規分布に対する共分散行列による偏微分の公式

( ∂ / ∂Vk ) ln N( xi | μj, Vj ) = ( -N / 2 )[ Vk-1 - Vk-1( xi - μk )( xi - μk )TVk-1 / N ] ≡ ( -N / 2 )( Vk-1 - Vk-1SiVk-1 )

を用いると ( 補足 1 )、上式は

Σi{1→N}( πkN( xi | μk, Vk )[ -( N / 2 )( Vk-1 - Vk-1SiVk-1 ) ] / Σj{1→K}( πjN( xi | μj, Vj ) ) )

となるので、これが 0 になる時、左右から Vk を掛けることで

Σi{1→N}( πkN( xi | μk, Vk )[ -( N / 2 )( Vk - Si ) ] / Σj{1→K}( πjN( xi | μj, Vj ) ) )
=-( N / 2 )Σi{1→N}( γik( Vk - Si ) )
=-( N / 2 )[ VkΣi{1→N}( γik ) - Σi{1→N}( γikSi ) ]
=-( N / 2 )[ NkVk - Σi{1→N}( γikSi ) ] = 0

より

Vk = Σi{1→N}( γikSi ) / Nk

となります。

K-平均法では、平均値 μk が最も近い k に特徴ベクトルを割り当てていたのに対し、混合ガウス分布を利用した場合は負担率 γik に応じて各クラスタに分配するような形になります。前者は「ハードな割り当て」であるのに対し、後者は「ソフトな割り当て」と呼ばれます。しかし、負担率を求めるためには πk, μk, Vk の値が必要になるため、まずはこれらの初期値を決めて負担率を計算します。こうして得られた負担率を使い、今度は πk, μk, Vk を求めるという操作を行います。これはちょうど、K-平均法でクラスタの平均を求める操作に該当します。K-平均法と同様にこの操作を繰り返すとやがて各パラメータは収束していきます。


混合ガウス分布を使ったクラスタリングのサンプル・プログラムを以下に示します。まずは多変量正規分布を構築するためのクラスを以下に示します。

/**
   多変量正規分布
**/
class MultiNormalDistribution
{
  // 型の定義
  typedef std::vector< double > array_type; // 確率変数の型

  array_type mu_; // 平均ベクトル
  SquareMatrix< double > inv_; // 共分散行列の逆行列
  double detV_; // 共分散行列の行列式

public:

  // 平均ベクトル mu と共分散行列 cov を指定して構築
  MultiNormalDistribution( const array_type& mu, const SymmetricMatrix< double >& cov );

  // 確率変数 a での確率密度を返す
  double operator[]( const array_type& a ) const;

  // 確率変数のサイズを返す
  array_type::size_type size() const
  { return( mu_.size() ); }
};

/*
  MultiNormalDistribution コンストラクタ : 平均ベクトル mu と共分散行列 cov を指定して構築
*/
MultiNormalDistribution::MultiNormalDistribution( const array_type& mu, const SymmetricMatrix< double >& cov )
  : mu_( mu )
{
  assert( mu.size() == cov.size() );

  // 共分散行列の逆行列を求める
  if ( ! InverseMatrix( cov, &inv_ ) )
    throw std::runtime_error( "Failed to get inverse matrix of covariance." );

  // 共分散行列の行列式を求める
  SquareMatrix< double > r;
  Householder_DoubleShiftQR< double >( cov, &r, 0, 1E-6, 10000 );
  detV_ = r[0][0];
  for ( SquareMatrix< double >::size_type i = 1 ; i < r.size() ; ++i )
    detV_ *= r[i][i];
}

/*
  MultiNormalDistribution::operator[] : 確率変数 a での確率密度を返す
*/
double MultiNormalDistribution::operator[]( const array_type& a ) const
{
  assert( a.size() == size() );

  // x = a - mu_
  array_type x( size() );
  std::transform( a.begin(), a.end(), mu_.begin(), x.begin(), std::minus< double >() );

  // xV = ( a - mu_ )V
  array_type xV( size() );
  for ( array_type::size_type i = 0 ; i < size() ; ++i )
    xV[i] = std::inner_product( x.begin(), x.end(), inv_.column( i ), double() );

  return( std::exp( -std::inner_product( xV.begin(), xV.end(), x.begin(), double() ) / 2.0 ) /
          ( std::pow( 2.0 * MathLib::Pi< double >(), size() / 2.0 ) * std::sqrt( detV_ ) )
          );
}

多変量正規分布の確率密度の計算には共分散行列の逆行列と行列式が必要になるので、コンストラクタの中で渡された共分散行列を使ってそれらを求めています。逆行列の計算には「数値計算法 (7) 連立方程式を解く -1-」の「3) 連立方程式による逆行列の計算」で紹介した関数 InverseMatrix を使っています。行列式の計算には、共分散行列が対称行列であることを利用して「固有値問題 (1) 対称行列の固有値」の「5) ハウスホルダー変換 ( Householder Transformation )」と「6) 三重対角行列とQR法」にある「ウィルキンソンの移動法 ( Wilkinson's Shift )」を利用した固有値算出用関数 Householder_DoubleShiftQR を使い、固有値をまず求めます。固有値分解は、対称行列 V を、固有値からなる対角行列 D とその固有ベクトルからなる直交行列 Q を使い V = QDQT の形に分解することを意味するのでした。Q の行列式は ±1 であり、D の行列式は固有値の積であること、行列の積の行列式が各行列の行列式の積に等しいことを利用すると、V の行列式は固有値の積に等しいことから、固有値がわかれば行列式の値を得ることができます (*3-1)。


次に、混合ガウス分布を利用したクラスタリングのサンプル・プログラムを示します。

/// 混合ガウス分布の型の宣言
typedef std::vector< std::pair< double, MultiNormalDistribution > > MixGaussian;

/*
  CalcResponsibility : 混合ガウス分布 gauss と特徴ベクトル x から負担率 r を計算する
*/
void CalcResponsibility( const MixGaussian& gauss, const Matrix< double >& x, Matrix< double >* r )
{
  r->assign( x.rows(), gauss.size() );

  for ( Matrix< double >::size_type i = 0 ; i < x.rows() ; ++i ) {
    // 特徴ベクトルをコピー
    std::vector< double > row( x.row( i ), x.row( i ).end() );
    // Σ πj・N( row | μj, Vj ) を求める
    double pSum = 0;
    for ( MixGaussian::size_type j = 0 ; j < gauss.size() ; ++j )
      pSum += gauss[j].first * gauss[j].second[row];
    // p(zk|row) = πk・N( row | μk, Vk ) / Σ πj・N( row | μj, Vj ) を求める
    if ( pSum != 0 )
      for ( MixGaussian::size_type k = 0 ; k < gauss.size() ; ++k )
        (*r)[i][k] = gauss[k].first * gauss[k].second[row] / pSum;
  }
}

/*
  CalcMixGaussian : 負担率 r と特徴ベクトル x から混合ガウス分布 gauss を計算する
*/
void CalcMixGaussian( const Matrix< double >& r, const Matrix< double >& x, MixGaussian* gauss )
{
  // Nk の計算
  std::vector< double > nk( r.cols() );
  for ( size_t k = 0 ; k < nk.size() ; ++k )
    nk[k] = std::accumulate( r.column( k ), r.column( k ).end(), double() );

  gauss->clear();
  for ( Matrix< double >::size_type k = 0 ; k < r.cols() ; ++k ) {
    // μk の計算
    std::vector< double > mu( x.cols(), double() );
    for ( Matrix< double >::size_type j = 0 ; j < x.cols() ; ++j ) {
      for ( Matrix< double >::size_type i = 0 ; i < x.rows() ; ++i )
        mu[j] += r[i][k] * x[i][j];
      if ( nk[k] != 0 ) mu[j] /= nk[k];
    }
    // Vk の計算
    SymmetricMatrix< double > cov( x.cols() );
    for ( Matrix< double >::size_type i = 0 ; i < r.rows() ; ++i ) {
      std::vector< double > xi( x.cols() ); // xi = x - mu
      std::transform( x.row( i ), x.row( i ).end(), mu.begin(), xi.begin(), std::minus< double >() );
      for ( Matrix< double >::size_type jr = 0 ; jr < cov.size() ; ++jr )
        for ( Matrix< double >::size_type jc = jr ; jc < cov.size() ; ++jc )
          cov[jr][jc] += r[i][k] * xi[jr] * xi[jc];
    }
    if ( nk[k] != 0 )
      for ( Matrix< double >::size_type jr = 0 ; jr < cov.size() ; ++jr )
        for ( Matrix< double >::size_type jc = jr ; jc < cov.size() ; ++jc )
          cov[jr][jc] /= nk[k];
    // 混合ガウス分布を求める
    gauss->push_back( std::pair< double, MultiNormalDistribution >( ( r.rows() != 0 ) ? nk[k] / r.rows() : double(), MultiNormalDistribution( mu, cov ) ) );
  }
}

/*
  EM_Algorithm_MixGaussian : 混合ガウス分布の EM アルゴリズム

  x : 特徴ベクトル
  r : 求めた負担率を保持する行列へのポインタ
  num : グループ数
  maxCnt : 最大処理回数

  戻り値 : 反復処理回数
*/
unsigned int EM_Algorithm_MixGaussian( const Matrix< double >& x, Matrix< double >* r, unsigned int num, unsigned int maxCnt )
{
  SymmetricMatrix< double > e( x.cols() ); // 単位行列(共分散行列の初期値)
  for( size_t i = 0 ; i < e.size() ; ++i )
    e[i][i] = 1;

  // μ の値をデータからランダムに抽出して混合ガウス分布を生成
  MixGaussian gauss;
  for ( unsigned int i = 0 ; i < num ; ++i ) {
    std::vector< double > mu( x.row( i * x.rows() / num ), x.row( i * x.rows() / num ).end() );
    gauss.push_back( std::pair< double, MultiNormalDistribution >( 1.0 / num, MultiNormalDistribution( mu, e ) ) );
  }

  unsigned int cnt = 1;    // 処理回数
  Matrix< double > buff_r; // 負担率用バッファ
  r->assign( x.rows(), num );
  for ( ; cnt <= maxCnt ; ++cnt ) {
    CalcResponsibility( gauss, x, &buff_r );
    CalcMixGaussian( buff_r, x, &gauss );
    // 負担率が収束したかをチェックしながらバッファからコピー
    bool isMatched = true;
    for ( Matrix< double >::size_type i = 0 ; i < r->rows() ; ++i ) {
      for ( Matrix< double >::size_type j = 0 ; j < r->cols() ; ++j ) {
        if ( std::abs( buff_r[i][j] - (*r)[i][j] ) >= 1E-6 )
          isMatched = false;
        (*r)[i][j] = buff_r[i][j];
      }
    }
    if ( isMatched ) break;
  }

  return( cnt - 1 );
}

処理内容は前に紹介したとおりです。CalcResponsibility で混合ガウス分布と特徴ベクトルを使って負担率を計算し、CalcMixGaussian で負担率と特徴ベクトルから混合ガウス分布を計算します。CalcMixGaussian の中では、まずはじめに平均 μj を計算してから、それを使って共分散行列 Vj を求めていることに注意して下さい。なお、混合ガウス分布は、混合率 πj と(多変量)正規分布 N( x | μj, Vj ) のペアからなる可変長配列 vector で表されています。

混合ガウス分布の初期値を作成する際に、μ の値は特徴ベクトルの中からランダムに抽出する方法を採用しています。これは、K-平均法で行っていたものと同じやり方です。混合ガウス分布を使ったやり方は K-平均法と比較して必要な繰り返し回数が多く、一回の処理にかかる時間も長くなるため、初期値の選び方は非常に重要になります。適切な初期値を得るためによく K-平均法が使われます。K-平均法でクラスタリングした結果から平均や共分散行列を求め、混合率は各クラスタの全体との比率から得ることができます。
初期値の選び方は結果にも影響し、K-平均法と同様に初期値によって結果は変化します。これは、尤度の極大値が一般には複数存在し、得られた結果が大域的な最大値であるとは限らないためです。


下図は、二次元データをサンプル・プログラムで二つのクラスタに分類したときの様子をいくつかの反復処理回数に対して出力した結果です。黒から白にかけて混合率で色分けしているので、例えば灰色は二つのクラスタの中間であることを表しています。14 回目で収束し、最終的には二つのクラスタにきれいに分割されていることがわかります。

図 3-1. 混合ガウス分布によるクラスタリングの推移
2 回5 回
反復処理2回 反復処理5回
10 回14 回 (収束)
反復処理10回 反復処理14回(収束)

混合ガウス分布を使ったクラスタリング処理は「EM アルゴリズム ( Expectaion-Miximization Algorithm )」と呼ばれるアルゴリズムの一種です。次の節で、一般化した EM アルゴリズムについて紹介します。


*3-1) 直交行列の行列式を含め、行列の積の行列式に関する説明は「確率・統計 (5) 正規分布」の「補足4) 行列の積の行列式」にて紹介しています。


4) EM アルゴリズム ( Expectaion-Miximization Algorithm )

観測された全データを行列 X で表し、その i 番目の特徴ベクトルは X の行ベクトル x'iT であるとします。X とは別に、本来は観測されるべきなのにそれができなかったデータを行列 Z で表し、その i 番目の特徴ベクトルは Z の行ベクトル z'iT であるとします。最後に、確率密度関数のパラメータを θ で表すと、同時確率は P( X, Z | θ ) で表せて、その対数尤度関数は ln P( X, Z | θ ) となります。特に、( x'iT, z'iT ) が各々独立であるならば、ln P( X, Z | θ ) は ln p( xiT, ziT | θ ) の和で表せて、通常の最尤推定法によって θ を求めることができます。

しかし、今 Z は未知の値となっています。そこで、同時確率を使う代わりに周辺確率 P( X | θ ) = ∫P( X, Z | θ ) dZ を利用することを検討します。ここで、∫ dZ は重積分 ∫...∫ dz1Tdz2T...dznT を意味しますが、Z が離散値であればこれは和の形になり、離散値と連続値を混合させた場合も考えることができます。このとき、対数尤度関数は

ln P( X | θ ) = ln ∫P( X, Z | θ ) dZ

となり、これを θ の一つ θj で偏微分すると、

ln P( X | θ ) / ∂θj=[ ∂P( X | θ ) / ∂θj ] / P( X | θ )
=[ ∂∫P( X, Z | θ ) dZ / ∂θj ] / P( X | θ )
=∫[ ∂P( X, Z | θ ) / ∂θj ] / P( X | θ ) dZ
=∫{ [ ∂P( X, Z | θ ) / ∂θj ] / P( X, Z | θ ) }[ P( X, Z | θ ) / P( X | θ ) ]dZ

と変形することができます。ここで、

[ ∂P( X, Z | θ ) / ∂θj ] / P( X, Z | θ ) = ∂ln P( X, Z | θ ) / ∂θj

であり、条件付確率の式から

P( X, Z | θ ) / P( X | θ ) = P( Z | X, θ )

なので、

ln P( X | θ ) / ∂θj = ∫[ ∂ln P( X, Z | θ ) / ∂θj ]P( Z | X, θ ) dZ

となりますが、このままでは式が複雑で解くことはできません。そこで、P( Z | X, θ ) の θ を定数 θ(k) として式変形すると、

ln P( X | θ ) / ∂θj=∫[ ∂ln P( X, Z | θ ) / ∂θj ]P( Z | X, θ(k) ) dZ
=∂∫ln P( X, Z | θ )P( Z | X, θ(k) ) dZ / ∂θj
∂Qk( θ, θ(k) ) / ∂θj

となります。最後の Qk( θ, θ(k) ) は、ln P( X, Z | θ ) に対する期待値を、P( Z | X, θ ) に θ(k) を代入して計算した値を表しています。偏微分が積分の外に出せるのは、P( Z | X, θ(k) ) が θ の関数ではなくなったからであることに注意して下さい。この関数を計算する部分を、期待値を計算する意味で「E ステップ」と呼びます。また、偏微分が 0 である、すなわち Qk( θ, θ(k) ) が極大値をとる θ を求める部分を極大値を求めるという意味で「M ステップ」と呼び、両者を併せて「EM アルゴリズム ( Expectaion-Miximization Algorithm )」といいます。M ステップで求められた θ は新たな値 θ(k+1) として E ステップで使われ、これを反復処理することで尤度は単調増加しやがて収束します。


混合分布は以下の式で表せるのでした。

p( x | Θ ) = Σj{1→K}( πjpj( x | θj ) )

今、特徴ベクトル x が 1 から K までのどの分布に属しているのかが既知であるとして潜在変数 z = { zj } を導入します。この変数は、属している分布の番号に対して 1 でありその他は 0 である二値変数であるとします。このとき、事前確率 p( z ) は

p( z ) = Πj{1→K}( πjzj )

z が与えられた下での x の条件付確率は

p( x | z, Θ ) = Πj{1→K}( pj( x | θj )zj )

とそれぞれ表すことができるので、同時確率 p( x, z | Θ ) は

p( x, z | Θ ) = p( x | z, Θ )p( z ) = Πj{1→K}( πjzj pj( x | θj )zj )

となります。特徴ベクトル { xi } の尤度は

L = Πi{1→N}( Πj{1→K}( πjzij pj( xi | θj )zij ) )

であり、対数尤度は

ln L = Σi{1→N}( Σj{1→K}( zij[ ln πj + ln pj( xi | θj ) ] ) )

ですが、πj に関しては総和が 1 になるという制約があるので、混合ガウス分布のところで行ったように「ラグランジュの未定乗数法 ( Lagrange Multiplier )」を用いて次の式を πk で偏微分します。

Σi{1→N}( Σj{1→K}( zij[ ln πj + ln pj( xi | θj ) ] ) ) - λ[ Σj( πj ) - 1 ]

上式を πk で偏微分すると

Σi{1→N}( zij / πk ) - λ

なので、これがゼロの時

λ = Σi{1→N}( zij ) / πk

となり、両辺に πk を掛けて k について和を取ると

Σk{1→K}( πk )λ = Σk{1→K}( Σi{1→N}( zij ) ) = N

より λ = N となります。従って、

πk = Σi{1→N}( zij ) / N

という結果が得られます。

しかし、実際には zij は未知の数です。そこで、同時確率 p( x, z ) の代わりに周辺確率 p( x ) = Σj{1→K}( p( x, z ) ) を利用すると、EM アルゴリズムにおいて、対数尤度関数は ln P( x, z | Θ ) の期待値を事後確率 P( z | x, Θ ) に関して計算した結果を意味するので

ln L=E[ Σi{1→N}( Σj{1→K}( zij[ ln πj + ln pj( x | θj ) ] ) ) ]
=Σi{1→N}( Σj{1→K}( E[zij][ ln πj + ln pj( x | θj ) ] ) )

となります。E[zij] は 事後確率 P( z | x, Θ ) に関する期待値であることに注意して

P( zij = 1 | xi, Θ ) = πjpj( xi | θj ) / Σj'{1→K}( πj'pj'( xi | θj' ) ) = γij

であり、zij が二値変数であることから ( zij = 0 は計算に寄与しないので ) 期待値 E[zij] は P( zij = 1 | xi, Θ ) に等しく、

ln L = Σi{1→N}( Σj{1→K}( γij[ ln πj + ln pj( x | θj ) ] ) )

となります。但し負担率 γij は前に求めたパラメータ { πj }, { θj } を使って計算し ( E ステップ )、次に負担率を定数として { πj }, { θj } を最大化します ( M ステップ )。このことから、混合ガウス分布を含む混合分布を使ったクラスタリングは EM アルゴリズムの特殊な場合であることがわかります。


EM アルゴリズムの反復処理は必ず収束します。これを証明するために、対数尤度 ln P( X | θ(k) ) が単調に増大することを示します。

ln P( X | θ(k+1) ) - ln P( X | θ(k) ) を以下のように変形します。

ln P( X | θ(k+1) ) - ln P( X | θ(k) )=ln ( P( X | θ(k+1) ) / P( X | θ(k) ) )
=ln ( ∫P( X, Z | θ(k+1) ) dZ / P( X | θ(k) ) )
=ln ∫[ P( X, Z | θ(k+1) ) / P( X | θ(k) ) ] dZ
=ln ∫{ [ P( X, Z | θ(k+1) ) / P( X, Z | θ(k) ) ][ P( X, Z | θ(k) ) / P( X | θ(k) ) ] } dZ
=ln ∫{ [ P( X, Z | θ(k+1) ) / P( X, Z | θ(k) ) ]P( Z | X, θ(k) ) } dZ
=ln E(k)[ P( X, Z | θ(k+1) ) / P( X, Z | θ(k) ) ]

但し、最後の期待値は事後確率 P( Z | X, θ(k) ) に対するものになります。ここで、「イェンセンの不等式 (Jensen's Inequality)」より ln E[x] ≥ E[ln x] が成り立つことから (補足 2)

ln P( X | θ(k+1) ) - ln P( X | θ(k) )E(k)[ ln [ P( X, Z | θ(k+1) ) / P( X, Z | θ(k) ) ] ]
=E(k)[ ln P( X, Z | θ(k+1) ) - ln P( X, Z | θ(k) ) ]
=E(k)[ ln P( X, Z | θ(k+1) ) ] - E[ ln P( X, Z | θ(k) ) ]
=Qk( θ, θ(k+1) ) - Qk( θ, θ(k) )

θ(k+1) は関数 Qk が最大になるように選んでいるので上式はゼロ以上となり、ln P( X | θ(k+1) ) ≥ ln P( X | θ(k) ) より対数尤度は減少することはないことが証明されます。これはもちろん、混合分布のクラスタリングでも成り立ちます。


「EM アルゴリズム」は応用範囲の広いアルゴリズムですが、特に混合ガウス分布によるそれがよく用いられることから、EM アルゴリズムというと混合ガウス分布によるクラスタリングを指すことが多いそうです。しかし、混合ガウス分布モデルは元々 EM アルゴリズムとは別に研究され、1977 年に EM アルゴリズムが発表された後で同じ構造を持つことが示されました。混合分布モデルは他にもベルヌーイ分布やポアソン分布など、離散確率密度関数に対しても使われているようです。


補足 1) 共分散行列による正規分布の偏微分

多変量正規分布の共分散行列による偏微分は結構複雑です。計算するためには行列の微分に関するいくつかの公式が必要になります。

関数 frc(x) を r 行 c 列目の要素に持つ行列を F としたとき、dfrc(x) / dx を r 行 c 列目の要素に持つ行列を dF / dx で表します。また、行列 X = { xrc } の要素を変数に持つ関数 f(X) に対し、偏微分 ∂f / ∂xrc を r 行 c 列目の要素に持つ行列を ∂f / ∂X で表します。

F, G はそれぞれ変数 x の関数 frc(x), grc(x) を r 行 c 列目の要素とする行列であるとします。このとき、行列の積 FG の r 行 c 列の成分は

Σi( fri(x)・gic(x) )

なので、その x による微分は

i( fri(x)・gic(x) ) / dx=Σi( fri(x)・g'ic(x) + f'ri(x)・gic(x) )
=Σi( fri(x)・g'ic(x) ) + Σi( f'ri(x)・gic(x) )

となります。従って、行列の微分に対する積の公式

dFG / dx = F( dG / dx ) + ( dF / dx )G

が得られます。FF-1 = E ( 単位行列 ) に積の公式を適用すると

dFF-1 / dx = F( dF-1 / dx ) + ( dF / dx )F-1 = 0 より

F( dF-1 / dx ) = -( dF / dx )F-1

であり、左側から F-1 を掛けることで

dF-1 / dx = -F-1( dF / dx )F-1

となって、逆行列の微分に対する公式が得られます。

行列式 |A| は、余因子展開により

|A| = Σi{1→n}( ( -1 )r+i・ari・cof Ari )

と変形できるのでした。cof Ari は「余因子 ( Cofactor )」といい、A から第 r 行と第 c 列を除いた行列の行列式を表します。この式を用いると、|A| の r 行 c 列の成分 arc による行列式の偏微分 ∂|A| / ∂arc

∂|A| / ∂arc = ( -1 )r+c・cof Arc

なので、

∂|A| / ∂A = A~T

が成り立ちます。ここで、A~ は「余因子行列 ( Adjugate Matrix )」であり、r 行 c 列の成分を ( -1 )c+r・cof Acr とする行列を表します。この結果を用いると、

ln |A| / ∂A=( ∂ln |A| / ∂|A| )( ∂|A| / ∂A )
=( 1 / |A| )A~T

となります。ところが、

AA~ = |A|E より

A-1 = A~ / |A|

なので、最終的に

ln |A| / ∂A = ( A-1 )T

となります。


共分散行列 V の r 行 c 列の要素を src、その逆行列 V-1 の r 行 c 列の要素を s'rc とした時、

V-1 / ∂src = -V-1( ∂V / ∂src )V-1
c
=-V-1|0,...0,...0|V-1
|:...:...:|
r|0,...1,...0|
|:...:...:|
|0,...0,...0|
=-|0,...s'1r,...0|V-1
|:...:...:|
|0,...s'cr,...0|
|:...:...:|
|0,...s'nr,...0|
=-|s'1rs'c1,...s'1rs'cc,...s'1rs'cD|
|:...:...:|
|s'rrs'c1,...s'rrs'cc,...s'rrs'cD|
|:...:...:|
|s'nrs'c1,...s'nrs'cc,...s'nrs'cD|

となるので、( xi - μ )T( ∂V-1 / ∂src )( xi - μ ) ≡ x'iT( ∂V-1 / ∂src )x'i

x'iT( ∂V-1 / ∂src )x'i
=-( x'i1, x'i2, ... x'iD )|s'1rs'c1,...s'1rs'cc,...s'1rs'cD||x'i1|
|:...:...:||x'i2|
|s'rrs'c1,...s'rrs'cc,...s'rrs'cD||:|
|:...:...:||:|
|s'nrs'c1,...s'nrs'cc,...s'nrs'cD||x'iD|
= -|s'1rs'c1x'i1,...s'1rs'ccx'i1,...s'1rs'cDx'i1||x'i1|
|:...:...:||x'i2|
|s'rrs'c1x'ir,...s'rrs'ccx'ir,...s'rrs'cDx'ir||:|
|:...:...:||:|
|s'nrs'c1x'iD,...s'nrs'ccx'iD,...s'nrs'cDx'iD||x'iD|
=-s'1rx'i1Σk{1→p}( s'ckx'ik ) + s'2rx'i2Σk{1→p}( s'ckx'ik ) + ... + s'nrx'iDΣk{1→p}( s'ckx'ik )
=j{1→p}( Σk{1→p}( s'jrx'ijx'iks'ck ) )
=-s'rTx'ix'iTs'c

と求めることができます。つまり、

x'iTV-1x'i / ∂V = -V-1x'ix'iTV-1

という結果が得られます。

単一の多変量正規分布の対数尤度は

Σi{1→N}( ln N( xi | μ, V ) ) = -( DN / 2 )ln 2π - ( N / 2 )ln |V| - ( 1 / 2 )Σi{1→N}( x'iTV-1x'i )

なので、V で微分すると、VV-1 が対称行列であることを考慮して

∂Σi{1→N}( ln N( xi | μ, V ) ) / ∂V=-( N / 2 )( V-1 )T + ( 1 / 2 )Σi{1→N}( V-1x'iTx'iV-1 )
=-( N / 2 )V-1 + ( 1 / 2 )V-1Σi{1→N}( x'iTx'i )V-1
-( N / 2 )( V-1 - V-1SV-1 )

となります。但し、

S ≡ Σi{1→N}( x'ix'iT ) / N = Σi{1→N}( ( xi - μ )( xi - μ )T ) / N

を表します。対数尤度が最大になるのは ∂Σi{1→N}( ln N( xi | μ, V ) ) / ∂V = 0 になるときなので、

-( N / 2 )( V-1 - V-1SV-1 ) = 0

より左右から V を掛けることで

V = S

となります。

補足 2) イェンセンの不等式 ( Jensen's Inequality )

f(x) が上に凸な関数であったとき、

f( E[x] ) ≥ E[ f(x) ]

が成り立ちます。これを「イェンセンの不等式 ( Jensen's Inequality )」といいます。pi を確率密度としたとき、任意の実数列 { xi } に対して

f( Σ( xipi ) ) ≥ Σ( f(xi)pi ) [ 但し Σ( pi ) = 1 ]

であり、p(x) が連続確率密度ならば

f( ∫xp(x) dx ) ≥ ∫f(x)p(x) dx

となります。

y = f(x) が上に凸ならば、m = E[x] としたとき ( m, f(m) ) を通る接線 y = f(m) + f'(m)( x - m ) に対して

f(x) ≤ f(m) + f'(m)( x - m )

となります。両辺の期待値を取れば

E[f(x)] ≤ E[f(m) + a( x - m )] = f(m)

となって、イェンセンの不等式が証明されます。


<参考文献>
  1. 「これなら分かる最適化数学」 金谷健一著 (共立出版)
  2. 「パターン認識と機械学習 下」 C.M.ビショップ著 (丸善出版)
  3. 「確率論」 伊藤清著 (岩波書店) ... 「イェンセンの不等式」の証明はここから拝借しました
  4. Wikipedia

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system