Salakhutdinov, Ruslan, and A. Mnih. “Bayesian probabilistic matrix factorization using markov chain monte carlo.” International Conference on Machine Learning 2008:880-887.

对PMF模型以及应用场景不熟悉的同学可以先阅读这篇基础PMF。
本论文的模型和前文类似，但在求解时，是从贝叶斯角度而不是传统概率角度出发：不再把系统参数当做一个固定值估计，而是作为一个服从某种分布的随机变量，转而估计该分布的参数。

问：明天什么天气？ 传统观点估计：呃，晴天吧。 贝叶斯观点估计：60%晴天，40%下雨。传统观点估计：你个心机婊。

模型

用 N×M 关系矩阵 R 来描述 N 个用户给 M 部电影的打分。这里 R 为观测值。

p(R)∼N(R|UTV,α−1)=∏ijN(Rij|UTiVj,α−1)

矩阵高斯函数等于其元素高斯函数的乘积。其中标量 α 为精度，等价于协方差的倒数（这样写是为了后续计算方便）。
R 的各个位置（即各次观测）的精度相同。

U,V 尺寸为 N×D,M×D ，也服从高斯分布。（在基础PMF中，U,V都是待求的固定值）

p(U)∼N(U|μU,Λ−1U)=∏iN(Ui|μU,Λ−1U)

p(V)∼N(V|μV,Λ−1V)=∏jN(Vj|μV,Λ−1V)

所有用户共享一组超参数 μU,ΛU ，尺寸为 D×1,D×D 。所有电影同样共享一组超参数 μV,ΛV 。
超参数 μ,Λ 服从Gaussian-Wishart分布，表示为一个均值的高斯分布，以及一个协方差的威沙特分布的乘积。

p(μ,Λ)∼N(μ|μ0,(β0Λ)−1)⋅W(Λ|W0,ν0)

其中 μ0=0,W0=I,ν0=D 。为书写简便，记为 ΘU={μU,ΛU} , ΘV={μV,ΛV} 。

预测

用 R∗ 表示预测出的关系矩阵。在贝叶斯观点下，推荐问题转化为：已知观测矩阵 R ，求预测矩阵的概率 p(R∗|R) 。
这个概率本身非常难求，所以退 而求其次：按照概率 p(R∗|R) ，生成预测 R∗ 即可，多输出几次在统计意义上就准确了。

在已知 U,V 的情况下， R∗ 的分布参数 α 已知。问题转化为：生成符合概率分布 p(U|R)，p(V|R) 的 U,V 。

核心变量是特征 U,V ,以及超参数 ΘU,ΘV 。设定初始值 U1,V1 。使用MCMC方法，按照如下步骤逐步逼近稳态分布
– 固定 U,V ，按照 p(ΘU|U),p(ΘV|V) 对 ΘU,ΘV 进行采样
– 固定 ΘU,V 以及观测 数据 R ，根据 p(U|R,V,ΘU) ，对 U 采样
– 固定 ΘV,U 以及观测 数据 R ，根据 p(V|R,U,ΘV) ，对 V 采样
每完成这样一轮，称为一个Epoch。随着epoch增加，生成的预测矩阵 R 的分布逐渐准确。

在迭代的第一个步骤中，由于已知 U 时， R 对 ΘU 无影响，所以条件中没有 R 。

根据特征，更新超参数

固定 U ，更新 ΘU ，要求的条件概率为 p(ΘU|U)
利用条件概率性质，凑出两项形式已知的条件概率：

p(ΘU|U)⋅p(U)=p(ΘU,U)=p(U|ΘU)⋅p(ΘU)

由于

U 已知，故p(U)为常量：

p(ΘU|U)∼p(U|ΘU)⋅p(ΘU)

已知p(A|B)，欲求p(B|A)。经典推导： p(A|B)p(B)=p(B|A)p(A)

右边两项均有明确的概率形式，但如果直接计算运算量太大，以下将其转换为简单的参数形式。
右侧第一项是一个高斯分布，第二项是高斯分布和威沙特分布的乘积。由于高斯乘积仍为高斯，两项乘积仍然是高斯分布和威沙特分布的乘积。需要求其参数 μ∗0,ν∗0,W∗0,β∗0 ，维度为 D 。

考察上式右侧的两个高斯分布和一个威沙特分布，写出和 μU,ΛU 有关的项：

A1=ΛN/2U⋅exp[−12∑i=1N(Ui−μU)TΛU(Ui−μU)]

A2=(β0ΛU)1/2⋅exp[−12(μU−μ0)Tβ0ΛU(μU−μ0)]

A3=|ΛU|(ν0−D−1)/2⋅exp[−12Tr(W−10ΛU)]

对比标准高斯威沙特分布：

B1=(β∗0Λ)1/2⋅exp[−12(μU−μ∗0)Tβ∗0ΛU(μU−μ∗0)]

B2=|ΛU|(ν∗0−D−1)/2⋅exp[−12Tr((W∗0)−1ΛU)]

应该找出这样的 μ∗0,ν∗0,W∗0 ，使得 A1⋅A2⋅A3=B1⋅B2

考察 ΛU 的指数：

N+1+ν0−D−1=1+ν∗0−D−1⇒ν∗0=ν0+N

考察指数中 μTUΛUμU 项：

β∗0=N+β0

考察指数中 ΛUμU 项：

∑i=1NUi+β0μ0=β∗0μ∗0⇒μ∗0=∑Ni=1Ui+β0μ0β0+N

W∗0 我实在导不出来了，哪位神指教一下论文里这是为啥：

根据超参数，更新特征

固定 V,R,ΘU ，更新 U 。要求的条件概率是 p(U|R,V,ΘU) 。
同样利用条件概率性质：

p(U|R,V,ΘU)⋅p(R|V,ΘU)=p(U,R|V,ΘU)

=p(R|U,V,ΘU)⋅p(U|V,ΘU)

在这一步中， R,V,ΘU 已知，所以等号左边第二项为常数； R 和 ΘU 无关， U,V 无关，等号右边相应删减：

p(U|R,V,ΘU)⋅∼p(R|U,V)⋅p(U|ΘU)

由于右侧两个高斯函数的乘积仍然是高斯函数，故待求概率也是高斯分布。需要求其参数 μ∗,Λ∗ 。

考察上式指数部分：

−∑ij[α(Rij−UTiVj)2]−∑i(UTi−μTU)ΛU(U−μU)+C

把指数部分和 U 有关的项提出来：

二次项系数=−12⎡⎣∑ijUTiVj⋅α⋅VTjUi+∑iUTiΛUUi⎤⎦

一次项系数=∑ijUTiαRijVj+∑iUTiΛUμU

对比一个经典高斯分布 N(μ,α) 的指数部分 −12(x−μ)Tα(x−μ)

二次项部分=−12xTΛx

一次项部分=xTΛμ

两相对比，可以推导出超参数的取值：

Λ∗=α⋅∑jVjVTj+ΛU

μ∗=(Λ∗)−1⋅⎛⎝α∑jRijVj+ΛUμU⎞⎠

从物理意义上来说：由于知道了 V,R ，对于 U 的了解更多了，所以精度增加了。

这种在指数部分“凑平方”(complete the square)的方法是处理高斯函数的经典方法。在《Pattern Recognition and Machine Learning》的2.3.1章有着重介绍。

p(V|R,U,ΘV) 用同样方法求得高斯参数。直接对高斯函数采样即可。

实验

实验数据和基础PMF相同。由于使用了贝叶斯估计，降低了过拟合现象，精度有所提高（1.7%）。
D = 10，每轮采样需要6.6分钟；D=300，则需要220分钟。