JP6453785B2

JP6453785B2 - 回帰分析装置、回帰分析方法および回帰分析プログラム

Info

Publication number: JP6453785B2
Application number: JP2016010117A
Authority: JP
Inventors: 靖宏藤原
Original assignee: Nippon Telegraph and Telephone Corp
Current assignee: Nippon Telegraph and Telephone Corp
Priority date: 2016-01-21
Filing date: 2016-01-21
Publication date: 2019-01-16
Anticipated expiration: 2036-01-21
Also published as: JP2017130116A

Description

本発明は、回帰分析装置、回帰分析方法および回帰分析プログラムに関する。

高次元データに対するｌ１正則化付き最小２乗法の代表的な手法としてlassoが知られている（例えば非特許文献１）。機械学習において、応答変数と関係を持つ予測変数を求める処理が重要である。このとき、良い予測変数には、応答変数と高い相関がありながら予測変数同士は相関が低いということが求められる。また、できる限り少ない予測変数で応答変数に対する高い予測精度を達成することが重要となる。

lassoは、ｌ１正則化項を２乗誤差項に足すことで疎な線形回帰を実現できることが優れている。lassoは、疎な予測変数の集合を求めるために、ｌ１正則化項による制約を効果的に利用し予測変数を削減しているのである。

lassoは１９９０年代の中頃に開発されたが、高い計算コストが必要であるため、２０００年前半までは大きな注目を集めなかった。計算コストが高かったのは、lassoが提案された論文では、大規模データの処理には適していない手法が予測変数を求めるのに使われていたからである。

そして、２００２年に、Tibshiraniらがlassoを処理するのにLARSというアルゴリズムを提案した。その結果、lassoを計算することが容易になり、lassoの関連手法であるelastic Netやgrouped Lasso等が開発された。

さらに、２００７年には、Tibshiraniらの研究グループが、LARSより高速な手法としてcoordinate descentに基づく手法を提案した（例えば非特許文献２）。coordinate descentに基づく手法は、予測変数の係数を１つずつ繰り返し更新し、lassoにおける解を求める。また、２０１０年および２０１２年における論文で、coordinate descentをさらに高速化する手法が提案された（例えば非特許文献３および非特許文献４）。lassoを用いるにあたり、coordinate descentに基づく手法は現在よく使われている。

Robert Tibshirani, Regression Shrinkage and Selection via the Lasso, Journal of the Royal Statistical Society, Series B, 1996 Jerome Friedman and Trevor Hastie and Holger Hofling and Robert Tibshirani, Pathwise coordinate optimization, Annals of Applied Statistics, 2007 Jerome H. Friedman and Trevor Hastie and Rob Tibshirani, Regularization Paths for Generalized Linear Models via Coordinate Descent, Journal of Statistical Software, 2010 Tibshirani, Robert and Bien, Jacob and Friedman, Jerome and Hastie, Trevor and Simon, Noah and Taylor, Jonathan and Tibshirani, Ryan J., Strong rules for discarding predictors in lasso-type problems, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2012

しかしながら、従来のlassoを用いた線形回帰の手法には、大きなサイズのデータを高速に処理することができない場合があるという課題があった。

例えば、lassoが提案された１９９０年の中頃、lassoは前立腺がんの解析に利用されたが、予測変数の数は高々１０程度であった。一方、最近の画像処理における予測変数の数は数千であり、さらに話題抽出のアプリケーションにおいては数万の数の予測変数が使われている。このように、近年lassoが扱うデータは非常に大きなサイズになっているため、全ての予測変数に対して係数の更新を行う従来の手法では、計算量が膨大になり、高速な処理を行うことができない場合があった。

本発明の回帰分析装置は、lassoを用いて、予測変数の持つ係数を繰り返し更新し回帰分析を行う回帰分析装置であって、前記lassoにおけるチューニングパラメータを順次設定するパラメータ設定部と、前記パラメータ設定部によって前記チューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する初期更新対象設定部と、前記パラメータ設定部によって前記チューニングパラメータが設定されるたびに、前記係数の初期値を設定する初期値設定部と、前記係数を基に参照ベクトルを設定する第１参照ベクトル設定部と、前記チューニングパラメータと、前記参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により必ず非零になるか否かを判定する非零係数判定部と、前記係数のうち、前記非零係数判定部によって必ず非零になると判定された係数を、収束するまで繰り返し更新する第１係数更新部と、前記第１係数更新部によって更新された前記係数を基に参照ベクトルをさらに設定する第２参照ベクトル設定部と、前記チューニングパラメータと、前記参照ベクトルを基にさらに計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により非零になる可能性があるか否かを判定する零係数判定部と、前記係数のうち、前記零係数判定部によって非零になる可能性があると判定された係数を、収束するまで繰り返し更新する第２係数更新部と、前記予測変数の集合に含まれる予測変数がＫＫＴ条件を満たすか否かを判定し、前記予測変数の全てがＫＫＴ条件を満たすまで係数の更新を繰り返し実行させるＫＫＴ条件判定部と、を有することを特徴とする。

