推荐系统初步窥探

学习篇

承接上篇，抽时间学习下svd在推荐系统中的简单应用，从谷歌上找到了一篇文章，大概1000+的引用，出自2000年,名字为Application of Dimensionality Reduction in Recommender System – A Case Study，虽然时间有点久远，但应该能学到很多东西。

• 文章主要介绍到随着发展的需求，用户和商品的数量增长，推荐系统遭遇到了各种问题，需要新的方法或者技术去提高效率

• 原文做了两个实验，第一个实验是使用不同的推荐系统来预测用户的对某个商品的评分，另外一个是用协议过滤来处理基于真实数据库中，用户购买数据的TOP-N商品。

• 原文章主要涉及到以下几点：

  1. 模型技术细节，以LSI/SVD用于推荐系统降维产生预测值
  2. 在低纬数据集中计算用来推荐的邻居组
  3. 作者的实验结果

当前（文章发表的时间）推荐系统用到的方法与其局限性

大多数CF模型主要是基于相同意愿的邻居，通过邻居的购买习惯来给新用户推荐产品。
邻居组的形成通常是通过皮尔森相关系数（Pearson correlation）或者余弦相似度（cosine similarity）来计算的。所以就有两种类型的推荐。

1. 预测一个用户C对一个产品P的喜爱程度。计算C在P上的喜爱程度，主要是通过计算一个C和C的所有邻居的求和权重，然后再加上C的平均评分。公式为：

$$C_{p_{pred}} = \overline{C}+\dfrac{\sum_{J\in{rates}} (J_p-\overline{J})}{\sum_{J}|r_{CJ}|}$$
这里， $r_{CJ}$代表着用户 $C$与其邻居组中的第$j$个用户的相互关系,比如皮尔森相关系数。 $J_p$是用户 $j$在产品P上的评分，$\overline{J}$和 $\overline{C}$代表着 $J$和$C$的平均评分。这个计算结果得到的是用户 $C$在产品$P$上的个性化预测值。然而也有一些非个性化预测方案，仅仅是简单的计算产品评分的均值，然后将这个值作为新用户的预测值。

2.给用户推产品，也就是所谓的top-N推荐，一旦一个邻居组计算出来了，推荐系统就可以计算这个邻居组的N个可以被用来推荐的产品。

这些系统在多个领域取得了不错的成绩，但是一些相关的算法已经开始表现出了局限性。例如：

• 稀疏矩阵（Sparsity）：最邻近算法主要依靠精确的用户匹配，一般的，不同用户之间至少有2个相同的评分才可以去计算。然后很多用户之间没有这种联系。在实践中，很多商业推荐系统中的活跃用户的评分数目小于产品的数目，有的甚至不到1%，相应的，皮尔森相关系数也许就不能很好的工作，这个问题被称作reduced coverage。由于评分数据不够多，就可以能导致很多很多问题的发生，比如推荐的精确度，还有就是相关性，比如用户A与用户B高度相关，用户B和用户C高度相关，由于评分数据的稀疏度，那么A和C可能并不相关，甚至会出现负相关情况。
• 可扩展性（Scalability）：主要是在随着日益增长的产品和用户数目，导致以网页为基础的推荐系统在邻居组算法上的开销问题。
• 同义词（Synonymy）：主要体现在，在一个真实的环境中，不同的产品名字可能指向相同的物品，但是基于相关性的推荐系统就不能找到这种潜在的联系，也就说认为本应该是同种物品，却因为不同的名字而认为他们是不同的物品。

在CF上应用SVD

由于皮尔森相关系数在大型稀疏矩阵的缺点，所以开始考虑其他的推荐算法。

预测产生

我们从一个用户产品评分的稀疏矩阵开始，把这矩阵成为R。首先得填充这个稀疏
矩阵R，尝试两种不同的方法，一个是以用户的平均评分，另外一种是产品的评价
评分。我们发现产品的评价评分效果要好些。同时考虑两种归一化处理方法：占位符

推荐产生

在第二个实验中，看看在低纬处理邻居组，然后推荐top-N。

实验

实验平台

数据集

