BPRで遊んだよ！ - 甲斐性なしのブログ

前回の記事の実装編です。

前回ちょっとだけ述べた、last.fmのデータがここに転がってたので、
こいつで、BPRアルゴリズムを試してみることにしました。
使用したのは、usersha1-artmbid-artname-plays.tsvファイルのデータのみです。
usersha1-profile.tsvは今回は使用してません。
usersha1-artmbid-artname-plays.tsvのフォーマットですが、
ユーザID、アーティストID、アーティスト名、再生数がタブ区切りで記述されています。
ただし、アーティストIDが付けれられていないアーティストもいるので、
その場合はアーティストIDの部分は空文字列です。

ただし、そのまま使うのは、色々と厳しいので、
以下のルールでデータを抽出することにしました。

・ユーザ $u$ が、アーティスト $i$ を100回以上再生している場合に、
$u$ は $i$ を好んでいるとみなす（ $S_{ui} = 1$ ）。
・100再生以上している場合、再生数に差があっても嗜好度に差をつけない。
・100再生未満のデータは学習に利用しない。
・アーティストIDがないアーティストは抽出しない。

従って、100回以上再生したアーティストがいないユーザや、
誰からも100回以上再生されていないアーティスト、
IDがないアーティストのデータは、今回の実験では除外され、
推薦されることはありません。

このルールによって抽出されたデータ以下の通り。
ユーザー数：7748
アーティスト数：28185

ただ、このデータをただアルゴリズムに突っ込んで、
精度出しましたって言うのは微妙な気がするので、
今回、次の3人の疑似ユーザを学習に紛れ込ませ、
どんなアーティストが推薦されるかをみることにしました。

・某アニメの元ネタくん（ユーザID：K-ON）
好きなアーティスト
p-model、電気グルーヴ、triceratops
soft ballet、筋肉少女帯、the pillows
計6組
（ニコニコ大百科の記事を参考にしました）

・暗いよ70年代後半〜80年代くん（ユーザID：gotham）
好きなアーティスト
joy division、virgin prunes、the birthday party、swans
sisters of mercy、theatre of hate、bauhaus、cristian death
killing joke、Siouxsie & The Banshees、Alien Sex Fiend
Mission、chameleons
計13組

・最近のジョギングのお供だよくん（ユーザID：joging）
好きなアーティスト
sonic youth、melt-banana、!!!、gang of four
echo & the bunnymen、the rapture、the stills、kula shaker
mars volta、doors、espers、nick drake、柴田淳
 光田康典、植松伸夫、L'Arc en ciel、凛として時雨
 king crimson、yes、screaming headless torsos
計20組

これで準備完了。
後は、このユーザ数7751（7748+3）、アーティスト数28185のデータを
BPRアルゴリズムで学習して、推薦結果を得るだけです。

てことで、ソースコードです。
今回もpythonで書きました。
以下は、行列 $S$ からBPRのMatrix Factorization（BPR MF）を行い、
WとHを取得する関数です。
引数ですが、Sは論文に示されている $\|U\| \times \|I\|$ の行列 $S$ です。
未観測の要素は0で補間します。
つまり、ユーザ $u$ がアーティスト $i$ を好んでいるなら $S_{ui}=1$ 、
そうじゃない（好みじゃない、または、知らない）なら $S_{ui}=0$ となります。
kは $W$ と $H$ の列数（隠れ変数の数ともいえます）、
alphaは学習係数 $\alpha$ 、Lamは正則化パラメータ $\lambda$ です。

def MatrixFactorizationBPR(S , k , alpha , Lam):
    #W , Hを初期化
    UserCount = S.shape[0]
    ItemCount = S.shape[1]
    W = numpy.random.randn(UserCount , k)
    H = numpy.random.randn(ItemCount , k)
    W_next = W.copy()
    H_next = H.copy()

    BPRopt_pre=0.0
    count = 0

    while(1):
        #ループ回数をカウント
        count = count + 1

        #u,i,jをランダムサンプリング
        (u,i,j) = Get_uijFromArray(S);

        #xuijを計算
        x_uij = numpy.dot(W[u , :] , H[i , :]) - numpy.dot(W[u , :] , H[j , :]);

        #次のWとHを計算
        G =  numpy.exp(-x_uij) / (1 + numpy.exp(-x_uij));
        W_next[u , :] = W[u , :] + alpha * (G * (H[i , :] - H[j , :]) - Lam * W[u , :]);
        H_next[i , :] = H[i , :] + alpha * (G * (W[u , :]) - Lam * H[i , :]);
        H_next[j , :] = H[j , :] + alpha * (G * (-W[u , :]) - Lam * H[j , :]);

        #WとHを更新
        W = W_next.copy();
        H = H_next.copy();

        #計算に時間がかかるので、50万回に1回のみ収束判定
        if count%500000 == 0:
            BPRopt = GetBprOpt(S , W , H , Lam);

            #BPRoptの値で収束判定
            if(abs((BPRopt - BPRopt_pre))<0.001):
                break;
            BPRopt_pre = BPRopt;

    return W , H

MatrixFactorizationBPR内にある、関数Get_uijFromArrayは、
(u , i , j)の三つ組をランダムサンプリングする関数、
また、GetBprOptは、 $BPR-Opt$ を計算する関数です。
これらの関数はそれぞれ、以下の通りです。

