Python
機械学習
OriginalAkatsukiDay 23

Alternating Direction Method of Multipliers (ADMM) で Lasso 回帰

この記事は Akatsuki Advent Calendar 2017 の 23 日目の記事です。
22日目: あなたのメモリはどこから来てる?malloc入門

株式会社アカツキでサーバーサイドをやっている新卒エンジニアです。
今回の記事では、機械学習ネタとしてスパースモデリングをやってみます。

深層学習とスパースモデリング

近年、深層学習が非常に注目を集めています。
深層学習では特徴量を自動的に抽出してくれるため、大量のデータさえあれば、モデリング(=数式化)をすることなく誰でも利用することが可能です。
しかしながら、データサイズが小さいと過学習を起こしてしまったり、あるいは充分な量のデータを使って適切な学習をした場合でも、入力と出力の関係性を説明することが難しいという問題があります。

一方でスパースモデリングは、適切なモデリングを要するものの、必ずしも大量のデータを必要とせず、しかも高次元データの本質的な部分を自動的に選択してくれるという性質があります。
これによって、例えば「ユーザー継続率の要因になりうる変数は大量にあるが、そのうち一体どれが本当に寄与しているのか?」といった分析が可能になります。

今回は、スパースモデリングの初歩的な手法である Lasso 回帰をやってみます。
Lasso 回帰に使えるアルゴリズムはいくつかあり、例えば scikit-learn の sklearn.linear_model.Lasso では Coordinate Descent が用いられていますが、今回は色々と応用の効く Alternating Direction Method of Multipliers (ADMM) という手法を使います。

解きたい問題

Lasso

Lasso (least absolute shrinkage and selection operator) 回帰は、最小二乗法の目的関数に L1 正則化項を加えた、以下の最適化問題を指します。

\mathrm{min} ~ \frac{1}{2N} \|Ax - b\|_2^2 + \lambda \|x\|_1

各変数の次元は $A \in \mathbb{R}^{n \times m}$, $b \in \mathbb{R}^{n}$, $x \in \mathbb{R}^{m}$ です。
それぞれの意味は

  • $A$: 計画行列。各行は説明変数ベクトル(を転置したもの)。
  • $b$: 観測ベクトル。各要素は1次元の観測値。
  • $x$: 回帰係数ベクトル。最適化問題を解くことで得ようとしているもの。

です。

L1 正則化項を用いることには、出力に影響のある入力変数が自動的に選択されるというメリットがあります。
上記の目的関数で $\lambda |x|_1$ を $\lambda |x|_2$ にすると Ridge 回帰になりますが、$\lambda |x|_2$ が $x$ の長さを小さくする働きをする一方、$\lambda |x|_1$ は $x$ の要素のうちできるだけ多くを0にする効果を持ちます。

より詳しくは、Wikipedia、もっともっと詳しくは Statistical Learning with Sparsity: The Lasso and Generalizations などの書籍を参照してください。

データセット

今回はボストンにおける住宅価格を推定してみます。データセットは scikit-learn のものを使用します。

入力の13次元ベクトルは「部屋数」のような住宅価格に影響のありそうなものや、「チャールズ川流域か否か」といった需要がニッチ(?)過ぎて影響の少なそうなものなど、様々な情報で構成されています。
sklearn.datasets.load_boston で読み込んだデータの DESCR 属性で表示される情報のうち、ベクトル要素に関する部分を以下に示します。

:Median Value (attribute 14) is usually the target
  :Attribute Information (in order):
      - CRIM     per capita crime rate by town
      - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
      - INDUS    proportion of non-retail business acres per town
      - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
      - NOX      nitric oxides concentration (parts per 10 million)
      - RM       average number of rooms per dwelling
      - AGE      proportion of owner-occupied units built prior to 1940
      - DIS      weighted distances to five Boston employment centres
      - RAD      index of accessibility to radial highways
      - TAX      full-value property-tax rate per $10,000
      - PTRATIO  pupil-teacher ratio by town
      - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
      - LSTAT    % lower status of the population
      - MEDV     Median value of owner-occupied homes in $1000's

これを Lasso の式に当てはめると、$A \in \mathbb{R}^{506 \times 13}$ の各行は1データの13属性を並べたもの、b の各要素は1データの MEDV を入れたものになります。

ADMM

Alternating Direction Method of Multipliers (ADMM) は拡張ラグランジュ法を改良した手法で、

  • 経験的に収束が速い。
  • 何らかの変換によってスパース性が期待できる基底が見つかる場合、特段の苦労なく適用可能。

という特長があります [今日からできるスパースモデリング]。

この記事では、「Alternating Direction Method of Multipliers」を参照しつつ、データサイズ $N (= 506)$ を考慮に入れて ADMM のアルゴリズムを使います。

一般的な ADMM

ADMM では以下の形の最適化問題を扱います。

\begin{array}{ll}
\mathrm{ min } & f(x) + g(z) \\
\mathrm{ s.t. } & Ax + Bz = c
\end{array}

この問題の拡張ラグランジュ関数は

L_\rho (x, z, y) = f(x) + g(z)  + y^{\mathrm{T}}(Ax + Bz - c) + \frac{\rho}{2} \|Ax + Bz -c \|_2^2

となります。ここで、ADMM は罰金項の係数 $\rho$ を動かすことなく、次のように $x, z, y$ を更新していきます。

