眾所周知，協同過濾（Collaboration Filtering）演算法是推薦系統中最常用的一種演算法。今天我們就以電影推薦為例，簡要論述基本原理，最終給出實現的python程式碼。

1. 問題定義

       假設現有一個二維表，記錄著每個使用者對所看電影的評分情況，如下圖所示：

       從圖中可知，該二維表記錄著4個使用者關於5部電影的評分情況，其中部分資料缺失。那麼我們能否根據表格已有資訊來對缺失值進行填充，實際意義就是能否根據使用者喜好和電影特徵來評估該使用者對該電影的評分情況。因此，該問題需要解決兩個小問題。一是使用者喜好如何來衡量？二是電影特徵如何確定？下面，我們將分別從這兩個方面進行解釋。

2. 已知電影特徵

       假設已知每部電影具有兩個特徵x1、x2，其中x1表示該電影浪漫的程度，x2表示該電影搞笑的程度。如果一部電影的特徵表示為（1,0），那麼就認為該電影是一部愛情浪漫片；如果特徵表示為（0,1），那麼就認為該電影是一部喜劇片；如果表示為（0.6,0.5），那麼就認為該電影是一部愛情浪漫喜劇片。

       如果現在已經通過手工標記，計算出所有電影的特徵資料，那麼我們就可以利用這些資料來學習每個使用者的喜好了。由於電影特徵x是二維的，使用者喜好θ也應該是二維的。如果學習到Alice的使用者喜好θ為（4.5,1），而《cute puppies of love》這部電影的特徵為（0.9,0），那麼就可以評估Alice對《cute puppies of love》的評分為4.05。該評分在現實中是合理的，Alice的使用者喜好特徵表明她喜歡愛情片，《cute puppies of love》的特徵表明它是一部愛情片，因此評分也就較高。下面將給出數學化描述。

      對於每一個使用者而言，其優化目標為：

      其中θj表示第j個使用者的愛好特徵，xi表示第i部電影的特徵，y(i,j)表示第j個使用者對第i部電影的評分，i:r(i,j)=1表示第j個使用者對第i部電影是評了分的（不是缺失值），mj表示第j個使用者評分電影的數目。右側那一項為正則項，是防止過擬合的（之前演算法經常見到，這裡不再贅述）。由於左右兩項具有mj,因此上式還可以寫成：

接著，利用梯度下降更新θj

最終所求θj即為第j個使用者的喜好特徵。

3. 已知使用者喜好

        如果是已知使用者喜好θ，那麼我們就可以學習電影的特徵x了。對於每一部電影，其優化函式為

接著，利用梯度下降更新xi

最終所求xi即為第i部電影的特徵。

4. 兩者結合

        如果現在電影特徵和使用者喜好特徵都未知，我們可以把上面兩個目標函式結合起來進行優化。如下圖所示

因此，可以得出如下協同過濾演算法：

最終，我們可以得到使用者喜好和電影特徵。根據訓練得到的這些引數，我們就可以預測那些缺失值了，然後就可以把那些預測的高分電影推薦給使用者了。

5. 程式碼實現

import numpy as np
import pandas as pd

# 定義初始化引數函式
# 輸入：特徵數量、使用者數量、產品數量
# 輸出：使用者特徵初始矩陣、產品特徵初始矩陣
def Initialize_Parameters(num_features, num_users, num_products):
    user_matrix = np.random.rand(num_users, num_features)
    product_matrix = np.random.rand(num_products, num_features)
    return user_matrix, product_matrix

# 計算當前的代價函式
# 輸入：當前的使用者矩陣、產品矩陣、缺失的二維表、懲罰因子lambdaa
# 輸出：當前代價
def get_cost(ori_data, user_matrix, product_matrix, lambdaa):
    nan_index = np.isnan(ori_data) # 記錄二維表中缺失的索引
    ori_data[nan_index] = 0 # 把缺失值填充為0
    predict_data = np.dot(user_matrix, product_matrix.T) # 計算預測的評分
    temp = predict_data - ori_data # 計算兩矩陣的差值
    temp[nan_index] = 0 # 缺失值不算在代價裡
    cost = 0.5*np.sum(temp**2) + 0.5*lambdaa*(np.sum(user_matrix**2)+np.sum(product_matrix**2))# 計算平方
    ori_data[nan_index] = np.nan # 恢復原資料
    return cost

