Python: scipy.sparseで疎行列計算する

疎行列計算用のモジュール scipy.sparse について解説します。

https://docs.scipy.org/doc/scipy/reference/sparse.html

疎行列クラス

疎行列とは要素のほとんど(90%以上)が0で構成される行列です。
レコメンドやソーシャルグラフなどの分野ではしばしば現れ、総じて巨大な行列(数千×数千〜)のため、そのままではCPUキャッシュやメモリに乗り切らずに非効率な計算になりやすく、扱いが難しいデータです。

こうした疎行列を効率的に扱うために、Pythonではscipy.sparseがよく使われます。

scipy.sparseで提供されている疎行列クラスはいくつかあります。
詳細は公式ドキュメントを参照ですが、代表的な3つのクラス(lil_matrix, csr_matrix, csc_matrix)について紹介します。

https://docs.scipy.org/doc/scipy/reference/sparse.html#sparse-matrix-classes

lil_matrix

lil_matrixは、リストやnumpyのクラス(ndarrayやmatrixなど)と相互変換するために使われることが多いクラスです。行単位のリンクトリストで実装されてます。

• lil_matrixのインスタンスをprintで表示すると、非ゼロ要素のインデックスと値だけしか持っていないことがわかります
• numpy.ndarrayに変換する時はtoarray()、numpy.matrixに変換する時はtodense()を使います
• csr_matrixに変換する時はtocsr()、csc_matrixに変換する時はtocsc()を使います

import numpy as np
from scipy import sparse


# リスト、numpy.ndarray、numpy.matrixのデータ
matrix1 = [[1, 0, 3], [0, 5, 0]]
matrix2 = np.array([[1, 0, 3], [0, 5, 0]])
matrix3 = np.matrix([[1, 0, 3], [0, 5, 0]])

# lil_matrixの生成
lil1 = sparse.lil_matrix(matrix1, dtype=np.float64) # リストから生成(dtypeで要素の型を指定できる)
lil2 = sparse.lil_matrix(matrix2) # numpy.ndarrayから生成
lil3 = sparse.lil_matrix(matrix3) # numpy.matrixから生成

# lil_matrixの表示
print(type(lil1), lil1)
# <class 'scipy.sparse.lil.lil_matrix'>
#  (0, 0)  1.0
#  (0, 2)  3.0
#  (1, 1)  5.0

# 要素にアクセス
print(lil1[0, 0], lil1[0, 1], lil1[1, 1]) # 1.0 0.0 5.0
# スライスで1〜2列目にアクセス
print(lil1[:, 0:2])
#  (0, 0)  1.0
#  (1, 1)  5.0

# lil_matrixからnumpy.ndarrayへの変換
lil1_ndarray = lil1.toarray()
print(type(lil1_ndarray), lil1_ndarray)
# <class 'numpy.ndarray'> [[ 1.  0.  3.] [ 0.  5.  0.]]

# lil_matrixからnumpy.matrixへの変換
lil1_matrix = lil1.todense()
print(type(lil1_matrix), lil1_matrix)
# <class 'numpy.matrixlib.defmatrix.matrix'> [[ 1.  0.  3.] [ 0.  5.  0.]]

# lil_matrixからcsr_matrixへの変換
lil1_csr = lil1.tocsr()
print(type(lil1_csr)) # <class 'scipy.sparse.csr.csr_matrix'>

# lil_matrixからcsc_matrixへの変換
lil1_csc = lil1.tocsc()
print(type(lil1_csc)) # <class 'scipy.sparse.csc.csc_matrix'>

csr_matrix

csr_matrixは、行ごとに圧縮されており、行単位なら高速にアクセスできます。また、csc_matrix同士の加減算やベクトル積・行列積も高速です。

# csr_matrixの生成
csr1 = sparse.csr_matrix(matrix1, dtype=np.float64) # リストから生成(dtypeで要素の型を指定できる)
csr2 = sparse.csr_matrix(matrix2) # numpy.ndarrayから生成
csr3 = sparse.csr_matrix(matrix3) # numpy.matrixから生成
csr4 = sparse.csr_matrix(lil1) # lil_matrixから生成

# csr_matrixの表示
print(type(csr1), csr1)
# <class 'scipy.sparse.csr.csr_matrix'>   
#  (0, 0)  1.0
#  (0, 2)  3.0
#  (1, 1)  5.0

# 要素にアクセス
print(csr1[0, 0], csr1[0, 1], csr1[1, 1]) # 1.0 0.0 5.0
# スライスで1〜2列目にアクセス
print(csr1[:, 0:2])
#  (0, 0)  1.0
#  (1, 1)  5.0

