上面的二部圖表示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的節點開始遊走，往前走一步。繼續遊走的概率是$\alpha$，停留在當前節點的概率是$1-\alpha$。

第一次遊走， 從A節點以各自50%的概率走到了a和c，這樣a和c就分得了A的部分重要度，$PR(a)=PR(c)=\alpha*PR(A)*0.5$。最後$PR(A)$變為$1-\alpha$。第一次遊走結束後PR不為0的節點有A a c。

第二次遊走，分別從節點A a c開始，往前走一步。這樣節點a分得A $\frac{1}{2}*\alpha$的重要度，節點c分得A $\frac{1}{2}*\alpha$的重要度，節點A分得a $\frac{1}{2}*\alpha$的重要度，節點A分得c $\frac{1}{3}*\alpha$的重要度，節點B分得a $\frac{1}{2}*\alpha$的重要度，節點B分得c $\frac{1}{3}*\alpha$的重要度，節點C分得c $\frac{1}{3}*\alpha$的重要度。最後$PR(A)$要加上$1-\alpha$

經過以上推演，可以概括成公式：

\begin{equation}PR(j)=\left\{\begin{matrix}\alpha*\sum_{i\in{in(j)}}\frac{PR(i)}{|out(i)|}\ \ \ \ if(j\ne{u})\\(1-\alpha)+\alpha*\sum_{i\in{in(j)}}\frac{PR(i)}{|out(i)|} \ \ \ \ if(j=u)\end{matrix}\right.\label{pr}\end{equation}

$u$是待推薦的使用者節點。

personalrank.py

#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經過多次的迭代遊走，使得各節點的重要度趨於穩定，實際上我們根據狀態轉移矩陣，經過一次矩陣運算就可以直接得到系統的穩態。公式\eqref{pr}的矩陣表示形式為：

\begin{equation}r=(1-\alpha)r_0+\alpha{M^T}\label{mat}\end{equation}

其中$r$是個n維向量，每個元素代表一個節點的PR重要度，$r_0$也是個n維向量，第i個位置上是1，其餘元素均為0，我們就是要為第i個節點進行推薦。$M$是n階轉移矩陣：

\begin{equation}M_{ij}=\left\{\begin{matrix}\frac{1}{|out(i)|}\ \ \ \ if(j\in{out(i)})\\0 \ \ \ \ else \end{matrix}\right.\end{equation}

按照矩陣乘法把\eqref{mat}展開就可以得到\eqref{pr}。

由\eqref{mat}可以得到兩種變形：

\begin{equation}(I-\alpha{M^T})r=(1-\alpha)r_0\label{eq}\end{equation}

\begin{equation}r=(I-\alpha{M^T})^{-1}(1-\alpha)r_0\label{inv}\end{equation}

 利用\eqref{eq}，解一次線性方程組就查以得到$r$，對$r$中的各元素降序排列取出前K個就得到了節點i的推薦列表。實踐中係數矩陣$(I-\alpha{M^T})$是一個高度稀疏的矩陣，解這種線性方程流行的方法是GMRES。另外在scipy中提供了多種稀疏矩陣的儲存方法：coo,lil,dia,dok,csr,csc等，各有各的優缺點，dok可以快速的按下標訪問元素，csr和csc適合做矩陣的加法、乘法運算，lil省記憶體且按下標訪問元素也很快。更多內容參見scipy中的稀疏矩陣。

由於我們只關心$r$中各元素的相對大小，並不關心其真實值，所以\eqref{inv}可以寫為

\begin{equation}r=(I-\alpha{M^T})^{-1}r_0\end{equation}

矩陣$(I-\alpha{M^T})^{-1}$乘以$r_0$相當於是取出矩陣的第i列，因此為節點i進行推薦時對矩陣$(I-\alpha{M^T})^{-1}$的第i列排序即可，所以求出矩陣$(I-\alpha{M^T})$的逆就得到了所有節點的推薦結果。

pr_matrix.py

# 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,