\begin{array}{rcl}
x^{k+1} &:=& \mathrm{argmin}_x ~ L_\rho (x, z^k, y^k) \\
z^{k+1} &:=& \mathrm{argmin}_z ~ L_\rho (x^{k+1}, z, y^k) \\
y^{k+1} &:=& y^k + \rho(Ax^{k+1} + Bz^{k+1} - c)
\end{array}

ADMM による Lasso 回帰

解きたい問題を再掲します。

\mathrm{min} ~ \frac{1}{2N} \|Ax - b\|_2^2 + \lambda \|x\|_1

これを ADMM の扱う問題形式に書き換えると以下のようになります。

\begin{array}{ll}
\mathrm{ min } & \frac{1}{2N} \|Ax - b\|_2^2 + \lambda \|z\|_1 \\
\mathrm{ s.t. } & x - z = 0
\end{array}

拡張ラグランジュ関数は次の通りです。

L_\rho (x, z, y) = \frac{1}{2N} \|Ax - b\|_2^2 + \lambda \|z\|_1  + y^{\mathrm{T}}(x - z) + \frac{\rho}{2} \|x - z \|_2^2

よって、ADMM を Lasso 回帰に当てはめると、以下の更新式が得られます。

\begin{array}{rcl}
x^{k+1} &:=& (\frac{1}{N} A^{\mathrm{T}}A + \rho I)^{-1} (\frac{1}{N} A^{\mathrm{T}} b + \rho z^k - y) \\
z^{k+1} &:=& S_{\lambda / \rho} (x^{k+1} + \frac{1}{\rho} y^k) \\
y^{k+1} &:=& y^k + \rho(x^{k+1} - z^{k+1})
\end{array}

ここに、$S_a(x)$ は以下のように定義される軟判定しきい値関数(をベクトルの要素単位で適用したもの)です。

S_a(x) = \begin{cases}
  x - a & (x > a) \\
  0 & (-a \leq x \leq a) \\
  -x + a & (x < -a)
\end{cases}

実験

以下の設定で実験を行いました。といっても、ちゃんとした実験ではなく、sckit-learn の結果と比べて「合ってるな」というぐらいのゆるふわな感じです。ご容赦ください。

  • 初期値
    • $x^0 := \frac{1}{N} A^{\mathrm{T}} b$
    • $z^0 := x^0$
    • $y^0 := 0$
  • パラメータ
    • $\lambda := 1$
    • $\rho := 1$

ソースコードの一部を以下に示します。 (全文はGitHubに置いてあります: https://github.com/s-capybara/lasso_by_admm)
コーディングには「Lassoの理論と実装 -スパースな解の推定アルゴリズム-」 を参考にしました。

class Admm(BaseEstimator, RegressorMixin):
    def __init__(self, lambd=1.0, rho=1.0, max_iter=1000):
        self.lambd = lambd
        self.rho = rho
        self.threshold = lambd / rho
        self.max_iter = max_iter
        self.coef_ = None
        self.intercept_ = 0.0

    def _soft_threshold(self, x):
        t = self.threshold

        positive_indexes = x >= t
        negative_indexes = x <= t
        zero_indexes = abs(x) <= t

        y = np.zeros(x.shape)
        y[positive_indexes] = x[positive_indexes] - t
        y[negative_indexes] = x[negative_indexes] + t
        y[zero_indexes] = 0.0

        return y

    def fit(self, A, b):
        N = A.shape[0]
        M = A.shape[1]
        inv_matrix = np.linalg.inv(np.dot(A.T, A) / N + self.rho * np.identity(M))

        x = np.dot(A.T, b) / N
        z = x.copy()
        y = np.zeros(len(x))

        for iteration in range(self.max_iter):
            x = np.dot(inv_matrix, np.dot(A.T, b) / N + self.rho * z - y)
            z = self._soft_threshold(x + y / self.rho)
            y += self.rho * (x - z)

        self.coef_ = x

        return self

    def predict(self, X):
        y = np.dot(X, self.coef_)
        return y
model = Admm(lambd=1.0, rho=1.0, max_iter=1000)
model.fit(A, b)
print("Lasso by ADMM")
print(model.coef_)
Lasso by ADMM
[-0.0000000000  0.0000000000 -0.0000000000 -0.0000000000  0.0000000000
  2.7131072809  0.0000000000  0.0000000000  0.0000000000  0.0000000000
 -1.3434986189  0.1807938799 -3.5436116588]

13個の回帰係数ベクトル要素のうち、9個が0になりました。スパースな解が得られています。
RM, PTRATIO, B, LSTAT がボストンの住宅価格の中央値に影響があるとの結果です。

比較として scikit-learn (Coordinate Descent) で解いたものが以下になります。

from sklearn import linear_model

model = linear_model.Lasso(alpha=1.0, max_iter=1000)
model.fit(A, b)
print("Lasso by sickit-learn (coordinate descent)")
print(model.coef_)
Lasso by sickit-learn (coordinate descent)
[-0.0000000000  0.0000000000 -0.0000000000  0.0000000000 -0.0000000000
  2.7133533427 -0.0000000000 -0.0000000000 -0.0000000000 -0.0000000000
 -1.3435484031  0.1809553656 -3.5433828848]

同じ結果が得られました。

応用

ADMM を使えば Fused Lasso や、はたまたガウスマルコフ確率場も扱えるらしく、色々楽しそうです。