また、本発明の回帰分析方法は、lassoを用いて、予測変数の持つ係数を繰り返し更新し回帰分析を行う回帰分析方法であって、前記lassoにおけるチューニングパラメータを順次設定するパラメータ設定工程と、前記パラメータ設定工程によって前記チューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する初期更新対象設定工程と、前記パラメータ設定工程によって前記チューニングパラメータが設定されるたびに、前記係数の初期値を設定する初期値設定工程と、前記係数を基に参照ベクトルを設定する第１参照ベクトル設定工程と、前記チューニングパラメータと、前記参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により必ず非零になるか否かを判定する非零係数判定工程と、前記係数のうち、前記非零係数判定工程によって必ず非零になると判定された係数を、収束するまで繰り返し更新する第１係数更新工程と、前記第１係数更新工程によって更新された前記係数を基に参照ベクトルをさらに設定する第２参照ベクトル設定工程と、前記チューニングパラメータと、前記参照ベクトルを基にさらに計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により非零になる可能性があるか否かを判定する零係数判定工程と、前記係数のうち、前記零係数判定工程によって非零になる可能性があると判定された係数を、収束するまで繰り返し更新する第２係数更新工程と、前記予測変数の集合に含まれる予測変数がＫＫＴ条件を満たすか否かを判定し、前記予測変数の全てがＫＫＴ条件を満たすまで係数の更新を繰り返し実行させるＫＫＴ条件判定工程と、を含んだことを特徴とする。

本発明によれば、大きなサイズのデータに対し、lassoを用いた線形回帰を高速に行うことができる。

図１は、第１の実施形態に係る回帰分析装置の構成の一例を示すブロック図である。図２は、第１の実施形態に係る回帰分析処理のアルゴリズムを示す図である。図３は、定義１を示す図である。図４は、定義２を示す図である。図５は、補題１を示す図である。図６は、補題２を示す図である。図７は、補題３を示す図である。図８は、補題４を示す図である。図９は、補題５を示す図である。図１０は、補題６を示す図である。図１１は、定理１を示す図である。図１２は、定理２を示す図である。図１３は、第１の実施形態に係る回帰分析処理の流れを示すフローチャートである。図１４は、プログラムが実行されることにより、回帰分析装置が実現されるコンピュータの一例を示す図である。

以下に、本願に係る回帰分析装置、回帰分析方法および回帰分析プログラムの実施形態を図面に基づいて詳細に説明する。なお、この実施形態により本発明が限定されるものではない。

［従来のlassoを用いた回帰分析について］
まず、lassoを用いた回帰分析を行う従来の手法について具体的に説明する。回帰分析において、予測変数の数はｐとし、観測値の数はｎとする。また、全てのベクトルは平均を０、分散を１に正規化されているものとする。また、応答変数は、式（１）のベクトルで表される。

また、予測変数は、式（２）の行列で表される。

lassoにおいては、式（２）の予測変数の行列における最小２乗誤差とｌ１正則化項の制限からなる以下の式を最小化する応答変数を予測するモデルを求める。

このとき、式（２）の予測変数の行列が列フルランク行列であれば、上記の最適化問題の解はただ１つになる。そうでない場合、解は１つに定まらない。なおｐ＞ｎである場合は、式（２）の予測変数の行列は、列フルランク行列にならない。

また、実際のアプリケーションにおいて最適なλの値は異なる。そのため、Ｋをチューニングパラメータの数としたときに、式（３）においてチューニングパラメータをλ_１＞λ_２＞...λ_Ｋと変化させる。そして、チューニングパラメータを変化させつつ最適化問題を解き、複数のチューニングパラメータを試した後に、事前に決定した赤池情報量等の基準を用いて最終的なlassoにおける解を決定する。

Tibshiraniらは、coordinate descentを用いて予測変数を１つずつ更新する手法を提案した。coordinate descentは他の係数はすでに更新したことを仮定し、係数ｗ［ｉ］に対して部分的に最適化を行う手法である。coordinate descentでは、まず式（４）が計算される。

さらに、coordinate descentにおいては、式（５）により係数が更新される。

式（５）のＳ（ｚ［ｉ］，λ）はsoft-thresholding operatorと呼ばれる処理であり、またｚ［ｉ］はｉ番目の予測変数に対応するパラメータとして式（６）により計算される。

