1 主要思路回顧

具體可見：論文筆記 Bayesian Probabilistic Matrix Factorizationusing Markov Chain Monte Carlo （ICML 2008）_UQI-LIUWJ的博客-CSDN博客

2? 導入庫

import numpy as np from numpy.linalg import inv as inv from scipy.linalg import khatri_rao as kr_prod from scipy.stats import wishart import matplotlib.pyplot as plt %matplotlib inline

3 數據預處理

import scipy.iotensor = scipy.io.loadmat('tensor.mat') tensor = tensor['tensor'] ##源數據集部分，214個路段,61 天，144個10分鐘 (214,61,144)random_tensor = scipy.io.loadmat('random_tensor.mat') random_tensor = random_tensor['random_tensor'] #(214, 61, 144)維度的[0,1]隨機矩陣dense_mat = tensor.reshape([tensor.shape[0], tensor.shape[1] * tensor.shape[2]]) #(214, 8784)missing_rate = 0.2binary_mat = (np.round(random_tensor + 0.5 - missing_rate).reshape([random_tensor.shape[0], random_tensor.shape[1] * random_tensor.shape[2]])) #(214, 8784)維的0、1矩陣，表示哪些元素保留哪些元素舍棄 #加0.5 表示四舍五入，-missin_rate表示多少的被舍棄sparse_mat = np.multiply(dense_mat, binary_mat) #(214, 8784) #保留有的部分

4 采樣U

分為兩步：

1）采樣θU

2） 采樣U

?

?

def sample_factor_U(alpha_sparse_mat, alpha_ind, U, V, alpha, beta0 = 1):"""Sampling N-by-R factor matrix W and its hyperparameters (mu_w, Lambda_w)."""'''alpha_sparse_mat, 觀測矩陣*alphaalpha_ind, 示性矩陣*alphaU, 路段特征矩陣 [dim1,rank]V, 時間特征矩陣 [dim2,rank]alpha, #WX和Y之間的精度（協方差矩陣的倒數）beta0 = 1, '''dim1, rank = U.shapeU_bar = np.mean(U, axis = 0)#每一列元素取一個均值，一共 rank個元素，U_bar是U矩陣各個行向量的均值向量temp = dim1 / (dim1 + beta0)#N/N+β0var_mu_hyper = temp * U_bar#N U_bar/N+β0，由于μ0為0，所以這個也可以看成μ0*var_U_hyper = inv(np.eye(rank) + cov_mat(U, U_bar) + temp * beta0 * np.outer(U_bar, U_bar))#np.eye(rank) 就是I，即我們的W0# cov_mat(U, U_bar) 就是N S_bar#np.outer(U_bar, U_bar) 相當于 U_bar.reshape(dim1,1) @ U_bar.reshape(1,dim)#temp * beta0 * np.outer(U_bar, U_bar)) 即求W0*的最后一項#所以這個整體可以看成W0*var_Lambda_hyper = wishart.rvs(df = dim1 + rank, scale = var_U_hyper)#df就是v0*#這一部分是從威沙特分布里面取一個[rank,rank]的矩陣樣本#也就是 Lambda_Uvar_mu_hyper = np.random.multivariate_normal(var_mu_hyper, np.linalg.inv((dim1 + beta0) * var_Lambda_hyper))#從N(μ0*，(β0* Lambda_U)……-1)中選取一個μU向量#[rank,1]################################################ ################################################ ## 上面的部分是在采樣ΘU，下面的部分在采樣U ## ################################################ ################################################var1 = V.T#[rank,dim2]var2 = kr_prod(var1, var1)#[ranl*rank,dim2]#alpha_ind [dim1,dim2]var3 = (var2 @ alpha_ind.T).reshape([rank, rank, dim1])+ var_Lambda_hyper[:, :, np.newaxis]#第一項相當于有值的部分相乘#var_Lambda_hyper——[rank,rank],這里給他強行添加了一個維度#第一項是12式的右邊，第二項是12式的左邊var4 = var1 @ alpha_sparse_mat.T + (var_Lambda_hyper @ var_mu_hyper)[:, np.newaxis]#第一項是13式的第一項，第二項是13式的第二項 【不含最前面的-1】for i in range(dim1):U[i, :] = np.random.multivariate_normal(np.linalg.inv(var3[:, :, i]) @ var4[:, i], var3[:, :, i])#采樣，每一行是rank維度的一個多元正態分布向量return U

?采樣V同理：

