2008-10-12

Recommender Systems 2007(http://recsys.acm.org/2007/)で発表された論文である，Bhaskar Mehta, Thomas Hofmann, and Wolfgang Nejdl, Robust Collaborative Filtering, In Proceedings of the 1st ACM Conference on Recommender Systems, ACM Press, October 2007, pp. 49–56. を読んだメモです．この論文では，ある種の攻撃に耐えられるような，頑強な協調フィルタリングの手法を提案していますが，その説明は後日行うことにして，今回は，関連研究に挙げられていた，特異値分解を用いたレコメンデーションアルゴリズムについて説明を行いたいと思います．

特異値分解

ある任意のm x n行列 $\mathbf{X}$ は下記のように分解可能であることが知られています．
$\mathbf{X} = \mathbf{U} \cdot \mathbf{\Sigma} \cdot \mathbf{V}^T$
ただし， $x_{i, j} \in \mathbb{R}$ であるとき， $\mathbf{U}$ と $\mathbf{V}$ はサイズm x mとn x nの正規直交行列になります．また， $\mathbf{\Sigma}$ はm x nの行列であり，対角成分 $\sigma_{i, i}$ は $\mathbf{X}$ のi番目に大きい固有値となり，対角成分以外の要素は0となります．

今， $\mathbf{u}_i$ および $\mathbf{v}_j$ を $\mathbf{U}$ と $\mathbf{V}$ のiとj番目の列ベクトルであるとすると， $\sigma_{1, 1}$ に対応する特異ベクトルは， $\mathbf{u}_1$ と $\mathbf{v}_1$ になります．ところで，この $\mathbf{u}_1$ および $\mathbf{v}_1$ のみで構成された行列 $\mathbf{U}' = [\mathbf{u}_1 \hspace{10px} 0 \hspace{10px} \cdots \hspace{10px} 0 ]$ と $\mathbf{V}' = [\mathbf{v}_1 \hspace{10px} 0 \hspace{10px} \cdots \hspace{10px} 0 ]$ の積は， $\mathbf{X}$ の第一主成分となります．すると， $\mathbf{X}$ から，この第一主成分を引いた行列，
$\mathbf{X}_{\mathrm{rem}} = \mathbf{X} - \mathbf{U}' \cdot \sigma_{1, 1} \cdot \mathbf{V}'^T$
の第一主成分は， $\mathbf{X}$ の第二主成分となることが明らかです．この特徴は，巨大な疎行列の固有値と固有ベクトルを求めるのに利用されており，Hotteling's Deflation Methodと呼ばれているそうです．

Hebbian learning

Hebbian learningとはニューラルネットワークで用いられる学習則の一つであり，ここでは，ニューラルネットワーク的な手法で，発見的に特異値分解を行うする方法を説明します．

ここで， $\mathbf{X}$ の第一主成分を表す特異ベクトルをそれぞれ， $\mathbf{u}$ ， $\mathbf{v}$ と表したとします．今，何らかの要因で， $\mathbf{X}$ の固有値 $\sigma_{1, 1}$ が， $\mathbf{u}$ ， $\mathbf{v}$ に吸収されたとすると， $\mathbf{X}$ の第一主成分のみで表されたi行j列目の要素は，
$\hat{x}_{i, j} = \hat{u}_i \cdot \hat{v}_j$
と表されます．ただし， $\hat{u}_i$ および， $\hat{v}_j$ はそれぞれ， $\sigma_{1, 1}$ を吸収した後の $\mathbf{u}$ ， $\mathbf{v}$ のi番目およびj番目の要素です．

いま，この $\hat{u}_i$ と $\hat{v}_j$ を何らかの方法で予測を行った値だと仮定すると，その残差は以下のように定義できます．
$r(x_{i, j}) = x_{i, j} - \hat{x}_{i, j} = x_{i, j} - \hat{u}_i \cdot \hat{v}_j$

Hebbian learningを用いた手法では，推測した特異行列から元の行列を推測し，その推測した行列と元の行列との残差を用いて，特異行列のアップデートを行い学習していきます．その学習規則は下記のようになります．
$\Delta \hat{u}_i = \lambda \cdot \hat{v}_j \cdot r(x_{i, j})$
$\Delta \hat{v}_i = \lambda \cdot \hat{u}_j \cdot r(x_{i, j})$
ここで， $\lambda$ は学習率となり，0.1 〜 0.001ぐらいの値が良いようです．第一主成分に対する推測が終わった後は，
$x_{i, j} \leftarrow x_{i, j} \hat{x}_{i, j} - \hat{u}_i \cdot \hat{v}_i$
として，第二主成分に対する予測を同じように行います．ここまでがオリジナルのHebbian learningを用いた手法ですが，Simon Frankはさらにこれを改良しています．

Simon's mehotd

Simonとは一時期，Netflix Prizeにチャレンジしていた人で，Netflix Prizeで結構上位の方に食い込んでいたそうです．彼は，そのアルゴリズムをブログで公開しています（http://sifter.org/~simon/journal/20061211.html）

Simonの方法は，Hebbian learningと若干学習規則が違い，以下のようになります．
$\Delta \hat{u}_i = \lambda(\hat{v}_j \cdot r(x_{i, j}) - k \cdot \hat{u}_i)$
$\Delta \hat{v}_i = \lambda(\hat{u}_j \cdot r(x_{i, j}) - k \cdot \hat{v}_i)$
$k \cdot \hat{u}_i$ および $k \cdot \hat{v}_i$ は，過剰適合を避けるように導入した罰則項であり，100以上の特異ベクトルを求める場合， $k = 0.02$ が良いそうです．

Simonはさらに，予測された値をNetflixのレーティングの範囲である1〜5に収めるために，値のクリップを行っています．クリップには単純に1〜5に収める方法と，シグモイド関数などを用いる方法が提案されていますが，ここではその説明は割愛します．

Pythonでのプログラム

Pythonでプログラムを書いてみると以下のようになります．slope one(http://d.hatena.ne.jp/ytakano/20081005/1223201692)やitem based(http://d.hatena.ne.jp/ytakano/20081002/1222970856)のプログラムと比較してみると分かりますが，ループの回数が格段に少なくなっているため，非常に高速に計算できます．

このプログラムではアイテム数が少ないので，learning rateを0.1としていますが，ユーザ数，アイテム数が多ければ，0.001とかの方が良いです．

# this program implements SVD based recommendation algorithms
#
# see section 3 of
# Bhaskar Mehta, Thomas Hofmann, and Wolfgang Nejdl, Robust Collaborative Filtering,
# In Proceedings of the 1st ACM Conference on Recommender Systems, ACM Press, October 2007,
# pp. 49-56
# and
# Simon's blog (http://sifter.org/~simon/journal/20061211.html)
from copy import deepcopy

class hebbian:
    def __init__(self, users, items):
        self.lrate = 0.1 # learning rate
        self.comps = 40  # the number of components

        self._users = deepcopy(users)
        self._items = items

        self._item2idx = {}
        i = 0
        for item in items:
            self._item2idx[item] = i
            i += 1

    def _init_uv(self):
        t1 = 0
        t2 = 0
        for user in self._users:
            for val in user.values():
                t1 += val
                t2 += 1

        ave = t1 / float(t2)
        p = (ave / self.comps) ** 0.5
        ulen = len(self._users)
        ilen = len(self._items)
        self._u = [[p for i in xrange(self.comps)] for j in xrange(ulen)]
        self._v = [[p for i in xrange(self.comps)] for j in xrange(ilen)]

    def _svd(self):
        for comp in xrange(self.comps):
            for i in xrange(len(self._users)):
                for item in self._users[i].keys():
                    j = self._item2idx[item]
                    u = self._u[i][comp]
                    v = self._v[j][comp]

                    x = self._users[i][item] # real
                    e = u * v # estimation
                    r = x - e # residual
                    self._u[i][comp] += self.lrate * v * r
                    self._v[j][comp] += self.lrate * u * r

            # remove
            for i in xrange(len(self._users)):
                for item in self._users[i].keys():
                    j = self._item2idx[item]
                    if comp > 0:
                        u = self._u[i][comp - 1]
                        v = self._v[j][comp - 1]
                        self._users[i][item] -= u * v

    def _predict(self, user, item):
        p = 0.0
        item = self._item2idx[item]
        for comp in xrange(self.comps):
            p += self._u[user][comp] * self._v[item][comp]

        return p

    def recommends(self, user):
        self._init_uv()
        self._svd()

        items = [x for x in self._items if not x in self._users[user]]
        result = []
        for item in items:
            p = self._predict(user, item)
            result.append((item, p))

        result.sort(lambda x, y: cmp(y[1], x[1]))
        return result


class simon(hebbian):
    def __init__(self, user, item):
        self.k = 0.02
        hebbian.__init__(self, user, item)

    def _svd(self):
        for comp in xrange(self.comps):
            for i in xrange(len(self._users)):
                for item in self._users[i].keys():
                    j = self._item2idx[item]
                    u = self._u[i][comp]
                    v = self._v[j][comp]

                    x = self._users[i][item] # real
                    e = u * v # estimation
                    r = x - e # residual
                    self._u[i][comp] += self.lrate * (v * r - self.k * u)
                    self._v[j][comp] += self.lrate * (u * r - self.k * v)

            # remove
            for i in xrange(len(self._users)):
                for item in self._users[i].keys():
                    j = self._item2idx[item]
                    if comp > 0:
                        u = self._u[i][comp - 1]
                        v = self._v[j][comp - 1]
                        self._users[i][item] -= u * v


users = [{'A': 4, 'B': 5, 'C': 2, 'D': 4,         'F': 5}, # user 0
         {'A': 2,         'C': 3, 'D': 4, 'E': 3        }, # user 1
         {'A': 1, 'B': 4,         'D': 5, 'E': 3, 'F': 4}, # user 2
         {        'B': 5,                 'E': 2, 'F': 4}, # user 3
         {        'B': 3, 'C': 1, 'D': 3,         'F': 3}] # user 4

items = ['A', 'B', 'C', 'D', 'E', 'F']

svd = hebbian(users, items)
print svd.recommends(3)

svd = simon(users, items)
print svd.recommends(3)

# outputs are
# [('D', 5.5970040690706924), ('A', 3.7606889484910573), ('C', 3.6752331888223218)]
# [('D', 5.5732517060635836), ('A', 3.7456618851422818), ('C', 3.6607220368100424)]