JP6662754B2 - Ｌ１グラフ計算装置、ｌ１グラフ計算方法及びｌ１グラフ計算プログラム - Google Patents

Description

本発明は、Ｌ_１グラフ計算装置、Ｌ_１グラフ計算方法及びＬ_１グラフ計算プログラムに関する。

グラフは基盤的なデータ構造であり、多次元データを解析する上で重要なデータ構造である。グラフ構造としてはｋ−近傍グラフが有名であるが、データのノイズに弱いという問題がある。この問題を解決するために、lassoにおけるＬ_１最適化を解くことで計算できるＬ_１グラフがMeinshausenらによって提案された（非特許文献１参照）。具体的に、非特許文献１記載の手法は、グラフの各ノードにlassoを適用し、lassoによって得られた疎な回帰結果をエッジの重みとするものである。このＬ_１グラフは、ｋ−近傍グラフよりノイズに強いという特徴がある。

N. Meinshausen and P. Buhlmann, "High Dimensional Graphs and Variable Selection with the Lasso", The Annals of Statistics, 34(3):1436−1462, June 2006.

しかしながら、従来の手法では、lassoを用いてＬ_１グラフを計算するには高い計算コストが必要であり、Ｌ_１グラフを高速に計算することが難しいという問題があった。これは、lassoでは座標降下法を用いて回帰を行うが、座標降下法は、繰り返し収束するまで重みを計算する必要があるためである。

本発明は、上記に鑑みてなされたものであって、lassoを用いてＬ_１グラフを高速に計算することができるＬ_１グラフ計算装置、Ｌ_１グラフ計算方法及びＬ_１グラフ計算プログラムを提供することを目的とする。

上述した課題を解決し、目的を達成するために、本発明に係るＬ_１グラフ計算装置は、lassoによるＬ_１グラフを計算するＬ_１グラフ計算装置であって、入力される多次元行列の特異値分解を計算するＳＶＤ計算部と、エッジの重みを計算する対象のノードを選択するノード選択部と、初めに重みを更新するエッジの集合を設定する第１エッジ集合設定部と、初めに重みを更新するエッジ集合に含まれるエッジの重みを収束するまで更新する第１重み更新部と、グラフの計算におけるパラメータの初期化を行うパラメータ初期化部と、重みを更新するエッジ集合として、特異値分解の計算結果を用いて、非零の重みを持ち得ないエッジを除外したエッジ集合を設定する第２エッジ集合設定部と、重みを更新するエッジ集合を、第２エッジ集合設定部が設定したエッジ集合に追加する追加エッジ計算部と、追加エッジ計算部が追加したエッジ集合から、非零の重みを有するエッジを、一つ一つ、重みを更新するエッジ集合に追加するエッジ追加部と、エッジ追加部がエッジを追加したエッジ集合に含まれるエッジの重みを、収束するまで更新する第２重み更新部と、選択されていないノードに対してエッジの重みを設定する重み設定部と、を有することを特徴する。

本発明によれば、lassoを用いてＬ_１グラフを高速に計算することができる。

図１は、本実施の形態１に係るＬ_１グラフ計算装置の要部構成の一例を示すブロック図である。 図２は、定義１を示す図である。 図３は、定義２を示す図である。 図４は、定義３を示す図である。 図５は、補題１を示す図である。 図６は、補題２を示す図である。 図７は、補題３を示す図である。 図８は、補題４を示す図である。 図９は、補題５を示す図である。 図１０は、行列ΧのＳＶＤを用いた近似方法を示す図である。 図１１は、図１０に示す定義４を示す図である。 図１２は、補題６を示す図である。 図１３は、補題７を示す図である。 図１４は、定義５を示す図である。 図１５は、補題８を示す図である。 図１６は、実施の形態１に係るＬ_１グラフ計算処理のアルゴリズムを示す図である。 図１７は、定理１を示す図である。 図１８は、定理２を示す図である。 図１９は、実施の形態１に係るＬ_１グラフ計算処理の流れを示すフローチャートである。 図２０は、プログラムが実行されることにより、Ｌ_１グラフ計算装置が実現されるコンピュータの一例を示す図である。

以下、図面を参照して、本発明の一実施の形態を詳細に説明する。なお、この実施の形態により本発明が限定されるものではない。また、図面の記載において、同一部分には同一の符号を付して示している。

［従来のlassoを用いたＬ_１グラフ計算手法］
まず、従来のlassoを用いたＬ_１グラフ計算手法について説明する。lassoによるＬ_１グラフでは、ノードが多次元データにおけるデータポイントに対応し、エッジが回帰における関係に対応する。ここで、ｐをグラフにおけるノードとし、Ｖをグラフにおけるノードの集合とすると、ｐ∈Ｖであるようなノードを一つ一つ選択し、lassoを計算することによって、エッジの重みを求める。lassoは、ほとんどのノードに零の重みを与えるため、グラフは疎な構造となる。なお、各集合を示すシンボルは、明細書中ではアルファベットの大文字で示し、図面中ではアルファベットの大文字を中抜きした字体で示す。

Χ∈Ｒ^Ｎ×Ｍを、Ｎ個のデータポイントと、Ｍ次元からなるデータとし、ｘ_ｐ＝（χ_ｐ［１］，χ_ｐ［２］，・・・，χ_ｐ［Ｍ］）を、行列Χにおけるｐ番目の行ベクトルとすると、行ベクトルｘ_ｐは、ｐ番目のデータポイント及びノードｐに対応する。本実施の形態１ではベクトルｘ_ｐは、平均が０、分散が１に正規化されているものとする。

