SVD/SVD++實現推薦算法

奇異值分解(Singular Value Decomposition，如下簡稱SVD)是在機器學習領域普遍應用的算法，它不只能夠用於降維算法中的特徵分解，還能夠用於推薦系統，以及天然語言處理等領域。

咱們首先回顧下特徵值和特徵向量的定義以下：

                                                                    Ax=λx

其中A是一個n×n的矩陣，x是一個n維向量，則咱們說λ是矩陣A的一個特徵值，而x是矩陣A的特徵值λ所對應的特徵向量。

　　求出特徵值和特徵向量有什麼好處呢？ 咱們能夠將矩陣A特徵分解。若是咱們求出了矩陣A的n個特徵值λ1≤λ2≤...≤λn,以及這n個特徵值所對應的特徵向量{w1,w2,...wn}，那麼矩陣A就能夠用下式的特徵分解表示：

A=WΣW⁻¹

其中W是這 n個特徵向量所張成的n×n維矩陣，而Σ爲這n個特徵值爲主對角線的n×n維矩陣。

　　通常咱們會把W的這n個特徵向量標準化，即知足||wi||2=1, 或者說w^T_iw_i=1，此時W的n個特徵向量爲標準正交基，知足W^TW=I，即W^T=W−1, 也就是說W爲酉矩陣。

　　　　這樣咱們的特徵分解表達式能夠寫成A=WΣW^T

　　　　注意到要進行特徵分解，矩陣A必須爲方陣。那麼若是A不是方陣，即行和列不相同時，咱們還能夠對矩陣進行分解嗎？答案是能夠，此時咱們的SVD登場了。

  SVD的定義

　　　　SVD也是對矩陣進行分解，可是和特徵分解不一樣，SVD並不要求要分解的矩陣爲方陣。假設咱們的矩陣A是一個m×nm×n的矩陣，那麼咱們定義矩陣A的SVD爲：

                                                       A=UΣV^T   

　　　　其中U是一個m×m的矩陣，Σ是一個m×n的矩陣，除了主對角線上的元素之外全爲0，主對角線上的每一個元素都稱爲奇異值，V是一個n×n的矩陣。U和V都是酉矩陣，即知足U^TU=I,V^TV=I。下圖能夠很形象的看出上面SVD的定義：

　　　　那麼咱們如何求出SVD分解後的U,Σ,V這三個矩陣呢？

　　　　若是咱們將A的轉置和A作矩陣乘法，那麼會獲得n×nn×n的一個方陣ATA。既然ATA是方陣，那麼咱們就能夠進行特徵分解，獲得的特徵值和特徵向量知足下式：

                                                                    (A^TA)vi=λivi

 

　　　　這樣咱們就能夠獲得矩陣ATA的n個特徵值和對應的n個特徵向量v了。將ATA的全部特徵向量張成一個n×n的矩陣V，就是咱們SVD公式裏面的V矩陣了。通常咱們將V中的每一個特徵向量叫作A的右奇異向量。

　　　　若是咱們將A和A的轉置作矩陣乘法，那麼會獲得m×m的一個方陣AAT。既然AAT是方陣，那麼咱們就能夠進行特徵分解，獲得的特徵值和特徵向量知足下式：

                                                               (AA^T)ui=λiui

 

　　　　這樣咱們就能夠獲得矩陣AAT的m個特徵值和對應的m個特徵向量u了。將AAT的全部特徵向量張成一個m×m的矩陣U，就是咱們SVD公式裏面的U矩陣了。通常咱們將U中的每一個特徵向量叫作A的左奇異向量。

　　　　U和V咱們都求出來了，如今就剩下奇異值矩陣Σ沒有求出了。因爲Σ除了對角線上是奇異值其餘位置都是0，那咱們只須要求出每一個奇異值σ就能夠了。

　　　　咱們注意到:

                                                              A=UΣV^T⇒AV=UΣV^TV⇒AV=UΣ⇒Av_i=σ_iu_i⇒σ_i=Av_i/u_i

 

 　　　 這樣咱們能夠求出咱們的每一個奇異值，進而求出奇異值矩陣Σ。

             SVD的一些性質　

　　　　上面幾節咱們對SVD的定義和計算作了詳細的描述，彷佛看不出咱們費這麼大的力氣作SVD有什麼好處。那麼SVD有什麼重要的性質值得咱們注意呢？

