読者です 読者をやめる 読者になる 読者になる

関係データ学習の実装 MovieLensの行列分解

概要

関係データ学習の後半のうち,行列分解を実装した. lossが下がっているので,たぶん合っていると思うが,分解したあとのベクトルを見るとよくわからない.

はじめに

以前こんなものをやった. ksknw.hatenablog.com

前半部分は終わったので(SBMのことは置いといて),後半部分の実装をやる. 今回は行列分解を実装する.

www.kspub.co.jp

MovieLensデータ

MovieLens 100K Dataset | GroupLens

これといって分解したい行列データもなかったので,本に書いてあったデータを使うことにする. MovieLensにはユーザが映画を評価したデータが含まれている.
データを作成した時期によってデータのサイズが違うらしいが,今回は一番小さいサイズのデータを使う.
このデータには943人のユーザの映画の評価(1から5)が含まれている.
映画は1682本で各ユーザはこれらのうち少なくとも5本を評価している.
とりあえずデータを見てプロットする.

%matplotlib inline
import pylab as plt
import seaborn
import numpy as np
import pandas as pd
data = pd.read_csv("./ml-100k/u.data", names=["uid", "mid", "rating", "timestamp"], sep="\t")
data
uid mid rating timestamp
0 196 242 3 881250949
1 186 302 3 891717742
2 22 377 1 878887116
3 244 51 2 880606923
4 166 346 1 886397596
29 160 234 5 876861185
... ... ... ... ...
99995 880 476 3 880175444
99996 716 204 5 879795543
99997 276 1090 1 874795795
99998 13 225 2 882399156
99999 12 203 3 879959583

100000 rows × 4 columns

X = np.zeros((data["uid"].max(), data["mid"].max()))
for i,item in data.iterrows():
    X[item["uid"]-1, item["mid"]-1] = item["rating"]

f, ax = plt.subplots(figsize=(7, 7))

import numpy.ma as ma
masked = ma.masked_where((X==0),X)


cmap = plt.cm.jet
cmap.set_bad('white',1.)
plt.pcolormesh(masked, cmap=cmap)
plt.xlabel("movie id")
plt.ylabel("user id")
plt.show()

f, ax = plt.subplots(figsize=(7, 7))

plt.pcolormesh(masked[:100, :100], cmap=cmap)
plt.xlabel("movie id")
plt.ylabel("user id")
plt.show()

f:id:ksknw:20170305215433p:plain

f:id:ksknw:20170305215506p:plain

人気の映画やよく映画を見ている人などがなんとなくわかる. 多くの人が好む映画,やたら低評価ばかりつけるユーザなど色々いる. なんらかソートされていることもわかる.
また大半は欠損データである.

まずは,(練習なので)欠損データは0とレーティングされているとして行列分解する.
その後,欠損値を除外したり,補完をするアルゴリズムを実装する.
また,本に載っているアルゴリズムでは収束判定をしていたが省略する.

import copy 
original_X = copy.deepcopy(X) 

# X = X+(X==0)*3
#X = X[:100, :100]
f, ax = plt.subplots(figsize=(7, 7))
plt.pcolormesh(X[:100, :100], cmap=cmap)

f:id:ksknw:20170305215522p:plain

1次交互勾配降下法

正則化つきの行列分解は特異値分解では求まらないので,勾配法を用いる必要がある.

それぞれ5つのベクトルでユーザ,映画が表現できると考えて,行列分解する.
動けばいいやの精神で書いているので,行列の計算はあれな感じだけど許してほしい. というかnumpy.matとarray関連の正しい使い方がわからない.

1次交互勾配法は不安定で,learning_rateを少しあげたら発散した.

nb_epoch = 100
reg = 0.01

R = 5
learning_rate = 0.001
I = X.shape[0]
J = X.shape[1]
U = np.mat(np.random.random((I, R)))
V = np.mat(np.random.random((J, R)))
dU = lambda U,V: -X*V + U * (V.T*V + reg*np.mat(np.identity(R)))
dV = lambda U,V: -X.T*U + V * (U.T*U + reg*np.mat(np.identity(R)))
error = lambda U,V: 0.5*np.sum(np.array(X-U*V.T)**2) + 0.5*reg*np.sum( np.sum(np.array(U)**2,axis=0) + np.sum(np.array(V)**2,axis=0))
err_hist = []
for i in range(nb_epoch):
    e = error(U,V)
    U -= learning_rate * dU(U,V)
    V -= learning_rate * dV(U,V)
    err_hist.append(e)