そして、ｗ_ｐを１×Ｎの重みベクトルとし、そのｕ番目の要素ｔω_ｐ［ｕ］は、ノードｕからｐのエッジの重みとする。lassoによるＬ_１グラフ計算では、ω_ｐ［ｐ］=０として、以下の式（１）に示す目的関数を最小化するように重みを計算する。

式（１）において最初の項は回帰のよる２乗誤差に対応し、２番目の項は重みのＬ_１制約に対応する。式（１）はチューニングパラメータλが大きくなるとグラフがより疎な構造を持つことを示している。もし行列Χがフル行ランクであれば式（１）は一つの解を持ち、そうでなければ解は一つにならない。

lassoを高速に解く方法として座標降下法がある。座標降下法は、収束するまでエッジの重みを一つ一つ更新するものである。具体的に、座標降下法では、以下の式（２）により重みの更新を行う。

ここでＳ［・，・］は、以下の式（３）を用いて計算する。

さらに、ｚ_ｐ［ｕ｜ｗ_ｐ］は、ノードｕのベクトルｗ_ｐのパラメータであり、以下の式（４）を用いて計算する。

上記の式にみられるように、各パラメータは、以下の性質を有する。

結果としてＴを収束までの繰り返し計算回数としたときに、式（２）によりノードｐへの重みを計算するにはＯ（ｄ_ｐＭＴ）の計算コストが必要になる。したがって、lassoによるＬ_１グラフを求めるためには、高い計算コストが必要になる。

［実施の形態１の概要］
本実施の形態１において用いられる手法によれば、より高速にlassoによるＬ_１グラフを計算することができる。まず、本実施の形態１の計算手法の概要を述べた後に、本実施の形態１による計算手法の詳細を述べる。

本実施の形態１では、重みを更新するエッジの集合を求めるために非零の重みを持ち得ない不要なエッジを枝狩りして、計算を高速化する。言い換えると、本実施の形態１では、必ず非零の重みをもつエッジに対してのみ、収束するまで重みの更新を行い、それ以外のエッジについては重みを更新する計算処理を実行しない。

続いて、本実施の形態１では、非零の重みを持つ可能性のあるエッジに対して更新を行う。この場合、本実施の形態１では、非零の重みを持つエッジを特定するために、ＫＫＴ（Karush-Kuhn-Tucker）条件のスコアの上限値と下限値とを計算する。ＫＫＴ条件は、誤って枝狩りされたエッジを見つけるために提案されたものである（詳細は、R. Tibshirani, J. Bien, J. Friedman, T. Hastie, N. Simon, J. Taylor and R. J. Tibshirani, “Strong Rules for Discarding Predictors in Lasso-type Problems”, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(2):245−266, 2012.参照）。

さらに、本実施の形態１では、高速に更新を行うために２つの計算方法を利用する。初めの計算方法は、全ての重みに対して共通であるlassoの残差を用いたものである（後述の式（８），式（９）参照）。この方法では、各重みの更新に異なる回帰結果を用いなくてよいため高速な計算が可能になる。本実施の形態１では、この計算方法を重みが零である場合の更新に用いる。

一方、本実施の形態１では、もし重みが非零であれば、既存研究で提案されている内積を用いる計算方法によって更新を行う（後述の式（１０）参照。詳細は、J. Friedman, T. Hastie, and R. Tibshirani, “Regularization Paths for Generalized Linear Models via Coordinate Descent”, Journal of Statistical Software, 33(1):1−22, 2 2010.参照）。

このように、本実施の形態１では、重み更新についての繰り返し計算を行う前に、非零の重みを持ち得ないエッジを枝狩りする。言い換えると、本実施の形態１では、非零の重みを持ち得ないエッジについては、重み更新についての繰り返し計算を行わない。

さらに、本実施の形態１では、各重み更新についての繰り返し計算において重みが零になるエッジを枝狩りする。言い換えると、本実施の形態１では、エッジが非零の重みを持たなければ、このエッジに対し、重み更新についての繰り返し計算を行わない。

したがって、本実施の形態１では、従来の手法よりも繰り返し計算を少ない回数とすることができるため、lassoを用いても高速にＬ_１グラフを計算することができる。

［実施の形態１の構成］
まず、図１を用いて、実施の形態１に係るＬ_１グラフ計算装置の構成について説明する。図１は、実施の形態１に係るＬ_１グラフ計算装置の要部構成の一例を示すブロック図である。図１に示すように、Ｌ_１グラフ計算装置１は、ＳＶＤ（Singular Value Decomposition：特異値分解）計算部１１、ノード選択部１２、第１エッジ集合設定部１３、第１重み更新部１４、パラメータ初期化部１５、第２エッジ集合設定部１６、追加エッジ計算部１７、エッジ追加部１８、第２重み更新部１９及び重み設定部２０を有する。

ＳＶＤ計算部１１は、入力される多次元行列の特異値分解を計算する。ノード選択部１２は、エッジの重みを計算する対象のノードを選択する。そして、ノード選択部１２の後段の各部は、この計算対象のノードに対して重み計算を繰り返し行う。第１エッジ集合設定部１３は、初めに重みを更新するエッジの集合を設定する。第１重み更新部１４は、第１エッジ集合設定部１３によって設定されたエッジ集合に含まれるエッジの重みを収束するまで更新する。パラメータ初期化部１５は、グラフの計算におけるパラメータの初期化を行う。