coordinate descentにおいては、係数が総当たり的に更新される。すなわち、式（５）を用いて収束するまで全ての予測変数の係数が繰り返し更新される。そこで、Tibshiraniらは、更新において計算されるパラメータｚ［ｉ］を式（７）のように計算することで高速に求める手法を提案した。

式（６）と式（７）は同じ計算結果となる。また、ｍを非零の係数を持つ予測変数の数とすると、係数を更新するのに必要な計算コストは、式（７）ではＯ（ｍ）であり、式（６）ではＯ（ｎ）となる。lassoでは、疎に予測変数が選択されるため、ｍ≪ｎとなる。そのため、式（７）を式（６）の代わりに用いることによって、高速に係数を更新することができる。

式（７）を用いるためには、繰り返し計算に入る前に応答変数と全ての予測変数との内積を計算する必要がある。さらに、ある予測変数が応答変数を予測するのに用いられると判断された場合、その予測変数とその他全ての応答変数との内積を計算する必要がある。すなわち、ある予測変数が解として非零の係数を持つ場合、式（７）を利用するために、
その予測変数における内積を計算しなければならない。

そこで、lassoを高速に計算するために、Tibshiraniらを含む多くの研究者によって、係数が０になる予測変数を特定するscreening techniqueが提案されている。特にsequential strong ruleにおいては、式（８）を用いた以下の方法によって予測変数の枝刈りが行われる。

sequential strong ruleは、誤って係数が０にならない予測誤差を枝刈りすることがあるため、収束後に全ての予測変数に対してＫＫＴ（Karush-Kuhn-Tucker）条件が満たされているかを確認する。ＫＫＴ条件は、式（９）を用いて以下のように表される。なお、ＫＫＴ条件は、式（７)を用いることにより、Ｏ（ｍ）のコストで計算できる。

［第１の実施形態の概要］
従来の手法では、screening techniqueで枝刈りされなかった全ての予測変数に対して収束するまで繰り返し係数の更新が行われる。このとき、従来の手法では、ある予測変数の係数の更新を行うのにその他の予測変数の係数で非零であるものを使う。そのため、もし繰り返し計算の中である予測変数が非零の係数を持てば、その係数が繰り返し計算の後０に収束したとしても、計算コストの増加につながってしまう。

提案手法、すなわち第１の実施形態において用いられる手法も、従来の手法と同様にcoordinate descentに基づく手法である。しかし、第１の実施形態においては、従来の手法と異なり、選択的に係数の更新を行う。すなわち、第１の実施形態における手法は、まず非零の係数を持つ予測変数のみを収束するまで更新し、その後に非零の係数を持つことが期待される予測変数を収束するまで更新する。この結果、第１の実施形態における手法は、繰り返し計算において係数が０になる予測変数に対する更新を枝刈りできるため、高速にlassoにおける解を計算することができる。

このように、screening techniqueが繰り返し計算に入る前に係数が０である予測変数を枝刈りするのに対し、第１の実施形態においては、繰り返し計算において係数が０である予測変数を枝刈りする。ｐを予測変数の数、ｎを観測値の数としたとき、与えられるデータはｎ×ｐの行列で表現することができる。もし、当該行列が列フルランクであれば、lassoにおける解は凸関数性から１つになるため、第１の実施形態では、従来の手法と同じ結果を出すことが保証される。一方、当該行列が列フルランクでない場合、第１の実施形態では勾配が０である予測変数の更新を行わないため、従来の手法より少ない予測変数で、より精度の高い予測を行うことができる。

［第１の実施形態の構成］
まず、図１を用いて、第１の実施形態に係る回帰分析装置の構成について説明する。図１は、第１の実施形態に係る回帰分析装置の構成の一例を示すブロック図である。図１に示すように、回帰分析装置１は、パラメータ設定部１０１、初期更新対象設定部１０２、初期値設定部１０３、第１参照ベクトル設定部１０４、非零係数判定部１０５、第１係数更新部１０６、第２参照ベクトル設定部１０７、零係数判定部１０８、第２係数更新部１０９、ＫＫＴ条件判定部１１０、および回帰分析出力部１１１を有する。

図２を用いて、回帰分析装置１の各部の処理について説明する。図２は、第１の実施形態に係る回帰分析処理のアルゴリズムを示す図である。図２に示すAlgorithm１において、チューニングパラメータの設定（例えばAlgorithm１の２行目）、式（７）を用いた係数の更新、ＫＫＴ条件の判定（例えばAlgorithm１の３１行目）等については、前述の従来の手法と同様の処理を適用することができる。一方、係数から計算される所定のパラメータの上限値および下限値を用いて予測変数の枝刈りを行う処理（例えばAlgorithm１の１５、２４行目）は、従来の手法には含まれない処理である。

