『ベイズ推論による機械学習入門』で解説されていた非負値行列因子分解 (Non-negative matrix factorization, NMF or NNMF) を R でやってみます.

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)
- 作者: 須山敦志,杉山将
- 出版社/メーカー: 講談社
- 発売日: 2017/10/21
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (1件) を見る
非負値行列因子分解では各要素が非負の値を持つ行列 V (D 行 N 列) を, 非負の値を持つ行列 W (D 行 M 列) と H (M 行 N 列) に分解します.
(『ベイズ推論による機械学習入門』では ,
となってるけど別に自然数である必要はないんじゃないかな? 非負の実数でいいと思う)
M は分析者が自分で決めます.
これによりデータの量が DN から DM + MH に削減できます. 場合によっては解釈が容易になることもあります.
Non-negative matrix factorization - Wikipedia
Rによる実装例
非負値行列因子分解を実行するコードはこんな風になりました.
アルゴリズムの詳細については『ベイズ推論による機械学習入門』を参照してください.
a_W, b_W, a_H, b_H はガンマ分布のパラメータ(ハイパーパラメータ)です.
NMF <- function(X,M=2,a_W=1,b_W=1,a_H=1,b_H=1,tol=1e-2,maxit=10000, seed=1234){ set.seed(seed) D <-nrow(X) N <-ncol(X) W<-matrix(rgamma(D*M,a_W,b_W),D,M) H<-matrix(rgamma(M*N,a_H,b_H),M,N) logW0 <-log(W) logH0 <-log(H) for(k in 1:maxit){ pi <-lapply(1:M,function(m)exp(outer(logW0[,m],logH0[m,],"+"))) den <-Reduce("+",pi) pi_norm <-lapply(1:M,function(i)pi[[i]]/den) S <- lapply(pi_norm,function(y)X*y) ahat_W <-sapply(S,function(y)apply(y,1,sum))+a_W bhat_W <-apply(H,1,sum)+b_W W <-sweep(ahat_W,2,bhat_W,"/") ahat_H <-t(sapply(S,function(y)apply(y,2,sum)))+a_H bhat_H <-apply(W,2,sum)+b_H H <-sweep(ahat_H,1,bhat_H,"/") logW <-sweep(digamma(ahat_W),2,log(bhat_W)) logH <-sweep(digamma(ahat_H),1,log(bhat_H)) if(all(abs(logH-logH0)<tol) & all(abs(logW-logW0)<tol)){ break } logW0<-logW logH0<-logH } return(list(W=W,H=H,iter=k)) }
正しく書けてるか自信がないので, 実データで試してみます.
尻目
分解する行列として, 今回は,
蕪村妖怪絵巻 - Wikipedia
より「京都帷子が辻のぬっぽり坊主」の画像データを選びました.
この妖怪は水木しげるが『尻目』という名前で紹介して有名になりました.

- 作者: 水木しげる
- 出版社/メーカー: 岩波書店
- 発売日: 1992/07/20
- メディア: 新書
- この商品を含むブログ (15件) を見る
jpeg パッケージで画像データを読み込むと, 400行772列の行列になりました.
一応, 今回はポアソン分布を使っているので, 100掛けて丸めて, 行列の値が非負の整数になるようにしています.
library(tidyverse) library(jpeg) mat4img <- function(mat1){ t(apply(mat1,2,rev)) } shirime <-readJPEG("Buson_Nopperabo.jpg") X <- round(shirime[,,1]*100) image(mat4img(X))
M の値をいくつかかえて, WH がどのくらい元の画像に近いか調べてみます.
res2<-NMF(X,M=2) res10<-NMF(X,M=10) res50<-NMF(X,M=50)
M=2 のとき:
image(mat4img(res2$W%*%res2$H),main="M=2")
がたがたですね.
M=10 のとき:
image(mat4img(res10$W%*%res10$H),main="M=10")
ちょっと尻目の輪郭らしきものが見えます.
M=50のとき:
image(mat4img(res50$W%*%res50$H),main="M=50")
尻目になった. これは尻目になったと言っていいでしょう.
感想
- 行列の分解なんかが確率モデルとして書けるというところがおもしろかった
- けっこう時間がかかる(M=50のときで15分くらい)
- 繰り返し計算が多いので R 向きのタスクではないかも