矩阵分解：从入门到高级入门

发表于 | 分类于

引言

矩阵分解可以视为无监督学习的一种形式，成名于Netflix Prize推荐比赛，由Yehuda Koren引申出诸多变种，在推荐系统中可谓是简单粗暴，非常有效。

矩阵分解，顾名思义，就是把已有的目标矩阵分解为两个矩阵相乘的形式，那么目标矩阵是什么？在推荐系统里一般我们可以理解为是用户（user）-商品（item）-评分（rating）所构成的矩阵。矩阵的每一行代表一个user，每一列代表item，矩阵元素$r_{ij}$代表第i个user为第j个item打得分数，可以是浏览数、加购数、购买数、评论数等等，如果第i个user对第j个item没有任何行为，我们令该矩阵元素为空。而我们的目标就是将这些空元素补齐，用以推断每个user对没有行为的这些item的兴趣度大小究竟是多少，从而找到用户可能最感兴趣的top-K个商品进行推荐。

目标函数

假设矩阵$R$维度为$M \times N$，也就是用户为$M$个，item为$N$个，目标是将$R$分解为$X \cdot Y$的形式，其中$X$维度为$M \times K$，$Y$维度为$K \times N$，$K << M,N$，这里的$K$称为隐因子维度。我们可以把user的隐因子解释为用户的兴趣维度，把item的隐因子解释为商品的属性维度，分解出的结果同样也可以基于cosine similarity计算用户之间的相似度和商品之间的相似度。其实传统的SVD同样也是这种思想，但是由于这里我们的场景有大量的缺失值SVD无法操作，盲目的填补缺失值也会造成模型的过拟合，因此也衍生出了大量诸如SVD++，RSVD等SVD加强版。事实上最基础和原始的矩阵分解优化目标形式非常简洁，思想就是尽可能的可以复原出已有矩阵中非空元素，使预测元素与原始元素保持一致，即

当然，还少不了regulation term，以防止overfitting，这样目标就变得复杂一点

Bias

有了这个基础版本的目标函数，其实就可以在上面进行修修改改，逐步完善。考虑到在实际推荐系统中user和item的多样性，我们还可以在此基础上加上user、item的bias项。举个直观例子来解释这个概念，假如现在要预测用户A对商品a的打分，已经知道的是在所有的商品里均值得分为4分，而a相比所有商品来说质量和口碑要相对好一些，因此均值要高出1分，但用户A又是一个比较较真和挑剔的用户，所以他对每个商品的打分要比所有人打分均值低0.5分，因此如果暂且不考虑兴趣因素，最终A对a的打分可以粗略估计为4+1-0.5=4.5，这里这些均值偏移数据就是所谓的bias项修正。如此目标就又复杂一些

Implicit Feedback

上面说了一大堆，但是有比较核心的一点被忽略了，那就是在现实的应用场景中，“rating”，也就是打分这个概念我们几乎是获取不到的。这相当于是用户给出最积极的反馈，比如豆瓣电影的主动打分，但是大多数用户在网站网行的行为也不过是点击、浏览，在电商场景中可能更多一些，比如关注、加购等，在一些音乐、影视平台上可能有收听、观看和评论等。但总的说来，这些都不能够让我们具体的概化为一个“打分”，或者称为显示反馈，所以需要引入隐式反馈这个概念来进一步刻画矩阵分解的优化目标。

首先引入一个二值变量，定义

这里的$r_{ui}$就不再是打分数据了，而是用户对商品有过任意行为的次数，而二值变量可以理解为一个用户对一个商品有过是/否有过任何行为。接下来优化目标也随之改变

，其中$c_{ui}$称为置信度，通常定义为

$\alpha$与$\gamma$同样都是hyper parameter，需要交叉验证来确定，对置信度定性的解释也就是用户对指定商品的行为越多，我们就越有理由认为用户对这个商品是真正感兴趣的，比较类似于统计检验中置信度的概念。设想我们在网购时买一个自己喜欢的东西基本上都要看来看去好多次甚至好多天，但如果是帮人买，多数情况下的行为是与这个商品只接触一次，买完即走。更generalized的想法是，把$c_{ui}$解释为用户对商品的行为权重，对电商而言，如果仅仅浏览过，权重最小，浏览次数大于一定阈值时，权重次之，加购过，权重再次加大，购买过，权重最大。