そして、第２エッジ集合設定部１６は、重みを更新するエッジの集合として、特異値分解の計算結果を用いて、非零の重みを持ち得ないエッジを除外したエッジ集合を設定する。追加エッジ計算部１７は、重みを更新するエッジ集合を、第２エッジ集合設定部１６が設定したエッジ集合に追加する。エッジ追加部１８は、追加エッジ計算部１７が追加したエッジ集合から、非零の重みを有するエッジを、一つ一つ、重みを更新するエッジ集合に追加する。第２重み更新部１９は、エッジ追加部１８がエッジを追加したエッジ集合に含まれるエッジの重みを、収束するまで更新する。重み設定部２０は、ノード選択部１２によって選択されていないノードに対してエッジの重みを設定する。これらの第２エッジ集合設定部１６、追加エッジ計算部１７、エッジ追加部１８、第２重み更新部１９及び重み設定部２０は、エッジの重み更新についての繰り返し計算を行う。

［エッジ集合設定処理］
ここで、実施の形態１に係るＬ_１グラフ計算方法のうち、エッジの重み更新についての繰り返し計算のうちのエッジ集合設定処理について説明する。まず、重みの更新のための計算を前に、第２エッジ集合設定部１６が、２つのエッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］とを計算する。エッジ集合Ｍ［ｗ_ｐ］は、重みベクトルｗ_ｐを用いて更新を行ったときに、重みが非零に必ずなるエッジの集合である。また、エッジ集合Ｃ［ｗ_ｐ］は、重みベクトルｗ_ｐを用いて更新したときに重みが非零になる可能性のあるエッジの集合である。

理論的に、エッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］とはＫＫＴ条件に基づいている。したがって、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］とを、ＫＫＴ条件のスコアの上限値と下限値とを用いて求める。第２エッジ集合設定部１６は、ＳＶＤ計算結果を用いて求めたＫＫＴ条件のスコアの上限値と下限値とを用いて、非零の重みを持ち得ない不要なエッジを特定し、このエッジを枝狩りする。

そこで、エッジ集合の定義を述べる前に、ＫＫＴ条件について述べる。ＫＫＴ条件は、図２に示す定義１のように定義される。図２は、定義１を示す図である。

ｕ∈Ｖ／ｐであるようなノードｕに対して、もし、ノードｕからｐのエッジの重みが重みベクトルｗ_ｐで更新されたとき、Ｋ［ｕ｜ｗ_ｐ］は、更新後の、ｗ_ｐ［ｕ］の値により、以下の式（７）のように計算できる。

直接ＫＫＴ条件のスコアを計算すれば重みが非零になるか否かが分かる。ただし、ＫＫＴ条件のスコアを計算するためには、最長でＯ（ＮＭ）の計算コストがかかる。これは、式（Ａ）における行ベクトルｘ_ｐ、重みベクトルｗ_ｐ、行列Χ、ｘ_ｕの転置ｘ_ｕ ^Ｔの大きさが、１×Ｍ、１×Ｎ、Ｎ×Ｍ、Ｍ×１であるためである。そこで、本実施の形態１では、第２エッジ集合設定部１６が、ＫＫＴ条件のスコアの上限値と下限値との計算を高速化できるように定義づけを行っている。なお、ＫＫＴ条件のスコアの上限値と下限値との計算方法は後述する。

続いて、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗ_ｐ］について、図３に示す定義２のように計算する。図３は、定義２を示す図である。第２エッジ集合設定部１６は、エッジ集合Ｃ［ｗ_ｐ］について、図４に示す定義３のように計算する。図４は、定義３を示す図である。そして、エッジ集合については、補題１，２が成り立つ。図５は、補題１を示す図である。また、図６は、補題２を示す図である。

エッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］との関係については、図７に示す補題３が成り立つ。図７は、補題３を示す図である。このように、第２エッジ集合設定部１６は、定義２，３を満たすエッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］を求める。言い換えると、第２エッジ集合設定部１６は、非零の重みを持ち得ない不要なエッジを除外したエッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］を求める。

［重み更新処理］
続いて、エッジの重み更新についての繰り返し計算のうちの重み更新処理における計算方法について述べる。Ｌ_１グラフ計算装置１では、二つの計算方法を用いる。まず、初めの計算方法（第１の計算方法）について説明する。初めの計算方法として、ノードｕのベクトルｗ_ｐのパラメータであるパラメータｚ_ｐ［ｕ｜ｗ_ｐ］を計算し、エッジごとに異なる解析結果を使わずともよい方法を用いる。具体的には、以下の式（８）を用いて計算を行う。

ここでｒ_ｐ［ｉ］は、要素χ_ｐ［ｉ］のlassoの回帰結果における残差であり、以下の式（９）を用いて計算される。

式（９）で使用される残差は、各重みに対して同じである。このため、もし回帰結果が更新後に変わることがなければ、パラメータｚ_ｐ［ｕ｜ｗ_ｐ］を高速に計算することができる。パラメータｚ_ｐ［ｕ｜ｗ_ｐ］を計算するために式（８）において必要になる計算コストは、Ｏ（Ｍ）である。この式（８）に対し、以下の補題４に示す性質が成り立つ。図８は、補題４を示す図である。

