浅谈贝叶斯张量分解(一):贝叶斯矩阵分解

贝叶斯张量分解(Bayesian Probabilistic Tensor Decomposition)是一种有效的张量分解方法,与其他张量分解方法不同的是,该方法能够很好地利用贝叶斯推断(Bayesian inference)进行参数估计,依据贝叶斯准则,只要写出参数(或超参数)的后验分布形式就能够对参数(或超参数)进行采样更新

由于矩阵分解往往作为张量分解的基础,因此,以下将以正态分布假设为基础,简单地介绍一下贝叶斯矩阵分解。限于本文篇幅较大,这里先列出了文章的大体结构:


1. 贝叶斯矩阵分解原理 (以矩阵的UV分解为例,介绍了一个高斯分布假设下的贝叶斯网络。)

2. 参数、超参数及待估计参数(介绍了贝叶斯网络中的参数划分。)

3. 参数和超参数的后验分布(详细介绍了如何依据贝叶斯准则推导出后验分布。)

3.1. 参数\boldsymbol{u}_i 的后验分布

3.2. 参数 \tau_{\epsilon} 的后验分布

3.3. 超参数 \boldsymbol{\lambda} 的后验分布

4. 算法实现和相关小实验(选取了一个公共开放的城市路网数据集,并采用贝叶斯矩阵分解测试缺失数据修复效果。)

4.1. 城市路网车速数据集

4.2. 用Matlab实现一个简单的贝叶斯矩阵分解


1. 贝叶斯矩阵分解原理

在矩阵分解中,同样可以利用贝叶斯推断进行分解结构的参数估计,即贝叶斯矩阵分解(Bayesian Probabilistic Matrix Factorization, BPMF),图1给出了一个矩阵分解的贝叶斯网络。

图1 矩阵分解的贝叶斯网络简图

对于一个稀疏矩阵 X\in\mathbb{R}^{m\times n} ,假设矩阵分解结构为

X\approx UV^T

任意第 i 行第 j 列的元素满足

x_{ij}\approx \boldsymbol{u}_{i}^T\boldsymbol{v}_j

其中,因子矩阵(factor matrix) U 的大小为 m\times r ,因子矩阵 V 的大小为 n\times r 。为方便后面的推导,用同样大小的矩阵 B 表示 X 内的元素是否发生缺失,若矩阵 X 的第i 行 第j 列元素缺失,则记作 b_{ij}=0 ,否则 b_{ij}=1

另外,需要特别说明的是,本文所提及的向量均为列向量,以图1为例,其中,\boldsymbol{u}_i\in\mathbb{R}^{r} 是取自矩阵 U\in\mathbb{R}^{m\times r}i 行的向量,为统一表示,将其记为列向量。

在图1中,观测值 x_{ij} 被假设为正态分布,其形式为

x_{ij}\sim\mathcal{N}\left(\boldsymbol{u}_i^T\boldsymbol{v}_j,\tau_{\epsilon}^{-1}\right)

其中, 符号\mathcal{N}(\cdot) 表示正态分布,其中, \tau_{\epsilon}^{-1} 表示方差, \tau_{\epsilon} 的意义是精度项(precision),根据共轭性质(conjugate prior),假设 \tau_{\epsilon}\sim\text{Gamma}\left(a_0,b_0\right) ,其中,伽马分布的形式为

p\left(\tau_{\epsilon}\mid a_0,b_0\right)=\frac{\left(b_0\right)^{a_0}}{\Gamma\left(a_0\right)}\left(\tau_{\epsilon}\right)^{a_0-1}\exp\left(-b_0\tau_{\epsilon}\right)

在伽马分布中, a_0 是形状参数(shape parameter), b_0 是比率参数(rate parameter).

在因子矩阵中,向量 \boldsymbol{u}_i\boldsymbol{v}_j 都被假设服从多元正态分布,如下:

\boldsymbol{u}_i\sim\mathcal{N}\left(\boldsymbol{0},\left[\text{diag}\left(\boldsymbol{\lambda}\right)\right]^{-1}\right)

\boldsymbol{v}_j\sim\mathcal{N}\left(\boldsymbol{0},\left[\text{diag}\left(\boldsymbol{\lambda}\right)\right]^{-1}\right)

其中,由于向量 \boldsymbol{u}_i,\boldsymbol{v}_j\in\mathbb{R}^{r} ,因此,为简化贝叶斯网络的结构,设定多元正态分布的均值向量为 \boldsymbol{0}\in\mathbb{R}^{r} ,给定一个向量 \boldsymbol{\lambda}\in\mathbb{R}^{r} ,其对角化后的矩阵便可转化为多元正态分布的协方差矩阵,更进一步,向量 \boldsymbol{\lambda} 任意元素被假设服从伽马分布,即 \lambda_k\sim\text{Gamma}\left(c_0,h_0\right) , k=1,2,...,r .