def hist_plot(err_hist):
    plt.plot(err_hist)
    plt.text(len(err_hist) -5, err_hist[-1]+500 ,"%.1f"%err_hist[-1])
hist_plot(err_hist)

f:id:ksknw:20170305215543p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)
(10.693971705841749, -2.4352853936451466)
def compare_plot(X, predict, cmap=plt.cm.Blues, vmin=None, vmax=None):
    f, ax = plt.subplots(figsize=(14, 7))
    plt.subplot(121)
    plt.pcolormesh(X[:100, :100], cmap=cmap, vmin=vmin, vmax=vmax)
    plt.xlabel("movie id")
    plt.ylabel("user id")
    
    plt.subplot(122)
    plt.pcolormesh(predict[:100, :100], cmap=cmap, vmin=vmin, vmax=vmax)
    plt.xlabel("movie id")
    plt.ylabel("user id")

    plt.show()
compare_plot(X, predict)

f:id:ksknw:20170305215555p:plain

最大値,最小値がずれているのでJetだと見にくかったので,青だけでプロットしている.
なんとなく出来ている気がする.

擬似2次交互勾配法

1次勾配だけだと収束が不安定で学習率をあまり大きくすることができず,収束が遅い.
2次勾配を使うことで収束を早くすることができる.

learning_rate = 1.

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.random.random((I, R)))
V = np.mat(np.random.random((J, R)))
err_hist2 = []

for i in range(nb_epoch):
    e = error(U,V)
    U = (1-learning_rate)*U + learning_rate * X*V*((V.T*V+reg*np.mat(np.identity(R))).I)
    V = (1-learning_rate)*V + learning_rate * X.T*U*((U.T*U+reg*np.mat(np.identity(R))).I)
    err_hist2.append(e)

hist_plot(err_hist)
hist_plot(err_hist2)

f:id:ksknw:20170305215610p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)
(10.414641991236321, -2.2410733139172789)
compare_plot(X, predict)

f:id:ksknw:20170305215631p:plain

学習率を1にしても発散せず,1次勾配(青)よりも2次勾配(緑)を利用したものは早く収束するようになった.

制約付きの最適化

ここまでの実装では,推定した評価値が負の値になっている. U,Vに制約をつけることで0以上の値にすることができる. (評価値は1から5までなので本当はもう少し制約しないといけないがどうするのかはわからない.)

具体的には,負の値をとった要素を0に射影しながら最適化する.

learning_rate = 1.

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.random.random((I, R)))
V = np.mat(np.random.random((J, R)))
err_hist_con = []

for i in range(nb_epoch):
    e = error(U,V)
    U = (1-learning_rate)*U + learning_rate * X*V*((V.T*V+reg*np.mat(np.identity(R))).I)
    U = np.maximum(U, 0, U)
    V = (1-learning_rate)*V + learning_rate * X.T*U*((U.T*U+reg*np.mat(np.identity(R))).I)
    V = np.maximum(V, 0, V)

    err_hist_con.append(e)

hist_plot(err_hist_con)

f:id:ksknw:20170305215653p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)
(8.9972633401464464, 0.0)

ちゃんと0以上に制約されている.

compare_plot(X, predict)

f:id:ksknw:20170305215708p:plain

欠損値をちゃんとやる

ここまでは欠損値を全て0だと仮定して実装した.

以下では欠損値を除外して行列分解する方法と補完する方法を行う.

除外する方法

lossの計算およびU,Vの勾配の計算のときに,欠損した値は利用しない.
実装が楽なので1次勾配を利用する方法で実装した.

M = np.array((X!=0)*1)

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.abs(np.random.random((I, R))))
V = np.mat(np.abs(np.random.random((J, R))))
dU = lambda U,V: -(np.array(X - U*V.T)*M)*V   + reg*U
dV = lambda U,V: -(np.array(X - U*V.T)*M).T*U + reg*V
error = lambda U,V: 0.5*np.sum((np.array(X-U*V.T)*M)**2) + 0.5*reg*np.sum( np.sum(np.array(U)**2,axis=0) + np.sum(np.array(V)**2,axis=0))
err_hist = []
for i in range(nb_epoch):
    e = error(U,V)
    U = (1-learning_rate)*U + learning_rate * X*V*((V.T*V+reg*np.mat(np.identity(R))).I)
    V = (1-learning_rate)*V + learning_rate * X.T*U*((U.T*U+reg*np.mat(np.identity(R))).I)
    err_hist.append(e)
