Patchouli的機器學習系列教程二：目標函數(Objective Functions)——『器』

01-30

之前的內容在我的個人專欄Pure Data入門教程里，日後『道』篇將發表在我的個人專欄里，『器』篇發到景略集智上去。

在形而上的部分里Patchouli的機器學習系列教程二：目標函數(Objective Functions)——『道』，我們介紹了目標函數，介紹了線性回歸，優化以及梯度下降演算法，在『器』篇我將用前面介紹的方法設計一個電影推薦系統。

數據源

import pylab as pltnimport numpy as npnimport pandas as pdnimport podsnmovies = pd.read_csv(movies.csv,encoding=latin-1).set_index(index)nuser_names = list(set(movies.columns)-set(movies.columns[:9]))nY = pd.melt(movies.reset_index(), id_vars=[Film, index], var_name=user, value_name=rating, value_vars=user_names)nY = Y.dropna(axis=0)n

推薦系統，就是曾經被人民日報連續批評過的系統，這個系統從社會學角度上是壓制人性中積極的一面，反對創新的；他把大量用戶喜歡的，偏愛的東西推薦給用戶，使得用戶的思想越來越極端，越來越偏激。媒體，本該是一條人們和外界交流的渠道，然而這種演算法的不當應用，卻成為了極端思潮的溫床。而我們的電影推薦系統同樣有著這樣的社會學缺陷，假設一個人看了《死神來了》，根據系統的推薦，他必定會連續看到死神來了全系列六部作品以及其他的宣揚宿命論的恐怖作品；假設他沒有其他的信息來源，那麼這個人很可能陷入可悲的基於宿命論的惴惴不安當中，從而對生活喪失一切希望，失去主觀能動性，不能成為新時代特色社會主義要求的好公民。因此，屏幕前的各位人工智慧科學家們，要小心自己的力量，你們的能量早已超乎你們的想像。人類絕不能只得到他們想要的東西，有些東西他們不想要，卻必須有，比如教育，比如審美，比如特色社會主義核心價值觀。

知識就是力量，但必須用於正途。

我想你們應該已經明白我的意思了，話不多說，正式進入本次課程。

拉普拉斯在其作品《概率理論的哲學原理》("Essai philosophique sur les probabilités")中聲稱一切存在都可以被數學抽象化，從而假想了拉普拉斯妖的存在；儘管拉普拉斯的理念在日後被從邏輯上和實踐上推翻，但是這種思想在數據科學和機器學習中確是十分有效的。電影，理念，觀點，興趣，無論哪個都無法被計算機所理解，我們面前的電腦，除了數學一無所知；因此我們必須想辦法把這些人類學的內容形式化，抽象化，從數學的角度來解釋和規範這些人類學的存在。一旦我們對於這些存在和概念的數學解釋成功了，那麼計算機也就能理解我們的想法，從而幫助我們進行數據分析了。我們假設幾十萬部電影被放到了一個平面圖書館當中，這個圖書館的正中間是最沒有特色的電影，我都想像不出來能有什麼電影這麼沒特色，總之有這麼個原點。從這個原點出發，往東走，適合情侶們，架子上的電影將變得越來越浪漫，越來越言情，如果我們能走到頭的話，將看到人類歷史上最浪漫的電影；單身狗們便會選擇往西走，越往西走，電影就越和愛情，和浪漫無關，走到頭就會是一部極端冷血的電影。往北走是歷史向，往南走是未來向。這樣一來便簡單明了了，排除電影的其他屬性不談，假設我要找一部極端未來言情向的作品，那麼我便要向著東南方前進，越往這個方向走，我找到的電影越符合我的要求。在我們為電影分好類之後，用戶也可以被這樣分類，通過他們對於電影的評分，計算出這個用戶在電影圖書館中的位置，然後把離他最近的電影推薦給他就好啦。接下來，我們要對用戶和電影做抽象化，由於我們的坐標系是二維的，所以每個對象的位置用一個二維向量就可以表示。值得一提的是，這個向量就是傳說中的特徵向量，我們日後的研究的對象會非常複雜，僅用一個二維向量是展現不了的，要用好長一串向量來表徵這個對象的特性，是為特徵向量。

1.向量的內積

剛剛我們為每個用戶和電影定義了他們在坐標系中的位置，這個位置將被如下表示：

$mathbf{v}_j = begin{bmatrix} v_{j,1} v_{j,2} end{bmatrix} mathbf{u}_i = begin{bmatrix} u_{i,1} u_{i,2} end{bmatrix}$

