スパースモデリングはなぜ生まれたか？ 代表的なアルゴリズム「LASSO」の登場

ITエンジニアのためのスパースモデリング入門　第2回

線形回帰とスパースモデリング

　ここからは前節で紹介したスパースモデリングの代表的な手法であるLASSOについて、その背景にある考え方も含めて具体例とともに紹介していきます。サンプルコードはGitHubに公開しています。

• hacarus/codezine-sparse-modeling（GitHub）

線形回帰モデルと最適化問題

　前回「未知数が二つだけれども方程式が一つしかない」ため答えを求められない場合にどうするかということを例に上げ、スパースモデリングの説明をしました。このようにスパースモデリングには、方程式よりも未知数、つまりデータ数よりも説明変数の数のほうが少ない場合にも使うことができるという側面がありますが、説明変数の多くがスパースであるという仮定をすることで説明変数、特徴量の選択をしているという側面もあります。実際の用途としてはこの特徴量選択のためにスパースモデリングを使うケースのほうが多いでしょう。

　具体例として店舗別の売上データ（sales.csv）を使って考えます。このデータには、以下の5つの説明変数と、売上のデータが含まれています。

• 店舗の面積
• 駅からの距離
• 周辺人口
• 競合店数
• 商品点数

　通常は基本統計データや相関係数などをみて必要に応じて変数選択をしたりしますが、今回は比較のためあえて特に何も考えず重回帰分析にかけてみます。

import pandas as pd
from sklearn.linear_model import LinearRegression

def print_model(model, X, y):
    print(f"Coefficients: {list(lr.coef_)}")
    print(f"Intercept: {lr.intercept_}")
    print(f"R square: {lr.score(x, y)}")

df = pd.read_csv('sales.csv')
X = df.iloc[:, 1:-1]
y = df["売上 (千円)"]
lr = LinearRegression()
lr.fit(X, y)
print_model(lr, X, y)

　これを実行すると以下の結果を得ます。

Coefficients: [2.274973606196981, -104.81754452236632, 3.070159435067054, -20.663131503414206, -0.4301726133789289]
Intercept: 857.3883792692491
R square: 0.9507107022338874

　売上を\(y\)、面積を\(x_1\)、駅からの距離を \(x_2\)、周辺人口を \(x_3\)、競合店数を\(x_4\)、商品点数を\(x_5\)とすると回帰式は以下のようになります。

　\(y = 2.27x_1 - 104.82x_2 + 3.07x_3 - 20.66x_4 - 0.43x_5 + 857.39\)

　学習に使ったデータから得られた値であるため少し高めの数値がでる傾向にありますが、決定係数（R square）も0.95で悪くない値です。ただ、一方で仮に真実の解がスパースだとすると、すべての変数が使われていることは、先に述べたオッカムの剃刀の考え方からは外れます。また、この例ですとまだ説明変数が少ないので実感しにくいかもしれませんが、入力する説明変数が多くなればなるほど、すべての変数が目的変数を求めるのに一定以上寄与することは考えづらいです。求める精度にもよりますが、大抵の場合はシンプルなモデルで十分実用に足ることが多いでしょう。

L₀ノルム最適化

　スパースモデリングを使うことで説明変数の選択をすると同時に係数を求めることができます。スパースモデリングにはさまざまな手法がありますが、真実の解が十分にスパースであれば、真実の解が得られることが知られている\(L_0\)最適化を最初に紹介します。

　売上のベクトルを\(y\)、説明変数となるデータの行列を\(\Phi\)、求める係数のベクトルを\(x\)とすると以下の式を解けばよいことになります。

　\(y = \Phi x\)

　\(L_0\)最適化は\(x\)の\(L_0\)ノルム^[1]を最小化するというものです。ベクトル\(x\)の\(L_0\)ノルムとは\(x\)の非ゼロ要素の個数です。つまり\(L_0\)最適化とは\(x\)の非ゼロ要素の個数ができるだけ少ないものを求めるというものです。この問題の最も直接的な解き方は以下で説明する総当たり法です。ここで、最もスパースな解を\(x^*\)とし、その時の\(L_0\)ノルム\(\|x^*\|_0\)を\(k\)とします。

1. もし\(y = 0\)なら\(x^* = 0\)として終了
2. \(k := 1\)
3. \(\|x\|_0 = k\)となる全ての\(x\)について\(y = \Phi x\)が解を持つか調べ、もし解を持つならそれを解として終了
4. \(k := k + 1\)として 3. に戻る

　このようにアルゴリズムそのものは非常に単純ではありますが、これは組合せ最適化でありベクトル\(x\)の次元数が多くなると組合せ爆発が発生し、現実的な時間内に計算が終わらない可能性が高くなります。

