论文辅助笔记（代码实现）：Bayesian Probabilistic Matrix Factorizationusing Markov Chain Monte Carlo

最新推荐文章于 2023-01-02 00:48:05 发布

UQI-LIUWJ

最新推荐文章于 2023-01-02 00:48:05 发布

阅读量358

点赞数

分类专栏： python库整理文章标签： python 开发语言后端

本文链接：https://blog.csdn.net/qq_40206371/article/details/120914479

版权

python库整理专栏收录该内容

321 篇文章 44 订阅

订阅专栏

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.io

tensor = 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.2


binary_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,  观测矩阵*alpha
    alpha_ind, 示性矩阵*alpha
    U,  路段特征矩阵 [dim1,rank]
    V,  时间特征矩阵 [dim2,rank]
    alpha, #WX和Y之间的精度（协方差矩阵的倒数）
    beta0 = 1, 
    '''
    dim1, rank = U.shape
    U_bar = np.mean(U, axis = 0)
    #每一列元素取一个均值，一共 rank个元素，U_bar是U矩阵各个行向量的均值向量

    temp = dim1 / (dim1 + beta0)
    #N/N+β0

    var_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_U

    var_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.shape
    V_bar = np.mean(V, axis = 0)
    #每一列元素去一个均值，一共 rank个元素
    temp = dim2 / (dim2 + beta0)
    #M/M+β0
    var_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_V
    var_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_bar
    return 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 = 200
    mat_hat_plus = np.zeros(sparse_mat.shape)
    #多次采样的预测矩阵R的和（在平稳状态之后）

    for it in range(burn_iter + gibbs_iter):
        alpha_ind = alpha * ind #有观测值的哪些点*alpha
        alpha_sparse_mat = alpha * sparse_mat #观测矩阵*alpha
        U = 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.T
        temp_hat += mat_hat
        if (it + 1) % show_iter == 0 and it < burn_iter:
            temp_hat = temp_hat / show_iter
            print('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 += U
            V_plus += V
            mat_hat_plus += mat_hat
    mat_hat = mat_hat_plus / gibbs_iter
    U = U_plus / gibbs_iter
    V = V_plus / gibbs_iter
    print('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))