Dynamical SimRank Search on Time-Varying Networks

In this article, we study the efficient dynamical computation of all-pairs SimRanks on time-varying graphs. Li {\em et al}.'s approach requires $O(r^4 n^2)$ time and $O(r^2 n^2)$ memory in a graph with $n$ nodes, where $r$ is the target rank of the low-rank SVD. (1) We first consider edge update that does not accompany new node insertions. We show that the SimRank update $\Delta S$ in response to every link update is expressible as a rank-one Sylvester matrix equation. This provides an incremental method requiring $O(Kn^2)$ time and $O(n^2)$ memory in the worst case to update all pairs of similarities for $K$ iterations. (2) To speed up the computation further, we propose a lossless pruning strategy that captures the"affected areas"of $\Delta S$ to eliminate unnecessary retrieval. This reduces the time of the incremental SimRank to $O(K(m+|AFF|))$, where $m$ is the number of edges in the old graph, and $|AFF| (\le n^2)$ is the size of"affected areas"in $\Delta S$, and in practice, $|AFF| \ll n^2$. (3) We also consider edge updates that accompany node insertions, and categorize them into three cases, according to which end of the inserted edge is a new node. For each case, we devise an efficient incremental algorithm that can support new node insertions. (4) We next design an efficient batch incremental method that handles"similar sink edges"simultaneously and eliminates redundant edge updates. (5) To achieve linear memory, we devise a memory-efficient strategy that dynamically updates all pairs of SimRanks column by column in just $O(Kn+m)$ memory, without the need to store all $(n^2)$ pairs of old SimRank scores. Experimental studies on various datasets demonstrate that our solution substantially outperforms the existing incremental SimRank methods, and is faster and more memory-efficient than its competitors on million-scale graphs.


Introduction
Recent rapid advances in web data management reveal that link analysis is becoming an important tool for similarity assessment. Due to the growing number of applications in e.g., social networks, recommender systems, citation analysis, and link prediction [9], a surge of graph-based similarity measures have surfaced over the past decade. For instance, Brin and Page [2] proposed a very successful relevance measure, called Google PageRank, to rank web pages. Jeh and Widom [9] devised SimRank, an appealing pair-wise similarity measure that quantifies the structural equivalence of two nodes based on link structure. Recently, Sun et al. [21] invented PathSim to retrieve nodes proximities in a heterogeneous graph. Among these emerging link based measures, SimRank has stood out as an attractive one in recent years, due to its simple and iterative philosophy that "two nodes are similar if they are pointed to by similar nodes", coupled with the base case that "every node is most similar to itself". This recursion not only allows SimRank to capture the global structure of a graph, but also equips SimRank with mathematical insights that attract many researchers. For example, Fogaras and Rácz [5] interpreted SimRank as the meeting time of the coalescing pair-wise random walks. Li et al. [13] harnessed an elegant matrix equation to formulate the closed form of SimRank.
Nevertheless, the batch computation of SimRank is costly: O(Kd ′ n 2 ) time for all node-pairs [24], where K is the total number of iterations, and d ′ ≤ d (d is the av-erage in-degree of a graph). Generally, many real graphs are large, with links constantly evolving with minor changes. This is especially apparent in e.g., co-citation networks, web graphs, and social networks. As a statistical example [17], there are 5%-10% links updated every week in a web graph. It is rather expensive to recompute similarities for all pairs of nodes from scratch when a graph is updated. Fortunately, we observe that when link updates are small, the affected areas for SimRank updates are often small as well. With this comes the need for incremental algorithms that compute changes to SimRank in response to link updates, to discard unnecessary recomputations. In this article, we investigate the following problem for SimRank evaluation: Problem (Incremental SimRank Computation) Given an old digraph G, old similarities in G, link changes ∆G 1 to G, and a damping factor C ∈ (0, 1). Retrieve the changes to the old similarities.
Our research for the above SimRank problem is motivated by the following real applications: need iteratively compute the SimRank matrix S of size (n×n) in a centralised way (by using a single machine). In contrast, our incremental approach can significantly improve the computational efficiency of all pairs of Sim-Ranks by retrieving S in a decentralised way as follows: We first employ a graph partitioning algorithm (e.g., METIS 2 ) that can decompose the large graph G into several small blocks such that the number of the edges with endpoints in different blocks is minimised. In this example, we partition G into 3 blocks, G 1 ∪ G 2 ∪ G 3 , along with 2 edges {(f, c), (f, k)} across the blocks, as depicted in the first row of Figure 1.
Let G old := G 1 ∪ G 2 ∪ G 3 and ∆G := {(f, c), (f, k)}. Then, G can be viewed as "G old perturbed by ∆G edge insertions". That is, Based on this decomposition, we can efficiently compute S over G by dividing S into two parts: where S old is obtained by using a batch SimRank algorithm over G old , and ∆S is derived from our proposed incremental method which takes S old and ∆G as input.
It is worth mentioning that this way of retrieving S is far more efficient than directly computing S over G via a batch algorithm. There are two reasons: 2 http://glaros.dtc.umn.edu/gkhome/views/metis Firstly, S old can be efficiently computed in a decentralised way. It is a block diagonal matrix with no need of n × n space to store S old . This is because G old is only comprised of several connected components (G 1 , G 2 , G 3 ), and there are no edges across distinct components. Thus, S old exhibits a block diagonal structure: To obtain S old , instead of applying the batch SimRank algorithm over the entire G old , we can apply the batch SimRank algorithm over each component G i (i = 1, 2, 3) independently to obtain the i-th diagonal block of S old , S Gi . Indeed, each S Gi is computable in parallel. Even if S old is computed using a single machine, only O(n 2 1 + n 2 2 + n 2 3 ) space is required to store its diagonal blocks, where n i is the number of nodes in each G i , rather than O(n 2 ) space to store the entire S old (see Figure 1).
Secondly, after graph partitioning, there are not many edges across components. Small size of ∆G often leads to sparseness of ∆S in general. Hence, ∆S is stored in a sparse format. In addition, our incremental SimRank method will greatly accelerate the computation of ∆S.
Hence, along with graph partitioning, our incremental SimRank research will significantly enhance the computational efficiency of SimRank on large graphs, using a decentralised fashion. ⊓ ⊔ Our research on the incremental SimRank problem not only can decentralise large-scale SimRank retrieval, but also will enable a substantial speedup on the batch computation of SimRank, as indicated below.

Example 2 (Accelerate Batch Computation of SimRank)
Consider the citation network G in Figure 2, where each node denotes a paper, and an edge a reference from one paper to another. We wish to assess all pairs of similarities between papers. Unlike existing batch computation that iteratively retrieves all-pairs SimRanks over the entire G, our incremental method significantly accelerates the batch computation of SimRank as follows: We first utilise BFS or DFS search to find a spanning tree (or arborescence) of G, denoted as G old . We observe that, due to the tree structure, all-pairs SimRank scores in G old are relatively easier to compute. For example, each entry of Li et al.'s SimRank matrix S old over G old can be obtained from a lightweight formula: if a and b are not on the same level of the tree G old ; where C is a damping factor, λ a,b is the number of edges from the lowest common ancestor of (a, b) to node a (or equivalently, to b) in the tree G old , and H is the number of edges from the root to node a (or equivalently, to b) in the tree G old .
Given S old , we next apply our incremental SimRank method that can significantly speed up the computation of new SimRank scores S in G. Specifically, we denote by ∆G the set of edges in G but not in G old . In Figure 2, Based on S old and ∆G, our incremental SimRank method can dynamically retrieve only the changes to S old in response to ∆G, whose result ∆S is a sparse matrix, as illustrated in Figure 2.
It is important to note that it does not require n× n memory to store S old because G old is a tree structure. If there are H levels in the tree G old with n l nodes on level l (l = 1, · · · , H), then we only need the space to store the nonzero diagonal blocks of S old . ⊓ ⊔ These examples show that our incremental SimRank is useful to (i) decentralise large-scale SimRank retrieval, and (ii) accelerate the batch computation of SimRank. Despite its usefulness, existing work on incremental Sim-Rank computation is rather limited. To the best of our knowledge, there is a relative paucity of work [10,13,20,25] on incremental SimRank problems. Shao et al. [20] proposed a novel two-stage random-walk sampling scheme, named TSF, which can support top-k SimRank search over dynamic graphs. In the preprocessing stage, TSF samples R g one-way graphs that serve as an index for querying process. At query stage, for each one-way graph, R q new random walks of node u are sampled. However, the dynamic SimRank problems studied in [20] and this work are different: This work focuses on all (n 2 ) pairs of SimRank retrieval, which requires O(K(m + |AFF|)) time to compute the entire matrix S in a deterministic style. In Section 7, we have proposed a memoryefficient version of our incremental method that updates all pairs of similarities in a column-by-column fashion within only O(Kn + m) memory. In comparison, [20] focuses on top-k SimRank dynamic search w.r.t. a given query u. It incrementally retrieves only k (≤ n) nodes with highest SimRank scores in a single column S ⋆,u , which requires O(K 2 R q R g ) average query time 3 to retrieve S ⋆,u along with O(n log k) time to return top-k results from S ⋆,u . Recently, Jiang et al. [10] pointed out that the probabilistic error guarantee of Shao et al.'s method is based on the assumption that no cycle in the given graph has a length shorter than K, and they proposed READS, a new efficient indexing scheme that improves precision and indexing space for dynamic SimRank search.  [13] on all-pairs dynamic search is exactly the same as ours: the goal is to retrieve changes ∆S to all-pairs SimRank scores S, given old graph G, link changes ∆G to G. To address this problem, the central idea of [13] is to factorize the backward transition matrix Q of the original graph into U · Σ · V T via a singular value decomposition (SVD) first, and then incrementally estimate the updated matrices of U, Σ, V T for link changes at the expense of exactness. Consequently, updating all pairs of similarities entails O(r 4 n 2 ) time and O(r 2 n 2 ) memory yet without guaranteed accuracy, where r (≤ n) is the target rank of the low-rank SVD approximation 4 . This method is efficient to graphs when r is extremely small, e.g., a star graph (r = 1). However, in general, r is not always negligibly small.
(Please refer to Appendix A for a discussion in detail, and Appendix C.1 for an example.)