# csc_matrixからlil_matrixへの変換
csr1_lil = csr1.tolil()
print(type(csr1_lil), csr1_lil)
# <class 'scipy.sparse.lil.lil_matrix'> 
#  (0, 0)  1.0
#  (0, 2)  3.0
#  (1, 1)  5.0

csc_matrix

csc_matrixは、列ごとに圧縮されており、列単位なら高速にアクセスできます。こちらも、csr_matrix同士の加減算やベクトル積・行列積も高速です。

• toarray()でnumpy.ndarray、todense()でnumpy.matrix、tolil()でlil_matrixへそれぞれ変換できます
• lil_matrixやcsr_matrixと異なり、csc_matrixクラスは列順で表示されていることがわかります((0, 0)要素の次に(1, 1)要素が表示されている)

# csc_matrixの生成
csc1 = sparse.csc_matrix(matrix1, dtype=np.float64) # リストから生成(dtypeで要素の型を指定できる)
csc2 = sparse.csc_matrix(matrix2) # numpy.ndarrayから生成
csc3 = sparse.csc_matrix(matrix3) # numpy.matrixから生成
csc4 = sparse.csc_matrix(lil1) # lil_matrixから生成

# csc_matrixの表示
print(type(csc1), csc1)
# <class 'scipy.sparse.csc.csc_matrix'>
#  (0, 0)  1.0
#  (1, 1)  5.0
#  (0, 2)  3.0

# 要素にアクセス
print(csc1[0, 0], csc1[0, 1], csc1[1, 1]) # 1.0 0.0 5.0
# スライスで1〜2列目にアクセス
print(csc1[:, 0:2])
#  (0, 0)  1.0
#  (1, 1)  5.0

# csc_matrixからlil_matrixへの変換
csc1_lil = csc1.tolil()
print(type(csc1_lil), csc1_lil)
# <class 'scipy.sparse.lil.lil_matrix'>
#  (0, 0)  1.0
#  (0, 2)  3.0
#  (1, 1)  5.0

計算

scipy.sparseの計算ですが、内積とアダマール積(要素ごとの積)を以外は、numpy.ndarrayと概ね同じです。

• lil_matrix、csr_matrix、csc_matrixを混合して計算できます
  • 演算子の前(オペランド)の型で評価されます
• numpyのndarrayやmatrixと混合して計算することもできます
  • その場合、matrixで評価されます
  • 性能的な恩恵も受けられにくくなります

加減算

演算子+と-で疎行列同士で加減算ができます。

csr1 = sparse.csr_matrix([[1, 2], [0, 4]])
csr2 = sparse.csr_matrix([[5, 0], [0, 8]])
csc1 = sparse.csc_matrix([[1, 2], [0, 4]])
ndarray1 = np.array([[1, 0], [3, 4]])

# 加算
result_add = csr1 + csr2
print(type(result_add), result_add)
# <class 'scipy.sparse.csr.csr_matrix'>
#  (0, 0)  6
#  (0, 1)  2
#  (1, 1)  12

result_add_ndarray = csr1 + ndarray1 # ndarrayとも加算可
print(type(result_add_ndarray), result_add_ndarray)
# <class 'numpy.matrixlib.defmatrix.matrix'>
# [[2 2]
#  [3 8]]

result_add_csr_csc = csr1 + csc1 # matrix_csrとmatrix_cscで加算可
print(type(result_add_csr_csc), result_add_csr_csc)
# <class 'scipy.sparse.csr.csr_matrix'>   
#  (0, 0)  2
#  (0, 1)  4
#  (1, 1)  8

# 減算
result_sub = csr1 - csr2
print(type(result_sub), result_sub)
# <class 'scipy.sparse.csr.csr_matrix'>
#  (0, 0)  -4
#  (0, 1)  2
#  (1, 1)  -4

result_sub_ndarray = csr1 - ndarray1 # ndarrayとも減算可
print(type(result_sub_ndarray), result_sub_ndarray)
# <class 'numpy.matrixlib.defmatrix.matrix'>
# [[ 0  2]
#  [-3  0]]

内積

numpy.ndarrayと異なり、scipy.sparseでは演算子*が内積の計算となります(numpyでは要素ごとの積)。また、dot()でも内積を計算できます。
numpyのクラスがオペランドの場合、要素ごとの積になることに注意してください。

# 内積
result_inner_product = csr1 * csr2
print(result_inner_product)
#  (0, 1)  16
#  (0, 0)  5
#  (1, 1)  32

# dot()でも可能
result_dot = csr1.dot(csr2)

# ndarrayがオペランドの場合、要素ごとの積
result_ndarray_csr = np.array([[1, 0], [3, 4]]) * csr1
print(type(result_ndarray_csr), result_ndarray_csr)
# <class 'numpy.ndarray'> 
# [[1 2]
#  [3 22]]