図２に示すように、Algorithm１においては、Ｕを繰り返し計算において更新を行う予測変数の集合、Ｐを全ての予測変数の集合、Ｐ_ｋをチューニングパラメータλ_ｋへの解として非零の係数となる予測変数の集合、Ｓ_ｋをチューニングパラメータλ_ｋにおけるsequential strong ruleで枝刈りされなかった予測変数の集合とする。

回帰分析装置１は、ｋを１，２，…，Ｋと設定し、以下の処理を行う（Algorithm１の１行目）。まず、パラメータ設定部１０１は、lassoにおけるチューニングパラメータを順次設定する。具体的に、パラメータ設定部１０１は、λ：＝λ_ｋとする（Algorithm１の２行目）。

そして、初期更新対象設定部１０２は、パラメータ設定部１０１によってチューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する。
具体的に、初期更新対象設定部１０２は、ｋ＝１であれば、Ｕを空集合とし、そうでなければＵ：＝Ｐ_ｋ−１とする（Algorithm１の３〜６行目）。

また、初期値設定部１０３は、パラメータ設定部１０１によってチューニングパラメータが設定されるたびに、係数の初期値を設定する。具体的に、初期値設定部１０３は、チューニングパラメータλ_ｋとλ_ｋ−１に対する解は似ているという知見に基づき、チューニングパラメータλ_ｋに対する係数の初期値を、式（１０）により決定する（Algorithm１の７行目）。

次に、第１参照ベクトル設定部１０４は、係数の初期値を基に参照ベクトルを設定する。具体的に、第１参照ベクトル設定部１０４は、繰り返し計算に入る前に、参照ベクトルを設定する（Algorithm１の１０行目）。

そして、非零係数判定部１０５は、チューニングパラメータと、参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、係数のそれぞれが更新により必ず非零になるか否かを判定する。

具体的に、非零係数判定部１０５は、上限値を図３に示す定義１のように定義する。図３は、定義１を示す図である。図３の式（Ａ）におけるベクトルおよびパラメータは以下のように説明することができる。

また、係数ベクトルが各繰り返しで更新される長さｐのベクトルであること等から、定義１は、以下の性質を持つ。

また、非零係数判定部１０５は、下限値を図４に示す定義２のように定義する。図４は、定義２を示す図である。上限値および下限値が定義１および定義２で定義されることは、図５および図６に示す補題１および補題２によって示される。図５は、補題１を示す図である。また、図６は、補題２を示す図である。

式（１３）は、Ｏ（１）で更新できることは明らかであるため、上限値と下限値の計算においては、図７に示す補題３が成り立つ。図７は、補題３を示す図である。ここで、必ず非零の係数を持つ予測変数は、図８の補題４に示す性質を持つ。図８は、補題４を示す図である。また、必ず係数が０になる予測変数は、図９の補題５に示す性質を持つ。図９は、補題５を示す図である。また、補題５より、非零の係数を持つ可能性のある予測変数は、図１０の補題６に示す性質を持つ。図１０は、補題６を示す図である。

予測変数が補題４の条件を満たすとき、その予測変数は必ず補題６の条件を満たす。すなわち、補題４による予測変数の集合は、必ず補題６による予測変数の集合に含まれる。これは、もし下限値がチューニングパラメータより小さければ、上限値もチューニングパラメータより小さく、上限値がチューニングパラメータを負にした値より小さければ、下限値もチューニングパラメータを負にした値より小さいからである。

そのため、もし予測変数において補題４より係数が非零となるとき、その予測変数は補題６より非零の係数を持つ可能性があるとされる。さらに、予測変数が補題６の条件を満たさないとき、補題６よりその予測変数の係数は０となる。第１の実施形態ではこの性質を用いて更新対象となる係数の枝刈りを行っている。

非零係数判定部１０５は、定義１および定義２により、Ｕに含まれる予測変数について、係数パラメータｚ［ｉ］、ｚ［ｉ］の上限値および下限値を計算し、上限値、下限値およびチューニングパラメータを用いてAlgorithm１の１５行目に示す条件が満たされるか否かを判定する（Algorithm１の１１〜１５行目）。

そして、第１係数更新部１０６は、係数のうち、非零係数判定部１０５によって必ず非零になると判定された係数を、収束するまで繰り返し更新する。具体的に、Algorithm１の１５行目に示す条件が満たされた場合、第１係数更新部１０６は、式（７）および式（１３）を用いて係数およびパラメータを更新する（Algorithm１の１６、１７行目）。なお、第１係数更新部１０６は、係数ベクトルが収束するまで更新処理を繰り返す（Algorithm１の１８行目）。