按照之前提到的2个实验，第一个实验用到MovieLens的数据去对评估以SVD为姜维基础的预测生产算法。
我们随机选取足够的用户吗，包含100K的评分记录，其中每条记录都是按照<用户,产品,评分>来记录的，我们按照一个比例来分割数据，

基础推荐系统

为了对SVD处理过的预测数据做对比，我们将训练数据输入到CF推荐系统中，在这个推荐今天中，用皮尔森最邻近算法作为衡量不同对象间的距离。
为了这个目的，我们在部署另外一个推荐系统，不同的是这个里面用余弦函数在高纬向量空间去计算不同对间的距离。（cosine-similarity），从而得到一个一个邻居社区，并且返回top-N推荐结果，也就是所谓的CF推荐

评估指标

评估一个推荐系统的好坏有多种方法，因为我们主要关注推荐系统的结果，所以我们只考虑推荐的质量（quality of prediction or recommendation）
事实上，在评估的过程中，有很多处可以阶段评估，比如针对每个用户的邻居社区质量的评估。
这里，我们主要用到讨论两种评估方法。

http://en.wikibooks.org/wiki/Statistics/Summary/Averages/Harmonic_Mean
http://argcv.com/articles/1036.c

我们选用平均绝对误差来作为衡量预测实验，因为平均绝对误差是最常用的，而且比较直观的了解。

推荐评价指标

为了评估top-N推荐系统，我们使用俩个在信息检索（IR）中常用到的指标，召回率（recall）和精确度（precision），我们把产品分成两组：测试组和top-N组。只有产品在俩组中只要出现过，就认为这个产品属于命中数据集（hit set），稍微修改下recall和precision的定义：

• 在推荐系统中，将recall定义为：
  $recall=\frac{size \space of\space hit \space set}{size \space of\space test \space set}=\frac{|test\bigcap topN|}{|test|}$

• 在推荐系统中，将precision定义为：
  $precision=\frac{size \space of\space hit \space set}{size \space of\space topN \space set}=\frac{|test\bigcap topN|}{|N|}$

在自然情况下，这俩个指标经常会出现冲突，例如，增加N会导致recall值变大，但是会导致precision值下降。我们认为这两个指标都很重要，所以用了一个调和均值（Harmonic Average）来连接着俩个指标：

$$F1=\frac{2*Recall*Precision}{Recall+Precision}$$

我们分别为每个用户计算F1的值，然后计算平均数作为最终指标。

实验步骤

预测实验

在矩阵R中，每个实体代表一个评分（范围在1.0-5.0之间），如果用户$i$没有对电影$j$做出评分，则$r_{ij}$为 空(null)。接下来开始实验。

针对每个用户，先计算每个用户的平均评分，每个电影的平均评分，然后将每个空值以其所在的列（也就是电影）的均值替代。然后我们归一化所有的元素$r_{ij}$,用$（r_{ij}-r_i）$来替代 $r_{ij}$，其中$r_{i}$是每个用户的平均评分，也就是行均值。
经过归一化处理后，利用工具，经行SVD分解处理，得到三个矩阵U、S、V。
S是一个包含奇异值的对角矩阵，对角元素按递减顺序排列，$S_{k}$是从矩阵S中，保留最大的K个奇异值，剩余奇异值以0替换。然后计算降维矩阵的平方根以 $U_kS_K^{1/2}和S^{1/2}V'$，然后相乘，得到一个矩阵943乘以1682的矩阵P。同时，行内积$U_k*S_k$ 和列内积 $S_k*V_k$ 给了我们一个预测分值，在新矩阵P中，每个元素 $p_{ij }$中保留每个用户–电影的评分，然后按照归一化的逆过程，即加上$r_i,p_{ij}+r_i$。然后将训练数据集加载到CF-Predict中，计算每个test数据集的评分。比较SVD的MAE和CF-Predict的MAE结果。
我们多次重复这个过程，控制变量为K, 取值K=2，5-21，25，50，100
然后发现当K=14的时候，有最优值。然后固定K=14，改变train/test数据比例x，从0.2变到0.95，步长为0.05.
对于每一个x值，我们重复10次实验，每次实验都选取不同的training/test集，然后计算平均数。

