JP6824872B2

JP6824872B2 - クラスタリング装置、クラスタリング方法及びプログラム

Info

Publication number: JP6824872B2
Application number: JP2017253902A
Authority: JP
Inventors: 一生青山; 正嗣服部
Original assignee: Nippon Telegraph and Telephone Corp
Current assignee: Nippon Telegraph and Telephone Corp
Priority date: 2017-12-28
Filing date: 2017-12-28
Publication date: 2021-02-03
Anticipated expiration: 2037-12-28
Also published as: JP2019121044A

Description

本発明は、クラスタリング装置、クラスタリング方法及びプログラムに関する。

膨大な情報を効率的かつ適切に処理する技術は、処理すべき情報量の増大に伴い、重要な技術になっている。効率的に情報を処理するための１つの方法として、各々の情報の間に類似度又は非類似度、距離等の情報の間の関係性を表す尺度を定義し、類似した情報を１つにまとめる又は与えられた情報群を類似した情報の集合に分割（クラスタリング）する方法等がある。

個々の情報が非構造化データである文書、文章、テキスト、ドキュメント等（以降、これらを総じて、単に「文書」又は「オブジェクト」と表す。）である場合は、文書（オブジェクト）をbag-of-wordsモデル（以降、「BoW」と表す。）を用いて、特徴空間であるユークリッド空間におけるベクトル表現された点として扱うことが行われている。BoWは、文書中に出現する単語（word）、語彙（vocabulary）、ターム（term）等（以降では、これらを総じて、単に「単語」と表す。）に代表される文字列に文書を分割し、これらの１つ１つの文字列を文書の特徴として扱って、１つの文書を、重みが与えられた各特徴を各要素とするベクトルで表現する。各特徴に与えられる重みは、当該文書中に特徴（単語）が出現する頻度（term frequency）と、クラスタリング対象となる文書集合の各文書に当該特徴が出現する頻度（document frequency）を入力変数とする関数値との積で表される。当該関数は、tf-idf（term frequency-inverse document frequency）であることが多く、素朴な関数は以下である。

ここで、Xはクラスタリング対象の文書の特徴ベクトル集合X={x₁, x₂,・・・,x_n}、x_i（i=1, 2,・・・,n）は特徴ベクトル、|X|は文書数（すなわち、特徴ベクトル数）であり、|X|=n、t_hは単語を表す特徴、hの最大値は文書集合の異なり単語数（すなわち、特徴数）、

は単語を表す特徴t_hを含む文書数である。なお、idfを表す関数のlogの底は特に限定されない。

また、文書の長さに関する依存性を除くために、ユークリッド距離（L2ノルム）で表された特徴ベクトル長で、当該特徴ベクトルを除算し、特徴ベクトル長を１に正規化し、特徴ベクトルを単位超球（unit hypersphere）上の点として扱うことも一般的である。このとき、文書集合における総特徴数は次元数（dimensionality：以降、「dim」と表す。）となる。文書数と、文書集合の文書に出現する異なり単語数との関係は、経験則であるヒープス則（Heaps' law）として知られており、dim=O(n^β), 0<β<1である。文書数に対して異なり単語数はサブリニア（sublinear）であるが、文書数が非常に大きくなると異なり単語数も大きくなる。一方で、或る文書に出現する単語は、その文書のトピック等によって或る程度限定される。この場合、単位超球は超高次元になり、特徴ベクトル（単位超球上の点）は非常に疎なベクトル、すなわち、０以外の要素数が非常に少ないベクトルとなる。

以上は文書を一例として説明したが、例えば、購入者（又はユーザ）と購入品（又はアイテム）との関係を表す購買履歴データ（又は関係性データ）や画像から抽出された特徴がbag-of-wordsで表される画像データ等は、上記の文書と同様に扱うことができる。

単位超球上の点をクラスタリングする方法の１つとして、spherical k-means clustering（球面k-means法：以降、「SPKM」とも表す。）がある（例えば、非特許文献１参照）。通常のk-means clusteringとSPKMとは、次の（１）及び（２）の２点が異なる。

（１）平均ベクトル（mean vector）長
通常のk-means clusteringでは、クラスタメンバ（すなわち、当該クラスタに属する特徴ベクトル）の平均をクラスタの代表ベクトル（代表ベクトルは、「平均ベクトル」とも称される。）とする。一方、SPKMでは、クラスタの代表ベクトルの長さを調整し、単位超球上の点（長さ１のベクトル）を新たな代表ベクトルとする。例えば、非特許文献１に開示されているSPKMを実施する方法は、クラスタの代表ベクトルとして、クラスタメンバである特徴ベクトルの平均ベクトルを求め、当該平均ベクトルの長さをL2ノルムで正規化し、代表ベクトル（すなわち、新たな平均ベクトル）とする。SPKMでは代表ベクトルも単位超球上のベクトルとするため、誤解を招くような場合を除き、このような新たな平均ベクトルのことを「平均ベクトル」と呼ぶ。

また、後述するように、２つのオブジェクト（前述の場合は特徴ベクトル）の間の関係性の尺度としては類似度が用いられるため、代表ベクトルとしては、平均ベクトルだけでなく、クラスタに所属する特徴ベクトルの代表となるベクトルであれば他のものでも良い。例えば、中央ベクトル（median vector）を単位超球上へ射影したベクトルやクラスタに所属する特徴ベクトルのメドイドベクトル（medoid vector）等を代表ベクトルとしても良い。

（２）２つのオブジェクト（文書）の関係性を表す尺度
通常のk-means clusteringでは関係性の尺度として、ユークリッド距離（又はユークリッド距離の二乗）が用いられるのに対して、SPKMでは、コサイン類似度（cosine similarity）、すなわち、２つの特徴ベクトルの内積が用いられる。或るオブジェクトの特徴ベクトルx_i（x_iは列ベクトル表現）と、或るクラスタの代表ベクトル（平均ベクトル）μ_j（μ_jは列ベクトル表現）とを用いると、２つの尺度の関係は、以下の式１の通りである。

ここで、||x_i||₂はベクトルx_iのL2ノルム、θはベクトルx_iとμ_jとの成す角を表す。

