推荐系统算法：矩阵分解_矩阵分解 python 基于rmse

---

基于模型的协同过滤推荐算法：隐语义模型之矩阵分解

在推荐系统算法中，协同过滤算法更适合给用户进行个性化的推荐，其中基于模型的协同过滤算法相较于基于领域的协同过滤算法（基于内存的算法），在大数据处理上，有更大的优势。

隐语义模型：

用户与物品之间存在着隐含的联系，基于聚类算法，将商品和顾客分别聚类，通过各个类别之间的关系，得到顾客与商品之间存在的隐含联系。隐含的特征就是类别的分类个数称为隐分类个数（隐含特征 Latent Factor），隐分类个数可多可少，相比较之下模型的可解释都不强，因为隐含特征是自动进行分类的。

一、稀疏矩阵分解

推荐系统中存在着用户行为的稀疏矩阵，由于稀疏矩阵无法给用户进行商品推荐，因此考虑补全这些缺失值，并根据缺失值TOP-N 给用户进行个性化推荐。依据打分商品推荐未打分商品。

可以对稀疏矩阵进行分解，假设有如下用户对电影的观影喜好：

假设隐分类个数：3

我们可以将评分矩阵分解为：User 矩阵和 Item 矩阵，User 矩阵代表用户对电影隐分类类别的喜好程度；item 矩阵代表各个电影都属于哪一类隐分类类别。如果可以分解，那么相乘必将得到完整的评分矩阵。

User 矩阵为 $X=[x_1,x_2,...,x_N]$隐分类个数 k=3 维列向量，
Item 矩阵为 $Y=[y_1,y_2,...y_M]$ 隐分类个数 k=3 维列向量，
$r_{ui}$ 在原始评分矩阵中若存在值，则为 1 ，else 0.

矩阵分解的目标函数：
$$Loss=\sum _{r_{ui}}(r_{ui}-x_u^T{y}_i{)}^2+\lambda (\sum _N||x_u|{|}_2^2+\sum _M||y_i|{|}_2^2)$$
L1 正则化与 L2 正则化
本身模型就是稀疏的，不需要使用 L1 正则化。接下来最小化 Loss 就可以了。

优化这个损失函数有两种方式：梯度下降法和 ALS 交替最小二乘法。对二维（用户、商品）矩阵的分解可以使用 “ALS交替最小二乘” 求解

二、 矩阵分解的优化方式

（一）、 ALS(Alternating Least Squares，交替最小二乘)

思想：固定一个求另一个最优化

1. 优化显式信息矩阵：

对 x 进行优化：
$$Loss=\sum _{r_{ui}}(r_{ui}-x_u^T{y}_i{)}^2+(\sum _N||x_u|{|}_2^2)$$

$$\frac{\partial L(x_u)}{x_u}=-2Y_u(R_u-Y_u{x}_u)+2\lambda x_u=0$$

$$x_u=(Y_u{Y}_u^T+\lambda I{)}^{-1}{Y}_u{R}_u$$

然后对 y 进行优化。

2. 对隐式信息矩阵进行优化：

显示矩阵转换：
$p_{ui}=1,r_{ui}>=1$
$p_{ui}=0,r_{ui}=0$

设置置信度：
$c_{ui}=1+\alpha r_{ui}$
$r_{ui}>0$，$c_{ui}$会线性的增长，$r_{ui}=0$，$c_{ui}=1$

目标函数：
$$Loss=\sum _{p_{ui}}{c}_{ui}(p_{ui}-x_u^T{y}_i{)}^2+\lambda (\sum _N||x_u|{|}_2^2+\sum _M||y_i|{|}_2^2)$$
使其达到最小即为目的。

1. 工具

• spark mllib库（spark3.0版本后废弃）
  支持常见的算法，包括 分类、回归、聚类和协同过滤
  from pyspark.mllib.recommendation import ALS, Rating, MatrixFactorizationModel
• spark ml库（官方推荐）
  功能更全面更灵活，ml在DataFrame上的抽象级别更高，数据和操作耦合度更低，使用起来像sklearn
  Spark安装很多坑，需要慢慢来
• Python代码
  https://github.com/tushushu/imylu/blob/master/imylu/recommend/als.py
• surprise

2. 示例：使用 MovieLens 数据集进行评分预测

方法 1 ： 使用 ALS 算法，采用分布式 spark 计算用户（userid == 100）的电影评分

• 模型目标函数：

$$Loss=\sum _{r_{ui}}(r_{ui}-x_u^T{y}_i{)}^2+\lambda (\sum _N||x_u|{|}_2^2+\sum _M||y_i|{|}_2^2)$$