L₁ノルム最適化とLASSO

　\(L_0\)最適化において組合せ爆発が生じるのは\(L_0\)ノルムを使っているからです。この\(L_0\)ノルムは以下の式で定義される\(L_1\)ノルムで近似することができます。

　\(\|x\|_1 = \sum_{i=1}^{n} |x_i|\)

　式から分かる通り\(L_1\)ノルムは\(x\)の要素の絶対値の総和であり、マンハッタン距離と呼ばれるものと同じです。

　近似により\(L_0\)ノルム最小化問題はn次元空間内の超平面\(y = \Phi x\)上の点のうち\(L_1\)ノルム最小となる点を求める問題となります。\(x\)が 2次元のとき図で表すと以下のようになります。

L1ノルム最適化

　\(L_1\)ノルムが一定となる点を結ぶと頂点が座標軸の上にあるひし形になります。このひし形を面積がゼロとなる中心点から徐々に大きくしていって最初に直線とぶつかる点が最適解となります。大抵の場合は座標軸上のひし形の頂点とぶつかるため、どちらかの座標がゼロとなりスパースな解が得られます。

　ここまではデータにノイズが混入していないことを前提に説明をしてきましたが、現実の世界においてはデータにノイズが混入していることが普通です。近似曲線を計算するのに最小二乗法を使うのはよく知られた方法です。しかしノイズの混入したデータにおいては冒頭で説明したように過学習が起こることがあります。詳細な説明は省きますが、この過学習を防ぐために正則化項として\(L_1\)ノルムを使うことができます。問題としては以下の式の値を最小化する最適化問題となります。

　\(\frac{1}{2}\|\Phi x - y\|^2_2 + \lambda \|x\|_1\)

　この最適化問題のことをLASSOといいます。LASSOは、スパースモデリングが広く知られるきっかけになったアルゴリズムです。LASSOの実装はscikit learnにありますので、重回帰分析にかけたデータをLASSOを使って線形回帰分析してみます。

from sklearn.linear_model import Lasso

clf = Lasso(alpha=5)
clf.fit(X, y)
print_model(clf, X, y)

　実行結果は以下のとおりです。

Coefficients: [2.1450202545748462, -103.13447263043422, 3.026543503672327, -18.8770989622993, -0.0]
Intercept: 838.5335314111367
R square: 0.9501869111042888

　5つめの変数の重みがゼロになっておりスパースモデリングの特徴の一つである変数選択されていることがわかります。決定係数も変数をすべて使っている重回帰の場合とほとんど変わりません。

L₀ノルム最適化を貪欲に解く Orthogonal Matching Pursuit

　\(L_0\)最適化問題で総当たり法を使うと組合せ爆発が発生して実用に耐えないことを説明しました。このような場合には貪欲法（greedy method）が有効です。貪欲法では問題をいくつかの部分問題に分割し、部分問題の解のうち最も良いものを選択するという方法です。

　\(L_0\)最適化問題を貪欲法で解く代表的な手法にOMP（Orthogonal Matching Pursuit）があります。OMPは非ゼロとなるベクトル\(x\)の要素を指すインデックス集合Sを見つけ出すアルゴリズムです。この非ゼロ係数のインデックス集合をサポートといいます。初めはサポートは空集合であるとして、\(x\)を基底\(\Phi\)の線型結合で近似した時の残差を最小にするように新たな基底をサポート集合に一つ一つ追加していき、サポートに含まれる基底のみで\(x\)を近似した時の残差がしきい値以下になったら停止するというものです。残差の低減に寄与する基底を順次選択していく貪欲法であり、解の最適性は保証されませんが、多くの場合優れた近似値が得られます。OMPを売上データに適用したコードを以下に示します。

from sklearn.linear_model import OrthogonalMatchingPursuit

omp = OrthogonalMatchingPursuit(n_nonzero_coefs=4)
omp.fit(X, y)
print_model(omp, X, y)

　これを実行すると以下の出力を得ます。

Coefficients: [2.2700090488584483, -104.80224608073459, 3.068194474356952, -20.666032177559337, 0.0]
Intercept: 842.4643015808057
R square: 0.950662082171787

　OrthogonalMatchingPursuitでは係数が非ゼロとなる特徴量の数は自動で計算されず、コンストラクタ引数にn_nonzero_coefs=4のようにして渡します。指定しない場合は特徴量の数の10％がデフォルト値として使われます。決定係数を見ながら調整していくのが良いでしょう。

[1] \(L_0\)ノルムは数学的には斉次性を満たさないため厳密な意味でのノルムではありませんが、慣例に従いここではノルムと呼びます。