Get_uijFromArray

def Get_uijFromArray(S):
    #ユーザ数とアイテム数をカウント
    UserCount = S.shape[0];
    ItemCount = S.shape[1];

    #ユーザをランダムに選択
    u = numpy.random.randint(UserCount);

    #選択したユーザの嗜好情報を取得
    S_u = S[u,:]
    #集合I+とI-を生成
    IplusNo = numpy.where(S_u == 1)[0];
    IMinusNo = numpy.where(S_u == 0)[0];

    #I+集合からランダムに一つサンプリング
    RandMax = IplusNo.shape[0];
    I_idx = numpy.random.randint(RandMax);
    i = IplusNo[I_idx]

    #I-集合からランダムに一つサンプリング
    RandMax = IMinusNo.shape[0];
    J_idx = numpy.random.randint(RandMax);
    j = IMinusNo[J_idx]

    return u , i , j

GetBprOpt
なお、論文通りのBPRoptを出力すると、
かなりでかい値になって、収束判定しづらくなるので、
最後に足し合わせた要素数（ソース内のZu）で割るようにしています。

def GetBprOpt(S , W , H , Lam):
    Bpropt = 0.0;
    Zu = 0.0;

    for u  in range(W.shape[0]):
        #ユーザuの全アーティストに対する嗜好度算出
        SC =numpy.dot(W[u , :] , H.T);

        #好みのアイテムの嗜好度のみ抽出
        IpScore = SC[ S[u,:] == 1 ];

        #未観測のアイテムの嗜好度のみ抽出
        ImScore = SC[ S[u,:] == 0 ];

        #(i , j)全ての組み合わせを計算し足しあわせる
        AA = numpy.tile(IpScore , (ImScore.shape[0] , 1))
        BB = numpy.tile(ImScore , (IpScore.shape[0] , 1))
        X=(AA-BB.T);
        lnSigX = numpy.log( 1/( 1 + numpy.exp(-X) ) );
        Bpropt = Bpropt + ( sum( sum(lnSigX) ) );

        #足し合わせた要素数も計算（最後にこれで割る）
        Zu = Zu + float(ImScore.shape[0] * IpScore.shape[0])

    Bpropt = (Bpropt - Lam * ((numpy.linalg.norm(W)**2)
                         + (numpy.linalg.norm(H)**2)))/Zu;
    return Bpropt

以上が今回組んだBPR MFのソースコードとなります。
BPROptの計算は全組み合わせを1度に計算するんじゃなくて、
$W$ と $H$ の更新時にBPROptも更新すれば速かったんじゃないかと、
いまさら思ってるんですが・・・まぁ、いいや

さて、いよいよ実験です。
まずパラメータ設定はk=50、alpha=0.7、Lam=0.1としました。
そして、このプログラムに上述の実験用データ突っ込み一晩寝かせます。
すると・・・

350万回繰り返したところで収束してました。
その結果ですが、先ほどの疑似ユーザ君3名に対する
推薦結果上位10アーティストは以下の通りです。

K-ON	gotham	joging
soft ballet	christian death	2 many dj's
ali project	coil	coil
the yellow monkey	the sisters of mercy	johann sebastian bach
perfume	the gathering	spiritualized
メリー	saltatio mortis	the fall
gackt	melvins	melvins
capsule	bauhaus	colour haze
電気グルーヴ	front 242	melt-banana
depapepe	ulver	guided by voices
lula	summoning	roxy music

遊びなんで定量的な評価はしませんが、
個人的には結構うまくいってるんじゃないかなと思います。というか、かなり満足。

以降、考察と言うより自分の趣味の観点からの感想です。

まず、ユーザK-ONに対する推薦結果ですが、
あの6組からcapsuleが推薦されたのは、個人的には結構嬉しかったりします。
イエモンが上位にきてるのも何となくわかる気が・・・ｗ
あと、20位くらいまで見るとsyrup16gやbuck-tickなんかもいて、なかなか面白い。
まぁ、これ見てると、テクノとかロックとかジャンル的な要素よりも、
日本のアーティストであるということが、推薦に最も寄与しているような気がします。
隠れ変数50個の中に、「日本」という属性を示す要素があるかもしれませんね。

つぎに、gothamくんの結果をみてみると、予想通りゴス、インダストリアル、
ブラックメタル系のアーティストが上位に来ています。うん、暗い。怖い。
このメンツが上位にいるなら、Throbbing GristleやPsychic TVなんかも
上位にいてもよさそうですが、残念ながら100位圏外でした。
20位まで見てみると、jesuやelectric wizardも挙がっていてました。
好みのアーティストとして挙げたのは70年代後半〜80年代の
いわゆるニューウェーブ期のバンドでしたが、
推薦結果には、それらに影響を受けた90〜00年代の
アーティストもいて、なかなか面白い。
年代、国関係なしに、その手の系統のアーティストが推薦されるといった感じでしょうか。

最後にユーザjoging。
これに関しては、本当にこの結果を利用して、
新たなジョギングソングを発掘しようと目論んでいました。
で、結果をみて最初に目についたのが、推薦結果の第3位。
まさかのバッハ。
確かにG線上のアリア聞きながら走るのも、
なかなか乙なものがあるのかもしれません。
他にも、monoという日本のポストロックバンドを知ることができたりと、
かなり使えそうです。