# 對使用者特徵進行偏導
# 輸入：當前的使用者矩陣、產品矩陣、缺失的二維表、懲罰因子lambdaa
# 輸出：使用者特徵偏導矩陣
def get_user_derivatives(ori_data, user_matrix, product_matrix, lambdaa=1):
    nan_index = np.isnan(ori_data) # 記錄二維表中缺失的索引
    ori_data[nan_index] = 0 # 把缺失值填充為0
    predict_data = np.dot(user_matrix, product_matrix.T) # 計算預測的評分
    temp = predict_data - ori_data # 計算兩矩陣的差值
    temp[nan_index] = 0 # 缺失值不算在代價裡
    ori_data[nan_index] = np.nan # 恢復原資料
    
    num_user = user_matrix.shape[0] # 計算使用者數目
    feature_user = user_matrix.shape[1] # 計算特徵數目
    user_dervatives = np.random.rand(num_user, feature_user) # 宣告使用者特徵偏導數矩陣
    
    for i in range(num_user):
        for j in range(feature_user):
            user_dervatives[i][j] = np.dot(temp[i], product_matrix[:,j]) + lambdaa * user_matrix[i][j]
    return user_dervatives

# 對產品特徵進行偏導
# 輸入：當前的使用者矩陣、產品矩陣、缺失的二維表、懲罰因子lambdaa
# 輸出：產品特徵偏導矩陣
def get_product_derivatives(ori_data, user_matrix, product_matrix, lambdaa=1):
    nan_index = np.isnan(ori_data) # 記錄二維表中缺失的索引
    ori_data[nan_index] = 0 # 把缺失值填充為0
    predict_data = np.dot(user_matrix, product_matrix.T) # 計算預測的評分
    temp = predict_data - ori_data # 計算兩矩陣的差值
    temp[nan_index] = 0 # 缺失值不算在代價裡
    ori_data[nan_index] = np.nan # 恢復原資料
    
    num_product = product_matrix.shape[0] # 計算產品數目
    feature_product = product_matrix.shape[1] # 計算特徵數目
    product_dervatives = np.random.rand(num_product, feature_product) # 宣告產品特徵偏導數矩陣
    
    for i in range(num_product):
        for j in range(feature_product):
            product_dervatives[i][j] = np.dot(temp[:,i], user_matrix[:,j]) + lambdaa * product_matrix[i][j]
    return product_dervatives
    

# 根據含有缺失值的二維表，學習相關引數
# 輸入：含有缺失值的二維表、使用者特徵初始矩陣、產品特徵初始矩陣、迭代次數、學習效率learning_rate、懲罰因子lambdaa
# 輸出：最優使用者特徵矩陣、最優產品特徵矩陣
def CF(ori_data, user_matrix, product_matrix, iterate_num=500, learning_rate=0.01, lambdaa=1):
    for i in range(iterate_num):
        cost = get_cost(ori_data, user_matrix, product_matrix, lambdaa) # 計算當前代價
        user_derivatives = get_user_derivatives(ori_data, user_matrix, product_matrix, lambdaa) # 對使用者特徵求偏導
        product_derivates = get_product_derivatives(ori_data, user_matrix, product_matrix, lambdaa) # 對產品特徵求偏導
        user_matrix = user_matrix - learning_rate * user_derivatives # 更新引數
        product_matrix = product_matrix - learning_rate * product_derivates
        print i,'th cost:', cost
        
    return user_matrix, product_matrix


# 根據學習的引數，進行評估 
# 輸入：使用者特徵矩陣、產品特徵矩陣
# 輸出：不含缺失值的二維表
def Evaluate_Score(user_matrix, product_matrix):
    return np.dot(user_matrix, product_matrix.T)