Main Contributions
Motivated by this, we propose an efficient and accurate scheme for incrementally computing all-pairs SimRanks on link-evolving graphs. Our main contributions consist of the following five ingredients: -To speed up the computation further, we also propose an effective pruning strategy that captures the "affected areas" of ∆S to discard unnecessary retrieval, without loss of accuracy. This reduces the time of incremental SimRank to O(K(m + |AFF|)), where |AFF| (≤ n 2 ) is the size of "affected areas" in ∆S, and in practice, |AFF| ≪ n 2 . (Section 4) -We also consider edge updates that accompany new node insertions, and distinguish them into three categories, according to which end of the inserted edge is a new node. For each case, we devise an efficient incremental SimRank algorithm that can support new nodes insertion and accurately update affected SimRank scores. (Section 5) 4 According to [13], using our notation, r ≤ rank(Σ + U T · ∆Q · V), where ∆Q is the changes to Q for link updates.
-We next investigate the batch updates of dynamical SimRank computation. Instead of dealing with each edge update one by one, we devise an efficient algorithm that can handle a sequence of edge insertions and deletions simultaneously, by merging "similar sink edges" and minimizing unnecessary updates. (Section 6) -To achieve linear memory efficiency, we also express ∆S as the sum of many rank-one tensor products, and devise a memory-efficient technique that updates all-pairs SimRanks in a column-by-column style in O(Kn + m) memory, without loss of exactness. (Section 7) -We conduct extensive experiments on real and synthetic datasets to demonstrate that our algorithm (a) is consistently faster than the existing incremental methods from several times to over one order of magnitude; (b) is faster than its batch counterparts especially when link updates are small; (c) for batch updates, runs faster than the repeated unit update algorithms; (d) entails linear memory and scales well on billion-edge graphs for all-pairs SimRank update; (e) is faster than L-TSF and its memory space is less than L-TSF; (f) entails more time on Cases (C0) and (C2) than Cases (C1) and (C3) for four edge types, and Case (C3) runs the fastest. (Section 8) This article is a substantial extension of our previous work [25]. We have made the following new updates: (1) In Section 5, we study three types of edge updates that accompany new node insertions. This solidly extends [25] and Li et al. 's incremental method [13] whose edge updates disallow node changes. (2) In Section 6, we also investigate batch updates for dynamic SimRank computation, and devise an efficient algorithm that can handle "similar sink edges" simultaneously and discard unnecessary unit updates further. (3) In Section 7, we propose a memory-efficient strategy that significantly reduces the memory from O(n 2 ) to O(Kn + m) for incrementally updating all pairs of SimRanks on millionscale graphs, without compromising running time and accuracy. (4) In Section 8, we conduct additional experiments on real and synthetic datasets to verify the high scalability and fast computational time of our memoryefficient methods, as compared with the L-TSF method. (5) In Section 9, we update the related work section by incorporating state-of-the-art SimRank research.

SimRank Background
In this section, we give a broad overview of SimRank. Intuitively, the central theme behind SimRank is that "two nodes are considered as similar if their incoming neighbors are themselves similar". Based on this idea, there have emerged two widely-used SimRank models: (1) Li et al.'s model (e.g., [6,8,13,18,27]) and (2) Jeh and Widom's model (e.g., [4,9,11,16,20]). Throughout this article, our focus is on Li et al. 's SimRank model, also known as Co-SimRank in [18], since the recent work [18] by Rothe and Schütze has showed that Co-SimRank is more accurate than Jeh and Widom's SimRank model in real applications such as bilingual lexicon extraction. (Please refer to Remark 1 for detailed explanations.)

Li et al.'s SimRank model
Given a directed graph G = (V, E) with a node set V and an edge set E, let Q be its backward transition matrix (that is, the transpose of the column-normalized adjacency matrix), whose entry [Q] i,j = 1/in-degree(i) if there is an edge from j to i, and 0 otherwise. Then, Li et al.'s SimRank matrix, denoted by S, is defined as where C ∈ (0, 1) is a damping factor, which is generally taken to be 0.6-0.8, and I n is an n × n identity matrix (n = |V |). The notation (⋆) T is the matrix transpose.
Recently, Rothe and Schütze [18] have introduced Co-SimRank, whose definition is Comparing Eqs. (1) and (2), we can readily verify that Li et al.'s SimRank scores equal Co-SimRank scores scaled by a constant factor (1 − C), i.e., S = (1 − C) ·S. Hence, the relative order of all Co-SimRank scores inS is exactly the same as that of Li et al.'s SimRank scores in S even though the entries inS can be larger than 1. That is, the ranking of Co-SimRankS( * , * ) is identical to the ranking of Li et al.'s SimRank S( * , * ).

Remark 1
The recent work by Kusumoto et al. [11] has showed that S and S ′ do not produce the same results.
Recently, Yu and McCann [27] have showed the subtle

old/new SimRank matrix
In n × n identity matrix  [13] for updating SimRank does not always obtain the correct solution S to Eq.(1). (Please refer to Appendix A for theoretical explanations). Table 1 lists the notations used in this article.

Edge Update without node insertions
In this section, we consider edge update that does not accompany new node insertions, i.e., the insertion of new edge (i, j) into G = (V, E) with i ∈ V and j ∈ V . In this case, the new SimRank matrixS and the old one S are of the same size. As such, it makes sense to denote the SimRank change ∆S asS − S. Below we first introduce the big picture of our main idea, and then present rigorous justifications and proofs.

The main idea
For each edge (i, j) insertion, we can show that ∆Q is a rank-one matrix, i.e., there exist two column vectors u, v ∈ R n×1 such that ∆Q ∈ R n×n can be decomposed into the outer product of u and v as follows: Based on Eq.(4), we then have an opportunity to efficiently compute ∆S, by characterizing it as where the auxiliary matrix M ∈ R n×n satisfies the following rank-one Sylvester equation: Here, u, w are two obtainable column vectors: u can be derived from Eq.(4), and w can be described by the old Q and S (we will provide their exact expressions later after some discussions); andQ = Q + ∆Q. Thus, computing ∆S boils down to solving M in Eq.(6). The main advantage of solving M via Eq.(6), as compared to directly computing the new scoresS via SimRank formulã is the high computational efficiency. More specifically, solvingS via Eq.(7) needs expensive matrix-matrix multiplications, whereas solving M via Eq.(6) involves only matrix-vector and vector-vector multiplications, which is a substantial improvement achieved by our observation that (C · uw T ) ∈ R n×n in Eq.(6) is a rank-one matrix, as opposed to the (full) rank-n matrix (1 − C) · I n in Eq. (7). To further elaborate on this, we can readily convert the recursive forms of Eqs.(6) and (7), respectively, into the series forms: To compute the sums in Eq.(8) for M, a conventional way is to memoize M 0 ← C · u · w T first (where the intermediate result M 0 is an n × n matrix), and then iterate as which involves expensive matrix-matrix multiplications (e.g.,Q · M k ). In contrast, our method takes advantage of the rank-one structure of u·w T to compute the sums in Eq.(8) for M, by converting the conventional matrixmatrix multiplicationsQ · (uw T ) ·Q T into only matrixvector and vector-vector multiplications (Qu) · (Qw) T . To be specific, we leverage two vectors ξ k , η k , and iteratively compute Eq.(8) as so that matrix-matrix multiplications are safely avoided.
3.2 Describing u, v, w in Eqs. (4) and (6) To obtain u and v in Eq.(4) at a low cost, we have the following theorem.
Theorem 1 Given an old digraph G = (V, E), if there is a new edge (i, j) with i ∈ V and j ∈ V to be added to G, then the change to Q is an n × n rank-one matrix, i.e., (Please refer to Appendix B.1 for the proof of Theorem 1, and Appendix C.2 for an example.) Theorem 1 suggests that the change ∆Q is an n × n rank-one matrix, which can be obtain in only constant time from d j and [Q] T j,⋆ . In light of this, we next describe w in Eq. (6) in terms of the old Q and S such that Eq. (6) is a rank-one Sylvester equation.
Theorem 2 Let (i, j) i∈V, j∈V be a new edge to be added to G (resp. an existing edge to be deleted from G). Let u and v be the rank-one decomposition of ∆Q = u·v T . Then, (i) there exists a vector w = y + λ such that Eq.(6) is the rank-one Sylvester equation.
(ii) Utilizing the solution M to Eq.(6), the SimRank update matrix ∆S can be represented by Eq. (5). ⊓ ⊔ (The proof of Theorem 2 is in Appendix B.2.) Theorem 2 provides an elegant expression of w in Eq. (6). To be precise, given Q and S in the old graph G, and an edge (i, j) inserted to G, one can find u and v via Theorem 1 first, and then resort to Theorem 2 to compute w from u, v, Q, S. Due to the existence of the vector w, it can be guaranteed that the Sylvester equation (6) is rank-one. Henceforth, our aforementioned method can be employed to iteratively compute M in Eq.(8), requiring no matrix-matrix multiplications.