• Loss 最小化的求解方法： ALS（通过先固定一个求解另一个的优化方式 ）

1. 初始化Items矩阵；
2. 固定Items矩阵，求Users矩阵使得Loss函数最小化；

当 $y_i$ 固定时，将其写为矩阵形式，可以对单个 $x_u$ 进行优化求解：
$$\frac{\partial L(x_u)}{x_u}=-2Y_u(R_u-Y_u{x}_u)+2\lambda x_u=0$$
$$x_u=(Y_u{Y}_u^T+\lambda I{)}^{-1}{Y}_u{R}_u$$
所以，可以将 $x_u$ 进行并行化操作，在每个分区执行不同的 $x_u$ ,在操作完成时，再合并到一起。

3. 固定Users矩阵，求Items矩阵使得Loss函数最小化；
4. 重复上述2，3步过程

• ALS 可以并行处理的原因：

当 $y_i$ 固定时，将其写为矩阵形式，可以对单个 $x_u$ 进行优化求解：
$$\frac{\partial L(x_u)}{x_u}=-2Y_u(R_u-Y_u{x}_u)+2\lambda x_u=0$$
$$x_u=(Y_u{Y}_u^T+\lambda I{)}^{-1}{Y}_u{R}_u$$
所以，可以将 $x_u$ 进行并行化操作，在每个分区执行不同的 $x_u$ ,在操作完成时，再合并到一起。

• spark中 ALS 参数设置：

• 显式行为：

numBlocks 是用于并行化计算的分块个数 (设置为-1为自动配置)。
maxTter: 最小迭代次数；
regParam：惩罚项系数$\lambda$（建议为0.01）
rank：隐分类个数
userCol= , 用户变量名称
itemCol= , 商品变量名称
ratingCol= ,评分变量名称
nonnegative:是否使用非负约束，默认不使用 false。
predictionCol:预测列的名字

• 隐式行为：

ImplicitPrefs = true此处表明rating矩阵是隐式评分

from pyspark.ml.recommendation import ALS
from pyspark import SparkContext
from pyspark.sql import SQLContext
from pyspark.sql.types import *
import pandas as pd
import os

#Spark 环境配置
sc = SparkContext()
sql_sc = SQLContext(sc)

#数据加载
pd_df_ratings = pd.read_csv('./ratings.csv')
#转换为 spark 下的 DataFrame 格式
pyspark_df_ratings = sql_sc.createDataFrame(pd_df_ratings)
pyspark_df_ratings = pyspark_df_ratings.drop('Timestamp')

print(pyspark_df_ratings.show(5, truncate=False))

# 创建ALS模型
als = ALS(rank = 8, 
		  maxIter = 15, 
		  regParam = 0.1, 
		  userCol= 'userId', 
		  itemCol='movieId', 
		  ratingCol='rating'
		  )
model = als.fit(pyspark_df_ratings)

# 对userId=100进行Top-N推荐
recommendations = model.recommendForAllUsers(5)

# 对user-id == 100 进行电影推荐
print(recommendations.where(recommendations.userId == 100).collect())
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

方法2：基于 surprise 的 Baseline 模型的 ALS 优化方法

• BaseLine的基本思想：

Baseline算法的思想就是设立基线 $\mu$，并引入用户的偏差 $b_u$ 以及item的偏差 $b_i$
$$\widehat{r_{ui}}=b_{ui}=\mu +b_u+b_i$$

• BaseLine的目标函数：

$$Loss=\sum _{{\hat{r}}_{ui}}({\hat{r}}_{ui}-\mu -b_u-b_i{)}^2+\lambda (\sum _N||b_u|{|}_2^2+\sum _M||b_i|{|}_2^2)$$

• 优化方法：ALS 交替最小二乘

• 验证方法：K 折交叉验证

#from surprise import SVD
from surprise import Dataset
from surprise import Reader
from surprise import BaselineOnly, KNNBasic, NormalPredictor
from surprise import accuracy
from surprise.model_selection import KFold
#import pandas as pd

# 数据读取
reader = Reader(line_format='user item rating timestamp', sep=',', skip_lines=1)
data = Dataset.load_from_file('./ratings.csv', reader=reader)
train_set = data.build_full_trainset()

# ALS优化
bsl_options = {
			   'method': 'als',
			   'n_epochs': 5,
			   'reg_u': 12, # user 的正则化参数，默认为 15
			   'reg_i': 5 # item 的正则化参数， 默认为 10
			   }
