对于一个大型的视频网站来说，我们需要“好”的内容的视频的权重被迅速提拔起来被更多的人看见，那么就需要一套完整的评分系统，预测评分显得尤为重要了。不过，并不是所有的协同过滤结果都是合理的，比如说，大部分电商平台面临的刷单问题，其实就是欺骗推荐系统获取权重的方法。

这一节，我们主要介绍 SVD，通过矩阵分解的方法来降低数据中的噪点，将数据最终转化成特征向量。

上文链接：

今晚的风儿很喧嚣：推荐算法入门（1）相似度计算方法大全​ zhuanlan.zhihu.com

数据来源：

以下题目和数据均来自于千里码，一个优质的程序员答题网站。

现在从豆瓣的用户中抽取了500左右个比较活跃的用户，这些用户都是忠实的电影迷，大部分人涉猎了上百部电影。

这里有个80多万行的
文本文件，文件的每行是三个数字，分别是 userid，movieid，rating。代表一个用户对一部电影的评分。rating代表评分的星级，如上图中的红框所示，星级从低到高依次是1-5。

接下来有个行数为10001的
文本文件（第一行为title），文件的每行为2个数字，分别代表userid和movieid，请你预测如果该用户观看了这部电影，会给该电影打多少分，你的预测值为1个大小为1-5的整数。

本题的答案是一个长度为1万的字符串，字符串的第k位代表你对第k行的预测结果。

如果你的预测结果和实际答案的差值的绝对值的和小于6000，通过该题。

或许你有更多的方法实现推荐，不来千里码试试吗？

答案提交地址：

http://www.qlcoder.com/task/7650​ www.qlcoder.com

一、矩阵奇异值分解（Singular Value Decomposition）

在很多情况下，数据的一小部分包含了数据的绝大部分信息，用线性代数的许多矩阵分解方法都可以将矩阵特征表现成新的易于理解的形式，在这一过程之中，实际是一个提取数据集合特征值的过程，有利于我们提高算法的精准度，减少数据运算量。

最常见的矩阵分解算法就是矩阵奇异值分解（SVD），SVD 在图像压缩、推荐系统、金融数学等领域都有应用，著名的主成成分分析算法（PCA）也是通过 SVD 实现的。

二、SVD 理论简介

假设数据 M 是一个 m*n 阶的样本矩阵，其中的元素全部属于域 K，那么矩阵分解可得：

描述成： U（m*m），sigma（m*n），VT（n*n）

矩阵 sigma 除了对角元素不为0以外，其他元素都为0，并且对角元素是从大到小顺序排列的，这些对角元素就是奇异值。

所以。对于任意的奇异值分解，矩阵 sigma 的对角线上的元素等于M的奇异值。U和V的列分别是奇异值中的左奇异向量和右奇异向量。

这一段是来源于百度百科的介绍，其实已经解释清楚了 SVD 是啥了，那么 SVD 的三个矩阵如何计算？如果你不想过多了解也没关系，因为单靠原始的 SVD 效率并不是很高，想要手撸 SVD的，这篇文章讲的比我好。

漫漫成长：奇异值分解（SVD）​ zhuanlan.zhihu.com

三、SVD 图像降噪

这是一个图像化说明 svd 的例子，虽然各种意义上面的意义不明。

test.jpg

test_70%.jpg

保留的奇异值越多，图片的特征保留的越明显，当奇异值减少时，图片中的像素间的差距逐渐减小（表现的模糊）。

import numpy as np
from PIL import Image

# 使用奇异值总和的百分比进行筛选
def svd(data,scale):
    # scale 代表你要保留的奇异值比例
    u,sigma,v = np.linalg.svd(data)
    svd_data = np.zeros(data.shape)
    total = sum(sigma)
    sum_data = 0
    for index,item in enumerate(sigma):  
        svd_data += item * np.dot(u[:,index].reshape(-1,1),v[index,:].reshape(1,-1))
        sum_data += item
        if sum_data >= scale * total:
            break
    return svd_data

def compress(data,scale):
    r = svd(data[:,:,0],scale)
    g = svd(data[:,:,1],scale)
    b = svd(data[:,:,2],scale)

    result = np.stack((r,g,b),2)
    result[result > 255] = 255
    result[result < 0] = 0
    result = result.astype(int)
    return result

if __name__ == '__main__':
    image = Image.open('test.jpg')
    width,height = image.size
    arr = np.zeros((width,height,3)) # RGB 
    for x in range(width):
        for y in range(height):
            arr[x][y] = image.getpixel((x, y))
    # 原生 range 不支持浮点数，所以用 np.arange 代替
    for c in np.arange(.1,.9,.2):
        result = compress(arr,c)
        for x in range(width):
            for y in range(height):
                image.putpixel((x, y),(result[x,y,0],result[x,y,1],result[x,y,2]))
        image.save('test_'+str(int(100 * c))+'%.jpg')

四、Netflix Prize

提到协同过滤算法，也提了 svd，就不得不提由美国 Netflix 公司举办的 Netflix Prize，这是一个旨在解决电影评分预测问题的竞赛，因为高昂的奖金，已经成为了协同过滤发展的风尚了。

也许你已经看过了很多的关于 svd 的文章了，也已经悉知了 svd 的原理了，决定向领导打包票开始写高阶协同过滤算法了。然后翻了翻往年的 Netflix Prize，纳尼。。。 我是谁？我在哪？然后可以收拾收拾删库跑路了。推荐算法水之深，我们都是新萌。

