(57)【特許請求の範囲】
【請求項１】
コンピュータの処理により、該コンピュータの記憶手段に記憶された複数のデータをクラタリングする方法であって、
前記コンピュータが、
前記複数のデータの間の類似度を与える分布に基づき複数のカーネル要素を計算するステップであって、該各カーネル要素には、非負の混合重みが割り当てられるステップと、
前記混合重みの添え字からなる活性要素の集合を用意するステップと、
以下のステップ(a) - (g)を前記活性要素の集合に適用するステップと、
(a) 前記複数のカーネル要素のうち１つのカーネルiを選ぶステップ、
(b) 重みが正で、カーネルiのあらわす分布に近い分布をもつ別のカーネルi'を選ぶステップ、
(c) カーネルiとカーネルi'の重みの和を計算するステップ、
(d) 前記混合重みの尤度関数について、負の該尤度関数の１次導関数を評価するステップ、
(e) もしその前記１次導関数が、前記カーネルiが重みゼロをもつ点で正であるなら、カーネルiとカーネルi'の重みの和を用いてカーネルi'の重みを更新し、前記カーネルiの重みをゼロにするとともに、前記活性要素の集合から要素iを刈り込むステップ、
(f) 一方、その１次導関数が、前記カーネルi'が重みゼロをもつ点で負であるなら、カーネルiとカーネルi'の重みの和を用いてカーネルiの重みを更新し、前記カーネルi'の重みをゼロにするとともに、前記活性要素の集合から要素i'を刈り込むステップ、及び
(g) もし前記尤度関数が単調でないなら、前記カーネルiにつき、混合重みの一方向最適化を実行するステップ。
前記混合重みの収束を判定し、まだ収束していないなら前記ステップ(a) - (g)を前記要素の集合に適用するステップに戻り、収束しているなら前記混合重みに基づき、前記複数のデータをクラタリングするステップを実行する、
クラタリング方法。
【請求項２】
前記クラタリングするステップが、前記収束した前記混合重みの非ゼロ成分に前記カーネル要素を掛けた値が最大値をとる添え字に基づき、所属するクラスタを決定する、請求項１に記載の方法。
【請求項３】
前記複数のデータの間の類似度を与える分布がガウス分布である、請求項１に記載の方法。
【請求項４】
前記複数のデータの間の類似度を与える分布がディリクレ複合多項分布である、請求項１に記載の方法。
【請求項５】
前記一方向最適化が、ニュートン・ラフソン法である、請求項１に記載の方法。
【請求項６】
コンピュータの処理により、該コンピュータの記憶手段に記憶された複数のデータをクラタリングするプログラムであって、
前記コンピュータに、
前記複数のデータの間の類似度を与える分布に基づき複数のカーネル要素を計算するステップであって、該各カーネル要素には、非負の混合重みが割り当てられるステップと、
前記混合重みの添え字からなる活性要素の集合を用意するステップと、
以下のステップ(a) - (g)を前記活性要素の集合に適用するステップと、
(a) 前記複数のカーネル要素のうち１つのカーネルiを選ぶステップ、
(b) 重みが正で、カーネルiのあらわす分布に近い分布をもつ別のカーネルi'を選ぶステップ、
(c) カーネルiとカーネルi'の重みの和を計算するステップ、
(d) 前記混合重みの尤度関数について、負の該尤度関数の１次導関数を評価するステップ、
(e) もしその前記１次導関数が、前記カーネルiが重みゼロをもつ点で正であるなら、カーネルiとカーネルi'の重みの和を用いてカーネルi'の重みを更新し、前記カーネルiの重みをゼロにするとともに、前記活性要素の集合から要素iを刈り込むステップ、
(f) 一方、その１次導関数が、前記カーネルi'が重みゼロをもつ点で負であるなら、カーネルiとカーネルi'の重みの和を用いてカーネルiの重みを更新し、前記カーネルi'の重みをゼロにするとともに、前記活性要素の集合から要素i'を刈り込むステップ、及び
(g) もし前記尤度関数が単調でないなら、前記カーネルiにつき、混合重みの一方向最適化を実行するステップ。
前記混合重みの収束を判定し、まだ収束していないなら前記ステップ(a) - (g)を前記要素の集合に適用する手段を実行し、収束しているなら前記混合重みに基づき、前記複数のデータをクラタリングするステップを実行させる、
クラタリング・プログラム。
【請求項７】
前記クラタリングするステップが、前記収束した前記混合重みの非ゼロ成分に前記カーネル要素を掛けた値が最大値をとる添え字に基づき、所属するクラスタを決定する、請求項６に記載のプログラム。
【請求項８】
前記複数のデータの間の類似度を与える分布がガウス分布である、請求項６に記載のプログラム。
【請求項９】
前記複数のデータの間の類似度を与える分布がディリクレ複合多項分布である、請求項６に記載のプログラム。
【請求項１０】
前記一方向最適化が、ニュートン・ラフソン法である、請求項６に記載のプログラム。
【請求項１１】
コンピュータの処理により、該コンピュータの記憶手段に記憶された複数のデータをクラタリングするシステムであって、
前記複数のデータの間の類似度を与える分布に基づき複数のカーネル要素を計算する手段であって、該各カーネル要素には、非負の混合重みが割り当てられる手段と、
前記混合重みの添え字からなる活性要素の集合を用意する手段と、
以下のステップ(a) - (g)を前記活性要素の集合に適用する手段と、
(a) 前記複数のカーネル要素のうち１つのカーネルiを選ぶステップ、
(b) 重みが正で、カーネルiのあらわす分布に近い分布をもつ別のカーネルi'を選ぶステップ、
(c) カーネルiとカーネルi'の重みの和を計算するステップ、
(d) 前記混合重みの尤度関数について、負の該尤度関数の１次導関数を評価するステップ、
(e) もしその前記１次導関数が、前記カーネルiが重みゼロをもつ点で正であるなら、カーネルiとカーネルi'の重みの和を用いてカーネルi'の重みを更新し、前記カーネルiの重みをゼロにするとともに、前記活性要素の集合から要素iを刈り込むステップ、
(f) 一方、その１次導関数が、前記カーネルi'が重みゼロをもつ点で負であるなら、カーネルiとカーネルi'の重みの和を用いてカーネルiの重みを更新し、前記カーネルi'の重みをゼロにするとともに、前記活性要素の集合から要素i'を刈り込むステップ、及び
(g) もし前記尤度関数が単調でないなら、前記カーネルiにつき、混合重みの一方向最適化を実行するステップ。
前記混合重みの収束を判定し、まだ収束していないなら前記ステップ(a) - (g)を前記要素の集合に適用するステップに戻り、収束しているなら前記混合重みに基づき、前記複数のデータをクラタリングする手段を有する、
クラタリング・システム。
【請求項１２】
前記クラタリングする手段が、前記収束した前記混合重みの非ゼロ成分に前記カーネル要素を掛けた値が最大値をとる添え字に基づき、所属するクラスタを決定する、請求項１１に記載のシステム。
【請求項１３】
前記複数のデータの間の類似度を与える分布がガウス分布である、請求項１１に記載のシステム。
【請求項１４】
前記複数のデータの間の類似度を与える分布がディリクレ複合多項分布である、請求項１１に記載のシステム。
【請求項１５】
前記一方向最適化が、ニュートン・ラフソン法である、請求項１１に記載のシステム。
【請求項１６】
コンピュータの処理により、該コンピュータの記憶手段に記憶された複数のデータをクラタリングする方法であって、
前記コンピュータが、
前記複数のデータの間の類似度を与える分布に基づき複数のカーネル要素を計算するステップであって、該各カーネル要素には、非負の混合重みが割り当てられるステップと、
前記混合重みの添え字からなる活性要素の集合を用意するステップと、
所与のカーネル要素と、該所与のカーネル要素のあらわす分布に近い分布をもつカーネル要素の添え字を前記活性要素の集合から選ぶステップと、
前記混合重みの尤度関数の単調性の判定に基づき、一方のカーネル要素に対応する活性配列要素からの刈り込み及び対応する混合重みを0とおくこと、他方のカーネル要素に対応する活性配列要素の刈り込み及び対応する混合重みを0とおくこと、または、一方のカーネル要素に対する一方向最適化を実行するステップと、
混合重みの収束を判定するステップと、
収束判定に応答して、混合重みに基づき、入力データ群のデータをクラスタリングするステップを実行する、
クラスタリング方法。
【請求項１７】
前記クラタリングするステップが、前記収束した前記混合重みの非ゼロ成分に前記カーネル要素を掛けた値が最大値をとる添え字に基づき、所属するクラスタを決定する、請求項１６に記載の方法。
【請求項１８】
前記複数のデータの間の類似度を与える分布がガウス分布である、請求項１６に記載の方法。
【請求項１９】
前記複数のデータの間の類似度を与える分布がディリクレ複合多項分布である、請求項１６に記載の方法。
【請求項２０】
前記一方向最適化が、ニュートン・ラフソン法である、請求項１６に記載の方法。
【請求項２１】
コンピュータの処理により、該コンピュータの記憶手段に記憶された複数のデータをクラタリングするプログラムであって、
前記コンピュータに、
前記複数のデータの間の類似度を与える分布に基づき複数のカーネル要素を計算するステップであって、該各カーネル要素には、非負の混合重みが割り当てられるステップと、
前記混合重みの添え字からなる活性要素の集合を用意するステップと、
所与のカーネル要素と、該所与のカーネル要素のあらわす分布に近い分布をもつカーネル要素の添え字を前記活性要素の集合から選ぶステップと、
前記混合重みの尤度関数の単調性の判定に基づき、一方のカーネル要素に対応する活性配列要素からの刈り込み及び対応する混合重みを0とおくこと、他方のカーネル要素に対応する活性配列要素の刈り込み及び対応する混合重みを0とおくこと、または、一方のカーネル要素に対する一方向最適化を実行するステップと、
混合重みの収束を判定するステップと、
収束判定に応答して、混合重みに基づき、入力データ群のデータをクラスタリングするステップを実行させる、
クラスタリング・プログラム。
【請求項２２】
前記クラタリングするステップが、前記収束した前記混合重みの非ゼロ成分に前記カーネル要素を掛けた値が最大値をとる添え字に基づき、所属するクラスタを決定する、請求項２１に記載のプログラム。
【請求項２３】
前記複数のデータの間の類似度を与える分布がガウス分布である、請求項２１に記載のプログラム。
【請求項２４】
前記複数のデータの間の類似度を与える分布がディリクレ複合多項分布である、請求項２１に記載のプログラム。
【請求項２５】
前記一方向最適化が、ニュートン・ラフソン法である、請求項２１に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、特徴量をもつ多数のデータの集合をクラスタリングするための技法に関するものである。
【背景技術】
【０００２】
従来より、統計解析、多変量解析、データ・マイニングなどの分野で利用される重要な技術として、クラスタリングがある。ある定義によれば、クラスタリングとは、分類対象の集合を、内的結合と外的分離が達成されるような部分集合に分割することである。
【０００３】
k-meansなどの従来の典型的なクラスタリング技法は、計算量的に簡易であるが、局所最適性(local optimality)に陥りやすいという欠点がある。また、結果の区分けは、ランダムな初期化に強く依存し、再現性に乏しい。
【０００４】
D Lashkari, P Golland, "Convex clustering with exemplar-based models", Advances in Neural Information Processing Systems 20, J. Patt, D. Koller, Y. Singer and S. Roweis, Eds, Cambridge, MA: MIT Press, 2008, pp.825-832は、ガウス混合モデルにおいて、制約されたカーネル分布で疎な混合重みを最適化する、凸クラスタリング(convex clustering)技法を開示する。この文献の開示における凸クラスタリング技法によれば、クラスタの大域最適性(global optimality)が保証されるが、そこで使用されているＥＭアルゴリズムが、極端に多数の反復計算を要し、計算時間的に妥当でないという問題がある。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】D Lashkari, P Golland, "Convex clustering with exemplar-based models", Advances in Neural Information Processing Systems 20, J. Patt, D. Koller, Y. Singer and S. Roweis, Eds, Cambridge, MA: MIT Press, 2008, pp.825-832
【発明の概要】
【発明が解決しようとする課題】
【０００６】
従って、本発明の目的は、クラスタの大域最適性を保証しつつ、計算量的に効率的なクラスタリング技法を提供することにある。
【０００７】
本発明の他の目的は、計算量的に効率的な密度推定技法を提供することにある。ここで、密度推定技法とは、限られた数の観察データを用いて、所与の分布の確率密度関数をあてはめることである。
【０００８】
本発明の更に他の目的は、データ・マイニングなどの解析に対して信頼性の高い出力を与えるクラスタリング技法を提供することにある。
【課題を解決するための手段】
【０００９】
本発明によれば、凸クラスタリング(convex clustering)における最適化手続きを高速化するために、以下のような手法が採用される：
-  先ず、入力データ群の各データ要素間の類似度を与える分布に基づき計算された固定中心、固定バンド幅の多数のカーネル要素を用意し、各カーネル要素には、非負の混合重みが割り当てられる。そして、混合重みが数値的に最適化される。
- カーネル要素の添え字に等しい初期値をもつ、活性要素(active component)の集合が用意される。
-  従来技術のＥＭアルゴリズムを、刈り込み(pruning)を行う、活性要素毎の(element-wise)最適化の反復に置き換える。
【００１０】
その反復ステップにおいて：
- １つのカーネルiを選ぶ。より大きい重みをもつカーネルには、選択のより高い優先順位が与えられる。
- 重みが正で、カーネルiのあらわす分布に近い分布をもつ別のカーネルi'を選ぶ。すなわち、
カーネルiの固定中心と固定バンド幅とカーネルi'の固定中心と固定バンド幅とが近くなるようにカーネルi'を選ぶ。
- カーネルiとカーネルi'の重みの和を計算する。
- 対数尤度関数の単調性を調べる。もし対数尤度関数が単調であるなら、カーネルiとカーネルi'の一方が刈り込まれ、他方のカーネルには、カーネルiとカーネルi'の重みの和が割り当てられる。
対数尤度関数の単調性を調べるには、例えば、次のようにする。すなわち先ず、負の対数尤度関数の１次導関数を評価する。
そして、もしその１次導関数が、カーネルiが重みゼロをもつ点で正であるなら、負の対数尤度関数はカーネルiの重みに対して単調増加であり、活性要素から要素iを刈り込む。
一方、その１次導関数が、カーネルi'が重みゼロをもつ点で負であるなら、負の対数尤度関数はカーネルiの重みに対して単調減少であり、活性要素から要素i'を刈り込む。
もし対数尤度関数が単調でないなら、カーネルiにつき、混合重みの一方向最適化を実行する。この一方向最適化とは例えば、要素毎のニュートン・ラフソン法による更新である。
- カーネルiの混合重みの現在の値において、負の対数尤度関数の１次導関数及び２次導関数を評価する。
そして、反復ステップにより混合重みが収束すると、処理は完了する。
【発明の効果】
【００１１】
この発明によれば、凸クラスタリングによりクラスタの大域最適性を保証しつつ、処理の高速化が達成される。例えば、本願発明者の実験によれば、所望の結果を得るために、本発明によれば、必要な反復ステップの回数は、ＥＭアルゴリズムを用いる凸クラスタリングの場合の１００分の１あるいは１０００分の１程度で済んだ。
【図面の簡単な説明】
【００１２】
【図１】本発明を実施するためのハードウェア構成のブロック図である。
【図２】本発明に係る機能論理ブロック図である。
【図３】本発明のクラスタリング処理のフローチャートを示す図である。
【発明を実施するための形態】
【００１３】
以下、図面に基づき、この発明の実施例を説明する。特に断わらない限り、同一の参照番号は、図面を通して、同一の対象を指すものとする。尚、以下で説明するのは、本発明の一実施形態であり、この発明を、この実施例で説明する内容に限定する意図はないことを理解されたい。
【００１４】
図１を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータ・ハードウェアのブロック図が示されている。図１において、システム・パス１０２には、ＣＰＵ１０４と、主記憶（ＲＡＭ）１０６と、ハードディスク・ドライブ（ＨＤＤ）１０８と、キーボード１１０と、マウス１１２と、ディスプレイ１１４が接続されている。ＣＰＵ１０４は、好適には、３２ビットまたは６４ビットのアーキテクチャに基づくものであり、例えば、インテル社のＰｅｎｔｉｕｍ（商標）  ４、Ｃｏｒｅ（商標）２  Ｄｕｏ、Ｃｏｒｅ（商標）２  Ｑｕａｄ、Ｘｅｏｎ（商標）、ＡＭＤ社のＡｔｈｌｏｎ（商標）などを使用することができる。主記憶１０６は、好適には、４ＧＢ以上の容量をもつものである。ハードディスク・ドライブ１０８は、クラスタリングすべき大量のデータを格納できるように、例えば、３２０ＧＢ以上の容量をもつものであることが望ましい。
【００１５】
ハードディスク・ドライブ１０８には、個々に図示しないが、オペレーティング・システムが、予め格納されている。オペレーティング・システムは、Ｌｉｎｕｘ（商標）、マイクロソフト社のＷｉｎｄｏｗｓ  ＸＰ（商標）、Ｗｉｎｄｏｗｓ（商標）２０００、アップルコンピュータのＭａｃ  ＯＳ（商標）などの、ＣＰＵ１０４に適合する任意のものでよい。
【００１６】
ハードディスク・ドライブ１０８にはまた、Ｃ、Ｃ＋＋、Ｃ＃、Ｊａｖａ（商標）などのプログラム言語処理系も格納されている。このプログラム言語処理系は、後で説明する、クラスタリング処理用のモジュールまたはツールを作成し、維持するために使用される。
【００１７】
ハードディスク・ドライブ１０８にはさらに、プログラム言語処理系でコンパイルするためのソースコードを書くためのテキスト・エディタ、及び、Ｅｃｌｉｐｓｅ（商標）などの開発環境を含んでいてもよい。
【００１８】
ハードディスク・ドライブ１０８にはさらに、クラスタリングすべきデータと、クラスタリングのための処理モジュールが保存されている。クラスタリングすべきデータ及び処理モジュールは、図２の機能ブロック図を参照して、後で説明する。
【００１９】
キーボード１１０及びマウス１１２は、オペレーティング・システムまたは、ハードディスク・ドライブ１０８から主記憶１０６にロードされ、ディスプレイ１１４に表示されたプログラム（図示しない）を起動したり、パラメータや文字を打ち込んだりするために使用される。
【００２０】
ディスプレイ１１４は、好適には、液晶ディスプレイであり、例えば、ＸＧＡ（１０２４×７６８の解像度）、またはＵＸＧＡ（１６００×１２００の解像度）などの任意の解像度のものを使用することができる。ディスプレイ１１４は、図示しないが、クラスタリングの途中経過や最終結果等を表示するために使用される。
【００２１】
図２は、本発明に係る処理モジュールの機能ブロック図である。これらのモジュールは、Ｃ、Ｃ＋＋、Ｃ＃、Ｊａｖａ（商標）など既存のプログラム言語で書かれ、実行可能バイナリ形式でハードディスク・ドライブ１０８に格納され、マウス１１２またはキーボード１１０の操作に応答して、オペレーティング・システム（図示しない）の働きで、主記憶１０６に呼び出されて、実行される。
【００２２】
図２において、ハードディスク・ドライブ１０８に格納されたデータ２０２は、クラスタリングするためのデータを含む。
【００２３】
データのクラスタリングを実行するため、本発明の一実施例のシステムは、データ取込みモジュール２０６、予備計算モジュール２０８、対数尤度関数の単調性判別モジュール２１０、刈り込みモジュール２１２、ニュートン・ラフソン法計算モジュール２１４及びクラスタリング・モジュール２１６、及びこれらのモジュールを適宜呼び出し、全体の処理を統合するメイン・ルーチン２０４を有する。
【００２４】
データ取込みモジュール２０６は、データ２０２からデータを取り込み、各々データ要素を、多次元ベクトルの形式に変換する。その際、必要に応じて、次元削減、正規化などの処理も行う。
【００２５】
予備計算モジュール２０８は、入力された各データ要素間の類似度を与える分布に基づき計算されたカーネル要素からなるカーネル行列を用意し、各カーネル要素には、非負の混合重みを割り当てるなどの処理を行う。また、活性インデックス配列、一時変数などの用意も行う。
【００２６】
対数尤度関数の単調性判別モジュール２１０は、カーネル要素についての対数尤度関数の単調性を判別する処理を行う。
【００２７】
刈り込みモジュール２１２は、活性インデックス配列から、要素を刈り込む処理を行う。
【００２８】
ニュートン・ラフソン法計算モジュール２１４は、対数尤度関数の単調性判別モジュール２１０の特定の判別条件に応じて、混合重みの値を、収束させるように更新する。
【００２９】
クラスタリング・モジュール２１６は、収束した混合重みの値に基づき、多次元ベクトルの形式のデータ要素の集合をクラスタリングする。
【００３０】
メイン・ルーチン２０４は、データ取込みモジュール２０６、予備計算モジュール２０８、対数尤度関数の単調性判別モジュール２１０、刈り込みモジュール２１２、ニュートン・ラフソン法計算モジュール２１４及びクラスタリング・モジュール２１６を適宜呼び出して、処理が進行するように制御を行う、
【００３１】
次に、図３のフローチャートを参照して、メイン・ルーチン２０４が、データ取込みモジュール２０６、予備計算モジュール２０８、対数尤度関数の単調性判別モジュール２１０、刈り込みモジュール２１２、ニュートン・ラフソン法計算モジュール２１４及びクラスタリング・モジュール２１６を適宜呼び出して実行するクラスタリング処理につい説明する。
【００３２】
ステップ３０２では、メイン・ルーチン２０４がデータ取込みモジュール２０６を呼び出し、データ２０２からデータを取り込むことにより、n個のベクトル・データx1,x2,...,xnを構成する。ここでnは、クラスタリングされるデータ要素の数である。
【００３３】
各ベクトル・データxi ( i = 1,...,n)は、d次元のベクトルであるとする。ここで、dは、各データ要素の特徴量の数である。すなわち、
xi = (xi1,xi2,...,xid)T
【００３４】
次にメイン・ルーチン２０４は、ステップ３０４で、予備計算モジュール２０８を呼び出す。すると、予備計算モジュール２０８は、m個のカーネル・ベクトルki (i=1,...,m)と、m個のカーネル・パラメータθi (i=1,...,m)を決定する。
【００３５】
このとき、nとmの大小関係は任意でよいが、ここでは便宜上、n = mと仮定する。
【００３６】
カーネル・ベクトルkiは、次の式に従い、データ要素間の類似度として定義される。
ki ≡ (p(x1|θi),p(x2|θi),...,p(xn|θi))T
すなわち、kij ≡ p(xj|θi)である。
【００３７】
一実施例では、θiは、データ要素xiに関連づけられたガウス分布の自然パラメータであり、すなわち、i=1,...,mについて、θi = (xi,σi2) 
すると、kij ≡ p(xj|xi,σi2)
【００３８】
ここでσiは、最近傍法またはパイロット・カーネル密度推定に基づく、局所適合等方性分散(locally-adaptive isotropic variances)であり、例えば次のような式で与えられる。
【数１】
0005552023.tif 000002
ここで、ε(i,j)は、iのj最近傍、すなわち、i番目のデータ要素の、j番目に近いデータ要素のインデックスをあらわす。また、||...||2は、ユークリッド・ノルムをあらわす。さらに、最近傍というときのデータ要素間のノルムも、データ要素のベクトル・データのユークリッド・ノルムであると考えてよい。
【００３９】
σiの値を用いて、kijをより詳細に書き下すと、次のとおりである。
【数２】
0005552023.tif 000003
【００４０】
次に、混合重みベクトルλの初期値を次のとおり与える。
λ ≡ (λ1 = 1/m, ..., λm = 1/m)
【００４１】
次に、活性インデックス配列、すなわち、活性要素の集合S初期値を次のとおり与える。
S = {1,2,...,m}
【００４２】
次に、各i = 1.,,,,mについて、ε(i,k)がiのk最近傍となるように、インデックス
(ε(i,1),...,ε(i,m-1))をソートし、キャッシュする。
【００４３】
さらに、次のとおり一時変数を割り当てる。
v = (v1,v2,...,vn)T
z = (z1,z2,...,zn)T
【００４４】
ところで、以下では、反復計算が行われるが、反復計算の回数をあらわす変数としてtを割り当てる。最初はまだ一度も反復計算が行われていないので、t ← 0とセットされる。
【００４５】
このtを以って、λ(t)を、t回目の反復計算でのλの値と定義する。すなわち、λの初期値は、λ(0)である。t回目の反復計算でのλのj番目の要素は、λj(t)と表記される。
【００４６】
一方、m個のカーネル・ベクトルki (i=1,...,m)を並べた行列、すなわち、
K = (k1,k2,...,km)をカーネル行列と呼ぶ。これは一般に、n×mの行列となる。
そこで、z ← Kλ(0)とセットする。
【００４７】
ここまでが、ステップ３０４で、予備計算モジュール２０８により行われる初期化処理である。次から反復計算に入る。
【００４８】
次のステップ３０６から、ステップ３２４までは、i ∈ Sについての、λiの昇順での反復計算である。
【００４９】
メイン・ルーチン２０４は、ステップ３０６であるiを選ぶと、ステップ３０８で、
i' ← mink ε(i,k)、但しε(i,k)∈ Sによって、iに基づきi'を選ぶ。
【００５０】
こうして、ステップ３０８での処理が完了すると、(i,i')というインデックスのペアが選ばれている。
【００５１】
メイン・ルーチン２０４は、ステップ３１０で、インデックスのペア(i,i')を以って、対数尤度関数の単調性判別モジュール２１０を呼び出す。対数尤度関数の単調性判別モジュール２１０は、具体的には、次のような計算を行う。
【数３】
0005552023.tif 000004
【００５２】
メイン・ルーチン２０４は、ステップ３１２で、その結果得られる、f'i0i'という値が正かどうか判断し、もしそうなら、ステップ３１４で、刈り込みモジュール２１２を呼び出して、Sからiを刈り込む。より具体的には、次のような処理を行う。
λi(t+1) ← 0
λi'(t+1) ← λi(t) + λi'(t)
z ← v
Sからiを取り除く。
【００５３】
次にステップ３２４に行って、次のiに進む。
【００５４】
ステップ３１２に戻って、f'i0i' > 0でないなら、メイン・ルーチン２０４は、処理をステップ３１６に進め、インデックスのペア(i,i')を以って、対数尤度関数の単調性判別モジュール２１０を呼び出す。対数尤度関数の単調性判別モジュール２１０は、具体的には、次のような計算を行う。
【数４】
0005552023.tif 000005
インデックスi,jの使い方が、ステップ３１０と少し異なることに留意されたい。メイン・ルーチン２０４は、ステップ３１８で、その結果得られる、f'ii'0という値が負かどうか判断し、もしそうなら、ステップ３２０で、刈り込みモジュール２１２を呼び出して、Sからiを刈り込む。より具体的には、次のような処理を行う。
λi'(t+1) ← 0
λi(t+1) ← λi(t) + λi'(t)
z ← v
Sからi'を取り除く。
【００５５】
そして、次にステップ３２４に行って、次のiに進む。
【００５６】
ステップ３１８で、f'ii'0 < 0でないなら、ステップ３２２で、メイン・ルーチン２０４は、ニュートン・ラフソン法計算モジュール２１４を呼び出して、次のような計算を行う。
【数５】
0005552023.tif 000006
【００５７】
そして、次にステップ３２４に行って、次のiに進む。
【００５８】
こうして、iについて、ステップ３０６からステップ３２４までのループを完了すると、メイン・ルーチン２０４はステップ３２６でtを1だけ増分し、ステップ３２８で、λ(t)が収束したかどうかの判断を行う。ここでの収束判定は、ある所定の正の閾値εをもちいて、
||λ(t) - λ(t-1)|| < εであることを以って収束とみなす。ここでのノルム||...||は、ユークリッド・ノルム、マンハッタン・ノルムなど任意のものでよい。
【００５９】
もしステップ３２６で、||λ(t) - λ(t-1)|| < εでないと判定されると、処理はステップ３０６に戻って、i ∈ Sについての、λiの昇順での反復計算を最初に戻って行う。
【００６０】
一方、ステップ３２６で、||λ(t) - λ(t-1)|| < εと判定されると、メイン・ルーチン２０４は、ステップ３３０に進み、クラスタリング・モジュール２１６を呼び出す。
【００６１】
ここで、λ(t) ≡ (λ1(t),λ2(t),...,λm(t))は、凸クラスタリングの性質により、疎なベクトル、すなわち、一部のλi(t)を除き、ほとんどの成分が0になる。そこで、クラスタリング・モジュール２１６は、各ベクトル・データxj (j=1,2,...,n)に対して、xjが所属するクラスターをλi(t)kijが最大となるiとする。このとき、非ゼロであるλi(t)のiだけがクラスター番号として選ばれうる。
【００６２】
なお、上記の計算では、予備計算モジュール２０８の計算で、データ要素間の類似度をガウス分布として仮定したが、これには限定されず、例えば、ディリクレ複合多項分布（ポリヤ分布とも呼ばれる）を使用してもよい。その場合、kijは、次の式で定義される。
【数６】
0005552023.tif 000007
【００６３】
この場合、θi = (μi1,μi2,...,μid,α)
そこで、μikは次のように与えられる。
加法的スムージングの場合：
【数７】
0005552023.tif 000008
減法的スムージングの場合：
【数８】
0005552023.tif 000009
これらの式で、α,β,δは、割引係数(discounting factor)であり、||...||1は、マンハッタン・ノルムをあらわす。
【００６４】
本発明で使用されるデータiとデータjの間の類似度を与える分布は、ガウス分布などの指数分布族、あるいはディリクレ複合多項分布に限定されず、クラスタリングされるデータの性質に応じた、任意の分布を使用することができる。
【００６５】
なお、上記の計算では、対数尤度関数の単調性が判別されるが、対数をとることは単調性の判別に影響を与えないので、単に尤度関数の単調性を判別することと等価であることを理解されたい。
【００６６】
また、上記の計算では、一方向最適化に、ニュートン・ラフソン法を用いたが、これには限定されず、解を含む区間の中間点を求める操作を繰り返すことによって方程式を解く求根アルゴリズムである二分法、あるいは、ニュートンラフソン法の接線の代わりに、2点を結ぶ直線(割線)を使い、この直線がx軸と交わる点を、次の近似解とする割線法などを使用することができる。
【００６７】
さらに、本発明は、コンピュータの任意のハードウェア、ソフトウェア及びプラットフォームで実施可能であるが、マルチコアあるいはマルチプロセッサの場合、対数尤度関数の単調性を判定するためのf'i0i'を計算する際に、複数ＣＰＵに処理を分割して割り当てることにより、処理を高速化することが可能である。
【符号の説明】
【００６８】
１０４    ＣＰＵ
１０６    ＲＡＭ
１０８    ハードディスク・ドライブ
２０２    データ
２０６    データ取込みモジュール
２０８    予備計算モジュール
２１０    対数尤度関数の単調性判別モジュール
２１２    刈り込みモジュール
２１４    ニュートン・ラフソン法計算モジュール
２１６    クラスタリング・モジュール
【図１】
0005552023.tif 000010
【図２】
0005552023.tif 000011
【図３】
0005552023.tif 000012