# csr_matrixがオペランドの場合、内積
result_csr_ndarray = csr1 * np.array([[1, 0], [3, 4]])
print(type(result_csr_ndarray), result_csr_ndarray)
# <class 'numpy.ndarray'> 
# [[7 8]
#  [12 16]]

アダマール積(要素ごとの積)

要素ごとの積が必要な場合は、multiply()で計算できます。

result_multiply = csr1.multiply(csr2)
print(result_multiply)
#  (0, 0)  5
#  (1, 1)  32

固有値・固有ベクトル

固有値・固有ベクトルを求める場合、scipy.sparse.linalgをインポートして、eigsで計算します。
引数kは求める固有値・固有ベクトルの数(デフォルトは6)ですが、行列の次元数-1未満である必要があります。全ての固有値・固有ベクトルを求めることができません。

import scipy.sparse.linalg


csr3 = sparse.csr_matrix([[1, 0, 3], [0, 5, 0], [7, 0, 9]], dtype=np.float64)

# 固有値・固有ベクトル
eigen_values, eigen_vectors = sparse.linalg.eigs(csr3, k=1)
print(eigen_values)
# [ 11.08276253+0.j]
print(eigen_vectors)
# [[ -2.85181797e-01+0.j]
#  [  7.54125462e-17+0.j]
#  [ -9.58473444e-01+0.j]]

特異値分解

特異値分解は、scipy.sparse.linalg.svdsで行います。
こちらも引数kで求める特異値の数(デフォルト6)で指定しますが、行列の低い方の次元数(2×3行列の場合は2)未満にする必要があります。

csr4 = sparse.csr_matrix([[1, 0, 3], [0, 5, 0]], dtype=np.float64)

# 特異値分解
u, s, v_t = sparse.linalg.svds(csr4, k=1)
print(u)
# [[  4.04373124e-17]
#  [  1.00000000e+00]]
print(s)
# [ 5.]
print(v_t)
# [[  8.08746247e-18   1.00000000e+00   2.42623874e-17]]

性能

最後にnumpyとscipy.sparseの各クラスで性能を比較してみたいと思います。

numpyとscipy.sparseの性能比較

疎行列を使った計算で一番多いのは内積かと思います。
そこで、1,000×1,000(1,000,000要素)の疎行列(非ゼロ要素が10%)を作り、それぞれで内積を計算してその時間を計測・比較します。
また、疎行列でない場合にndarrayよりも性能劣化するのかどうかを検証するために、同じサイズの密行列(非ゼロ要素が100%)も作って内積を行います。

以下の結果となりました。カッコ内は試行回数で、10回以上の場合は上位3件の平均値です。

• 疎行列の場合、ndarrayと比較すると、lil_matrixでもcsr_matrixでも高速
  • 特にcsr_matrixは、ndarrayの約3,000倍、lil_matrixの約5倍、高速
• 密行列の場合でも、ndarrayと比較して性能が低下する、ということはなかった
  • 疎行列になる"かも"しれないケースにおいても、scipy.sparseのクラスに変換・利用する、というのも1つの手段ということが示唆されます

クラス	疎行列	密行列
ndarray	2.70s (1)	2.66s (1)
lil_matrix	4.57ms (100)	2.49s (1)
csr_matrix	863µs (1000)	2.30s (1)

# 1000×1000の疎行列の作成
dimensions = 1000
sparse_matrix = np.zeros((dimensions, dimensions), dtype=np.int64)

non_zero_elements = 10000 # 非ゼロ要素は全体の1%
count = 0
rows = np.random.randint(0, dimensions, 20000)
columns = np.random.randint(0, dimensions, 20000)
values = np.random.randint(0, 10, 20000)
for i in range(0, 20000):
    if sparse_matrix[rows[i], columns[i]] == 0:
        sparse_matrix[rows[i], columns[i]] = values[i]
        count += 1
        if count == non_zero_elements:
            break

sparse_lil = sparse.lil_matrix(sparse_matrix)
sparse_csr = sparse.csr_matrix(sparse_matrix)

# 疎行列の内積速度の計測
%timeit np.dot(sparse_matrix, sparse_matrix)
%timeit sparse_lil * sparse_lil
%timeit sparse_csr * sparse_csr

# 1000×1000の密行列の作成
dense_matrix = np.random.randint(1, 10, (dimensions, dimensions))
dense_lil = sparse.lil_matrix(dense_matrix)
dense_csr = sparse.csr_matrix(dense_matrix)

# 密行列の実行速度の計測
%timeit np.dot(dense_matrix, dense_matrix)
%timeit dense_lil * dense_lil
%timeit dense_csr * dense_csr