Netflix Prize – Wikipedia​ en.wikipedia.org

五、Netflix Prize 大放异彩的 Funk SVD 算法

1. 传统 SVD 分解在元素缺失上面的问题

历史上对缺失值的研究有很多，对于一个没有被打分的物品来说，到底是应该给它补一个 0 值，还是应该给它补一个平均值呢？由于在实际过程中，元素缺失值是非常多的，这就导致了早期的 SVD 不论通过以上哪种方法进行补全在实际的应用之中都是不可以被接受的。

2.LFM （Latent Factor Model）

直到 2006年 Netflix Prize 中 Simon Funk 在博客公开的算法。将评分矩阵分解成两个低维矩阵相乘，Simon Funk的思想很简单：可以直接通过训练集中的观察值利用最小化均方根学习P，Q矩阵。这种模型也被称作是 LFM （隐语义模型），下面是他当年的博客，有兴趣可以了解一下。

Try This at Home​ sifter.org

简单的来说就是将原本的 SVD 的思想上加上了线性回归，也就是说，我们可以用均方差作为损失函数，来寻找 P 和 q 的最终值，线性回归和均方差对于机器学习的同学们来说一定不陌生了，如果你还没有了解过，可能一下子理解不了下面的公式，那么我建议还是先从线性回归学起，便于理解。不过，线性回归也就是一句话 —— 线性函数参数调优。

今晚的风儿很喧嚣：斯坦福大学机器学习课程（介绍和线性回归）​ zhuanlan.zhihu.com

2.加入偏移项后的 Funk-SVD

在 Funk-SVD 获得巨大成功之后，很多著名的模型都是对 Funk-SVD 进行缝缝补补得到的（详情可参见 Netflix Prize `Koren:2009` `Ricci:2010`），于是就有了在预测模型中添加三项偏移的模型，被称为 BaisSVD。

• Biased Item
• Biased User
• Biased Mean

Biased Item（物品偏移），表示了物品接受的评分和用户没有多大关系，物品本身质量决定了的偏移。

Biased User（用户偏移），有些用户喜欢打高分，有些用户喜欢打低分，用户决定的偏移。

Biased Mean（全局平均值偏移），根据网站全局打分设置的偏移，可能和整体用户群和物品质量有相对应的关系。

3.公式说明

符号含义：

rates 学习率

regularization 正则化项

b biased 偏移

矩阵 M 经过Funk SVD 分解之后，只分解出两个矩阵P和Q：

对于某一个用户的评分使用 Funk SVD 进行矩阵分解得到的结果就是：

那么我们需要最小化的损失函数就是：

损失就是：

使用随机梯度下降调参：

还有一种情况就是，抛弃biased，也就是：

Python 代码：

以下是算法核心的随机梯度下降，注意，代码是不完整的，无法正常运行。只是为想要自己动手写写的同学提供一个思路，完整的代码在 github 仓库。

def sgd(self, trainset):
	"""  这里是算法的核心部分，由于需要效率，这里使用了 cpython，你可以找到根目录下的setup.py 进行编译  """
	cdef np.ndarray[np.double_t] bu
	# item biases
	cdef np.ndarray[np.double_t] bi
	# user factors
	cdef np.ndarray[np.double_t, ndim=2] pu
	# item factors
	cdef np.ndarray[np.double_t, ndim=2] qi

	cdef int u, i, f
	cdef double r, err, dot, puf, qif
	cdef double global_mean = self.trainset.global_mean()

	cdef double lr_bu = self.lr_bu
	cdef double lr_bi = self.lr_bi
	cdef double lr_pu = self.lr_pu
	cdef double lr_qi = self.lr_qi

	cdef double reg_bu = self.reg_bu
	cdef double reg_bi = self.reg_bi
	cdef double reg_pu = self.reg_pu
	cdef double reg_qi = self.reg_qi

	# random samples from a normal (Gaussian) distribution.
	rng = np.random.mtrand._rand

	bu = np.zeros(trainset.n_users, np.double)
	bi = np.zeros(trainset.n_items, np.double)
	pu = rng.normal(self.init_mean, self.init_std_dev,
					(trainset.n_users, self.n_factors))
	qi = rng.normal(self.init_mean, self.init_std_dev,
					(trainset.n_items, self.n_factors))

	if not self.biased:
		global_mean = 0

	for current_epoch in range(self.n_epochs):
		if self.verbose:
			print("Processing epoch {}".format(current_epoch))

		for u, i, r in trainset.all_ratings():

			# compute current error
			dot = 0  # <q_i, p_u>
			for f in range(self.n_factors):
				dot += qi[i, f] * pu[u, f]
			err = r - (global_mean + bu[u] + bi[i] + dot)

			# update biases
			if self.biased:
				bu[u] += lr_bu * (err - reg_bu * bu[u])
				bi[i] += lr_bi * (err - reg_bi * bi[i])

			# update factors
			for f in range(self.n_factors):
				puf = pu[u, f]
				qif = qi[i, f]
				pu[u, f] += lr_pu * (err * qif - reg_pu * puf)
				qi[i, f] += lr_qi * (err * puf - reg_qi * qif)

	self.bu = bu
	self.bi = bi
	self.pu = pu
	self.qi = qi

如有错误欢迎指出。关爱新萌，人人有责。

Github仓库链接：

NoisyWinds/Recommend​ github.com

欢迎点赞，欢迎评论一起交流学习，感谢支持。