エッジ追加部１８は、重み更新の前に、もしエッジが零の重みを持てば、式（８）を用いて重みベクトルに関するパラメータｚ_ｐ［ｕ｜ｗ_ｐ］を計算する。この零の重みを持つエッジは、更新後も、再度、零の重みを持つことが期待されるため、残差は、更新の前後で同じ値を持つことが期待できる。このため、次の繰り返し計算においても、効果的に残差ｒ_ｐ［ｉ］を利用することができる。言い換えると、この零の重みを持つエッジは、更新後も、再度、零の重みを持つため、次の重み更新についても繰り返し計算の対象から除外することができる。

一方、エッジが非零の重みを更新の後で持つ場合であっても、第２重み更新部１９は、式（９）を用いることで、逐次的に、残差ｒ_ｐ［ｉ］をＯ（１）の計算コストで更新することができる。

そして、第２重み更新部１９は、もしエッジが更新の前に非零の重みを持つ場合には、既存研究で提案されている以下の内積を用いる式（１０）を用いて、重みを更新する。

次に、重み更新についての繰り返し計算において、更新対象となるエッジ集合について説明する。前述したように、Ｌ_１グラフ計算装置１では、重み更新についての繰り返し計算の前に、エッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］とを計算する。Ｕを、更新を行うエッジの集合としたときに、素朴にはＵ＝Ｍ［ｗ_ｐ］またはＵ＝Ｃ［ｗ_ｐ］とすることにより更新を行うエッジを減らすことができる。

ただし、Ｌ_１グラフ計算装置１では、ＫＫＴ条件の上限値と下限値とを用いてエッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］とを求めているため、エッジ集合には重みが零になるエッジが含まれる。そして、Ｌ_１グラフ計算装置１では、重みが零になるエッジに対して繰り返し計算が行われることがあり得る。

そこで、高速に重みの更新を行うために、Ｌ_１グラフ計算装置１では、エッジ集合Ａからエッジ集合Ｕへ一つ一つエッジを足す処理を行う。Ｌ_１グラフ計算装置１では、まず、第１重み更新部１４が、エッジ集合Ｕを非零の重みを持つエッジから求め、そのエッジ集合Ｕに対して収束するまで重みの更新を行う。収束後、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗ_ｐ］を計算し、エッジを加える集合ＡをＡ＝Ｍ［ｗ_ｐ］＼Ｕとして計算する。

このエッジ集合Ａのエッジに対して、追加エッジ計算部１７、エッジ追加部１８及び第２重み更新部１９が、一つ一つ重みの更新を行う。そして、もし、あるエッジが非零の重みを持つ場合には、エッジ追加部１８が、そのエッジをエッジ集合Ｕへ追加し、第２重み更新部１９が、集合Ｕについて、エッジの重みが収束するまで重み更新について繰り返し計算を行う。

一方、あるエッジが非零の重みを持たない場合、エッジ追加部１８は、このエッジに対して、繰り返し計算をすることなく枝狩りを行う。同様に、第２エッジ集合設定部１６は、Ａ＝Ｃ［ｗ_ｐ］＼Ｕを計算し、追加エッジ計算部１７、エッジ追加部１８及び第２重み更新部１９が、このエッジ集合Ａのエッジに対して重みの更新を行う。

Ｌ_１グラフ計算装置１は、もしエッジが非零の重みを持たなければ、そのエッジの重み更新のための繰り返し計算を行わない。このため、Ｌ_１グラフ計算装置１では、全てのエッジについてエッジの重み更新のための繰り返し計算を行っていた従来装置と比較し、高速な計算処理が可能になる。理論的には、この手法は、以下の補題５に示す性質に基づいている。図９は、補題５を示す図である。

補題５にあるように、ノードｕからｐへのエッジの重みの更新処理において、その重み更新についての繰り返し計算において、零であり座標降下法により収束する場合には、この結果は、座標降下法による解析結果に影響しない。したがって、重み更新についての繰り返し計算において零であるとともに座標降下法により収束するエッジについては、このエッジを枝狩りできることを示している。上記のように、Ｌ_１グラフ計算装置１は、エッジが非零の重みを持たない場合には、そのエッジの繰り返し計算を行わないため、全てのエッジについてエッジの重み更新のための繰り返し計算を行っていた従来装置と比較して高速な処理が可能になる。

［ＫＫＴスコアの上限値と下限値とを求めるための計算処理］
次に、ＫＫＴスコアの上限値と下限値とを求めるための計算処理について説明する。上限値と下限値を計算する方法として、Ｌ_１グラフ計算装置１は、ＳＶＤを用いる方法と、逐次的に更新する方法とを用いる。

まず、ＳＶＤを用いる方法について説明する。この方法では、行列Χを、ＳＶＤを用いて、図１０に示すように近似する。図１０は、行列ΧのＳＶＤを用いた近似方法を示す図である。そして、図１１は、図１０に示す定義４を示す図である。そして、ＫＫＴ条件のスコアの上限値と下限値とに対し、以下の補題６に示す性質が成り立つ。図１２は、補題６を示す図である。