如果考虑得再复杂一点，可以把$x_u^T y_i$更加具体化，比如加上用户的属性特征$U$，有过行为的商品embedding归一化特征$I$，扩充用户的隐因子矩阵$x$，就变为

对于$y_i$，如果有必要的话也可以同样进行这种扩充方式处理，从而解决一些数据不足或冷启动的问题。

学习算法

最小化目标函数的方法主要有梯度下降法和交替最小二乘法（Alternative Least Square, ALS），接下来分别简要说明一下。

梯度下降

其实严格的来说，矩阵分解的梯度下降应该叫做随机梯度下降（逐个元素更新）。先看最基础的，定义损失函数

，待更新参数为$x_{ik}$和$y_{kj}$，每一步梯度更新后两个参数分别为

加上正则项后，形式类似，

这里有一点需要注意，并不是$X$和$Y$中的所有元素都需要按顺序依次更新，因为目标矩阵的稀疏性导致大量rating值缺失，所以只需要更新那些$e_{ij}$不为空对应的user和item矩阵元素。

python的简单实现
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
import numpy
def matrix_factorization(R, X, Y, rank, iter, lr = 0.01, gamma =0.02):
    Y = Y.T
    # update user-item latent factor
    for step in xrange(iter):
        for i in xrange(R.shape[0]):
            for j in xrange(R.shape[1]):
                if R[i][j] > 0:
                    eij = R[i][j] - numpy.dot(X[i,:],Y[:,j])
                    for k in xrange(rank):
                        X[i][k] = X[i][k] + 2 * lr * (eij * Y[k][j] - gamma * X[i][k])
                        Y[k][j] = Y[k][j] + 2 * lr * (eij * X[i][k] - gamma * Y[k][j])
    return X, Y.T
R = numpy.array([[5,3,0,1],
     [4,0,0,1],
     [1,1,0,5],
     [1,0,0,4],
     [0,1,5,4]])
rank = 2 
X = numpy.random.rand(R.shape[0], rank)
Y = numpy.random.rand(R.shape[1], rank)
X_new, Y_new = matrix_factorization(R, X, Y, rank, 1000)
R_new = numpy.dot(X_new, Y_new.T)
print "predict matrix:\n", R_new	
print "raw matrix:\n", R
R = R.reshape(-1)
R_new = R_new.reshape(-1)
mse = 0
for i in xrange(R.size):
    if R[i] > 0:
        mse += pow(R[i] - R_new[i], 2) 
print "MSE: ", mse/R.size

输出结果

1
2
3
4
5
6
7
8
9
10
11
12
13
predict matrix:
[[ 4.95397969  2.96469808  4.39008599  1.00365562]
 [ 3.96264341  2.38721232  3.71827269  1.00055567]
 [ 1.00646214  0.98090633  5.85147223  4.94890212]
 [ 0.99598939  0.8966756   4.82086649  3.96964276]
 [ 1.2039566   1.01871901  4.97353667  3.98151943]]
raw matrix:
[[5 3 0 1]
 [4 0 0 1]
 [1 1 0 5]
 [1 0 0 4]
 [0 1 5 4]]
MSE:  0.000506024533095

ALS

从目标函数$e_{ij}^2$的形式上可以看出来，这是个non-convex function，但却是bi-convex的。也就是说，在已知$x$的情况下，就是个关于$y$的quadratic形式的function，同理在已知$y$的情况下，就是个关于$x$的quadratic形式的function，而quadratic形式是一定能找到全局最优的。针对这种问题，我们就可以固定其中一个变量来更新另一个变量，然后做一次相反的操作，迭代多次，就可以获得最下二乘下的优化解，所以这种方法就称之为交替最小二乘。正是由于在每一次迭代都可以获得真正意义的全局最优，而非梯度下降方法每次迭代只在梯度方向移动一小步，因此ALS在理论上要比gradient descent方法收敛快很多。