csr_matrixとcsc_matrixの性能比較

最後に、csr_matrixとcsc_matrixについて、行単位と列単位でアクセスした場合の速度について計測・比較してみます。

scipy.sparseの疎行列クラスで行・列にアクセスするには、スライスを使う方法とメソッドを使う方法の2つがあります。それぞれの方法で計測します。

• スライスを使う方法
  • i行目にアクセスするなら、matrix[i]
  • i列目にアクセスするなら、matrix[:, i]
• メソッドを使う方法
  • i行目にアクセスするなら、matrix.getrow(i)
  • i列目にアクセスするなら、matrix.getcol(i)

結果は以下のとおりです。

• 総じて、スライスよりも、メソッドを使う方が高速
  • 行単位・列単位でアクセスしたい場合は、極力メソッドを使う方が良さそうです
• 行アクセスでは、csr_matrix > lil_matrix > csc_matrixの順で、予想通りでした
• 列アクセス(メソッドの場合)では、csc_matrix > csc_matrix >> lil_matrixの順でした
  • lil_matrixは遥かに低速で、何らか計算を行う際にはlil_matrixから別のクラスに変換した方が良さそうです
  • スライスの場合は、csc_matrixよりもcsr_matrixの方が高速になっており、予想に反した結果となりました

クラス	行アクセス(スライス)	列アクセス(スライス)	行アクセス(getrow())	列アクセス(getcol())
lil_matrix	498ms (1)	7.40s (1)	469ms (1)	2.41s (1)
csr_matrix	330ms (1)	386ms (1)	304ms (1)	367ms (1)
csc_matrix	504ms (1)	442ms (1)	425ms (1)	305ms (1)

sparse_csc = sparse.csc_matrix(sparse_matrix)

# 行単位(スライス)で平均値を算出
%timeit [sparse_lil[i].mean() for i in range(dimensions)]
%timeit [sparse_csr[i].mean() for i in range(dimensions)]
%timeit [sparse_csc[i].mean() for i in range(dimensions)]

# 列単位(スライス)で平均値を算出
%timeit [sparse_lil[:, i].mean() for i in range(dimensions)]
%timeit [sparse_csr[:, i].mean() for i in range(dimensions)]
%timeit [sparse_csc[:, i].mean() for i in range(dimensions)]

# 行単位(getrow())で平均値を算出
%timeit [sparse_lil.getrow(i).mean() for i in range(dimensions)]
%timeit [sparse_csr.getrow(i).mean() for i in range(dimensions)]
%timeit [sparse_csc.getrow(i).mean() for i in range(dimensions)]

# 列単位(getcol())で平均値を算出
%timeit [sparse_lil.getcol(i).mean() for i in range(dimensions)]
%timeit [sparse_csr.getcol(i).mean() for i in range(dimensions)]
%timeit [sparse_csc.getcol(i).mean() for i in range(dimensions)]

まとめ

今回はscipy.sparseを使った疎行列計算について、調べて紹介しました。
特にlil_matrix、csr_matrix、csc_matrixに焦点を当て、計算方法や性能比較を行っています。

非ゼロ要素の割合が数%以下の疎行列はしばしば現れますので、計算ツールとして覚えておくと便利かと思います。

参考

参考にさせていただいた投稿です。

scipy.sparseの概要の把握に役立ちます。
計算高速化のための強力な武器、NumPyのブロードキャスティングとSciPyの疎行列処理：Pythonで始める機械学習入門（4） - ＠IT

csr_matrixとcsc_matrixの内部構造の違いについて解説されています。
scipy.sparseの内部データ構造 – はむかず！

Python: レコメンドの行列分解を確率的勾配降下法で実装する

過去3回の投稿で、行列分解(SVDとNMF)によるレコメンドを実装してきました。

• Pythonで特異値分解(SVD)を理解する - け日記
• SVDでMovieLensのレコメンドを実装する - け日記
• NMFでMovieLensのレコメンドを実装する - け日記

ですが、いずれも欠測値やユーザごと・アイテムごとのバイアスに対応していませんでした。

今回の投稿では、欠測値やバイアスの概念を導入した行列分解について理解し、確率的勾配降下法で実装してみます。

欠測値とバイアスを考慮した行列分解(SVD)

欠測値とバイアスを考慮した行列分解について、Simon Funkのブログで投稿されています。元々はNetflix Prize(wikipedia)で考えられたアルゴリズムだそうです。

Netflix Update: Try This at Home

解説としてはこちらの論文が参考になります。

• Y.Koren、R.Bell、C.Volinsky ： MATRIX FACTORIZATION TECHNIQUES FOR RECOMMENDER SYSTEM (2009) (PDF)
• Y.Koren、R.Bell ： Advances in Collaborative Filtering (2011) (PDF)