上記の式１の左辺が２つのベクトルのユークリッド距離の二乗である。ベクトルが単位超球上の点である場合、ベクトル長は１に正規化されているため、上記の式１は、以下の式２の通りとなる。

ここで、x_i ^Tはベクトルx_iの転置、x_i ^Tμ_jはベクトルx_iとμ_jとの内積を表す。このベクトルの内積がコサイン類似度に相当する。

上記の（１）及び（２）の相違により、通常のk-means clusteringとSPKMとでは、クラスタリングを実行する際の目的関数J(C;X,k)が異なる。ただし、Cは、オブジェクトが重複せずに（すなわち、同一のオブジェクトが複数のクラスタに属することなく）、Xをk個のクラスタに分割した際のクラスタの集合C={C₁, C₂,・・・,C_k}を表す。

このとき、通常のk-means clusteringの目的関数J_KM(C;X,k)と、SPKMの目的関数J_SPKM(C;X,k)とは、それぞれ以下の式３及び式４で表される。

ここで、μ_jはクラスタC_jの代表ベクトル（平均ベクトル）を表す。

通常のk-means clusteringの目的は、上記の式３の目的関数J_KM(C;X,k)を最小化するクラスタ集合C^* _KMを求めることであり、SPKMの目的は、上記の式４の目的関数J_SPKM(C;X,k)を最大化するクラスタ集合C^* _SPKMを求めることである。すなわち、通常のk-means clusteringの目的と、SPKMの目的とは、それぞれ以下の式５及び式６を求めることである。

ただし、特徴ベクトルと平均ベクトルとの双方が単位超球上の点となる場合、すなわち、SPKMの問題設定時には、上記の式５及び式６は等価になる。

大規模データに対して、上記の式５を厳密に解くことは膨大な計算量を要するため、発見的方法（heuristic algorithm）であるLloyd-type algorithm（以降、「Lloyd法」と表す。）が知られている。Lloyd法の概略は次の通りである。

Ｓ１）初期状態生成：与えられたオブジェクト集合の各オブジェクトをランダムにk個のクラスタに割り当てた初期状態を生成し、各クラスタの平均ベクトルを求める。

Ｓ２）割当ステップ（assignment step）：各オブジェクトと、全ての平均ベクトルとの距離を計算し、最近傍の平均ベクトルのクラスタに当該オブジェクトを割り当てる。

Ｓ３）更新ステップ（update step）：新たに割り当てられたクラスタのオブジェクトを用いて、各クラスタの平均ベクトルを計算する。

Ｓ４）終了又は繰り返し：予め設定された終了条件を満たした場合、アルゴリズムを終了し、クラスタ割当又は平均ベクトルを出力する。終了条件を満たさない場合は、上記のＳ２）及びＳ３）を繰り返す。この繰り返す操作を、iteration又は繰り返しと呼び、現在のiteration（繰り返し）、１つ前のiteration（繰り返し）と呼ぶことで、繰り返しを特定するものとする。なお、終了条件としては、全ての平均ベクトルが更新されなくなる場合（すなわち、オブジェクトの割当が不変となる場合）、目的関数値の変化が設定した閾値未満に場合等、種々の条件を用いることができる。

Lloyd法により問題を解く場合、上記の式５及び式６の等価性を考慮すると、コサイン類似度を用いた場合とユークリッド距離（L2ノルム）を用いた場合とは、初期状態が同じであれば同じ解を得ることになる。

一般的には、上記のＳ２）の割当ステップにおける距離計算のコストがアルゴリズムの計算量を支配する。このため、距離計算コストを削減し、アルゴリズムの高速化を図る方法が知られている（例えば、非特許文献２乃至４参照）。これらに共通し、最も距離計算コストを削減する効果のある方法が、平均ベクトルの繰り返し間における距離（１つ前のiteration時の平均ベクトルと、現在のiteration時の平均ベクトルとの距離）と、オブジェクトと平均ベクトルとの距離の上限値(upper bound)及び下限値(lower bound)とにより、不要な距離計算を省く方法である。なお、オブジェクトと平均ベクトルとの距離の上下限値は、距離公理の１つである三角不等式を用いて求められる。

一方、上記の式６に対しては、その発見的方法をアルゴリズム的側面から高速にする方法はあまり知られていない。

I. S. Dhillon and D. S. Modha, "Concept decomposition for large sparse text data using clustering," Machine Learning, vol. 42, pp. 143-175, 2001. C. Elkan, "Using the triangle inequality to accelerate k-means," Proc. Int. Conf. Machine Learning (ICML), 2003. G. Hamerly, "Making k-means even faster" Proc. SIAM Int. Conf. Data Mining (SDM), pp. 130-140, 2010. Y. Ding, Y. Zhao, X. Shen, M. Musuvathi, and T. Mytkowicz, "Yinyang k-means: A drop-in replacement of the classic k-means with consistent speedup," Proc. Int. Conf. Machine Learning (ICML), 2015.

このため、SPKMを大規模データに適用すると膨大な計算量を要するという問題がある。これに対して、上述した等価性を考慮して、通常のk-means clusteringの高速化に利用されている方法（上述した、オブジェクトと平均ベクトルとの距離の上限値及び下限値を利用する方法）を、平均ベクトルを単位超球上の点となるように正規化しユークリッド距離を用いた方法に適用することも考えられる。しかしながら、この方法の計算効率は必ずしも良いものではない。

そこで、本発明は、球面上のベクトルをクラスタリングする場合の計算コストを削減することを目的とする。

そこで、本発明の実施の形態では、球面上にそれぞれ配置された複数の特徴ベクトルを繰り返し手続きにより複数のクラスタに割り当てるクラスタリング装置であって、前記特徴ベクトル毎に、前記特徴ベクトルと、前記複数のクラスタそれぞれの代表ベクトルとの類似度を計算する類似度計算手段と、繰り返しの回数がｔ回目における前記類似度の計算を行うか否かを判定する判定手段と、を有し、前記判定手段は、前記特徴ベクトルと前記代表ベクトルとが成す角に基づいて前記類似度に対して算出される上限値であって、前記繰り返しの回数がｔ−１回目における上限値を少なくとも用いて、前記繰り返しの回数がｔ回目における前記類似度の計算を行うか否かを判定し、前記類似度計算手段は、前記判定手段により前記類似度の計算を行わないと判定された場合に、前記ｔ回目における前記類似度の計算を行わない。