次に、逐次的に上限値と下限値とを計算する方法について説明する。前述したように、Ｌ_１グラフ計算装置１では、第２エッジ集合設定部１６、追加エッジ計算部１７、エッジ追加部１８及び第２重み更新部１９が、エッジ集合Ａからエッジを一つ一つエッジ集合Ｕに足しながらエッジの重みを更新する。この時、ＳＶＤによる上限値と下限値とを用いることにより高速にエッジ集合Ｍ［ｗ_ｐ］とＣ［ｗ_ｐ］を計算することができる。ここで、さらに、高速にエッジ集合を計算するために、第２エッジ集合設定部１６は、逐次的にＫＫＴ条件のスコアの上限値と下限値とを計算する方法を用いる。

逐次的に上限値と下限値を更新する方法を述べる前に、まず、ＫＫＴ条件とパラメータｚ_ｐ［ｕ｜ｗ_ｐ］との関係について説明する。図１３は、補題７を示す図である。図１３の補題７は、ＫＫＴ条件とパラメータｚ_ｐ［ｕ｜ｗ_ｐ］との関係を示すものである。

補題７は、もし（１）あるエッジがエッジ集合に含まれ、（２）そのエッジのパラメータｚ_ｐ［ｕ｜ｗ_ｐ］が既に計算済みであれば、ＫＫＴ条件のスコアはＯ（１）の計算コストで求めることができることを示している。第２エッジ集合設定部１６は、エッジを追加した後、上限値と下限値とを、以下の定義５に示すように逐次的に計算する。図１４は、定義５を示す図である。

そして、ＫＫＴ条件のスコアの上限値及び下限値については、補題８に示す性質が成り立つ。図１５は、補題８を示す図である。

前述したように、Ｌ_１グラフ計算装置１は、高速に重みを更新するために追加されたエッジが回帰結果を変えることがあるたびに回帰の残差を計算する。座標降下法は、重みを一つずつ更新するため、式（Ｈ）（図１４参照）から、δ_ｐ ^ｉとΔ_ｐをＯ（１）の計算コストで更新することができる。Ｌ_１グラフ計算装置１は、結果として定義５を用いることによって、ＫＫＴ条件のスコアの上限値と下限値とを逐次的にＯ（１）の計算コストで更新することができる。

なお、Ｌ_１グラフ計算装置１において、ＫＫＴ条件のスコアの上限値と下限値を計算するという観点では、定義５で与えられる逐次的な方法は、定義４で与えられるＳＶＤを用いる方法と同じである。ただし、これらは、それぞれの計算方法がまったく異なるため、異なる上限値と下限値となる。さらに、逐次的な方法は、ＳＶＤによる方法より高速に上限値と下限値を計算することができる。このため、Ｌ_１グラフ計算装置１は、まず逐次的な方法でＫＫＴ条件のスコアの上限値と下限値を計算し、それでも枝狩りできない場合にＳＶＤによる方法を用いて上限値と下限値を計算する。

［アルゴリズム］
図１６を参照して、Ｌ_１グラフ計算装置１の各部の処理について説明する。図１６は、実施の形態１に係るＬ_１グラフ計算処理のアルゴリズムを示す図である。図１６に示すAlgorithm１では、行列Χ、チューニングパラメータλ、ＳＶＤのランクｍを入力とし、行列Ｗを出力とする。Ｗは、そのｐ番目の行が重みベクトルｗ_ｐに対応するＮ×Ｎの行列である。Ｐは、重みを計算するために選択されたノードの集合である。実施の形態１では、Ｌ_１グラフにおいてノードｕ ｔｏ ｖの重みはノードｖ to ｕの重みと似ているという性質を用いて重みを初期化する。

ＳＶＤ計算部１１は、まずノード集合Ｐを初期化し、上限値と下限値を計算するために用いる行列ΧのＳＶＤを計算する（Algorithm１の１，２行目）。そして、ノード選択部１２は、重みベクトルに対して最大のＬ_１を持つノードを選択する（Algorithm１の４行目）。続いて、第１エッジ集合設定部１３は、初めに重みを更新するエッジ集合Ｕを設定する（Algorithm１の５行目）。そして、第１重み更新部１４は、式（１０）を用いて非零の重みを持つエッジに対して更新を行い（Algorithm１の６行目）、パラメータ初期化部１５は、ＫＫＴ条件のスコアの上限値及び下限値を初期化する（Algorithm１の８行目）。

ここで、本実施の形態１では、２つのステップによって選択されたノードの重みを計算する。第２エッジ集合設定部１６、追加エッジ計算部１７、エッジ追加部１８及び第２重み更新部１９は、ノードの重みの更新についての計算における初めのステップにおいて、エッジ集合Ｍ［ｗ_ｐ］を用いる。続いて、第２エッジ集合設定部１６、追加エッジ計算部１７、エッジ追加部１８及び第２重み更新部１９は、ノードの重みの計算における次のステップにおいて、エッジ集合Ｃ［ｗ_ｐ］を用いる。

第２エッジ集合設定部１６は、不要なエッジを枝刈りして更新処理回数を低減するために、逐次的な方法である定義５の方法を用いて、ＫＫＴ条件のスコアの上限値と下限値を計算する（Algorithm１の１１−１２行目）。そして、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗ_ｐ］またはエッジ集合Ｃ［ｗ_ｐ］との条件を確認する（Algorithm１の１６行目及び２４行目）。

続いて、エッジがエッジ集合の条件（定義２及び定義３）を満たす可能性がある場合には、第２エッジ集合設定部１６が、ＳＶＤによって、すなわち、定義４の方法を用いて、ノードの上限値と下限値とを計算し（Algorithm１の１７−１９行目と２５−２７行目）、エッジ集合を決定する（Algorithm１の２０行目と２８行目）。