Characterizing ∆S
Leveraging Theorems 1 and 2, we next characterize the SimRank change ∆S. Theorem 3 If there is a new edge (i, j) with i ∈ V and j ∈ V to be inserted to G, then the SimRank change ∆S can be characterized as where the auxiliary vector γ is obtained as follows: (ii) when d j > 0, (15) and the scalar λ can be derived from Theorem 3 provides an efficient method to compute the incremental SimRank matrix ∆S, by utilizing the previous information of Q and S, as opposed to [13] that requires to maintain the incremental SVD.

Deleting an edge
For an edge deletion, we next propose a Theorem 3-like technique that can efficiently update SimRanks.
Theorem 4 When an edge (i, j) i∈V, j∈V is deleted from G = (V, E), the changes to Q is a rank-one matrix, which can be described as The changes ∆S to SimRank can be characterized as where the auxiliary vector γ :=

Inc-uSR Algorithm
We present our efficient incremental approach, denoted as Inc-uSR (in Appendix D.1), that supports the edge insertion without accompanying new node insertions. The complexity of Inc-uSR is bounded by O(Kn 2 ) time and O(n 2 ) memory 5 in the worst case for updating all n 2 pairs of similarities.
(Please refer to Appendix D.1 for a detailed description of Inc-uSR, and Appendix C.3 for an example.)

Pruning Unnecessary Node-Pairs in ∆S
After the SimRank update matrix ∆S has been characterized as a rank-one Sylvester equation, pruning techniques can further skip node-pairs with unchanged Sim-Ranks in ∆S (called "unaffected areas").

Affected Areas in ∆S
We next reinterpret the series M in Theorem 3, aiming to identify "affected areas" in ∆S. Due to space limitations, we mainly focus on the edge insertion case of d j > 0. Other cases have the similar results.
By substituting Eq.(15) back into Eq.(13), we can readily split the series form of M into three parts: tallies the weighted sum of the following new paths for node-pair (a, b): Such paths are the concatenation of four types of sub-paths (as depicted above) associated with four ma- the inserted edge j ⇐ i. When such entire concatenated paths exist in the new graph, they should be accommodated for assessing the new SimRank [S] a,b in response to the edge insertion (i, j) because our reinterpretation of SimRank indicates that SimRank counts all the symmetric in-link paths, and the entire concatenated paths can prove to be symmetric in-link paths. Likewise, Parts 2 and 3 of [M] a,b , respectively, tally the weighted sum of the following paths for pair (a, b): Indeed, when edge (i, j) is inserted, only these three kinds of paths have extra contributions for M (therefore for ∆S). As incremental updates in SimRank merely tally these paths, node-pairs without having such paths could be safely pruned. In other words, for those pruned node-pairs, the three kinds of paths will have "zero contributions" to the changes in M in response to edge insertion. Thus, after pruning, the remaining node-pairs in G constitute the "affected areas" of M.
We next identify "affected areas" of M, by pruning redundant node-pairs in G, based on the following.
Then, for every iteration k = 0, 1, · · · , the matrix M k has the following sparse property: For the edge (i, j) deletion case, all the above results hold except that, in Eq.(21), the conditions d j = 0 and d j > 0 are, respectively, replaced by d j = 1 and d j > 1. ⊓ ⊔ (Please refer to Appendix B.5 for the proof and intuition of Theorem 5, and Appendix C.4 for an example.) Theorem 5 provides a pruning strategy to iteratively eliminate node-pairs with a-priori zero values in M k (thus in ∆S). Hence, by Theorem 5, when edge (i, j) is updated, we just need to consider node-pairs in

Inc-SR Algorithm with Pruning
Based on Theorem 5, we provide a complete incremental algorithm, referred to as Inc-SR, by incorporating our pruning strategy into Inc-uSR. The total time of (22), being the average size of "affected areas" in M k for K iterations.
(Please refer to Appendix D.2 for Inc-SR algorithm description and its complexity analysis.)

Edge Update with node insertions
In this section, we focus on the edge update that accompanies new node insertions. Specifically, given a new edge (i, j) to be inserted into the old graph G = (V, E), we consider the following cases when For each case, we devise an efficient incremental algorithm that can support new node insertions and can accurately update only "affected areas" of SimRanks.
Remark 2 Let n = |V |, without loss of generality, it can be tacitly assumed that a) in case (C1), new node j / ∈ V is indexed by (n + 1); b) in case (C2), new node i / ∈ V is indexed by (n + 1); c) in case (C3), new nodes i / ∈ V and j / ∈ V are indexed by (n + 1) and (n + 2), respectively.

Inserting an edge
In this case, the inserted new edge (i, j) accompanies the insertion of a new node j. Thus, the size of the new SimRank matrixS is different from that of the old S. As a result, we cannot simply evaluate the changes to S by adoptingS − S as we did in Section 3.
To resolve this problem, we introduce the block matrix representation of new matrices for edge insertion. Firstly, when a new edge (i, j) i∈V,j / ∈V is inserted to G, the new transition matrixQ can be described as Intuitively,Q is formed by bordering the old Q by 0s except [Q] j,i = 1. Utilizing this block structure ofQ, we can obtain the new SimRank matrix, which exhibits a similar block structure, as shown below: where S ∈ R n×n is the old SimRank matrix of G. ⊓ ⊔ Proof We substitute the newQ in Eq.(23) back into the SimRank equationS = C ·Q ·S ·Q T + (1 − C) · I n+1 : By expanding the right-hand side, we can obtain The above block matrix equation implies that Due to the uniqueness of S in Eq.(1), it follows that Thus, we havẽ Combining all blocks ofS together yields Eq. (24). ⊓ ⊔ Theorem 6 provides an efficient incremental way of computing the new SimRank matrixS for unit insertion of the case (C1). Precisely, the newS is formed by bordering the old S by the auxiliary vector y. To obtain y (and therebyS), we just need use the i-th column of S with one matrix-vector multiplication (Q[S] ⋆,i ). Thus, the total cost of computing newS requires O(m) time, as illustrated in Algorithm 1. Fig. 4. If the new edge (i, p) with new node p is inserted to G, the newS can be updated from the old S as follows:

Example 3 Consider the citation digraph G in
According to Theorem 6, since C = 0.8 and We now focus on the case (C2), the insertion of an edge (i, j) with i / ∈ V and j ∈ V . Similar to the case (C1), the new edge accompanies the insertion of a new node i. Hence,S − S makes no sense.
However, in this case, the dynamic computation of SimRank is far more complicated than that of the case (C1), in that such an edge insertion not only increases the dimension of the old transition matrix Q by one, but also changes several original elements of Q, which may recursively influence SimRank similarities. Specifically, the following theorem shows, in the case (C2), how Q changes with the insertion of an edge (i, j) i / ∈V,j∈V .
Theorem 7 Given an old digraph G = (V, E), if there is a new edge (i, j) with i / ∈ V and j ∈ V to be added to G, then the new transition matrix can be expressed as where Q is the old transition matrix of G. ⊓ ⊔ Proof When edge (i, j) with i / ∈ V and j ∈ V is added, there will be two changes to the old Q: (ii) The size of the old Q is added by 1, with new entry [Q] j,i = 1 dj+1 in the bordered areas and 0s elsewhere: Combining Eqs. (26) and (27) yields (25). ⊓ ⊔ Theorem 7 exhibits a special structure of the newQ: it is formed by borderingQ by 0s except [Q] j,i = 1 dj+1 , whereQ is a rank-one update of the old Q. The block structure ofQ inspires us to partition the new SimRank matrixS conformably into the similar block structure: To determine each block ofS with respect to the old S, we next present the following theorem. (28) such that the new SimRank matrixS is expressible as where S is the old SimRank of G, and ∆S 11 satisfies the rank-two Sylvester equation: withQ being defined by Theorem 7. ⊓ ⊔ Proof We plugQ of Eq.(25) into the SimRank formula: In 0 0 1 By using block matrix multiplications, the above equation can be simplified as Block-wise comparison of both sides of Eq.(31) yields Combing the above equations with Eq.(32) produces ApplyingS 11 = S+∆S 11 and S = CQSQ T +(1−C)I n to Eq.(33) and rearranging the terms, we have with α and y being defined by Eq. (28). ⊓ ⊔ Theorem 8 implies that, in the case (C2), after a new edge (i, j) is inserted, the new SimRank matrixS takes an elegant diagonal block structure: the upper-left block ofS is perturbed by ∆S 11 which is the solution to the rank-two Sylvester equation (30); the lower-right block ofS is a constant (1 − C). This structure ofS suggests that the inserted edge (i, j) i / ∈V,j∈V only has a recursive impact on the SimRanks with pairs (x, y) ∈ V × V , but with no impacts on pairs ( Thus, our incremental way of computing the newS will focus on the efficiency of obtaining ∆S 11 from Eq.(30). Fortunately, we notice that ∆S 11 satisfies the rank-two Sylvester equation, whose algebraic structure is similar to that of ∆S in Eqs. (5) and (6) (in Section 3). Hence, our previous techniques to compute ∆S in Eqs. (5) and (6) can be analogously applied to compute ∆S 11 in Eq.(30), thus eliminating costly matrix-matrix multiplications, as will be illustrated in Algorithm 2.
One disadvantage of Theorem 8 is that, in order to get the auxiliary vector z for evaluatingS, one has to memorize the entire old matrix S in Eq. (28). In fact, we can utilize the technique of rearranging the terms of the SimRank Eq.(1) to characterize QS[Q] T j,⋆ in terms of only one vector [S] ⋆,j so as to avoid memoizing the entire S, as shown below.

Theorem 9
The auxiliary matrix ∆S 11 in Theorem 8 can be represented as whereQ is defined by Theorem 7 and and S is the old SimRank matrix of G. ⊓ ⊔ Proof We multiply the SimRank equation by e j to get Combining this with y = QS[Q] T j,⋆ in Eq.(28) produces Plugging these results into Eq.(28), we can get Eq.(35). Also, the recursive form of ∆S 11 in Eq.(30) can be converted into the following series: For edge insertion of the case (C2), Theorem 9 gives an efficient method to compute the update matrix ∆S 11 . We note that the form of ∆S 11 in Eq.(34) is similar to that of ∆S in Eq. (13). Thus, similar to Theorem 3, the follow method can be applied to compute M so as to avoid matrix-matrix multiplications.
In Algorithm 2, we present the edge insertion of our method for the case (C2) to incrementally update new SimRank scores. The total complexity of Algorithm 2 is O(Kn 2 ) time and O(n 2 ) memory in the worst case for retrieving all n 2 pairs of scores, which is dominated by Line 8. To reduce its computational time further, the similar pruning techniques in Section 4 can be applied to Algorithm 2. This can speed up the computational time to O(K(m + |AFF|)), where |AFF| is the size of "affected areas" in ∆S 11 .
Example 4 Consider the citation digraph G in Fig.5. If the new edge (p, j) with new node p is inserted to G, the newS can be incrementally derived from the old S as follows: First, we obtain ∆S 11 according to Theorem 9. Note that C = 0.8, d j = 2, and the old SimRank scores 0, · · · , 0, 0.2064, 0, 0, 0.3104, 0.5104, 0, · · · , 0 T It follows from Eq.(35) that the auxiliary vector Utilizing z, we can obtain M from Eq.(34). Thus, ∆S 11 can be computed from M as Next, by Theorem 8, we obtain the new SimRank Inserting an edge (i, j) with i / ∈ V and j / ∈ V We next focus on the case (C3), the insertion of an edge (i, j) with i / ∈ V and j / ∈ V . Without loss of generality, it can be tacitly assumed that nodes i and j are indexed by n+1 and n+2, respectively. In this case, the inserted edge (i, j) accompanies the insertion of two new nodes, which can form another independent component in the new graph.
In this case, the new transition matrixQ can be characterized as a block diagonal matrix With this structure, we can infer that the new SimRank matrixS takes the block diagonal form as This is because, after a new edge (i, j) i / ∈V,j / ∈V is added, all node-pairs (x, y) ∈ (V × {i, j} ∪ {i, j} × V ) have zero SimRank scores since there are no connections between nodes x and y. Besides, the inserted edge (i, j) is an independent component that has no impact on s(x, y) Hence, the submatrixŜ of the new SimRank matrix can be derived by solving the equation: This suggests that, for unit insertion of the case (C3), the new SimRank matrix becomes Algorithm 3 presents our incremental method to obtain the new SimRank matrixS for edge insertion of the case (C3), which requires just O(1) time.

Batch Updates
In this section, we consider the batch updates problem for incremental SimRank, i.e., given an old graph G = (V, E) and a sequence of edges ∆G to be updated to G, the retrieval of new SimRank scores in G ⊕ ∆G. Here, the set ∆G can be mixed with insertions and deletions: The straightforward approach to this problem is to update each edge of ∆G one by one, by running a unit update algorithm for |∆G| times. However, this would produce many unnecessary intermediate results and redundant updates that may cancel out each other.
Example 5 Consider the old citation graph G in Fig. 6, and a sequence of edge updates ∆G to G: Example 5 suggests that a portion of redundancy in ∆G arises from the insertion and deletion of the same edge that may cancel out each other. After cancellation, it is easy to verify that To obtain ∆G net from ∆G, we can readily use hashing techniques to count occurrences of updates in ∆G. More specifically, we use each edge of ∆G as a hash key, and initialize each key with zero count. Then, we scan each edge of ∆G once, and increment (resp. decrement) its count by one each time an edge insertion (resp. deletion) appears in ∆G. After all edges in ∆G are scanned, the edges whose counts are nonzeros make a net update ∆G net . All edges in ∆G net with +1 (resp. −1) counts make a net insertion update ∆G + net (resp. a net deletion update ∆G − net ). Clearly, we have Having reduced ∆G to the net edge updates ∆G net , we next merge the updates of "similar sink edges" in ∆G net to speedup the batch updates further.
We first introduce the notion of "similar sink edges".
Definition 1 Two distinct edges (a, c) and (b, c) are called "similar sink edges" w.r.t. node c if they have a common end node c that both a and b point to. ⊓ ⊔ "Similar sink edges" is introduced to partition ∆G net . To be specific, we first sort all the edges {(i p , j p )} of ∆G + net (resp. ∆G − net ) according to its end node j p . Then, the "similar sink edges" w.r.t. node j p form a partition of ∆G + net (resp. ∆G − net ). For each block {(i p k , j p )} in ∆G + net , we next split it further into two sub-blocks according to whether its end node i p k is in the old V . Thus, after partitioning, each block in ∆G + net (resp. ∆G − net ), denoted as {(i 1 , j), (i 2 , j), · · · , (i δ , j)}, falls into one of the following cases: Example 6 Let us recall ∆G net derived by Example 5, We first partition ∆G + net by "similar sink edges" into In the first block of ∆G + net , since the nodes q / ∈ V , j ∈ V , and k ∈ V , we will partition this block further into {(q, i, +)} ∪ {(j, i, +), (k, i, +)}. Eventually,   The main advantage of our partitioning approach is that, after partition, all the edge updates in each block can be processed simultaneously, instead of one by one. To elaborate on this, we use case (C0) as an example, i.e., the insertion of δ edges {(i 1 , j), (i 2 , j), · · · , (i δ , j)} into G = (V, E) when i 1 ∈ V, · · · , i δ ∈ V , and j ∈ V . Analogous to Theorem 1, one can readily prove that, after such δ edges are inserted, the changes ∆Q to the old transition matrix is still a rank-one matrix that can be decomposed asQ where e I is an n×1 vector with its entry [e I ] x = 1 if x ∈ I {i 1 , i 2 , · · · , i δ }, and [e I ] x = 0 if x / ∈ V . Since the rank-one structure of ∆Q is preserved for updating δ edges, Theorem 2 still holds under the new settings of u and v for batch updates. Therefore, the changes ∆S to the SimRank matrix in response to δ edges insertion can be represented as a similar formulation to Theorem 3, as illustrated in the first row of Table 2. Similarly, we can also extend Theorems 6-9 in Section 5 to support batch updates of δ edges for other cases (C1)-(C3) that accompany new node insertions. Table 2 summarizes the new Q and S in response to such batch edge updates of all the cases. When δ = 1, these batch update results in Table 2 can be reduced to the unit update results of Theorems 1-9.
Algorithm 4 presents an efficient batch updates algorithm, Inc-bSR, for dynamical SimRank computation. The actual computational time of Inc-bSR depends on the input parameter ∆G since different update types in Table 2 would result in different computational time. However, we can readily show that Inc-bSR is superior to the |∆G| executions of the unit update algorithm, because Inc-bSR can process the "similar sink updates" of each block simultaneously and can cancel out redundant updates. To clarify this, let us assume that |∆G net | can be partitioned into |B| blocks, with δ t denoting the number of edge updates in t-th block. In the worst case,  we assume that all edge updates happen to be the most time-consuming case (C0) or (C2). Then, the total time for handling |∆G| updates is bounded by Note that |B| ≤ |∆G net |, in general |B| ≪ |∆G net |. Thus, Inc-bSR is typically much faster than the |∆G| executions of the unit update algorithm that is bounded by O |∆G|K(nd + ∆G + |AFF|) .
Example 7 Recall from Example 5 that a sequence of edge updates ∆G to the graph G = (V, E) in Fig. 6. We want to compute new SimRank scores in G ⊕ ∆G.
First, we can use hashing method to obtain the net update ∆G net from ∆G, as shown in Example 5.
Next, by Example 6, we can partition ∆G net into . } δ rows → row j  Then, for each block, we can apply the formulae in Table 2 to update all edges simultaneously in a batch fashion. The results are partially depicted as follows: Node The column '(q, i, +)' represents the updated SimRank scores after the edge (q, i) is added to The last column is the new SimRanks in G ⊕ ∆G. ⊓ ⊔

Memory Efficiency
In previous sections, our main focus was devoted to speeding up the computational time of incremental Sim-Rank. However, for updating all pairs of SimRank scores, the memory requirement for Algorithms 1-4 remains at O(n 2 ) since they need to store all (n 2 ) pairs of old Sim-Rank S into memory, which hinders its scalability on large graphs. We call Algorithms 1-4 in-memory algorithms.

Line Description Required Elements from old S
all elements of old S and newS Table 3: Lines of Inc-uSR (in Appendix D.1) that require to get elements from old S (highlighted in red color) In this section, we propose a novel scalable method based on Algorithms 1-4 for dynamical SimRank search, which updates all pairs of SimRanks column by column using only O(Kn + m) memory, with no need to store all (n 2 ) pairs of old SimRank S into memory, and with no loss of accuracy.
Let us first analyze the O(n 2 ) memory requirement for Algorithms 1-4 in Sections 3-5. We notice that there are two factors dominating the original O(n 2 ) memory: (1) the storage of the entire n × n old SimRank matrix S, and (2) the computation of M k from one outer product. For example, in Inc-uSR (in Appendix D.1), Lines 3, 4, 6, 9, 15 need to get elements from old S (see Table 3); Lines 10, 14, 15 require to store n × n entries of matrix M k (see Table 4). Indeed, the storage of S and M k are the main obstacles to the scalability of our in-memory algorithms on large graphs, resulting in O(n 2 ) memory space. Apart from these lines, the memory required for the remaining steps of Inc-uSR is O(m), dominated by (a) the storage of sparse matrix Q and (b) sparse matrix-vector products.
To overcome the bottleneck of the O(n 2 ) memory, our main idea is to update all pairs of S in a columnby-column style, with no need to store the entire S and M k . Specifically, we update S by updating each column [S] ⋆,x (∀x = 1, 2, · · · ) of S individually. Let us rewrite Line 15 of Table 3 into the column-wise style: Applying the following facts This implies that, to compute one column of ∆S, we only need prepare one row and one column of M K . To compute only the x-th row and x-th column of M K , there are two challenges: (1) From Line 10 of Table 3, we notice that M K is derived from the auxiliary vector γ, and γ depends on the i-th and j-th column of old S according to Lines 3, 4, 6, 9 of Table 3. Since the update edge (i, j) can be arbitrary, it is hard to determine which columns of old S will be used in future. Thus, all our in-memory algorithms in Section 5 prepare n × n elements of S into memory, leading to O(n 2 ) memory.
(2) According to Lines 10, 14, 15 of Table 4, it also requires O(n 2 ) memory to iteratively compute M K . It is not easy to use just linear memory for iteratively computing only one row and one column of M K . In the next two subsections, we will address these two challenges, respectively.