2. 参数、超参数及待估计参数

在图1中,所涉及到的变量可以归为三类,第一类是输入的观测值 x_{ij} (observations);第二类是模型的参数(model parameters),即 \left\{\boldsymbol{u}_i,\boldsymbol{v}_j,\tau_{\epsilon}\right\} ;第三类则是超参数(hyper-parameters),向量 \boldsymbol{\lambda} 以及需要在迭代过程中预先设置的 \left\{a_0,b_0,c_0,h_0\right\} 都属于该类。

在这里,依据贝叶斯准则(后验正比于似然与先验的乘积),我们需要推导并写出 \left\{\boldsymbol{u}_i,\boldsymbol{v}_j,\tau_{\epsilon},\boldsymbol{\lambda}\right\} (可视为待估计参数)的后验分布,然后再进行采样和更新。

3. 参数和超参数的后验分布

在这部分,我们将讨论如何利用Gibbs采样(不太熟悉Gibbs采样原理的读者可参考浅谈「Gibbs采样」一文)对模型参数和超参数进行更新。

------------***------------

3.1. 参数\boldsymbol{u}_i 的后验分布

从图1可以看到,与 \boldsymbol{u}_i 关联的变量有观测值 x_{ij} 和超参数 \boldsymbol{\lambda} ,根据“箭头”指向,我们能够很容易地写出似然函数(likelihood)

\mathcal{L}\left(\boldsymbol{x}_i\mid\boldsymbol{u}_i,V,\tau_{\epsilon}\right)\propto\exp\left\{-\frac{\tau_{\epsilon}}{2}\sum_{j=1}^{n}b_{ij}(x_{ij}-\boldsymbol{u}_i^T\boldsymbol{v}_j)^2\right\}

\propto\exp\left\{-\frac{\tau_{\epsilon}}{2}\left[\boldsymbol{b}_i^T\odot\left(\boldsymbol{x}_i^T-\boldsymbol{u}_i^TV^T\right)\right]\left[\boldsymbol{b}_i^T\odot\left(\boldsymbol{x}_i^T-\boldsymbol{u}_i^TV^T\right)\right]^T\right\}

\propto\exp\left\{-\frac{\tau_{\epsilon}}{2}\boldsymbol{u}_i^TA_iA_i^T\boldsymbol{u}_i+\frac{\tau_{\epsilon}}{2}\boldsymbol{u}_i^TA_i\boldsymbol{x}_i\right\}

在这里,为方便书写,记作 A_i=\boldsymbol{b}_i^T\odot V^T\in\mathbb{R}^{r\times n} ,符号\odot 表示Khatri-Rao积(运算规则可参考浅谈张量分解(二):张量分解的数学基础一文), \boldsymbol{x}_i,\boldsymbol{b}_i\in\mathbb{R}^{n} 均为列向量。

与此同时,我们也可以写出先验分布的形式为

p\left(\boldsymbol{u}_i\mid\boldsymbol{\lambda}\right)\propto\exp\left\{-\frac{1}{2}\boldsymbol{u}_i^T\text{diag}\left(\boldsymbol{\lambda}\right)\boldsymbol{u}_i\right\}

由此,依据贝叶斯准则,后验分布为

p\left(\boldsymbol{u}_i\mid\boldsymbol{x}_i,V,\tau_{\epsilon},\boldsymbol{\lambda}\right)\propto\mathcal{L}\left(\boldsymbol{x}_i\mid\boldsymbol{u}_i,V,\tau_{\epsilon}\right)p\left(\boldsymbol{u}_i\mid\boldsymbol{\lambda}\right)

\propto\exp\left\{-\frac{1}{2}\boldsymbol{u}_i^T\left[\tau_{\epsilon}A_iA_i^T+\text{diag}\left(\boldsymbol{\lambda}\right)\right]\boldsymbol{u}_i+\frac{\tau_{\epsilon}}{2}\boldsymbol{u}_i^TA_i\boldsymbol{x}_i\right\}

并且满足 \boldsymbol{u}_i\sim\mathcal{N}\left(\tilde{\boldsymbol{\mu}}_i^{(1)},\tilde{\Lambda}_i^{(1)}\right) ,其中,

\tilde{\Lambda}_i^{(1)}=\left[\tau_{\epsilon}A_iA_i^T+\text{diag}\left(\boldsymbol{\lambda}\right)\right]^{-1}

\tilde{\boldsymbol{\mu}}_i^{(1)}=\tau_{\epsilon}\tilde{\Lambda}_i^{(1)}A_i\boldsymbol{x}_i

是更新后的多元正态分布参数,有了这些参数即可进行采样,并可将采样后得到的向量用于更新原来的 \boldsymbol{u}_i 。另外,参数 \boldsymbol{v}_j 的后验分布与 \boldsymbol{u}_i 类似,感兴趣的读者可以自行推导。

