確率・統計

(23) 正則化最小二乗法 ( Regularized Least Square Method )

最小二乗法は 19 世紀初めに公開されて以来、様々な場面で利用されている重要なデータ解析手法です。最も単純なものはデータの直線への当てはめであり、それから発展して様々な手法が利用されています。今回はその中から「正則化最小二乗法 ( Regularized Least Square Method )」を紹介したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) リッジ回帰 ( Ridge Regression )

線形重回帰モデルでは、独立変数に対して従属変数が以下のようなモデル式で表されると仮定していました。

E[y]=a0 + a1x1 + a2x2 + ... + apxp
=a0 + Σk{1→p}( akxk )
=( a, x ) = aTx

上式の左辺 E[y] は、従属変数 y の期待値を表します。ak ( k = 0, 1, ... p ) は各独立変数にかかる係数であり、E[y] に対して作用する重みを表すと考えることができます。また、a0 は定数項なので、右辺の大きさを従属変数に揃える役割を持っています。独立変数は p 個あり、それに定数項にかかる要素 1 を加えたベクトルが独立変数ベクトル x = ( 1, x1, x2, ... xp )T です。また、係数ベクトルは a = ( a0, a1, a2, ... ap )T で表され、モデル式はこの二つのベクトルの内積として表されます。「重回帰分析」「分散分析法」「共分散分析」は全て y が正規分布に従うとして最尤推定法により係数を求める手法でした。また、E[y] に連結関数を作用させ、さらに確率分布を正規分布以外のものまで利用するようにした手法が「一般化線形モデル」でした。

線形という言葉は独立変数ではなくその係数に対して用いられています。つまり、独立変数が非線形であっても係数に対して線形式になっていれば問題ありません。従って、線形重回帰モデルはもっと一般に

E[y]=a0 + a1f1(x) + a2f2(x) + ... + apfp(x)
=a0 + Σk{1→p}( akfk(x) )
=( a, f ) = aTf

という形で表すこともできます。f = ( 1, f1(x), f2(x), ... fp(x) ) は x を変数とする任意の関数からなるベクトルになります。例えば、あるデータを p 次の多項式で近似したいとすれば、f = ( 1, x, x2, ... xp ) として係数 a を解けばいいわけです。

y の従う確率分布に、分散が一定の正規分布 N( E[y], σ2 ) を仮定した場合、尤度を最大にするためには以下の正規方程式を解く必要があります。

XTXa = XTy

上式で、X は独立変数ベクトル x を行とする行列で、y は従属変数の観測値です。X の行数を N とした場合、y のサイズも当然 N でなければなりません。この式は、p + 1 個の未知数 ak を持った p + 1 個の線形式からなる連立方程式であり、それを解くことで係数を得ることができます。この式を変形すれば

a = (XTX)-1XTy

となって、右辺によって a の値が定まります。

下の表は、参考文献に従って、sin( 2πx ) の値に正規分布 N( 0, 0.12 ) に従う乱数をノイズとして加算して測定値を人工的に作成したデータです。

図表 1-1. 最尤推定用テストデータ
xsin(2πx)乱数 ey = sin(2πx) + eテストデータのグラフ
0.0000.000-0.054-0.054
0.1110.643-0.1480.495
0.2220.9850.0150.999
0.3330.8660.0160.882
0.4440.3420.0320.374
0.556-0.3420.073-0.269
0.667-0.866-0.041-0.907
0.778-0.9850.173-0.812
0.889-0.643-0.267-0.910
1.0000.000-0.041-0.041

右側にあるグラフは、曲線が実際に推定したい関数 ( sin 2πx ) で、プロットが乱数を加えた擬似測定値です。このデータに対して、一次(直線)から 9 次までの多項式として最小二乗法により係数を推定した結果が下表になります。

図表 1-2. 最小二乗法による係数の推定結果
多項式の次数推定した多項式のグラフ
1次2次3次4次5次6次7次8次9次
係数定数項6.52E-015.64E-01-1.55E-01-1.17E-01-7.35E-02-5.00E-02-5.67E-02-5.36E-02-5.40E-02
x-1.35E+00-7.58E-011.11E+019.64E+005.72E+00-5.55E-016.60E+00-1.52E+012.08E+01
x2-5.94E-01-3.18E+01-2.45E+018.17E+008.55E+01-3.26E+014.08E+02-4.23E+02
x32.08E+019.20E+00-8.33E+01-4.16E+022.84E+02-3.04E+034.45E+03
x45.79E+001.12E+027.55E+02-1.22E+031.13E+04-2.39E+04
x5-4.24E+01-6.15E+022.27E+03-2.38E+047.20E+04
x61.91E+02-1.90E+032.85E+04-1.28E+05
x75.99E+02-1.80E+041.34E+05
x84.66E+03-7.59E+04
x91.79E+04