　　　　對於奇異值,它跟咱們特徵分解中的特徵值相似，在奇異值矩陣中也是按照從大到小排列，並且奇異值的減小特別的快，在不少狀況下，前10%甚至1%的奇異值的和就佔了所有的奇異值之和的99%以上的比例。也就是說，咱們也能夠用最大的k個的奇異值和對應的左右奇異向量來近似描述矩陣。也就是說：

                                                   Am×n=Um×mΣm×nV^Tn×n≈Um×kΣk×kV^Tk×n

　　　　其中k要比n小不少，也就是一個大的矩陣A能夠用三個小的矩陣Um×k,Σk×k,V^Tk×n來表示。以下圖所示，如今咱們的矩陣A只須要灰色的部分的三個小矩陣就能夠近似描述了。

　　　　因爲這個重要的性質，SVD能夠用於PCA降維，來作數據壓縮和去噪。也能夠用於推薦算法，將用戶和喜愛對應的矩陣作特徵分解，進而獲得隱含的用戶需求來作推薦。同時也能夠用於NLP中的算法，好比潛在語義索引（LSI）。

以上轉自：http://www.cnblogs.com/pinard/p/6251584.html

SVD協同過濾：

假設存在如下user和item的數據矩陣：

這是一個極其稀疏的矩陣，這裏把這個評分矩陣記爲R，其中的元素表示user對item的打分，「？」表示未知的，也就是要你去預測的，如今問題來了：如何去預測未知的評分值呢？從上面的SVD的性質:  Am×n=Um×mΣm×nVTn×n≈Um×kΣk×kVTk×n，能夠獲得：

一個m*n的打分矩陣R能夠由分解的兩個小矩陣U（m*k）和V（k*n）的乘積來近似，即 R=UV^T,k<=m,n 

將這種分解方式體現協同過濾中，即有：

           (matrix factorization model，MF模型 )

 在這樣的分解模型中,Pu表明用戶隱因子矩陣（表示用戶u對因子k的喜愛程度),Qi表示電影隱因子矩陣（表示電影i在因子k上的程度）。

 

SVD推薦算法公式以下：

    

這裏須要解釋一下各個參數的含義：

對於電影評分實例，首先獲得訓練數據集 user_id,movie_id和rating，u表示打分矩陣中全部評分值的平均值，bi在這個公式中應該是一個參數值，而不是向量，能夠這樣理解，首先初始化一個表明item的向量bi，向量維度是item的個數，公式中的bi是指bi[movie_id],同理，bu表明bu[user_id],r_{ui_}hat 表示預測的評分值.

加入防止過擬合的 λ 參數，能夠獲得下面的優化函數：

     

利用隨機梯度降低算法更新參數：

 代碼體現：

# -*- coding: utf-8 -*-
"""
Created on Thu Apr 19 16:28:10 2018

@author: 
"""
import numpy as np  
import random  
import os 
  
class SVD:  
    def __init__(self,mat,K=20):  
        self.mat=np.array(mat)  
        self.K=K  
        self.bi={}  
        self.bu={}  
        self.qi={}  
        self.pu={}  
        self.avg=np.mean(self.mat[:,2])  
        for i in range(self.mat.shape[0]):  
            uid=self.mat[i,0]  
            iid=self.mat[i,1]  
            self.bi.setdefault(iid,0)  
            self.bu.setdefault(uid,0)  
            self.qi.setdefault(iid,np.random.random((self.K,1))/10*np.sqrt(self.K))  
            self.pu.setdefault(uid,np.random.random((self.K,1))/10*np.sqrt(self.K))  
    def predict(self,uid,iid):  #預測評分的函數  
        #setdefault的做用是當該用戶或者物品未出現過期，新建它的bi,bu,qi,pu，並設置初始值爲0  
        self.bi.setdefault(iid,0)  
        self.bu.setdefault(uid,0)  
        self.qi.setdefault(iid,np.zeros((self.K,1)))  
        self.pu.setdefault(uid,np.zeros((self.K,1)))  
        rating=self.avg+self.bi[iid]+self.bu[uid]+np.sum(self.qi[iid]*self.pu[uid]) #預測評分公式  
        #因爲評分範圍在1到5，因此當分數大於5或小於1時，返回5,1.  
        if rating>5:  
            rating=5  
        if rating<1:  
            rating=1  
        return rating  
      
    def train(self,steps=30,gamma=0.04,Lambda=0.15):    #訓練函數，step爲迭代次數。  
        print('train data size',self.mat.shape)  
        for step in range(steps):  
            print('step',step+1,'is running')  
            KK=np.random.permutation(self.mat.shape[0]) #隨機梯度降低算法，kk爲對矩陣進行隨機洗牌  
            rmse=0.0;mae=0
            for i in range(self.mat.shape[0]):  
                j=KK[i]  
                uid=self.mat[j,0]  
                iid=self.mat[j,1]  
                rating=self.mat[j,2]  
                eui=rating-self.predict(uid, iid)  
                rmse+=eui**2  
                mae+=abs(eui)
                self.bu[uid]+=gamma*(eui-Lambda*self.bu[uid])    
                self.bi[iid]+=gamma*(eui-Lambda*self.bi[iid])  
                tmp=self.qi[iid]  
                self.qi[iid]+=gamma*(eui*self.pu[uid]-Lambda*self.qi[iid])  
                self.pu[uid]+=gamma*(eui*tmp-Lambda*self.pu[uid])  
            gamma=0.93*gamma  #gamma以0.93的學習率遞減 
            print('rmse is {0:3f}, ase is {1:3f}'.format(np.sqrt(rmse/self.mat.shape[0]),mae/self.mat.shape[0]))  
      
    def test(self,test_data):   
          
        test_data=np.array(test_data)  
        print('test data size',test_data.shape)  
        rmse=0.0;mae=0  
        for i in range(test_data.shape[0]):  
            uid=test_data[i,0]  
            iid=test_data[i,1]  
            rating=test_data[i,2]  
            eui=rating-self.predict(uid, iid)  
            rmse+=eui**2  
            mae+=abs(eui)
        print('rmse is {0:3f}, ase is {1:3f}'.format(np.sqrt(rmse/self.mat.shape[0]),mae/self.mat.shape[0]))  
              
              