第２参照ベクトル設定部１０７は、第１係数更新部１０６によって更新された係数を基に参照ベクトルをさらに設定する。具体的に、第２参照ベクトル設定部１０７は、第１係数更新部１０６による処理の終了後、次の繰り返し計算に入る前に、参照ベクトルを式（１１）のように設定する（Algorithm１の１９行目）。

零係数判定部１０８は、チューニングパラメータと、参照ベクトルを基にさらに計算された上限値および下限値と、が所定の条件を満たすか否かによって、係数のそれぞれが更新により非零になる可能性があるか否かを判定する。

具体的に、零係数判定部１０８は、定義１および定義２により、Ｕに含まれる予測変数について、係数パラメータｚ［ｉ］およびｚ［ｉ］の上限値および下限値を計算し、上限値、下限値およびチューニングパラメータを用いてAlgorithm１の２４行目に示す条件が満たされるか否かを判定する（Algorithm１の２０〜２４行目）。

そして、第２係数更新部１０９は、係数のうち、零係数判定部１０８によって非零になる可能性があると判定された係数を、収束するまで繰り返し更新する。具体的に、Algorithm１の２４行目に示す条件が満たされた場合、第２係数更新部１０９は、式（７）を用いて係数を更新する（Algorithm１の２５行目）。また、Algorithm１の２４行目に示す条件が満たされなかった場合、第２係数更新部１０９は、係数を０に更新する（Algorithm１の２７行目）。その後、第２係数更新部１０９は、式（１３）を用いてパラメータを更新する（Algorithm１の２８行目）。なお、第２係数更新部１０９は、係数ベクトルが収束するまで更新処理を繰り返す（Algorithm１の２９行目）。

ＫＫＴ条件判定部１１０は、予測変数の集合に含まれる予測変数がＫＫＴ条件を満たすか否かを判定し、予測変数の全てがＫＫＴ条件を満たすまで係数の更新を繰り返し実行させる。

具体的に、ＫＫＴ条件判定部１１０は、第２係数更新部１０９における処理が完了した後、各係数に対応する予測変数のうち、ｐ_ｉ∈Ｓ_ｋであるような予測変数、およびｐ_ｉ∈Ｐである予測変数に対してそれぞれＫＫＴ条件を用いて係数が非零になる予測変数がないか判定を行う。そして、ＫＫＴ条件判定部１１０は、ＫＫＴ条件を満たさない予測変数をさらに更新対象の予測変数の集合に追加し、ＫＫＴ条件を満たさない予測変数がなくなるまでAlgorithm１の８行目からの処理を繰り返し実行させる（Algorithm１の３０〜３７行目）。

回帰分析出力部１１１は、ＫＫＴ条件判定部１１０によって、チューニングパラメータにおいて全ての予測変数がＫＫＴ条件を満たすと判定された場合、パラメータ設定部１０１に、チューニングパラメータにおける予測変数を基に次のチューニングパラメータを設定させる。具体的に、回帰分析出力部１１１は、チューニングパラメータλ_ｋについて、ＫＫＴ条件判定部１１０によって全ての予測変数がＫＫＴ条件を満たすと判定された場合、予測変数を出力する。そして、パラメータ設定部１０１は、回帰分析出力部１１１によって出力された予測変数を基に、チューニングパラメータλ_ｋを設定する。

そして、回帰分析装置１は、最後に、それぞれのチューニングパラメータに対して解を求めた後、事前に決めた基準に対して最適な解を設定する。なお、パラメータ設定部１０１はチューニングパラメータとして正の値を設定した場合、非零係数判定部１０５および零係数判定部１０８における条件判定は、次のように説明することができる。すなわち、非零係数判定部１０５は、下限値がチューニングパラメータより大きい、または上限値がチューニングパラメータを負にした値より小さい場合に、係数が更新により必ず非零になると判定する。また、零係数判定部１０８は、上限値がチューニングパラメータより大きい、または下限値がチューニングパラメータを負にした値より小さい場合に、係数が更新により非零になる可能性があると判定する。

第１の実施形態における手法の計算コストについては、図１１の定理１に示すことがいえる。図１１は定理１を示す図である。定理１ではｋ＝０のとき、ｍ_ｋ＝０とする。第１の実施形態における手法では、係数が０となる予測変数の更新を枝刈りできるため、従来の手法より高速な処理が可能となる。

