yt coffee

Study hard, play harder.

Matrix Factorization のための ALS アルゴリズムの数理と実装

  • 2020-01-06 01:09
  • #julia

Apache Spark の mllib には情報推薦でよく使われる Matrix Factorization (以下 MF)を計算するための ALS アルゴリズムが実装されています。 MF は直感的に理解しやすいモデルですし、 Spark を使えば大規模な user-item 行列にも適用できるのですが、実際に ALS がどういうアルゴリズムなのか知って使うのか、それともブラックボックスとして使うのかでは大きな違いがある、というのが個人的な考えです。

というわけでこの記事では Large-scale Parallel Collaborative Filtering for the Netflix Prize を参考に ALS が内部で何をしているのかを確認してから、実際に実装してみたいと思います。

MF の数理

MF ではレーティング行列2つの行列 PPQQ で近似します。

RP×QT=R^R \approx P \times Q^\mathrm{T} = \hat{R}

ユーザの総数を U|U|、 商品の総数を D|D| とすると PPQQ はそれぞれ U×K|U| \times K 行列と D×K|D| \times K 行列となります。このとき KK はモデルに与えられるパラメータで潜在因子の数を表します。

このときユーザ uu の商品 ii に対するレーティングの推定値 r^ui\hat{r}_{ui}

r^ui=puTqi\hat{r}_{ui} = \mathbf{p}_u^\mathrm{T} \mathbf{q}_i

で与えられます。 PPQQ を計算するには、 IIRR の中で実際に値が与えられている uuii の組を表すとして、二乗誤差に L1 正則化を加えた

L=(u,i)I(ruir^ui)2+αP+βQL = \sum_{(u, i) \in I} (r_{ui} - \hat{r}_{ui})^2 + \alpha \|P\| + \beta \|Q\|

を損失関数として最小化する方法があります。このアルゴリズムは 2006 年の Netflix Prize で Simon Funk によって提案されたので、一般的に Funk SVD と呼ばれています。

Funk SVD は単一のコンピュータ上で計算するには適していますが、スケーラビリティに問題がありました。今日ではテラバイトあるいはペタバイトという大きなのデータを扱うこともあるので、そのような規模のデータを一つのマシンで処理するのは大変です。

Alternating Least Square (ALS)

スケーラビリティを改善するために提案されたのが ALS です。 ALS では L2 正則化を行います。

L=(u,i)I(ruipuTqi)2+λ(unpupu2+inqiqi2)L = \sum_{(u,i) \in I} (r_{ui} - \mathbf{p}_u^\mathrm{T}\mathbf{q}_i)^2 + \lambda \Biggl(\sum_u n_{p_u}\|\mathbf{p}_u\|^2 + \sum_i n_{q_i}\|\mathbf{q}_i\|^2\Biggr)

ここで npun_{p_u}nqin_{q_i} はそれぞれユーザ uu と商品 ii のレーティングの数を表します。また IuI_u がユーザ uu が評価した商品の集合です。この損失関数を pukp_{uk} で偏微分してみると:

12Lpuk=iIu(puTqirui)qik+λnpupku=iIuqikqiTpuiIuqikrui+λnpupku\begin{aligned} \frac{1}{2} \frac{\partial L}{\partial p_{uk}} &= \sum_{i \in I_u}(\mathbf{p}_u^\mathrm{T}\mathbf{q}_i - r_{ui})q_{ik} + \lambda n_{p_u} p_{ku} \\ &= \sum_{i \in I_u} q_{ik} \mathbf{q}_i^\mathrm{T}\mathbf{p}_u - \sum_{i \in I_u} q_{ik}r_{ui} + \lambda n_{p_u} p_{ku} \end{aligned}

11 から KK までの潜在因子 kk を並べてベクトルとしてまとめると:

12Lpu=QIuQIuTpuQIuRT(u,In)+λnpupu=(QIuQIuT+λnpE)puQIuRT(u,In)\begin{aligned} \frac{1}{2} \frac{\partial L}{\partial \mathbf{p}_u} &= Q_{I_u} Q_{I_u}^\mathrm{T} \mathbf{p}_u - Q_{I_u} R^\mathrm{T} (u, I_n) + \lambda n_{p_u} \mathbf{p}_u \\ &= (Q_{I_u} Q_{I_u}^\mathrm{T} + \lambda n_p E) \mathbf{p}_u - Q_{I_u} R^\mathrm{T} (u, I_n) \end{aligned}