Avoid storing n × n elements of old S
Our above analysis imply that, to compute each column [∆S] ⋆,x , we only need prepare two columns information (i-th and j-th) from old S. Since the update edge (i, j) can be arbitrary, there are no prior knowledge which i-th and j-th columns in old S will be used. As opposed to Algorithms 1-4 that memoize all (n 2 ) pairs of old S, we use the following scalable method to compute only the i-th and j-th columns of old S on demand in linear memory. Specifically, based on our previous work [26] on partial-pairs SimRank retrieval, we can readily verify that the following iterations will yield [S] ⋆,i and [S] ⋆,j in just O(Kn + m) memory.  Table 3 to derive the vector γ in linear memory.

Algorithm 5: Inc-SR-All-P (G, ∆G, [S] ⋆,x , K, C)
Input : an old digraph G = (V, E), a collection of edges ∆G inserted into G, x-th column [S]⋆,x of old SimRank in G, number of iterations K, damping factor C. In addition, since Line 15 of Table 3 can be computed column-wisely via Eq.(37). Throughout all lines in Table 3, we do not need store n 2 pairs of old S in memory. However, O(n 2 ) memory is still required to store M k . In the next subsection, we will show how to avoid O(n 2 ) memory to compute M k .

Compute [M K ] ⋆,x and [M K ] x,⋆ in linear memory
Using γ, we next devise our method based on Table 4, aiming to use linear memory to compute each column [M K ] ⋆,x and each row [M K ] x,⋆ for Eq.(37). In Table 4, our key observation is that M k is the summation of   Table 4 in a column-wise style for [M K ] ⋆,x as follows: and in a row-wise style for [M K ] x,⋆ as follows: The main advantage of our method is that, throughout the entire updating process, we need not store n×n pairs of M k and S, and thereby, significantly reduce the memory usage from O(n 2 ) to O(Kn + m).
In addition to the insertion case (C0), our memoryefficient methods are applicable to other insertion cases in Subsection 5.1. The complete algorithm, denoted as Inc-SR-All-P, is described in Algorithm 5. Inc-SR-All-P is a memory-efficient version of Algorithms 1-4. It includes a procedure PartialSim that allows us to compute two columns information of old S on demand in linear memory, rather than store n 2 pairs of old S in memory. In response to each edge update (i, j), once the two old columns S ⋆,i and S ⋆,j are computed via PartialSim for updating the x-th column [∆S] ⋆,x , they can be memoized in only O(n) memory and reused later to compute another y-th column [∆S] ⋆,y in response to the edge update (i, j).
Correctness. Inc-SR-All-P correctly returns similarity. It consists of four update cases: lines 6-22 for Case (C0), lines 23-30 for Case (C1), lines 31-45 for Case (C2), and lines 46-54 for Case (C3). The correctness of each case can be verified by Theorems 3, 6, 8, and 9, respectively. For instance, to verify the correctness for Case (C0), we apply successive substitution to for-loop in lines 14-21, which produces the following result: This is consistent with Eq.(36), implying that our memoryefficient method does not compromise any accuracy for scalability. It is worth mentioning that Inc-SR-All-P can be also combined with our batch updating method in Section 6. This will speed up the dynamical update of Procedure 1: PartialSim(Q, q, K, C) Input : transition matrix Q in G, query node q, number of iterations K, damping factor C. Output: q-th column [S]⋆,q of SimRank scores in G.
SimRank further, with O(n(max |B| t=1 δ t )+m+Kn) memory. Here O(nδ t ) memory is needed to store δ t columns of S when [S] ⋆,I is required for processing the t-th block.

Experimental Evaluation
In this section, we present a comprehensive experimental study on real and synthetic datasets, to demonstrate (i) the fast computational time of Inc-SR to incrementally update SimRanks on large time-varying networks, (ii) the pruning power of Inc-SR that can discard unnecessary incremental updates outside "affected areas"; (iii) the exactness of Inc-SR, as compared with Inc-SVD; (iv) the high efficiency of our complete scheme that integrates Inc-SR with Inc-uSR-C1, Inc-uSR-C2, Inc-uSR-C3 to support link updates that allow new node insertions; (v) the fast computation time of our batch update algorithm Inc-bSR against the unit update method Inc-SR; (vi) the scalability of our memory-efficient algorithm Inc-SR-All-P in Section 7 on million-scale large graphs for dynamical updates; (vii) the performance comparison between Inc-SR-All-P and L-TSF in terms of computational time, memory space, and top-k exactness; (viii) the average updating time and memory usage of Inc-SR-All-P for each case of edge updates.

Experimental Settings
Datasets. We adopt both real and synthetic datasets. The real datasets include small-scale (DBLP and CitH), medium-scale (YouTu, WebB and WebG), and largescale graphs (CitP, SocL, UK05, and IT04). Table 5 summarises the description of these datasets.
(Please refer to Appendix E for details.) To generate synthetic graphs and updates, we adopt GraphGen 6 generation engine. The graphs are controlled  We produce a sequence of graphs that follow the linkage generation model [7]. To control graph updates, we use two parameters simulating real evolution: (a) update type (edge/node insertion or deletion), and (b) the size of updates |∆G|.
Algorithms. We implement the following algorithms: (a) Inc-SVD, the SVD-based link-update algorithm [13]; (b) Inc-uSR, our incremental method without pruning; (c) Batch, the batch SimRank method via fine-grained memoization [24]; (d) Inc-SR, our incremental method with pruning power but not supporting node insertions; (e) Inc-SR-All, our complete enhanced version of Inc-SR that allows node insertions by incorporating Inc-uSR-C1, Inc-uSR-C2, and Inc-uSR-C3; (f) Inc-bSR, our batch incremental update version of Inc-SR; (g) Inc-SR-All-P, our memory-efficient version of Inc-SR-All that dynamically computes the SimRank matrix column by column without the need to store all pairs of old similarities; (h) L-TSF, the log-based implementation of the existing competitor, TSF [20], which supports dynamic Sim-Rank updates for top-k querying.
Parameters. We set the damping factor C = 0.6, as used in [9]. By default, the total number of iterations is set to K = 15 to guarantee accuracy C K ≤ 0.0005 [16]. On CitH and YouTu, we set K = 10; On large graphs (CitP, SocL, UK05, and IT04), we set K = 5. The target rank r for Inc-SVD is a speed-accuracy trade-off, we set r = 5 in our time evaluation since, as shown in the experiments of [13], the highest speedup is achieved when r = 5. In our exactness evaluation, we shall tune this value. For L-TSF algorithm, we set the number of one-way graphs R g = 100, and the number of samples at query time R q = 20, as suggested in [20].
All the algorithms are implemented in Visual C++ and Matlab. For small-scale graphs, we use a machine with an Intel Core 2.80GHz CPU and 8GB RAM. For medium-and large-scale graphs, we use a processor with Intel Core i7-6700 3.40GHz CPU and 64GB RAM.