推定した係数を見ると、次数を増やすに従って大きくなってゆく傾向にあることがわかります。また、右側のグラフは次数が 1, 3, 7, 9 の場合の多項式のグラフを示したもので、黒の点線で示したのが sin 2πx のグラフ、プロットは測定値をそれぞれ示しています。9 次の多項式のグラフを見ると、測定値のプロットを正確にたどっていることから誤差は完全にゼロとなっているにもかかわらず、得られた多項式は特に測定誤差(乱数の値)の絶対値が大きな場所で激しく振れています。最小二乗法は、一定の分散を持った正規分布に従った誤差をデータが持っていることを仮定して係数を推定する手法でした。推定したい係数に対してデータ数が充分にない場合、誤差成分は推定値に影響を及ぼすようになることをこの結果は示しています。9 次の係数での結果はその顕著な例であり、データ数と推定する係数の数が等しくなったため、誤差成分はゼロとなって、その分が係数に取り込まれた形になっています。このような現象を「過剰適合 (Overfitting)」といい、機械学習の分野では「過学習 (Overtraining)」とも呼ばれています。


直感的に考えると、曲線が大きく振れてしまうのは得られた推定値(係数)の変動が大きくなっているからです。そこで、係数が大きくなった場合はそれを抑制するような新たなパラメータを追加することを考えます。最小二乗法は、以下の二乗誤差を最小にすることを目的とした手法でした。

J = (1/2)Σi{1→N}( ( xiTa - yi )2 )

xi は独立変数ベクトルで yi は対応する従属変数、a が推定したい係数ベクトル (未知数) になります。a の係数である xi が小さくなると、得られる解が大きくなるので、それを抑制するよう左辺に (λ/2)||a||2 を加算して係数を補正します。

J = (1/2)Σi{1→N}( ( xiTa - yi )2 ) + (λ/2)||a||2

λ は第二項の影響度を調整するためのパラメータになります。この第二項は「ペナルティー項 (Penalty Term)」と呼ばれ、ペナルティー項を使った過剰適合への対応法を「正則化 (Regularization)」といいます。k 番目の係数 ak による偏微分 ∂J / ∂ak

∂J / ∂ak = Σi{1→N}( ( xiTa - yi )xi,k ) + λak

であり、これがゼロになるときは

Σi{1→N}( xiTaxi,k ) + λak = Σi{1→N}( yixi,k )

が成り立ちます。左辺は

Σi{1→N}( xiTaxi,k ) + λak=Σi{1→N}( Σj{0→p}( xi,jajxi,k ) ) + λak
=Σj{0→p}( Σi{1→N}( xi,jxi,k )aj ) + λak
=Σj{0→p}( x'jTx'kaj ) + λak

となって、正規方程式は

|x'0Tx'0,x'1Tx'0,x'2Tx'0,...x'pTx'0||a0|+|λa0|=|x'0Ty|
|x'1Tx'0,x'1Tx'1,x'2Tx'1,...x'pTx'1||a1||λa1||x'1Ty|
|x'2Tx'0,x'2Tx'1,x'2Tx'2,...x'pTx'2||a2||λa2||x'2Ty|
|:::...:||:||:||:|
|x'pTx'0,x'pTx'1,x'pTx'2,...x'pTx'p||ap||λap||x'pTy|

で表されます。但し、x'k は独立変数ベクトルの k 番目の要素を N 個並べたベクトルになります。x'k を k 番目の行ベクトルとした (つまり xi を i 番目の列ベクトルとした ) N 行 p + 1 列の行列は「デザイン行列 (Design Matrix)」といい、これを X で表せば、正規方程式の第一項は XTXa になるので、通常の最小二乗法の左辺と一致します。第二項は λa = λEa ( E は単位行列 ) で表され、右辺は XTy なので、この正規方程式は最終的に

( XTX + λE )a = XTy

となります。重回帰分析用のプログラムがあれば、連立方程式の係数行列を求める時に対角成分に対して λ を加算するだけで対応することができるので、実装は非常に簡単に行うことができます。但し、独立変数ベクトルの要素の中に定数項を持たせている場合、それには λ は適用しないか、他の補正項を適用するのが一般的なようです。このように、係数の二乗和をペナルティー項として利用する回帰分析は「リッジ回帰 ( Ridge Regression )」と呼ばれています。


独立変数ベクトルの一番目の要素が 1 ならば、a0 が定数項を意味します。このとき、正規方程式の左辺の係数 XTX の一番目の対角成分はデータ数 N であり、一列目の非対角成分 ( XTX が対称行列であることから一行目の要素でもあります ) は独立変数ベクトルの各要素の和 Σi{1→N}( xi,j ) になります。従って、定数項に対してペナルティー項を加えない場合、正規方程式の第一式は

N・a0 + Σj{1→p}( Σi{1→N}( xi,j )aj ) = Σi{1→N}( yi )

となり、両辺を N で割ると

a0 + mxTa = my

となります。但し、mx = ( mx1, mx2, ... mxp ) は定数項を除く独立変数ベクトルの標本平均であり、my は従属変数の標本平均です。また、係数ベクトル a は定数項を除いた形で再定義されます。この結果から