Simonのアルゴリズムでは、ユーザuのアイテムiに対する評価値は、以下式で予測します。

• μは全評価値の平均値
• はアイテムiのバイアスを表しており、アイテムごとの高評価・低評価の差が反映されています
• はユーザuのバイアスを表しており、ユーザごとの評価の甘さ・辛さの差が反映されています
• qとuが分解された評価値行列となります
  • qは、アイテムごとの特徴がk次元で表現された行列で、k行×M列(Mはアイテム数)
    • したがって、は、アイテムiのk次元の特徴
  • pは、ユーザごとの特徴がk次元で表現された行列で、k行×N列(Nはユーザ数)
    • したがって、は、ユーザuのk次元の特徴

なお、推薦システムの界隈ではこの方法が一般的にSVDと呼ばれているようです。固有値から導出する数学的な特異値分解(SVD)と異なりますが、同じ"SVD"という用語が使われていますので、文脈で読み分ける必要があります。この投稿でも、SVDと呼称します。

確率的勾配降下法によるSVD

上の式では、バイアス(、)と分解後の行列(q, p)が未知のため、何らかの方法で導出する必要があります。

Simonの投稿では、確率的勾配降下法(SGD)が使われています。
SGDは、データ(ここでは評価値)をランダムに1つずつ取り出し、目的関数の値が小さくなるように変数を更新していく方法です。SGDについてはこちらの投稿が参考になります。

postd.cc

SGDで最小化を目指す目的関数は以下式となります。

• 第1項は、予測値と実際値の誤差()の2乗和です
• 第2項は、正則化項です
  • ハイパーパラメータλは正則化制約の強さを表し、値が大きいほど強くなります
• 値がある評価値()のみで計算され、欠測値は除外されていることに注目してください

各値は以下の式で更新します。

• 
• 
• 
• 
• は誤差で、
• ハイパーパラメータγは学習率を表し、値が大きいほど1回の学習で値が大きく更新されます

実装

それではPythonで実装していきます。

例として、5本の映画に対して、4人のユーザが5段階(1〜5)で評価した、というケースを想定します(前回の投稿と同じです)。
未評価の値(欠測値)は0と置いておきます。この後の計算では除外されます。

import numpy as np
import pandas as pd
np.set_printoptions(precision=2, suppress=True)


# 評価値行列
x = np.matrix([[5, 5, 0, 3], [4, 4, 0, 0], [0, 3, 4, 0], [0, 1, 4, 5], [1, 0, 5, 4]])
# matrix([[5, 5, 0, 3],
#         [4, 4, 0, 0],
#         [0, 3, 4, 0],
#         [0, 1, 4, 5],
#         [1, 0, 5, 4]])

変数の初期化

q、p、、を、-0.05〜+0.05の小さな乱数で初期化します(初期値は0でも良いかもしれません)。

def get_initial_values(x, k):
    # 各行列またはベクトルの値を、-0.05〜+0.05の乱数で初期化
    mu = np.sum(x) / np.count_nonzero(x)
    bu = (np.random.rand(x.shape[1]) - 0.5) * 0.1
    bi = (np.random.rand(x.shape[0]) - 0.5) * 0.1
    q = (np.random.rand(k, x.shape[0]) - 0.5) * 0.1
    p = (np.random.rand(k, x.shape[1]) - 0.5) * 0.1
    
    return (mu, bu, bi, q, p)


k = 2
get_initial_values(x, k)
# mu: 3.69
# bu:  array([-0.03,  0.  ,  0.  ,  0.04])
# bi:  array([-0.01, -0.02, -0.01, -0.05, -0.  ])
# q:
#  array([[-0.04, -0.01, -0.01,  0.02, -0.03],
#        [ 0.03, -0.04, -0.04, -0.02,  0.04]]),
# p:
#  array([[-0.04, -0.03,  0.  , -0.03],
#         [ 0.02,  0.02,  0.04,  0.02]]))

誤差の計算

誤差と目的関数の計算を実装します。

• get_error_matrix関数は、1つの誤差の値が1要素となる行列にして返します
• error_function関数は、目的関数の実装で、欠測値(0)はここで除外することで、学習へ影響を与えないようにします

def get_error_matrix(x, mu, q, p, bu, bi):
    # 誤差eの計算
    # 計算の簡略化のため、評価値0は除外していない
    return x - (mu + np.dot(q.T, p) + bu + np.matrix(bi).T) 


