CodeZine(コードジン)

特集ページ一覧

最先端のスパースモデリング~HMLassoとPliable Lasso~

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

  • ブックマーク
  • LINEで送る
  • このエントリーをはてなブックマークに追加
2020/09/24 11:00

 本連載は「これから機械学習に取り組みたい」「ディープラーニングや機械学習を使った経験がある」といったエンジニアに向けて、データ量が少なくても分析が実現できる「スパースモデリング」という手法を紹介します。 前回はスパースモデリングの画像処理への発展的な応用として、画像の欠損補間、異常検知、超解像の3つを紹介しました。今回はスパースモデリングの最近の学術分野におけるスパースモデリングの発展の様子、最新の手法をご紹介します。

目次

スパースモデリングにおける数理モデル研究の重要性

 データが少なく説明責任が求められる状況に力を発揮するスパースモデリングの中には、様々な手法があります。例えば、第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)

 学習された回帰係数を出力してみると、半数ものデータが欠損していたにも関わらず、欠損のないときと結果が大きくは変わっていないことがわかると思います。

 モデルを工夫することで、このように強力なアルゴリズムを生み出すことができたことが見て取れます。


  • ブックマーク
  • LINEで送る
  • このエントリーをはてなブックマークに追加

著者プロフィール

  • 田辺 広樹(タナベヒロキ)

     京都大学工学部情報学科に首席入学・首席卒業した後に進学した京都大学大学院情報学研究科修士課程を1.5年で早期修了し、現在は同研究科の後期課程に在籍中。日本学術振興会特別研究員DC1。専門は数理最適化だが、Webや機械学習など、幅広いエンジニアリングスキルを持つ。「ざるご」という名義でYouTube...

バックナンバー

連載:ITエンジニアのためのスパースモデリング入門
All contents copyright © 2005-2020 Shoeisha Co., Ltd. All rights reserved. ver.1.5