また、第１の実施形態における手法の回帰の精度については、図１２の定理２に示すことがいえる。図１２は定理２を示す図である。ｐ＜ｎである場合、予測変数の行列は列フルランクにならないため１つの解に必ずしも収束しない。しかし、発明者は、第１の実施形態における手法によれば、従来の手法より少ない予測変数で高い精度の回帰を行えることを実験的に確かめている。

［第１の実施形態の処理］
図１３を用いて第１の実施形態の処理について説明する。図１３は、第１の実施形態に係る回帰分析処理の流れを示すフローチャートである。図１３に示すように、まず、回帰分析装置１は、分析対象となる予測変数の行列が入力されると（ステップＳ１０１）、ｋが１からＫに至るまで、以下の処理を実行する（ステップＳ１０２〜Ｓ１３９）。

まず、パラメータ設定部１０１は、λ＝λ_ｋとする（ステップＳ１０３）。ｋ＝１である場合（ステップＳ１０４、Ｙｅｓ）、初期更新対象設定部１０２は、Ｕを空集合とする（ステップＳ１０５）。ｋ＝１でない場合（ステップＳ１０４、Ｎｏ）、初期更新対象設定部１０２は、ＵをＰ_ｋ−１とする（ステップＳ１０６）。次に、初期値設定部１０３は、係数の初期値を計算する（ステップＳ１０７）。

そして、第１参照ベクトル設定部１０４は、参照ベクトルを設定する（ステップＳ１０８）。ここで、非零係数判定部１０５は、全てのｐ_ｉについて係数パラメータを計算する（ステップＳ１０９〜Ｓ１１１）。さらに、非零係数判定部１０５は、全てのｐ_ｉについて係数パラメータの上限値および下限値を計算し、上限値、下限値およびチューニングパラメータが所定の条件を満たすか否かを判定する（ステップＳ１１２、Ｓ１１３）ことにより、係数が必ず非零になるか否かを判定する。

所定の条件が満たされる場合（ステップＳ１１３、Ｙｅｓ）、第１係数更新部１０６は、係数およびパラメータを更新する（ステップＳ１１４、Ｓ１１５）。また、所定の条件が満たされない場合（ステップＳ１１３、Ｎｏ）、第１係数更新部１０６は、更新を行わない。そして、非零係数判定部１０５および第１係数更新部１０６は、係数ベクトルが収束するまでこれらの処理を繰り返す（ステップＳ１１２、Ｓ１１６、Ｓ１１７）。

そして、第２参照ベクトル設定部１０７は、更新された係数ベクトルを基に参照ベクトルを設定する（ステップＳ１１８）。ここで、零係数判定部１０８は、全てのｐ_ｉについて係数パラメータを計算する（ステップＳ１１９〜Ｓ１２１）。さらに、零係数判定部１０８は、全てのｐ_ｉについて係数パラメータの上限値および下限値を計算し、上限値、下限値およびチューニングパラメータが所定の条件を満たすか否かを判定する（ステップＳ１２２、Ｓ１２３）ことにより、係数が更新により非零になる可能性があるか否かを判定する。

所定の条件が満たされる場合（ステップＳ１２３、Ｙｅｓ）、第２係数更新部１０９は、係数を更新する（ステップＳ１２４）。また、所定の条件が満たされない場合（ステップＳ１２３、Ｎｏ）、第２係数更新部１０９は、係数を０に更新する（ステップＳ１２５）。そして、第２係数更新部１０９は、パラメータを更新する（ステップＳ１２６）。零係数判定部１０８および第２係数更新部１０９は、係数ベクトルが収束するまでこれらの処理を繰り返す（ステップＳ１２２、Ｓ１２７、Ｓ１２８）。

そして、ＫＫＴ条件判定部１１０は、Ｓ_ｋに含まれる全てのｐ_ｉについてＫＫＴ条件を満たすか否かを判定し、ＫＫＴ条件を満たさないｐ_ｉをＵに追加する（ステップＳ１２９〜Ｓ１３２）。そして、ＫＫＴ条件を満たさないｐ_ｉが存在する場合（ステップＳ１３３、Ｎｏ）は、ステップＳ１０８以降の処理をさらに実行させる。

また、Ｓ_ｋに含まれる全てのｐ_ｉがＫＫＴ条件を満たす場合（ステップＳ１３３、Ｙｅｓ）、ＫＫＴ条件判定部１１０は、Ｐに含まれる全てのｐ_ｉについてＫＫＴ条件を満たすか否かを判定し、ＫＫＴ条件を満たさないｐ_ｉをＵに追加する（ステップＳ１３４〜Ｓ１３７）。そして、ＫＫＴ条件を満たさないｐ_ｉが存在する場合（ステップＳ１３８、Ｎｏ）は、ステップＳ１０８以降の処理をさらに実行させる。そして、全てのｋについて処理が完了すると、線形回帰の結果が得られる（ステップＳ１４０）。