a0 = my - mxTa

を正規方程式の二式め以降に代入すると、

Σi{1→N}( xi,j )( my - mxTa ) + Σk{1→p}( x'kTx'jak ) + λf( aj ) = x'jy

より左辺は

(左辺)=Σi{1→N}( xi,j )my - Σi{1→N}( xi,jk{1→p}( mxkak ) + Σk{1→p}( Σi{1→N}( xi,kxi,j )ak ) + λf( aj )
=Σi{1→N}( xi,jmy ) - Σk{1→p}( Σi{1→N}( xi,jmxkak ) ) + Σk{1→p}( Σi{1→N}( xi,kxi,jak ) ) + λf( aj )
=Σi{1→N}( xi,jmy ) + Σk{1→p}( Σi{1→N}( xi,kxi,jak - xi,jmxkak ) ) + λf( aj )
=Σi{1→N}( xi,jmy ) + Σk{1→p}( Σi{1→N}( xi,j( xi,k - mxk )ak ) ) + λf( aj )

となって、第一項を右辺に移行して整理することで式は

Σk{1→p}( Σi{1→N}( xi,j( xi,k - mxk )ak ) ) + λf( aj ) = Σi{1→N}( xi,j( yi - my ) )

となります。ここで、f( aj ) はペナルティー項を表し、リッジ回帰ならば f( aj ) = aj です。

正規方程式の第一式に a0 を再代入すると

N( my - mxTa ) + Σk{1→p}( Σi{1→N}( xi,kak ) ) = Σi{1→N}( yi )

より

Σk{1→p}( Σi{1→N}( ( xi,k - mxk )ak ) ) = Σi{1→N}( yi - my )

となるので、両辺に mxj を掛けて二式め以降から辺々引くと

Σk{1→p}( Σi{1→N}( ( xi,j - mxj )( xi,k - mxk )ak ) ) + λf( aj ) = Σi{1→N}( ( xi,j - mxj )( yi - my ) )

であり、正規方程式の式・未知数の数は定数項の分を除くことができます。独立変数ベクトルの(標本)共分散行列を V、 xi,j と yi との標本共分散からなるベクトルを sy = ( sy1, sy2, ... syp )、ペナルティー項からなるベクトルを f とすれば、正規方程式は

Va + ( λ / N )f = sy

になります。ペナルティー項 ( λ / N )f を除いた式は、「確率・統計 (13) 回帰分析法」の中で導いたものと全く同じなので、定数項に対してペナルティー項を使わなければ、共分散行列を利用した解法においても係数行列の対角成分を変化させるだけで対応できることになります (*1-1)。


リッジ回帰を行うためのサンプル・プログラムを以下に示します。

/*
  線形モデルのチェック

  x : 独立変数ベクトル
  y : 従属変数
  n : 標本数を取得する変数へのポインタ
  p : 独立変数ベクトルのサイズを取得する変数へのポインタ
*/
void CheckLinearModel( const vector< vector< double > >& x, const vector< double >& y, size_t* n, size_t* p )
{
  *p = x.size();

  if ( *p < 1 )
    throw domain_error( "The size of x must be greater than zero." );

  *n = x[0].size();
  if ( y.size() != *n )
    throw domain_error( "The size of sample data ( x, y ) must be the same size." );
  if ( *n <= *p )
    throw domain_error( "The size of data must be greater than the size of independent variable." );
  for ( size_t i = 1 ; i < *p ; ++i ) {
    if ( x[i].size() != *n )
      throw domain_error( "All the variable size of data x must be same." );
  }
}

/*
  リッジ回帰モデル

  x : 独立変数ベクトル
  y : 従属変数
  lambda : 補正の重みパラメータ
  a : 推定した係数を取得するための配列へのポインタ
*/
void RidgeRegression( const vector< vector< double > >& x, const vector< double >& y, double lambda,
                      vector< double >* a )
{
  assert( &x != 0 && &y != 0 && a != 0 );

  size_t n; // 標本数
  size_t p; // 独立変数ベクトルのサイズ
  CheckLinearModel( x, y, &n, &p );

  // x, y の平均
  vector< double > mx( p ); // x[i] の平均
  for ( size_t j = 0 ; j < p ; ++j )
    mx[j] = accumulate( x[j].begin(), x[j].end(), 0.0 ) / n;
  double my = accumulate( y.begin(), y.end(), 0.0 ) / n; // y の平均

  // 共分散行列の作成
  LinearEquationSystem< double > s( p );
  for ( size_t r = 0 ; r < p ; ++r ) {
    s[r][r] = sampleVariance( x[r] ) + lambda / n;
    s.ans( r ) = sampleCovariance( x[r], y );
    for ( size_t c = r + 1 ; c < p ; ++c )
      s[r][c] = s[c][r] = sampleCovariance( x[r], x[c] );
  }

  // 回帰係数の推定
  if ( ! MathLib::GaussianElimination( s ) )
    throw runtime_error( "Failed to determine an estimator of regression coefficient." );

  a->resize( p + 1 );
  for ( size_t i = 0 ; i < p ; ++i )
    (*a)[i + 1] = s.ans( i );
  a->front() = my;
  for ( size_t i = 0 ; i < p ; ++i )
    a->front() -= mx[i] * (*a)[i + 1];
}

処理の内容は、以前紹介した重回帰モデルの回帰係数推定プログラムとほとんど同じで、唯一の違いは係数行列の対角成分となる独立変数ベクトルの分散から重みパラメータの λ を加算している個所のみになります。また、重回帰モデルではクラスの形をとりましたが、今回は推定した係数を返すだけに用途を絞り込んだので関数の形としています。
引数の x は独立変数ベクトルからなる二次元配列 (行列) を STL(Standard Template Library) のコンテナ・クラス vector で表現していますが、x の要素となる各 vector はあるデータ種の全標本からなる配列とみなして処理しているため、x のサイズはデータ種の数を表し、各要素となる配列のサイズが標本数を表すことに注意して下さい。各データ種の集計を行うため、このような構成にした方がプログラムが書きやすいことが理由ですが、x[i][j] の形で要素にアクセスした場合は j 番目の標本の i 番目のデータ種を取得することになり、行列としてみた場合は通常 i が標本番号、j がデータ種の番号となるのが一般的なのでそれとは異なることになります。


サンプル・プログラムを利用して、先に示した sin( 2πx ) の多項式近似を次数 9 次に対していくつかの重みパラメータで行った結果を以下に示します。

図表 1-3. リッジ回帰による係数の推定結果
重みパラメータ推定した多項式のグラフ(Ridge)
ln λ=-∞ (λ=0)ln λ=-20ln λ=-10ln λ=0 (λ=1)
係数定数項-5.40E-02-5.75E-02-8.76E-023.83E-01
x2.07E+014.61E+008.16E+00-4.04E-01
x2-4.23E+021.07E+01-1.55E+01-4.30E-01
x34.45E+03-5.48E+01-6.15E+00-3.02E-01
x4-2.39E+044.53E+017.12E+00-1.67E-01
x57.20E+04-1.75E+028.72E+00-5.36E-02
x6-1.28E+053.48E+023.32E+003.57E-02
x71.34E+051.28E+02-2.26E+001.05E-01
x8-7.59E+04-6.32E+02-3.79E+001.60E-01
x91.79E+043.27E+023.92E-012.03E-01

重みパラメータ λ は表の右になるほど大きくなり、それに対して係数の振れ方は小さくなることがこの結果からわかります。グラフを見ると、特に ln λ = -10 の場合の近似曲線は真の曲線にかなり近い状態となっています。

リッジ回帰を利用する必要がある場面は、推定しなければならないパラメータ(係数)に対してサンプルが少ないようなときが挙げられます。今回は例としてそのような場面を想定した条件で処理したわけですが、わざわざサンプルと同程度の次数にせずに、次数を減らして近似しても充分に精度のよい結果は得られたので、パラメータを減らすかサンプルを増やせばリッジ回帰を利用しなくても済みます。しかし、実際にはたくさんのパラメータを推定したいのにサンプルはそれほど多く取得できないような場合もあり得るため、そのようなときには発揮します。


リッジ回帰には利用する上で大きな問題があり、それは重みパラメータをどの程度の値にすればよいのか試行錯誤しなければわからないという点です。それでも、経験則として重みパラメータがある程度既知の状態となるなら、リッジ回帰は処理も非常に単純なのでよい選択肢となります。

以前、ベイスの定理を利用して、事後確率が尤度と事前確率の積で表せることを示しました。式に表すと以下のようになります。

P( θ | x ) ∝ P( θ )・P( x | θ )

パラメータ θ となる確率を事前分布、その条件下で x となる確率(尤度)が P( x | θ ) であり、その積が、x となったときにパラメータが θ である確率 (事後確率) に比例するというものです。この時の事前分布は「共役事前分布 (Conjugate Prior Distribution)」と呼ばれ、この確率分布が正規分布なら事後確率も正規分布となるのでした (*1-2)。

今、独立変数ベクトル xi = ( xi1, xi2, ... xip ) と従属変数 yi との関係が線形式 yi = xiα + ε で表され、誤差 ε は正規分布 N( 0, σ2 ) に従うと仮定します。事前確率の平均は未知の状態なので、パラメータの期待値は仮に 0 であり、分散も全て一定で独立であると考えて正規分布 N( 0, τ2E ) ( E は単位行列 ) に従うとすれば、事後確率による尤度はこの二つの分布の積で計算できるので、対数尤度は

-ln P( μ | y ) = Σi{1→N}( ( yi - xiα )2 ) / 2σ2 + αTα / 2τ2 + (定数)

となります。σ2 / τ2 = λ とすれば、これはリッジ回帰の式と一致することから、リッジ回帰はベイス推定における事後確率を使った最尤推定法と考えることができます。当然、各分布のパラメータは未知であり、パラメータ λ を決めることはこの式からもできません ( σ はデータから推定することは可能ですが、τ については推定を複数回行わない限り不可能です )。また、このことは当てはめにベイズ推定を用いることが非常に有効であることも表しています。このような手法は「最大事後確率推定 (Maximum a Posteriori Estimation)」またはその略称より「MAP 推定」と呼ばれます。


*1-1) 共分散行列を利用した式において、最初はペナルティー項を λf としていました。しかし、元の式から変形した場合は ( λ / N )f が正しいことになるので修正しました。どちらも影響度が変わるだけで処理自体は正常にできますが、N で割った場合は標本数が増えるに従ってペナルティー項の影響度が小さくなります。標本数が増えれば精度は上がるので影響が小さくなることをよしとするか、標本数に関係なく λ の値に完全に依存させるかについては選択の余地があると思います。