さらに、第２エッジ集合設定部１６が、エッジ集合Ｕに一つ一つエッジを追加する。そのため、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］を計算した後に、エッジを加えるエッジ集合ＡをＡ＝Ｍ［ｗ_ｐ］＼ＵまたはＡ＝Ｃ［ｗ_ｐ］＼Ｕとして計算する（Algorithm１の２０行目と２８行目）。もしエッジの重みが零であれば、補題５に示すように回帰結果に影響がなく、第２重み更新部１９は、式（８）を用いて、その重みを高速に計算することができる（Algorithm１の３０行目）。もし追加するエッジがなければ、重みは明らかに収束する。そのため、Ｌ_１グラフ計算装置１では、エッジ集合Ａが、Ａ＝０であれば、グラフを計算するステップを進める（Algorithm１の３４−３５行目）。

そして、重みの更新についての繰り返し計算の後、重み設定部２０は、選択されなかったノードの重みを設定する（Algorithm１の３７−３８行目)。

なお、Algorithm１では、一つのスレッドで処理することを前提にしていたが、本実施の形態では、複数のスレッドで処理することも可能である。複数のスレッドで処理するため、ノードを複数のグループに分割し、並列に各グループに対してAlgorithm１を実行する。各グループにおいてエッジはすべてのノードから計算するために、並列処理を行ってもひとつのスレッドで処理するのと同じグラフを計算することができる。ノードを複数のグループに分割するためには、k-means法を用いる。k-means法は、一つのスレッドで実行するが、ＳＶＤを用いて次元数を削減しており、また、k-means++を用いて初期のシードを決定するため、高速にノードを分割することができる（詳細は、D. Arthur and S. Vassilvitskii, “k-means++: The Advantages of Careful Seeding”, In SODA, pages 1027−1035, 2007.参照）。

そして、このAlgorithm１に対して、定理１，２に示す性質が成り立つ。図１７は、定理１を示す図である。図１８は、定理２を示す図である。

［実施の形態１の処理］
次に、図１９は、実施の形態１に係るＬ_１グラフ計算処理の流れを示すフローチャートである。

図１９に示すように、まず、行列Χ、チューニングパラメータλ、ＳＶＤのランクｍが入力されると（ステップＳ１）、ＳＶＤ計算部１１は、ノード集合Ｐを初期化し（ステップＳ２）、行列ΧのＳＶＤを計算する（ステップＳ３）。そして、ノード選択部１２は、次元数ｉ、データポイント数Ｎを初期化し（ステップＳ４）、重みベクトルに対して最大のＬ_１を持つノードを選択する（ステップＳ５）。第１エッジ集合設定部１３は、初めに重みを更新するエッジ集合Ｕを設定する（ステップＳ６）。続いて、第１重み更新部１４は、式（１０）を用いて非零の重みを持つエッジに対して更新を行い（ステップＳ７）、パラメータ初期化部１５は、ＫＫＴ条件のスコアの上限値及び下限値を初期化する（ステップＳ８〜ステップＳ１０）。

続いて、第２エッジ集合設定部１６は、ステップを初期化し（ステップＳ１１）、逐次的な方法である定義５の方法を用いて、ＫＫＴ条件のスコアの上限値と下限値を計算する（ステップＳ１２〜ステップＳ１４）。そして、第２エッジ集合設定部１６は、ノードの重みの計算におけるステップが最初のステップであるか否かを判断する（ステップＳ１５）。

第２エッジ集合設定部１６は、ノードの重みの計算におけるステップが最初のステップであると判断した場合（ステップＳ１５：Ｙｅｓ）、エッジ集合Ｍ［ｗ_ｐ］を初期化する（ステップＳ１６）。第２エッジ集合設定部１６は、各ノードｕについて、エッジ集合Ｍ［ｗ_ｐ］の条件を確認し（ステップＳ１７，ステップＳ１８）、エッジ集合Ｍ［ｗ_ｐ］の条件を満たす場合には（ステップＳ１８：Ｙｅｓ）、不要な更新処理を枝狩りするために、逐次的な方法である定義４の方法を用いて、ＫＫＴ条件のスコアの上限値と下限値を計算する（ステップＳ１９）。

第２エッジ集合設定部１６は、計算したＫＫＴ条件のスコアの上限値と下限値とが定義２を満たすか否かを判断する（ステップＳ２０）。第２エッジ集合設定部１６は、計算した上限値と下限値とが定義２を満たすと判断した場合には（ステップＳ２０：Ｙｅｓ）、このエッジをエッジ集合Ｍ［ｗ_ｐ］に追加する（ステップＳ２１）。第２エッジ集合設定部１６は、計算した上限値と下限値とが定義２を満たさないと判断した場合には（ステップＳ２０：Ｎｏ）、このエッジを、エッジ集合Ｍ［ｗ_ｐ］に追加せず、枝刈りする。第２エッジ集合設定部１６は、各ノードｕについて、ステップＳ１８〜ステップＳ２１の処理を行う。