------------***------------

3.2. 参数 \tau_{\epsilon} 的后验分布

与参数 \boldsymbol{u}_i 求后验分布类似, \tau_{\epsilon} 的后验分布为

p\left(\tau_{\epsilon}\mid-\right)\propto\prod_{i=1}^{m}\prod_{j=1}^{n}p\left(x_{ij}\mid\boldsymbol{u}_i,\boldsymbol{v}_j,\tau_{\epsilon}\right)p\left(\tau_{\epsilon}\mid a_0,b_0\right)

其中,为简化书写,后验分布 p\left(\tau_{\epsilon}\mid x_{ij},\boldsymbol{u}_i,\boldsymbol{v}_j,\tau_{\epsilon},a_0,b_0\right) 被记作 p\left(\tau_{\epsilon}\mid-\right) ,进一步,后验分布可以被写成如下形式

p\left(\tau_{\epsilon}\mid-\right)

\propto\prod_{i=1}^{m}\prod_{j=1}^{n}b_{ij}\left(\tau_{\epsilon}\right)^{1/2}\exp\left\{-\frac{\tau_{\epsilon}}{2}\left(x_{ij}-\boldsymbol{u}_i^T\boldsymbol{v}_j\right)^2\right\}\left(\tau_{\epsilon}\right)^{a_0-1}e^{-b_0\tau_{\epsilon}}

\propto\left(\tau_{\epsilon}\right)^{a_0+\frac{1}{2}\sum_{i,j}b_{ij}-1}\exp\left\{-\left(b_0+\frac{1}{2}\sum_{i,j}b_{ij}\left(x_{ij}-\boldsymbol{u}_i^T\boldsymbol{v}_j\right)^2\right)\tau_{\epsilon}\right\}

\tau_{\epsilon}\sim\text{Gamma}\left(\tilde{a},\tilde{b}\right) ,其中

\tilde{a}=a_0+\frac{1}{2}\sum_{i,j}b_{ij}

\tilde{b}=b_0+\frac{1}{2}\sum_{i,j}b_{ij}\left(x_{ij}-\boldsymbol{u}_i^T\boldsymbol{v}_j\right)^2 .

------------***------------

3.3. 超参数 \boldsymbol{\lambda} 的后验分布

尽管超参数 \boldsymbol{\lambda} 与参数 \tau_{\epsilon} 都被假设服从伽马分布,但不同的是,参数 \boldsymbol{\lambda} 作为一个向量,对应着多元正态分布中的协方差矩阵,在这里,不妨以 \boldsymbol{u}_i 为例,先写一下多元正态分布的形式

p\left(\boldsymbol{u}_i\mid\boldsymbol{\lambda}\right)=\frac{|\text{diag}\left(\boldsymbol{\lambda}\right)|^{1/2}}{\left(2\pi\right)^{r/2}}\exp\left\{-\frac{1}{2}\boldsymbol{u}_i^T\text{diag}\left(\boldsymbol{\lambda}\right)\boldsymbol{u}_i\right\}

从这条公式中,对于任意 \lambda_k,k=1,2,...,r ,有

p\left({u}_{ik}\mid\lambda_k\right)\propto\left(\lambda_k\right)^{1/2}\exp\left\{-\frac{1}{2}\lambda_k{u}_{ik}^2\right\}

因此,关于 \lambda_k 的似然函数为

\mathcal{L}\left(U,V\mid\lambda_k\right)=\prod_{i=1}^{m}p\left({u}_{ik}\mid\lambda_k\right)\prod_{j=1}^{n}p\left(v_{jk}\mid\lambda_k\right)

\propto\left(\lambda_k\right)^{\left(m+n\right)/2}\exp\left\{-\frac{1}{2}\lambda_k\left(\sum_{i=1}^{m}{u}_{ik}^2+\sum_{j=1}^{n}{v}_{jk}^2\right)\right\}

从而,我们可以推导出后验分布 p\left(\lambda_k\mid U,V,c_0,h_0\right) ,即

p\left(\lambda_k\mid U,V,c_0,h_0\right)\propto \mathcal{L}\left(U,V\mid\lambda_k\right) p\left(\lambda_k\mid c_0,h_0\right)

\propto\left(\lambda_{k}\right)^{c_0+\left(m+n\right)/2-1}\exp\left\{-\left[h_0+\frac{1}{2}\left(\sum_{i=1}^{m}{u}_{ik}^2+\sum_{j=1}^{n}{v}_{jk}^2\right)\right]\lambda_{k}\right\}

更新后的超参数 \tilde{c},\tilde{h} 分别为

\tilde{c}=c_0+\frac{1}{2}\left(m+n\right)

\tilde{h}=h_0+\frac{1}{2}\left(\sum_{i=1}^{m}{u}_{ik}^2+\sum_{j=1}^{n}{v}_{jk}^2\right)