作者选用的数据是100K大小的那个，我选择的是1M大小的，结果跟原文上有很大出入，而且我不能保证自己是充分理解的作者的想法以及思考，只能是按照自己的对原文的理解加以实验。

---

数据的获取

数据的获取，我是直接从grouplens.org上获取的，文章中，作者用的100K的那个数据，我选择的是1M大小的数据。

实验测试

从解压后得到的README.txt可以直观的了解到，在“ratings.dat”文件中，有我的目标数据，里面的文件格式为：

UserID::MovieID::Rating::Timestamp

• UserIDs range between 1 and 6040
• MovieIDs range between 1 and 3952
• Ratings are made on a 5-star scale (whole-star ratings only)
• Timestamp is represented in seconds since the epoch as returned by time(2)
• Each user has at least 20 ratings

这里，我主要关心前三者也就是，UserID，MovieID，Rating。
原始部分数据：

1::1193::5::978300760
1::661::3::978302109
1::914::3::978301968
1::3408::4::978300275
1::2355::5::978824291
1::1197::3::978302268

用python脚本先来处理一下数据。

    # coding=utf-8
    import numpy as np

    dat={} #创建一个字典用来保存数据
    movie_num=3952
    user_num=6040
    m=open('ratings.dat','r')
    for line in m.readlines(): 
        record=line.split('::')
        userid=int(record[0])
        user_ever_movieid=int(record[1])
        the_rate=float(record[2])
        if not dat.has_key(userid):
            dat[userid]=[0]*movie_num #没有评分的作品用0来填充
            dat[userid][user_ever_movieid-1]=the_rate
        else:
            dat[userid][user_ever_movieid-1]=the_rate
    m.close()
    longlist=[]
    for user in xrange(user_num):
        longlist.extend(dat[user+1])
    sparse_tag=np.array(longlist,dtype=float).reshape(user_num,movie_num) #将超长array重构成二维数组，当然也可以转化成一个稀疏矩阵。
    sparse_tag.dump('output.mat')

这样当前目录下就有了一个output.mat的文件了，但是新文件的大小却比原始的数据集合大了很多。

为了每次读取数据的方便，暂时就先这样。
有了数据之后，先自己搭建一个以皮尔森相关系数为核心的简易CF系统。
命名为：benchmark_system.py。
主要思路比较简单，是以user-item为基础作为评分的。
直接利用scipy包中的皮尔森函数。可以自己定义函数做计算。
随机选取trian和test的集合：

    def choice_train_and_test_database(self,ratios=0.95): #按照比例从原始数据中分割train集合和test集合
            n=self.data_set.shape[0]
            total=int(n*ratios)
            len_list=range(n)
            np.random.shuffle(len_list)
            self.train_set_index=len_list[:total]
            self.test_set_index=len_list[total:]
            self.train_test,self.test_set=self.data_set
            [self.train_set_index,:],self.data_set[self.test_set_index,:]
            return True

从test集合中选取用户，然后从train集合中选取符合要求的邻居,这里面遇到这样的一个问题，比如针对用户C，找到了C的邻居A和B，然后只有B在对物品p上有过评分，而A没有在物品p上评分，尽管A和C用户的相关系数高度相关。

    def comput_neighbors(self,currentursrid,sim=sim):
        '利用皮尔森相关系数计算邻近邻居，返回一个以索引为key的字典'
        useridlist={}
        current_data=self.data_set[currentursrid]  #取出当前用户的历史评分
        # index_list=np.where([~np.isnan(self.data_set[currentursrid])])[1]
        index_list=np.where(current_data!=0)[0] #用当前用户的哪些有效列做为就按邻居的依据
        self.index_list=index_list
        nn_matrix=self.train_test[:,index_list] #从train数据中构建全局新的矩阵
        trian_list=self.get_train_set_index()
        for x in range(len(trian_list)):
            temp=self.sim(current_data[index_list],nn_matrix[x])
            if temp>self.k:
                useridlist[trian_list[x]]=temp
        self.global_neighbor[currentursrid]=sorted(useridlist.iteritems(),key=lambda x:x[1],reverse=True)#当前userID有哪些满足条件的邻居，保存到实例的全局字典中
        return sorted(useridlist.iteritems(),key=lambda x:x[1],reverse=True)

