スパースモデリングにおける数理モデル研究の重要性
データが少なく説明責任が求められる状況に力を発揮するスパースモデリングの中には、様々な手法があります。例えば、第2、3回で紹介したLASSOや、第4、5回で紹介した辞書学習などはその中の代表的なものとして知られています。スパースモデリングを扱うエンジニアには、問題の質に応じて、これらの中から最も適切な手法を選択することが求められています。
一方、機械学習の社会実装が急速に進む昨今では、個別の問題に対してまだ適切な手法が存在しないということもあります。そこで重要となるのは、そのような問題を捉えることのできる新しい数理モデルを提案し、そのモデルを解く手法を開発することです。国内外の研究機関では、より広い範囲の問題を解決することができるよう、日夜新しいスパースモデリングの数理モデルが開発されています。
連載の最後となる今回の記事では、スパースモデリングの代表的な手法であるLASSOに対して、近年開発された発展的なモデルを紹介していきます。
新しい「LASSO」たち（1）
ここでは、近年開発されたLASSOの発展形を2つ紹介します。従来のLASSOでは、入力信号\(y\)、計画行列\(X\)に対して、次の最適化問題を解くことにより、回帰係数\(\hat{\beta}\)を推定していました。
\( \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \hat{\beta} = \argmin_{\beta} \left[ \frac{1}{2} \|y - X\beta\|_2^2 + \lambda \|\beta\|_1 \right] \)
ここで、第2項は正則化項、\(\lambda\)は正則化パラメータと呼ばれます。この最適化問題を解くことにより、いくつかのパラメータがゼロになるような回帰係数\(\hat{\beta}\)を推定することができました（詳しくは、第3回の記事を参照してください）。
欠損値を含むデータにも使える「HMLasso」
現実の問題では、収集したデータに欠損値が含まれているということがしばしばあります。そこで昨年、東芝と統計数理研究所によって開発されたのが、欠損値を含む問題に対する新しいLASSOである「HMLasso」です[1]。ここでは、HMLassoの簡単なアルゴリズムと、実際のプログラムを解説していきます。
注
[1] Takada, Masaaki, Hironori Fujisawa, and Takeichiro Nishikawa. "HMLasso: lasso with high missing rate." arXiv preprint arXiv:1811.00255 (2018).
従来のLASSOは、欠損値を含む問題に対してそのまま適用することはできません。そこで、LASSOの式を次のように変形して考えます。
\( \hat{\beta} = \argmin_{\beta} \left[ \frac{1}{2} \beta^\top S \beta - \rho^\top \beta + \lambda \|\beta\|_1 \right] \)
ただし、\(S = \frac{1}{2} X^\top X, \rho = X^\top y\)とします。ここで、\(S\)は計画行列\(X\)の分散共分散行列で、半正定値であるため、この最適化問題は凸最適化問題であり、簡単に解くことができることに注意します。
いま、計画行列\(X\)に欠損が含まれているとすると、上の式をそのまま計算することはできません。そこで、観測できた値のみを用いて、\( S = S^\text{pair} = \left( S^\text{pair}_{jk} \right) = \sum_{i \in I_{jk}} X_{ij} X_{ik}, \rho = \rho^\text{pair} = \left( \rho^\text{pair}_j \right) = \sum_{i \in I_{jj}} X_{ij} y_i \)と、\(S, \rho\)を推定することが考えられます（ただし、\(I_{jk}\)は\(X_{ij}\)と\(X_{ik}\)がともに観測されているような\(i\)の集合とします）。
しかし、このようにして推定した\(S\)は半正定値行列ではないため、上の最適化問題は非凸最適化問題になってしまい、解くことが難しくなってしまいます。HMLassoでは、この問題を解決するために、\(S^\text{pair}\)に近い半正定値行列を、重み行列\(W\)を用いた次の最適化問題を解くことによって求めます。
\( S = \argmin_{\Sigma \succeq O} \| W \circ (\Sigma - S^\text{pair})\|_F^2 \)
ただし、\(\circ\)は要素ごとの積を表します。このようにして得られた半正定値行列\(S\)を分散共分散行列の推定値として用いることで、回帰係数\(\hat{\beta}\)を推定することができます。
次に、実際にプログラムを用いてアルゴリズムを動かしてみましょう。同じ内容のコードを載せたJupyter Notebookを以下のリポジトリで公開しているので、興味ある方は実際に触ってみてください。
ここでは、R言語を用いて実装を行っていくことにします。従来のLASSOについてはglmnet、HMLassoについてはhmlassoというパッケージがそれぞれCRANで公開されているので、これらをインストールしてみましょう。
# Install package for Lasso
if(!require("glmnet")){
install.packages("glmnet")
}
# Install package for HMLasso
if(!require("hmlasso")){
install.packages("hmlasso")
}
まずは、欠損値を含まないデータに対して両アルゴリズムを適用してみます。データには、glmnetに含まれているQuickStartExampleデータセットを用います。このデータに対して、Cross Validationで推定したパラメータを用いて学習を行ってみましょう。
library(glmnet) library(hmlasso) data(QuickStartExample) # Cross Validation for Lasso cv_fit_glmnet <- cv.glmnet(x, y) fit_glmnet <- glmnet(x, y, lambda=cv_fit_glmnet$lambda.min) # Cross Validation for HMLasso cv_fit_hmlasso <- cv.hmlasso(x, y) fit_hmlasso <- hmlasso(x, y, lambda=cv_fit_hmlasso$lambda.min)
学習した回帰係数をそれぞれ出力してみると、どちらのアルゴリズムでもほぼ同じ回帰係数が出力されていることがわかると思います。
次に、欠損値を含む問題を考えてみます。先ほどと同じQuickStartExampleのデータのうち、50%ほどを無作為に欠損させてみます。
library(glmnet) data(QuickStartExample) x_miss <- x len <- length(x) misses <- sample(1:len, len%/%2) x_miss[misses] <- NA
半数が欠損値であるこのデータに、普通のLASSOを適用しても、エラーが出てしまいます。
fit_glmnet <- glmnet(x_miss, y) # Error in elnet(x, is.sparse, ix, jx, y, weights, offset, type.gaussian, : NA/NaN/Inf in foreign function call (arg 5) # Traceback: # # 1. glmnet(x_miss, y) # 2. elnet(x, is.sparse, ix, jx, y, weights, offset, type.gaussian, # . alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, flmin, ulam, # . thresh, isd, intr, vnames, maxit)
しかし、同じデータにHMLassoを適用すると、問題なく学習を行うことができます。先ほどと同様に、Cross Validationで推定したパラメータを学習に使ってみます。
# Cross Validation for HMLasso cv_fit_hmlasso <-cv.hmlasso(x_miss, y) fit_hmlasso <- hmlasso(x_miss, y, lambda=cv_fit_hmlasso$lambda.min)
学習された回帰係数を出力してみると、半数ものデータが欠損していたにも関わらず、欠損のないときと結果が大きくは変わっていないことがわかると思います。
モデルを工夫することで、このように強力なアルゴリズムを生み出すことができたことが見て取れます。