# SGD优化
#bsl_options = {'method': 'sgd','n_epochs': 5}
algo = BaselineOnly(bsl_options = bsl_options)
#algo = BaselineOnly()
#algo = NormalPredictor()

# 定义K折交叉验证迭代器，K=3
kf = KFold(n_splits=3)
for trainset, testset in kf.split(data):
    # 训练并预测
    algo.fit(trainset)
    predictions = algo.test(testset)
    # 计算RMSE
    accuracy.rmse(predictions, verbose = True)

uid = str(196)
iid = str(302)

# 输出uid对iid的预测结果
pred = algo.predict(uid, iid, r_ui=4, verbose=True)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

3. ALS 的缺点：计算时所需要的时间与空间复杂度太高

一般来讲：

ALS 的方法只能优化较少的特征，对多特征时运行较慢；
但 SVD 的方法可以同时优化多种特征。

（二）、梯度下降法

三种梯度下降法：BGD、SGD、Mini-batch 梯度下降。

三、SVD 矩阵分解：

1. 原理：

一、对称方阵矩阵的分解：

对于对称方阵 $A$ 可以分解为：$A=VD{V}^T$ 的形式，
若矩阵 $A_{n,m}$ 不是对称方阵，可以将 $A_{n,m}$ 转换为 $A{A}_{n,n}^T$ 进行矩阵分解。

二、$A{A}^T$ 与 $A^{T}A$ 的特征值是相同的（$A_{m,n}$ ）:

若 $A^{T}A$ 存在特征值、特征向量 $\lambda ,v_i$ 则$$A^{T}A\cdot v_i=\lambda \cdot v_i$$ $$A{A}^T(A{v}_i)=A\lambda v_i$$ $$A{A}^T\cdot u_i=\lambda (u_i)$$所以两者有着相同的特征值$$A^T{A}_{n,n}=V_{n,n}D{V}^T,A{A}_{m,m}^T=U_{m,m}B{U}^T$$ $D,B$ 中的特征值是相同的 $\lambda$

三、非方阵的矩阵分解：

若 $A^{T}A$ 存在特征值、特征向量 $\lambda ,v_i$ 则$$A^{T}A\cdot v_i=\lambda \cdot v_i$$
由 二 可得：相应的$u_i=A{v}_i$也是特征值 $\lambda$ 的特征向量。
对 $u_i$ 进行标准化可以得到：$$w_i=\frac{A{v}_i}{|A{v}_i|}=\frac{A{v}_i}{\sqrt{\lambda }}$$ $$A{v}_i=\sqrt{\lambda }\cdot w_i$$
$$AV=W\Sigma$$ $$A=W\Sigma V^T$$

四、特性：

1. SVD中的奇异值 $\Sigma$ 为$D$ 的特征值开方；

2. 右奇异值向量为$A^{T}A$ 的特征向量；
3. 左奇异值向量为 $A{v}_i$的标准化向量。
4. $W_{m,m}$ $V_{n,n}$ 分别代表 $A$ 中的行特征矩阵与列特征矩阵，即：$User$ 矩阵与 $Item$ 矩阵

参考

2. 示例：

from scipy.linalg import svd
import numpy as np

from scipy.linalg import svd
A = np.array([[1,2],
	    [1,1],
	    [0,0]])
p,s,q = svd(A,full_matrices=False)

print('P=', p)
print('S=', s)
print('Q=', q)
1
2
3
4
5
6
7
8
9
10
11
12

在推荐系统中，矩阵通常是稀疏的，而 SVD 要求矩阵是稠密的。
事实上，SVD 更适合进行降维

SVD降维

GitHub：图像降维


可以看到，减少特征数量很少时，并不影响图像的表达，基于此，应用于推荐系统。

3. SVD 在推荐系统中的应用

可以看到使用两阶特征矩阵亦可以还原评分矩阵
这实际上也是一个隐语义模型，其隐特征维度为 2

Funk SVD、Bias SVD、SVD++

使用 SVD 构造的这几个隐语义模型其实本质上就是损失函数的不同：

• 由于 SVD 方法不能解决稀疏矩阵的问题，因此使用机器学习方法，在 SVD 降维的基础上(设置超参数 构造特征数量： $k$)。
• 只关注有值部分，构造损失函数，使用梯度下降的优化方式，得到最小Loss的两个矩阵，以达到矩阵分解的效果。
• 三者的区别在于损失函数的构造方式不同。本质与前述稀疏矩阵的 Loss 损失函数类似

