主要参考论文《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 rated
double 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);
else
BIASU[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);
else
BIASI[j] = 0;
}
}
void train()
{
//read rating matrix
memset(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 decomposation
for (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;
}
效果图如下：
可以看到，rmse能非常快的收敛，训练数据中的rmse能很快收敛到0.8左右，然后拿测试集的数据去测试，rmse为0.949，也是蛮不错的预测结果了，当然这里可以调各种参数来获得更优的实验结果。。。就是所谓的黑科技？。。。：）