球面上のベクトルをクラスタリングする場合の計算コストを削減することができる。

本発明の実施の形態におけるクラスタリング装置の機能構成の一例を示す図である。本発明の実施の形態におけるクラスタリング装置のハードウェア構成の一例を示す図である。本発明の実施の形態におけるクラスタリング装置が実行する全体処理の一例を示すフローチャートである。コサイン類似度計算の一例を説明するための図である。本発明の実施の形態におけるクラスタリング装置が実行する更新処理の一例を示すフローチャートである。扇形回転体が第１象限に含まれる場合における上下限値計算の一例を説明するための図である。扇形回転体が第１象限に含まれない場合における上限値計算の一例を説明するための図（その１）である。扇形回転体が第１象限に含まれない場合における上限値計算の一例を説明するための図（その２）である。扇形回転体が第１象限に含まれない場合における下限値計算の一例を説明するための図である。本発明と従来法との比較例を示す図である。コサイン類似度とL2ノルム（ユークリッド距離）との比較例を示す図である。コサイン類似度の代わりにL2ノルム（ユークリッド距離）を用いた計算の一例を説明するための図である。 L2ノルム（ユークリッド距離）を用いた厳密距離計算回数の一例を説明するための図である。

以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。以降では、SPKMを用いたクラスタリングを行うクラスタリング装置１０について説明する。

＜機能構成＞
まず、本発明の実施の形態におけるクラスタリング装置１０の機能構成について、図１を参照しながら説明する。図１は、本発明の実施の形態におけるクラスタリング装置１０の機能構成の一例を示す図である。

図１に示すように、本発明の実施の形態におけるクラスタリング装置１０は、入力部１０１と、初期設定部１０２と、割当部１０３と、更新部１０４と、終了判定部１０５と、出力部１０６とを有する。これら各機能部は、クラスタリング装置１０にインストールされた１以上のプログラムCPU（Central Processing Unit）等に実行させる処理により実現される。

また、本発明の実施の形態におけるクラスタリング装置１０は、記憶部１０７を有する。記憶部１０７は、例えばHDD（Hard Disk Drive）やSSD（Solid State Drive）等の記憶装置を用いて実現可能である。なお、記憶部１０７は、例えば、クラスタリング装置１０とネットワークを介して接続される記憶装置等を用いて実現されていても良い。

入力部１０１は、オブジェクトから抽出された特徴ベクトルの集合X={x₁, x₂,・・・,x_n}と、クラスタ数k（<n）とを入力する。ただし、||x_i||₂=1（すなわち、各特徴ベクトルが示す点は、単位超球上の点）であるものとする。なお、|X|=nである。

初期設定部１０２は、初期状態を設定する。すなわち、初期設定部１０２は、初期状態として、各特徴ベクトルを各クラスタにランダムに割り当てると共に、初期平均ベクトルとして各クラスタから特徴ベクトルを無作為に選択する。

割当部１０３は、特徴ベクトルと平均ベクトルとの上限値（以降、「ＵＢ」とも表す。）及び下限値（以降、「ＬＢ」とも表す。）を用いて、コサイン類似度（以降、単に「類似度」とも表す。）の計算（厳密な類似度計算）を行うか否かを判定する。或る特徴ベクトルのＬＢとは、当該特徴ベクトルが属するクラスタの平均ベクトルと、当該特徴ベクトルとの類似度の下限値のことである。一方で、或る特徴ベクトルのＵＢとは、各クラスタの平均ベクトルのうち２番目に類似する平均ベクトル（すなわち、２番目に類似度が高い平均ベクトル。これを「第二類似平均ベクトル」とも表す。）と、当該特徴ベクトルとの類似度の上限値とのことである。

そして、割当部１０３は、コサイン類似度の計算を行うと判定した場合に、特徴ベクトルと、各クラスタの平均ベクトルとの類似度を計算（厳密な類似度計算）し、最も類似する平均ベクトル（これを「最類似平均ベクトル」とも表す。）及び第二類似平均ベクトルを求める。このとき、割当部１０３は、最類似平均ベクトルのクラスタを識別するクラスタＩＤと、最類似平均ベクトルとの類似度と、第二類似平均ベクトルのクラスタを識別するクラスタＩＤと、第二類似平均ベクトルとの類似度とを記憶部１０７に記憶させる。

ここで、割当部１０３がコサイン類似度の計算を行うと判定される場合は、ＬＢがＵＢ未満であって、かつ、ｓｉｍがＵＢ未満である場合である。ｓｉｍとは、特徴ベクトルが属するクラスタの平均ベクトルと、当該特徴ベクトルとの類似度である。一方で、割当部１０３がコサイン類似度の計算を行うと判定されない場合は、ＬＢがＵＢ以上である場合、ＬＢがＵＢ未満であって、かつ、ｓｉｍがＵＢ以上である場合である。

また、割当部１０３は、記憶部１０７に記憶されているクラスタＩＤを用いて、特徴ベクトルが属するクラスタを変更する。すなわち、割当部１０３は、特徴ベクトルが属するクラスタの平均ベクトルと、最類似平均ベクトルとが同一でない場合、当該特徴ベクトルを、当該最類似平均ベクトルのクラスタに割り当てる。

更新部１０４は、平均ベクトルの計算及び更新と、各特徴ベクトルのＵＢ及びＬＢの計算とを行う。

以上の割当部１０３による処理と、更新部１０４による処理とは、予め設定された終了条件を満たすまで繰り返し実行される。終了条件としては、上述したように、例えば、全ての平均ベクトルが更新されなくなる場合、目的関数値の変化が設定した閾値未満に場合等、種々の条件を用いることができる。

終了判定部１０５は、終了条件を満たしたか否かを判定する。出力部１０６は、終了判定部１０５により終了条件を満たしたと判定された場合、特徴ベクトルのクラスタへの割当結果（これは、オブジェクトのクラスタへの割当結果と同一視できる。）を出力する。なお、出力部１０６は、各クラスタの平均ベクトルを出力しても良い。出力部１０６による出力先は、限定されない。出力先の具体例としては、例えば、ディスプレイ、補助記憶装置、外部記録媒体、ネットワークを介して接続されるサーバ装置等が挙げられる。