def getData(file_name): 
    """
    獲取訓練集和測試集的函數 
    """  
    data=[] 
    with open(os.path.expanduser(file_name)) as f:                          
        for line in f.readlines():  
            list=line.split('::')  
            data.append([int(i) for i in list[:3]])  
    random.shuffle(data)  
    train_data=data[:int(len(data)*7/10)]  
    test_data=data[int(len(data)*7/10):]  
    print('load data finished')  
    print('total data ',len(data))  
    return train_data,test_data  
           
if __name__=='__main__':  
    train_data,test_data=getData('D:/Downloads/ml-1m/ratings.dat')  
    a=SVD(train_data,30)    
    a.train()  
    a.test(test_data)  

測試結果

在訓練集上
rmse is 0.869038, ase is 0.690794
在測試集上
rmse is 0.583027, ase is 0.303116

 

 

 SVD++算法：

SVD算法是指在SVD的基礎上引入隱式反饋，使用用戶的歷史瀏覽數據、用戶歷史評分數據等做爲新的參數。

這裏的N(u)表示用戶u行爲記錄（包括瀏覽的和評過度的商品集合），yj爲隱藏的「評價了電影 j」反映出的我的喜愛偏置。其餘參數同SVD中的參數含義一致。

 利用隨機梯度降低算法更新參數：

代碼體現：

# -*- coding: utf-8 -*-
"""
Created on Sat Apr 21 14:44:25 2018

@author: fanchao3
"""
import numpy as np  
import random
import os