Time Efficiency of Inc-SR and Inc-uSR
We first evaluate the computational time of Inc-SR and Inc-uSR against Inc-SVD and Batch on real datasets.
Note that, to favor Inc-SVD that only works on small graphs (due to memory crash for high-dimension SVD n > 10 5 ), we just use Inc-SVD on DBLP and CitH. Fig.8 depicts the results when edges are added to DBLP, CitH, YouTu, respectively. For each dataset, we fix |V | and increase |E| by |∆E|, as shown in the x-axis. Here, the edge updates are the differences between snapshots w.r.t. the "year" (resp. "video age") attribute of DBLP, CitH (resp. YouTu), reflecting their real-world evolution. We observe the following. (1) Inc-SR always outperforms Inc-SVD and Inc-uSR when edges are increased. For example, on DBLP, when the edge changes are 10.7%, the time for Inc-SR (83.7s) is 11.2x faster than Inc-SVD (937.4s), and 4.2x faster than Inc-uSR (348.7s). This is because Inc-SR deploys a rankone matrix method to update the similarities, with an effective pruning strategy to skip unnecessary recomputations, as opposed to Inc-SVD that entails rather expensive costs to incrementally update the SVD. The results on CitH are more pronounced, e.g., Inc-SR is 30x better than Inc-SVD when |E| is increased to 401K. (2) Inc-SR is consistently better than Batch when the edge changes are fewer than 19.7% on DBLP, and 7.2% on CitH. When link updates are 5.3% on DBLP (resp. 3.9% on CitH), Inc-SR improves Batch by 10.2x (resp. 4.9x). This is because (i) Inc-SR can exploit the sparse structure of ∆S for incremental update, and (ii) small link perturbations may keep ∆S sparsity. Hence, Inc-SR is highly efficient when link updates are small. (3) The computational time of Inc-SR, Inc-uSR, Inc-SVD, unlike Batch, is sensitive to the edge updates |∆E|. The reason is that Batch needs to reassess all similarities from scratch in response to link updates, whereas Inc-SR and Inc-uSR can reuse the old information in SimRank for incremental updates. In addition, Inc-SVD is too sensi-  tive to |∆E|, as it entails expensive tensor products to compute SimRank from the updated SVD matrices. Fig.9 shows the target rank r required for the Li et al.'s lossless SVD approach w.r.t. the edge changes |∆E| on DBLP and CitH. The y-axis is r n × 100%. On each dataset, when increasing |∆E| from 6K to 18K, we see that r n is 95% on DBLP (resp. 80% on CitH), Thus, r is not negligibly smaller than n in real graphs. Due to the quartic time w.r.t. r, Inc-SVD may be slow in practice to get a high accuracy.
On synthetic data, we fix |V | = 79, 483 and vary |E| from 485K to 560K (resp. 560K to 485K) in 15K increments (resp. decrements). The results are shown in Fig.10. We can see that, when 6.4% edges are increased, Inc-SR runs 8.4x faster than Inc-SVD, 4.7x faster than Batch, and 2.7x faster than Inc-uSR. When 8.8% edges are deleted, Inc-SR outperforms Inc-SVD by 10.4x, Batch by 5.5x, and Inc-uSR by 2.9x. This justifies our complexity analysis of Inc-SR and Inc-uSR. Fig.11 shows the pruning power of Inc-SR as compared with Inc-uSR on DBLP, CitH, and YouTu, in which the percentage of the pruned node-pairs of each graph is depicted on the black bar. The y-axis is in a log scale. It can be discerned that, on every dataset, Inc-SR constantly outperforms Inc-uSR by nearly 0.5 order of magnitude. For instance, the running time of Inc-SR (64.9s) improves that of Inc-uSR (314.2s) by 4.8x on CitH, with approximately 82.1% node-pairs being pruned. That is, our pruning strategy is powerful to discard unnecessary node-pairs on graphs with different link distributions.

Effectiveness of Pruning
Since our pruning strategy hinges mainly on the size of the "affected areas" of the SimRank update matrix, Fig.12 illustrates the percentage of the "affected areas" of SimRank scores w.r.t. link updates |∆E| on DBLP, CitH, and YouTu. We find the following. (1) When |∆E| is varied from 6K to 18K on every real dataset, the "affected areas" of SimRank scores are fairly small. For instance, when |∆E| = 12K, the percentage of the "affected areas" is only 23.9% on DBLP, 27.5% on CitH, and 24.8% on YouTu, respectively. This highlights the effectiveness of our pruning method in real applications, where a larger number of elements of the SimRank update matrix with zero scores can be discarded. (2) For each dataset, the size of the "affect areas" mildly grows when |∆E| is increased. For example, on YouTu, the percentage of |AFF| increases from 19.0% to 24.8% when |∆E| is changed from 6K to 12K. This agrees with our time efficiency analysis, where the speedup of Inc-SR is more obvious for smaller |∆E|.

Time Efficiency of Inc-SR-All and Inc-bSR
We next compare the computational time of Inc-SR-All with Inc-SVD and Batch on DBLP, CitH, and YouTu. For each dataset, we increase |E| by |∆E| that might accompany new node insertions. Note that Inc-SR cannot deal with such incremental updates as ∆S does not make any sense in such situations. To enable Inc-SVD to handle new node insertions, we view new inserted nodes   Fig. 13 depicts the results. We can discern that (1) on every dataset, Inc-SR-All runs substantially faster than Inc-SVD and Batch when |∆E| is small. For example, as |∆E| = 6K on CitH, Inc-SR-All (186s) is 30.6x faster than Inc-SVD (5692s) and 15.1x faster than Batch (2809s). The reason is that Inc-SR-All can integrate the merits of Inc-SR with Inc-uSR-C1, Inc-uSR-C2, Inc-uSR-C3 to dynamically update SimRank scores in a rank-one style with no need to do costly matrix-matrix multiplications. Moreover, the complete framework of Inc-SR-All allows itself to support link updates that enables new node insertions.
(2) When |∆E| grows larger on each dataset, the time of Inc-SVD increases significantly faster than Inc-SR-All. This larger increase is due to the SVD tensor products used by Inc-SVD. In contrast, Inc-SR-All can effectively reuse the old SimRank scores to compute changes even if such changes may accompany new node insertions. Fig. 14 compares the computational time of Inc-bSR with Inc-SR-All. From the results, we can notice that, on each graph, Inc-bSR is consistently faster than Inc-SR-All. The last column "(%)" denotes the percentage of Inc-bSR improvement on speedup. On each dataset, the speedup of Inc-bSR is more apparent when |∆E| grows larger. For example, on DBLP, the improvement of Inc-bSR over Inc-SR-All is 8.8% when |E| = 75K, and 14.0% when |E| = 83K. On CitH (resp. YouTu), the highest speedup of Inc-bSR over Inc-SR-All is 20.7% for |E| = 419K (resp. 16.4% for |E| = 901K). This is because the large size of |∆E| may increase the number of the new inserted edges with one endpoint overlapped. Hence, more edges can be handled simultaneously by Inc-bSR, highlighting its high efficiency over Inc-SR-All. Fig. 15 evaluates the total memory usage of Inc-SR-All and Inc-bSR against Inc-SVD on real datasets. Note that the total memory usage includes the storage of the old SimRanks required for all-pairs dynamical evaluation. For Inc-SR-All, we test its three versions: (a) We first switch off our methods of "pruning" and "column-wise partitioning", denoted as "No Optimization"; (b) next turn on "pruning" only; and (c) finally turn on both. For Inc-SVD, we also tune the default target rank r = 5 larger to see how the memory space is affected by r.

Total Memory Usage
The results indicate that (1) on each dataset when the memory of Inc-SVD and Inc-bSR does not explode, the total spaces of Inc-SR-All and Inc-bSR are consistently much smaller Inc-SVD whatever target rank r is. This is because, unlike Inc-SVD, Inc-SR-All and Inc-bSR need not memorize the results of SVD tensor products.
(2) When the "pruning" switch is open, the space of Inc-SR-All can be reduced by ∼ 4x further on real data, due to our pruning method that discards many zeros in auxiliary vectors and final SimRanks. (3) When the "column-wise partitioning" switch is open, the space of Inc-SR-All can be saved by ∼ 100x further. The reason is that, as all pairs of SimRanks can be computed in a column-by-column style, there is no need to memorize the entire old SimRank S and auxiliary M. This improvement agrees with our space analysis in Section 7. (4) The space of Inc-bSR is 8-11x larger than Inc-SR-All, but is still acceptable. This is because batch updates require more space to memoize several columns from the old S to handle a subset of edge updates simultaneously. (5) For Inc-SVD, when the target rank r is varied from Datasets Inc-SR-All-P L-TSF  Fig. 18: Avg Time for Each Insertion Case 5 to 25, its total space increases from 1.36G to 3.86G on DBLP, but crashes on CitH and YouTu. This implies that r has a huge impact on the space of Inc-SVD, and is not negligible in the big-O analysis of [13].

Exactness
We next evaluate the exactness of Inc-SR-All, Inc-bSR, and Inc-SVD on real datasets. We leverage the NDCG metrics [13] to assess the top-100 most similar pairs. We adopt the results of the batch algorithm [6] on each dataset as the NDCG 100 baselines, due to its exactness. For Inc-SR-All, we evaluate its two enhanced versions: "with column-wise partitioning" and "with pruning"; for Inc-SVD, we tune its target rank r from 5 to 25. Fig. 16 depicts the results, showing the following. (1) On each dataset, the NDCG 100 s of Inc-SR-All and Inc-bSR are 1, which are better than Inc-SVD (< 0.62). This agrees with our observation that Inc-SVD may loss eigen-information in real graphs. In contrast, Inc-SR-All and Inc-bSR guarantee the exactness. (2) The NDCG 100 s for the two versions of Inc-SR-All are exactly the same, implying that both our pruning and column-wise partitioning methods are lossless while achieving high speedup.

Scalability on Large Graphs
To evaluate the scalability of our incremental techniques, we run Inc-SR-All-P, a memory-efficient version of Inc-SR, on six real graphs (WebB, WebG, CitP, SocL, UK05, and IT04), and compare its performance with L-TSF. Both Inc-SR-All-P and L-TSF can compute any single column, S ⋆,u , of S with no need to memoize all n 2 pairs of the old S. To choose the query node u, we randomly pick up 10,000 queries from each medium-sized graph (WebB and WebG), and 1,000 queries from each largesized graph (CitP, SocL, UK05, and IT04). To ensure every selected u can cover a board range of any possible queries, for each dataset, we first sort all nodes in V in descending order based on their importance that is measured by PageRank (PR), and then split all nodes into 10 buckets: nodes with PR ∈ [0.9, 1] are in the first bucket; nodes with PR ∈ [0.8, 0.9) the second, etc. For every medium-(resp. large-) sized graph, we randomly select 1,000 (resp. 100) queries from each bucket, such that u contains a wide range of various types of queries. To generate dynamical updates, we follow the settings in [20], randomly choosing 1,000 edges, and considering 80% of them as insertions and 20% deletions. It can be discerned that, on each dataset, Inc-SR-All-P is scalable well over large graphs, and runs consistently 4-7x faster than log-based L-TSF per edge update. On one-billion edge graphs (IT04), for every edge update, the updating time of Inc-SR-All-P (69.301s) is 7.3x faster than that of L-TSF (505.794s). This is because the time of L-TSF is dominated by its cost of merging R g = 100 one-way graphs' log buffers for updating the index. For example, on large IT04, almost 99.92% time required by L-TSF is due to its merge operations. In comparison, our memory-efficient method for Inc-SR-All-P takes advantage of the rank-one Sylvester equation which computes the updates to S ⋆,u in a column-by-column style on demand, without the need to merge one-way graphs and memoize all pairs of old S in advance. Fig. 18 shows the time complexities of Inc-SR-All-P for four cases of edge insertions on each real dataset. For every graph, we randomly select 1,000 edges {(i, j)} for insertion updates, with nodes i and j respectively having the probability 1/2 to be picked up from the old vertex set V . Hence, each case of edge insertion occurs at 1/4 probability. For each insertion case, we sum all the time spent in this case, and divide it by the total number of edge insertions counted for this case. Fig. 18 reports the average time per edge update for each case, together with the preprocessing time over each dataset (including the cost of loading the graph and preparing its transition matrix Q). From the results, we see that, on each dataset, the time spent for Cases (C0) and (C2) is moderately higher than that for Case (C1); Case (C0) is slightly slower than Case (C2); Case (C3) entails the lowest time cost. These results are consistent with our intuition and mathematical formulation of ∆S for each case. Case (C0) has the most expensive time cost as it