ここで QIuQ_{I_u}QQ の中からユーザ uu が評価した商品 IuI_u に関する行だけを抜き出したもので、RT(u,In)R^\mathrm{T} (u, I_n)RR からユーザ uu の行の中から商品 IuI_u に関するものだけを抜き出したベクトルを表します。

Lpu=0  pu=(QIuQIuT+λnpE)1QIuRT(u,In)  pu=Au1Vu\begin{aligned} \frac{\partial L}{\partial \mathbf{p}_u} &= 0 \\ \Rightarrow \; \mathbf{p}_u &= (Q_{I_u} Q_{I_u}^\mathrm{T} + \lambda n_p E)^{-1} Q_{I_u} R^\mathrm{T} (u, I_n) \\ \Rightarrow \; \mathbf{p}_u &= A_u^{-1}V_u \end{aligned}

このようにユーザ因子行列 PPQQRR だけに依存して更新できることが分かります。 qikq_{ik} で偏微分して整理していけば同様に商品因子行列 QQPPRR だけを使って更新できることが示されます。

ALS ではこの式を用いて PPQQ を交互に更新し、損失関数が変化しなくなるまで反復します。問題は AuA_u の逆行列が存在するのか定かでないことですが、この点についてもとの論文にも特に言及がなかったのでよくわかっていません。

分散処理

巨大なデータセットに対して ALS を実行する際にもっとも影響が大きな影響をもつのが RR です。 RRU×D|U| \times |D| 行列であり、ユーザ数と商品数が増えると猛烈な勢いで大きくなっていきますが、一般に一人のユーザが評価するのは商品のごく一部でしかないので、極めて疎に行列です。データセットが巨大になったとき RR をそのままメモリにもつことは困難です。

前述のように ALS では PPQQRR だけを使って更新します。そのため、 PPQQ の計算は各反復において互いに独立して並列に行うことができます。

さらに RR は更新しているユーザが評価した商品の評価値だけを取り出せばいいので、 CSC のように列方向に圧縮された形式で保持することで、データサイズを小さくしつつ効率的に計算をすることが可能です。

ただこの記事では単純に大変であることに加えて、単一の Jupyter Notebook 上で実装する方法がよく分からないので、分散処理は実装しません。

ALS を実装する

ALS を実装する際にもっとも紛らわしいのが QIuQ_{I_u} を求める部分でしょう。

仮にレーティング行列として次の RR が与えられたとします。

R =  [
    0   4.5 2.0 0;
    4.0 0   3.5 0;
    0   5.0 0   2.0;
    0   3.5 4.0 1.0
]
4×4 Array{Float64,2}:
 0.0  4.5  2.0  0.0
 4.0  0.0  3.5  0.0
 0.0  5.0  0.0  2.0
 0.0  3.5  4.0  1.0

行列の行が上からユーザ 1, 2, 3, 4 に対応し、列が左から商品 1, 2, 3, 4 に対応しているとすると、この行列はユーザ 1 は商品 2 を 4.5、商品 3 を 2.0 と評価し、商品 1 と 4 を評価していないとみなすことができます。なので K=2K=2

Q=(q11q12q21q22q31q32q41q42)Q = \left( \begin{array}{ccc} q_{11} & q_{12} \\ q_{21} & q_{22} \\ q_{31} & q_{32} \\ q_{41} & q_{42} \end{array} \right)

とすると

QI1=(01000010)Q=(q21q22q31q32)Q_{I_1} = \left( \begin{array}{ccc} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{array} \right) Q = \left( \begin{array}{ccc} q_{21} & q_{22} \\ q_{31} & q_{32} \end{array} \right)

です。この QQ の左からかけているフィルターとなる行列を作る関数を用意します。

"""
ユーザあるいは商品のレーティングベクトルを受け取って、フィルターとなる行列を返す
"""
function filtermatrix(rating::Array)
    result = zeros(Int, (count(rating .!= 0), length(rating)))
    counter = 1
    for (index, value) in enumerate(rating)
        if value != 0
            result[counter, index] = 1
            counter += 1
        end
    end
    result
end
filtermatrix

ユーザ 1 のレーティングを入力として与えると、たしかにさきほど QQ にかけたものと同じ行列が返ってきます。qi\mathbf{q}_i を更新する時に必要になる PIiP_{I_i} も同様に求めることができます:

filtermatrix(R[1:1, :])
2×4 Array{Int64,2}:
 0  1  0  0
 0  0  1  0

本筋からは離れますが、 ALS の中で使うために MSE を計算する関数を用意しておきます:

"""
ベクトル x と y の平均二乗誤差を計算する

x: ベクトル
y: ベクトル
"""
function mse(x::Vector, y::Vector)::Float64
    @assert length(x) == length(y)
    sum((x .- y).^2) / length(x)