开始推导，这里我们只针对带正则项的explicit feedback的目标函数形式，

先固定$y_i$，求$x_u$，求导有

令上式等于0，直接求最优，可得$x_u$表达式

简单解释下各个符号含义。假设用户维度为$N$维，商品维度为$M$维，隐因子的rank为$K$维，用户$u$打过分的商品为$m$个，那么$Y_{I_i}$代表用户$u$打过分的item集合的向量堆叠，即$K \times m$维，$R(u, I_i)$用户$u$打过分的商品得分，为一个$m*1$列向量，$I$为单位向量，维度为$K \times K$，$x_u$即为用户$u$的向量表示，是一个$k \times 1$的列向量。在更新完$x_u$后，同理可以更新$y_i$，由于目标函数是对称的，就不推导了，直接上公式

含义和上面的刚好相反。如果商品$i$被$n$个用户打过分，那么$X_{I_u}$代表商品i被打过分的user集合的向量堆叠，即$K \times n$维，$R(i, I_u)$为商品$i$被打过的得分，为一个$n \times 1$的列向量，$I$为单位向量，维度为$K*K$，$y_i$即为商品$i$的向量表示，是一个$k \times 1$的列向量。在实际应用里，由于矩阵的求逆操作计算量很大，通常就通过解线性方程组的方法来求$x_u$和$y_i$。

python的简单实现
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
41
42
43
import numpy as np
def matrix_factorization(R, rank, iter, gamma =0.02):
	#initialization
	X = np.random.rand(R.shape[0], rank)
	Y = np.random.rand(R.shape[1], rank)
	I = np.eye(rank)
	for step in xrange(iter):
		for user, rating in enumerate(R):
			is_rating_index = np.nonzero(rating)
			A = np.dot(Y[is_rating_index].T, Y[is_rating_index]) + gamma * I
			B = np.dot(Y[is_rating_index].T, rating[is_rating_index])
			# using solve method instead of calculating inverse matrix
			X[user,:] = np.linalg.solve(A, B)
		for item, rating in enumerate(R.T):
			is_rating_index = np.nonzero(rating)
			A = np.dot(X[is_rating_index].T, X[is_rating_index]) + gamma * I
			B = np.dot(X[is_rating_index].T, rating[is_rating_index])
			Y[item,:] = np.linalg.solve(A, B)
	return X, Y
R = np.array([[5,3,0,1],
     [4,0,0,1],
     [1,1,0,5],
     [1,0,0,4],
     [0,1,5,4]])
rank = 2 
X, Y = matrix_factorization(R, rank, 10)
R_new = X.dot(Y.T)
print "predict matrix:\n", 	R_new
print "raw matrix:\n", R
R = R.reshape(-1)
R_new = R_new.reshape(-1)
mse = 0
for i in xrange(R.size):
    if R[i] > 0:
        mse += pow(R[i] - R_new[i], 2) 
print "MSE: ", mse/R.size

输出结果

1
2
3
4
5
6
7
8
9
10
11
12
13
predict matrix:
[[ 5.0161259   2.9859823  -3.84990048  1.01292524]
 [ 3.9720155   2.38887852 -2.78909332  0.9972788 ]
 [ 0.96012708  1.16133008  5.5264977   4.90613676]
 [ 1.04836492  1.09127496  4.15699704  3.94456546]
 [ 0.50838622  0.809576    4.99343708  4.15301915]]
raw matrix:
[[5 3 0 1]
 [4 0 0 1]
 [1 1 0 5]
 [1 0 0 4]
 [0 1 5 4]]
MSE:  0.00514865377649

至此，矩阵分解的基本原理和求解方法都已经啰嗦完了，但是现实应用中其实远没有这么简单，比如关于隐反馈的定义（用户-商品打分的量化），大数据量的内存问题（spark ALS处理亿量级的用户）等。总的来说坑还是很多的，所以这篇文章也只能算是从入门到高级入门。

参考资料

Matrix factorization techniques for recommender systems
Large-scale parallel collaborative filtering for the netflix prize
Collaborative filtering for implicit feedback datasets