以上のように、本発明の実施の形態のクラスタリング装置１０は、当該特徴ベクトルが属するクラスタの平均ベクトルと、当該特徴ベクトルとの類似度の下限値（ＬＢ）が、第二類似平均ベクトルと、当該特徴ベクトルとの類似度の上限値（ＵＢ）未満である場合にのみ、平均ベクトルとの厳密な類似度を計算する。これにより、本発明の実施の形態のクラスタリング装置１０では、類似度の計算回数を低減（すなわり、計算コストを削減）することができ、球面k-means法（SPKM）を用いたクラスタリングの計算時間を短縮することができる。

＜ハードウェア構成＞
次に、本発明の実施の形態におけるクラスタリング装置１０のハードウェア構成について、図２を参照しながら説明する。図２は、本発明の実施の形態におけるクラスタリング装置１０のハードウェア構成の一例を示す図である。

図２に示すように、本発明の実施の形態におけるクラスタリング装置１０は、入力装置２０１と、表示装置２０２と、外部I/F２０３と、RAM（Random Access Memory）２０４と、ROM（Read Only Memory）２０５と、CPU（Central Processing Unit）２０６と、通信I/F２０７と、補助記憶装置２０８とを有する。これら各ハードウェアは、それぞれがバス２０９を介して通信可能に接続されている。

入力装置２０１は、例えばキーボードやマウス、タッチパネル等であり、ユーザが各種装置を入力するのに用いられる。表示装置２０２は、例えばディスプレイ等であり、クラスタリング装置１０の処理結果等を表示する。なお、クラスタリング装置１０は、入力装置２０１及び表示装置２０２のうちの少なくとも一方を有していなくても良い。

外部I/F２０３は、外部装置とのインタフェースである。外部装置には、記録媒体２０３ａ等がある。クラスタリング装置１０は、外部I/F２０３を介して、記録媒体２０３ａ等の読み取りや書き込みを行うことができる。記録媒体２０３ａには、例えば、本発明の実施の形態におけるクラスタリング装置１０が有する各機能部を実現するためのプログラムが格納されていても良い。

記録媒体２０３ａには、例えば、フレキシブルディスク、CD（Compact Disc）、DVD（Digital Versatile Disk）、SDメモリカード（Secure Digital memory card）、USB（Universal Serial Bus）メモリカード等がある。

RAM２０４は、プログラムやデータを一時保持する揮発性の半導体メモリである。ROM２０５は、電源を切ってもプログラムやデータを保持することができる不揮発性の半導体メモリである。ROM２０５には、例えば、OS（Operating System）設定やネットワーク設定等が格納されている。CPU２０６は、ROM２０５や補助記憶装置２０８等からプログラムやデータをRAM２０４上に読み出して処理を実行する演算装置である。

通信I/F２０７は、クラスタリング装置１０が他の装置等と通信を行うためのインタフェースである。本発明の実施の形態におけるクラスタリング装置１０が有する各機能部を実現するためのプログラムは、例えば、通信I/F２０７を介して、所定のサーバ等から取得（ダウンロード）されても良い。また、本発明の実施の形態におけるクラスタリング装置１０は、例えば、通信I/F２０７を介して、これら各機能部を実現するためのプログラムを他の装置に提供しても良い。

補助記憶装置２０８は、例えばHDDやSSD等であり、プログラムやデータを格納している不揮発性の記憶装置である。補助記憶装置２０８に格納されているプログラムやデータには、例えば、OS、当該OS上において各種機能を実現するアプリケーションプログラム、本発明の実施の形態におけるクラスタリング装置１０が有する各機能部を実現するためのプログラム等がある。

本発明の実施の形態におけるクラスタリング装置１０は、上記の各ハードウェア構成を有することにより、後述する各種処理を実現することができる。

＜全体処理＞
次に、本発明の実施の形態におけるクラスタリング装置１０が実行する全体処理について、図３を参照しながら説明する。図３は、本発明の実施の形態におけるクラスタリング装置１０が実行する全体処理の一例を示すフローチャートである。

ステップＳ１０１：入力部１０１は、特徴ベクトル集合Xと、クラスタ数kとを入力する。

ステップＳ１０２：初期設定部１０２は、初期状態を設定する。すなわち、初期設定部１０２は、初期状態として、各特徴ベクトルを各クラスタにランダムに割り当てると共に、初期平均ベクトルとして各クラスタから特徴ベクトルを無作為に選択する。

ここで、初期設定部１０２により初期状態が設定された直後は、以降のステップＳ１０３及びステップＳ１０４の処理は実行されない。言い換えれば、割当部１０３による処理と更新部１０４による処理とが繰り返し実行された回数（すなわち、iteration（繰り返し又は反復）回数）を表すインデックスをｔとして、ｔ＝０（すなわち、繰り返し回数が０回）のときは、以降のステップＳ１０３及びステップＳ１０４の処理は実行されない。

以降のステップＳ１０３〜ステップＳ１０８の処理は、特徴ベクトル集合Xに含まれる特徴ベクトル毎に実行される。

ステップＳ１０３：割当部１０３は、１つの特徴ベクトルについて、現在の繰り返し回数をｔとした場合にｔ−１のときに計算されたＬＢ及びＵＢを用いて、ＬＢ≧ＵＢであるか否かを判定する。

ＬＢ≧ＵＢであると判定された場合、当該特徴ベクトルの類似度の計算は行われず、次の特徴ベクトルについてステップＳ１０３の処理が実行される（次の特徴ベクトルが無い場合は、ステップＳ１０９の処理に進む。）。この場合、割当部１０３は、当該特徴ベクトルのＵＢ及びＬＢを更新部１０４に渡す。

一方で、ＬＢ≧ＵＢであると判定されなかった場合（すなわち、ＬＢ＜ＵＢであると判定された場合）、ステップＳ１０４の処理に進む。ただし、ステップＳ１０４の処理を行わずに、ステップＳ１０５の処理に進んで類似度計算を実行しても良い。