end
mse

では ALS を実装します。さきほどの更新式をそのままコードにおとしていくだけです。

using LinearAlgebra

"""
R: レーティング行列
rank: 潜在因子の数
iter: 反復回数
λ: 正規化パラメータ
"""
function als(
        R::Matrix;
        rank::Int=10,
        iter::Int=10,
        λ::Float64=0.001
    )
    usersnum = size(R, 1) # ユーザ数
    itemsnum = size(R, 2)  # 商品数
    
    P = rand(usersnum, rank)  # ユーザ因子行列
    Q = rand(itemsnum, rank)  # 商品因子行列
    
    # 各ユーザと商品の評価数
    np = [count(R[u:u, :] .> 0) for u in 1:usersnum]
    nq = [count(R[:, i:i] .> 0) for i in 1:itemsnum]

    # 各ユーザと商品のフィルター行列
    Iu = [filtermatrix(R[u:u, :]) for u in 1:usersnum]
    Ii = [filtermatrix(R[:, i:i]) for i in 1:itemsnum]
    
    x = R[R .!= 0] # レーティング行列のうち 0 ではないもの
    losses = []  # 各イテレーションごとの mse の値を入れる
    
    for _ in 1:iter
        # ユーザ因子を更新する
        for u in 1:usersnum
            QIu = (Iu[u] * Q)'
            Au = QIu * QIu'  + λ * np[u] * I
            Vu = QIu * Iu[u] * (R[u:u, :])'
            P[u:u, :] = inv(Au) * Vu
        end
        
        # 商品因子を更新する
        for i in 1:itemsnum
            PIi = (Ii[i] * P)'
            Ai = PIi * PIi' + λ * nq[i] * I
            Vi = PIi * Ii[i] * R[:, i:i]
            Q[i:i, :] = inv(Ai) * Vi
        end
        
        # 推定された R
        = P * Q'
        #  R̂ のうちレーティング行列で 0 ではない要素に対応するもの
        y =[R .!= 0]
        push!(losses, mse(x, y))
    end
    
    P, Q, losses
end
als

動作確認

実際に計算してみます:

using Random
Random.seed!(255)
P, Q, losses = als(R)
= P * Q'
4×4 Array{Float64,2}:
 1.95614  4.49856  2.00126  1.93523 
 3.99916  2.49496  3.5      0.577751
 2.59633  5.00034  2.96279  1.99935 
 2.92678  3.49851  3.99614  1.00025 

RRR^\hat{R} から II に含まれるものを取り出してみると、かなり近い値が並んでいることが分かります:

println("R = ", R[R .!= 0])
println("R̂ = ", round.([R .!= 0], digits=3))
R = [4.0, 4.5, 5.0, 3.5, 2.0, 3.5, 4.0, 2.0, 1.0]
R̂ = [3.999, 4.499, 5.0, 3.499, 2.001, 3.5, 3.996, 1.999, 1.0]

誤差の遷移をグラフに表示してみると、反復すると急激に値が小さくなっていくことが分かります。今回使ったレーティング行列はとても小さいので数回の反復で収束しています:

using Plots
Plots.gr()
plot(1:10, losses, label="Losses", yscale=:log10)
24681010-5.50 10-5.25 10-5.00 10-4.75 Losses

ユーザ数と商品数ともに 100 くらいのもっと大きなデータを使ってみます:

R2 = zeros(100, 100)
for u in 1:100, p in 1:100
    if rand() < 0.1
        R2[u, p] = rand() * 5
    end
end

収束が緩やかになっていることが分かります:

iter = 100
_, _, losses2 = als(R2, iter=iter)
Plots.gr()
plot(1:iter, losses2, label="Losses", yscale=:log10)
025507510010-5 10-4 10-3 10-2 Losses

おわりに

協調フィルタリングで広く使われる MF の実装の一つである ALS の数学的な背景とその実装を確認しました。実際に MF する際は Apache Spark に実装された ALS を用いることが多いと思いますが、こうして中身が分かると少しだけ安心感も増す感じがします。

ALS は MF のためのものではありますが、交互に二乗誤差を修正する、という考え方自体は他にも適用できるものです。実際、 Factorization Marchines(FM)の論文では FM が ALS で解けると書かれています。

実はこのブログは当初 FM の勉強のために書こうとしたものなのです。しかしいかんせん、長くなってきたのでここで一旦終わります。機会があれば FM を ALS (または MCMC)で解いてみたいと思います。