(57)【要約】【課題】本発明の目的は、一般化された公式を用いて大きなデータセットにおけるクラスタリングを実行する新規の方法を提供することである。【解決手段】我々は、（１）最適輸送距離に基づいてヒストグラムのセットの重心（つまり平均）を計算する新規な方法、および（２）第１の方法の拡張として、クラスタのサイズに任意の制約を持つRｄにおけるベクトルのデータセットをクラスタリングする効果的かつ新規な方法を提供する。【選択図】図７
2015203946.tif 000002
【特許請求の範囲】
【請求項１】
近似した双対および主最適変数の合計を計算する方法であって、
ｎ個の変数のシンプレックスΣｎにおけるＮ個のヒストグラムのデータセット{c1, c2, …, cN}をｎ行Ｎ列の行列Ｃに格納して、収束許容誤差TOLを設定するステップと、
ｄ行ｄ列の行列ＫおよびＱを計算するステップであって、その要素(i,j)は、kij = exp(-λ mij ); qij = mij exp(-λ mij )に等しいステップと、
ｄ行Ｎ列の行列Ｕを初期化するステップであって、Ｕの各要素は１／ｄに等しいステップと、
行列Ｌ＝diag(1./r)Ｋを計算して、ｚ＝無限大に設定するステップと、
ｚ<TOLになるまでａ）およびｂ）を反復するステップであって、ａ）およびｂ）は、
ａ）U= 1./ (L (C./(K’U)))を計算するステップ、
ｂ）i. V=C./(K’U); U=1./ (LV)を求め、ii. 出口条件z =||V.*(K’U)-C)||を更新することを所定の回数反復するステップ、
であるステップと、
合計された近似目的、双対最適条件ρ★λおよび主最適条件Ｔ★λを計算するステップであって、
【数１】
2015203946.tif 000034
【数２】
2015203946.tif 000035
【数３】
2015203946.tif 000036
であるステップと、
を有する方法。
【請求項２】
Ｎ個のヒストグラムのワッサーステイン重心を計算する方法であって、
ｎ個の変数のシンプレックスΣｎにおけるＮ個のヒストグラムのデータセット{c1, c2, …, cN}、およびｎ行ｎ列の行列Ｍを収集するステップと、
Σｎの関連部分集合Θを当該部分集合上の投影PΘとともに定義するステップあって、該投影は、所与の任意のベクトルａを、前記ポイントのΘ内の最も近いポイントに戻す関数であるステップと、
ｒをベクトル[1/n, 1/n, ..., 1/n]に初期化するステップと、
所望に収束するまでｃ）〜ｆ）を反復するステップと、
を有し、ｃ）〜ｆ）は、
ｃ）Ｎ個の双対問題{dλM(r,c1), dλM(r,c2), …, dλM(r,cN)}を解いて、後述のサブルーチンを用いて、Ｎ個の距離diとＮ個の双対最適変数ρi★λを回復するステップ、
ｄ）請求項１に記載の方法を用いて、近似目的および近似勾配を形成するステップ、
【数４】
2015203946.tif 000037
ｅ）現在の変数r ← PΘ (r - ερ)を更新するステップ、
ｆ）２つの連続した反復間の目的における絶対的な差があらかじめ定義された許容値を下回った場合に停止するステップ、
である方法。
【請求項３】
Σkの部分集合Θに制約されている重み付けを有するRdにおける経験的計測のワッサーステイン重心を計算する方法であって、
ｎ個の変数のシンプレックスΣnの重みベクトルｂを有する、Rdにおけるn個の重み付けされたポイント群{x1, x2, …, xn}を合計するステップであって、前記ポイントは、ｄ行ｎ列の行列Ｘとして表すことができるステップと、
Σｋの関連する部分集合Θを該部分集合上の投影PΘとともに定義するステップであって、該投影は、所与の任意のベクトルａを、前記ポイントのΘ内の最も近いポイントに戻す関数であるステップと、 
ＹをＤ行ｋ列の行列に初期化するステップであって、各列が、Ｘの列の中から無作為に抽出されるステップと、
ａをベクトル[1/k, 1/k, ..., 1/k]に設定するステップと、
所望に収束するまでｇ）〜ｊ）を反復するステップと、
を有し、ｇ）〜ｊ）は、
ｇ）距離行列
【数５】
2015203946.tif 000038
を形成するステップ、
ｈ）請求項２に記載の方法を用いて最適な重みａを計算するステップであって、Ｍを距離行列パラメータとして、ｂを入力ヒストグラム（Ｎ＝１）として用いるステップ、
ｉ）請求項１に記載の方法を用いて近似的な最適輸送Ｔ★λを計算するステップ、
ｊ）勾配ステップ：Ｙ← Ｘ Ｔ★T diag(b-1)を更新するステップ、
である方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ヒストグラムのセットの重心を計算する方法に関し、特に、Ｎのヒストグラムのワッサースタイン重心、近似した双対および主最適変数の合計、および経験的計測のワッサーステイン重心を計算する方法に関する。
【背景技術】
【０００２】
非常に大きなデータセットを効率的に簡素化し要約することは、データ解析アルゴリズムの基本的な処理である。そのような処理を実行するための２つの強力なツールが、前処理ステップとして日常的に使用されている。前処理ステップとは、データセットの平均を計算するステップ、および、ｋ平均問題（ｋは１より大きい整数）として知られている、データセットの複数の平均を計算するステップである。
【０００３】
（データセットの平均）
データセットの特徴の基本的な概要は、当該データセットの平均によって与えられる。データセットの平均は、当該データセット内のすべてのポイントへの距離の合計が最小となるポイントである。例えば、任意の集合Ωから得られたＮ個のポイントの集合をＸ ={x1、x2、...、xN}、１以上の整数をpとし、Ωにおける２つの対象を比較する距離関数をdistとすると、Ｘの平均は、Ωにおける次の量を最小化する対象ｙである。
【０００４】
【数１】
2015203946.tif 000003
【０００５】
データセット内のポイントが、ｄ次元のベクトル空間Rdにおけるベクトルであり、これらのベクトル間のユークリッド距離が考慮され、さらに、指数ｐ＝２とすると、Rdにおけるそれらの平均は、単純に、その座標がデータセット内のすべてのポイントの座標の平均値となるポイントである。すなわち、
【０００６】
【数２】
2015203946.tif 000004
【０００７】
である。
【０００８】
（データセットのｋ平均）
１つの平均の概念の自然な拡張は、１の平均要素だけでなく、多数の平均要素を用いてデータセットをまとめることである。ｋが１以上の整数であり、正確なｋ平均を求めることは、ｋ平均問題として知られている。各ΩにおけるＮポイントのデータセットを{x1、x2、...、xN}とすると、標準的なk平均アルゴリズムは、Ωにおける、以下の基準を最小化するｋ平均ポイント{y1、y2、...、yk}を求めようとする。
【０００９】
【数３】
2015203946.tif 000005
【００１０】
これらのｋ平均も重心と呼ばれている。ｋ重心{y1、y2、...、yk}のそれぞれは、Ωにおけるポイントである。最適なｋ重心を求めることは、集合{y1、y2、...、yk}内の各ポイントxiとその最短の隣接ポイントとの間の平均残留誤差を最小限にする重心を求めるという問題である。最適化問題に現れる、{1、...、k}におけるサイズＮのベクトルσは、属性ベクトルと呼ばれる。そのベクトルのｉ番目の値であるσiは、ｋ重心ポイントxiがそれに属すべきであることを示している。
【００１１】
図１（ａ）は、元データセットを示す図であり、図１（ｂ）〜（ｄ）はそれぞれ、単一平均の計算結果、ｋ＝２のｋ平均の計算結果、およびｋ＝３のｋ平均の計算結果を示す図である。この例では、Ω＝R2であり、距離関数はユークリッド距離であり、ｐ＝２である。図１（ａ）〜（ｄ）において、交点は平面R2におけるポイントを表している。図１（ｂ）の四角は、これらの交点の通常の平均値を表している。図１（ｃ）の２つの四角は、ｋが２に設定されている場合のｋ平均法の結果を表している。図１（ｄ）の３つの四角は、ｋが３に設定されている場合のｋ平均法の結果を表している。１番目の図１（ｂ）の四角を計算することは簡単であり、すべての点の座標を合計して、その合計をポイントの個数によって割るだけでよい。ｋ＞１である場合のｋ平均を計算することは、ＮＰ困難である問題であることが知られている。
【００１２】
そこで、開示のシステムは、依然として有用な設定を試みるクラスタリングを行うために、二つのアルゴリズムによる解決法を提案する。第１の方法は、最適輸送距離に基づいてヒストグラムのセットの単一平均を計算する。第１の方法を拡張して、我々は、クラスターサイズ、または等価的には属性ベクトルにおいて追加的な制約を考慮した場合の、ｋ平均アルゴリズムの解法を計算するための第２の方法を提案する。
【００１３】
（最適輸送距離におけるヒストグラムのための平均の高速計算）
我々が提案する第１の方法は、最適輸送距離（ワッサースタイン距離としても知られる）を距離関数distとして用いて、ヒストグラムのセットの平均を計算することができる。このような平均はワッサー重心として文献で知られている。図１（ｂ）に示すように、ベクトルのセットの単一平均を計算することは簡単な処理であるが、ここでは我々は挑戦的な事例を想定する。すなわち、観察点（図１（ｂ）の交点）が、ユークリッド空間における単純なベクトルではなく、ヒストグラムである事例である。さらに、我々は次の事例も想定する。すなわち、距離関数distが、適切な距離行列Ｍによってパラメータ化されたヒストグラム間の最適輸送距離である事例である。標準的なユークリッド平均と最適輸送距離によって導出された平均との違いを、図２および図３に示す。図２において、我々は、３０の画像セットのワッサースタイン重心（それだけでなく、異なる距離を使用する他の重心）を計算する。３０の画像は、画像内の100×100 = 10,000の可能なピクセル位置の各々における輝度のヒストグラムを表している。
【００１４】
具体的には、図２は二重楕円の３０の画像を示している。図３（ａ）〜（ｄ）は、図２に示す３０の画像の、異なる距離および前処理アプローチを用いた異なる平均を示している。図３（ａ）は、ユークリッド距離を用いた通常の算術による平均を示している。図３（ａ）に示すように、ベクトルのセットの単一平均を計算することは簡単な処理であるが、ここでは我々は、各観察点（図３（ａ）における交点）がユークリッド空間における単純なベクトルではなく、ヒストグラムである挑戦的な事例を想定する。図３（ｂ）は、各画像をセンタリングした後のユークリッド距離を用いたものを示している。図３（ｃ）は、ジェフリー重心を示しており、図３（ｄ）は、ＲＫＨＳ平均を示している。その設定に使用される距離パラメータは、100×100ピクセルグリッドのピクセル間の標準的なユークリッド距離である。
【００１５】
（制約されたｋ平均法）
我々が想定する第２の問題は、ｋ平均法の一般化であり、それは、各クラスタによって行われる総質量の制約を考慮することができる。実際には、ｋ平均アルゴリズムの結果は、属性ベクトルσが不平衡分布の値をとるという意味で、不均衡であり得る。例えば、ある人は、ｋ平均アルゴリズムの結果において、殆どのポイントがｋ個の重心の１つのみに属することができ、残りのｋ−１の重心は、少数の残りのポイントのために用いられる、と想像するかもしれない。このような状況は、例えば図１（ｄ）に観察される。図１（ｄ）では、３のポイントのみが最も小さい四角で示す重心に属する。これに対し、図１（ａ）に示す元データセットにおける２５のポイントは、最も大きい四角で示す重心に属する。
【００１６】
図４は、米国の４８の隣接州における平均所得の分布を示している。５７６４７の交点のそれぞれは、データセットにおけるデータポイントを表す。交点の大きさ（および明暗度）は、その地点で記録された平均所得の値に比例している（明るい交点であるほど所得が高いことを示す）。図５において、標準的なｋ平均法（ここではｋ＝４８）の結果が、逆三角形を用いて示されている。
【００１７】
〔数学的な定義〕
（ヒストグラム）
ベクトルuは、任意の次元dを有し、合計１になる非負の座標を有する。次元ｄのヒストグラムΣｄのセットは、次のように定義される。
【００１８】
【数４】
2015203946.tif 000006
【００１９】
図６は、３成分のヒストグラムのセットΣ３を、２つの例示的なベクトルｒおよびｃに沿って示している。ベクトルｒおよびｃは、それらの合計が１であり、非負の要素を有するため、Σ３内に存在する。
【００２０】
（重み付けされたポイント群）
空間Ωにおけるポイントの有限数ｎのファミリーは、すべての重みの合計が１に等しくなるように、非負の重みが互いに関連付けられている。すなわち、
【００２１】
【数５】
2015203946.tif 000007
【００２２】
また、重み付けされたポイント群は、しばしば、ディラック記法δｘを用いた確率測度として表現される。ディラック記法δｘは、他の位置における{Ｘ}と０のセットにおいて、質量１を有する確率測度を表す。例えば、ポイント群
【００２３】
【数６】
2015203946.tif 000008
【００２４】
は、測度表記で等価的には以下のように表すことができる。
【００２５】
【数７】
2015203946.tif 000009
【００２６】
（ヒストグラムのための最適輸送距離）
Σｎ におけるｒ、および、Σｍ におけるｃの２つのヒストグラム、並びに、ｎ行ｍ列のコスト行列Ｍ＝[ｍｉｊ]が所与であり、ｒとｃの間の最適輸送距離（つまりワッサースタイン距離）ｄＭが、ｎ行ｍ列の非負の行列である変数Ｔ＝[ｔｉｊ]ｉｊを有する以下の線形プログラムの結果として定義される。
【００２７】
【数８】
2015203946.tif 000010
【００２８】
上記の式を単純化するために、我々はまた、同じサイズの二つの行列のために以下の表記も使用する。
【００２９】
【数９】
2015203946.tif 000011
【００３０】
（Ωにおける重み付けされたポイント群のための２つのワッサーステイン距離）
２つの重み付けされたポイント群μ,νを、
【００３１】
【数１０】
2015203946.tif 000012
【００３２】
とする。また、１より大きい数をｐとすると、μとνの間のワッサーステイン距離Ｗｐ(μ,ν)は、ヒストグラムのための上記式によって、次のように直接的に規定される。
【００３３】
【数１１】
2015203946.tif 000013
【００３４】
この式において、ｎ行ｍ列の行列ＭｐＸＹは、ＸおよびＹで表されるポイント間の距離によって、次のように与えられる。
【００３５】
【数１２】
2015203946.tif 000014
【００３６】
（ヒストグラムのワッサースタイン重心）
ｎ個の変数のシンプレックスΣｎにおけるヒストグラムのファミリー{c1, c2, …, cN}、および、各ヒストグラムの全てのｎ個のビン間の対の距離を計算するためのｎ行ｎ列の行列Ｍが所与である場合、これらのヒストグラムのワッサースタイン重心の計算問題は、
【００３７】
【数１３】
2015203946.tif 000015
【００３８】
を最小化するヒストグラムｒを求める問題である。
【００３９】
なお、この定義は、任意の距離関数のための実施形態１において我々が提案する一般的な定義と整合する。
【００４０】
〔先行研究〕
（ヒストグラムのためのワッサースタイン重心の計算）
AguehとCarlier（２０１１年）は、任意の確率測度のためのワッサースタイン重心の枠組みを初めて導入し、その理論的な特性の一部を説明した（非特許文献１）。しかし彼らは、ポイント群、ヒストグラム又は観察によって立証できる測度でそれらを計算するための実用的なアプローチは提案していなかった。その後、Rabinら（２０１２年）は、ワッサースタイン距離の近似の最小化に基づいたアルゴリズムのアプローチを提案した。彼らは、それをスライスされたワッサースタイン距離と名付けた（非特許文献２）。彼らの方法は、距離がユークリッドであり、セットΩが低次元のユークリッド空間であると仮定している。我々の方法では、これらの制約はなく、それゆえ、ＢｏＷ表現またはＢｏＶＦ表現といった、構造化されたデータタイプを表現するヒストグラムに適用することができる。
【００４１】
（クラスタのサイズに任意の制約を持つ一般化されたｋ平均）
クラスターサイズに制約のないｋ平均法は、早くとも５０年代初頭にLloydが提案していたが、公に出版されたのは１９８２年である。 Ng（２０００年）は、図４に示したように、重みに均一な制約のあるｋ平均問題を研究することを提案した（非特許文献４）。Ngのアルゴリズムのアプローチは、最適輸送の系統だった演算に依存しており、計算上集約的である。また、重みが一意に固定されている場合にしか適用することができない。
【先行技術文献】
【非特許文献】
【００４２】
【非特許文献１】Agueh, M., Carlier, G., 「Barycenters in the Wasserstein space」, SIAM Journal on Mathematical Analysis, 2011年, 第43巻第2号, 904-924頁
【非特許文献２】Rabin, J., Peyre, G., Delon, J., & Bernot, M., 「Wasserstein barycenter and its application to texture mixing」, Scale Space and Variational Methods in Computer Vision, Springer Berlin Heidelberg, 2012年, 435-446頁
【非特許文献３】Lloyd, Stuart P., 「Least squares quantization in PCM」, IEEE Transactions on Information Theory, 1982年, 第28巻第2号, 129-137頁
【非特許文献４】Ng, M. K., 「A note on constrained k-means algorithms」, Pattern Recognition, 2000年, 第33巻第3号, 515-519頁
【発明の概要】
【発明が解決しようとする課題】
【００４３】
本発明は、上記課題に鑑みてなされたものであり、ヒストグラムの平均要素を容易に計算するために最適化された装置を提供することを目的としている。さらに拡張として、本発明は、各クラスタの重みに凸面制約があるクラスタリングを実行するためのアルゴリズムを提案する。
【課題を解決するための手段】
【００４４】
上記課題を解決するため、我々は、請求項に係る方法、すなわち、後述のアルゴリズム１〜３を提案する。
【発明の効果】
【００４５】
これらの方法によれば、最適輸送距離におけるヒストグラムのセットの重心（平均）を容易に計算することができ、クラスタのサイズに制約があるRdにおけるベクトルのデータセットをクラスタリングすることができる。
【図面の簡単な説明】
【００４６】
【図１】（ａ）は、元データセットを示す図であり、（ｂ）〜（ｄ）はそれぞれ、単一平均の計算結果、ｋ＝２のｋ平均の計算結果、およびｋ＝３のｋ平均の計算結果を示す図である。
【図２】二重楕円の３０の画像を示す図であり、各画像は、ピクセルの合計輝度が１に標準化された100×100 ピクセルイメージである。
【図３】（ａ）〜（ｄ）は、図２に示す３０の画像の、異なる距離および前処理アプローチを用いた異なる平均を示す図であり、（ａ）は、ユークリッド距離を用いた通常の算術による平均を示しており、（ｂ）は、各画像をセンタリングした後のユークリッド距離を用いたものを示しており、（ｃ）は、ジェフリー重心を示しており、（ｄ）は、ＲＫＨＳ平均を示している。
【図４】米国の４８の隣接州における平均所得の分布を示す図である。
【図５】標準的なｋ平均法（ここではｋ＝４８）の結果が逆三角形を用いて示された図である。
【図６】ｄ＝３におけるヒストグラムのセットΣｄを示す図である。
【図７】本発明に係る距離および前処理アプローチを用いた、図２に示す３０の画像の平均を示す図である。
【図８】円を用いて示された本発明に係るクラスタリング方法の結果である。
【発明を実施するための形態】
【００４７】
本発明は、２つの部分から構成されている。まず、実施形態１において、Ｎ個のヒストグラムのセットのワッサースタイン重心を効率的に計算するための方法を詳細に説明する。続いて、実施形態２において、各クラスタの重み分布がシンプレックスの部分集合に存在するように制約されているクラスタリングを計算するための本発明の方法を説明する。
【００４８】
［実施形態１：重心の重みの制約を伴う特徴間の距離が所与であるＮ個のヒストグラムのワッサースタイン重心の高速計算］
この方法は次のように構築される。すなわち、Σｎ におけるｒ、および、Σｍ におけるｃの２つのヒストグラム、並びに、ｎ行ｍ列のコスト行列Ｍ＝[ｍｉｊ]が所与である場合、ｒとｃの間の最適輸送距離（つまりワッサースタイン距離）dＭは、以下のように最初に定義される。
【００４９】
【数１４】
2015203946.tif 000016
【００５０】
dＭ(r,c)は、全く同じ最適値を有する、双対最適輸送問題として知られている、最適化の双対公式を用いて以下のように計算することができる。
【００５１】
【数１５】
2015203946.tif 000017
【００５２】
ｃとＭが固定パラメータであることを考慮するならば、この双対式は、ｒとｃとの間の輸送距離が、ｒの凸状で区分的な線形関数であることを明確に示している。さらに、その解が唯一である場合、上記の式を最大化する最適な双対ベクトルρ★は、dＭ(r,c)の勾配であり、その解が唯一でない場合、劣勾配であることが知られている。
【００５３】
この公式は、重要かつ実用的な結果をもたらす。すなわち、あらゆるヒストグラムｒについて、ベクトルρ★ の小さな正の分数ε（εは０よりもはるかに大きい）をｒから減算した場合のベクトル(r - ερ★)が、（射影演算子Pｎを用いて）再度Σｎにおいて予測することにより、Σｎに存在することを確認する。そして、Pｎ(r - ερ★) がｃよりもｒに近くなるくらい十分に小さいεが、元来、
【００５４】
【数１６】
2015203946.tif 000018
【００５５】
であることを見積もることができる。
【００５６】
要約すると、ｃとＭが所与である場合、dＭ(r,c) が最小であるようなヒストグラムｒを、単にr ← Pｎ(r - ερ★)の演算を反復することによって、（劣）勾配降下アプローチを用いて求めることができる。当然、ρ★は各反復において再計算される。従って、Ｎ個のヒストグラムのファミリー{c1, c2, …, cN}のワッサーステイン重心を導入するために定義した目的関数である
【００５７】
【数１７】
2015203946.tif 000019
【００５８】
もまた、凸状で区分的な線形であり、（劣）勾配降下アルゴリズムを用いて最小化することができる。（劣）勾配降下アルゴリズムは、各反復において、r ← Pｎ(r - εΣiρi★)の演算を単純に実行するだけである。ここで、ρi★ の各々は、双対最適輸送式によってdＭ(r,ci)を計算して得られる最適双対変数である。
【００５９】
このアルゴリズムは非常に単純である。しかし、実際に実行するのは非常に高コストである。このアルゴリズムは、勾配降下の各段階で、Ｎ個の最適双対変数の計算に依存している。つまり、勾配降下の各反復ごとに、データセット内の各ヒストグラムciのための最適変数ρi★のベクトルを計算する必要がある。比較されるヒストグラムｒおよびciと行列Ｍが任意である一般的な事例では、最も効率的な最適輸送解は、操作の超立方体番号O(n3log(n)) を必要とする。そのため、ヒストグラムの次元ｎが大きい場合には、そのような最適双対変数の計算は非常に高コストになり得る。
【００６０】
この問題を解決するために、我々は、双対問題に現れる制約の平滑化された公式を用いて、双対輸送問題を近似する方法を提案する。インデックス（i、j）のすべてのペアごとにmij - ρi -γj ≧ 0であることを要求するのではなく、我々は、数λ>0が非常に大きい場合は常に、非常に急速に負になる急峻な負のペナルティ-exp(-λ (mij - ρi+γj )) を対象に追加することを選択する。
【００６１】
【数１８】
2015203946.tif 000020
【００６２】
問題dλM(r,c)を解くことは、従来の問題dM(r,c)を解くことよりもはるかに単純である。つまり、平滑化された問題に対する解ρ★λは、シンクホーン行列スケーリングアルゴリズムを用いて回復可能であることを示すことができる。この方法は、近似的な双対最適解を計算し、それを例えば降下方向として用いることにより、シンプレックスΣｎまたはΣｎの関連する部分集合Θにおいてそれを投影する前に、現在の反復で変数ｒを修正することを提案する。
【００６３】
本発明のワッサースタイン重心アルゴリズムは、１つの距離だけでなく、最小にされる必要があるＮ個の距離の合計も考慮するので、本発明に係る方法の本質的な貢献の１つは、全てのＮ個の近似的な双対最適解の合計を計算して、それらを変数ρ★λに格納するための効率的な方法を、アルゴリズム２において提供することである。この近似的な双対最適解の合計を計算することによって、さらに、同じアルゴリズムで、そして、追加コストなしに、輸送問題のすべての第一の最適解の合計の近似を高速に回復することができることを、我々は示すことができる。これはまた、アルゴリズム３を用いる場合に有用である。
【００６４】
（特徴点１）
アルゴリズム２の重要な側面は、このアルゴリズムが、線形代数、及び、より正確な行列積演算（要素的な重複と逆転だけでなく）を含むだけであるということである。したがって、アルゴリズム２の計算は、グラフィック処理装置（ＧＰＵ）で簡単に実行することができ、その結果、グラフィックスカードの安価な演算力を利用することができる。
【００６５】
（特徴点２）
アルゴリズム１の重要な特徴は、合計が同一ではないヒストグラムの処理を容易に一般化できるということである。輸送距離の合計が同一ではない２つの可能な測定値を比較する標準的な方法は、Ωにおける他の全てのポイントに対して固定距離Δ＞０を有する、Ωにおける仮想点ωを設定し、さらに、他の全てのポイント間の絶対的な差に等しいωの重みの最小の合計を測定に加えることである。以下では、我々の方法をこの事例に一般化するアイデアについて説明する。
【００６６】
実際には、異なる合計を有するｎ個のビンの２つのヒストグラムｒ、ｃを、ｎ行ｎ列の距離行列Ｍを用いて比較する場合に、この一般化は、以下のアプローチを適用することによって実現される。（最適輸送距離の対称性による）一般化を欠くことなく、ヒストグラムｒの総合計がヒストグラムｃの総合計より小さいと仮定すると、Ｍによってパラメータ化されたｒからｃまでの最適輸送距離は、Ｍ’によってパラメータ化されたｒ’とｃとの間の最適輸送距離として定義される。なお、（１）Ｍ’はｎ＋１行ｎ列の行列であり、一様にΔに等しい長さｎの一定の列ベクトルがその底に付加された行列Ｍと等しい。（２）ｒ’は、ｎ＋１に等しい構成要素を有するヒストグラムであり、ｒ’の最初のｎ個の要素は、ｒのそれらと等しく、一方、最後の要素は、ｒのｃ個の負の要素の合計と異なる。ｒ’の定義によって、ｒ’の合計とｃの合計とが同一になることは、容易に認識できる。標準化されていないヒストグラムに適用されるこの定義を用いることにより、アルゴリズムの全てのステップにおいて必要な場合はいつでも、それと比較されるヒストグラムに応じて、考慮されたヒストグラムをその標準化されたバージョンに置き換えることによって、必ずしも合計が同一でないＮ個のヒストグラムの重心を、アルゴリズム１によって計算することは容易になる。
【００６７】
以下、我々のアルゴリズムを段階的に説明する。
【００６８】
（アルゴリズム１：Ｎ個のヒストグラムのワッサーステイン重心）
１．ｎ個の変数のシンプレックスΣｎにおけるＮ個のヒストグラムのデータセット{c1, c2, …, cN}、およびｎ行ｎ列の行列Ｍを収集する。
２．Σｎの関連部分集合Θを当該部分集合上の投影PΘとともに定義する。投影は、任意のベクトルyが所与である場合、Θにおいて最も近い点を返す関数である。
３．ｒをベクトル[1/n, 1/n, ..., 1/n]に初期化する。
４．所望に収束するまでａ〜ｄを反復する。
ａ．Ｎ個の双対問題{dλM(r,c1), dλM(r,c2), …, dλM(r,cN)}を解いて、後述のサブルーチンを用いて、Ｎ個の距離diとＮ個の双対最適変数ρi★λを回復する。
ｂ．アルゴリズム２を用いて、近似目的および近似勾配を形成する。
【００６９】
【数１９】
2015203946.tif 000021
【００７０】
ｃ．現在の変数r ← PΘ (r - ερ)を更新する。
ｄ．２つの連続した反復間の目的における絶対的な差があらかじめ定義された許容値を下回った場合に停止する。
【００７１】
（アルゴリズム２：近似した双対および主最適変数の合計の計算）
１．ｎ個の変数のシンプレックスΣｎにおけるＮ個のヒストグラムのデータセット{c1, c2, …, cN}をｎ行Ｎ列の行列Ｃに格納する。収束許容誤差TOLを設定する。
２．ｄ行ｄ列の行列ＫおよびＱを計算する。その要素(i,j)は、kij = exp(-λ mij ); qij = mij exp(-λ mij )に等しい。
３．ｄ行Ｎ列の行列Ｕを初期化する。Ｕの各要素は１／ｄに等しい。
４．行列L=diag(1./r) Kを計算する。ｚ＝無限大に設定。
５．ｚ<TOLになるまでａおよびｂを反復する。
ａ．U= 1./ (L (C./(K’U)))
ｂ．iおよびiiを１０回程度反復する。
i.    V=C./(K’U); U=1./ (LV)を求める。
ii.   出口条件z =||V.*(K’U)-C)||を更新する。
６．合計された近似目的、双対最適条件ρ★λおよび主最適条件T★λを計算する。
【００７２】
【数２０】
2015203946.tif 000022
【００７３】
【数２１】
2015203946.tif 000023
【００７４】
【数２２】
2015203946.tif 000024
【００７５】
［実施例１］
図７は、本発明に係る距離および前処理アプローチを用いた、図２に示す３０の画像の平均を示している。すなわち、図７は、アルゴリズム１を用いた最適輸送距離重心（つまり、ワッサーステイン重心）を表している。
【００７６】
［実施形態２：各クラスタの重みに制約のある、重み付けされたポイント群のｋ平均法の演算］
この方法の出発点は、以下の考察から得られる。それぞれがΩに存在する点{x1, x2, … ,xn}の１つの集合において、ｋ平均で最小化されている目的は、重み付けされた２つポイント群の間のワッサースタイン距離の最小化によって、等価的に書き換えることができる。
【００７７】
【数２３】
2015203946.tif 000025
【００７８】
上記式において、ベクトルbはΣｎにおける均一ベクトルb=[1/n,…,1/n]である。上記式の右辺の式は、依然として不均一な重みに対して有効であり、我々は、ポイント{x1, x2, … ,xn}が、ｎ個の変数のシンプレックスΣｎにおける任意のベクトルbによって重み付けされている、さらに一般的な事例を考慮する。
【００７９】
ｋ平均目的が位置(y1,…,yk)とそれらの位置の重みａとの両方を最適化することを示唆していることを、この再定式化は示している。Lloydの独自のアルゴリズムは、ａの値の制約を何ら考慮していないので、我々のアプローチよりも実装は簡単である。しかし、我々は、ａのみの最適化の上記問題の制限、すなわち、
【００８０】
【数２４】
2015203946.tif 000026
【００８１】
の演算が、aがシンプレックスにある場合だけでなく、次の場合においても可能であることを、アルゴリズム１において示している。すなわち、我々が有効な投影PΘを有するためのΣｋのどの凸状の部分集合Θにも存在するという制約をaが受けている場合であっても、アルゴリズム１において詳細に説明された近似的な（劣）勾配降下法を用いることにより、上記演算は可能である。
【００８２】
【数２５】
2015203946.tif 000027
【００８３】
ここで、ａ, ｂ, Ｘが固定されていると仮定して、Ｙのみの関数としての式を最小化したいと考える。Ｙの現在の推定値が利用できるとする。最適条件において、最適輸送変数Ｔ★を計算してワッサースタイン距離Ｗｐを用意することは可能であるので、上記式を以下の式に置き換えることができる。
【００８４】
【数２６】
2015203946.tif 000028
【００８５】
ここで、Ｔ★は最適輸送であるため、上記式の後半のみ有効である。
【００８６】
我々が本項の残りの部分で行うように、Ωがユークリッド空間Rdであると仮定すると、各ポイントyiに関する目的fの（勾配）導関数が、以下の距離関数の知識を活用して簡単な方法で計算できることを達成するための式に多変数微積分を適用することができる。
【００８７】
【数２７】
2015203946.tif 000029
【００８８】
そして、ｄ行（次元）ｋ列（それぞれが上記式で得られる）の勾配行列∇を形成する。当該行列は、各ポイントyiのために、これらの個々の勾配の全てを組み合わせることができる。
【００８９】
要約すると、２つのポイント群に関連する最適輸送Ｔ★、位置Ｘおよび重みｂの一方、位置Ｙおよび重みａの他方、および、各ポイントyiに関する距離dist(xi,yi)の各々の勾配情報が既知であれば、十分に小さいステップεおよび勾配∇, Ｙ←Ｙ- ε∇によって、行列Ｙを更新することができるので、f(Y - ε∇) < f(Y)であることを保証することができる。
【００９０】
以下に示すアルゴリズムでは、簡略化するために、２点間の距離がユークリッド距離であり、ｐ＝２、つまり、任意の２つのベクトルｕとｖが
【００９１】
【数２８】
2015203946.tif 000030
【００９２】
であると仮定している。
【００９３】
この選択は、より単純な式や単純なアルゴリズム記述に変換される。具体的には、この距離を選択する場合、Ｙに近似する第１オーダー近似値ｆの最小値を、閉形式において、得ることができる。ここで、最適輸送Ｔ★は全てのＹについて同一であると仮定する。この場合、このｆの近似値がマトリックスＹの二次関数であり、その解が、Ｘ Ｔ★T diag(b-1)であることが、基本的な微積分法によって示される。ここで、diag(b-1)は、サイズｎの対角行列であり、サイズｎの対角係数は、ベクトルb-1、すなわち、その値がｂの各値の逆数であるベクトルによって構成される。
【００９４】
［アルゴリズム３：Σｋの部分集合Θに制約された重み付けを有するRdにおける経験的計測のワッサーステイン重心］
１．ｎ個の変数のシンプレックスΣnの重みベクトルｂを有する、Rdにおけるｎ個の重み付けされたポイント群{x1, x2, …, xn}を合計する。前記ポイントは、ｄ行ｎ列の行列Ｘとして表すことができる。 
２．Σｋの関連する部分集合Θを該部分集合上の投影PΘとともに定義する。該投影は、所与の任意のベクトルａを、前記ポイントのΘ内の最も近いポイントに戻す関数である。  ３．ＹをＤ行ｋ列の行列に初期化する。各列は、Ｘの列の中から無作為に抽出される。ａをベクトル[1/k, 1/k, ..., 1/k]に設定する。
４．所望に収束するまでａ〜ｄを反復する。
ａ．距離行列
【００９５】
【数２９】
2015203946.tif 000031
【００９６】
を形成する。
ｂ．アルゴリズム１を用いて最適な重みａを計算する工程。Ｍを距離行列パラメータとして、ｂを入力ヒストグラム（Ｎ＝１）として用いる。
ｃ．アルゴリズム２を用いて近似的な最適輸送Ｔ★λを計算する。
ｄ．勾配ステップ：Ｙ← Ｘ Ｔ★T diag(b-1)を更新する。
【００９７】
［実施例２］
本発明のアルゴリズムは、属性ベクトル上の明示的な配分の制約を考慮することを提案する。そして、必要であれば、所望の平滑性を有することができる属性を適用することを提案する。例えば、各クラスタの重心の質量が等しいことが必要な場合、本発明の方法は、図８に示す、除去できない誤差の総和を最小にするような米国国勢調査データのクラスタリングを得ることができる。当然、各重心は、その地図に示された米国の総所得の均一な配分を捕捉していると考えられる。この例では、Ω＝R2であり、距離関数はユークリッド距離であり、ｐ＝２である。本発明のクラスタリング方法（すなわち、アルゴリズム３）の結果は、円を用いて示されており、これらの円は、各重心の均一な重みを表している。
【産業上の利用可能性】
【００９８】
本発明の方法に実用的に関連する適用例を、以下説明する。
【００９９】
［ヒストグラムデータのクラスタリングへの適用（アルゴリズム１）］
ワッサースタイン重心の計算（アルゴリズム１）は、次の場合にヒストグラムのデータセットを要約するために用いることができる。すなわち、当該ヒストグラムに記載されている特徴間の距離が所与である場合（アルゴリズム１のプレゼンテーションの行列Ｍが所与である場合）である。例えば、次の適用例が考えられる。
【０１００】
１．各画像が（例えばＳＩＦＴ特徴を使用して得られる）特徴のヒストグラムとして表現されている、画像のデータベースを考える。これらの特徴間の距離が所与である場合、それらのヒストグラムのワッサースタイン距離に基づいて平均ヒストグラムを計算することができる。
【０１０１】
２．テキストのデータベースが所与である場合、各テキストのＢｏＷ表現（すなわち、単語表現のヒストグラム）、および、辞書内のすべての単語間の適当な距離行列を用いて、ワッサースタイン距離に基づき独自の平均ヒストグラムを生成することができる。このヒストグラムでは、すべてのテキストに共通のキーワードが強調される。そして、通常使用されている単純な等差中項よりも疎性を持ち有益であることが期待できる。
【０１０２】
３．オーディオ録音のデータベースが所与である場合、そのデータセットを備えた辞書抽出ツールを使用することにより、特徴のセットが得られる。当該セットは、各オーディオ録音をそのような特徴のヒストグラムとして特徴付けて表現することができる。我々のアルゴリズムは、最適輸送の意味での平均録音を得るため使用することができ、及び/又は、ｋ平均法の実行時における中央揃え工程として使用することができる。
【０１０３】
４．離散化された集合体（不規則な高さを有するマップ上の人口密度、脳内の活性化、グラフのノードでの動き）における濃度パターンのデータベースを考える。データベースにおける各ポイントは、その集合体上の位置と同数のビンを有するヒストグラムと解釈することができる。集合体は、必ずしもユークリッドではない。すなわち、２つのノード間の距離は、必ずしもユークリッド距離として計算されるのではなく、その代わりに、局所的な制約を考慮する測地距離として計算される。この測地距離は、例えば、不規則な起伏のあるマップ上の２点間の高低差を考慮した徒歩距離、および、変化する端質量を有する不完全な連続グラフのノードであり、いずれの場合でも、全ての組み合わせの最短経路アルゴリズムの結果として、距離が計算される。アルゴリズム１を用いることにより、グラフ中の全てのノード対の間の距離を表す距離行列Ｍに加えて提供されるデータセットを用いることにより、ワッサースタインの意味での平均濃度パターンを計算することができる。
【０１０４】
［制限されたクラスタリング、容量制限を伴うアンテナ中継配置への適用（アルゴリズム２）］
各点ｘにおける劣勾配を認める距離関数Ｄ（ｘ、ｙ）を用いて、ユークリッド空間に不規則に分布した集合が所与であると仮定する。実際には、その集団に適用される標準的なｋ平均のコスト関数は、重心Ｘの組と、ａのエントロピーが非常に小さくなるような重みベクトルａとを用いて最小化することができる。非専門用語では、データセットにおける大半の元アイテムは、重心の小さな部分集合≪ kに帰することができる。これに対し、他の全ての重心は、集合全体の非常にわずかな量を捕捉する。
【０１０５】
これは、より規則正しい属性が要求されるｋ平均の適用において望ましくない可能性がある。例えば、センサの配置において、各重心（センサ）が機能することができるデータポイント（利用者）の数が限られている場合、我々は、属性がそれらの制限を受け入れることを保証したいと考える。従来のｋ平均はそのような制限を考慮することができない。これに対し我々は、アルゴリズム２を用いて、log(k)-α（α≧0はエントロピーの閾値を規定する）の閾値よりも大きいエントロピーを有する、ΣｋにおけるヒストグラムのセットΘのために設定することにより、それらを保証することができる。その場合、カルバック・ライブラー情報量を用いて設定された非負ベクトルの投影は、べき指数ｔを求めるように、簡単に実行することができる。ここで、構成要素的に、
【０１０６】
【数３０】
2015203946.tif 000032
【０１０７】
と定義されるΣｋにおけるエントロピー
【０１０８】
【数３１】
2015203946.tif 000033
【０１０９】
は、log(k)-αに等しい。α=0の場合、セットΘは、重み1/kの均一なヒストグラムに減少する。
【図１】
2015203946.tif 000039
【図２】
2015203946.tif 000040
【図３】
2015203946.tif 000041
【図４】
2015203946.tif 000042
【図５】
2015203946.tif 000043
【図６】
2015203946.tif 000044
【図７】
2015203946.tif 000045
【図８】
2015203946.tif 000046
【外国語明細書】
2015203946000001.pdf
2015203946000002.pdf
2015203946000003.pdf
2015203946000004.pdf