Precision
To compare the precision of Inc-SR-All-P and L-TSF, we define the precision measure [10] for top-k querying: The original batch algorithm in [9] (resp. [13]) serves as the exact solution to obtain SimRank results for L-TSF (resp. Inc-SR-All-P). We evaluate the precision of both algorithms on several real datasets. Fig. 19 reports the result on YouTu; the tendencies on other datasets are similar. We see that, when top-k varies from 10 to 10 5 , the precision of L-TSF remains high (> 84%) for small top-k (< 1000), but is lower (68%-75%) for large top-k (> 10 4 ). This is because the probabilistic guarantee for the error bound of L-TSF is based on the assumption that no cycle in the given graph has a length shorter than K (the total number of steps). Hence, L-TSF is highly efficient for top-k single-source querying, where k is not large. In contrast, the precision of Inc-SR-All-P is stable at 1, meaning that it will produce exact Sim-Rank results of [13], regardless of top-k values. Thus, Inc-SR-All-P is better for non top-k query.

Memory of
Inc-SR-All-P Fig. 20 evaluates the memory usage of Inc-SR-All-P and L-TSF over six real datasets. We observe that both algorithms scale well on large graphs. On WebB, IT04, and UK05, the memory space of Inc-SR-All-P is almost the same as L-TSF; On WebG, CitP, and SocL, the memory usage of Inc-SR-All-P is 5-8x less than L-TSF. This is because, unlike L-TSF that need load a one-way graph to memory, Inc-SR-All-P only requires to prepare the vector information of ξ k , η k , old S ⋆,i , and old S ⋆,j to assess the changes to each column of S in response to edge update (i, j). The memory space of these auxiliary vectors can sometimes be comparable to the size of the one-way graph, and sometimes be much smaller. However, such memory space is linear to n as we do not need n 2 space to store the entire old S. Note that the old S ⋆,j and S ⋆,i can be computed on demand with only linear memory by our partial-pairs SimRank approach [26]. Moreover, we see that, with the growing scale of the real datasets, the memory space of Inc-SR-All-P is increasing linearly, highlighting its scalability on large graphs. Fig. 21 depicts further the average memory usage of Inc-SR-All-P for each case of edge insertion. We randomly pick up 1,000 edges {(i, j)} for insertion updates on each dataset, with nodes i and j respectively having the probability 1/2 to be chosen from the old vertex set V . The average memory space of Inc-SR-All-P for each case is reported in Fig. 21. We see that, on each dataset, the memory required for Cases (C0), (C1), and (C2) are similar, whereas the memory space of Case (C3) is much smaller than the other cases. The reason is that, for Cases (C0), (C1), and (C2), Inc-SR-All-P needs linear memory to store some auxiliary vectors (e.g., ξ k , η k , y, old S ⋆,i , and old S ⋆,j ) for updating SimRank scores, whereas for Case (C3), no auxiliary vectors are required for precomputation, thus saving much memory space.

Incremental SimRank
Li et al. [13] devised an interesting matrix representation of SimRank, and was the first to show a SVD method for incrementally updating all pairs of Sim-Ranks, which requires O(r 4 n 2 ) time and O(r 2 n 2 ) memory, where r (≤ n) is the target rank of the low-rank approximation. However, their incremental techniques are inherently inexact, with no guaranteed accuracy.
Recently, Shao et al. [20] provided an excellent exposition of a two-stage random sampling framework, TSF, for top-k SimRank dynamic search w.r.t. one query u. In the preprocessing stage, they sampled a collection of one-way graphs to index random walks in a scalable manner. In the query stage, they retrieved similar nodes by pruning unqualified nodes based on the connectivity of one-way graph. To retrieve top-k nodes with highest SimRank scores in a single column S ⋆,u , [20] requires O(K 2 R q R g ) average query time to retrieve S ⋆,u along with O(n log k) time to return top-k results from S ⋆,u . The recent work of Jiang et al. [10] has argued that, to retrieve S ⋆,u , the querying time of [20] is O(KnR q R g ). The n factor is due to the time to traverse the reversed one-way graph; in the worst case, all n nodes are visited. Moreover, Jiang et al. [10] observed that the probabilistic error guarantee of Shao et al.'s method is based on the assumption that no cycle in the given graph has a length shorter than K, and they proposed READS, a new efficient indexing scheme that improves precision and indexing space for dynamic SimRank search. The query time of READS is O(rn) to retrieve one column S ⋆,u , where r is the number of sets of random walks. Hence, TSF and READS are highly efficient for top-k single-source SimRank search. In comparison, our dynamical method focuses on all (n 2 )-pairs SimRank search in O(K(m+|AFF|)) time. Optimization methods in this work are based on a rank-one Sylvester matrix equation characterising changes to n 2 pairs of SimRank scores, which is fundamentally different from [10,20]'s methods that maintain one-way graphs (or SA forests) updating. It is important to note that, for large-scale graphs, our incremental methods do not need to memoize all (n 2 ) pairs of old SimRank scores, and can dynamically update S column-wisely in only O(Kn + m) memory. For updating each column of S, our experiments in Section 8 verify that our memory-efficient incremental method is scalable on large real graphs while running 4-7 times faster than the dynamical TSF [20] per edge update, due to the high cost of [20] merging one-way graph's log buffers for TSF indexing.
There has also been a body of work on incremental computation of other graph-based relevance measures. Banhmani et al. [1] utilized the Monte Carlo method for incrementally computing Personalized PageRank. Desikan et al. [3] proposed an excellent incremental PageRank algorithm for node updating. Their central idea revolves around the first-order Markov chain. Sarma et al. [19] presented an excellent exposition of randomly sampling random walks of short length, and merging them together to estimate PageRank on graph streams.

Batch SimRank
In comparison to incremental algorithms, the batch Sim-Rank computation has been well-studied on static graphs. For deterministic methods, Jeh and Widom [9] were the first to propose an iterative paradigm for SimRank, entailing O(Kd 2 n 2 ) time for K iterations, where d is the average in-degree. Later, Lizorkin et al. [16] utilized the partial sums memoization to speed up SimRank computation to O(Kdn 2 ). Yu et al. [24] have also improved SimRank computation to O(Kd ′ n 2 ) time (with d ′ ≤ d) via a fine-grained memoization to share the common parts among different partial sums. Fujiwara et al. [6] exploited the matrix form of [13] to find the top-k similar nodes in O(n) time w.r.t. a given query node. All these methods require O(n 2 ) memory to output all pairs of SimRanks. Recently, Kusumoto et al. [11] proposed a linearized method that requires only O(dn) memory and O(K 2 dn 2 ) time to compute all pairs of SimRanks. The recent work of [26] proposes an efficient aggregate method for computing partial pairs of SimRank scores. The main ideas of partial-pairs SimRank search are also incorporated into the incremental model of our work, achieving linear memory to update n 2 -pairs similarities.
For parallel SimRank computing, Li et al. [15] proposed a highly parallelizable algorithm, called Cloud-Walker, for large-scale SimRank search on Spark with ten machines. Their method consists of offline and online phases. For offline processing, an indexing vector is derived by solving a linear system in parallel. For online querying, similarities are computed instantly from the index vector. Throughout, the Monte Carlo method is used to maximally reduce time and space.
The recent work of Zhang et al. [28] conducted extensive experiments and discussed in depth many existing SimRank algorithms in a unified environment using different metrics, encompassing efficiency, effectiveness, robustness, and scalability. The empirical study for 10 algorithms from 2002 to 2015 shows that, despite many recent research efforts, the running time and precision of known algorithms have still much space for improvement. This work makes a further step towards this goal.
Fogaras and Rácz [5] proposed P-SimRank in linear time to estimate a single-pair SimRank s(a, b) from the probability that two random surfers, starting from a and b, will finally meet at a node. Li et al. [14] harnessed the random walks to compute local SimRank for a single node-pair. Later, Lee et al. [12] employed the Monte Carlo method to find top-k SimRank node-pairs. Tao et al. [22] proposed an excellent two-stage way for the top-k SimRank-based similarity join.

