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

上面的二部图表示user A对item a和c感兴趣，B对a b c d都感兴趣，C对c和d感兴趣。本文假设每条边代表的感兴趣程度是一样的。

现在我们要为user A推荐item，实际上就是计算A对所有item的感兴趣程度。在personal rank算法中不区分user节点和item节点，这样一来问题就转化成：对节点A来说，节点A B C a b c d的重要度各是多少。重要度用PR来表示。

初始赋予 PR(A)=1,PR(B)=PR(C)=PR(a)=PR(b)=PR(c)=PR(d)=0，即对于A来说，他自身的重要度为满分，其他节点的重要度均为0。

然后开始在图上游走。每次都是从PR不为0的节点开始游走，往前走一步。继续游走的概率是α，停留在当前节点的概率是1−α。

第一次游走， 从A节点以各自50%的概率走到了a和c，这样a和c就分得了A的部分重要度，PR(a)=PR(c)=α∗PR(A)∗0.5。最后PR(A)变为1−α。第一次游走结束后PR不为0的节点有A a c。

第二次游走，分别从节点A a c开始，往前走一步。这样节点a分得A 12∗α的重要度，节点c分得A 12∗α的重要度，节点A分得a 12∗α的重要度，节点A分得c 13∗α的重要度，节点B分得a 12∗α的重要度，节点B分得c 13∗α的重要度，节点C分得c 13∗α的重要度。最后PR(A)要加上1−α

经过以上推演，可以概括成公式：