ステップＳ１０４：割当部１０３は、当該特徴ベクトルについて、現在の繰り返し回数をｔとした場合に、ｔにおけるｓｉｍを計算する。そして、割当部１０３は、ｔにおけるｓｉｍと、ｔ−１のときに計算されたＵＢとを用いて、ｓｉｍ≧ＵＢであるか否かを判定する。

ｓｉｍ≧ＵＢであると判定された場合、当該特徴ベクトルの類似度の計算は行われず、次の特徴ベクトルについてステップＳ１０３の処理が実行される（次の特徴ベクトルが無い場合は、ステップＳ１０９の処理に進む。）。この場合、割当部１０３は、当該特徴ベクトルのｓｉｍ及びＵＢを更新部１０４に渡す。

一方で、ｓｉｍ≧ＵＢであると判定されなかった場合（すなわち、ｓｉｍ＜ＵＢであると判定された場合）、ステップＳ１０５の処理に進む。

なお、上記のステップＳ１０３の処理を行わずに、上記のステップＳ１０４の処理を実行するようにしても良い。この場合、上記のステップＳ１０２で初期状態の設定が完了した後、割当部１０３は、上記のステップＳ１０４の処理を行えば良い。また、この場合、後述するステップＳ１０９では、更新部１０９は、ＵＢのみを計算すれば良い。

本発明は、上記のステップＳ１０３の処理又は上記のステップＳ１０４の処理のいずれか一方のみが行われる構成であってもその効果を奏する。このため、本発明は、上記のステップＳ１０３の処理及び上記のステップＳ１０４の処理の少なくとも一方が行われれば良い。

ステップＳ１０５：割当部１０３は、当該特徴ベクトルと、各クラスタの平均ベクトルとの類似度を計算（厳密な類似度計算）する。

ここで、例えば文書の特徴ベクトルの各要素はtf-idf（非負値）であるため、図４に示すように、各特徴ベクトルx_iは、超高次元空間における単位超球の第１象限表面上の点となる。このため、特徴ベクトルx_iと平均ベクトルμ_jとのコサイン類似度は、cosθ=x_i ^Tμ_jで表される。なお、jはクラスタＩＤを表す。

ステップＳ１０６：次に、割当部１０３は、上記のステップＳ１０５における類似度の計算結果を用いて、最類似平均ベクトル及び第二類似平均ベクトルを求める。また、割当部１０３は、最類似平均ベクトルのクラスタを識別するクラスタＩＤと、最類似平均ベクトルとの類似度と、第二類似平均ベクトルのクラスタを識別するクラスタＩＤと、第二類似平均ベクトルとの類似度とを記憶部１０７に記憶させる。

更に、このとき、割当部１０３は、当該特徴ベクトルと最類似平均ベクトルとの類似度と、当該特徴ベクトルと第二類似平均ベクトルとの類似度とを更新部１０４に渡す。

ステップＳ１０７：割当部１０３は、当該特徴ベクトルが属するクラスタ（すなわち、現在の繰り返し回数をｔとした場合に、ｔ−１回目の繰り返し終了時点における当該特徴ベクトルが属するクラスタ）の平均ベクトルと、最類似平均ベクトルとが同一であるか否かを判定する。

当該特徴ベクトルが属するクラスタの平均ベクトルと、最類似平均ベクトルとが同一であると判定された場合、次の特徴ベクトルについてステップＳ１０３の処理が実行される（次の特徴ベクトルが無い場合は、ステップＳ１０９の処理に進む。）。この場合、割当部１０３は、当該特徴ベクトルの属するクラスタは変更されない（言い換えれば、割当部１０３は、当該特徴ベクトルの属するクラスタと同一のクラスタに当該特徴ベクトルを割り当てる。）。

一方で、当該特徴ベクトルが属するクラスタの平均ベクトルと、最類似平均ベクトルとが同一でないと判定された場合、ステップＳ１０８の処理に進む。

ステップＳ１０８：割当部１０３は、当該特徴ベクトルを、当該最類似平均ベクトルのクラスタ（このクラスタのクラスタＩＤは記憶部１０７に記憶されている。）に割り当てる。これにより、当該特徴ベクトルが属するクラスタが、当該最類似平均ベクトルのクラスタに変更される。

ステップＳ１０９：更新部１０４は、現在の繰り返し回数がｔにおける平均ベクトルの計算及び更新と、各特徴ベクトルのＵＢ及びＬＢの計算とを行う処理（更新処理）を実行する。このステップＳ１０９の処理の詳細については後述する。

ステップＳ１１０：終了判定部１０５は、終了条件を満たしたか否かを判定する。終了条件を満たしたと判定された場合、ステップＳ１１１の処理に進む。一方で、終了条件を満たしていないと判定された場合、ステップＳ０１０３の処理に戻る。これにより、終了条件を満たすまで、各特徴ベクトルについて、上記のステップＳ１０３〜ステップＳ１０８の処理が実行され、ステップＳ１０９で平均ベクトルの計算及び更新とＵＢ及びＬＢの計算とが実行される。

ステップＳ１１１：出力部１０６は、各特徴ベクトルのクラスタへの割当結果を出力する。なお、出力部１０６は、各クラスタの平均ベクトルを出力しても良い。

以上により、本発明の実施の形態のクラスタリング装置１０は、球面k-means法（SPKM）を用いて、特徴ベクトルの集合（すなわち、オブジェクトの集合）をクラスタリングすることができる。しかも、本発明の実施の形態のクラスタリング装置１０は、各特徴ベクトルについて、ＬＢがＵＢ未満である場合にのみ平均ベクトルとの厳密な類似度を計算する。これにより、類似度の計算回数を低減することができ、SPKMを用いたクラスタリングの計算時間を短縮することができる。

次に、上記のステップＳ１０９における更新処理について、図５を参照しながら説明する。図５は、本発明の実施の形態におけるクラスタリング装置１０が実行する更新処理の一例を示すフローチャートである。