一方、第２エッジ集合設定部１６は、ノードの重みの計算におけるステップが最初のステップでないと判断した場合（ステップＳ１５：Ｎｏ）、エッジ集合Ｃ［ｗ_ｐ］を初期化する（ステップＳ２４）。第２エッジ集合設定部１６は、各ノードｕについて、エッジ集合Ｃ［ｗ_ｐ］の条件を確認し（ステップＳ２５，ステップＳ２６）、エッジ集合Ｃ［ｗ_ｐ］の条件を満たす場合には（ステップＳ２６：Ｙｅｓ）、不要な更新処理を枝狩りするために、逐次的な方法である定義４の方法を用いて、ＫＫＴ条件のスコアの上限値と下限値を計算する（ステップＳ２７）。

第２エッジ集合設定部１６は、計算したＫＫＴ条件のスコアの上限値と下限値とが定義３を満たすか否かを判断する（ステップＳ２８）。第２エッジ集合設定部１６は、計算した上限値と下限値とが定義３を満たすと判断した場合には（ステップＳ２８：Ｙｅｓ）、このエッジをエッジ集合Ｃ［ｗ_ｐ］に追加する（ステップＳ２９）。第２エッジ集合設定部１６は、計算した上限値と下限値とが定義３を満たさないと判断した場合には（ステップＳ２８：Ｎｏ）、このエッジを、エッジ集合Ｃ［ｗ_ｐ］に追加せず、枝刈りする。第２エッジ集合設定部１６は、各ノードｕについて、ステップＳ２６〜ステップＳ２９の処理を行う。

そして、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗ_ｐ］とエッジ集合Ｃ［ｗ_ｐ］を計算した後に、エッジを加えるエッジ集合ＡをＡ＝Ｍ［ｗ_ｐ］＼ＵまたはＡ＝Ｃ［ｗ_ｐ］＼Ｕとして計算する（ステップＳ２３及びステップＳ３１）。

続いて、第２重み更新部１９は、各ノードｕについて、式（８）式から重みを計算し（ステップＳ３２，ステップＳ３３）、計算した重みが非零であるか否かを判断する（ステップＳ３４）。第２重み更新部１９は、計算した重みが非零であると判断した場合には（ステップＳ３４：Ｙｅｓ）、このエッジをエッジ集合Ｕに追加して式（１０）を用いて収束するまで重みを更新する（ステップＳ３５）、一方、第２重み更新部１９は、計算した重みが非零でないと判断した場合には（ステップＳ３４：Ｎｏ）、このエッジをエッジ集合Ｕに追加せず、枝刈りする。

そして、第２重み更新部１９は、エッジ集合がＡ＝０であるか否かを判断する。第２重み更新部１９は、エッジ集合がＡ＝０であれば（ステップＳ３７：Ｙｅｓ）、グラフを計算するステップを進める（ステップＳ３８）。エッジ集合がＡ＝０でないと第２重み更新部１９が判断した場合（ステップＳ３７：Ｎｏ）、或いは、ステップＳ３８終了後、Ｌ_１グラフ計算装置１は、ステップが２以下であるか否かを判断する（ステップＳ３９）。Ｌ_１グラフ計算装置１は、ステップが２以下であると判断した場合（ステップＳ３９：Ｙｅｓ）。ステップＳ１２に戻る。

一方、ステップが２以下でないとＬ_１グラフ計算装置１が判断した場合（ステップＳ３９：Ｎｏ）、選択されなかったノードをノード集合Ｐに追加しする（ステップＳ４０）。そして、重み設定部２０は、選択されなかったノードごとに、該ノードの重みを設定する（ステップＳ４１〜ステップＳ４３）。Ｌ_１グラフ計算装置１は、各ｉ、各Ｎについて、ステップＳ３以降の処理を終了した場合、行列Ｗをを出力する（ステップＳ４５）。

［実施の形態１の効果］
このように、実施の形態１では、重み更新についての繰り返し計算を行う前に非零の重みを持ち得ないエッジを枝狩りしている。言い換えると、実施の形態１では、非零の重みを持ち得ないエッジについては重みを更新するための繰り返し計算を実行しない。さらに、本実施の形態１では、各重み更新についての繰り返し計算において重みが零になるエッジを枝狩りしている。したがって、本実施の形態１では、全てのエッジについてエッジの重み更新のための繰り返し計算を行っていた従来装置と比較し、重み更新についての繰り返し計算を少ない回数とすることができる。この結果、本実施の形態１によれば、lassoを用いて高速にＬ_１グラフを計算することができる。

［実施の形態１のシステム構成について］
図１に示したＬ_１グラフ計算装置１の各構成要素は機能概念的なものであり、必ずしも物理的に図示のように構成されていることを要しない。すなわち、Ｌ_１グラフ計算装置１の機能の分散及び統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散または統合して構成することができる。

また、Ｌ_１グラフ計算装置１においておこなわれる各処理は、全部または任意の一部が、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）、及び、ＣＰＵ，ＧＰＵにより解析実行されるプログラムにて実現されてもよい。また、Ｌ_１グラフ計算装置１においておこなわれる各処理は、ワイヤードロジックによるハードウェアとして実現されてもよい。

また、実施の形態において説明した各処理のうち、自動的におこなわれるものとして説明した処理の全部または一部を手動的に行うこともできる。もしくは、手動的におこなわれるものとして説明した処理の全部または一部を公知の方法で自動的に行うこともできる。この他、上述及び図示の処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて適宜変更することができる。

［プログラム］
図２０は、プログラムが実行されることにより、Ｌ_１グラフ計算装置１が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。これらの各部は、バス１０８０によって接続される。