Conclusions
In this article, we study the problem of incrementally updating SimRank scores on time-varying graphs. Our complete scheme, Inc-SR-All, consists of five ingredients: (1) For edge updates that do not accompany new node insertions, we characterize the SimRank update matrix ∆S via a rank-one Sylvester equation. Based on this, a novel efficient algorithm is devised, which reduces the incremental computation of SimRank from O(r 4 n 2 ) to O(Kn 2 ) for each link update. (2) To eliminate unnecessary SimRank updates further, we also devise an effective pruning strategy that can improve the incremental computation of SimRank to O(K(m + |AFF|)), where |AFF| (≪ n 2 ) is the size of the "affected areas" in the SimRank update matrix. (3) For edge updates that accompany new node insertions, we consider three insertion cases, according to which end of the inserted edge is a new node. For each case, we devise an efficient incremental SimRank algorithm that can support new node insertions and accurately update the affected similarities. (4) For batch updates, we also propose efficient batch incremental methods that can handle "similar sink edges" simultaneously and eliminate redundant edge updates. (5) To optimize the memory for all-pairs SimRank updates, we also devise a columnwise memory-efficient technique that significantly reduces the storage from O(n 2 ) to O(Kn + m), without the need to memoize n 2 pairs of SimRank scores. Our experimental evaluations on real and synthetic datasets demonstrate that (a) our incremental scheme is consistently 5-10 times faster than Li et al. 's SVD based method; (b) our pruning strategy can speed up the incremental SimRank further by 3-6 times; (c) the batch update algorithm enables an extra 5-15% speedup, with just a little compromise in memory; (d) our memoryefficient incremental method is scalable on billion-edge graphs; for every edge update, its computational time can be 4-7 times faster than L-TSF and its memory space can be 5-8 times less than L-TSF; (e) for different cases of edge updates, Cases (C0) and (C2) entail more time than Case (C1), and Case (C3) runs the fastest.

Appendix A Limitation of Li et al.'s SVD [13]
We rigorously explain the reason why Li et al.'s incremental method may miss some eigen-information even if a lossless SVD is utilized for SimRank computation.
Let us first revisit the main idea of Li et al.'s incremental method [13]. Briefly, [13] characterizes SimRank matrix S in Eq.(1) in terms of three matrices U, Σ, V, where U, Σ, V are derived by the SVD of Q, i.e., Then, when links are changed, [13] incrementally computes the new SimRank matrixS by updating the old matrices U, Σ, V respectively as where U C , Σ C , V C are derived from the SVD of the auxiliary matrix C Σ + U T · ∆Q · V, i.e., and ∆Q is the changes to Q in response to link updates. However, the main problem is that the derivation of Eq.(39) rests on the assumption that Unfortunately, Eq.(41) does not hold (unless Q is a fullrank matrix, i.e., rank(Q) = n) because in the case of rank(Q) < n, even a "perfect" (lossless) SVD of Q via Eq.(38) would produce n × α rectangular matrices U and V with α = rank(Q) < n. Thus, rank(U · U T ) = α < n = rank(I n ), which implies that U · U T = I n . Similarly, V · V T = I n when rank(Q) < n. Hence, Eq.(41) is not always true, as visualized in Fig. 22.
Example 8 Consider a graph with the matrix Q = [ 0 1 0 0 ], and its lossless SVD: One can readily verify that  Step 1. Initially, when links are changed, the old Q is updated to newQ = Q + ∆Q. By replacing Q with Eq.(38), it follows that Step 2. Premultiply by U T and postmultiply by V on both sides of Eq.(42), and then apply the property Step 3. Let C be the right-hand side of Eq.(43). Applying Eq.(40) to Eq.(43) yields Step 4. Li et al. [13] attempted to premultiply by U and postmultiply by V T on both sides of Eq.(44) first, and then rested on the assumption of Eq.(41) to obtain which is the result of Eq.(39).
However, the problem lies in Step 4. As mentioned before, Eq.(41) does not hold when rank(Q) < n, which means thatQ =Ũ ·Σ ·Ṽ T in Eq.(45). Consequently, updating the old U, Σ, V via Eq.(39) may produce an error (up to I n −U·U T 2 = 1, which is not practically small) for incrementally "approximating" S.

Example 9
Recall the old Q and its SVD in Example 8. Suppose there is a new edge insertion, associated with ∆Q = [ 0 0 1 0 ]. [13] first computes auxiliary matrix C as Then, the matrix C is decomposed via Eq.(40) into Finally, [13] updates the new SVD ofQ via Eq.(39) as However, one can readily verify that In comparison, a "true" SVD ofQ should bẽ Besides, the approximation error is not small in practice Our analysis suggests that, only when (i) Q is fullrank, and (ii) the SVD of Q is lossless (n = rank(Q) = α), Li et al.'s incremental way [13] can produce the exact S, but the time complexity of [13], O(r 4 n 2 ), would become O(n 6 ), which is prohibitively expensive. In practice, as evidenced by our statistical experiments in Fig.9 on Stanford Large Network Datasets (SNAP), most real graphs are not full-rank, highlighting our need to devise an efficient method for dynamic SimRank computation.
T = u · w T + w · u T , with w = y + λ 2 u.
(b) We next convert the recursive form of ∆S into the series form. One can readily verify that Thus, based on Eq.(51), the recursive definition of ∆S in Eq.(48) naturally leads itself to the series form: Combining this with Eq.(50) yields ∆S = ∞ k=0 C k+1 ·Q k · u · w T + w · u T · (Q T ) To simplify Q · S · [Q] T j,⋆ in y, and [Q] j,⋆ · S · [Q] T j,⋆ in λ, we postmultiply both sides of Eq.(1) by e j to obtain We also premultiply both sides of Eq.(53) by e T j to get Combining this with Theorem 2 completes the proof of the case when d j = 1.
(ii) If d j > 1, then all nonzeros in old [Q] j,⋆ are 1 dj . The deleted edge (i, j) will update [Q] j,⋆ via 2 steps: first, all nonzeros in [Q] j,⋆ are changed from 1 dj to 1 dj −1 ; then, the entry [Q] j,i is changed from 1 dj to 0.
For k > 0, one can readily prove that the k-th iterative M k in Line 14 of Algorithm 6 is the first k-th partial sum of M in Eq. (13). Thus, M k+1 can be derived from M k as follows: Thus, the (a, b)-entry form of the above equation is To show that [M k ] a,b = 0 for (a, b) / ∈ A 0 ×B 0 ∪A k ×B k , we follow the 2 steps: (i) For (a, b) / ∈ A 0 × B 0 , as proved in the case k = 0, the term C · [e j ] a [γ] b in the above equation is obviously 0. (ii) For each node x ∈ T , we then find all out-neighbors of x in G, which produces F 1 , i.e., F 1 = x∈T O(x). Analogously, the set F 2 in Eq. (21), in the case of d j > 0, consists of the nodes "⋆" depicted in (18). When d j = 0, F 2 = ∅ as the term [S] ⋆,i is not in the expression of γ in Eq. (14), in contrast to the case d j > 0.
After obtaining F 1 and F 2 , we can readily find A 0 × B 0 , according to Eq. (22). For k > 0, to iteratively derive the node-pair set A k × B k , we take the following two steps: (i) we first construct a node-pair set T 1 × T 2 := {(x, y)|[M k−1 ] x,y = 0}. (ii) For every node x ∈ T 1 (resp. y ∈ T 2 ), we then find all out-neighbors of x (resp. y) in G ∪ {(i, j)}, which yields A k (resp. B k ), i.e., A k = x∈T1Õ (x) and B k = y∈T2Õ (y).
The node selectivity of Theorem 5 hinges on ∆S sparsity. Since real graphs are constantly updated with minor changes, ∆S is often sparse in general. Hence, many node-pairs with zero scores in ∆S can be discarded. As demonstrated by our experiments in Fig.11, 76.3% paper-pairs on DBLP can be pruned, significantly reducing unnecessary similarity recomputations.

C.1 Li et al.s SVD incremental approach
Example 10 Figure 3 depicts a citation graph G, a tiny fraction of DBLP, where each node is a paper, and an edge represents a reference from one paper to another. Suppose G is updated by adding an edge (i, j), denoted by ∆G (see the dash arrow). Using the damping factor C = 0.8, we would like to compute SimRank scores in the new graph G ∪ ∆G.
The results are compared in the table of Figure 3, where Column 'sim Li et al. ' denotes the approximation of SimRank scores returned by Li et al.'s Algorithm 3 [13], and Column 'sim true ' denotes the "true" SimRank scores returned by a batch algorithm [6] that runs in G ∪ ∆G from scratch. It can be noticed that for some nodepairs (not highlighted in gray), the similarities obtained by Li et al. 's incremental method are different from the "true" SimRank scores even if lossless SVD is used 10 during the process of updating U, Σ, V T . This suggests that Li et al.'s incremental approach [13] is inherently approximate. In fact, their incremental strategy would neglect some useful eigen-information whenever rank(Q) < n.
We also notice that the target rank r for the SVD of the matrix C 11 is not always negligibly smaller than n. For example, in Column 'sim Li et al. ' of Figure 3, r is chosen to be rank(C) = 9 to get a lossless SVD of C. Although r = 9 appears not negligibly smaller than n = 15, the accuracy of 'sim Li et al. ' is still undesirable as compared with 'sim true ', not to mention using r < 9. ⊓ ⊔ Example 10 implies that Li et al.'s incremental approach [13] is approximate and may produce high computational overheads since r is not always much smaller.
10 A rank-α SVD of the matrix X ∈ R n×n is a factorization of the form Xα = U · Σ · V T , where U, V ∈ R n×α are columnorthonormal matrices, and Σ ∈ R α×α is a diagonal matrix, α is called the target rank of the SVD, as specified by the user. If α = rank(X), then Xα = X, and we call it the lossless SVD. If α < rank(X), then X − Xα 2 gives the least square estimate error, and we call it the low-rank SVD. 11 As defined in [13], r is the target rank for the SVD of the auxiliary matrix C Σ + U T · ∆Q · V, where ∆Q is the changes to Q for link updates.