ステップＳ２０１：更新部１０４は、現在の繰り返し回数ｔにおける平均ベクトルの計算及び更新を行う。すなわち、更新部１０４は、各クラスタについて、代表ベクトルを計算し、計算した代表ベクトルの長さが単位超球上の点となるようにL2ノルムで正規化することで、平均ベクトルを計算する。

ステップＳ２０２：更新部１０４は、現在の繰り返し回数ｔにおける各特徴ベクトルのＬＢ及びＵＢを計算する。すなわち、更新部１０４は、各特徴ベクトルについて、当該特徴ベクトルが属するクラスタの平均ベクトルと、当該特徴ベクトルとの類似度の下限値（ＬＢ）を計算する。また、更新部１０４は、各特徴ベクトルについて、当該特徴ベクトルと、第二類似平均ベクトルとの類似度の上限値（ＵＢ）を計算する。これにより、繰り返し回数ｔにおけるＵＢ及びＬＢが特徴ベクトル毎に得られる。

ここで、上限値及び下限値の計算方法ついて説明する。繰り返し回数がｔ−１回目の平均ベクトルをμ_j ^[t-1]、繰り返し回数がｔ回目の平均ベクトルをμ_j ^[t]とした場合に、ｔ−１回目の平均ベクトルμ_j ^[t-1]を主軸とし、ｔ回目の平均ベクトルμ_j ^[t]とｔ−１回目の平均ベクトルμ_j ^[t-1]との成す角を頂点の最大角とする円錐状の形状（扇形回転体）を考える。

図６に示すように、この扇形回転体が第１象限に含まれる場合、当該扇形回転体と特徴ベクトルx_iとが最近傍及び最遠となる単位超球面上の点は、平均ベクトルμ_j ^[t-1]と特徴ベクトルx_iとを含む超平面と当該扇形回転体との交点である。

ｔ−１回目の平均ベクトルをμ_j ^[t-1]が、繰り返し回数がｔのときにδだけ移動したとすれば（すなわち、μ_j ^[t]とμ_j ^[t-1]との成す角の角度をδとすれば）、類似度の上限値は、平均ベクトルμ_j ^[t]が特徴ベクトルx_iに仮想的に最も近付いたときの角度であり、cos(θ-δ)である。同様に、類似度の下限値は、平均ベクトルμ_j ^[t]が特徴ベクトルx_iに最も遠ざかったときの角度であり、cos(θ+δ)である。なお、θは平均ベクトルμ_j ^[t-1]と特徴ベクトルx_iとの成す角の角度である。

したがって、この場合、cosθ=x_i ^Tμ_j ^[t-1]、sinθ=√(1-cos²θ)、cosδ=μ_j ^[t]Tμ_j ^[t-1]、sinθ=√(1-cos²δ)を用いて、類似度の上限値はcos(θ-δ)=cosθcosδ+sinθsinδ、類似度の下限値はcos(θ+δ)=cosθcosδ-sinθsinδにより計算することができる。なお、jとしては、類似度の上限値を計算する場合は第二類似平均ベクトルのクラスタＩＤ、類似度の下限値を計算する場合は当該特徴ベクトルが属するクラスタのクラスタＩＤを用いる。

図７に示すように、扇形回転体が第１象限に含まれない場合も有り得る。この場合、類似度の上限値は、特徴ベクトルx_iと、ｔ回目の平均ベクトルをμ_j ^[t]との成す角の角度が最小となるときである。そこで、|θ-δ|を考えると、δ-θ＞０の場合は、図８（ａ）に示すように、扇形回転体の境界という意味での角度は、特徴ベクトルx_iと、平均ベクトルμ_j ^[t-1]及び特徴ベクトルx_iを含む超平面と交わる座標軸との成す角の角度εである。平均ベクトルは、この扇形回転体の内部の何れかに存在することになるため、特徴ベクトルx_iと一致する可能性がある。従って、類似度の上限値は１となる。

一方で、δ-θ≦０の場合は、図８（ｂ）に示すように、扇形回転体が第１象限に含まれる。したがって、この場合は、図６を用いて説明した通りである。

また、図９に示すように、δが最大角となった場合に、扇形回転体が第１象限に含まれないことが有り得る。この場合、第１象限内における最小の類似度は、最大角δにおける類似度の下限値より大きくなる。すなわち、第１象限内という制約を外した場合の下限値は、実際に起こり得る最小類似度より小さくなくため、緩い下限値（loose lower bound）であるが、下限としての役割を果たすことができるため、この緩い下限値を、類似度の下限値として用いることができる。したがって、この場合、類似度の下限値はcos(θ+δ)とすれば良い。

＜本発明の効果＞
以降では、本発明の効果について説明する。

≪本発明による手法と従来法との比較≫
本発明による手法（これを「BND法（similarity-based bounding method（類似度に基づく上下限法）」と呼ぶ。）と従来法との比較例を図１０に示す。図１０は、同一の或る計算環境下で、SPKMを用いたクラスタリングを行った場合における繰り返し回数と経過時間との関係をプロットした図である。

図１０示す例では、クラスタリングの対象は或る新聞記事データであり、特徴ベクトル数n=1,285,944（記事数）、ユークリッド空間の次元数dim=495,126（使用されている語彙数）、クラスタ数k=10,000である。

図１０に示すように、本発明による手法（BND法）では、従来法と比べて、経過時間が大幅に低減されているのがわかる。

また、SPKMが収束するまでにおける１回の繰り返しあたりの平均経過時間が、従来法では2195.1（秒）、BND法では1402.4であった。このため、経過時間低減率を1.0-(BND法の平均経過時間)/(従来法の平均経過時間)により計算すると、経過時間低減率=0.36であった。

≪コサイン類似度を用いた場合とL2ノルムを用いた場合との比較≫
コサイン類似度を用いた場合（すなわち、本発明による手法（BND法）を用いた場合）と、コサイン類似度の代わりにL2ノルム（ユークリッド距離）を用いた場合との比較例を図１１に示す。図１１は、同一の或る計算環境下で、SPKMを用いたクラスタリングを行った場合における初期状態（Seed(initial state)）と、計算回数の平均削減率（Average reduction rate）との関係をプロットした図である。なお、初期状態(Seed)には、０から２０までの２１通りを用いた。