# 主函式
if __name__ == '__main__':
    ori_data = np.array([[5,5,np.nan,0,0],[5,np.nan,4,0,0],[0,np.nan,0,5,5],[0,0,np.nan,4,np.nan]])
    #user_matrix = np.array([[5,0.1],[5,0.1],[0.1,5],[0.1,5]])
    #product_matrix = np.array([[0.9,0.1],[1.0,0.01],[0.99,0.01],[0.1,1.0],[0.1,0.9]])         
    user_matrix, product_matrix = Initialize_Parameters(2,ori_data.shape[0],ori_data.shape[1])
    user_matrix, product_matrix = CF(ori_data, user_matrix, product_matrix, iterate_num=100, learning_rate=0.1, lambdaa=0)
    score = Evaluate_Score(user_matrix, product_matrix)
    print score

6. 實驗結果

      以問題定義中的例子作為測試，經過100次迭代，最終預測的資料如下：

[[  4.99997719e+00   4.99998362e+00   3.99998939e+00   2.70274725e-03    2.70273945e-03]
 [  4.99997832e+00   4.99998474e+00   3.99999029e+00   2.70274079e-03    2.70273299e-03]
 [ -1.55079571e-07  -9.98770235e-08  -7.57981453e-08   5.00001356e+00    5.00000909e+00]
 [ -1.18565589e-07  -7.44035873e-08  -5.62400915e-08   4.00000695e+00    4.00000337e+00]]

由此可見，預測的資料和原始資料很接近，因此該演算法很有效。在測試的過程中發現一個問題：人為初始化引數有時無法進行梯度下降，我認為由於該演算法的特殊性致使在前幾次迭代中代價函式突然變大，從而無法繼續迭代（有更好理解的歡迎評論交流）。在實驗中，直接使用隨機初始化引數即可。

附改進的程式碼（改用RMSprop梯度下降，原理）

import numpy as np
import pandas as pd

# 定義初始化引數函式
# 輸入：特徵數量、使用者數量、產品數量
# 輸出：使用者特徵初始矩陣、產品特徵初始矩陣
def Initialize_Parameters(num_features, num_users, num_products):
    user_matrix = np.random.rand(num_users, num_features)
    product_matrix = np.random.rand(num_products, num_features)
    return user_matrix, product_matrix

# 計算當前的代價函式
# 輸入：當前的使用者矩陣、產品矩陣、缺失的二維表、懲罰因子lambdaa
# 輸出：當前代價
def get_cost(ori_data, user_matrix, product_matrix, lambdaa):
    nan_index = np.isnan(ori_data) # 記錄二維表中缺失的索引
    ori_data[nan_index] = 0 # 把缺失值填充為0
    predict_data = np.dot(user_matrix, product_matrix.T) # 計算預測的評分
    temp = predict_data - ori_data # 計算兩矩陣的差值
    temp[nan_index] = 0 # 缺失值不算在代價裡
    cost = 0.5*np.sum(temp**2) + 0.5*lambdaa*(np.sum(user_matrix**2)+np.sum(product_matrix**2))# 計算平方
    ori_data[nan_index] = np.nan # 恢復原資料
    
    return cost

# 對使用者特徵進行偏導
# 輸入：當前的使用者矩陣、產品矩陣、缺失的二維表、懲罰因子lambdaa、加權平均引數
# 輸出：使用者特徵偏導矩陣、加權平均矩陣
def get_user_derivatives(ori_data, user_matrix, product_matrix, weight_average_matrix, lambdaa=1, weight_average_para = 0):
    
    nan_index = np.isnan(ori_data) # 記錄二維表中缺失的索引
    ori_data[nan_index] = 0 # 把缺失值填充為0
    predict_data = np.dot(user_matrix, product_matrix.T) # 計算預測的評分
    temp = predict_data - ori_data # 計算兩矩陣的差值
    temp[nan_index] = 0 # 缺失值不算在代價裡
    ori_data[nan_index] = np.nan # 恢復原資料
    
    num_user = user_matrix.shape[0] # 計算使用者數目
    feature_user = user_matrix.shape[1] # 計算特徵數目
    user_dervatives = np.random.rand(num_user, feature_user) # 宣告使用者特徵偏導數矩陣
    
    for i in range(num_user):
        for j in range(feature_user):
            user_dervatives[i][j] = np.dot(temp[i], product_matrix[:,j]) + lambdaa * user_matrix[i][j]
    
    weight_average_matrix = weight_average_para * weight_average_matrix + (1 - weight_average_para) * (user_dervatives ** 2) # 計算加權平均
    user_dervatives = user_dervatives / (weight_average_matrix**0.5)  # 計算變換的偏導
    return user_dervatives, weight_average_matrix 

