主要参考论文《A Guide to Singular Value Decomp osition for Collab orative Filtering》
其实一开始是比较疑惑的，因为一开始没有查看论文，只是网上搜了一下svd的概念和用法，搜到的很多都是如下的公式：
其中假设C是m*n的话，那么可以得到三个分解后的矩阵，分别为m*r,r*r,r*n,这样的话就可以大大降低存储代价，但是这里特别需要注意的是：这个概念一开始是用于信息检索方面的，它的C矩阵式完整的，故他们可以直接把这个矩阵应用svd分解，但是在推荐系统中，用户和物品的那个评分矩阵是不完整的，是个稀疏矩阵，故不能直接分解成上述样子，也有的文章说把空缺的补成平均值等方法，但效果都不咋滴。。。在看了上述论文之后，才明白，在协同过滤中应用svd，其实是个最优化问题，我们假设用户和物品之间没有直接关系，但是定义了一个维度，称为feature，feature是用来刻画特征的，比如描述这个电影是喜剧还是悲剧，是动作片还是爱情片，而用户和feature之间是有关系的，比如某个用户必选看爱情片，另外一个用户喜欢看动作片，物品和feature之间也是有关系的，比如某个电影是喜剧，某个电影是悲剧，那么通过和feature之间的联系，我们可以把一个评分矩阵rating = m*n(m代表用户数，n代表物品数)分解成两个矩阵的相乘：user_feature*T(item_feature), T()表示转置，其中user_feature是m*k的（k是feature维度，可以随便定），item_feature是n*k的，那么我们要做的就是求出这两者矩阵中的值，使得两矩阵相乘的结果和原来的评分矩阵越接近越好。
这里所谓的越接近越好，指的是期望越小越好，期望的式子如下：
其中n表示用户数目，m表示物品数目，I[i][j]是用来表示用户i有没有对物品j评过分，因为我们只需要评过分的那些越接近越好，没评过的就不需要考虑，Vij表示训练数据中给出的评分，也就是实际评分，p(Ui,Mj)表示我们对用户i对物品j的评分的预测，结果根据两向量点乘得到，两面的两项主要是为了防止过拟合，之所以都加了系数1/2是为了等会求导方便。
这里我们的目标是使得期望E越小越好，其实就是个期望最小的问题，故我们可以用随机梯度下降来实现。随机梯度下降说到底就是个求导问题，处于某个点的时候，在这个点上进行求导，然后往梯度最大的反方向走，就能快速走到局部最小值。故我们对上述式子求导后得：
所以其实这个算法的流程就是如下过程：
实现起来还是比较方便快捷的，这里rmse是用来评测效果的，后面会再讲。
上述算法其实被称为批处理式学习算法，之所以叫批处理是因为它的期望是计算整个矩阵的期望（so big batch），其实还存在增量式学习算法，批处理和增量式的区别就在于前者计算期望是计算整个矩阵的，后者只计算矩阵中的一行或者一个点的期望，其中计算一行的期望被称为不完全增量式学习，计算一个点的期望被称为完全增量式学习。
不完全增量式学习期望如下（针对矩阵中的一行的期望，也就是针对一个用户i的期望）：
那么求导后的式子如下：
算法的大致思想如下：
完全增量式学习算法是对每一个评分进行期望计算，期望如下：
求导后如下：
所以整个算法流程是这样的：
上述都是svd的变种，只不过实现方式不一样，根据论文所说，其中第三种完全增量式学习算法效果最好，收敛速度非常快。
当然还有更优的变种，考虑了每个用户，每个物品的bias,这里所谓的bias就是每个人的偏差，比如一个电影a,b两人都认为不错，但是a评分方便比较保守，不错给3分，，b评分比较宽松，不错给4分，故一下的评分方式考虑到了每个用户，每个物品的bias，要比上述算法更加精准。原来评分的话是直接计算user_feature*T(item_feature), T()表示转置，但现在要考虑各种bias，如下：
其中a表示所有评分的平均数，ai表示用户i的bias，Bj表示物品j的偏差，相乘的矩阵还是和上面一样的意思。
故这时候的期望式子和求导的式子如下（这里只写了bias的求导，矩阵求导还是和上面一样）：
当然了，光说不练假把式，我们选择了最后一种算法，及考虑bias的算法来实现了一把，数据源是来自movielens的100k的数据，其中包含了1000个用户对2000件物品的评分（当然，我这里是直接开的数组，要是数据量再大的话，就不这么实现了，主要是为了验证一把梯度下降的效果），用其中的base数据集来训练模型，用test数据集来测试数据，效果评测用一下式子来衡量：
说白了就是误差平方和。。。同时我们也记录下每一次迭代后训练数据集的rmse.
代码如下：
#include <iostream>#include <string>#include <fstream>#include <math.h>using namespace std;const int USERMAX = 1000;const int ITEMMAX = 2000;const int FEATURE = 50;const int ITERMAX = 20;double rating[USERMAX][ITEMMAX];int I[USERMAX][ITEMMAX];//indicate if the item is rateddouble UserF[USERMAX][FEATURE];double ItemF[ITEMMAX][FEATURE];double BIASU[USERMAX];double BIASI[ITEMMAX];double lamda = 0.15;double gamma = 0.04;double mean;double predict(int i, int j){double rate = mean + BIASU[i] + BIASI[j];for (int f = 0; f < FEATURE; f++)rate += UserF[i][f] * ItemF[j][f];if (rate < 1)rate = 1;else if (rate>5)rate = 5;return rate;}double calRMSE(){int cnt = 0;double total = 0;for (int i = 0; i < USERMAX; i++){for (int j = 0; j < ITEMMAX; j++){double rate = predict(i, j);total += I[i][j] * (rating[i][j] – rate)*(rating[i][j] – rate);cnt += I[i][j];}}double rmse = pow(total / cnt, 0.5);return rmse;}double calMean(){double total = 0;int cnt = 0;for (int i = 0; i < USERMAX; i++)for (int j = 0; j < ITEMMAX; j++){total += I[i][j] * rating[i][j];cnt += I[i][j];}return total / cnt;}void initBias(){memset(BIASU, 0, sizeof(BIASU));memset(BIASI, 0, sizeof(BIASI));mean = calMean();for (int i = 0; i < USERMAX; i++){double total = 0;int cnt = 0;for (int j = 0; j < ITEMMAX; j++){if (I[i][j]){total += rating[i][j] – mean;cnt++;}}if (cnt > 0)BIASU[i] = total / (cnt);elseBIASU[i] = 0;}for (int j = 0; j < ITEMMAX; j++){double total = 0;int cnt = 0;for (int i = 0; i < USERMAX; i++){if (I[i][j]){total += rating[i][j] – mean;cnt++;}}if (cnt > 0)BIASI[j] = total / (cnt);elseBIASI[j] = 0;}}void train(){//read rating matrixmemset(rating, 0, sizeof(rating));memset(I, 0, sizeof(I));ifstream in("ua.base");if (!in){cout << "file not exist" << endl;exit(1);}int userId, itemId, rate;string timeStamp;while (in >> userId >> itemId >> rate >> timeStamp){rating[userId][itemId] = rate;I[userId][itemId] = 1;}initBias();//train matrix decomposationfor (int i = 0; i < USERMAX; i++)for (int f = 0; f < FEATURE; f++)UserF[i][f] = (rand() % 10)/10.0 ;for (int j = 0; j < ITEMMAX; j++)for (int f = 0; f < FEATURE; f++)ItemF[j][f] = (rand() % 10)/10.0 ;int iterCnt = 0;while (iterCnt < ITERMAX){for (int i = 0; i < USERMAX; i++){for (int j = 0; j < ITEMMAX; j++){if (I[i][j]){double predictRate = predict(i, j);double eui = rating[i][j] – predictRate;BIASU[i] += gamma*(eui – lamda*BIASU[i]);BIASI[j] += gamma*(eui – lamda*BIASI[j]);for (int f = 0; f < FEATURE; f++){UserF[i][f] += gamma*(eui*ItemF[j][f] – lamda*UserF[i][f]);ItemF[j][f] += gamma*(eui*UserF[i][f] – lamda*ItemF[j][f]);}}}}double rmse = calRMSE();cout << "Loop " << iterCnt << " : rmse is " << rmse << endl;iterCnt++;}}void test(){ifstream in("ua.test");if (!in){cout << "file not exist" << endl;exit(1);}int userId, itemId, rate;string timeStamp;double total = 0;double cnt = 0;while (in >> userId >> itemId >> rate >> timeStamp){double r = predict(userId, itemId);total += (r – rate)*(r – rate);cnt += 1;}cout << "test rmse is " << pow(total / cnt, 0.5) << endl;}int main(){train();test();return 0;}效果图如下：
命运如同手中的掌纹，无论多曲折，终掌握在自己手中。
