(57)【特許請求の範囲】
【請求項１】
lassoを用いて、予測変数の持つ係数を繰り返し更新し回帰分析を行う回帰分析装置であって、
前記lassoにおけるチューニングパラメータを順次設定するパラメータ設定部と、
前記パラメータ設定部によって前記チューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する初期更新対象設定部と、
前記パラメータ設定部によって前記チューニングパラメータが設定されるたびに、前記係数の初期値を設定する初期値設定部と、
前記係数を基に参照ベクトルを設定する第１参照ベクトル設定部と、
前記チューニングパラメータと、前記参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により必ず非零になるか否かを判定する非零係数判定部と、
前記係数のうち、前記非零係数判定部によって必ず非零になると判定された係数を、収束するまで繰り返し更新する第１係数更新部と、
前記第１係数更新部によって更新された前記係数を基に参照ベクトルをさらに設定する第２参照ベクトル設定部と、
前記チューニングパラメータと、前記参照ベクトルを基にさらに計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により非零になる可能性があるか否かを判定する零係数判定部と、
前記係数のうち、前記零係数判定部によって非零になる可能性があると判定された係数を、収束するまで繰り返し更新する第２係数更新部と、
前記予測変数の集合に含まれる予測変数がＫＫＴ条件を満たすか否かを判定し、前記予測変数の全てがＫＫＴ条件を満たすまで係数の更新を繰り返し実行させるＫＫＴ条件判定部と、
を有することを特徴とする回帰分析装置。
【請求項２】
前記ＫＫＴ条件判定部によって、前記チューニングパラメータにおける全ての予測変数がＫＫＴ条件を満たすと判定された場合、前記パラメータ設定部に、前記チューニングパラメータにおける前記予測変数を基に次のチューニングパラメータを設定させる回帰分析出力部をさらに有することを特徴とする請求項１に記載の回帰分析装置。
【請求項３】
前記パラメータ設定部は前記チューニングパラメータとして正の値を設定し、
前記非零係数判定部は、前記下限値が前記チューニングパラメータより大きい、または前記上限値が前記チューニングパラメータを負にした値より小さい場合に、前記係数が更新により必ず非零になると判定することを特徴とする請求項１または２に記載の回帰分析装置。
【請求項４】
前記パラメータ設定部は前記チューニングパラメータとして正の値を設定し、
前記零係数判定部は、前記上限値が前記チューニングパラメータより大きい、または前記下限値が前記チューニングパラメータを負にした値より小さい場合に、前記係数が更新により非零になる可能性があると判定することを特徴とする請求項１から３のいずれか１項に記載の回帰分析装置。
【請求項５】
前記非零係数判定部は、前記係数のうちのi番目の係数について、上限値を、式（１）により計算し、
前記零係数判定部は、前記係数のうちのi番目の係数について、下限値を、式（２）により計算することを特徴とする請求項１から４のいずれか１項に記載の回帰分析装置。
【数１】
0006453785.tif 000015
【請求項６】
回帰分析装置が、lassoを用いて、予測変数の持つ係数を繰り返し更新し回帰分析を行う回帰分析方法であって、
前記lassoにおけるチューニングパラメータを順次設定するパラメータ設定工程と、
前記パラメータ設定工程によって前記チューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する初期更新対象設定工程と、
前記パラメータ設定工程によって前記チューニングパラメータが設定されるたびに、前記係数の初期値を設定する初期値設定工程と、
前記係数を基に参照ベクトルを設定する第１参照ベクトル設定工程と、
前記チューニングパラメータと、前記参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により必ず非零になるか否かを判定する非零係数判定工程と、
前記係数のうち、前記非零係数判定工程によって必ず非零になると判定された係数を、収束するまで繰り返し更新する第１係数更新工程と、
前記第１係数更新工程によって更新された前記係数を基に参照ベクトルをさらに設定する第２参照ベクトル設定工程と、
前記チューニングパラメータと、前記参照ベクトルを基にさらに計算された上限値および下限値と、が所定の条件を満たすか否かによって、前記係数のそれぞれが更新により非零になる可能性があるか否かを判定する零係数判定工程と、
前記係数のうち、前記零係数判定工程によって非零になる可能性があると判定された係数を、収束するまで繰り返し更新する第２係数更新工程と、
前記予測変数の集合に含まれる予測変数がＫＫＴ条件を満たすか否かを判定し、前記予測変数の全てがＫＫＴ条件を満たすまで係数の更新を繰り返し実行させるＫＫＴ条件判定工程と、
を含んだことを特徴とする回帰分析方法。
【請求項７】
コンピュータを、請求項１から５のいずれか１項に記載の回帰分析装置として機能させるための回帰分析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、回帰分析装置、回帰分析方法および回帰分析プログラムに関する。
【背景技術】
【０００２】
高次元データに対するｌ１正則化付き最小２乗法の代表的な手法として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を用いた回帰分析を行う従来の手法について具体的に説明する。回帰分析において、予測変数の数はｐとし、観測値の数はｎとする。また、全てのベクトルは平均を０、分散を１に正規化されているものとする。また、応答変数は、式（１）のベクトルで表される。
【００１６】
【数１】
0006453785.tif 000002
【００１７】
また、予測変数は、式（２）の行列で表される。
【数２】
0006453785.tif 000003
【００１８】
lassoにおいては、式（２）の予測変数の行列における最小２乗誤差とｌ１正則化項の制限からなる以下の式を最小化する応答変数を予測するモデルを求める。
【００１９】
【数３】
0006453785.tif 000004
【００２０】
このとき、式（２）の予測変数の行列が列フルランク行列であれば、上記の最適化問題の解はただ１つになる。そうでない場合、解は１つに定まらない。なおｐ＞ｎである場合は、式（２）の予測変数の行列は、列フルランク行列にならない。
【００２１】
また、実際のアプリケーションにおいて最適なλの値は異なる。そのため、Ｋをチューニングパラメータの数としたときに、式（３）においてチューニングパラメータをλ１＞λ２＞...λＫと変化させる。そして、チューニングパラメータを変化させつつ最適化問題を解き、複数のチューニングパラメータを試した後に、事前に決定した赤池情報量等の基準を用いて最終的なlassoにおける解を決定する。
【００２２】
Tibshiraniらは、coordinate  descentを用いて予測変数を１つずつ更新する手法を提案した。coordinate  descentは他の係数はすでに更新したことを仮定し、係数ｗ［ｉ］に対して部分的に最適化を行う手法である。coordinate  descentでは、まず式（４）が計算される。
【００２３】
【数４】
0006453785.tif 000005
【００２４】
さらに、coordinate  descentにおいては、式（５）により係数が更新される。
【００２５】
【数５】
0006453785.tif 000006
【００２６】
式（５）のＳ（ｚ［ｉ］，λ）はsoft-thresholding  operatorと呼ばれる処理であり、またｚ［ｉ］はｉ番目の予測変数に対応するパラメータとして式（６）により計算される。
【００２７】
【数６】
0006453785.tif 000007
【００２８】
coordinate  descentにおいては、係数が総当たり的に更新される。すなわち、式（５）を用いて収束するまで全ての予測変数の係数が繰り返し更新される。そこで、Tibshiraniらは、更新において計算されるパラメータｚ［ｉ］を式（７）のように計算することで高速に求める手法を提案した。
【００２９】
【数７】
0006453785.tif 000008
【００３０】
式（６）と式（７）は同じ計算結果となる。また、ｍを非零の係数を持つ予測変数の数とすると、係数を更新するのに必要な計算コストは、式（７）ではＯ（ｍ）であり、式（６）ではＯ（ｎ）となる。lassoでは、疎に予測変数が選択されるため、ｍ≪ｎとなる。そのため、式（７）を式（６）の代わりに用いることによって、高速に係数を更新することができる。
【００３１】
式（７）を用いるためには、繰り返し計算に入る前に応答変数と全ての予測変数との内積を計算する必要がある。さらに、ある予測変数が応答変数を予測するのに用いられると判断された場合、その予測変数とその他全ての応答変数との内積を計算する必要がある。すなわち、ある予測変数が解として非零の係数を持つ場合、式（７）を利用するために、
その予測変数における内積を計算しなければならない。
【００３２】
そこで、lassoを高速に計算するために、Tibshiraniらを含む多くの研究者によって、係数が０になる予測変数を特定するscreening  techniqueが提案されている。特にsequential  strong  ruleにおいては、式（８）を用いた以下の方法によって予測変数の枝刈りが行われる。
【００３３】
【数８】
0006453785.tif 000009
【００３４】
sequential  strong  ruleは、誤って係数が０にならない予測誤差を枝刈りすることがあるため、収束後に全ての予測変数に対してＫＫＴ（Karush-Kuhn-Tucker）条件が満たされているかを確認する。ＫＫＴ条件は、式（９）を用いて以下のように表される。なお、ＫＫＴ条件は、式（７)を用いることにより、Ｏ（ｍ）のコストで計算できる。
【００３５】
【数９】
0006453785.tif 000010
【００３６】
［第１の実施形態の概要］
従来の手法では、screening  techniqueで枝刈りされなかった全ての予測変数に対して収束するまで繰り返し係数の更新が行われる。このとき、従来の手法では、ある予測変数の係数の更新を行うのにその他の予測変数の係数で非零であるものを使う。そのため、もし繰り返し計算の中である予測変数が非零の係数を持てば、その係数が繰り返し計算の後０に収束したとしても、計算コストの増加につながってしまう。
【００３７】
提案手法、すなわち第１の実施形態において用いられる手法も、従来の手法と同様にcoordinate  descentに基づく手法である。しかし、第１の実施形態においては、従来の手法と異なり、選択的に係数の更新を行う。すなわち、第１の実施形態における手法は、まず非零の係数を持つ予測変数のみを収束するまで更新し、その後に非零の係数を持つことが期待される予測変数を収束するまで更新する。この結果、第１の実施形態における手法は、繰り返し計算において係数が０になる予測変数に対する更新を枝刈りできるため、高速にlassoにおける解を計算することができる。
【００３８】
このように、screening  techniqueが繰り返し計算に入る前に係数が０である予測変数を枝刈りするのに対し、第１の実施形態においては、繰り返し計算において係数が０である予測変数を枝刈りする。ｐを予測変数の数、ｎを観測値の数としたとき、与えられるデータはｎ×ｐの行列で表現することができる。もし、当該行列が列フルランクであれば、lassoにおける解は凸関数性から１つになるため、第１の実施形態では、従来の手法と同じ結果を出すことが保証される。一方、当該行列が列フルランクでない場合、第１の実施形態では勾配が０である予測変数の更新を行わないため、従来の手法より少ない予測変数で、より精度の高い予測を行うことができる。
【００３９】
［第１の実施形態の構成］
まず、図１を用いて、第１の実施形態に係る回帰分析装置の構成について説明する。図１は、第１の実施形態に係る回帰分析装置の構成の一例を示すブロック図である。図１に示すように、回帰分析装置１は、パラメータ設定部１０１、初期更新対象設定部１０２、初期値設定部１０３、第１参照ベクトル設定部１０４、非零係数判定部１０５、第１係数更新部１０６、第２参照ベクトル設定部１０７、零係数判定部１０８、第２係数更新部１０９、ＫＫＴ条件判定部１１０、および回帰分析出力部１１１を有する。
【００４０】
図２を用いて、回帰分析装置１の各部の処理について説明する。図２は、第１の実施形態に係る回帰分析処理のアルゴリズムを示す図である。図２に示すAlgorithm１において、チューニングパラメータの設定（例えばAlgorithm１の２行目）、式（７）を用いた係数の更新、ＫＫＴ条件の判定（例えばAlgorithm１の３１行目）等については、前述の従来の手法と同様の処理を適用することができる。一方、係数から計算される所定のパラメータの上限値および下限値を用いて予測変数の枝刈りを行う処理（例えばAlgorithm１の１５、２４行目）は、従来の手法には含まれない処理である。
【００４１】
図２に示すように、Algorithm１においては、Ｕを繰り返し計算において更新を行う予測変数の集合、Ｐを全ての予測変数の集合、Ｐｋをチューニングパラメータλｋへの解として非零の係数となる予測変数の集合、Ｓｋをチューニングパラメータλｋにおけるsequential  strong  ruleで枝刈りされなかった予測変数の集合とする。
【００４２】
回帰分析装置１は、ｋを１，２，…，Ｋと設定し、以下の処理を行う（Algorithm１の１行目）。まず、パラメータ設定部１０１は、lassoにおけるチューニングパラメータを順次設定する。具体的に、パラメータ設定部１０１は、λ：＝λｋとする（Algorithm１の２行目）。
【００４３】
そして、初期更新対象設定部１０２は、パラメータ設定部１０１によってチューニングパラメータが設定されるたびに、更新対象である係数を持つ予測変数の集合を設定する。
具体的に、初期更新対象設定部１０２は、ｋ＝１であれば、Ｕを空集合とし、そうでなければＵ：＝Ｐｋ−１とする（Algorithm１の３〜６行目）。
【００４４】
また、初期値設定部１０３は、パラメータ設定部１０１によってチューニングパラメータが設定されるたびに、係数の初期値を設定する。具体的に、初期値設定部１０３は、チューニングパラメータλｋとλｋ−１に対する解は似ているという知見に基づき、チューニングパラメータλｋに対する係数の初期値を、式（１０）により決定する（Algorithm１の７行目）。
【００４５】
【数１０】
0006453785.tif 000011
【００４６】
次に、第１参照ベクトル設定部１０４は、係数の初期値を基に参照ベクトルを設定する。具体的に、第１参照ベクトル設定部１０４は、繰り返し計算に入る前に、参照ベクトルを設定する（Algorithm１の１０行目）。
【００４７】
そして、非零係数判定部１０５は、チューニングパラメータと、参照ベクトルを基に計算された上限値および下限値と、が所定の条件を満たすか否かによって、係数のそれぞれが更新により必ず非零になるか否かを判定する。
【００４８】
具体的に、非零係数判定部１０５は、上限値を図３に示す定義１のように定義する。図３は、定義１を示す図である。図３の式（Ａ）におけるベクトルおよびパラメータは以下のように説明することができる。
【００４９】
【数１１】
0006453785.tif 000012
【００５０】
また、係数ベクトルが各繰り返しで更新される長さｐのベクトルであること等から、定義１は、以下の性質を持つ。
【００５１】
【数１２】
0006453785.tif 000013
【００５２】
【数１３】
0006453785.tif 000014
【００５３】
また、非零係数判定部１０５は、下限値を図４に示す定義２のように定義する。図４は、定義２を示す図である。上限値および下限値が定義１および定義２で定義されることは、図５および図６に示す補題１および補題２によって示される。図５は、補題１を示す図である。また、図６は、補題２を示す図である。
【００５４】
式（１３）は、Ｏ（１）で更新できることは明らかであるため、上限値と下限値の計算においては、図７に示す補題３が成り立つ。図７は、補題３を示す図である。ここで、必ず非零の係数を持つ予測変数は、図８の補題４に示す性質を持つ。図８は、補題４を示す図である。また、必ず係数が０になる予測変数は、図９の補題５に示す性質を持つ。図９は、補題５を示す図である。また、補題５より、非零の係数を持つ可能性のある予測変数は、図１０の補題６に示す性質を持つ。図１０は、補題６を示す図である。
【００５５】
予測変数が補題４の条件を満たすとき、その予測変数は必ず補題６の条件を満たす。すなわち、補題４による予測変数の集合は、必ず補題６による予測変数の集合に含まれる。これは、もし下限値がチューニングパラメータより小さければ、上限値もチューニングパラメータより小さく、上限値がチューニングパラメータを負にした値より小さければ、下限値もチューニングパラメータを負にした値より小さいからである。
【００５６】
そのため、もし予測変数において補題４より係数が非零となるとき、その予測変数は補題６より非零の係数を持つ可能性があるとされる。さらに、予測変数が補題６の条件を満たさないとき、補題６よりその予測変数の係数は０となる。第１の実施形態ではこの性質を用いて更新対象となる係数の枝刈りを行っている。
【００５７】
非零係数判定部１０５は、定義１および定義２により、Ｕに含まれる予測変数について、係数パラメータｚ［ｉ］、ｚ［ｉ］の上限値および下限値を計算し、上限値、下限値およびチューニングパラメータを用いて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）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３およびプログラムデータ１０９４は、他のコンピュータから、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【符号の説明】
【００９４】
１  回帰分析装置
１０１  パラメータ設定部
１０２  初期更新対象設定部
１０３  初期値設定部
１０４  第１参照ベクトル設定部
１０５  非零係数判定部
１０６  第１係数更新部
１０７  第２参照ベクトル設定部
１０８  零係数判定部
１０９  第２係数更新部
１１０  ＫＫＴ条件判定部
１１１  回帰分析出力部
【図１】
0006453785.tif 000016
【図２】
0006453785.tif 000017
【図３】
0006453785.tif 000018
【図４】
0006453785.tif 000019
【図５】
0006453785.tif 000020
【図６】
0006453785.tif 000021
【図７】
0006453785.tif 000022
【図８】
0006453785.tif 000023
【図９】
0006453785.tif 000024
【図１０】
0006453785.tif 000025
【図１１】
0006453785.tif 000026
【図１２】
0006453785.tif 000027
【図１３】
0006453785.tif 000028
【図１４】
0006453785.tif 000029