有了邻居的评分后就可以尝试预测了。

        def predict_scores(self,user_id,product_id):#预测评分，预测user_id的用户在产品product_id上的评分
            '''
            user_id代表需要经行评分的用户id，
            nn_dict包含前N个邻居,key为id号，value为皮尔森相关系数
            index代表需要计算的是哪些物品，
            data_matrix代表完整的矩阵，
            product_id代表需要用户C在产品P上预测评分
            '''
            if self.global_neighbor.has_key(user_id):
                self.all_nn_dict=self.global_neighbor[user_id]
            else:
                self.all_nn_dict=self.comput_neighbors(user_id)
            self.nn_dict=dict(self.all_nn_dict[:self.neihbors])#选取满足条件的前N个邻
            C=None
            if self.global_C_average.has_key(user_id):
                C=self.global_C_average[user_id]
            else:
                self.global_C_average[user_id]=self.r_i(self.data_set[[user_id]])  #用户C的平均评分
                C=self.global_C_average[user_id]
            J=self.r_i(self.data_set[self.nn_dict.keys()])#满足条件邻居的平均评分
            J_p=self.data_set[self.nn_dict.keys(),product_id] #满足条件邻居在P上的评分 
            valist=[]
            for m in self.nn_dict.keys(): #选取真正对评分有贡献的邻居
                if self.data_set[m,product_id]!=0:
                    valist.append(m)
            if len(valist):
                J=self.r_i(self.data_set[valist])  #用户邻居组的平均 评分
                J_p=self.data_set[valist,product_id] #用户邻居组的在P上的评分
                r_cj=np.array([self.nn_dict[score] for score in valist],dtype=float)  #邻居与用户C之间的皮尔森相关系数
                C_p=C+((J_p-J)*r_cj).sum()/linalg.norm(r_cj,1) #用户C在未知物品P上的评分预测
            else:#如果所有邻居对物品的评分没有贡献，则返回一个中间值
                C_p=3.0

            return C_p