*1-2) ベイズ推定については「確率・統計 (17) ベイズの定理(Bayes' Theorem)」にて紹介しています。


1) Least Absolute Shrinkage and Selection Operator ( LASSO )

リッジ回帰では、ペナルティー項が係数の二乗和で表されていましたが、これとは異なるペナルティー項を検討することも可能です。例えば、リッジ回帰での係数の二乗和を一般化して

J = (1/2)Σi{1→N}( ( xiTa - yi )2 ) + (λ/q)Σj{0→p}( |aj|q )

を利用することが考えられ、特に q = 1 のときは「Lasso (Least Absolute Shrinkage and Selection Operator)」という名が付けられています。なお、係数の λ / q は様々な定義があり、分母の q がなかったり、リッジ回帰と同じ 2 が使われていたりしていました。これはそれほど重要なものではなく、q をつければ単に偏微分するときに分母がなくなるというだけなので、あまり気にすることはないと思います。また、リッジ回帰と Lasso のハイブリット型で「Elastic Net」というものもあり、その定義式は以下のようになります。

J = (1/2)Σi{1→N}( ( xiTa - yi )2 ) + λΣj{0→p}( (p/2)aj2 + ( 1 - p )|aj| )

p は 0 から 1 までの値をとる任意の定数で、p = 0 なら Lasso、p = 1 ならリッジ回帰そのものになります。