def error_function(x, mu, q, p, bu, bi, l):
    # 誤差に正則化制約を加えた目的関数を定義
    error_matrix = get_error_matrix(x, mu, q, p, bu, bi)
    error = np.sum(np.square(error_matrix[x > 0])) # 0を除外して誤差を計算
    
    regularization = l * (np.sum(np.square(bu)) + np.sum(np.square(bi)) + np.sum(np.square(q)) + np.sum(np.square(p)))
    
    return error + regularization


error_function(x, mu, q, p, bu, bi, 0.02)
# 23.03739901928402

SGDによる値の更新

最後にSGDで値を更新する実装です。

各エポックの最初に評価値をランダムに並び替えています(ここでも0は除外しています)。
そして、評価値を先頭から1つずつ取り出し、目的関数で値を求めて、各値を更新しています。

次元数kは2、エポック数は100の条件下で、目的関数の値は22.7から0.87まで下がります。

def sgd(x, k, epochs=100, l=0.02, gamma=0.02):
    # 確率的勾配法で分解した行列とバイアスを求める
    mu, bu, bi, q, p = get_initial_values(x, k)
    errors = [] # エポック毎の目的関数の値
    
    for epoch in range(epochs):
        error = error_function(x, mu, q, p, bu, bi, l)
        errors.append(error)
        
        # 欠測値を除いた評価値のインデックスを取得してランダムに並び替える
        x_i, x_u = np.where(x > 0)
        targets = np.arange(len(x_i))
        np.random.shuffle(targets)
        
        for target in targets:
            # 誤差行列を計算して、そのうち1つを選択
            error_matrix = get_error_matrix(x, mu, q, p, bu, bi)
            
            i = x_i[target]
            u = x_u[target]
            
            e_ui = error_matrix[i, u]
            
            # 分解した行列とバイアスを更新
            bu[u] = bu[u] + gamma * (e_ui - l * bu[u])
            bi[i] = bi[i] + gamma * (e_ui - l * bi[i])
            q[:, i] = q[:, i] + gamma * (e_ui * p[:, u] - l * q[:, i])
            p[:, u] = p[:, u] + gamma * (e_ui * q[:, i] - l * p[:, u])

    error = error_function(x, mu, q, p, bu, bi, l)
    errors.append(error)
    print('Error: {}'.format(error))
    # Error: 0.8715459314064704

    expected = mu + bu + np.matrix(bi).T + np.dot(q.T, p)
    
    return expected, errors
    

expected, errors = sgd(x, 2)

得られた値は以下になります。アイテムやユーザの特徴・嗜好というよりは、アイテムやユーザごとのバイアスからの影響が強そうです。

• では、ユーザA(1列目)のユーザが辛め、ユーザC(3列目)のユーザが甘めの傾向が反映されています
• では、タイタニック(2行目)のアイテムが高評価、サマーウォーズ(4行目)のアイテムが低評価となっています
• qは、1列目を見ると、ローマの休日(1行目)とサマーウォーズ(4行目)、AKIRA(5行目)が近くなっているようです
• pは、1行目を見るとユーザA(1列目)とB(2列目)、また、2行目ではBとC(3列目)が近いことを示しているようです

0で置き換えていた欠測値にも、予測された評価値が入っており、例えばユーザAとBのアニメーション映画に対する評価は低い、といった嗜好が反映されています。

$$ \hat{R}=\mu + b_u + b_i + q^Tp = 3.7+\left( \begin{matrix} -0.8 & -0.6 & 1.0 & 0.2 \\ -0.8 & -0.6 & 1.0 & 0.2 \\ -0.8 & -0.6 & 1.0 & 0.2 \\ -0.8 & -0.6 & 1.0 & 0.2 \\ -0.8 & -0.6 & 1.0 & 0.2 \end{matrix} \right)+\left( \begin{matrix} 0.7 & 0.7 & 0.7 & 0.7 \\ 0.8 & 0.8 & 0.8 & 0.8 \\ -0.4 & -0.4 & -0.4 & -0.4 \\ -0.7 & -0.7 & -0.7 & -0.7 \\ -0.6 & -0.6 & -0.6 & -0.6 \end{matrix} \right)+\left( \begin{matrix} 1.2 & 0.4 \\ -0.3 & 0.0 \\ -0.2 & 0.0 \\ 1.1 & 0.7 \\ 0.9 & 0.2 \end{matrix} \right) \left( \begin{matrix} -1.2 & -0.9 & 0.4 & 1.1 \\ 0.0 & 0.5 & 0.2 & -0.6 \end{matrix} \right)= \left( \begin{matrix} 4.9 & 5.1 & 5.4 & 3.1 \\ 3.9 & 4.0 & 5.6 & 4.6 \\ 2.7 & 2.8 & 4.3 & 3.3 \\ 1.2 & 1.1 & 4.1 & 4.9 \\ 1.4 & 2.0 & 4.5 & 4.1 \end{matrix} \right) $$

