(57)【要約】【課題】lassoを用いてＬ１グラフを高速に計算することができる。【解決手段】Ｌ１グラフ計算装置１は、lassoによるＬ１グラフを計算するＬ１グラフ計算装置であって、入力される多次元行列の特異値分解を計算するＳＶＤ計算部１１と、エッジの重みを計算する対象のノードを選択するノード選択部１２と、グラフの計算におけるパラメータの初期化を行うパラメータ初期化部１５と、ＳＶＤ計算結果を用いて、非零の重みを持ち得ないエッジを除外したエッジ集合を設定する第２エッジ集合設定部１６と、重みを更新するエッジ集合をエッジ集合に追加する追加エッジ計算部１７と、追加したエッジ集合から、非零の重みを有するエッジを、一つ一つエッジ集合に追加するエッジ追加部１８と、エッジ集合に含まれるエッジの重みを収束するまで更新する第２重み更新部１９と、選択されていないノードに対してエッジの重みを設定する重み設定部２０と、を有する。【選択図】図１
2018073285.tif 000002
【特許請求の範囲】
【請求項１】
lassoによるＬ１グラフを計算するＬ１グラフ計算装置であって、
入力される多次元行列の特異値分解を計算するＳＶＤ計算部と、
エッジの重みを計算する対象のノードを選択するノード選択部と、
初めに重みを更新するエッジの集合を設定する第１エッジ集合設定部と、
前記初めに重みを更新するエッジ集合に含まれるエッジの重みを収束するまで更新する第１重み更新部と、
グラフの計算におけるパラメータの初期化を行うパラメータ初期化部と、
重みを更新するエッジ集合として、前記特異値分解の計算結果を用いて、非零の重みを持ち得ないエッジを除外したエッジ集合を設定する第２エッジ集合設定部と、
重みを更新するエッジ集合を、前記第２エッジ集合設定部が設定したエッジ集合に追加する追加エッジ計算部と、
前記追加エッジ計算部が追加したエッジ集合から、非零の重みを有するエッジを、一つ一つ、前記重みを更新するエッジ集合に追加するエッジ追加部と、
前記エッジ追加部がエッジを追加したエッジ集合に含まれるエッジの重みを、収束するまで更新する第２重み更新部と、
前記ノード選択部によって選択されていないノードに対してエッジの重みを設定する重み設定部と、
を有することを特徴するＬ１グラフ計算装置。
【請求項２】
前記第２エッジ集合設定部は、前記特異値分解の計算結果を用いて求めたＫＫＴ（Karush-Kuhn-Tucker）条件のスコアの上限値及び下限値を用いて、前記非零の重みを持ち得ないエッジを特定することを特徴とする請求項１に記載のＬ１グラフ計算装置。
【請求項３】
前記エッジ追加部は、追加したエッジ集合のエッジに対し、lassoの回帰結果における残差を用いた演算式により重みを計算し、前記計算した重みが非零であるエッジを、前記重みを更新するエッジ集合に追加し、前記計算した重みが零であるエッジを、前記重みを更新するエッジ集合から除外することを特徴とする請求項１または２に記載のＬ１グラフ計算装置。
【請求項４】
lassoによるＬ１グラフを計算するＬ１グラフ計算装置が行うＬ１グラフ計算方法であって、
入力される多次元行列の特異値分解を計算するＳＶＤ計算工程と、
エッジの重みを計算する対象のノードを選択するノード選択工程と、
初めに重みを更新するエッジの集合を設定する第１エッジ集合設定工程と、
前記初めに重みを更新するエッジ集合に含まれるエッジの重みを収束するまで更新する第１重み更新工程と、
グラフの計算におけるパラメータの初期化を行うパラメータ初期化工程と、
重みを更新するエッジ集合として、前記特異値分解の計算結果を用いて、非零の重みを持ち得ないエッジを除外したエッジ集合を設定する第２エッジ集合設定工程と、
重みを更新するエッジ集合を、前記第２エッジ集合設定工程において設定されたエッジ集合に追加する追加エッジ計算工程と、
前記追加エッジ計算工程において追加したエッジ集合から、非零の重みを有するエッジを、一つ一つ、前記重みを更新するエッジ集合に追加するエッジ追加工程と、
前記エッジ追加工程においてエッジを追加されたエッジ集合に含まれるエッジの重みを、収束するまで更新する第２重み更新工程と、
前記ノード選択工程において選択されていないノードに対してエッジの重みを設定する重み設定工程と、
を含んだことを特徴とするＬ１グラフ計算方法。
【請求項５】
コンピュータを、請求項１〜３のいずれか一つに記載のＬ１グラフ計算装置として機能させるためのＬ１グラフ計算プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、Ｌ１グラフ計算装置、Ｌ１グラフ計算方法及びＬ１グラフ計算プログラムに関する。
【背景技術】
【０００２】
グラフは基盤的なデータ構造であり、多次元データを解析する上で重要なデータ構造である。グラフ構造としてはｋ−近傍グラフが有名であるが、データのノイズに弱いという問題がある。この問題を解決するために、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によるＬ１グラフ計算では、ωｐ［ｐ］=０として、以下の式（１）に示す目的関数を最小化するように重みを計算する。
【００１３】
【数１】
2018073285.tif 000003
【００１４】
式（１）において最初の項は回帰のよる２乗誤差に対応し、２番目の項は重みのＬ１制約に対応する。式（１）はチューニングパラメータλが大きくなるとグラフがより疎な構造を持つことを示している。もし行列Χがフル行ランクであれば式（１）は一つの解を持ち、そうでなければ解は一つにならない。
【００１５】
lassoを高速に解く方法として座標降下法がある。座標降下法は、収束するまでエッジの重みを一つ一つ更新するものである。具体的に、座標降下法では、以下の式（２）により重みの更新を行う。
【００１６】
【数２】
2018073285.tif 000004
【００１７】
ここでＳ［・，・］は、以下の式（３）を用いて計算する。
【００１８】
【数３】
2018073285.tif 000005
【００１９】
さらに、ｚｐ［ｕ｜ｗｐ］は、ノードｕのベクトルｗｐのパラメータであり、以下の式（４）を用いて計算する。
【００２０】
【数４】
2018073285.tif 000006
【００２１】
【数５】
2018073285.tif 000007
【００２２】
上記の式にみられるように、各パラメータは、以下の性質を有する。
【００２３】
【数６】
2018073285.tif 000008
【００２４】
結果としてＴを収束までの繰り返し計算回数としたときに、式（２）によりノードｐへの重みを計算するにはＯ（ｄｐＭＴ）の計算コストが必要になる。したがって、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：特異値分解）計算部１１、ノード選択部１２、第１エッジ集合設定部１３、第１重み更新部１４、パラメータ初期化部１５、第２エッジ集合設定部１６、追加エッジ計算部１７、エッジ追加部１８、第２重み更新部１９及び重み設定部２０を有する。
【００３４】
ＳＶＤ計算部１１は、入力される多次元行列の特異値分解を計算する。ノード選択部１２は、エッジの重みを計算する対象のノードを選択する。そして、ノード選択部１２の後段の各部は、この計算対象のノードに対して重み計算を繰り返し行う。第１エッジ集合設定部１３は、初めに重みを更新するエッジの集合を設定する。第１重み更新部１４は、第１エッジ集合設定部１３によって設定されたエッジ集合に含まれるエッジの重みを収束するまで更新する。パラメータ初期化部１５は、グラフの計算におけるパラメータの初期化を行う。
【００３５】
そして、第２エッジ集合設定部１６は、重みを更新するエッジの集合として、特異値分解の計算結果を用いて、非零の重みを持ち得ないエッジを除外したエッジ集合を設定する。追加エッジ計算部１７は、重みを更新するエッジ集合を、第２エッジ集合設定部１６が設定したエッジ集合に追加する。エッジ追加部１８は、追加エッジ計算部１７が追加したエッジ集合から、非零の重みを有するエッジを、一つ一つ、重みを更新するエッジ集合に追加する。第２重み更新部１９は、エッジ追加部１８がエッジを追加したエッジ集合に含まれるエッジの重みを、収束するまで更新する。重み設定部２０は、ノード選択部１２によって選択されていないノードに対してエッジの重みを設定する。これらの第２エッジ集合設定部１６、追加エッジ計算部１７、エッジ追加部１８、第２重み更新部１９及び重み設定部２０は、エッジの重み更新についての繰り返し計算を行う。
【００３６】
［エッジ集合設定処理］
ここで、実施の形態１に係るＬ１グラフ計算方法のうち、エッジの重み更新についての繰り返し計算のうちのエッジ集合設定処理について説明する。まず、重みの更新のための計算を前に、第２エッジ集合設定部１６が、２つのエッジ集合Ｍ［ｗｐ］とエッジ集合Ｃ［ｗｐ］とを計算する。エッジ集合Ｍ［ｗｐ］は、重みベクトルｗｐを用いて更新を行ったときに、重みが非零に必ずなるエッジの集合である。また、エッジ集合Ｃ［ｗｐ］は、重みベクトルｗｐを用いて更新したときに重みが非零になる可能性のあるエッジの集合である。
【００３７】
理論的に、エッジ集合Ｍ［ｗｐ］とエッジ集合Ｃ［ｗｐ］とはＫＫＴ条件に基づいている。したがって、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗｐ］とエッジ集合Ｃ［ｗｐ］とを、ＫＫＴ条件のスコアの上限値と下限値とを用いて求める。第２エッジ集合設定部１６は、ＳＶＤ計算結果を用いて求めたＫＫＴ条件のスコアの上限値と下限値とを用いて、非零の重みを持ち得ない不要なエッジを特定し、このエッジを枝狩りする。
【００３８】
そこで、エッジ集合の定義を述べる前に、ＫＫＴ条件について述べる。ＫＫＴ条件は、図２に示す定義１のように定義される。図２は、定義１を示す図である。
【００３９】
ｕ∈Ｖ／ｐであるようなノードｕに対して、もし、ノードｕからｐのエッジの重みが重みベクトルｗｐで更新されたとき、Ｋ［ｕ｜ｗｐ］は、更新後の、ｗｐ［ｕ］の値により、以下の式（７）のように計算できる。
【００４０】
【数７】
2018073285.tif 000009
【００４１】
直接ＫＫＴ条件のスコアを計算すれば重みが非零になるか否かが分かる。ただし、ＫＫＴ条件のスコアを計算するためには、最長でＯ（ＮＭ）の計算コストがかかる。これは、式（Ａ）における行ベクトルｘｐ、重みベクトルｗｐ、行列Χ、ｘｕの転置ｘｕＴの大きさが、１×Ｍ、１×Ｎ、Ｎ×Ｍ、Ｍ×１であるためである。そこで、本実施の形態１では、第２エッジ集合設定部１６が、ＫＫＴ条件のスコアの上限値と下限値との計算を高速化できるように定義づけを行っている。なお、ＫＫＴ条件のスコアの上限値と下限値との計算方法は後述する。
【００４２】
続いて、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗｐ］について、図３に示す定義２のように計算する。図３は、定義２を示す図である。第２エッジ集合設定部１６は、エッジ集合Ｃ［ｗｐ］について、図４に示す定義３のように計算する。図４は、定義３を示す図である。そして、エッジ集合については、補題１，２が成り立つ。図５は、補題１を示す図である。また、図６は、補題２を示す図である。
【００４３】
エッジ集合Ｍ［ｗｐ］とエッジ集合Ｃ［ｗｐ］との関係については、図７に示す補題３が成り立つ。図７は、補題３を示す図である。このように、第２エッジ集合設定部１６は、定義２，３を満たすエッジ集合Ｍ［ｗｐ］とエッジ集合Ｃ［ｗｐ］を求める。言い換えると、第２エッジ集合設定部１６は、非零の重みを持ち得ない不要なエッジを除外したエッジ集合Ｍ［ｗｐ］とエッジ集合Ｃ［ｗｐ］を求める。
【００４４】
［重み更新処理］
続いて、エッジの重み更新についての繰り返し計算のうちの重み更新処理における計算方法について述べる。Ｌ１グラフ計算装置１では、二つの計算方法を用いる。まず、初めの計算方法（第１の計算方法）について説明する。初めの計算方法として、ノードｕのベクトルｗｐのパラメータであるパラメータｚｐ［ｕ｜ｗｐ］を計算し、エッジごとに異なる解析結果を使わずともよい方法を用いる。具体的には、以下の式（８）を用いて計算を行う。
【００４５】
【数８】
2018073285.tif 000010
【００４６】
ここでｒｐ［ｉ］は、要素χｐ［ｉ］のlassoの回帰結果における残差であり、以下の式（９）を用いて計算される。
【００４７】
【数９】
2018073285.tif 000011
【００４８】
式（９）で使用される残差は、各重みに対して同じである。このため、もし回帰結果が更新後に変わることがなければ、パラメータｚｐ［ｕ｜ｗｐ］を高速に計算することができる。パラメータｚｐ［ｕ｜ｗｐ］を計算するために式（８）において必要になる計算コストは、Ｏ（Ｍ）である。この式（８）に対し、以下の補題４に示す性質が成り立つ。図８は、補題４を示す図である。
【００４９】
エッジ追加部１８は、重み更新の前に、もしエッジが零の重みを持てば、式（８）を用いて重みベクトルに関するパラメータｚｐ［ｕ｜ｗｐ］を計算する。この零の重みを持つエッジは、更新後も、再度、零の重みを持つことが期待されるため、残差は、更新の前後で同じ値を持つことが期待できる。このため、次の繰り返し計算においても、効果的に残差ｒｐ［ｉ］を利用することができる。言い換えると、この零の重みを持つエッジは、更新後も、再度、零の重みを持つため、次の重み更新についても繰り返し計算の対象から除外することができる。
【００５０】
一方、エッジが非零の重みを更新の後で持つ場合であっても、第２重み更新部１９は、式（９）を用いることで、逐次的に、残差ｒｐ［ｉ］をＯ（１）の計算コストで更新することができる。
【００５１】
そして、第２重み更新部１９は、もしエッジが更新の前に非零の重みを持つ場合には、既存研究で提案されている以下の内積を用いる式（１０）を用いて、重みを更新する。
【００５２】
【数１０】
2018073285.tif 000012
【００５３】
次に、重み更新についての繰り返し計算において、更新対象となるエッジ集合について説明する。前述したように、Ｌ１グラフ計算装置１では、重み更新についての繰り返し計算の前に、エッジ集合Ｍ［ｗｐ］とエッジ集合Ｃ［ｗｐ］とを計算する。Ｕを、更新を行うエッジの集合としたときに、素朴にはＵ＝Ｍ［ｗｐ］またはＵ＝Ｃ［ｗｐ］とすることにより更新を行うエッジを減らすことができる。
【００５４】
ただし、Ｌ１グラフ計算装置１では、ＫＫＴ条件の上限値と下限値とを用いてエッジ集合Ｍ［ｗｐ］とエッジ集合Ｃ［ｗｐ］とを求めているため、エッジ集合には重みが零になるエッジが含まれる。そして、Ｌ１グラフ計算装置１では、重みが零になるエッジに対して繰り返し計算が行われることがあり得る。
【００５５】
そこで、高速に重みの更新を行うために、Ｌ１グラフ計算装置１では、エッジ集合Ａからエッジ集合Ｕへ一つ一つエッジを足す処理を行う。Ｌ１グラフ計算装置１では、まず、第１重み更新部１４が、エッジ集合Ｕを非零の重みを持つエッジから求め、そのエッジ集合Ｕに対して収束するまで重みの更新を行う。収束後、第２エッジ集合設定部１６は、エッジ集合Ｍ［ｗｐ］を計算し、エッジを加える集合ＡをＡ＝Ｍ［ｗｐ］＼Ｕとして計算する。
【００５６】
このエッジ集合Ａのエッジに対して、追加エッジ計算部１７、エッジ追加部１８及び第２重み更新部１９が、一つ一つ重みの更新を行う。そして、もし、あるエッジが非零の重みを持つ場合には、エッジ追加部１８が、そのエッジをエッジ集合Ｕへ追加し、第２重み更新部１９が、集合Ｕについて、エッジの重みが収束するまで重み更新について繰り返し計算を行う。
【００５７】
一方、あるエッジが非零の重みを持たない場合、エッジ追加部１８は、このエッジに対して、繰り返し計算をすることなく枝狩りを行う。同様に、第２エッジ集合設定部１６は、Ａ＝Ｃ［ｗｐ］＼Ｕを計算し、追加エッジ計算部１７、エッジ追加部１８及び第２重み更新部１９が、このエッジ集合Ａのエッジに対して重みの更新を行う。
【００５８】
Ｌ１グラフ計算装置１は、もしエッジが非零の重みを持たなければ、そのエッジの重み更新のための繰り返し計算を行わない。このため、Ｌ１グラフ計算装置１では、全てのエッジについてエッジの重み更新のための繰り返し計算を行っていた従来装置と比較し、高速な計算処理が可能になる。理論的には、この手法は、以下の補題５に示す性質に基づいている。図９は、補題５を示す図である。
【００５９】
補題５にあるように、ノードｕからｐへのエッジの重みの更新処理において、その重み更新についての繰り返し計算において、零であり座標降下法により収束する場合には、この結果は、座標降下法による解析結果に影響しない。したがって、重み更新についての繰り返し計算において零であるとともに座標降下法により収束するエッジについては、このエッジを枝狩りできることを示している。上記のように、Ｌ１グラフ計算装置１は、エッジが非零の重みを持たない場合には、そのエッジの繰り返し計算を行わないため、全てのエッジについてエッジの重み更新のための繰り返し計算を行っていた従来装置と比較して高速な処理が可能になる。
【００６０】
［ＫＫＴスコアの上限値と下限値とを求めるための計算処理］
次に、ＫＫＴスコアの上限値と下限値とを求めるための計算処理について説明する。上限値と下限値を計算する方法として、Ｌ１グラフ計算装置１は、ＳＶＤを用いる方法と、逐次的に更新する方法とを用いる。
【００６１】
まず、ＳＶＤを用いる方法について説明する。この方法では、行列Χを、ＳＶＤを用いて、図１０に示すように近似する。図１０は、行列ΧのＳＶＤを用いた近似方法を示す図である。そして、図１１は、図１０に示す定義４を示す図である。そして、ＫＫＴ条件のスコアの上限値と下限値とに対し、以下の補題６に示す性質が成り立つ。図１２は、補題６を示す図である。
【００６２】
次に、逐次的に上限値と下限値とを計算する方法について説明する。前述したように、Ｌ１グラフ計算装置１では、第２エッジ集合設定部１６、追加エッジ計算部１７、エッジ追加部１８及び第２重み更新部１９が、エッジ集合Ａからエッジを一つ一つエッジ集合Ｕに足しながらエッジの重みを更新する。この時、ＳＶＤによる上限値と下限値とを用いることにより高速にエッジ集合Ｍ［ｗｐ］とＣ［ｗｐ］を計算することができる。ここで、さらに、高速にエッジ集合を計算するために、第２エッジ集合設定部１６は、逐次的にＫＫＴ条件のスコアの上限値と下限値とを計算する方法を用いる。
【００６３】
逐次的に上限値と下限値を更新する方法を述べる前に、まず、ＫＫＴ条件とパラメータｚｐ［ｕ｜ｗｐ］との関係について説明する。図１３は、補題７を示す図である。図１３の補題７は、ＫＫＴ条件とパラメータｚｐ［ｕ｜ｗｐ］との関係を示すものである。
【００６４】
補題７は、もし（１）あるエッジがエッジ集合に含まれ、（２）そのエッジのパラメータｚｐ［ｕ｜ｗｐ］が既に計算済みであれば、ＫＫＴ条件のスコアはＯ（１）の計算コストで求めることができることを示している。第２エッジ集合設定部１６は、エッジを追加した後、上限値と下限値とを、以下の定義５に示すように逐次的に計算する。図１４は、定義５を示す図である。
【００６５】
そして、ＫＫＴ条件のスコアの上限値及び下限値については、補題８に示す性質が成り立つ。図１５は、補題８を示す図である。
【００６６】
前述したように、Ｌ１グラフ計算装置１は、高速に重みを更新するために追加されたエッジが回帰結果を変えることがあるたびに回帰の残差を計算する。座標降下法は、重みを一つずつ更新するため、式（Ｈ）（図１４参照）から、δｐｉとΔｐをＯ（１）の計算コストで更新することができる。Ｌ１グラフ計算装置１は、結果として定義５を用いることによって、ＫＫＴ条件のスコアの上限値と下限値とを逐次的にＯ（１）の計算コストで更新することができる。
【００６７】
なお、Ｌ１グラフ計算装置１において、ＫＫＴ条件のスコアの上限値と下限値を計算するという観点では、定義５で与えられる逐次的な方法は、定義４で与えられるＳＶＤを用いる方法と同じである。ただし、これらは、それぞれの計算方法がまったく異なるため、異なる上限値と下限値となる。さらに、逐次的な方法は、ＳＶＤによる方法より高速に上限値と下限値を計算することができる。このため、Ｌ１グラフ計算装置１は、まず逐次的な方法でＫＫＴ条件のスコアの上限値と下限値を計算し、それでも枝狩りできない場合にＳＶＤによる方法を用いて上限値と下限値を計算する。
【００６８】
［アルゴリズム］
図１６を参照して、Ｌ１グラフ計算装置１の各部の処理について説明する。図１６は、実施の形態１に係るＬ１グラフ計算処理のアルゴリズムを示す図である。図１６に示す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）により代替されてもよい。
【００９５】
また、上述した実施の形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０９０に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０９０に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
【００９６】
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０９０に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１１００等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３及びプログラムデータ１０９４は、ネットワーク（ＬＡＮ、ＷＡＮ等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３及びプログラムデータ１０９４は、他のコンピュータから、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【００９７】
以上、本発明者によってなされた発明を適用した実施の形態について説明したが、本実施の形態による本発明の開示の一部をなす記述及び図面により本発明は限定されることはない。すなわち、本実施の形態に基づいて当業者等によりなされる他の実施の形態、実施例及び運用技術等は全て本発明の範疇に含まれる。
【符号の説明】
【００９８】
１  Ｌ１グラフ計算装置
１１  ＳＶＤ計算部
１２  ノード選択部
１３  第１エッジ集合設定部
１４  第１重み更新部
１５  パラメータ初期化部
１６  第２エッジ集合設定部
１７  追加エッジ計算部
１８  エッジ追加部
１９  第２重み更新部
２０  重み設定部
【図１】
2018073285.tif 000013
【図２】
2018073285.tif 000014
【図３】
2018073285.tif 000015
【図４】
2018073285.tif 000016
【図５】
2018073285.tif 000017
【図６】
2018073285.tif 000018
【図７】
2018073285.tif 000019
【図８】
2018073285.tif 000020
【図９】
2018073285.tif 000021
【図１０】
2018073285.tif 000022
【図１１】
2018073285.tif 000023
【図１２】
2018073285.tif 000024
【図１３】
2018073285.tif 000025
【図１４】
2018073285.tif 000026
【図１５】
2018073285.tif 000027
【図１６】
2018073285.tif 000028
【図１７】
2018073285.tif 000029
【図１８】
2018073285.tif 000030
【図１９】
2018073285.tif 000031
【図２０】
2018073285.tif 000032