1. FunkSvd： $$\underset{p_i,q_j}{\mathrm{argmin}}\sum _{i,j\in K}(m_{ij}-q_j^T{p}_i{)}^2+\lambda (||p_i|{|}_2^2+||q_j|{|}_2^2)$$
2. BiasSVD: 通过 baseline 的方式设置损失函数：
   $$\underset{p_i,q_j}{\mathrm{argmin}}\sum _{i,j\in K}(m_{ij}-\mu -b_u-b_i-q_j^T{p}_i{)}^2+\lambda (||p_i|{|}_2^2+||q_j|{|}_2^2+||b_u|{|}_2^2+||q{b}_i|{|}_2^2)$$
3. SVD ++ 加入了用户的隐式行为：
   $$\underset{p_i,q_j}{\mathrm{argmin}}\sum _{i,j\in K}(m_{ij}-\mu -b_u-b_i-q_j^T(p_i+|I_i{|}^{-\frac{1}{2}}+\sum _{j\in I_i}{y}_j){)}^2+\lambda (||p_i|{|}_2^2+||q_j|{|}_2^2+||b_u|{|}_2^2+||b_i|{|}_2^2+||\sum _{j\in I_i}{y}_j|{|}_2^2)$$

• 优化方式：梯度下降法

工具

准备工作:

上传Google云盘 movielens 数据集

使用Google colab

#安装必要的库和包
from google.colab import drive
drive.mount('/content/drive')

#进入安装环境
!ls "/content/drive/My Drive/"

#安装 surprise 包
!pip install surprise
1
2
3
4
5
6
7
8
9
10

数据读取：

from surprise import SVD
from surprise import Dataset
from surprise.model_selection import cross_validate
from surprise import Reader
from surprise import accuracy
from surprise.model_selection import KFold
import pandas as pd
import time


# from google.colab import drive
# drive.mount('/content/drive')
import os
os.chdir("/content/drive/My Drive/Colab Notebooks/")

time1=time.time()

# 数据读取
reader = Reader(line_format='user item rating timestamp', sep=',', skip_lines=1)
data = Dataset.load_from_file('ratings.csv', reader=reader)
train_set = data.build_full_trainset()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

(1) 使用Funk SVD方法：

time1=time.time()

# 使用funkSVD
algo = SVD(biased=False)# True 即为 Bias SVD

# 定义K折交叉验证迭代器，K=3
kf = KFold(n_splits=3)
for trainset, testset in kf.split(data):
    # 训练并预测
    algo.fit(trainset)
    predictions = algo.test(testset)
    # 计算RMSE
    accuracy.rmse(predictions, verbose=True)

uid = str(196)
iid = str(302)
# 输出uid对iid的预测结果
pred = algo.predict(uid, iid, r_ui=4, verbose=True)
time2=time.time()
print(time2-time1)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

输出结果：

RMSE: 0.8724
RMSE: 0.8737
RMSE: 0.8722
user: 196 item: 302 r_ui = 4.00 est = 4.57 {‘was_impossible’: False}
145.8252305984497

(2) 使用biased SVD：

time1=time.time()
# 使用 Biased SVD
algo = SVD(biased = True)

# 定义K折交叉验证迭代器，K=3
kf = KFold(n_splits=3)
for trainset, testset in kf.split(data):
    # 训练并预测
    algo.fit(trainset)
    predictions = algo.test(testset)
    # 计算RMSE
    accuracy.rmse(predictions, verbose=True)

uid = str(196)
iid = str(302)
# 输出uid对iid的预测结果
pred = algo.predict(uid, iid, r_ui=4, verbose=True)
time2=time.time()
print(time2-time1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

输出结果：

RMSE: 0.8465
RMSE: 0.8432
RMSE: 0.8462
user: 196 item: 302 r_ui = 4.00 est = 3.97 {‘was_impossible’: False}
141.18379878997803

(3) 使用 SVD++ ：

from surprise import SVDpp

time1=time.time()
# 使用 Biased SVD
algo = SVDpp(n_factors = 5)

# 定义K折交叉验证迭代器，K=3
kf = KFold(n_splits=3)
for trainset, testset in kf.split(data):
    # 训练并预测
    algo.fit(trainset)
    predictions = algo.test(testset)
    # 计算RMSE
    accuracy.rmse(predictions, verbose=True)

uid = str(196)
iid = str(302)
# 输出uid对iid的预测结果
pred = algo.predict(uid, iid, r_ui=4, verbose=True)
time2=time.time()
print(time2-time1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

输出结果：

RMSE: 0.8308
RMSE: 0.8293
RMSE: 0.8309
user: 196 item: 302 r_ui = 4.00 est = 3.90 {‘was_impossible’: False}
4811.0583329200745