非負値行列因子分解をRで（ベイズ推論による機械学習入門）

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

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

非負値行列因子分解では各要素が非負の値を持つ行列 V (D 行 N 列) を, 非負の値を持つ行列 W (D 行 M 列) と H (M 行 N 列)　に分解します.

(『ベイズ推論による機械学習入門』では , となってるけど別に自然数である必要はないんじゃないかな？　非負の実数でいいと思う)

M は分析者が自分で決めます.

これによりデータの量が DN から DM + MH に削減できます. 場合によっては解釈が容易になることもあります.

Non-negative matrix factorization - Wikipedia

モデル

となるような潜在変数 S を導入します.

はパラメータ のポアソン分布に従い, W, H はそれぞれガンマ分布に従うことにします.

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
より「京都帷子が辻のぬっぽり坊主」の画像データを選びました.

この妖怪は水木しげるが『尻目』という名前で紹介して有名になりました.

カラー版 妖怪画談 (岩波新書)

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 向きのタスクではないかも