def sample_factor_V(alpha_sparse_mat, alpha_ind, U, V, beta0 = 1):dim2, rank = V.shapeV_bar = np.mean(V, axis = 0)#每一列元素去一個均值，一共 rank個元素temp = dim2 / (dim2 + beta0)#M/M+β0var_mu_hyper = temp * V_bar#M X_bar/M+β0，由于μ0為0，所以這個也可以看成μ0*var_V_hyper = inv(np.eye(rank) + cov_mat(V, V) + temp * beta0 * np.outer(V_bar, V_bar))#np.eye(rank) 就是I，即我們的X0# cov_mat(V, V_bar) 就是N S_bar#np.outer(V_bar, V_bar) 相當于 V.reshape(dim2,1) @ V_bar.reshape(1,dim2)#temp * beta0 * np.outer(V_bar, V_bar)) 即求X0*的最后一項#所以這個可以看成X0*var_Lambda_hyper = wishart.rvs(df = dim2 + rank, scale = var_V_hyper)#df就是v0*#這一部分是從威沙特分布里面去一個[rank,rank]的樣本#也就是 Lambda_Vvar_mu_hyper = np.random.multivariate_normal(var_mu_hyper, (dim2 + beta0) * var_Lambda_hyper)#從N(μ0*，(β0* Lambda_X)……-1)中選取一個μV#直接multivariate_normal即可#[rank,1]var1 = U.T#[rank,dim1]var2 = kr_prod(var1, var1)#[ranl*rank,dim2]#alpha_ind [dim1,dim2]var3 = (var2 @ alpha_ind).reshape([rank, rank, dim2]) + var_Lambda_hyper[:, :, np.newaxis]#第一項相當于有值的部分相乘#var_Lambda_hyper——[rank,rank],這里給他強行添加了一個維度#第一項是12式的右邊，第二項是12式的左邊var4 = var1 @ alpha_sparse_mat + (var_Lambda_hyper @ var_mu_hyper)[:, np.newaxis]#采樣，每一行是一個多元正態分布for t in range(dim2):V[t, :] = np.random.multivariate_normal(np.linalg.inv(var3[:, :, t]) @ var4[:, t], var3[:, :, t])return V

?5 輔助函數

def compute_mape(var, var_hat):return np.sum(np.abs(var - var_hat) / var) / var.shape[0]def compute_rmse(var, var_hat):return np.sqrt(np.sum((var - var_hat) ** 2) / var.shape[0])def cov_mat(mat, mat_bar):mat = mat - mat_barreturn mat.T @ mat

6 BPMF

def BPMF(dense_mat, sparse_mat, init, rank, burn_iter, gibbs_iter):"""Bayesian Probabilistic Matrix Factorization, BPMF."""#dense_mat, ————沒有丟失值得矩陣#sparse_mat, ——有了丟失值得矩陣#init, ——初始的W和X矩陣#rank, ——特征矩陣的秩#burn_iter, ——表示經過了這么多步之后，MCMC代表的馬氏鏈趨近于平穩狀態#gibbs_iter——吉布斯采樣的步數dim1, dim2 = sparse_mat.shape#(214, 8784)U = init["U"]V = init["V"]ind = sparse_mat != 0#相當于ind = （sparse_mat != 0）pos_obs = np.where(ind)#輸出滿足條件 (即這個位置路段的速度非0) 元素的坐標。pos_test = np.where((dense_mat != 0) & (sparse_mat == 0))#在沒有丟失數據的矩陣上有值，但是我們人為丟掉的那些點的坐標dense_test = dense_mat[pos_test]#測試集，在沒有丟失數據的矩陣上有值，但是我們人為丟掉的那些點的坐標alpha = 2 #WX和Y之間的精度（協方差矩陣的倒數）U_plus = np.zeros((dim1, rank))V_plus = np.zeros((dim2, rank))#多次采樣的U和V的和temp_hat = np.zeros(sparse_mat.shape)#多次采樣的預測矩陣R的和（在平穩狀態之前）show_iter = 200mat_hat_plus = np.zeros(sparse_mat.shape)#多次采樣的預測矩陣R的和（在平穩狀態之后）for it in range(burn_iter + gibbs_iter):alpha_ind = alpha * ind #有觀測值的哪些點*alphaalpha_sparse_mat = alpha * sparse_mat #觀測矩陣*alphaU = sample_factor_U(alpha_sparse_mat, alpha_ind, U, V, alpha)#采樣一個W矩陣V = sample_factor_V(alpha_sparse_mat, alpha_ind, U, V)#采樣一個X矩陣mat_hat = U @ V.Ttemp_hat += mat_hatif (it + 1) % show_iter == 0 and it < burn_iter:temp_hat = temp_hat / show_iterprint('Iter: {}'.format(it + 1))print('MAPE: {:.6}'.format(compute_mape(dense_test, temp_hat[pos_test])))print('RMSE: {:.6}'.format(compute_rmse(dense_test, temp_hat[pos_test])))temp_hat = np.zeros(sparse_mat.shape)print()if it + 1 > burn_iter:U_plus += UV_plus += Vmat_hat_plus += mat_hatmat_hat = mat_hat_plus / gibbs_iterU = U_plus / gibbs_iterV = V_plus / gibbs_iterprint('Imputation MAPE: {:.6}'.format(compute_mape(dense_test, mat_hat[pos_test])))print('Imputation RMSE: {:.6}'.format(compute_rmse(dense_test, mat_hat[pos_test])))print()return mat_hat, U, V

7 求解

import time start = time.time() dim1, dim2 = sparse_mat.shape #(214, 8784) rank = 10 #低維特征矩陣的維度 init = {"U": 0.01 * np.random.randn(dim1, rank), "V": 0.01 * np.random.randn(dim2, rank)} burn_iter = 1 gibbs_iter = 2 mat_hat, U, V = BPMF(dense_mat, sparse_mat, init, rank, burn_iter, gibbs_iter) end = time.time() print('Running time: %d seconds'%(end - start))

《新程序員》：云原生和全面數字化實踐50位技術專家共同創作，文字、視頻、音頻交互閱讀

總結

以上是生活随笔為你收集整理的论文辅助笔记（代码实现）：Bayesian Probabilistic Matrix Factorizationusing Markov Chain Monte Carlo的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。