図１１に示す例では、クラスタリングの対象は或る新聞記事データであり、特徴ベクトル数n=116,905（記事数）、ユークリッド空間の次元数dim=220,114（使用されている語彙数）、クラスタ数k=1,000である。

図１１に示すように、全ての初期状態において、コサイン類似度を用いた場合の方がL2ノルムを用いた場合よりも計算回数の平均削減率が高いことがわかる。平均削減率とは、コサイン類似度の場合は類似度の計算回数削減率の平均値であり、L2ノルムの場合は距離の計算回数削減率の平均値である。計算回数削減率は、以下の式７で表される。

ここで、#itrはSPKMが収束するまでの繰り返し回数（この繰り返し回数は初期状態で異なる。）、#calcは繰り返し回数１回あたりの実際の類似度又は距離の計算回数である。

更に、コサイン類似度を用いた場合に、２１通りの初期状態で平均削減率を平均すると0.258であった。一方で、L2ノルムを用いた場合に、２１通りの初期状態で平均削減率を平均すると0.224であった。これにより、コサイン類似度を用いる場合の方がL2ノルムを用いる場合よりも計算回数の削減率が高い（すなわち、計算効率が良い）ことがわかる。

上述したように、計算効率は劣るものの、本発明は、コサイン類似度の代わりにL2ノルム（ユークリッド距離）を用いることもできる。そこで、コサイン類似度の代わりにL2ノルム（ユークリッド距離）を用いる場合について説明する。

ユークリッド距離を用いる場合、最類似平均ベクトルの代わりに、距離が最小となる最近傍平均ベクトルを求める。同様に、ユークリッド距離を用いる場合は、第二類似平均ベクトルの代わりに、距離が２番目に小さい第二近傍平均ベクトルを求める。

単位超球上にある特徴ベクトルをx_i、平均ベクトルをμ_j、x_iとμ_jとの成す角の角度をθとすると、コサイン類似度は、cosθ=x_i ^Tμ_jとなる。このため、図１２に示すように、コサイン類似度は、x_iとμ_jとの間の弧（arc）の長さに相当する。

一方で、ユークリッド距離は、√(2(1-cosθ))= √(2(1- x_i ^Tμ_j))となる。このため、図１２に示すように、ユークリッド距離は、x_iとμ_jとの間の弦（chord）の長さに相当する。なお、最近傍平均ベクトルを求めるだけの場合は、２乗距離（squared Euclidean distance）で代用することができるため、2(1-x_i ^Tμ_j)が最小になる平均ベクトル（すなわち、x_i ^Tμ_jが最大になる平均ベクトル）を求めることとなり、コサイン類似度を用いる場合と同じ解を得ることとなる。

また、μ_jが最も近づく下限値の場合であって第１象限内に存在する場合を考える。このとき、μ_a ^[t]をｔ回目における平均ベクトルであってx_iが属するクラスタの平均ベクトル、μ_j ^[t]（j≠a）を他のクラスタの平均ベクトル、Δ_a=d(μ_a ^[t-1], μ_a ^[t])、Δ_max=Δ_j=d(μ_j ^[t-1], μ_j ^[t])とする。ユークリッド距離を用いた場合の上限値（ＵＢ）及び下限値（ＬＢ）は、ＵＢ=d(x_i, μ_a ^[t-1])+Δ_a=√(2(1- x_i ^Tμ_a ^[t-1]))-√(2(1-μ_a ^[t]Tμ_a ^[t-1]))、ＬＢ=d(x_i, μ_j ^[t-1])-Δ_maxと表される。

ユークリッド距離を用いた場合の距離計算（厳密な距離計算）の計算回数については、例えば、図１３に示すように、或る領域Ｄを考えた場合に、この領域Ｄに含まれる単位超球面上に平均ベクトルμ_j ^[t]が存在する場合、ユークリッド距離を用いた場合（すなわち、弦を用いた場合）はＵＢ＞ＬＢとなり、厳密な距離計算が必要となる。一方で、コサイン類似度（弧を用いた場合）はＬＢ≧ＵＢとなり、厳密な類似度計算は不要となる。このように、弧と弦との違いによって、ユークリッド距離を用いた場合の厳密な距離計算の回数は、コサイン類似度を用いた場合の厳密な類似度計算の回数と比べて多くなる。このため、計算回数の削減率も、ユークリッド距離を用いた場合はコサイン類似度を用いた場合よりも低くなる。

ところで、非特許文献３に開示されている方法と同様にユークリッド距離における上下限値を用いる方法であって、下限値をより厳しく（tighter）する方法として、非特許文献４に開示されている方法や以下の参考文献１に開示されている方法がある。

［参考文献１］
J. Drake and G. Hamerly, "Accelerated k-means with adaptive distance bounds," Proc. 5^th NIPS Workshop on Optimization for Machine Learning, 2012.
非特許文献３に開示されている方法は、ユークリッド距離を用いた一般的なk-meansクラスタリング法において、或るデータ点を表す特徴ベクトルが所属するクラスタ以外のクラスタの平均ベクトルのうち、最も近い平均ベクトル（すなわち、２番目に近い平均ベクトル）までの距離から、(ｔ−１)回目からｔ回目のiteration（繰り返し）で最も移動した平均ベクトルの最大変化距離を減算することにより、下限値を算出する。

参考文献１に開示されている方法は、２番目に近い平均ベクトルではなく、ｂ＋１番目に近い平均ベクトルから同様に最大変化距離を減算し、２番目からb番目までに近い平均ベクトルに対しては、自身の移動距離を減算する。下限値を確実なものとするために、２番目からｂ番目までに近い平均ベクトルに基づいて算出された下限値が、ｂ＋１番目に近い平均ベクトルに基づいて算出された下限値よりも大きい距離となっているものは、ｂ＋１番目に近い平均ベクトルに基づいて算出された下限値に置換される。この方法のユークリッド距離を用いて算出された下限値を本発明の類似度に基づく下限値に置き換えることもできる。