class SVDPP:  
    def __init__(self,mat,K=20):  
        self.mat=np.array(mat)  
        self.K=K  
        self.bi={}  
        self.bu={}  
        self.qi={}  
        self.pu={}  
        self.avg=np.mean(self.mat[:,2])  
        self.y={}  
        self.u_dict={}  
        for i in range(self.mat.shape[0]):  
              
            uid=self.mat[i,0]  
            iid=self.mat[i,1]  
            self.u_dict.setdefault(uid,[])  
            self.u_dict[uid].append(iid)  
            self.bi.setdefault(iid,0)  
            self.bu.setdefault(uid,0)  
            self.qi.setdefault(iid,np.random.random((self.K,1))/10*np.sqrt(self.K))  
            self.pu.setdefault(uid,np.random.random((self.K,1))/10*np.sqrt(self.K))  
            self.y.setdefault(iid,np.zeros((self.K,1))+.1)  
    def predict(self,uid,iid):  #預測評分的函數  
        #setdefault的做用是當該用戶或者物品未出現過期，新建它的bi,bu,qi,pu及用戶評價過的物品u_dict，並設置初始值爲0  
        self.bi.setdefault(iid,0)  
        self.bu.setdefault(uid,0)  
        self.qi.setdefault(iid,np.zeros((self.K,1)))  
        self.pu.setdefault(uid,np.zeros((self.K,1)))  
        self.y.setdefault(uid,np.zeros((self.K,1)))  
        self.u_dict.setdefault(uid,[])  
        u_impl_prf,sqrt_Nu=self.getY(uid, iid)  
        rating=self.avg+self.bi[iid]+self.bu[uid]+np.sum(self.qi[iid]*(self.pu[uid]+u_impl_prf)) #預測評分公式  
        #因爲評分範圍在1到5，因此當分數大於5或小於1時，返回5,1.  
        if rating>5:  
            rating=5  
        if rating<1:  
            rating=1  
        return rating  
      
    #計算sqrt_Nu和∑yj  
    def getY(self,uid,iid):  
        Nu=self.u_dict[uid]  
        I_Nu=len(Nu)  
        sqrt_Nu=np.sqrt(I_Nu)  
        y_u=np.zeros((self.K,1))  
        if I_Nu==0:  
            u_impl_prf=y_u  
        else:  
            for i in Nu:  
                y_u+=self.y[i]  
            u_impl_prf = y_u / sqrt_Nu  
          
        return u_impl_prf,sqrt_Nu  
      
    def train(self,steps=30,gamma=0.04,Lambda=0.15):    #訓練函數，step爲迭代次數。  
        print('train data size',self.mat.shape)  
        for step in range(steps):  
            print('step',step+1,'is running')  
            KK=np.random.permutation(self.mat.shape[0]) #隨機梯度降低算法，kk爲對矩陣進行隨機洗牌  
            rmse=0.0  
            for i in range(self.mat.shape[0]):  
                j=KK[i]  
                uid=self.mat[j,0]  
                iid=self.mat[j,1]  
                rating=self.mat[j,2]  
                predict=self.predict(uid, iid)  
                u_impl_prf,sqrt_Nu=self.getY(uid, iid)  
                eui=rating-predict  
                rmse+=eui**2  
                self.bu[uid]+=gamma*(eui-Lambda*self.bu[uid])    
                self.bi[iid]+=gamma*(eui-Lambda*self.bi[iid])  
                self.pu[uid]+=gamma*(eui*self.qi[iid]-Lambda*self.pu[uid])  
                self.qi[iid]+=gamma*(eui*(self.pu[uid]+u_impl_prf)-Lambda*self.qi[iid])  
                for j in self.u_dict[uid]:  
                    self.y[j]+=gamma*(eui*self.qi[j]/sqrt_Nu-Lambda*self.y[j])  
                                      
            gamma=0.93*gamma  
            print('rmse is',np.sqrt(rmse/self.mat.shape[0]))  
      
    def test(self,test_data):  #gamma以0.93的學習率遞減  
          
        test_data=np.array(test_data)  
        print('test data size',test_data.shape)  
        rmse=0.0  
        for i in range(test_data.shape[0]):  
            uid=test_data[i,0]  
            iid=test_data[i,1]  
            rating=test_data[i,2]  
            eui=rating-self.predict(uid, iid)  
            rmse+=eui**2  
        print('rmse of test data is',np.sqrt(rmse/test_data.shape[0]))  
      
      
def getData(file_name): 
    """
    獲取訓練集和測試集的函數 
    """  
    data=[] 
    with open(os.path.expanduser(file_name)) as f:                          
        for line in f.readlines():  
            list=line.split('::')  
            data.append([int(i) for i in list[:3]])  
    random.shuffle(data)  
    train_data=data[:int(len(data)*7/10)]  
    test_data=data[int(len(data)*7/10):]  
    print('load data finished')  
    print('total data ',len(data))  
    return train_data,test_data  
      
if __name__=='__main__':  
    train_data,test_data=getData('D:/Downloads/ml-1m/ratings.dat')  
    a=SVDPP(train_data,30)    
    a.train()  
    a.test(test_data)  

 

# -*- coding: utf-8 -*-
"""
Created on Thu Apr 19 17:53:34 2018

@author: 
"""

import numpy as np  
import random  
import os 
  