完整代码如下：

    # coding=utf-8
    import numpy as np
    from scipy.stats.stats import pearsonr
    from scipy import linalg
    from multiprocessing import Pool
    import time
    class CF_Predict():
        '''
        all_nn_dict包含所有邻居
        nn_dict包含前N个满足条件的邻居,key为id号，value为皮尔森相关系数
        neihbors为N的大小
        index_list 表示需要在哪些产品上去计算all_nn_dict。
        '''
        all_nn_dict=None
        nn_dict=None
        neihbors=0
        index_list=None
        k=0
        train_test=None
        test_set=None
        test_set_index=[]
        train_set_index=[]
        global_neighbor={}
        global_C_average={}
        def __init__(self,arg=None,k=0.1,neihbors=30):#实例初始化的时候，设置邻居大小，读取数据到内存，以及相关系数K的选择
            self.all_nn_dict=None
            self.neihbors=neihbors                    
            self.data_set=self.load_data() 
            self.k=k

        def load_data(self,path=None):#读取样本数据
            return np.load('output.mat')

        def sim(self,a,b):              #计算皮尔森相关系数
            return pearsonr(a,b)[0]

        def r_i(self,R=None):           #计算数组中有效的行平均数
            ret=[]
            if R is None:
                r_i=self.data_set.mean(1)
                return r_i
            else:
                for row in R:
                    ret.append(row[np.where(row!=0)[0]].mean())  
                return  np.array(ret)

        def comput_neighbors(self,currentursrid,sim=sim):
            '利用皮尔森相关系数计算邻近邻居，返回一个以索引为key的字典'
            useridlist={}
            current_data=self.data_set[currentursrid]  #取出当前用户的历史评分
            index_list=np.where(current_data!=0)[0] #用当前用户的哪些有效列做为就按邻居的依据
            self.index_list=index_list
            nn_matrix=self.train_test[:,index_list] #从train数据中构建全局新的矩阵
            trian_list=self.get_train_set_index()
            for x in range(len(trian_list)):
                temp=self.sim(current_data[index_list],nn_matrix[x])
                if temp>self.k:
                    useridlist[trian_list[x]]=temp
            self.global_neighbor[currentursrid]=sorted(useridlist.iteritems(),key=lambda x:x[1],reverse=True)#当前userID有哪些满足条件的邻居，保存到实例的全局字典中
            return sorted(useridlist.iteritems(),key=lambda x:x[1],reverse=True)

        def predict_scores(self,user_id,product_id):#预测评分，预测user_id的用户在产品product_id上的评分
            '''
            user_id代表需要经行评分的用户id，
            nn_dict包含前N个邻居,key为id号，value为皮尔森相关系数
            index代表需要计算的是哪些物品，
            data_matrix代表完整的矩阵，
            product_id代表需要用户C在产品P上预测评分
            '''
            if self.global_neighbor.has_key(user_id):
                self.all_nn_dict=self.global_neighbor[user_id]
            else:
                self.all_nn_dict=self.comput_neighbors(user_id)
            self.nn_dict=dict(self.all_nn_dict[:self.neihbors])#选取满足条件的前N个邻
            C=None
            if self.global_C_average.has_key(user_id):
                C=self.global_C_average[user_id]
            else:
                self.global_C_average[user_id]=self.r_i(self.data_set[[user_id]])  #用户C的平均评分
                C=self.global_C_average[user_id]
            J=self.r_i(self.data_set[self.nn_dict.keys()])#满足条件邻居的平均评分
            J_p=self.data_set[self.nn_dict.keys(),product_id] #满足条件邻居在P上的评分

            valist=[]
            for m in self.nn_dict.keys():
                if self.data_set[m,product_id]!=0:
                    valist.append(m)
            if len(valist):
                J=self.r_i(self.data_set[valist])  #用户邻居组的平均 评分
                J_p=self.data_set[valist,product_id] #用户邻居组的在P上的评分
                r_cj=np.array([self.nn_dict[score] for score in valist],dtype=float)  #邻居与用户C之间的皮尔森相关系数
                C_p=C+((J_p-J)*r_cj).sum()/linalg.norm(r_cj,1) #用户C在未知物品P上的评分预测
            else:
                C_p=3.0
            return C_p

        def choice_train_and_test_database(self,ratios=0.95): #按照比例从原始数据中分割train集合和test集合
            n=self.data_set.shape[0]
            total=int(n*ratios)
            len_list=range(n)
            np.random.shuffle(len_list)
            self.train_set_index=len_list[:total]
            self.test_set_index=len_list[total:]
            self.train_test,self.test_set=self.data_set[self.train_set_index,:],self.data_set[self.test_set_index,:]
            return True

        def get_test_set(self):
            return self.test_set

        def get_train_set(self):
            return  self.train_test

        def get_test_set_index(self):
            return self.test_set_index

        def get_train_set_index(self): 
            return self.train_set_index

        def auto_compute(self):

            test_set_list=self.get_test_set_index()
            test_set=self.get_test_set()
            predict_scores=[]
            count_0=0 #本次实验中，一共做了多少次的评分


            for n in test_set_list:
                # print "global n:",n
                valid_item=np.where(self.data_set[n]!=0)#计算当前user在哪些产品上有评分

                for valid_id in valid_item[0]:
                    count_0+=1
                    # print "item:",valid_id
                    temp_score=self.predict_scores(user_id=n,product_id=valid_id)#预测用户n在产品上valid_id上的一次评分
                    predict_scores.append(temp_score-self.data_set[n][valid_id])#计算预测评分与真实评分的差值，用一个列表保存结果


            print 'count:',count_0
            return np.linalg.norm(predict_scores,1)/count_0
    def main():
        cf=CF_Predict(neihbors=45,k=0.1)
        cf.choice_train_and_test_database()
        print cf.auto_compute() 

    if __name__ == "__main__":
        #start = time.clock()
        elapsed = (time.clock() - start)
        #print("Time used:",elapsed)