兩個向量 $mathbf{a}$ 和 $mathbf{b}$ 的內積一般寫作 $mathbf{a}cdotmathbf{b}$ ，在涉及到矩陣的乘法時記作 $mathbf{a}^top mathbf{b}$ 。矩陣的乘法記不清的話就記住口訣：橫乘豎等於數，豎乘橫等於陣；內積是個數，向量是個列。所以內積是把第一個向量轉置乘第二個向量，這在以後的計算中至關重要，千萬不要搞混。而內積的計算有如下公式：

$mathbf{a}cdotmathbf{b}= mathbf{a}^top mathbf{b} = sum_i a_i b_i= |mathbf{a}| |mathbf{b}| costheta$

其中 $|mathbf{x}|$ 是向量的長度， $theta$ 是向量之間的夾角。這個公式可以用來計算向量之間的相似度，放回到我們的電影推薦系統，就是電影和用戶之間的相性。用戶和電影的相性越高，用戶給電影的評分就越高，正好符合了我們擁有的數據集——用戶為電影的評分。評分和內積是正相關的，這對我們來說是個絕好的消息。

根據我們在最一開始導入的數據集，我們使用melt函數將一個用戶和電影的評價矩陣拆分成了電影用戶和評分的三列數據集，並排除了空評分的數據，使得數據更容易被操作。現在我們需要構建兩個包含用戶和坐標的大向量，也就是特徵向量集：

$mathbf{U} = begin{bmatrix} mathbf{u}_1 dots mathbf{u}_n end{bmatrix}^top$ 和 $mathbf{V} = begin{bmatrix} mathbf{v}_1 dots mathbf{v}_m end{bmatrix}^top$

向量中的每個元素都是一個二維向量，代表用戶和電影在坐標系中的位置。

據此我們可以構建出目標函數：

$E(mathbf{U},mathbf{V}) = sum_{i,j} s_{i,j} (y_{i,j} - mathbf{u}_i^top mathbf{v}_j)^2$

其中 $s_{i,j}$ 是一個布爾變數，1則是用戶為此電影打過分，0則是用戶未給此電影打過分。

2.目標函數

與『道』篇中，這個函數只有兩個自變數，我們對其中一個自變數求偏導可以得到：

$frac{partial E(mathbf{U},mathbf{V})}{partial u_{k,ell}} = -2 sum_j s_{k,j} v_{j,ell} (y_{k,j} - mathbf{u}_k^top mathbf{v}_j)$

用這個函數來不斷優化：

$u_{k,ell} leftarrow u_{k,ell} - eta frac{partial E(mathbf{U},mathbf{V})}{partial u_{k,ell}}$

對上面這段推理不熟悉的就回去補習高數吧，以下過程都是該一眼看出來的，如果沒有的話回去做一萬道求偏導數的題再來看這個教程：

$begin{align} frac{partial E(mathbf{U},mathbf{V})}{partial u_{k,ell}} &= frac{partial}{partial u_{k,ell}} sum_j s_{k,j}(y_{k,j} - mathbf{u}_k^top mathbf{v}_j)^2 &= sum_j s_{k,j} left( frac{partial}{partial u_{k,ell}} (y_{k,j} - mathbf{u}_k^top mathbf{v}_j)^2 right) &= sum_j s_{k,j} left( 2(y_{k,j} - mathbf{u}_k^top mathbf{v}_j) frac{partial}{partial u_{k,ell}} (y_{k,j} - mathbf{u}_k^top mathbf{v}_j) right) &= -2sum_j s_{k,j} (y_{k,j} - mathbf{u}_k^top mathbf{v}_j) left( frac{partial}{partial u_{k,ell}} (mathbf{u}_k^top mathbf{v}_j) right) end{align} Rightarrow frac{partial E(mathbf{U},mathbf{V})}{partial u_{k,ell}} = -2 sum_j s_{k,j} v_{j,ell} (y_{k,j} - mathbf{u}_k^top mathbf{v}_j)$

與之相對的， $v_{k,ell}$ 的偏導數如下：

$frac{partial E(mathbf{U},mathbf{V})}{partial v_{k,ell}} = -2sum_i s_{i,k} u_{i,ell} (y_{i,k} - mathbf{u}_i^top mathbf{v}_k)$

就像『道』篇里所講的那樣，我們不知道電影和用戶在坐標系中的位置，那麼就像我們當初設定 $m$ 和 $c$ 那樣隨機一個出來。

q = 2nlearn_rate = 0.01nU = pd.DataFrame(np.random.normal(size=(len(user_names), q))*0.001, index=user_names)nV = pd.DataFrame(np.random.normal(size=(len(movies.index), q))*0.001, index=movies.index)nY[rating] -= Y[rating].mean()n