エポック毎に目的関数の値をプロットしました。エポック毎に0に近づくことがわかります。

まとめ

レコメンドの欠測値とバイアスを考慮した行列分解(SVD)を確率的勾配降下法(SGD)で実装し、小さなデータで概ね納得できる結果が得られました。
次回以降の投稿で、MovieLens100kに適用してみたいと思います。

NMFでMovieLensのレコメンドを実装する

前回の投稿( SVDでMovieLensのレコメンドを実装する - け日記 )に引き続き、今回はNMFで映画のレコメンドを実装します。

NMF

非負値行列因子分解(Non-negative Matrix Factorization: NMF)は、ある行列XをW・Hに近似的に分解する方法の一つです。その名の通り、各行列の要素は全て0以上となります。

SVDなどの行列分解では、負の値が現れることがあります。
しかし、実世界で使われる値(画素、頻度、評価値など)の多くは非負値で、分解した結果に解釈を与えるなど非負の性質を保ちたいケースもあります(負の画素値などは存在しませんし)。 そうした場合にNMFが使われるようです。

NMFでは、M×Nの行列Xを、N×kの行列Wとk×Mの行列Hに分解します。kはM,Nよりも小さな値が使われます。

未知のWとHを求めるために、Xとの距離を定義します。
ここでは一番直感的なフロベニウスノルムを使います。

学習では、この距離(つまり誤差)が小さくなるように、HとWを更新します。

上述以外にもいくつか距離の定義と更新式が提案されていますが、割愛します。こちらの方の投稿がよくまとまっています。

非負値行列因子分解（NMF : Non-negative Matrix Factorization）を改めてやり直してみた - Qiita

NMFを推薦システムに適用する

以下のように、5本の映画に対して、4人のユーザが5段階(1〜5)で評価した、というケースを想定します。
恣意的に選択しているのでそうなっているのですが、これら5本の映画にはロマンス要素とアニメーション要素の2つがあり、ユーザもこの2つのグループに分けられそうです(AとBはロマンス要素が強い映画、CとDはアニメーション要素が強い映画をそれぞれ好んでいる)。

この表は5×4の行列(ここではX)で表現できます。なお、未評価となっている要素は0で埋めています。

XをNMFでWとHに分解します。ここではk=2としています。

分解されたWとHを見ると、映画の持つ要素やユーザの嗜好が反映された行列になっていることがわかります。

• Wは映画数×kの行列
  • 上2行(ローマの休日、タイタニック)は1列目、下2行(サマーウォーズ、AKIRA)は2列目の値がそれぞれ大きくなっており、映画の持つ要素が各列に反映されているように見えます
• Hはk×ユーザ数の行列
  • 左2列は1行目、右2列は2行目の値がそれぞれ大きくなっており、ユーザごとの嗜好がグループ分けされているように見えます
• kによって集約された特徴と捉えることができます

レコメンドの実装

それではNMFを使ったレコメンドを実装します。データセットの準備とランキングの作成については前回の投稿( SVDでMovieLensのレコメンドを実装する - け日記 )とほぼ同じです。

データセットの準備

今回もMovieLens 100K Datasetを使います。

以下からダウンロードして、今回作成するPythonファイルと同じディレクトリに解凍しておきます。以前の投稿( word2vecでitem2vecを実装して映画を推薦する - け日記 )も参考にしてみてください。

MovieLens | GroupLens

評価値をnumpy行列にして持ちます。欠測値は0で補完しています。あわせてアイテムIDとタイトルの一覧もDataFrameにロードしておきます。

import numpy as np
import pandas as pd
import codecs


# 評価値を持つファイルはu.data。楽に行列化するために、まずはDataFrameでロードする
ratings = pd.read_csv('ml-100k/u.data', delimiter='\t', header=None).iloc[:, :3]
ratings.rename(columns={0: 'user_id', 1: 'item_id', 2: 'rating'}, inplace=True)
print('ratings.shape: {}'.format(ratings.shape)) # ratings.shape: (100000, 3)

# item_idを行、user_idを列とする1682×943の行列(欠測値は0埋め)
rating_matrix = ratings.pivot(index='item_id', columns='user_id', values='rating').fillna(0).as_matrix()
print('rating_matrix.shape: {}'.format(rating_matrix.shape)) # rating_matrix.shape: (1682, 943)

# pd.read_csvで読み込もうとすると文字コードエラーになるため、
# codecs.openでエラーを無視してファイルを読み込んだ後にDataFrameにする
with codecs.open('ml-100k/u.item', 'r', 'utf-8', errors='ignore') as f:
    item_df = pd.read_table(f, delimiter='|', header=None).ix[:, 0:1]
    item_df.rename(columns={0: 'item_id', 1: 'item_title'}, inplace=True)
    item_df.set_index('item_id', inplace=True)