最小二乗法では、誤差の二乗和が推定する係数を変数とした楕円体になることを利用して、その極値を求める手法であるとも考えることができます。最も単純な例として、直線への当てはめ

J = (1/2)Σi{1→N}( ( axi + b - yi )2 )

を変形すると

J=(1/2)[ Σi{1→N}( a2xi2 + b2 + yi2 + 2abxi - 2byi - 2axiyi ) ]
=(1/2)[ ( a2Sxx - 2aSxy ) + 2abSx + ( Nb2 - 2bSy ) + Syy ]
=(1/2)[ Sxx( a - Sxy / Sxx )2 - Sxy / Sxx + 2abSx + N( b - Sy / N )2 - Sy2 / N + Syy ]
=(1/2)[ Sxx( a - Sxy / Sxx )2 + 2abSx + N( b - Sy / N )2 - Sxy / Sxx - Sy2 / N + Syy ]

となります。但し、Sxx, Sxy, Syy, Sx, Sy はそれぞれ xi2, xiyi, yi2, xi, yi の和を表します。これは 3 次元上の楕円体を表し、特に ab 項の係数 2Sx = 0 なら極値は ( a, b ) = ( Sxy / Sxx, Sy / N ) の時と簡単に求めることができます ( これは、回帰分析において x の標本平均がゼロの場合の値に一致します )。

ペナルティー項を加えた場合、楕円体はその分だけ値が大きくなりますが、ペナルティー項は原点を極値とする曲面 ( 特にリッジ回帰なら楕円体 ) を描くので、例えばペナルティー項による曲面が誤差項を完全に覆ってしまったとすれば、加算した結果の極値は原点側にシフトして係数は小さくなります。LASSO の場合は原点を頂点とした錐体になるため、重みパラメータ λ が大きくなるとゼロに収束してしまう現象が見られます。以下に、その様子を示したグラフを示します。

ペナルティー項による誤差関数の推移
図 2-1. ペナルティー項による誤差関数の推移

このグラフは一次元上の係数(横軸)に対する放物線 y = 3( x - 0.5 )2 + 0.5 にペナルティー項 y = λ|x| を λ = 0.1, 0.5, 1, 2, 5 の場合において加算した時の結果を表しています。λ が大きくなるに従い、極値は原点側にシフトし、特に λ = 5 のときは x = 0 にて最小値となっている様子が確認できます。


誤差関数部分を f(a)、ペナルティー項を g(a) として、g(a) = 0 の制約のもとで f(a) を最小にする問題を考えます。a の次元が p のとき、g(a) = 0 は p 次元の空間上にある超平面になります。例えば g(a) = K1a12 + K2a22 - r2 ( K1 > 0, K2 > 0 なら楕円体 ) のとき g(a) = 0 は空間上の楕円 K1a1 + K2a2 = r2 を意味します。f(a) が極値をとる点では、f(a) = (定数) と g(a) = 0 は接しています。なぜなら、もし接していなければ、g(a) = 0 上の点を f(a) の値が変化するように移動させることができるからです。つまり、f(a) と g(a) の各変数の微分係数からなるベクトルは平行になり、このベクトルをそれぞれ ∇f, ∇g とすれば