メモリ１０１０は、ＲＯＭ（Read Only Memory）１０１１及びＲＡＭ１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic Input Output System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０９０に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１１００に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１１００に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１１１０、キーボード１１２０に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１１３０に接続される。

ハードディスクドライブ１０９０は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、Ｌ_１グラフ計算装置１の各処理を規定するプログラムは、コンピュータ１０００により実行可能なコードが記述されたプログラムモジュール１０９３として実装される。プログラムモジュール１０９３は、例えばハードディスクドライブ１０９０に記憶される。例えば、Ｌ_１グラフ計算装置１における機能構成と同様の処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０９０に記憶される。なお、ハードディスクドライブ１０９０は、ＳＳＤ（Solid State Drive）により代替されてもよい。

また、上述した実施の形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０９０に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０９０に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。

なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０９０に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１１００等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３及びプログラムデータ１０９４は、ネットワーク（ＬＡＮ、ＷＡＮ等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３及びプログラムデータ１０９４は、他のコンピュータから、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。

以上、本発明者によってなされた発明を適用した実施の形態について説明したが、本実施の形態による本発明の開示の一部をなす記述及び図面により本発明は限定されることはない。すなわち、本実施の形態に基づいて当業者等によりなされる他の実施の形態、実施例及び運用技術等は全て本発明の範疇に含まれる。

１ Ｌ_１グラフ計算装置
１１ ＳＶＤ計算部
１２ ノード選択部
１３ 第１エッジ集合設定部
１４ 第１重み更新部
１５ パラメータ初期化部
１６ 第２エッジ集合設定部
１７ 追加エッジ計算部
１８ エッジ追加部
１９ 第２重み更新部
２０ 重み設定部

Claims (5)

1. lassoによるＬ_１グラフを計算するＬ_１グラフ計算装置であって、
   入力される多次元行列の特異値分解を計算するＳＶＤ計算部と、
   エッジの重みを計算する対象のノードを選択するノード選択部と、
   初めに重みを更新するエッジの集合を設定する第１エッジ集合設定部と、
   前記初めに重みを更新するエッジ集合に含まれるエッジの重みを収束するまで更新する第１重み更新部と、
   グラフの計算におけるパラメータの初期化を行うパラメータ初期化部と、
   重みを更新するエッジ集合として、前記特異値分解の計算結果を用いて、非零の重みを持ち得ないエッジを除外したエッジ集合を設定する第２エッジ集合設定部と、
   重みを更新するエッジ集合を、前記第２エッジ集合設定部が設定したエッジ集合に追加する追加エッジ計算部と、
   前記追加エッジ計算部が追加したエッジ集合から、非零の重みを有するエッジを、一つ一つ、前記重みを更新するエッジ集合に追加するエッジ追加部と、
   前記エッジ追加部がエッジを追加したエッジ集合に含まれるエッジの重みを、収束するまで更新する第２重み更新部と、
   前記ノード選択部によって選択されていないノードに対してエッジの重みを設定する重み設定部と、
   を有することを特徴するＬ_１グラフ計算装置。
2. 前記第２エッジ集合設定部は、前記特異値分解の計算結果を用いて求めたＫＫＴ（Karush-Kuhn-Tucker）条件のスコアの上限値及び下限値を用いて、前記非零の重みを持ち得ないエッジを特定することを特徴とする請求項１に記載のＬ_１グラフ計算装置。
3. 前記エッジ追加部は、追加したエッジ集合のエッジに対し、lassoの回帰結果における残差を用いた演算式により重みを計算し、前記計算した重みが非零であるエッジを、前記重みを更新するエッジ集合に追加し、前記計算した重みが零であるエッジを、前記重みを更新するエッジ集合から除外することを特徴とする請求項１または２に記載のＬ_１グラフ計算装置。
4. lassoによるＬ_１グラフを計算するＬ_１グラフ計算装置が行うＬ_１グラフ計算方法であって、
   入力される多次元行列の特異値分解を計算するＳＶＤ計算工程と、
   エッジの重みを計算する対象のノードを選択するノード選択工程と、
   初めに重みを更新するエッジの集合を設定する第１エッジ集合設定工程と、
   前記初めに重みを更新するエッジ集合に含まれるエッジの重みを収束するまで更新する第１重み更新工程と、
   グラフの計算におけるパラメータの初期化を行うパラメータ初期化工程と、
   重みを更新するエッジ集合として、前記特異値分解の計算結果を用いて、非零の重みを持ち得ないエッジを除外したエッジ集合を設定する第２エッジ集合設定工程と、
   重みを更新するエッジ集合を、前記第２エッジ集合設定工程において設定されたエッジ集合に追加する追加エッジ計算工程と、
   前記追加エッジ計算工程において追加したエッジ集合から、非零の重みを有するエッジを、一つ一つ、前記重みを更新するエッジ集合に追加するエッジ追加工程と、
   前記エッジ追加工程においてエッジを追加されたエッジ集合に含まれるエッジの重みを、収束するまで更新する第２重み更新工程と、
   前記ノード選択工程において選択されていないノードに対してエッジの重みを設定する重み設定工程と、
   を含んだことを特徴とするＬ_１グラフ計算方法。
5. コンピュータを、請求項１〜３のいずれか一つに記載のＬ_１グラフ計算装置として機能させるためのＬ_１グラフ計算プログラム。

Images