Personal Rank算法的原理及实现

1.简介

PageRank[1]是Google创始人拉里·佩奇和谢尔盖·布林于1997提出的用于客观评价网页的重要度的方法，最初的PageRank算法是主题无关的，它不依赖于任何特定的搜索查询。为了得到主题相关的搜索结果，Haveliwala提出主题敏感的PageRank方法，称为PersonalRank[2]，该方法用于在二分图中为用户进行推荐。二分图又称为二部图，是图论中的一种特殊模型，设G=(V,E) 是一个图，如果顶点V可分割为两个互不相交的子集(A, B)，并且图中的每条边(i, j)所关联的两个顶点i和j分别属于这两个不同的顶点集(i in A, j in B) ，则称图G为一个二分图，如图所示：

假设给“任小牛”进行个性化推荐，从节点“任小牛”开始游走，游走到一个节点时，首先alpha概率决定继续游走，或者以(1-alpha)的概率停止这次游走并从“任小牛”节点开始重新游走；如果决定继续游走，那么就从当前节点指向的节点中按照权重随机选择一个节点作为下次经过的节点，这样经过很多次的随机游走后，每个节点被访问到的概率就会收敛。最终推荐列表中节点的权重就是节点的访问概率，PersonalRank方法公式如下1-1所示，其中PR(j)表示物品j的访问概率，out(i) 是物品节点i的出度，alpha决定继续访问的概率。 

但是，迭代形式的Personal Rank算法计算复杂度较高，它需要经过多次的迭代游走，才能使得各节点的重要度趋于稳定，其改进方案是经过一次矩阵运算直接得到系统的稳态，公式1-2为上式的矩阵表示形式： 

​​

其中r是n维向量，每个元素表示一个节点的PR重要度， 也是个n维向量，第i个位置上是1，其余位置均为0，表示对第i个节点进行推荐，M 是n阶转移矩阵，定义见公式1-3：

由公式1-2变形可以得到公式1-4和公式1-5：

基于公式1-5，解一次线性方程组即可以得到r的值，对r中的各元素降序排列即为节点i的推荐列表。

2.实现

基于Python语言设计并实现了PersonalRank类，它的逻辑较为简单，见流程图2-1所示：

在Personal Rank类中，不同实现的输入和输出都是一致的，输入都是networkx库DiGraph（有向图），输出则是对应节点的推荐列表，只不过有着三种不同的实现：

1. train()方法基于迭代游走；
2. train_matrix()方法基于numpy库进行矩阵运算；
3. train_csr_matrix()基于scipy的稀疏矩阵进行矩阵运算。

对于Personal Rank算法来说，迭代游走是最为直观的方式，也是最耗时的实现方式。

首先构造图，即PersonalRank类的输入：

if __name__ == '__main__':graph = networkx.DiGraph()graph.add_edge('任小牛', '笔记本电脑', weight=1)graph.add_edge('任小牛', '风扇', weight=0.1)graph.add_edge('任小牛', '键盘', weight=0.1)graph.add_edge('卡洛斯', '笔记本电脑', weight=0.2)graph.add_edge('卡洛斯', '风扇', weight=0.3)graph.add_edge('詹姆斯', '风扇', weight=0.4)graph.add_edge('詹姆斯', '键盘', weight=0.5)graph.add_edge('卡尔', '笔记本电脑', weight=0.7)graph.add_edge('卡尔', '键盘', weight=0.9)rank = PersonalRank(alpha=0.85)target = '任小牛'rs = rank.train(graph, target, 2000)rs = sorted(rs.items(), key=lambda x: x[1], reverse=True)print(rs)# 另一个rs = rank.train_matrix(graph, target)rs = sorted(rs.items(), key=lambda x: x[1], reverse=True)print(rs)# gmresrs = rank.train_csr_matrix(graph, target)rs = sorted(rs.items(), key=lambda x: x[1], reverse=True)print(rs)

本文实现的PersonalRank类的输入是有向有权图，有向是为了权重服务的，而权重则是根据需求自己定义的，一般情况下，权重应该小于等于1。

接着是构造函数，alpha的作用表示转移概率：

import time
import networkx
import functools
import numpy as np
from numpy.linalg import solve
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import gmres, lgmresclass PersonalRank(object):def __init__(self, alpha):self.alpha = alpha

然后是train方法的实现：

    @log_timedef train(self, graph, root, iterations):rank = {x: 0 for x in graph.nodes}rank[root] = 1count = 0# 迭代while True:tmp = {x: 0 for x in graph.nodes}# 节点i和它的出度节点集合rifor node in graph.nodes:out_nodes = list(graph.neighbors(node))out_degree = len(out_nodes)# 节点j和边权重for j in out_nodes:data = graph.get_edge_data(node, j)tmp[j] += self.alpha * rank[node] / out_degree * data['weight']# 每次游走都是从root出发，因此root节点的权重需要加上1- alphatmp[root] += (1 - self.alpha)rank = tmpcount += 1if count >= iterations:# print('PersonalRank:%d' % count)breakreturn rank