超参数 \lambda_k\sim\text{Gamma}\left(\tilde{c},\tilde{h}\right),k=1,2,...,r .

------------***------------

4. 算法实现和相关小实验

到这里,我们已经了解到贝叶斯矩阵分解的基本原理,并学会了如何推导参数(和待估计超参数)的后验分布,但在实际应用中,我们该如何实现贝叶斯矩阵分解呢?

4.1. 城市路网车速数据集

下面将以路网车速数据为例来具体看看贝叶斯矩阵分解如何用于修复缺失数据,感兴趣的读者可以下载该城市路网车速数据集(Urban Traffic Speed Dataset of Guangzhou, China)。


图2 广州城市路网车速数据集


数据集中的张量(tensor.mat)可在Matlab中直接打开,张量的大小为 214\times 61\times 144 ,其中,214表示214条匿名路段,61对应着2016年8月1日至9月30日的61天,144对应着144个时间窗,即每天的1440分钟按每10分钟进行切分恰好可以得到144个时间窗口。需要注意的是,该车速数据集的观测样本为185.56万,原始缺失率为1.29%.

4.2. 用Matlab实现一个简单的贝叶斯矩阵分解

将该张量展开成大小为 214\times 8784 的矩阵,其中,张量展开的函数在之前的文章(浅谈张量填补算法——HaLRTC)中已经给出。

clear;
load tensor.mat;
OriMat = unfolding(tensor,214,61,144,1); % original matrix

紧接着,生成一个0-1之间均匀分布的随机矩阵,调整该随机矩阵元素数值的大小可用以控制缺失率,以下设置的缺失率为20%

dim = size(OriMat);
RandMat = rand(dim); % uniform distributed random numbers
RandMat0 = RandMat+0.5-0.2;
SparseMat = round(RandMat0).*OriMat; % sparse matrix with 20% missing ratio
pos1 = find(SparseMat>0);
pos2 = find(OriMat>0 & SparseMat==0);
BinMat = SparseMat;
BinMat(pos1) = 1; % binary matrix

设置初始化参数的大小。

r = 40;
maxiter = 200;
FactMat = cell(2,1);
for k = 1:2
    FactMat{k} = 0.1*randn(dim(k),r);
end
tau = 1;
a0 = 1e-6;
b0 = 1e-6;
lambda = ones(r,1);
c0 = 1e-6;
h0 = 1e-6;

根据前面已经推导的后验分布进行迭代采样。

第一步:对因子矩阵 U,V 进行采样更新。

rmse = zeros(maxiter,1);
for iter = 1:maxiter
    % Update factor matrices
    for i = 1:dim(1)
        pos = find(SparseMat(i,:)>0);
        var = FactMat{2}(pos,:);
        Lambda = (tau*var'*var+diag(lambda))^(-1);
        Lambda = (Lambda+Lambda')./2;
        mu = tau*Lambda*var'*SparseMat(i,pos)';
        FactMat{1}(i,:) = mvnrnd(mu,Lambda);
    end
    for i = 1:dim(2)
        pos = find(SparseMat(:,i)>0);
        var = FactMat{1}(pos,:);
        Lambda = (tau*var'*var+diag(lambda))^(-1);
        Lambda = (Lambda+Lambda')./2;
        mu = tau*Lambda*var'*SparseMat(pos,i);
        FactMat{2}(i,:) = mvnrnd(mu,Lambda);
    end
    MatHat = FactMat{1}*FactMat{2}';
    rmse(iter,1) = sqrt(sum((OriMat(pos2)-MatHat(pos2)).^2)./length(pos2));

第二步:对超参数 \boldsymbol{\lambda} 进行采样更新。

    % Update lambda
    c = (0.5*sum(dim)+c0)*ones(r,1);
    h = 0;
    for k = 1:2
        h = h+diag(FactMat{k}'*FactMat{k});
    end
    h = h0+0.5*h;
    lambda = gamrnd(c,1./h);

第三步:对精度项 \tau_{\epsilon} 进行采样更新,同时输出当前迭代的误差。

    % Update precision
    a = a0+0.5*length(pos1);
    error = SparseMat-MatHat;
    b = b0+0.5*sum(error(pos1).^2);
    tau = gamrnd(a,1./b);
    
    % Print the results
    fprintf('iteration = %g, RMSE = %g km/h.\n',iter,rmse(iter));
end


到这里,我们已经完成了一个贝叶斯矩阵分解的实验。需要注意的是,为了提高程序的可读性,这里所给出的代码在运算效率方面表现一般(迭代200代大约需要40分钟)。以下截取了迭代过程中的误差(RMSE=4.49km/h,如图3)用于表征缺失数据修复的效果。


图3 贝叶斯矩阵分解计算结果截图

编辑于 2018-05-22