U和V便是兩個隨機的數據集，為每個用戶和電影隨機了一個位置出來，讓他們不斷優化以尋找自己在坐標系中的位置。由於我們的坐標系是以零為原點，所以將整體的評分減去平均值，使得觀眾對電影的評分均勻地分布在坐標系零點的兩側。我們繼續設計目標函數：

import sysndef objective_gradient(Y, U, V):n gU = pd.DataFrame(np.zeros((U.shape)), index=U.index)n gV = pd.DataFrame(np.zeros((V.shape)), index=V.index)n obj = 0.n for ind, series in Y.iterrows():n film = series[index]n user = series[user]n rating = series[rating]n prediction = np.dot(U.loc[user], V.loc[film])n diff = prediction - rating n obj += diff*diffn gU.loc[user] += 2*diff*V.loc[film]n gV.loc[film] += 2*diff*U.loc[user]n return obj, gU, gVnlearn_rate = 0.01niterations = 100nfor i in range(iterations):n obj, gU, gV = objective_gradient(Y, U, V)n print("迭代次數為", i+1, "偏差為", obj)n U -= learn_rate*gUn V -= learn_rate*gVn

運行時我們就發現了，哪怕就這幾十個用戶，幾百個數據，這樣巨大的優化次數對計算機來說也是吃力的，在我的電腦要跑一兩分鐘才好，當我們面對數以十萬計的數據時，梯度下降演算法必須要進行改進，我們不能每次都把U和V兩個坐標集向量全迭代一遍，而換做隨機取點，對每個點進行迭代，因此我們要使用隨機梯度下降法來處理這份數據：

def stochastic_descent(obj, Y, U, V, learn_rate):n Y = Y.iloc[np.random.permutation(len(Y))] n obj = 0.n for ind, series in Y.iterrows():n film = series[index]n user = series[user]n rating = series[rating]n prediction = np.dot(U.loc[user], V.loc[film]) n diff = prediction - rating n obj += diff**2n U.loc[user] -= 2*learn_rate*diff*V.loc[film]n V.loc[film] -= 2*learn_rate*diff*U.loc[user]n return obj, U, VnU = pd.DataFrame(np.random.normal(size=(len(user_names), q))*0.001, index=user_names)nV = pd.DataFrame(np.random.normal(size=(len(movies.index), q))*0.001, index=movies.index)nlearn_rate = 0.05niterations = 100nobj = 0.nfor i in range(iterations):n obj, U, V = stochastic_descent(obj, Y, U, V, learn_rate)n print("update {} objective function: {}".format(i+1, obj))nfig, ax = plt.subplots()nplt.plot(U[0], U[1], rx, label = users)nplt.plot(V[0], V[1], go, label = movies)nfor name in U.index:n ax.text(U[0][name], U[1][name], name)nlegend = ax.legend(loc=upper left)n

結果：

前面說的70萬數據集我回去測了又測，覺得還是會拖垮集智的網站，於是用to_csv函數閹了份10萬版的給大家測試演算法，接下來我們就要用隨機梯度下降演算法來處理這份數據：

100K

def sgd(obj, Y, U, V, learn_rate):n Y = Y.iloc[np.random.permutation(len(Y))] n obj = 0.n for ind, series in Y.iterrows():n item = series[item]n user = series[user]n rating = series[rating]n prediction = np.dot(U.loc[user], V.loc[item]) n diff = prediction - rating n obj += diff**2n U.loc[user] -= 2*learn_rate*diff*V.loc[item]n V.loc[item] -= 2*learn_rate*diff*U.loc[user]n return obj, U, VnY = pd.read_csv(100k.csv)nY = Y.reset_index()nY[rating] -= Y[rating].mean()nq = 2nU = pd.DataFrame(np.random.normal(size=(len(Y.user.unique()), q))*0.001, index = Y.user.unique())nV = pd.DataFrame(np.random.normal(size=(len(Y.item.unique()), q))*0.001, index = Y.item.unique())nplt.figure(figsize=(11, 7))nplt.scatter(V[0], V[1], c=r, marker=x)nplt.show()nlearn_rate = 0.008niterations = 2nobj = 0.nfor i in range(iterations):n t0 = time.time()n obj, U, V = sgd(obj, Y, U, V, learn_rate)n t1 = time.time()n print("Update {}, objective function = {}, computation time = {} minutes".format(i, obj, (t1 - t0)/60.))n n plt.figure(figsize=(11, 7))n plt.scatter(V[0], V[1], c=r, marker=x)n plt.show()n

迭代前（隨機點）

第一次迭代後的結果

第二次迭代後的結果

以上便是隨機梯度演算法在電影推薦系統中的應用。