train_matrix方法和train_csr_matrix方法基于矩阵计算得到推荐结果，前者使用numpy，后者使用了scipy的稀疏矩阵，它们比迭代游走多出的一个获取权重矩阵的步骤：

def _get_weight_matrix(graph, is_using_out=True):nodes = list(graph.nodes)matrix = networkx.adjacency_matrix(graph).Aif is_using_out:degrees = [1.0/graph.degree(node) for node in nodes]degrees = np.tile(degrees, (len(nodes), 1))matrix = csr_matrix(np.multiply(matrix, degrees))return matrix

_get_weight_matrix()方法较为简单，因为使用了networkx的缘故，因此直接调用networkx.adjacency_matrix方法得到权重矩阵，之后的判断表示是否使用出度均分权重。

class PersonalRank(object):...@log_timedef train_matrix(self, graph, root):"""使用矩阵求解:param graph: networkx.DiGraph:param alpha: 系数，一般为0.8左右:param root: 从哪出发:return:"""# 生成矩阵nodes = list(graph.nodes)matrix = _get_weight_matrix(graph)# 除了边的权重外，还需要为每个点乘以出度的倒数r0 = np.matrix([[1 if node == root else 0] for node in nodes])n = matrix.shape[0]# 求解A = np.eye(n) - self.alpha * matrix.Tb = (1 - self.alpha) * r0r = solve(A, b)rank = {}for j in range(n):rank[nodes[j]] = r[j, 0]return rank@log_timedef train_csr_matrix(self, graph, root):# 生成矩阵nodes = list(graph.nodes)matrix = _get_weight_matrix(graph)n = matrix.shape[0]r0 = np.matrix([[1 if node == root else 0] for node in nodes])A = np.eye(n) - self.alpha * matrix.Tb = (1 - self.alpha) * r0r = lgmres(A, b, tol=1e-8, atol=1e-8, maxiter=1)[0]rank = {}for j in range(n):rank[nodes[j]] = r[j]return rank

train方法、train_matrix方法和train_csr_matrix方法都使用了log_time装饰器来计算执行的时间，实现如下：

def log_time(func):@functools.wraps(func)def wrapper(*args, **kwargs):start = time.time()rs = func(*args, **kwargs)print('%s:duration:%.6f' % (func.__name__, time.time() - start))return rsreturn wrapper

最后则是查看执行结果：

train:duration:0.014966
[('任小牛', 0.15000000000000002), ('笔记本电脑', 0.0425), ('风扇', 0.00425), ('键盘', 0.00425), ('卡洛斯', 0), ('詹姆斯', 0), ('卡尔', 0)]
train_matrix:duration:0.000994
[('任小牛', 0.15000000000000002), ('笔记本电脑', 0.0425), ('风扇', 0.00425), ('键盘', 0.00425), ('卡洛斯', 0.0), ('詹姆斯', 0.0), ('卡尔', 0.0)]
train_csr_matrix:duration:0.000996
[('任小牛', 0.15000000000000002), ('笔记本电脑', 0.042499999999999996), ('风扇', 0.0042499999999999994), ('键盘', 0.0042499999999999994), ('卡洛斯', 0.0), ('詹姆斯', 0.0), ('卡尔', 0.0)]

注：本示例中只有从“人”到“物品”的边，因此对于没有“人”->“物品”边的推荐分数为0。

从执行结果来看，三个方法的推荐结果基本一致；从执行时间上来看，迭代游走花费的时间是其他两个方法的15倍；而目前对于后两个方法来看，执行结果基本一致，但是当网络变大时，第三个方法的执行时间要小于第二个方法。

3.实例验证

本文的实例验证各个实现所花费的时间，其场景为基于《战略性新兴产业》为地区推荐高校，它们之间的关系见图3-1所示：

暂且不提单位和战略性新兴产业的边及边的权重的由来，若选择近10年全国发表专利最多的150个高校，为地区推荐高校，那么在构造好有向有权图后，得到推荐列表见图3-2：

在本次示例中，地区、高校和战略性新兴产业共1460个节点，其花费时间见表3-1：

4.参考文献

1. Page, Lawrence, Brin, et al. Page, L. et al.: The PageRank citation ranking: Bringing order to the web. stanford digital libraries working paper, 1998.
2. Haveliwala T H. Topic-sensitive PageRank[C]// International world wide web conference. ACM, 2002.

5.参考帖子

个性化推荐召回算法——Personal Rank

基于图的推荐算法及Python实现（PersonalRank）

PersonalRank：一种基于图的推荐算法