［第１の実施形態の効果］
lassoを高速に処理するため多くの研究者はscreening techniqueという手法の研究に注力している。screening techniqueは、係数の繰り返し更新に入る前に係数の値が０になる予測係数を特定し、coordinate descentに基づく手法を高速にする手法であるが、大規模のデータを高速に扱うためにはcoordinate descentに基づく手法自体を高速にする必要がある。

そこで、第１の実施形態においては、パラメータ設定部１０１は、lassoにおけるチューニングパラメータを順次設定する。そして、初期更新対象設定部１０２は、パラメータ設定部１０１によってチューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する。また、初期値設定部１０３は、パラメータ設定部１０１によってチューニングパラメータが設定されるたびに、係数の初期値を設定する。

第１参照ベクトル設定部１０４は、係数を基に参照ベクトルを設定する。次に、非零係数判定部１０５は、チューニングパラメータと、参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、係数のそれぞれが更新により必ず非零になるか否かを判定する。そして、第１係数更新部１０６は、係数のうち、非零係数判定部１０５によって必ず非零になると判定された係数を、収束するまで繰り返し更新する。

第２参照ベクトル設定部１０７は、第１係数更新部１０６によって更新された係数を基に参照ベクトルをさらに設定する。次に、零係数判定部１０８は、チューニングパラメータと、参照ベクトルを基にさらに計算された上限値および下限値と、が所定の条件を満たすか否かによって、係数のそれぞれが更新により非零になる可能性があるか否かを判定する。そして、第２係数更新部１０９は、係数のうち、零係数判定部１０８によって非零になる可能性があると判定された係数を、収束するまで繰り返し更新する。

これにより、繰り返し計算において係数が０になる予測変数に対する更新を枝刈りできるため、coordinate descentに基づく手法自体を高速化することができる。また、第１の実施形態によれば、定理２等に示すように、従来の手法と同等以上の精度で線形回帰を行うことができる。また、第１の実施形態によれば、各パラメータが自動的に計算されるため、ユーザによるパラメータ調整の必要がない。

また、回帰分析出力部１１１は、ＫＫＴ条件判定部１１０によって、チューニングパラメータにおいて全ての予測変数がＫＫＴ条件を満たすと判定された場合、パラメータ設定部１０１に、チューニングパラメータにおける予測変数を基に次のチューニングパラメータを設定させる。これにより、効率的にチューニングパラメータの設定を行うことができる。

また、例えば、パラメータ設定部１０１はチューニングパラメータとして正の値を設定する。このとき、非零係数判定部１０５は、下限値がチューニングパラメータより大きい、または上限値がチューニングパラメータを負にした値より小さい場合に、係数が更新により必ず非零になると判定する。また、零係数判定部１０８は、上限値がチューニングパラメータより大きい、または下限値がチューニングパラメータを負にした値より小さい場合に、係数が更新により非零になる可能性があると判定する。このように、条件を設定し判定を行うことで、更新対象となる係数の枝刈りを効率良く行うことができる。

［システム構成等］
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。さらに、各装置にて行なわれる各処理機能は、その全部または任意の一部が、ＣＰＵ（Central Processing Unit）および当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。

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

［プログラム］
一実施形態として、回帰分析装置は、パッケージソフトウェアやオンラインソフトウェアとして上記の回帰分析を実行する回帰分析プログラムを所望のコンピュータにインストールさせることによって実装できる。例えば、上記の回帰分析プログラムを情報処理装置に実行させることにより、情報処理装置を回帰分析装置として機能させることができる。ここで言う情報処理装置には、デスクトップ型またはノート型のパーソナルコンピュータが含まれる。また、その他にも、情報処理装置にはスマートフォン、携帯電話機やＰＨＳ（Personal Handyphone System）等の移動体通信端末、さらには、ＰＤＡ（Personal Digital Assistant）等のスレート端末等がその範疇に含まれる。

また、回帰分析装置は、ユーザが使用する端末装置をクライアントとし、当該クライアントに上記の回帰分析に関するサービスを提供するサーバ装置として実装することもできる。例えば、回帰分析装置は、予測変数および応答変数を入力とし、回帰式を出力とする回帰分析サービスを提供するサーバ装置として実装される。この場合、回帰分析装置は、Ｗｅｂサーバとして実装することとしてもよいし、アウトソーシングによって上記の回帰分析に関するサービスを提供するクラウドとして実装することとしてもかまわない。

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

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

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

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

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