hist_plot(err_hist)

f:id:ksknw:20170305215722p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)
(10.414841938424392, -2.2413625198086238)
import numpy.ma as ma
masked_predict = ma.masked_where((X==0),predict)

compare_plot(masked[:100, :100], masked_predict, cmap=cmap, vmin=1, vmax=5)
compare_plot(masked[:100, :100], predict, cmap=cmap, vmin=1, vmax=5)

f:id:ksknw:20170305215738p:plain

f:id:ksknw:20170305215755p:plain

全体的に評価が低めな感じになってしまっている. また,欠損していた部分も低めの値で埋まっている感じ.

補完する方法

こちらの方法では欠損値を毎回,予測値で補完する. 欠損値を除外する方法と等価らしいので,こちらは2次勾配を使って,制約もつける方法で実装する.

M_hat = np.array((X==0)*1)

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.abs(np.random.random((I, R))))
V = np.mat(np.abs(np.random.random((J, R))))
error = lambda U,V: 0.5*np.sum((np.array(X-U*V.T)*M)**2) + 0.5*reg*np.sum( np.sum(np.array(U)**2,axis=0) + np.sum(np.array(V)**2,axis=0))

err_hist_comp = []
X_hat = X + np.mat(np.array(U*V.T)*M_hat)

for i in range(nb_epoch):
    e = error(U,V)
    U = (1-learning_rate)*U + learning_rate * X_hat*V*((V.T*V+reg*np.mat(np.identity(R))).I)
    U = np.maximum(U, 0, U)
    X_hat = X + np.mat(np.array(U*V.T)*M_hat)
    V = (1-learning_rate)*V + learning_rate * X_hat.T*U*((U.T*U+reg*np.mat(np.identity(R))).I)
    V = np.maximum(V, 0, V)
    X_hat = X + np.mat(np.array(U*V.T)*M_hat)
    err_hist_comp.append(e)
hist_plot(err_hist_comp)

f:id:ksknw:20170305215813p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)
(7.7494405681324903, 0.0082254500224862269)
import numpy.ma as ma
masked_predict = ma.masked_where((X==0),predict)

compare_plot(masked[:100, :100], masked_predict, cmap=cmap, vmin=1, vmax=5)
compare_plot(masked[:100, :100], predict, cmap=cmap, vmin=1, vmax=5)

f:id:ksknw:20170305215827p:plain

f:id:ksknw:20170305215841p:plain

だいたいできているような気がする.

ベクトルをみる

行列分解で得られたU,Vでは近いユーザもしくは映画が近いベクトルで表現される.
MovieLensには評価だけでなく,どんな映画かやどんなユーザかといった情報も含まれている.
これらとU,Vを比較することで,正しい行列分解ができたかを評価する.

普段あまり映画を見ないのと,古い映画のデータしかないので,Toy StoryとGolden Eyeぐらいしかわからなかった.
ユーザの属性ならなんとなくわかるので,そちらを使う. ユーザの情報は以下.

user_info = pd.read_csv('ml-100k/u.user', sep='|',
                         encoding='latin-1',header=None)
user_info
0 1 2 3 4
0 1 24 M technician 85711
1 2 53 F other 94043
2 3 23 M writer 32067
3 4 24 M technician 43537
... ... ... ... ... ...
937 938 38 F technician 55038
938 939 26 F student 33319
939 940 32 M administrator 02215
940 941 20 M student 97229
941 942 48 F librarian 78209
942 943 22 M student 77841

943 rows × 5 columns

職業がわかりやすそうなので,プロットしてみる.

from sklearn.manifold import TSNE
tsne = TSNE()
tsned = tsne.fit_transform(U) 
plt.scatter(tsned[:,0], tsned[:,1], c=(user_info[3]=="scientist"), cmap=plt.cm.Accent, marker=".")
plt.show()
plt.scatter(tsned[:,0], tsned[:,1], c=(user_info[3]=="engineer"), cmap=plt.cm.Accent, marker=".")
plt.show()
plt.scatter(tsned[:,0], tsned[:,1], c=(user_info[3]=="student"), cmap=plt.cm.Accent, marker=".")
plt.show()

f:id:ksknw:20170305215854p:plain

f:id:ksknw:20170305215903p:plain

f:id:ksknw:20170305215911p:plain

そうでもない感じだったので,うまくいっているのかよくわからない.

まとめ

勉強のために行列分解の実装をした. 収束はしたけど,うまくいっているのかはよくわからない.

参考