items = pd.Series(item_df['item_title'])
print('items.shape: {}'.format(items.shape)) # items.shape: (1682, 1)

評価ランキングの予測と表示

指定したユーザに対して、未評価のアイテムの評価値を予測し、ランキングを返す関数を作ります。

未評価の映画については0で補完しており、また映画やユーザごとの評価値の平均なども考慮していないため、絶対値の予測精度は高くないと予想されます。そのため、ランキングにしています。

def predict_ranking(user_index, reducted_matrix, original_matrix, n):
    # 対象ユーザのSVD評価値
    reducted_vector = reducted_matrix[:, user_index]
    
    # 評価済みのアイテムの値は0にする
    filter_vector = original_matrix[:, user_index] == 0
    predicted_vector = reducted_vector * filter_vector

    # 上位n個のアイテムのインデックスと予測した評価値を返す
    return [(i, predicted_vector[i]) for i in np.argsort(predicted_vector)[::-1][:n]]


def print_ranking(user_ids, items, reducted_matrix):
    for user_id in user_ids:
        predicted_ranking = predict_ranking(user_id - 1, reducted_matrix, rating_matrix, 10)
        print('User: {}:'.format(user_id))
        for item_id, rating in predicted_ranking:
            # アイテムID, 映画タイトル, 予測した評価値を表示
            print('{}: {} [{}]'.format(item_id, items[item_id + 1], rating))

NMFによる行列分解

NMFによる行列分解ですが、sklearn.decomposition.NMFでお手軽にできますので、これを使います。

• n_componentsにk=30を指定しています(前回のSVDと同じ)
• initで初期値の設定方法を指定できます(今回はランダム)
• MAE(絶対平均誤差)は約0.26
  • XとWHの各要素は平均0.26の差がある、ということです

from sklearn.decomposition import NMF


def reduct_matrix_nmf(matrix, k):
    # NMFで分解
    model = NMF(n_components=k, init='random', random_state=1)
    w = model.fit_transform(matrix)
    h = model.components_
    
    return np.dot(w, h)


reducted_matrix_nmf = reduct_matrix_nmf(rating_matrix, k=30)

print('MAE: {}'.format(np.average(np.absolute(rating_matrix - reducted_matrix_nmf))))
# MAE: 0.25617962729902044

評価

それでは実際に推薦をさせてみます。

前回のSVDの結果と比較するために、前回と同じ2人のユーザ(ユーザID 267と868)を対象にランキングを生成しました。
概ねSVDの時と同じようなタイトルが得られましたが、負値を含まない分、予測評価値は全体的にNMFの方が高いことがわかります。

• この2人は"Akira(1988)"と"Ghost in the Shell(Kokaku kidotai)(1995)"を5で評価しています(つまりSF+アニメーション好きです)

User: 267:
95: Terminator 2: Judgment Day (1991) [4.355731842820824]
10: Seven (Se7en) (1995) [3.5465670297733514]
231: Young Guns (1988) [3.3581447476132875]
78: Fugitive, The (1993) [3.3281046778126173]
153: Monty Python's Life of Brian (1979) [3.3028868626155874]
41: Clerks (1994) [3.1889273065852977]
172: Princess Bride, The (1987) [3.158916000181793]
356: One Flew Over the Cuckoo's Nest (1975) [3.1348566836324423]
0: Toy Story (1995) [3.128039827666597]
3: Get Shorty (1995) [2.8694586381163]

User: 868:
174: Brazil (1985) [4.270610564574103]
257: Contact (1997) [2.965456387373195]
356: One Flew Over the Cuckoo's Nest (1975) [2.914023120253882]
143: Die Hard (1988) [2.5656501189688483]
430: Highlander (1986) [2.5194608160565846]
195: Dead Poets Society (1989) [2.319622366933584]
21: Braveheart (1995) [2.2844329680205355]
41: Clerks (1994) [2.212755642335326]
379: Star Trek: Generations (1994) [2.212123031118453]
27: Apollo 13 (1995) [2.1517770562746565]

まとめ

今回はNMFを使った映画のレコメンドを実装しました。

SVDが含んでいた負値を取り除くことができましたが、欠測値や、ユーザごとのバイアス(甘め・辛め)が考慮されたモデルには、まだなっていません。

次は、欠測値やバイアスを考慮した推薦システムに取り組みたいと思います。

参考にした文献・記事

岩波データサイエンス Vol.5

qiita.com

http://lab.adn-mobasia.net/?p=1220lab.adn-mobasia.net

非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみる | かものはしの分析ブログ