class SVDPP:  
    def __init__(self,mat,K=20):  
        self.mat=np.array(mat)  
        self.K=K 
        self.avg=np.mean(self.mat[:,2]) 
        self.user_num = len(set(self.mat[:,0]))
        self.item_num = len(set(self.mat[:,1]))
        #print("item_num:",self.item_num )
        #user bias
        self.bu = np.zeros(self.user_num, np.double)

        #item bias
        self.bi = np.zeros(self.item_num, np.double)

        #user factor
        self.p = np.zeros((self.user_num, self.K), np.double) + .1

        #item factor
        self.q = np.zeros((self.item_num, self.K), np.double) + .1

        #item preference facotor
        self.y = np.zeros((self.item_num, self.K), np.double) + .1
        self.u_items={}
        for i in range(self.mat.shape[0]):
            uid=self.mat[i,0]
            iid=self.mat[i,1]
            if uid not in self.u_items.keys():                
                self.u_items[uid]=[iid]
            else:
                self.u_items[uid].append(iid)
                
    def train(self,steps=30,gamma=0.04,Lambda=0.15):    #訓練函數，step爲迭代次數。  
        #print('train data size',self.mat.shape) 
        
        for step in range(steps):  
            print('step',step+1,'is running')  
            KK=np.random.permutation(self.mat.shape[0]) #隨機梯度降低算法，kk爲對矩陣進行隨機洗牌  
            rmse=0.0;mae=0
            for i in range(self.mat.shape[0]):  
                j=KK[i]  
                uid=self.mat[j,0]  
                iid=self.mat[j,1]  
                rating=self.mat[j,2]
                Nu=self.u_items[uid]
                I_Nu = len(Nu)
                sqrt_N_u = np.sqrt(I_Nu)
                #基於用戶u點評的item集推測u的implicit偏好
                y_u = np.sum(self.y[Nu], axis=0)
                u_impl_prf = y_u / sqrt_N_u
                #預測值
                rp = self.avg + self.bu[uid] + self.bi[iid] + np.dot(self.q[iid], self.p[uid] + u_impl_prf)
                eui=rating- rp
                rmse+=eui**2  
                mae+=abs(eui)
                #sgd
                self.bu[uid] += gamma * (eui - Lambda * self.bu[uid])
                self.bi[iid] += gamma * (eui - Lambda * self.bi[iid])
                self.p[uid] += gamma * (eui * self.q[iid] - Lambda * self.p[uid])
                self.q[iid] += gamma * (eui * (self.p[uid] + u_impl_prf) - Lambda * self.q[iid])
                for j in Nu:
                    self.y[j] += gamma * (eui * self.q[j] / sqrt_N_u - Lambda * self.y[j])

            gamma=0.93*gamma  #gamma以0.93的學習率遞減 
            print('rmse is {0:3f}, ase is {1:3f}'.format(np.sqrt(rmse/self.mat.shape[0]),mae/self.mat.shape[0]))  
      
    def test(self,test_data):   
          
        test_data=np.array(test_data)  
        print('test data size',test_data.shape)  
        rmse=0.0;mae=0  
        for i in range(test_data.shape[0]):  
            uid=test_data[i,0]  
            iid=test_data[i,1]  
            Nu=self.u_items[uid]
            I_Nu = len(Nu)
            sqrt_N_u = np.sqrt(I_Nu)
            y_u = np.sum(self.y[Nu], axis=0) / sqrt_N_u
            est = self.avg + self.bu[uid] + self.bi[iid] + np.dot(self.q[iid], self.p[uid] + y_u)
            rating=test_data[i,2]  
            eui=rating-est  
            rmse+=eui**2  
            mae+=abs(eui)
        print('rmse is {0:3f}, ase is {1:3f}'.format(np.sqrt(rmse/self.mat.shape[0]),mae/self.mat.shape[0]))  
              
              
def getData(file_name): 
    """
    獲取訓練集和測試集的函數 
    """  
    data=[] 
    with open(os.path.expanduser(file_name)) as f:                          
        for line in f.readlines():
            List=line.split('::')  
            data.append([int(i) for i in List[:3]])  

    random.shuffle(data)  
    train_data=data[:int(len(data)*7/10)]  
    test_data=data[int(len(data)*7/10):]  
    new_train_data=mapping(train_data)
    new_test_data=mapping(test_data)
    print('load data finished')  
    return new_train_data,new_test_data  
def mapping(data):
    """
    將原始的uid，iid映射爲從0開始的編號
    """
    data=np.array(data)
    users=list(set(data[:,0]))
    u_dict={}
    for i in range(len(users)):
        u_dict[users[i]]=i
    items=list(set(data[:,1]))
    i_dict={}
    for j in range(len(items)):
        i_dict[items[j]]=j
    new_data=[]
    for l in data:
        uid=u_dict[l[0]]
        iid=i_dict[l[1]]
        r=l[2]
        new_data.append([uid,iid,r])
    return new_data
     
if __name__=='__main__':  
    train_data,test_data=getData('D:/Downloads/ml-1m/ratings.dat')  
    a=SVDPP(train_data,30)    
    a.train()  
    a.test(test_data)