∇f - λ∇g = 0

が成り立ちます。この微分係数からなるベクトルは「勾配 (Gradient)」と呼ばれ、g(a) = 0 に対してどちらも垂直になります。上式に制約条件 g(a) = 0 を含めた p + 1 個の式を使い、p 個の変数 a と λ を求める手法は「ラグランジュの未定乗数法(Lagrange Multiplier)」といいます (*2-1)。このことから、正則化最小二乗法はラグランジュの未定乗数法を使って係数を求める操作と同じであることもわかります。但し、λ の値は最初から与えられ、それに従って制約条件の右辺の値も変化します。制約条件がない場合、J が最小となるのは明らかに λ = 0 のときです。これでは意味がないので、制約条件を付加してその中で最小値を決めるというのが正則化最小二乗法の考え方であるともいえます。

正則化最小二乗法の与える解の様子
図 2-2. 正則化最小二乗法の与える解の様子

リッジ回帰の場合とは異なり、LASSO では誤差関数との接し方に特徴があります。もし直線上で接すればリッジ回帰とよく似た結果になりますが、上図のようにちょうど頂点のところで接した場合は傾きの値にゼロが含まれることになります。これは、先に示した一次元での例において、λ = 5 の場合に発生したのと同じ現象です。このような解は「疎な解 (Sparse Solutions)」と呼ばれ、無関係な項に対する係数を完全に取り除くという効果があります。一次元での例で示したように、この効果は λ の値が大きいほど顕著になります。


ここで問題となるのが係数の求め方です。ペナルティー項は微分不可な個所があり、リッジ回帰のように解析的に求めることができません。そのため、特別なアルゴリズムが必要になります。高速なアルゴリズムとして知られているものに「Coordinate Descend Algorithm ; CDA」というものがあります。

最小二乗法では、全ての係数を連立方程式を使って解いていました。しかし、LASSO の場合は絶対値があるため、各係数に対して符号が正負のどちらになるかをあらかじめ決めておくことができません。組み合わせは 2 のべき乗となり、しかも係数の符号が変化すれば他の係数にも影響を及ぼすので、全ての組み合わせに対して計算を行うのは不可能です。そこで、一つの係数を除き他の係数は固定として計算を行い、全係数が収束するまで順番に処理を繰り返すという手法をとります。こうすると、個々の係数を計算するときは連立方程式を計算する必要がなく、処理を高速に行うこともできます。これが CDA の原理になります。

LASSO での正規方程式は次のようになります。

Σi{1→N}( ( xiTa - yi )xi,j ) + sign( aj )λ = 0

但し、sign( x ) は x の符号を表し、x ≥ 0 なら 1、x < 0 なら -1 となります。aj 以外を既知の数とした場合、

Σi{1→N}( ( xi,-jTa-j + xi,jaj - yi )xi,j ) + sign( aj )λ = 0 より

aj = -[ Σi{1→N}( ( xi,-jTa-j - yi )xi,j ) + sign( aj )λ ] / Σi{1→N}( xi,j2 )

で解を得ることができます。但し、xi,-ja-j はそれぞれ xia から j 番目の要素を除いたベクトルとします。右辺に aj の符号があるため 2 種類の解が得られることになり、左辺の aj が、sign( aj ) = 1 のときに非負、あるいは sign( aj ) = -1 のときに負となれば結果は正しいことになります。分母は必ず非負となることから、もし sign( aj ) = 1 のときに非負であれば、sign( aj ) = -1 のときも必ず非負であり、sign( aj ) = -1 のときに負であれば、sign( aj ) = 1 のときも必ず負数となるので、このときはどちらも前者が正しく後者が誤りとなります。もう一つの可能性は、sign( aj ) = 1 のときに負で sign( aj ) = -1 のときに正となる場合で、このときはどちらも正しくないことになります。これは | Σi{1→N}( ( xiTa - yi )xi,j ) | < λ のときに発生します。

r 回目の反復処理で得られる解を aj(r) としたとき、1 ≤ k < j の範囲の ak はすでに r + 1 回目の処理が完了しているので、上式は

aj(r+1) = -[ Σi{1→N}( ( Σk{0→j-1}( xi,kak(r+1) ) + Σk{j+1→p}( xi,kak(r) ) - yi )xi,j ) + sign( aj(r) )λ ] / Σi{1→N}( xi,j2 )

と変形することができて、さらに

Σi{1→N}( Σk{0→j-1}( xi,kak(r) )xi,j )
=Σk{0→j-1}( ak(r+1)Σi{1→N}( xi,kxi,j ) )
=Σk{0→j-1}( ak(r+1)Sj,k )