整个代码跑起来比较慢，并且跑起来的时候，CPU平没有完全利用，跑下来一看，只有25%的利用率，单进程，谷歌搜了下，决定引入multiprocessing模块，并行执行，将CPU完全利用起来，在python2.7的版本中，需要并行执行的函数，得单独放到类外面，放到类里面的话有问题（可以放到类里面，但是传参数的时候，只能穿一个参数。而且需要重写__call__函数，我没折腾出来）。因为另一台机器有8核，所以决定并发8个进程。需要处理函数定义如下,放到类外面:

    def sig_core(self):
        predict_scores=[]
        count_0=0
        valid_item=np.where(self[0].data_set[self[1]]!=0)
        for valid_id in valid_item[0]:
            count_0+=1
            temp_score=self[0].predict_scores(user_id=self[1],product_id=valid_id)
            predict_scores.append(temp_score-self[0].data_set[self[1]][valid_id])
        return float(np.linalg.norm(predict_scores,1)),float(count_0)

self是传进一个类实例，因为需要用到同一个实例中的属性和方法。
同时在mian（）中重新写,因为ratios的长度是16,一共要跑16轮，每一轮中在计算user的MAE时候用8个进程并行处理计算，以提高速度：

    def main():    
        cf=CF_Predict(neihbors=45,k=0.1)
        pool=Pool(processes=8)
        ratios=np.arange(0.2,1,0.05)
        for x in ratios:
            cf.choice_train_and_test_database(ratios=x)
            arg=cf.auto_compute()
            result=pool.map(sig_core,arg)
            a,b= np.array(result).sum(0)
            f=open('result.txt','a')
            f.write(str(a)+' '+str(b)+' '+str(a/b)+'\n')
            f.close()

同时类中的auto_compute()也改写下:

    def auto_compute(self):
        test_set_list=self.get_test_set_index()
        arg=[[self,x] for x in test_set_list]
        return arg

跑完大概用了90分钟左右。
结果如下：

train/test:0.2  MAE:0.733111839491
train/test:0.25  MAE:0.737595316004
train/test:0.3  MAE:0.730551987283
train/test:0.35  MAE:0.733572669895
train/test:0.4  MAE:0.73522774335
train/test:0.45  MAE:0.732973165653
train/test:0.5  MAE:0.728804514563
train/test:0.55  MAE:0.726490568005
train/test:0.6  MAE:0.728009338323
train/test:0.65  MAE:0.726796155946
train/test:0.7  MAE:0.731139866798
train/test:0.75  MAE:0.71820201418
train/test:0.8  MAE:0.722951850975
train/test:0.85  MAE:0.730886284679
train/test:0.9  MAE:0.73905248
train/test:0.95  MAE:0.710399727915

画个建议图，方便视觉：

另一方面，用SVD处理train数据，将train中的用户数据投射到一个低维空间，然后将test数据也投射到这个空间中，然后去寻找邻居，由于从高维降到了低维，所以在时间和空间上，计算开销得到了很明显的节约(在K值已经确定的情况下)。
不过在跑程序之前，先把皮尔森相关系数计算的函数得手动写一遍，因为按照scipy包中提供的pearsonr函数，整个程序跑下来，发现一半的时间都花费在这个函数里面了，应该是完整的函数考虑的东西比较健壮和健全，所以在里面执行的计算比较多，为了使计算的时间更见快点，可以按照公式，写个简易的皮尔森系数计算：

    def my_pearsonr(a,b):
        l=len(a)
        a_sum=sum(a)
        b_sum=sum(b)
        a_average=a_sum/l
        b_average=b_sum/l

        up=0.0  #分子计算
        d_l=0.0
        d_r=0.0
        for x in xrange(l):    
            m=a[x]-a_average
            n=b[x]-b_average
            up+=m*n
            d_l+=m*m    #分母计算
            d_r+=n*n    #分母计算
        return up/(math.sqrt(d_l*d_r))

加快计算还有其他的办法，比如可以利用GPU对计算做加速。nvida官方也发布了支持Python的cuda编程库。）:我折腾了一段时间，没折腾明白，等有时间了，在整理整理GPU上跑python代码。
有了邻居后，就可以根据相关的数学关系式去计算其他数据了。