非特許文献４に開示されている方法は、平均ベクトルをＴ個のグループに分割し（平均的にはｋ／Ｔ個が１つのグループになる）、或る特徴ベクトルに対して、各グループの中で最も近い平均ベクトルからそのグループ内での最大変化距離を減算することで各グループにおける下限値を算出する。ただし、特徴ベクトルが所属するクラスタに関しては最も近い平均ベクトルではなく、２番目に近いものを基に下限値が算出される。このユークリッド距離を用いた方法の上下限値算出を、本発明の類似度に基づく上下限値を算出する方法に置き換えることもできる。

本発明は、具体的に開示された上記の実施形態に限定されるものではなく、特許請求の範囲から逸脱することなく、種々の変形や変更が可能である。

１０クラスタリング装置
１０１入力部
１０２初期設定部
１０３割当部
１０４更新部
１０５終了判定部
１０６出力部
１０７記憶部

Claims

球面上にそれぞれ配置された複数の特徴ベクトルを繰り返し手続きにより複数のクラスタに割り当てるクラスタリング装置であって、
繰り返しの回数を表す変数をｔとして、ｔ回目において、前記特徴ベクトル毎に、前記特徴ベクトルと、前記複数のクラスタそれぞれの代表ベクトルとの類似度を計算する類似度計算手段と、
ｔ回目において各特徴ベクトルに対して前記類似度の計算を行うか否かをそれぞれ判定する判定手段と、
ｔ回目において各特徴ベクトルに対して前記類似度の上限値をそれぞれ算出する上限算出手段と、
を有し、
前記上限算出手段は、
前記代表ベクトルのうち、前記特徴ベクトルと２番目に類似する代表ベクトルを第二類似代表ベクトルとして、
前記第二類似代表ベクトルと同一クラスタのｔ−１回目における代表ベクトルと前記特徴ベクトルとの第一の成す角と、前記ｔ−１回目における代表ベクトルと前記第二類似代表ベクトルとの第二の成す角との２つの成す角を用いて、前記第二類似代表ベクトルが前記特徴ベクトルに仮想的に最も近づいたときの類似度を、ｔ回目における上限値として算出し、
前記判定手段は、
ｔ−１回目における前記上限値を少なくとも用いて、ｔ回目において各特徴ベクトルに対して前記類似度の計算を行うか否かを判定し、
前記類似度計算手段は、
前記判定手段により前記類似度の計算を行わないと判定された特徴ベクトルについては、ｔ回目において前記類似度の計算を行わない、ことを特徴とするクラスタリング装置。
前記判定手段によりｔ回目において類似度の計算を行うと判定された特徴ベクトルについて、前記類似度計算手段により計算された類似度に基づいて、前記特徴ベクトルを前記複数のクラスタのうちのいずれかのクラスタに割り当てる割当手段と、
ｔ回目において、前記クラスタ毎に、前記割当手段により前記クラスタに割り当てられた特徴ベクトルから、ｔ＋１回目における各クラスタの代表ベクトルを計算する代表ベクトル計算手段と、
を有することを特徴とする請求項１に記載のクラスタリング装置。
ｔ回目において各特徴ベクトルに対して前記類似度の下限値をそれぞれ算出する下限算出手段を有し、
前記下限算出手段は、
前記代表ベクトルのうち、前記特徴ベクトルが属するクラスタの代表ベクトルを所属代表ベクトルとして、
前記所属代表ベクトルと同一クラスタのｔ−１回目における代表ベクトルと前記特徴ベクトルとの第一の成す角と、前記ｔ−１回目における代表ベクトルと前記所属代表ベクトルとの第二の成す角との２つの成す角を用いて、前記所属代表ベクトルが前記特徴ベクトルに仮想的に最も遠ざかったときの類似度を、ｔ回目における下限値として算出し、
前記判定手段は、
ｔ−１回目における前記上限値が、ｔ−１回目における前記下限値以下でなく、かつ、ｔ−１回目における前記上限値が、前記特徴ベクトルと該特徴ベクトルが所属する代表ベクトルとの類似度以下でない場合、ｔ回目において前記類似度の計算を行うと判定し、
ｔ−１回目における前記上限値が、ｔ−１回目における前記下限値以下の場合、又は、ｔ−１回目における前記上限値が、前記特徴ベクトルと該特徴ベクトルが所属する代表ベクトルとの類似度以下の場合、のいずれか一方を満たす場合、ｔ回目において前記類似度の計算を行わないと判定する、ことを特徴とする請求項１又は２に記載のクラスタリング装置。
前記類似度計算手段は、
前記類似度として、コサイン類似度又はＬ２ノルムを計算する、ことを特徴とする請求項１乃至３の何れか一項に記載のクラスタリング装置。
球面上にそれぞれ配置された複数の特徴ベクトルを繰り返し手続きにより複数のクラスタに割り当てるコンピュータが、
繰り返しの回数を表す変数をｔとして、ｔ回目において、前記特徴ベクトル毎に、前記特徴ベクトルと、前記複数のクラスタそれぞれの代表ベクトルとの類似度を計算する類似度計算手順と、
ｔ回目において各特徴ベクトルに対して前記類似度の計算を行うか否かをそれぞれ判定する判定手順と、
ｔ回目において各特徴ベクトルに対して前記類似度の上限値をそれぞれ算出する上限算出手順と、
を実行し、
前記上限算出手順は、
前記代表ベクトルのうち、前記特徴ベクトルと２番目に類似する代表ベクトルを第二類似代表ベクトルとして、
前記第二類似代表ベクトルと同一クラスタのｔ−１回目における代表ベクトルと前記特徴ベクトルとの第一の成す角と、前記ｔ−１回目における代表ベクトルと前記第二類似代表ベクトルとの第二の成す角との２つの成す角を用いて、前記第二類似代表ベクトルが前記特徴ベクトルに仮想的に最も近づいたときの類似度を、ｔ回目における上限値として算出し、
前記判定手順は、
ｔ−１回目における前記上限値を少なくとも用いて、ｔ回目において各特徴ベクトルに対して前記類似度の計算を行うか否かを判定し、
前記類似度計算手順は、
前記判定手順により前記類似度の計算を行わないと判定された特徴ベクトルについては、ｔ回目において前記類似度の計算を行わない、ことを特徴とするクラスタリング方法。
コンピュータを、請求項１乃至４の何れか一項に記載のクラスタリング装置における各手段として機能させるためのプログラム。