# 對產品特徵進行偏導
# 輸入：當前的使用者矩陣、產品矩陣、缺失的二維表、懲罰因子lambdaa
# 輸出：產品特徵偏導矩陣
def get_product_derivatives(ori_data, user_matrix, product_matrix, weight_average_matrix, lambdaa=1, weight_average_para=0):
    nan_index = np.isnan(ori_data) # 記錄二維表中缺失的索引
    ori_data[nan_index] = 0 # 把缺失值填充為0
    predict_data = np.dot(user_matrix, product_matrix.T) # 計算預測的評分
    temp = predict_data - ori_data # 計算兩矩陣的差值
    temp[nan_index] = 0 # 缺失值不算在代價裡
    ori_data[nan_index] = np.nan # 恢復原資料
    
    num_product = product_matrix.shape[0] # 計算產品數目
    feature_product = product_matrix.shape[1] # 計算特徵數目
    product_dervatives = np.random.rand(num_product, feature_product) # 宣告產品特徵偏導數矩陣
    
    for i in range(num_product):
        for j in range(feature_product):
            product_dervatives[i][j] = np.dot(temp[:,i], user_matrix[:,j]) + lambdaa * product_matrix[i][j]
    
    weight_average_matrix = weight_average_para * weight_average_matrix + (1 - weight_average_para) * (product_dervatives ** 2) # 計算加權平均
    product_dervatives = product_dervatives / (weight_average_matrix**0.5)  # 計算變換的偏導

    return product_dervatives, weight_average_matrix 
    

# 根據含有缺失值的二維表，學習相關引數
# 輸入：含有缺失值的二維表、使用者特徵初始矩陣、產品特徵初始矩陣、迭代次數、學習效率learning_rate、懲罰因子lambdaa
# 輸出：最優使用者特徵矩陣、最優產品特徵矩陣
def CF(ori_data, user_matrix, product_matrix, iterate_num=500, learning_rate=0.01, lambdaa=1, weight_average_para=0.5):
    user_weight_average_matrix = np.zeros(user_matrix.shape) # 初始化使用者偏導加權平均為0
    product_weight_average_matrix = np.zeros(product_matrix.shape) # 初始化產品偏導加權平均為0
    for i in range(iterate_num):
        cost = get_cost(ori_data, user_matrix, product_matrix, lambdaa) # 計算當前代價
        user_derivatives, user_weight_average_matrix = get_user_derivatives(ori_data, user_matrix, product_matrix, user_weight_average_matrix, lambdaa, weight_average_para) # 對使用者特徵求偏導
        product_derivates, product_weight_average_matrix = get_product_derivatives(ori_data, user_matrix, product_matrix, product_weight_average_matrix, lambdaa, weight_average_para) # 對產品特徵求偏導
        user_matrix = user_matrix - learning_rate * user_derivatives # 更新引數
        product_matrix = product_matrix - learning_rate * product_derivates
        print i,'th cost:', cost
        
    return user_matrix, product_matrix


# 根據學習的引數，進行評估 
# 輸入：使用者特徵矩陣、產品特徵矩陣
# 輸出：不含缺失值的二維表
def Evaluate_Score(user_matrix, product_matrix):
    return np.dot(user_matrix, product_matrix.T)



# 主函式
if __name__ == '__main__':
    ori_data = pd.read_csv('cf_data.csv')
    columns = ori_data.columns
    ori_data = np.array(ori_data)
    #ori_data = np.array([[5,5,np.nan,0,0],[5,np.nan,4,0,0],[0,np.nan,0,5,5],[0,0,np.nan,4,np.nan]])
    
    user_matrix, product_matrix = Initialize_Parameters(20,ori_data.shape[0],ori_data.shape[1])
    user_matrix, product_matrix = CF(ori_data, user_matrix, product_matrix, iterate_num=100, learning_rate=0.01, lambdaa=0)
    score = Evaluate_Score(user_matrix, product_matrix)
    predict_cf_data = pd.DataFrame(score, columns=columns)
    predict_cf_data.to_csv('predict_cf_data.csv', index=False)