(1)PR(j)={α∗∑i∈in(j)PR(i)|out(i)|    if(j≠u)(1−α)+α∗∑i∈in(j)PR(i)|out(i)|    if(j=u)

u是待推荐的用户节点。

personalrank.py

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

#coding=utf-8

__author__ = 'orisun'

 

import time

 

 

def PersonalRank(G, alpha, root, max_step):

    rank = dict()

    rank = {x:0 for x in G.keys()}

    rank[root] = 1

    #开始迭代

    begin=time.time()

    for k in range(max_step):

        tmp = {x:0 for x in G.keys()}

        #取节点i和它的出边尾节点集合ri

        for i, ri in G.items():

            #取节点i的出边的尾节点j以及边E(i,j)的权重wij, 边的权重都为1，归一化之后就上1/len(ri)

            for j, wij in ri.items():

                #i是j的其中一条入边的首节点，因此需要遍历图找到j的入边的首节点，

                #这个遍历过程就是此处的2层for循环，一次遍历就是一次游走

                tmp[j] += alpha * rank[i] / (1.0 * len(ri))

        #我们每次游走都是从root节点出发，因此root节点的权重需要加上(1 - alpha)

        tmp[root] += (1 - alpha)

        rank = tmp

    end=time.time()

    print 'use time', end - begin

 

    li = sorted(rank.items(), cmp=lambda x, y: cmp(x[1], y[1]), reverse=True)

    for ele in li:

        print "%s:%.3f, \t"%(ele[0], ele[1]),

    print

     

    return rank

 

 

if __name__ == '__main__' :

    alpha = 0.8

    G = {'A' : {'a' : 1, 'c' : 1},

         'B' : {'a' : 1, 'b' : 1, 'c':1, 'd':1},

         'C' : {'c' : 1, 'd' : 1},

         'a' : {'A' : 1, 'B' : 1},

         'b' : {'B' : 1},

         'c' : {'A' : 1, 'B' : 1, 'C':1},

         'd' : {'B' : 1, 'C' : 1}}

 

    PersonalRank(G, alpha, 'b', 50)     #从'b'节点开始游走

输出：

use time 0.00107598304749
B:0.312, b:0.262, c:0.115, a:0.089, d:0.089, A:0.066, C:0.066,

注意这里有一个现象，上述代码从节点b开始游走，最终算得的PR重要度最高的不是b，而是B。

personalrank经过多次的迭代游走，使得各节点的重要度趋于稳定，实际上我们根据状态转移矩阵，经过一次矩阵运算就可以直接得到系统的稳态。公式(1)的矩阵表示形式为：

(2)r=(1−α)r0+αMT

其中r是个n维向量，每个元素代表一个节点的PR重要度，r0也是个n维向量，第i个位置上是1，其余元素均为0，我们就是要为第i个节点进行推荐。M是n阶转移矩阵：

(3)Mij={1|out(i)|    if(j∈out(i))0    else

按照矩阵乘法把(2)展开就可以得到(1)。

由(2)可以得到两种变形：

(4)(I−αMT)r=(1−α)r0

(5)r=(I−αMT)−1(1−α)r0

 利用(4)，解一次线性方程组就查以得到r，对r中的各元素降序排列取出前K个就得到了节点i的推荐列表。实践中系数矩阵(I−αMT)是一个高度稀疏的矩阵，解这种线性方程流行的方法是GMRES。另外在scipy中提供了多种稀疏矩阵的存储方法：coo,lil,dia,dok,csr,csc等，各有各的优缺点，dok可以快速的按下标访问元素，csr和csc适合做矩阵的加法、乘法运算，lil省内存且按下标访问元素也很快。更多内容参见scipy中的稀疏矩阵。

由于我们只关心r中各元素的相对大小，并不关心其真实值，所以(5)可以写为

(6)r=(I−αMT)−1r0

矩阵(I−αMT)−1乘以r0相当于是取出矩阵的第i列，因此为节点i进行推荐时对矩阵(I−αMT)−1的第i列排序即可，所以求出矩阵(I−αMT)的逆就得到了所有节点的推荐结果。

pr_matrix.py

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

# coding=utf-8

__author__ = 'orisun'

 

import numpy as np

from numpy.linalg import solve

import time

from scipy.sparse.linalg import gmres, lgmres

from scipy.sparse import csr_matrix

 

if __name__ == '__main__':

    alpha = 0.8

    vertex = ['A', 'B', 'C', 'a', 'b', 'c', 'd']

    M = np.matrix([[0,        0,        0,        0.5,      0,        0.5,      0],

                   [0,        0,        0,        0.25,     0.25,     0.25,     0.25],

                   [0,        0,        0,        0,        0,        0.5,      0.5],

                   [0.5,      0.5,      0,        0,        0,        0,        0],

                   [0,        1.0,      0,        0,        0,        0,        0],

                   [0.333,    0.333,    0.333,    0,        0,        0,        0],

                   [0,        0.5,      0.5,      0,        0,        0,        0]])

    # print np.eye(n) - alpha * M.T

    r0 = np.matrix([[0], [0], [0], [0], [1], [0], [0]])  # 从'b'节点开始游走

    n = M.shape[0]

 

    # 直接解线性方程法

    A = np.eye(n) - alpha * M.T

    b = (1 - alpha) * r0

    begin = time.time()

    r = solve(A, b)

    end = time.time()

    print 'use time', end - begin

    rank = {}

    for j in xrange(n):

        rank[vertex[j]] = r[j]

    li = sorted(rank.items(), cmp=lambda x, y: cmp(x[1], y[1]), reverse=True)

    for ele in li:

        print "%s:%.3f, \t" % (ele[0], ele[1]),

    print

 

    # 采用CSR法对稀疏矩阵进行压缩存储，然后解线性方程

    A =np.eye(n) - alpha * M.T

    b = (1 - alpha) * r0

    data = list()

    row_ind = list()

    col_ind = list()

    for row in xrange(n):

        for col in xrange(n):

            if(A[row, col] != 0):

                data.append(A[row, col])

                row_ind.append(row)

                col_ind.append(col)

    AA = csr_matrix((data, (row_ind, col_ind)), shape=(n, n))

    begin = time.time()

    # 系数矩阵很稀疏时采用gmres方法求解。解方程的速度比上面快了许多

    r = gmres(AA, b, tol=1e-08, maxiter=1)[0]

    # r = lgmres(AA, (1 - alpha) * r0, tol=1e-08,maxiter=1)[0]  #lgmres用来克服gmres有时候不收敛的问题，会在更少的迭代次数内收敛

    end = time.time()

    print 'use time', end - begin

    rank = {}

    for j in xrange(n):

        rank[vertex[j]] = r[j]

    li = sorted(rank.items(), cmp=lambda x, y: cmp(x[1], y[1]), reverse=True)

    for ele in li:

        print "%s:%.3f, \t" % (ele[0], ele[1]),

    print

 

    # 求逆矩阵法。跟gmres解方程的速度相当

    A = np.eye(n) - alpha * M.T

    b = (1 - alpha) * r0

    begin = time.time()

    r = A.I * b

    end = time.time()

    print 'use time', end - begin

    rank = {}

    for j in xrange(n):

        rank[vertex[j]] = r[j, 0]

    li = sorted(rank.items(), cmp=lambda x, y: cmp(x[1], y[1]), reverse=True)

    for ele in li:

        print "%s:%.3f, \t" % (ele[0], ele[1]),

    print

 

    # 实际上可以一次性计算出从任意节点开始游走的PersonalRank结果。从总体上看，这种方法是最快的

    A = np.eye(n) - alpha * M.T

    begin = time.time()

    D = A.I

    end = time.time()

    print 'use time', end - begin

    for j in xrange(n):

        print vertex[j] + "\t",

        score = {}

        total = 0.0  # 用于归一化

        for i in xrange(n):

            score[vertex[i]] = D[i, j]

            total += D[i, j]

        li = sorted(score.items(), cmp=lambda x,

                    y: cmp(x[1], y[1]), reverse=True)

        for ele in li:

            print "%s:%.3f, \t" % (ele[0], ele[1] / total),

        print

输出：

use time 7.60555267334e-05
B:0.312, b:0.262, c:0.115, d:0.089, a:0.089, C:0.066, A:0.066, 
use time 0.000385999679565
B:0.312, b:0.262, c:0.115, a:0.089, d:0.089, A:0.066, C:0.066, 
use time 0.000133991241455
B:0.312, b:0.262, c:0.115, d:0.089, a:0.089, C:0.066, A:0.066, 
use time 0.000274181365967
A A:0.314, c:0.189, B:0.166, a:0.159, C:0.076, d:0.063, b:0.033, 
B B:0.390, c:0.144, d:0.111, a:0.111, C:0.083, A:0.083, b:0.078, 
C C:0.314, c:0.189, B:0.166, d:0.159, A:0.076, a:0.063, b:0.033, 
a a:0.308, B:0.222, A:0.159, c:0.133, d:0.070, C:0.063, b:0.044, 
b B:0.312, b:0.262, c:0.115, d:0.089, a:0.089, C:0.066, A:0.066, 
c c:0.340, B:0.192, C:0.126, A:0.126, d:0.089, a:0.089, b:0.038, 
d d:0.308, B:0.222, C:0.159, c:0.133, a:0.070, A:0.063, b:0.044,