１回帰分析装置
１０１パラメータ設定部
１０２初期更新対象設定部
１０３初期値設定部
１０４第１参照ベクトル設定部
１０５非零係数判定部
１０６第１係数更新部
１０７第２参照ベクトル設定部
１０８零係数判定部
１０９第２係数更新部
１１０ＫＫＴ条件判定部
１１１回帰分析出力部

Claims

lassoを用いて、予測変数の持つ係数を繰り返し更新し回帰分析を行う回帰分析装置であって、
前記lassoにおけるチューニングパラメータを順次設定するパラメータ設定部と、
前記パラメータ設定部によって前記チューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する初期更新対象設定部と、
前記パラメータ設定部によって前記チューニングパラメータが設定されるたびに、前記係数の初期値を設定する初期値設定部と、
前記係数を基に参照ベクトルを設定する第１参照ベクトル設定部と、
前記チューニングパラメータと、前記参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により必ず非零になるか否かを判定する非零係数判定部と、
前記係数のうち、前記非零係数判定部によって必ず非零になると判定された係数を、収束するまで繰り返し更新する第１係数更新部と、
前記第１係数更新部によって更新された前記係数を基に参照ベクトルをさらに設定する第２参照ベクトル設定部と、
前記チューニングパラメータと、前記参照ベクトルを基にさらに計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により非零になる可能性があるか否かを判定する零係数判定部と、
前記係数のうち、前記零係数判定部によって非零になる可能性があると判定された係数を、収束するまで繰り返し更新する第２係数更新部と、
前記予測変数の集合に含まれる予測変数がＫＫＴ条件を満たすか否かを判定し、前記予測変数の全てがＫＫＴ条件を満たすまで係数の更新を繰り返し実行させるＫＫＴ条件判定部と、
を有することを特徴とする回帰分析装置。
前記ＫＫＴ条件判定部によって、前記チューニングパラメータにおける全ての予測変数がＫＫＴ条件を満たすと判定された場合、前記パラメータ設定部に、前記チューニングパラメータにおける前記予測変数を基に次のチューニングパラメータを設定させる回帰分析出力部をさらに有することを特徴とする請求項１に記載の回帰分析装置。
前記パラメータ設定部は前記チューニングパラメータとして正の値を設定し、
前記非零係数判定部は、前記下限値が前記チューニングパラメータより大きい、または前記上限値が前記チューニングパラメータを負にした値より小さい場合に、前記係数が更新により必ず非零になると判定することを特徴とする請求項１または２に記載の回帰分析装置。
前記パラメータ設定部は前記チューニングパラメータとして正の値を設定し、
前記零係数判定部は、前記上限値が前記チューニングパラメータより大きい、または前記下限値が前記チューニングパラメータを負にした値より小さい場合に、前記係数が更新により非零になる可能性があると判定することを特徴とする請求項１から３のいずれか１項に記載の回帰分析装置。
前記非零係数判定部は、前記係数のうちのi番目の係数について、上限値を、式（１）により計算し、
前記零係数判定部は、前記係数のうちのi番目の係数について、下限値を、式（２）により計算することを特徴とする請求項１から４のいずれか１項に記載の回帰分析装置。
回帰分析装置が、lassoを用いて、予測変数の持つ係数を繰り返し更新し回帰分析を行う回帰分析方法であって、
前記lassoにおけるチューニングパラメータを順次設定するパラメータ設定工程と、
前記パラメータ設定工程によって前記チューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する初期更新対象設定工程と、
前記パラメータ設定工程によって前記チューニングパラメータが設定されるたびに、前記係数の初期値を設定する初期値設定工程と、
前記係数を基に参照ベクトルを設定する第１参照ベクトル設定工程と、
前記チューニングパラメータと、前記参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により必ず非零になるか否かを判定する非零係数判定工程と、
前記係数のうち、前記非零係数判定工程によって必ず非零になると判定された係数を、収束するまで繰り返し更新する第１係数更新工程と、
前記第１係数更新工程によって更新された前記係数を基に参照ベクトルをさらに設定する第２参照ベクトル設定工程と、
前記チューニングパラメータと、前記参照ベクトルを基にさらに計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により非零になる可能性があるか否かを判定する零係数判定工程と、
前記係数のうち、前記零係数判定工程によって非零になる可能性があると判定された係数を、収束するまで繰り返し更新する第２係数更新工程と、
前記予測変数の集合に含まれる予測変数がＫＫＴ条件を満たすか否かを判定し、前記予測変数の全てがＫＫＴ条件を満たすまで係数の更新を繰り返し実行させるＫＫＴ条件判定工程と、
を含んだことを特徴とする回帰分析方法。
コンピュータを、請求項１から５のいずれか１項に記載の回帰分析装置として機能させるための回帰分析プログラム。