より

aj(r+1) = [ Sy,j - Σk{0→j-1}( ak(r+1)Sj,k ) - Σk{j+1→p}( ak(r)Sj,k ) - sign( aj(r) )λ ] / Sj,j

という結果が得られます。ここで、Sj,k は xi,j と xi,k の、Sy,j は yi と xi,j の積和をそれぞれ表します。よって、直前に求めた解 ak と Sj,k の j 以外の積和を求め、Sy,j - sign( aj(r) )λ から減算して Sj,j で割る操作を繰り返すことで、CDA を実現することができます。これは、連立方程式を「ガウス=ザイデル法 (Gauss-Seidel Method)」で求める操作と同等です。

S = Sy,j - Σk{0→j-1}( ak(r+1)Sj,k ) - Σk{j+1→p}( ak(r)Sj,k )

としたとき、S ≥ λ なら aj(r) の正負にかかわらず aj(r+1) は正なので、sign( aj(r) ) = 1 として結果を求めればよいことになります。逆に S ≤ -λ なら aj(r+1) は必ず負なので、sign( aj(r) ) = -1 とします。| S | < λ の場合はどちらも正しくない状態となりますが、先の二つの条件において | S | = λ のときはゼロになるので、| S | < λ のときは aj(r+1) = 0 とすると連続な関数が得られます。以上をまとめると

aj(r+1)=( S - λ ) / Sj,j[ S > λ ]
=0[ | S | ≤ λ ]
=( S + λ ) / Sj,j[ S < -λ ]

という結果になります。この関数は「Soft-Thresholding Operator」と呼ばれます。

先ほど、| S | ≤ λ の領域については 0 とすることで連続となると説明しましたが、その根拠については特に触れていませんでした。元々、最小二乗法での正規方程式は、二次式を偏微分することから各未知数に対して一次式になります。よって、一つの未知数のみを変数とすれば、元の式に対する傾きが線形に増大することを意味します。一方、ペナルティー項の偏微分結果はゼロを境界に符号が反転する不連続な階段関数です。従って、両者を加算した結果はちょうどゼロのところで分断された直線となります。二つの式は

y = Sj,jaj(r+1) - S

y = sign( aj(r)

であり、二つを合わせると

y = Sj,jaj(r+1) - S + sign( aj(r)

となって、そのグラフは以下のようになります。

S > λ| S | ≤ λS < -λ
係数が正になる場合係数がゼロになる場合係数が負になる場合
図 2-3. y = Sj,jaj(r+1) - S + sign( aj(r) )λ のグラフの推移

求めたいのは直線が aj 軸と交差するときの aj の値であり、S > λ ではその値は正の側、S < -λ では負の側にあるのがグラフからわかります。そして、| S | ≤ λ については直線が aj 軸に交わらず、その値は常にゼロであると考えることができます。このため、| S | ≤ λ の範囲では 0 としているわけです。

CDA は「最急降下法 (Gradient Descent Method)」にもよく似た手法です。最急降下法では勾配の方向に最小値となる点を探索し、その点で勾配を求めて次の最小値を探索するという処理を繰り返しますが、CDA の場合は勾配ではなく各軸の方向へ探索を行うところが異なります。また、係数が収束して誤差関数が小さくなると求める未知数がゼロになる傾向が強まるため、一気にゼロへ収束しやすいといった特徴もあります。下図は、係数が 2 つの場合の探索経路を比較したものになります。

最急降下法CDA
最急降下法の探索経路CDAの探索経路
図 2-4. 最急降下法と CDA の探索経路の違い

リッジ回帰ときのように定数項に対してペナルティー項を加えないのであれば、LASSO での正規方程式は

Σk{1→p}( Σi{1→N}( ( xi,j - mxj )( xi,k - mxk )ak ) ) + sign( aj )λ = Σi{1→N}( ( xi,j - mxj )( yi - my ) )

なので、

aj=[ Σi{1→N}( ( xi,j - mxj )( yi - my ) ) - Σk{0→p;k≠j}( Σi{1→N}( ( xi,j - mxj )( xi,k - mxk ) )ak ) - sign( aj )λ ] / Σi{1→N}( ( xi,j - mxj )2 )
=[ Syj - Σk{0→p;k≠j}( Sj,kak ) - sign( aj )λ ] / Sjj

で求めることができます。但し、ここでの Syj と Sjk はそれぞれ平均との差を計算した上で積和を計算した結果であることに注意して下さい。なお、共分散 Vyj = Sjk / N, Vjk = Sjk / N を使う場合は

aj = [ Vyj - Σk{0→p;k≠j}( Vj,kak ) - sign( aj )λ / N ] / Vjj

となります。


LASSO を行うためのサンプル・プログラムを以下に示します。

/*
  LASSO 回帰モデル

  x : 独立変数ベクトル
  y : 従属変数
  lambda : 補正の重みパラメータ
  a : 推定した係数を取得するための配列へのポインタ
  threshold : 収束判定のためのしきい値
  maxCnt : 反復処理の最大処理回数
*/
void LASSORegression( const vector< vector< double > >& x, const vector< double >& y, double lambda,
                      vector< double >* a, double threshold, unsigned int maxCnt )
{
  assert( &x != 0 && &y != 0 && a != 0 );
  assert( threshold >= 0 );

  size_t n; // 標本数
  size_t p; // 独立変数ベクトルのサイズ
  CheckLinearModel( x, y, &n, &p );

  // x, y の平均
  vector< double > mx( p ); // x[i] の平均
  for ( size_t j = 0 ; j < p ; ++j )
    mx[j] = accumulate( x[j].begin(), x[j].end(), 0.0 ) / n;
  double my = accumulate( y.begin(), y.end(), 0.0 ) / n; // y の平均

  // 共分散行列の作成
  vector< vector< double > > sx( p, vector< double >( p ) );
  vector< double > sy( p );
  for ( size_t r = 0 ; r < p ; ++r ) {
    sx[r][r] = sampleVariance( x[r] );
    sy[r] = sampleCovariance( x[r], y );
    for ( size_t c = r + 1 ; c < p ; ++c )
      sx[r][c] = sx[c][r] = sampleCovariance( x[r], x[c] );
  }

  a->assign( p + 1, 0 );

  unsigned int cnt = 0; // 現在の反復回数
  for ( ; cnt < maxCnt ; ++cnt ) {
    bool isConv = true;
    for ( size_t j = 0 ; j < p ; ++j ) {
      // 推定係数を求める
      double s = sy[j] - inner_product( sx[j].begin(), sx[j].end(), a->begin() + 1, 0.0 ) + sx[j][j] * (*a)[j + 1];
      double d =
        ( s > lambda / n ) ? ( s - lambda / n ) / sx[j][j] :
        ( ( s < -lambda / n ) ? ( s + lambda / n ) / sx[j][j] : 0 );
      if ( fabs( (*a)[j + 1] - d ) > threshold ) isConv = false;
      (*a)[j + 1] = d;
    }
    if ( isConv ) break;
  }
  // 定数項の計算
  a->front() = my - inner_product( mx.begin(), mx.end(), a->begin() + 1, 0.0 );

  if ( cnt >= maxCnt )
    throw runtime_error( "Iteration count exceeds maximum." );
}

リッジ回帰の場合では連立方程式を作成するときにペナルティー項を係数行列へ加算して解いていたのに対し、LASSO では通常の最小二乗法で使う連立方程式を作成し、これを解くのではなく要素を利用して上述した反復処理を行う流れになっています。係数を取得する配列 a の最初の要素は定数項を保持するために使うので、反復処理には 2 番目の要素 ( 添字 1 の要素 ) からのデータを利用することに注意して下さい。


サンプル・プログラムを利用して、sin( 2πx ) の多項式近似を次数 9 次に対していくつかの重みパラメータで行った結果を以下に示します。しきい値は 10-6、最大反復回数は 100,000 回としてあり、λ = 0 の場合だけ処理回数が最大値をオーバーしたため途中で処理を中断しています。

図表 2-1. LASSO による係数の推定結果
重みパラメータ推定した多項式のグラフ(LASSO)
λ=0λ=0.001λ=0.01λ=0.1λ=1
係数定数項-1.09E-01-7.60E-021.11E-015.91E-011.12E-01
x7.98E+008.16E+005.31E+00--
x2-9.30E+00-1.70E+01-1.10E+01-2.44E+00-3.87E-01
x3-3.17E+01----
x43.47E+01----
x51.70E+011.14E+012.58E+00--
x6-7.93E+00-2.97E+00--
x7-1.54E+01----
x8-6.64E+00----
x91.13E+01-2.56E+00-1.63E+00-

LASSO では、重みパラメータ λ がゼロのときすでに振れ幅が小さくなっています。λ = 0 なら LASSO はガウス=ザイデル法と全く同じ処理になります。計算は反復法で行われ、前の処理との変動がしきい値以下になれば処理を中断することから、ある程度の誤差は許容範囲内とされます。この結果は、ペナルティー項を用いなくても反復法による処理で精度を上げることができる ( 誤差による過度な影響が抑制できる ) ことを意味します。表の中で "-" で示されてる個所は全てゼロの値です。 λ が大きくなるに従って、ゼロになる個所が増えていく様子がこの結果から読み取れます。無関係な係数をふるい落とすことが目的であれば、やはりペナルティー項は必要となります。


*2-1) 「ラグランジュの未定乗数法」については「確率・統計 (11) 二標本の解析 -1-」の中の「補足 1」でも紹介しています。


<参考文献>
1. 「確率・統計入門」 小針あき宏著 (岩波書店)
2. 「統計数学入門」 本間鶴千代著 (森北出版)
3. 「統計のための行列代数(上)」 D.A.ハーヴィル 著 (シュプリンガー・ジャパン社)
4. 我楽多頓陳館 --- 統計学入門
5. Wikipedia

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system