Algorithmica (2008) 51: 387–427 DOI 10.1007/s00453-007-9051-4

Mantaining Dynamic Matrices for Fully Dynamic Transitive Closure Camil Demetrescu · Giuseppe F. Italiano

Received: 22 November 2005 / Accepted: 20 February 2007 / Published online: 13 October 2007 © Springer Science+Business Media, LLC 2007

Abstract In this paper we introduce a general framework for casting fully dynamic transitive closure into the problem of reevaluating polynomials over matrices. With this technique, we improve the best known bounds for fully dynamic transitive closure. In particular, we devise a deterministic algorithm for general directed graphs that achieves O(n2 ) amortized time for updates, while preserving unit worst-case cost for queries. In case of deletions only, our algorithm performs updates faster in O(n) amortized time. We observe that fully dynamic transitive closure algorithms with O(1) query time maintain explicitly the transitive closure of the input graph, in order to answer each query with exactly one lookup (on its adjacency matrix). Since an update may change as many as (n2 ) entries of this matrix, no better bounds are possible for this class of algorithms. Keywords Dynamic graph algorithms · Transitive closure

This work has been partially supported by the Sixth Framework Programme of the EU under contract number 507613 (Network of Excellence “EuroNGI: Designing and Engineering of the Next Generation Internet”), and number 001907 (“DELIS: Dynamically Evolving, Large Scale Information Systems”), and by the Italian Ministry of University and Research (Project “ALGO-NEXT: Algorithms for the Next Generation Internet and Web: Methodologies, Design and Experiments”). Portions of this paper have been presented at the 41st Annual Symp. on Foundations of Computer Science, 2000. C. Demetrescu () Dipartimento di Informatica e Sistemistica, Università di Roma “La Sapienza”, Rome, Italy e-mail: [email protected] G.F. Italiano Dipartimento di Informatica, Sistemi e Produzione, Università di Roma “Tor Vergata”, Rome, Italy e-mail: [email protected]

388

Algorithmica (2008) 51: 387–427

1 Introduction In this paper we present fully dynamic algorithms for maintaining the transitive closure of a directed graph. A dynamic graph algorithm maintains a given property on a graph subject to dynamic updates, such as edge insertions and edge deletions. We say that an algorithm is fully dynamic if it can handle both edge insertions and edge deletions. A partially dynamic algorithm can handle either edge insertions or edge deletions, but not both: we say that it is incremental if it supports insertions only, and decremental if it supports deletions only. In the fully dynamic transitive closure problem we wish to maintain a directed graph G = (V , E) under an intermixed sequence of the following operations: Insert(x, y): insert an edge from x to y in G; Delete(x, y): delete the edge from x to y in G; Query(x, y): report yes if there is a path from x to y in G, and no otherwise. Throughout the paper, we denote by m and by n the number of edges and vertices in G, respectively. Previous Work Research on dynamic transitive closure spans over two decades. Before describing the results known, we list the bounds obtainable with simple-minded methods. If we do nothing during each update, then we have to explore the whole graph in order to answer reachability queries: this gives O(n2 ) time per query and O(1) time per update in the worst case. On the other extreme, we could recompute the transitive closure from scratch after each update; as this task can be accomplished via matrix multiplication [1, 18], this approach yields O(1) time per query and O(nω ) time per update in the worst case, where ω is the best known exponent for matrix multiplication (currently ω < 2.38 [2]). For the incremental version of the problem, the first algorithm was proposed by Ibaraki and Katoh [10] in 1983: its running time was O(n3 ) over any sequence of insertions. This bound was later improved to O(n) amortized time per insertion by Italiano [11] and also by La Poutré and van Leeuwen [16]. Yellin [22] gave an O(m∗ δmax ) algorithm for m edge insertions, where m∗ is the number of edges in the final transitive closure and δmax is the maximum out-degree of the final graph. All these algorithms maintain explicitly the transitive closure, and so their query time is O(1). The first decremental algorithm was again given by Ibaraki and Katoh [10], with a running time of O(n2 ) per deletion. This was improved to O(m) per deletion by La Poutré and van Leeuwen [16]. Italiano [12] presented an algorithm that achieves O(n) amortized time per deletion on directed acyclic graphs. Yellin [22] gave an O(m∗ δmax ) algorithm for m edge deletions, where m∗ is the initial number of edges in the transitive closure and δmax is the maximum out-degree of the initial graph. Again, the query time of all these algorithms is O(1). More recently, Henzinger and King [9] gave a randomized decremental transitive closure algorithm for general directed graphs with a query time of O(n/ log n) and an amortized update time of O(n log2 n).

Algorithmica (2008) 51: 387–427

389

The first fully dynamic transitive closure algorithm was devised by Henzinger and King [9] in 1995: they gave a randomized Monte Carlo algorithm with onesided error supporting a query time of O(n/ log n) and an amortized update time of O(nm ˆ 0.58 log2 n), where m ˆ is the average number of edges in the graph throughout the whole update sequence. Since m ˆ can be as large as O(n2 ), their update time 2 2.16 is O(n log n). Khanna, Motwani and Wilson [13] proved that, when a lookahead of (n0.18 ) in the updates is permitted, a deterministic update bound of O(n2.18 ) can be achieved. Recently, King and Sagert [15] showed how to support queries in O(1) time and updates in O(n2.26 ) time for general directed graphs and O(n2 ) time for directed acyclic graphs; their algorithm is randomized with one-sided error. In a major breakthrough [14], King improved the bounds of King and Sagert, by exhibiting a deterministic algorithm on general digraphs with O(1) query time and O(n2 log n) amortized time per update operations, where updates are insertions of a set of edges incident to the same vertex and deletions of an arbitrary subset of edges. We remark that all these algorithms (except [14]) use fast matrix multiplication as a subroutine. We observe that fully dynamic transitive closure algorithms with O(1) query time maintain explicitly the transitive closure of the input graph, in order to answer each query with exactly one lookup (on its adjacency matrix). Since an update may change as many as (n2 ) entries of this matrix, O(n2 ) seems to be the best update bound that one could hope for this class of algorithms. It is thus quite natural to ask whether the O(n2 ) update bound can be actually realized for fully dynamic transitive closure on general directed graphs while maintaining one lookup per query. Our Results In this paper, we affirmatively answer this question. We exhibit a deterministic algorithm for fully dynamic transitive closure on general digraphs that does exactly one matrix look-up per query and supports updates in O(n2 ) amortized time, thus improving over [14]. Our algorithm can also support within the same time bounds the generalized updates of [14], i.e., insertion of a set of edges incident to the same vertex and deletion of an arbitrary subset of edges. In the special case of deletions only, our algorithm achieves O(n) amortized time for deletions and O(1) time for queries: this generalizes to directed graphs the bounds of [12], and improves over [9]. In [5], we show how to break through the O(n2 ) barrier in the case of directed acyclic graphs. Our results are based on a novel technique: we introduce a general framework for maintaining polynomials defined over matrices, and we cast fully dynamic transitive closure into this framework. In particular, our algorithm hinges upon the equivalence between transitive closure and matrix multiplication on a closed semiring; this relation has been known for over 30 years (see, e.g., the results of Munro [18], Furman [8] and Fischer and Meyer [7]) and yields the fastest known static algorithm for transitive closure. Surprisingly, no one before seems to have exploited this equivalence in the dynamic setting: some recent algorithms [9, 13, 15] make use of fast matrix multiplication, but only as a subroutine for fast updates. Differently from other approaches, the crux of our method is to use dynamic reevaluation of products of Boolean matrices as the kernel for solving dynamic transitive closure. The remainder of this paper is organized as follows. In Sect. 2 we first define the fully dynamic transitive closure problem and show that it is equivalent to a fully

390

Algorithmica (2008) 51: 387–427

dynamic Boolean matrix closure problem. In Sect. 3 we address the problem of maintaining polynomials defined over Boolean matrices during updates of their variables. In Sect. 4 we present a first method for casting fully dynamic transitive closure into the dynamic problem on matrices considered in Sect. 3: as a warm-up, we revisit within this framework the algorithm of King [14] from a completely different perspective. In Sect. 5 we show how the algebraic approach of Sect. 3 can be further refined: this yields our improved algorithm for fully dynamic transitive closure. Finally, in Sect. 6 we list some concluding remarks.

2 Fully Dynamic Transitive Closure In this section we give a more formal definition of the fully dynamic transitive closure problem considered in this paper. We assume the reader to be familiar with standard graph and algebraic terminology as contained for instance in [1, 3]. Definition 1 Let G = (V , E) be a directed graph and let TC(G) = (V , E  ) be its transitive closure. The F ULLY DYNAMIC T RANSITIVE C LOSURE P ROBLEM consists of maintaining a data structure G for graph G under an intermixed sequence of I NITIALIZATION, U PDATE, and Q UERY operations. Each operation can be either one of the following: • Init(A): perform the initialization operation E ← A, where A ⊆ V × V . • Insert(v, I ): perform the update E ← E ∪ {(x, y) ∈ I | x = v ∨ y = v }, where I ⊆ V × V and v ∈ V . We call this update a v-C ENTERED insertion in G. • Delete(D): perform the update E ← E − D, where D ⊆ E. • Query(x, y): perform a query operation on TC(G) by returning 1 if (x, y) ∈ E  and 0 otherwise. A few remarks are in order at this point. First, the generalized Insert and Delete updates considered here have been first introduced by King in [14]. With just one operation, they are able to change the graph by adding or removing a whole set of edges, rather than a single edge, as illustrated in Fig. 1. Second, we consider explicitly initializations of the graph G and, more generally than in the traditional definitions of dynamic problems, we allow them to appear everywhere in the sequence of operations. This gives more generality to the problem, and allows for more powerful data structures, i.e., data structures that can be restarted at run time on a completely different input graph. Differently from other variants of the problem, we do not address the issue of returning actual paths between nodes, and we just consider the problem of answering reachability queries. It is well known that, if G = (V , E) is a directed graph and XG is its adjacency ∗ of X is equivalent to computing the (rematrix, computing the Kleene closure XG G flexive) transitive closure TC(G) of G. We remind the reader that the Kleene closure of XG is defined as follows: ∗ XG

=

∞  i=0

(XG )i .

Algorithmica (2008) 51: 387–427

391

Fig. 1 a Insert operation; b Delete operation as in Definition 1 ∗ , instead of considIn this paper, because of the equivalence between TC(G) and XG ering directly the problem introduced in Definition 1, we study an equivalent problem on matrices. Before defining it formally, we need some preliminary notation.

Definition 2 If X is a matrix, we denote by IX,i and JX,j the matrices equal to X in the i-th row and j -th column, respectively, and zero elsewhere:  X[x, y] if x = i, IX,i [x, y] = 0 otherwise,  X[x, y] if y = j, JX,j [x, y] = 0 otherwise. Definition 3 Let X and Y be two n × n Boolean matrices. We say that X ⊆ Y if and only if X[x, y] = 1 ⇒ Y [x, y] = 1 for any x, y ∈ {1, . . . , n}. Definition 4 Let X and Y be two n × n Boolean matrices. We denote by X + Y the bitwise or of X and Y , by X ⊕ Y the bitwise xor of X and Y , and by X · Y the Boolean matrix multiplication of X and Y . We are now ready to define a dynamic version of the problem of computing the Kleene closure of a Boolean matrix. Definition 5 Let X be an n × n Boolean matrix and let X ∗ be its Kleene closure. We define the F ULLY DYNAMIC B OOLEAN M ATRIX C LOSURE P ROBLEM as the problem of maintaining a data structure X for matrix X under an intermixed sequence of initialization, update, and query operations. Each operation can be either one of the following: • Init∗ (Y ): perform the initialization operation X ← Y , where Y is an n × n Boolean matrix.

392

Algorithmica (2008) 51: 387–427

• Set∗ (i, X): perform the update X ← X + IX,i + JX,i , where X is an n × n Boolean update matrix and i ∈ {1, . . . , n}. We call this kind of update an i-C ENTERED set operation on X. • Reset∗ (X): perform the update X ← X ⊕ X, where X ⊆ X is an n × n Boolean update matrix. • Lookup∗ (x, y): return the value of X ∗ [x, y], where x, y ∈ {1, . . . , n}. Notice that Set∗ is allowed to modify only the i-th row and the i-th column of X, while Reset∗ and Init∗ can modify any entries of X. We stress the strong correlation between Definition 1 and Definition 5: if G is a graph and X is its adjacency matrix, operations G.Init, G.Insert, G.Delete, and G.Query are equivalent to operations X.Init∗ , X.Set∗ , X.Reset∗ , and X.Lookup∗ , respectively. We remark that for this to be true it is crucial that X.Reset∗ is defined only for X ⊆ X.

3 Dynamic Matrices In this section we address the problem of reevaluating polynomials over Boolean matrices under modifications of their variables. We propose a data structure for maintaining efficiently the special class of polynomials of degree 2 consisting of single products of Boolean matrices. We show then how to use this data structure for solving the more general problem on arbitrary polynomials. This problem will be central to designing efficient algorithms for the fully dynamic Boolean matrix closure problem. 3.1 Dynamic Reevaluation of Polynomials over Boolean Matrices Let Bn be the set of n × n Boolean matrices and let P=

h 

T a = T 1 + T2 + · · · + Th

a=1

be a polynomial1 with h terms defined over Bn , where each Ta =

k 

Xba = X1a · X2a · · · Xka

b=1

has degree exactly k and variables Xba ∈ Bn are distinct. We now study the problem of maintaining the value of polynomial P under updates of variables Xba . We will denote by Y the maintained value of polynomial P , assuming that Y can be different 1 In the following, we omit specifying explicitly the dependence of a polynomial on its variables, and when-

ever the correct interpretation is clear from the context, we denote by P both the function P (X1 , . . . , Xk ) and the value of this function for fixed values of X1 , . . . , Xk .

Algorithmica (2008) 51: 387–427

393

from P : as we will see later on, tolerating errors in Y will allow us to support updates for our original problem of dynamic transitive closure much faster than by using an exact implementation. In this setting, we will consider two kinds of updates on our data structures: L AZY U PDATES and S TRICT U PDATES, defined as follows. Definition 6 We say that an update is S TRICT if it satisfies the following properties for any x, y: 1. if P [x, y] flips from 0 to 1 because of the update, then Y [x, y] flips from 0 to 1 as well; 2. if P [x, y] = 0 after the update then Y [x, y] = 0 after the update; We also say that an update is L AZY if it satisfies Property 2, but not Property 1 above. The above definition implies that strict updates affect immediately the maintained value Y of the polynomial P , i.e., any change in P is immediately reflected in Y , while lazy updates may cause Y to be different from P , in the sense that 1’s that appear in P due to a lazy update may not necessarily appear immediately in Y . If initially Y = P , then both strict and lazy updates maintain the invariant that Y ⊆ P . We now need some preliminary definitions. Definition 7 We denote by 0n the n × n null matrix such that 0n [x, y] = 0 for each x, y and by In the n × n unit matrix such that In [x, y] = 1 for x = y and In [x, y] = 0 if x = y. Definition 8 Let X be a Boolean matrix. We denote by X(t) the value of X at time t, i.e., the value of X after the t-th operation in a sequence of operations that modify X. By convention, we assume that at time 0 any value in X is zero, i.e., X(0) = 0n . Definition 9 Let Bn be the set of n × n Boolean matrices and let T = X1 · · · Xk be a product of matrices in Bn . For each entry T [x, y] = 1, there must be at least one sequence of indices πxy = x = u0 , u1 , . . . , uk−1 , uk = y such that Xj [uj −1 , uj ] = 1, 1 ≤ j ≤ k. We define any such sequence πxy to be a W ITNESS PATH for the entry T [x, y] = 1. Note that an entry T [x, y] = 1 may have several witness paths. We now characterize better how the maintained value of the polynomial Y can be different from its exact value P , by introducing the following notation. Definition 10 Let X be a Boolean matrix. We denote by L AST S ET (X[x, y]) the time of the latest strict update that set entry X[x, y] of matrix X to 1, including updates that overwrite a value already set to 1. Definition 11 Let X be a Boolean matrix. We denote by L AST F LIP (X[x, y]) the time the entry X[x, y] of matrix X last flipped from 0 to 1. We observe that, while LastSet can only be affected by strict updates, LastFlip may be changed by both strict and lazy updates. Note that:

394

Algorithmica (2008) 51: 387–427

• LastSet(X[x, y]) < LastFlip(X[x, y]) when the latest update on X[x, y] was a lazy update. • LastSet(X[x, y]) = LastFlip(X[x, y]) when the latest update on X[x, y] was a strict update that flipped X[x, y] from 0 to 1. • LastSet(X[x, y]) > LastFlip(X[x, y]) when the latest update on X[x, y] was a strict update that overwrote X[x, y] again to 1. Definition 12 Let Bn be the set of n × n Boolean matrices and let T = X1 · · · Xk be a product of matrices in Bn subject to a sequence of updates such that at each time exactly one matrix is updated. Let πxy = x = u0 , u1 , . . . , uk−1 , uk = y be a witness path for T [x, y] = 1. We define S TRICT W ITNESS FOR πxy the unique pair of indices (uq−1 , uq ), if any, such that both the following conditions hold: LastSet(Xq [uq−1 , uq ]) = max LastSet(Xj [uj −1 , uj ]), 1≤j ≤k

LastSet(Xq [uq−1 , uq ]) ≥ max LastFlip(Xj [uj −1 , uj ]). 1≤j ≤k

In this case, we also say that Xq [uq−1 , uq ] is a S TRICT W ITNESS FOR T [x, y] = 1. Notice that, if T [x, y] = 1 has a strict witness, then necessarily Y [x, y] = P [x, y] = 1. On the other hand, if T [x, y] = 1 has no strict witness, it may be Y [x, y] = 0 = P [x, y] = 1: we claim that this can only happen if the witness paths of T [x, y] = 1 appeared due to lazy updates. Indeed, let πxy = x = u0 , u1 , . . . , uk−1 , uk = y be a witness path for T [x, y] = 1. If there is no strict witness for πxy , then by Definition 12: max LastSet(Xj [uj −1 , uj ]) < max LastFlip(Xj [uj −1 , uj ]),

1≤j ≤k

1≤j ≤k

which implies that πxy appeared due to a lazy update. We now formally introduce our problem on polynomials. Definition 13 Let Bn be the set of n × n Boolean matrices and let P=

h 

Ta

a=1

be a polynomial with h terms defined over Bn , where each Ta =

k 

Xba

b=1

has degree exactly k and variables Xba ∈ Bn are distinct. We consider the problem of maintaining a data structure P for the polynomial P under an intermixed sequence of initialization, update, and query operations. Each operation can be either one of the following:

Algorithmica (2008) 51: 387–427

395

• Init(Z11 , . . . , Zkh ): clean up the data structure and perform the initialization Xba ← Zba of the variables of polynomial P , where each Zba is an n × n Boolean matrix. • SetRow(i, X, a, b): perform the row update operation Xba ← Xba + IX,i , where X is an n × n Boolean update matrix. In other words, the operation sets to 1 the entries in the i-th row of variable Xba of polynomial P as specified by matrix X. • SetCol(i, X, a, b): perform the column update operation Xba ← Xba + JX,i , where X is an n × n Boolean update matrix. In other words, the operation sets to 1 the entries in the i-th column of variable Xba of polynomial P as specified by matrix X. • LazySet(X, a, b): perform the update operation Xba ← Xba + X, where X is an n × n Boolean update matrix. In other words, the operation sets to 1 the entries of variable Xba of polynomial P as specified by matrix X. This operation does not affect immediately the maintained value Y of the polynomial. • Reset(X, a, b): perform the update operation Xba ← Xba ⊕ X, where X is an n × n Boolean update matrix such that X ⊆ Xba . In other words, the operation resets to 0 the entries of variable Xba of polynomial P as specified by matrix X. • Lookup(): return a matrix Y such that Y [x, y] = 1 if and only if there exists an entry Ta [x, y] with at least one strict witness (see Definition 12). Init, SetRow, SetCol and Reset are strict updates, while LazySet is a lazy update. We note that SetRow and SetCol are allowed to modify only the i-th row and the i-th column of variable Xba , respectively, while LazySet, Reset and Init can modify all entries of Xba . Note that only Init, SetRow and SetCol may affect LastSet, while LastFlip may be affected by Init, SetRow, SetCol, and LazySet. It is crucial to observe that the definition of Lookup allows one-sided errors in answering queries on the value of P . In particular, let Y be the matrix returned by Lookup. If Y [x, y] = 1 then necessarily P [x, y] = 1. However, if P [x, y] = 1, but no entry Ta [x, y] has a strict witness, then Y [x, y] = 0. We notice that, if no LazySet operation is ever performed, then every entry Ta [x, y] = 1 will have at least one strict witness: in this case Y = P and Lookup makes no errors by always returning the correct value of polynomial P . Thus, the presence of errors is related to the presence of LazySet operations in the sequence. It is not difficult to extend the results of this section to the general class of polynomials with terms of different degrees and multiple occurrences of the same variable. In the following, for the sake of simplicity, we will sometimes abuse the notation: for instance, we will use SetRow(i, X, X) instead of SetRow(i, X, a, b) to denote a SetRow operation on each occurrence of variable X in the polynomial. Before considering the general case where polynomials have arbitrary degree k, we focus on the special class of polynomials of degree 2. 3.1.1 Data Structure for Polynomials of Degree k = 2 In this section we define a data structure for a polynomial P of degree 2 that allows us to maintain explicitly the value Y (t) of the matrix Y returned by Lookup at any

396

Algorithmica (2008) 51: 387–427

Table 1 Data structures for polynomials of degree 2. Times are amortized

Simple-minded

Witness sets [4]

This paper

Init

O(hnω )

O(hnω )

O(hn3 )

SetRow

O(hn2 )

O(hn2 )

O(n2 )

SetCol

O(hn2 )

O(hn2 )

O(n2 )

LazySet

O(n2 )

O(n2 )

O(n2 )

Reset

O(nω )

O(n2 )

O(1)

Lookup

O(n2 )

O(n2 )

O(n2 )

Space

O(hn2 )

O(hn3 )

O(hn2 )

time t during a sequence of operations. This makes it possible to perform Lookup operations in optimal quadratic time. We propose efficient techniques for propagating to Y the effects of changes of variables Xba due to SetRow, SetCol and Reset operations. In case of LazySet, we only need to update the maintained value of the affected variables, leaving the other portions of the data structure untouched. This, implies that after a LazySet at time t, the maintained value Y (t) will not be synchronized with the correct value P (t) of the polynomial. Most technical difficulties of this section come specifically from this lazy maintenance of Y (t). Before giving the details of our data structure, we briefly list the bounds that can be obtained by simple-minded algorithms. Let P = X11 · X21 + · · · + X1h · X2h be the polynomial to be maintained under the sequence of operations. We keep explicitly P and the value of each term Ta = X1a ·X2a , for a = 1, 2, . . . , h, thus requiring a total of O(hn2 ) space. An Init operation can be simply implemented by h Boolean matrix products and (h − 1) matrix sums, and thus requires O(hn2 + hnω ) = O(hnω ) time. SetRow (respectively SetCol) requires multiplying a row (respectively a column) by a matrix, and performing (h − 1) matrix sums: this can be clearly done in O(hn2 ) time. LazySet can be implemented in O(n2 ) time by simply updating the matrix Xab , without propagating the changes. Reset can be supported by updating the matrix Xab , by recomputing Ta via fast matrix multiplication, and by performing (h − 1) matrix sums to rebuild the value of P : this requires O(nω + hn2 ) time. Finally, Lookup simply returns the maintained value of P in O(n2 ) time. Note that with this approach, Reset operations are much slower than SetRow, SetCol, LazySet and Lookup. In [4] we have shown how to achieve O(n2 ) amortized time per Reset operation at the expense of increasing to O(hn3 ) the space required by the data structure. The main idea is to keep explicitly for each a = 1, 2, . . . , h all the witness paths x, y, z of the Boolean product Ta = X1a ·X2a that contain a strict witness. In this section, we show how to further reduce the amortized bound of Reset and to reduce the space to O(hn2 ) by using a technique that we call lazy counting. The bounds obtained in this paper are illustrated in Table 1. Our data structure for representing a polynomial of degree 2 is presented below. Data Structure 1 We maintain the following elementary data structures with O(h · n2 ) space:

Algorithmica (2008) 51: 387–427

397

1. a counter T ime of the number of performed operations; 2. 2h Boolean matrices X1a and X2a for 1 ≤ a ≤ h; 3. h integer matrices Prod1 , . . . , Prodh such that Proda [x, y] ≤ n maintains a count of the number of strict witnesses for term entry Ta [x, y] = (X1a · X2a )[x, y] (see Definition 12). 4. an integer matrix S such that S[x, y] = |{a : Proda [x, y] > 0}|. Note that Y [x, y] = 1 ⇔ S[x, y] > 0, and therefore we do not need to maintain explicitly matrix Y . 5. 2h integer matrices LastFlipX , one for each matrix X = Xba . For each entry X[x, y] = 1, LastFlipX [x, y] stores the time of the most recent operation that caused X[x, y] to flip from 0 to 1 (i.e., LastFlipX [x, y] = LastFlip(X[x, y]) as in Definition 11). 6. 3h integer vectors LastLev1a , LastLev2a , LastLev3a such that: • LastLev1a [i] stores the time of the latest Init operation, or SetRow operation on the i-th row of X1a ; • LastLev2a [i] stores the time of the latest Init operation, or SetCol operation on the i-th column of X1a , or SetRow operation on the i-th row of X2a ; • LastLev3a [i] stores the time of the latest Init operation, or SetCol operation on the i-th column of X2a . Note that LastLev1a , LastLev2a , LastLev3a encode information about LastSet as in Definition 10. Namely, for each entry X1a [x, y] = 1 we have: LastSet(X1a [x, y]) = max{LastLev1a [x], LastLev2a [y]} and for each entry X2a [y, z] = 1 it holds: LastSet(X2a [y, z]) = max{LastLev2a [y], LastLev3a [z]}. Before getting into the full details of our implementation of operations, we give an overview of the main ideas. We consider how the various operations should affect the data structure. In particular, we suppose that an operation changes some entries of variable X1a in a term Ta = X1a · X2a , and we define what our implementation should do on matrix Proda : SetRow/SetCol: if some entry X1a [x, y] is set to 1, then it becomes a strict witness in the product X1a · X2a for any pair x, z such that X2a [y, z] = 1. Then we should increase the count Proda [x, z] by 1, if neither X1a [x, y] nor X2a [y, z] was already counted for. LazySet: this operation never creates strict witnesses, and thus we do not need to update Proda . Reset: if some entry X1a [x, y] is set to 0, then we should decrease the count Proda [x, z] by 1, if either X1a [x, y] or X2a [y, z] was previously counted for. We proceed similarly if any entry X2a [y, z] is set to 0. Note that after performing LazySet there may be triples (x, y, z) such that both X1a [x, y] = 1 and X2a [y, z] = 1, but neither of them is counted in Proda [x, z]. To

398

Algorithmica (2008) 51: 387–427

tackle this case, we introduce a predicate Pa (x, y, z), 1 ≤ x, y, z ≤ n, such that Pa (x, y, z) is true if and only if the witness path x, y, z has a strict witness, i.e.: Pa (x, y, z) := max{LastFlipXa [x, y], LastFlipXa [y, z]} 1

2

≤ max{LastSet(X1a [x, y]), LastSet(X2a [y, z])}. Since for each entry X1a [x, y] = 1 it holds LastSet(X1a [x, y]) = max{LastLev1a [x], LastLev2a [y]} and for each entry X2a [y, z] = 1 we have LastSet(X2a [y, z]) = max{LastLev2a [y], LastLev3a [z]}, this is equivalent to: Pa (x, y, z) := max{LastFlipXa [x, y], LastFlipXa [y, z]} 1

2

≤ max{LastLev1a [x], LastLev2a [y], LastLev3a [z]}. We remark that we do not need to maintain Pa explicitly in our data structure as it can be computed on demand in constant time by simply accessing LastFlip, LastLev1, LastLev2, and LastLev3. Property Pa allows it to define the following invariant that we maintain in our data structure. Invariant 1 For any term Ta = X1a · X2a in polynomial P , at any time during a sequence of operations, the following invariant holds for any pair of indices x, z: Proda [x, z] = |{y : X1a [x, y] = 1 ∧ X2a [y, z] = 1 ∧ Pa (x, y, z)}|. In accordance with Invariant 1, it should be clear that the value of each entry Proda [x, z] counts the number of strict witnesses of the Boolean matrix product Ta [x, z] = (X1a · X2a )[x, z]. We implement the operations introduced in Definition 13 as described next, assuming that the operation Time ← Time + 1 is performed just before each operation:

Init Init cleans up the data structure, initializes variables X1a and X2a , and builds the other portions of Data Structure 1. In particular, LastFlipX [x, y] is set to Time for any X[x, y] that becomes 1 and LastLev1a [i], LastLev2a [i], and LastLev3a [i] are set to Time, for any i and a. Proda is initialized by computing integer matrix multiplication X1a · X2a , i.e., by looking at Xba as integer 0/1 matrices.

Algorithmica (2008) 51: 387–427

399

SetRow procedure SetRow(i, X, a, b) 1. begin 2. Xba ← Xba + IX,i 3. if b = 1 then 4. for each x s.t. X1a [i, x] = 1 do 5. for each y s.t. X2a [x, y] = 1 do 6. if not Pa (i, x, y) then 7. Proda [i, y] ← Proda [i, y] + 1 8. if Proda [i, y] = 1 then S[i, y] ← S[i, y] + 1 9. LastLev1a [i] ← Time 10. else 11. for each x : X1a [x, i] = 1 do 12. for each y : X2a [i, y] = 1 do 13. if not Pa (x, i, y) then 14. Proda [x, y] ← Proda [x, y] + 1 15. if Proda [x, y] = 1 then S[x, y] ← S[x, y] + 1 16. LastLev2a [i] ← Time 17. update LastFlipXa b 18. end

{b=2}

After a SetRow(i, X, a, 1) operation (b = 1), each entry X1a [i, x] = 1 is a strict witness for witness paths of the form i, x, y. Clearly, Pa (i, x, y) is true after the update. If Pa (i, x, y) was false before the update, then Proda [i, y] (and possibly S[i, y]) is increased by one. Similarly, after a SetRow(i, X, a, 2) operation (b = 2), each entry X2a [i, y] = 1 becomes a strict witness for witness paths of the form x, i, y. Clearly, Pa (x, i, y) is true after the update. If Pa (x, i, y) was false before the update, then Proda [x, y] (and possibly S[x, y]) is increased by one.

SetCol The operation is completely symmetrical to SetRow.

LazySet procedure LazySet(X, a, b) 1. begin 2. Xba ← Xba + X {LastLev1a , LastLev2a and LastLev3a are not changed} 3. update LastFlipXa b 4. end

LazySet simply sets to 1 any entries in Xba as specified by X, and updates LastFlipXa . We remark that no other portion of the data structure is changed. b

400

Algorithmica (2008) 51: 387–427

Reset procedure Reset(X, a, b) 1. begin 2. if b = 1 then 3. for each x, y s.t. X[x, y] = 1 do 4. if max{LastLev1a [x], LastLev2a [y]} ≥ LastFlipXa [x, y] then 1 5. for each z s.t. X2a [y, z] = 1 do 6. if Pa (x, y, z) then 7. Proda [x, z] ← Proda [x, z] − 1 8. if Proda [x, z] = 0 then S[x, z] ← S[x, z] − 1 9. else {here max{LastLev1a [x], LastLev2a [y]} < LastFlipXa [x, y]} 1 10. for each z s.t. X2a [y, z] = 1 ∧ LastLev3a [z] ≥ LastFlipXa [x, y] do 1 11. if Pa (x, y, z) then 12. Proda [x, z] ← Proda [x, z] − 1 13. if Proda [x, z] = 0 then S[x, z] ← S[x, z] − 1 14. else b = 2 similar to b = 1 15. Xba ← Xba ⊕ X 16. end

The objective of the operation is to reset the entries of Xba as specified by X (line 15). Before doing so in line 15, in lines 2–14 it updates Proda and S so as to maintain Invariant 1, using information encoded in LastLev1a , LastLev2a , LastLev3a , and LastFlipXa . We only describe the case b = 1, since the case b = 2 is completely b analogous. For each pair (x, y) such that X[x, y] = 1 (line 3), Reset looks for triples (x, y, z) such that Pa (x, y, z) is invalidated by resetting Xba [x, y] (lines 5–6 and lines 10–11). Proda and S are adjusted accordingly (lines 7–8 and lines 12–13). The distinction between the two cases max{LastLev1a [x], LastLev2a [y]} ≥ LastFlipXa [x,y] and max{LastLev1a [x], LastLev2a [y]} < LastFlipXa [x,y] in line 4 and 1 1 in line 9, respectively, is important to achieve O(n2 ) amortized running times as it will be discussed in the proof of Theorem 2. Here we only point out that if the test in line 4 succeeds, then we can scan any z s.t. X2a [y, z] = 1 without affecting the claimed running time. If this is not the case, then we process only indices z such that the condition LastLev3a [z] ≥ LastFlipXa [x, y] is satisfied, and avoid scanning 1 other indices. For this reason line 10 must be implemented very carefully by maintaining indices z in a list and by using a move-to-front strategy that brings index z to the front of the list as any operation Init(. . .), SetRow(z, . . .) or SetCol(z, . . .) is performed on z. In this way indices are sorted in accordance with the dates of operations on them. Notice that we can safely skip lines 11–13 if the test at line 10 fails since Pa (x, y, z) cannot be true in that case. Indeed, the test at line 10 fails if it is either X2a [y, z] = 0 or LastLev3a [z] < LastFlipXa [x, y]. If X2a [y, z] = 0 then clearly 1 Pa (x, y, z) is false. On the other side, if LastLev3a [z] < LastFlipXa [x, y], then 1

max{LastLev1a [x], LastLev2a [y], LastLev3a [z]} < LastFlipXa [x, y] ≤ max{LastFlipXa [x, y], LastFlipXa [y, z]}, 1

1

2

Algorithmica (2008) 51: 387–427

401

and thus Pa (x, y, z) is false. This corresponds to the case where the witness path

x, y, z has been previously created by a LazySet, but this information has not been propagated yet to the maintained value of the polynomial, and thus Reset has no other work to do besides resetting Xba [x, y] to 0.

Lookup function Lookup() 1. begin 2. return matrix Y s.t. Y [x, y] = 1 ⇔ S[x, y] > 0 3. end

Lookup simply returns a binarized version Y of matrix S defined in Data Structure 1. Analysis The correctness of our implementation is discussed in the following theorem. Theorem 1 Operations SetRow, SetCol, LazySet, Reset and Init maintain a matrix S such that S[x, y] > 0 if and only if there exists some entry Ta [x, y] = 1 with at least one strict witness (see Definition 13). Proof The implementation of the operations on our data structure maintains Invariant 1 for each a, x, y by means of simple bookkeeping operations: Proda [x, y] = |{z : X1a [x, z] = 1 ∧ X2a [z, y] = 1 ∧ Pa (x, z, y)}|. It is also easy to see that the following invariant is maintained as well: S[x, y] > 0 ⇔ ∃a such that Proda [x, y] > 0. Thus, S[x, y] > 0 if and only if there exist two indices a and z such that X1a [x, z] = 1 ∧ X2a [z, y] = 1 ∧ Pa (x, z, y), which is equivalent to saying that there exists some entry Ta [x, y] = 1 with at least one strict witness.  The previous theorem implies that Lookup correctly returns a matrix Y such that Y [x, y] = 1 if and only if there exists some entry Ta [x, y] = 1 with at least one strict witness. We now analyze the running time of our implementation of the operations on polynomials. Theorem 2 Lookup, SetRow, SetCol and LazySet operations require O(n2 ) time. Init requires O(h·n3 ) time, while Reset requires O(1) time. All time bounds are amortized. The space required is O(h · n2 ). Proof The space bound derives from the fact that we have to maintain O(h) matrices, each of size O(n2 ). To prove the time bounds, we use a credit-based argument. Each SetRow (respectively SetCol) deposits n credits on each entry in the row (respectively column) of

402

Algorithmica (2008) 51: 387–427

Xba affected by the operation, while Init deposits n credits on each entry of a variable Xba . This gives a total of O(n2 ) credits for SetRow and SetCol, and O(h · n3 ) credits for Init. Since Lookup just returns Y , it requires O(n2 ) actual time. It is straightforward to see from the pseudocode of the operations that any SetRow, SetCol and LazySet operation requires O(n2 ) actual time. Since SetRow and SetCol requires O(n2 ) credits and Lookup and LazySet require no credits, their amortized cost is O(n2 ). Init takes O(h · nω ) actual time: in more detail, each Proda can be directly computed via matrix multiplication and any other initialization step requires no more than O(n2 ) worst-case time. The O(h · n3 ) amortized bound derives from the deposited credits. To complete the proof, we show that the actual cost of each Reset can be entirely payed for by the credits deposited by the other operations. Consider the distinction between the two cases max{LastLev1a [x], LastLev2a [y]} ≥ LastFlipXa [x, y] in line 4 1 and max{LastLev1a [x], LastLev2a [y]} < LastFlipXa [x, y] in line 9. In the first case, 1 we can pay the cost of scanning triples (x, y, z) with the credits on X1a [x, y]. In the second case, we only scan those triples (x, y, z) for which some Init or SetCol was performed on the z-th column of X2a [y, z] after both X1a [x, y] and X2a [y, z] were set to 1. This cost can be payed for with the credits on X2a [y, z].  The amortized analysis of Theorem 2 allows Reset operations to charge up to a O(h · n3 ) cost to a single Init operation. Thus, in an arbitrary mixed sequence with any number of Init, Init takes O(h · n3 ) amortized time per update. If, however, we allow Init operations to appear in the sequence only every (n) Reset operations, the bound for Init drops down to O(h · n2 ) amortized time per operation. As a consequence of Theorem 2, we have the following corollaries that refine the analysis of the running time of Init operations. Corollary 1 If we perform just one Init operation in a sequence of (n) operations, or more generally one Init operation every (n) Reset operations, then the amortized cost of Init is O(h · n2 ) per operation. Corollary 2 If we perform just one Init operation in a sequence of (n2 ) operations, or more generally one Init operation every (n2 ) Reset operations, and no operations SetRow and SetCol, then the amortized cost of Init is O(h · n) per operation. In the following, we show how the general case of polynomials of degree k > 2 can be reduced to the case k = 2 considered in this section. 3.1.2 Data Structure for Polynomials of Degree k > 2  of To support terms of degree k > 2 in P , we consider an equivalent representation P P such that the degree of each term is 2. This allows us to maintain a data structure  with the operations defined in the previous section. for P

Algorithmica (2008) 51: 387–427

403

Lemma 1 Consider a polynomial P=

h 

Ta =

a=1

h 

X1a · · · Xka

a=1

with h terms where each term Ta has degree exactly k and variables Xba are Boolean  be the polynomial over Boolean matrices of degree 2 defined as matrices. Let P = P

h  k 

a Lab,b−1 · Rb,k−b−1 ,

a=1 b=0 a are polynomials over Boolean matrices of degree ≤ 2 defined as where Lab,j and Rb,j

 Lab,j =  a Rb,j =

a · Lab,j −1 Xb−j In a a Rb,j −1 · Xb+1+j In

if j ∈ [0, b − 1], if j = −1, if j ∈ [0, k − b − 1], if j = −1.

. Then P = P Proof To prove the claim, it suffices to check that Ta =

k 

a Lab,b−1 · Rb,k−b−1 .

b=0

Unwinding the recursion for Lab,b−1 , we obtain: Lab,b−1 = X1a · Lab,b−2 = X1a · X2a · Lab,b−3 = · · · = X1a · X2a · · · Xba · In . a a Likewise, Rb,k−b−1 = In · Xb+1 · · · Xka holds. Thus, by idempotence of the closed semiring of Boolean matrices, we finally have: k  b=0

a Lab,b−1 · Rb,k−b−1 =

k  b=0

a X1a · · · Xba · Xb+1 · · · Xka = X1a · · · Xka = Ta .



a = I · Xa a  a Note that Lab,0 = Xba · In = Xba and Rb,0 n b+1 = Xb+1 . Since P , Lb,j and a are all polynomials of degree ≤ 2, they can be represented and maintained effiRb,j ciently with instances of Data Structure 1 presented in Sect. 3.1.1. Our data structure for maintaining polynomials of degree > 2 is presented below: a Data Structure 2 We maintain explicitly the O(h · k 2 ) polynomials Lab,j and Rb,j  exactly, we for j > 0 with instances of Data Structure 1. Rather than maintaining P tolerate one-sided errors. Namely, we maintain a matrix Y such that: if Y [x, y] = 1 [x, y] = 1; if P [x, y] = 1, but no term Ta [x, y] of P  has a strict then necessarily P witness, then Y [x, y] = 0. Y is maintained with an instance of Data Structure 1.

404

Algorithmica (2008) 51: 387–427

Fig. 2 Updating a term T of degree k with a SetCol operation on variable Xb , or with a SetRow operation on variable Xb+1

We now consider how to support SetRow, SetCol, LazySet, Reset, Init and Lookup in the case of arbitrary degree. We denote by SetRow2 , SetCol2 , and LazySet2 the versions of SetRow, SetCol, and LazySet implemented for polynomials of degree 2.

SetCol 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.

procedure SetCol(i, X, Xba ) begin Lab,1 .SetCol2 (1, X, a, 1) for j ← 2 to b − 1 do a {it holds Lab,j = Xb−j · Lab,j −1 } Lab,j .SetCol2 (i, Lab,j −1 .Lookup(), a, 2) for j ← 2 to k − b − 1 do a a a a Ra b,j .SetRow2 (i, Rb,j −1 .Lookup(), a, 1) { it holds Rb,j = Rb,j −1 · Xb+1+j } a Y.SetCol2 (i, Lb,b−1 .Lookup(), a, 1) a .Lookup(), a, 2) Y.SetRow2 (i, Rb,k−b−1 for j ← 1 to k − b do Lab+j,j .LazySet2 (X, a, 2) a for j ← 1 to b − 2 do Rb−j −1,j .LazySet2 (X, a, 1) end

The main idea behind SetCol (see Fig. 2) is to exploit the associativity of Boolean matrix multiplication in order to propagate changes of intermediate polynomials that are always limited to a row or to a column: this can be efficiently handled by means of operations SetRow2 and SetCol2 .

Algorithmica (2008) 51: 387–427

405

In lines 3–4 SetCol propagates via SetCol2 the changes of the i-th column of Xba to Lab,1 , then the changes of the i-th column of Lab,1 to Lab,2 , and so on throughout the recursive decomposition: = Xba Lab,0 = Xba · In a a a a Lb,1 = Xb−1 · Lb,0 = Xb−1 · Xba a a a a a Lb,2 = Xb−2 · Lb,1 = Xb−2 · Xb−1 · Xba .. .. .. . . . a a · Xb−1 · Xba Lab,b−1 = X1a · Lab,b−2 = X1a · · · Xb−2 Likewise, in lines 5–6 it propagates via SetRow2 a null matrix of changes of the a a , then the changes (possibly none) of the i-th row of R a to Rb,1 i-th row of Xb+1 b,1 a , and so on through the (due to the late effects of some previous LazySet) to Rb,2 recursive decomposition: a = In · Xba Rb,0 a a · Xa Rb,1 = Rb,0 b+1 a a · Xa Rb,2 = Rb,1 b+2 .. .. . . a a = Rb,k−b−2 · Xka Rb,k−b−1

a = Xb+1 a a = Xb+1 · Xb+2 a a a = Xb+1 · Xb+2 · Xb+3 .. . a a a = Xb+1 · Xb+2 · Xb+3 · · · Xka

We remark that both loops in lines 3–4 and in lines 5–6 reveal, gather and propagate any 1’s that appear in the intermediate polynomials due to the late effects of some previous LazySet. In particular, even if the presence of lines 5–6 may seem a = 0n , these lines are executed just to propagate the apparently strange since Xb+1 effects of some LazySet. a are propagated to Y , and Finally, in lines 7–8 changes of Lab,b−1 and Rb,k−b−1 in lines 9–10 new 1’s are lazily inserted in any other polynomials that feature Xba as a variable. The pseudocode for SetRow is similar to SetCol, and hence it is omitted.

Reset, LazySet, Init, Lookup Reset(X, Xba ) can be supported by propagating via Reset2 any changes of Xba a that contains it, then changes of such to any intermediate polynomial Lau,v and Ru,v polynomials to any polynomials that depend on them, and so on up to Y . LazySet(X, Xba ) can be supported by performing LazySet2 operations on each a that contains X a . polynomial Lau,v and Ru,v b Init(Z11 , . . . , Zkh ) can be supported by calling Init2 on each polynomial Lw u,v , w and by propagating the intermediate results up to Y . Ru,v . Lookup() can be realized by returning the maintained value Y of P Analysis To conclude this section, we discuss the correctness and the running time of our operations in the case of polynomials of arbitrary degree. The correctness of our implementation (see Definition 13) hinges on the following theorem.

406

Algorithmica (2008) 51: 387–427

Theorem 3 Operations SetRow, SetCol, LazySet, Reset and Init maintain  matrix Y such that Y [x, y] = 1 if and only if there exists some entry Ta [x, y] in P with at least one strict witness (see Definition 12). Proof We first observe that Init operations rebuild from scratch the data structure, Reset calls Reset2 on each maintained polynomial, and set operations (SetRow, SetCol, and LazySet) never add more 1’s than required. This implies that Y ⊆ P . We now prove that, if there is an entry Ta [x, y] with a strict witness, then it must be Y [x, y] = 1. Let (ub−1 , ub ) be the strict witness, and t = LastSet(Xba [ub−1 , ub ]) be the time of the latest strict update that set entry Xba [ub−1 , ub ] to 1. The strict set update at time t can only be one of the following: • SetCol(ub , Xba , Xba ). Since (ub−1 , ub ) is a strict witness for Ta [x, y] at time t, then there must be a witness path x = u0 , u1 , u2 , . . . , ub−1 , ub , . . . , uk = y, i.e., Xja [uj −1 , uj ] = 1 for every j ∈ [1, k]. After the first iteration of the first loop in SetCol, Lab,1 .Lookup()[ub−2 , ub ] = 1; after the second iteration, Lab,2 .Lookup()[ub−3 , ub ] = 1, etc. At the end of the loop, Lab,b−1 .Lookup() a [x, ub ] = 1 (see Fig. 2). Similarly, at the end of the second loop, Rb,k−b−1 .Lookup()[ub , y] = 1. Thus, after the set operations on Y in lines 7–8, Y [x, y] = a [ub , y] = 1. Lab,b−1 [x, ub ] · Rb,k−b−1 • SetRow(ub−1 , Xba , Xba ). Analogous to SetCol. To conclude the proof, we show that, if Y [x, y] = 1 then there must be a strict witness for some Ta [x, y]. Assume by contradiction that no Ta [x, y] has a strict witness. Then either P [x, y] = 0, which implies Y [x, y] = 0, or P [x, y] = 1 and all witness paths for each Ta [x, y] appeared due to LazySet operation, which do not affect Y . Again, it must be Y [x, y] = 0, clearly a contradiction.  Theorem 4 Lookup requires O(n2 ) time, while SetRow, SetCol and LazySet require O(k · n2 ) time. Init requires O(h · k 2 · n3 ) time, while Reset requires O(1) time. All time bounds are amortized. The space required is O(h · k 2 · n2 ). Proof The data structure maintains O(h · k 2 ) polynomials of degree 2 of the form L · R, to which all bounds of Theorem 2 apply with h = 1. In particular, each polynomial requires O(n2 ) space, and thus the total space requirement is O(h · k 2 · n2 ). We now discuss the time bounds. Since Lookup just returns Y , it requires O(n2 ) time. SetCol performs O(b) SetCol2 operations, O(k − b) SetRow2 operations, O(k) LazySet operations, and O(k) Lookup; since b = O(k), the O(k · n2 ) bound follows directly from Theorem 2. A similar argument holds for SetRow. LazySet performs O(k) LazySet2 , while Init performs O(h · k 2 ) Init2 operations, and thus they require O(k · n2 ) and O(h · k 2 · n3 ) amortized time, respectively. Finally, the cost of any Reset operation can be charged against previous SetRow, SetCol, and Init operations, exactly as in the proof of Theorem 2.  As in the previous section, we have the following corollaries.

Algorithmica (2008) 51: 387–427

407

Corollary 3 If we perform just one Init operation in a sequence of (n) operations, or more generally one Init operation every (n) Reset operations, then the amortized cost of Init is O(h · k · n2 ) per operation. Corollary 4 If we perform just one Init operation in a sequence of (n2 ) operations, or more generally one Init operation every (n2 ) Reset operations, and we perform no operations SetRow and SetCol, then the amortized cost of Init is O(h · k · n) per operation.

4 Transitive Closure Updates in O(n2 log n) Time In this section we show a first method for casting fully dynamic transitive closure into the problem of reevaluating polynomials over Boolean matrices presented in Sect. 3.1. Based on the technique developed in Sect. 3.1, we revisit the dynamic graph algorithm of King in [14] in terms of dynamic matrices and we present a matrixbased variant of it which features better initialization time while maintaining the same bounds on the running time of update and query operations, i.e., O(n2 · log n) time per update and O(1) time per query. The space requirement of our algorithm is M(n) · log n, where M(n) is the space used for representing a polynomial over Boolean matrices. As stated in Theorem 4, M(n) is O(n2 ) if h and k are constant. This will be a warmup for our improved algorithm of Sect. 5. In the remainder of this section we first describe our data structure and then we show how to support efficiently operations introduced in Definition 5 for the equivalent problem of fully dynamic Boolean matrix closure. 4.1 Data Structure As it is well known, the Kleene closure of a Boolean matrix X can be computed from scratch via matrix multiplication by computing log2 n polynomials Pk = 2 , 1 ≤ k ≤ log n. In the static case where X ∗ has to be computed only Pk−1 + Pk−1 2 once, intermediate results can be thrown away as only the final value X ∗ = Plog2 n is required. In the dynamic case, instead, intermediate results provide useful information for updating efficiently X ∗ whenever X gets modified. In this section we consider a slightly different definition of polynomials P1 , . . . , Plog2 n with the property that each of them has degree ≤ 3: Definition 14 Let X be an n × n Boolean matrix. We define the sequence of log2 n + 1 polynomials over Boolean matrices P0 , . . . , Plog2 n as:  Pk =

X 2 +P3 Pk−1 + Pk−1 k−1

if k = 0, if k > 0.

Before describing our data structure for maintaining the Kleene closure of X, we discuss some useful properties.

408

Algorithmica (2008) 51: 387–427

Lemma 2 Let X be an n ×n Boolean matrix and let Pk be formed as in Definition 14. Then for any 1 ≤ u, v ≤ n, Pk [u, v] = 1 if and only if there is a path u  v of length at most 3k in X. Proof The proof is by induction on k. The base (k = 0) is trivial. We assume by induction that the claim is satisfied for Pk−1 and we prove that it is satisfied for Pk as well. Sufficient condition: Any path of length up to 3k between u and v in X is either of length up to 3k−1 or it can be obtained as concatenation of three paths of length up to 3k−1 in X. Since all these paths are correctly reported in Pk−1 by the inductive 2 [u, v] = 1 or P 3 [u, v] = 1. Thus hypothesis, it follows that Pk−1 [u, v] = 1 or Pk−1 k−1 3 2 Pk [u, v] = Pk−1 [u, v] + Pk−1 [u, v] + Pk−1 [u, v] = 1. Necessary condition: If Pk [u, v] = 1 then at least one among Pk−1 [u, v], 2 [u, v] and P 3 [u, v] must be 1. If P Pk−1 k−1 [u, v] = 1, then by the inductive hyk−1 2 [u, v] = 1, then there are pothesis there is a path of length up to 3k−1 < 3k . If Pk−1 two paths of length up to 3k−1 whose concatenation yields a path no longer than 3 [u, v] = 1, then there are three paths of length up to 3k−1 whose 3k . Finally, if Pk−1 concatenation yields a path no longer than 3k .  Lemma 3 Let X be an n ×n Boolean matrix and let Pk be formed as in Definition 14. Then X ∗ = In + Plog2 n . Proof The proof easily follows from Lemma 2 and from the observation that the length of the longest simple path in X is no longer than n − 1 < 3log3 n ≤ 3log2 n . In is required to guarantee the reflexivity of X ∗ .  Our data structure for maintaining X ∗ is the following: Data Structure 3 We maintain an n × n Boolean matrix X and we maintain the log2 n polynomials P1 . . . Plog2 n of degree 3 given in Definition 14 with instances of Data Structure 2 presented in Sect. 3.1.2. 3 As we will see in Sect. 4.2, the reason for considering the extra term Pk−1 in our data structure is that polynomials need to be maintained using not only SetRow/SetCol, but also LazySet. As stated in Definition 13, using LazySet yields a weaker representation of polynomials, and this forces us to increase the degree if complete information about X ∗ has to be maintained. This aspect will be discussed in more depth in the proof of Theorem 5.

4.2 Implementation of Operations In this section we show that operations Init∗ , Set∗ , Reset∗ and Lookup∗ introduced in Definition 5 can all be implemented in terms of operations Init, LazySet, SetRow, and SetCol (described in Sect. 3.1.2) on polynomials P1 . . . Plog2 n .

Algorithmica (2008) 51: 387–427

409

In the following pseudocode we use SetRow(i, X, X) instead of SetRow (i, X, a, b) to denote a SetRow operation on each occurrence of variable X in the polynomial. We use similar notation for the other operations.

Init∗ 1. 2. 3. 4. 5. 6.

procedure Init∗ (X) begin Z←X for k = 1 to log2 n do Pk .Init(Z) Z ←Pk .Lookup() end

Init∗ performs Pk .Init operations on each Pk by propagating intermediate results from X to P1 , then from P1 to P2 , and so on up to Plog2 n .

Lookup∗ procedure Lookup∗ (x, y) 1. begin 2. Y ←Plog2 n .Lookup() 3. return In + Y [x, y] 4. end

Lookup∗ returns the value of (In + Plog2 n )[x, y] = X ∗ [x, y].

Set∗ 1. 2. 3. 4. 5. 6. 7. 8.

procedure Set∗ (i, X) begin Z ← X for k = 1 to log2 n do Pk .LazySet(Z, Pk−1 ) Pk .SetRow(i, Z, Pk−1 ) Pk .SetCol(i, Z, Pk−1 ) Z ←Pk .Lookup() end

Set∗ propagates changes of Pk−1 to Pk for any k = 1 to log2 n. Notice that any new 1’s that appear in Pk−1 are inserted in Pk via LazySet, but only the i-th row and the i-th row column of Pk−1 are taken into account by SetRow and SetCol in order to determine changes of Pk . As re-inserting 1’s already present in a variable is allowed by our operations on polynomials, for the sake of simplicity in line 7 we assign the update matrix Z with Pk and not with the variation of Pk .

410

Algorithmica (2008) 51: 387–427

Reset∗ 1. 2. 3. 4. 5. 6. 7.

procedure Reset∗ (X) begin Z ← X for k = 1 to log2 n do Z ←Pk .Lookup() Pk .Reset(Z, Pk−1 ) Z ← Y −Pk .Lookup() end

Reset∗ performs Pk .Reset operations on each Pk by propagating changes specified by X to P1 , then changes of P1 to P2 , and so on up to Plog2 n . Notice that we use an auxiliary matrix Z to compute the difference between the value of Pk before and after the update and that the computation of Z in line 6 always yields a Boolean matrix. 4.3 Analysis In the following, we discuss the correctness and the complexity of our implementation of operations Init∗ , Set∗ , Reset∗ , and Lookup∗ presented in Sect. 4.2. We recall that X is an n × n Boolean matrix and Pk , 0 ≤ k ≤ log2 n, are the polynomials introduced in Definition 14. Theorem 5 If at any time during a sequence of operations there is a path of length up to 2k between x and y in X, then Pk [x, y] = 1. Proof We proceed by induction. The base is trivial. We assume that the claim holds inductively for Pk−1 , and we show that, after any operation, the claim holds also for Pk . • Init∗ : since any Init∗ operation rebuilds from scratch Pk , the claim holds from Lemma 2. • Set∗ : let us assume that a Set∗ operation is performed on the i-th row and column of X and a new path π of length up to 2k , say π = x, . . . , i, . . . , y, appears in X due to this operation. We prove that Pk [x, y] = 1 after the operation. Observe that Pk .LazySet(Pk−1 , Pk−1 ) puts in place any new 1’s in any occurrence of the variable Pk−1 in data structure Pk . We remark that, although the maintained value of Pk in data structure Pk is not updated by LazySet and therefore the correctness of the current operation is not affected, this step is very important: indeed, new 1’s corresponding to new paths of length up to 2k−1 that appear in X will be useful in future Set∗ operations for detecting the appearance of new paths of length up to 2k . If both the portions x  i and i  y of π have length up to 2k−1 , then π gets 2 , and therefore in P , thanks to one of P .SetRow(i, P recorded in Pk−1 k k k−1 , Pk−1 ) or Pk .SetCol(i, Pk−1 , Pk−1 ). On the other hand, if i is close to (but does not coincide with) one endpoint of π , the appearance of π may be recorded

Algorithmica (2008) 51: 387–427

411

3 , but not in P 2 . This is the reason why degree 2 does not suffice for P in Pk−1 k k−1 in this dynamic setting. • Reset∗ : by inductive hypothesis, we assume that Pk−1 [x, y] flips to zero after a Reset∗ operation only if no path of length up to 2k−1 remains in X between x and y. Since any Pk .Reset operation on Pk leaves it as if cleared 1’s in Pk−1 were never set to 1, Pk [x, y] flips to zero only if no path of length up to 2k remains  in X.

We remark that the condition stated in Theorem 5 is only sufficient because Pk may keep track of paths having length strictly more than 2k , though no longer than 3k . However, for k = log2 n the condition is also necessary as no shortest path can be longer than n = 2k . Thus, it is straightforward to see that a path of any length between x and y exists at any time in X if and only if Plog2 n [x, y] = 1. The following theorem addresses the running time and space requirements of operations Init∗ , Set∗ and Reset∗ . Theorem 6 Lookup∗ requires O(n2 ) time and Set∗ requires O(n2 log n) time. Init∗ requires O(n3 log n) time, while Reset∗ requires O(1) time. All time bounds are amortized. The space required is O(n2 log n). Proof Since we maintain O(log n) polynomials with constant degree k = 3 and constant number of terms h = 3, it follows from Theorem 4 that the total space usage is O(n2 · log n). We now discuss the time bounds. Lookup∗ requires one Lookup operation, while Init∗ requires O(log n) Init operations. Set∗ can be implemented via O(log n) LazySet, SetRow, and SetCol operations. Thus their amortized bound follows from Theorem 4. As in Theorem 4, the cost of Reset∗ can be amortized  against previous Init∗ and Set∗ operations. Corollary 5 If we perform just one Init∗ operation in a sequence of (n) operations, or more generally one Init operation every (n) Reset operations, then the amortized cost of Init is O(n2 · log n) per operation. Corollary 6 If we perform just one Init∗ operation in a sequence of (n2 ) operations, or more generally one Init operation every (n2 ) Reset operations, and we perform no operations SetRow and SetCol, then the amortized cost of Init is O(n · log n) per operation. In the case where Init∗ is performed just once at the beginning of the sequence of operations, previous corollaries state that the amortized time per operation is at most O(n2 · log n). In the decremental case where only Reset∗ operations are performed, the amortized time per operation is at most O(n · log n). The algorithm that we have presented in this section can be viewed as a variant that features very different data structures of the fully dynamic transitive closure algorithm presented by King in [14]. King’s algorithm is based on a data structure for a graph G = (V , E) that maintains a logarithmic number of edge subsets E0 , . . . , Elog2 n with the property that

412

Algorithmica (2008) 51: 387–427

E0 = E and (x, y) ∈ Ei if there is a path x  y of length up to 2i in G. Moreover, if y is not reachable from x in G, then (x, y) ∈ Ei for all 0 ≤ i ≤ log2 n. The maintained values of our polynomials P0 , . . . , Plog2 n here correspond to the sets E0 , . . . , Elog2 n . The algorithm by King also maintains log2 n forests F0 , . . . , Flog2 n−1 such that Fi uses edges in Ei and includes 2n trees Outi (v) and Ini (v), two for each node v ∈ V , such that Outi (v) contains all nodes reachable from v using at most 2 edges in Ei , and Ini (v) contains all nodes that reach v using at most 2 edges in Ei . For each pair of nodes, also a table Counti is maintained, where Counti [x, y] is the number of nodes v such that x ∈ Ini (v) and y ∈ Outi (v). Now, Ei is maintained so as to contain edges (x, y) such that Counti−1 [x, y] > 0. Trees Ini (v) and Outi (v) are maintained for any node v by means of deletions-only data structures [6] which are rebuilt from scratch after each v-centered insertion of edges. Our data structures for polynomials over Boolean matrices Pi play the same role as King’s forests Fi of Ini and Outi trees and of counters Counti . While King’s data structures require O(n3 · log n) worst-case initialization time on dense graphs, the algebraic properties of Boolean matrices allow us to exploit fast matrix multiplication subroutines for initializing more efficiently our data structures in O(nω · log n) time in the worst case, where ω ≤ 2.38.

5 Transitive Closure Updates in O(n2 ) Time In this section we show a more powerful method for casting fully dynamic transitive closure into the problem of reevaluating polynomials over Boolean matrices presented in Sect. 3.1. This method hinges upon the well-known equivalence between transitive closure and matrix multiplication on a closed semiring and yields a new deterministic algorithm that improves the best known bounds for fully dynamic transitive closure. Our algorithm supports each update operation in O(n2 ) amortized time and answers each reachability query with just one matrix lookup. The space used is O(n2 ). 5.1 Data Structure Let X be a Boolean matrix and let X ∗ be its Kleene closure. Before discussing the dynamic case, we recall the main ideas behind the algorithm for computing statically X ∗ . Definition 15 Let Bn be the set of n × n Boolean matrices and let X ∈ Bn . Without loss of generality, we assume that n is a power of 2. Define a mapping F : Bn → Bn by means of the following equations: ⎧ E = (A + BD ∗ C)∗ , ⎪ ⎪ ⎨ F = EBD ∗ , (1) G = D ∗ CE, ⎪ ⎪ ⎩ ∗ ∗ ∗ H = D + D CEBD ,

Algorithmica (2008) 51: 387–427

413

where A, B, C, D and E, F, G, H are obtained by partitioning X and Y = F(X) into sub-matrices of dimension n2 × n2 as follows: A C

X=

B , D

Y=

E G

F . H

The following fact is well known [18]: if X is an n × n Boolean matrix, then F(X) = X ∗ . Another equivalent approach is given below: Definition 16 Let Bn be the set of n × n Boolean matrices, let X ∈ Bn and let G : Bn → Bn be the mapping defined by means of the following equations: ⎧ E = A∗ + A∗ BH CA∗ , ⎪ ⎪ ⎨ F = A∗ BH, G = H CA∗ , ⎪ ⎪ ⎩ H = (D + CA∗ B)∗ ,

(2)

where X and Y = G(X) are defined as: X=

A C

B , D

Y=

E G

F . H

It is easy to show that, for any X ∈ Bn , G(X) = F(X). Both F(X) and G(X) can be computed in O(nω ) worst-case time [18], where ω is the exponent of Boolean matrix multiplication. We now define another function H such that H(X) = X ∗ , based on a new set of equations obtained by combining (1) and (2). Our goal is to define H in such a way that it is well-suited for efficient reevaluation in a fully dynamic setting. Lemma 4 Let Bn be the set of n × n Boolean matrices, let X ∈ Bn and let H : Bn → Bn be the mapping defined by means of the following equations: ⎧ P = D∗, ⎪ ⎪ ⎪ ⎪ ⎨ E1 = (A + BP 2 C)∗ , F1 = E12 BP , ⎪ ⎪ ⎪ G1 = P CE12 , ⎪ ⎩ H1 = P CE12 BP ,

E2 = E1 BH22 CE1 , F2 = E1 BH22 , G2 = H22 CE1 , H2 = (D + CE12 B)∗ ,

where X and Y = H(X) are defined as: X=

A C

B , D

Then, for any X ∈ Bn , H(X) = X ∗ .

Y=

E G

F . H

E = E1 + E2 , F = F1 + F2 , G = G1 + G2 , H = H1 + H2 ,

(3)

414

Algorithmica (2008) 51: 387–427

Proof We prove that E1 + E2 , F1 + F2 , G1 + G2 and H1 + H2 are the sub-matrices of X ∗ : X∗ =

E1 + E2 G1 + G2

F1 + F2 . H1 + H2

We first observe that, by definition of Kleene closure, (X ∗ )2 = X ∗ . Thus, since E1 = (A + BP 2 C)∗ , H2 = (D + CE12 B)∗ and P = D ∗ are all closures, then we can replace E12 with E1 , H22 with H2 and P 2 with P . This implies that E1 = (A + BP C)∗ = (A + BD ∗ C)∗ and then E1 = E by 1. Now, E is a sub-matrix of X ∗ and encodes explicitly all paths in X with both end-points in V1 = {1, . . . , n2 }, and since E2 = EB(D + CEB)∗ CE, then E2 ⊆ E. It follows that E1 + E2 = E + E2 = E. With a similar argument, we can prove that F1 + F2 , G1 + G2 and H1 + H2 are sub-matrices of X ∗ . In particular, for H = H1 + H2 we also need to observe that  D ∗ ⊆ H2 . Note that the definition of H suggests a method for computing the Kleene closure of an n × n Boolean matrix, provided that we are able to compute Kleene closures of Boolean matrices of size n2 × n2 . The reason of using E12 , H22 and P 2 instead of E1 , H2 and P in (3) will be clear in Lemma 6 after presenting a fully dynamic version of the algorithm that defines H. In the next lemma we show that a divide and conquer algorithm that recursively uses H to solve sub-problems of smaller size requires asymptotically the same time as computing the product of two Boolean matrices. Theorem 7 Let X be an n × n Boolean matrix and let T (n) be the time required to compute recursively H(X). Then T (n) = O(nω ), where O(nω ) is the time required to multiply two Boolean matrices. Proof It is possible to compute E, F , G and H with three recursive calls of H, a constant number cm of multiplications, and a constant number cs of additions of n n 2 × 2 matrices. Thus: T (n) ≤ 3 T

n 2

+ cm M

n 2

+ cs

 n 2 2

,

where M(n) = O(nω ) is the time required to multiply two n × n Boolean matrices. Solving the recurrence relation, since log2 3 < max{ω, 2} = ω, we obtain that  T (n) = O(nω ). The previous theorem shows that, even if for H we need to compute one more closure than for F and for G, asymptotically we get the same running times. In the following, we study how to reevaluate efficiently H(X) = X ∗ under changes of X. Our data structure for maintaining the Kleene closure X ∗ is the following:

Algorithmica (2008) 51: 387–427

415

Data Structure 4 We maintain two n × n Boolean matrices X and Y decomposed in sub-matrices A, B, C, D, and E, F , G, H : X=

A C

B , D

Y=

E G

F . H

We also maintain the following 12 polynomials over n × n Boolean matrices with the data structure presented in Sect. 3.1.2: Q = A + BP 2 C, F1 = E12 BP , G1 = P CE12 , H1 = P CE12 BP ,

E2 = E1 BH22 CE1 , F2 = E1 BH22 , G2 = H22 CE1 , R = D + CE12 B,

E = E1 + E2 , F = F1 + F2 , G = G1 + G2 , H = H1 + H2

and we recursively maintain 3 Kleene closures P , E1 and H2 : P = D∗, with instances of size

n 2

×

n 2

E1 = Q∗ ,

H2 = R ∗

of Data Structure 2 presented in Sect. 3.1.2.

We note that Data Structure 4 is defined recursively: P , E1 and H2 are Kleene closures of n2 × n2 matrices. Also observe that the polynomials Q, F1 , G1 , H1 , E2 , F2 , G2 , R, E, F , G and H that we maintain have all constant degree ≤ 6. In Fig. 3 we show the acyclic graph of dependencies between portions of our data structure: there is an arc from node u to node v if the polynomial associated to u is a variable of the polynomial associated to v. For the sake of readability, we do not report nodes for the final polynomials E, F , G, H . A topological sort of this graph, e.g., τ = A, B, C, D, P , Q, E1 , R, H2 , F1 , G1 , H1 , E2 , F2 , G2 , E, F , G, H , yields a correct evaluation order for different portions of the data structure and thus gives a method for computing H(X). We remark that our data structure keeps memory of all the intermediate values produced when computing H(X) from scratch and maintains such values upon updates of X. As it was already observed in Sect. 4, maintaining intermediate results of some static algorithm for computing X ∗ is a fundamental idea for updating efficiently X ∗ whenever X gets modified. Since our data structure reflects the way H(X) is computed, it basically represents X ∗ as the sum of two Boolean matrices: the first, say X1∗ , is defined by submatrices E1 , F1 , G1 , H1 , and the second, say X2∗ , by submatrices E2 , F2 , G2 , H2 : X1∗ =

E1 G1

F1 , H1

X2∗ =

E2 G2

F2 . H2

In the next section we show how to implement operations Init∗ , Set∗ , Reset∗ and Lookup∗ introduced in Definition 5 in terms of operations Init, LazySet, SetRow and SetCol (see Sect. 3.1) on the polynomials of Data Structure 4.

416

Algorithmica (2008) 51: 387–427

Fig. 3 Data dependencies between polynomials and closures

5.2 Implementation of Operations From a high-level point of view, our approach is the following. We maintain X1∗ and X2∗ in tandem (see Fig. 4): whenever a Set∗ operation is performed on X, we update X ∗ by computing how either X1∗ or X2∗ are affected by this change. Such updates are lazily performed so that neither X1∗ nor X2∗ encode complete information about X ∗ , but their sum does. On the other side, Reset∗ operations update both X1∗ and X2∗ and leave the data structures as if any reset entry was never set to 1. We now describe in detail our implementation. To keep pseudocodes shorter and more readable, we assume that implicit Lookup and Lookup∗ operations are performed in order to retrieve the current value of data structures so as to use them in subsequent steps. Furthermore, we do not deal explicitly with base recursion steps.

Init∗ procedure Init∗ (Z) 1. begin 2. X←Z 3. P.Init∗ (D) 4. Q.Init(A, B, P , C) 5. E1 .Init∗ (Q) 6. R.Init(D, C, E1 , B)

Algorithmica (2008) 51: 387–427

417

Fig. 4 Overview of operations Init∗ , Set∗ and Reset∗

7. H2 .Init∗ (R) 8. F1 .Init(E1 , B, P ) 9. similarly for G1 , H1 , E2 , F2 , G2 , and then for E, F , G, H 10. end

Init∗ sets the initial value of X (line 2) and initializes the different portions of Data Structure 4 in accordance with the topological order τ of the graph of dependencies as explained in the previous subsection (lines 3–9).

Set∗ Before describing our implementation of Set∗ , we first define a useful shortcut for performing simultaneous SetRow and SetCol operations with the same i on more than one variable in a polynomial P : procedure P.Set(i, X1 , . . . , Xq ) 1. begin 2. P.SetRow(i, X1 , X1 ) 3. P.SetCol(i, X1 , X1 ) .. 4. . 5. P.SetRow(i, Xq , Xq ) 6. P.SetCol(i, Xq , Xq ) 7. end

Similarly, we give a shortcut2 for performing simultaneous LazySet operations on more than one variable in a polynomial P : 2 For the sake of simplicity, we use the same identifier LazySet for both the shortcut and the native operation on polynomials, assuming to use the shortcut in defining Set∗ .

418

Algorithmica (2008) 51: 387–427

procedure P.LazySet(X1 , . . . , Xq ) 1. begin 2. P.LazySet(X1 , X1 ) .. 3. . 4. P.LazySet(Xq , Xq ) 5. end

We also define an auxiliary operation LazySet∗ on closures that performs LazySet operations for variables A, B, C and D on the polynomials Q, R, F1 , G1 , H1 , E2 , F2 , and G2 and recurses on the closure P which depend directly on them. We assume that, if M is a variable of a polynomial maintained in our data structure, M = Mcurr − Mold is the difference between the current value Mcurr of M (i.e., the value after an update) and the old value Mold of M (i.e., the value before the update). procedure LazySet∗ (X) 1. begin 2. X ← X + X 3. Q.LazySet(A, B, C) 4. R.LazySet(B, C, D) 5. similarly for F1 , G1 , H1 , E2 , F2 , and G2 6. P.LazySet∗ (D) 7. end

Using the shortcuts Set and LazySet and the new operation LazySet∗ , we are now ready to define Set∗ . procedure Set∗ (i, X) 1. begin 2. X ← X + IX,i + JX,i 3. if 1 ≤ i ≤ n2 then 4. Q.Set(i, A, B, C) 5. E1 .Set∗ (i, Q) 6. F1 .Set(i, E1 , B) 7. G1 .Set(i, C, E1 ) 8. H1 .Set(i, C, E1 , B) 9. R.Set(i, C, E1 , B) 10. H2 .LazySet∗ (R) 11. G2 .LazySet(C, E1 ) 12. F2 .LazySet(E1 , B) 13. E2 .LazySet(E1 , B, C) 14. else 15. i ← i − n2 16. P.Set∗ (i, D) 17. R.Set(i, B, C, D) 18. H2 .Set∗ (i, R) 19. G2 .Set(i, H2 , C) 20. F2 .Set(i, B, H2 ) 21. E2 .Set(i, B, H2 , C) 22. Q.Set(i, B, P , C)

{ n2 + 1 ≤ i ≤ n }

Algorithmica (2008) 51: 387–427

419

23. E1 .LazySet∗ (Q) 24. F1 .LazySet(B, P ) 25. G1 .LazySet(P , C) 26. H1 .LazySet(B, P , C) 27. E.Init(E1 , E2 ) 28. F.Init(F1 , F2 ) 29. G.Init(G1 , G2 ) 30. H.Init(H1 , H2 ) 31. end

Set∗ performs an i-centered update in X and runs through the closures and the polynomials of Data Structure 4 to propagate any changes of A, B, C, D to E, F , G, H . The propagation order is Q, E1 , F1 , G1 , H1 , R, H2 , G2 , F2 , E2 , E, F , G, H  if 1 ≤ i ≤ n2 and P , R, H2 , G2 , F2 , E2 , Q, E1 , F1 , G1 , H1  if n2 + 1 ≤ i ≤ n and is defined in accordance with a topological sort of the graph of dependencies between portions in Data Structure 4 shown in Fig. 3. Roughly speaking, Set∗ updates the different portions of the data structure in accordance with the value of i as follows: 1. If 1 ≤ i ≤ n2 , Set∗ fully updates Q, R, E1 , F1 , G1 , H1 (lines 4–9) and lazily updates E2 , F2 , G2 , H2 (lines 10–13). See Fig. 5(a). 2. If n2 + 1 ≤ i ≤ n, Set∗ fully updates P , Q, R, E2 , F2 , G2 , H2 (lines 16–22) and lazily updates E1 , F1 , G1 , H1 (lines 23–26). See Fig. 5(b). We highlight that it is not always possible to perform efficiently full updates of all the portions of Data Structure 4. Actually, some matrices may change everywhere, and not only in a row or in a column. Such unstructured changes imply that we can only perform lazy updates on such data structures, as they cannot be efficiently handled by means of i-centered SetRow and SetCol operations. We now explain in detail the operations performed by Set∗ in accordance with the two cases 1 ≤ i ≤ n2 and n2 + 1 ≤ i ≤ n. Case 1: 1 ≤ i ≤ n2 . In this case an i-centered update of X may affect the i-th row and the i-th column of A, the i-th row of B and the i-th column of C, while D is not affected at all by this kind of update (see Fig. 4). The operations performed by Set∗ when 1 ≤ i ≤ n2 are therefore the following: Line 2: an i-centered set operation is performed on X. Line 4: Q = A + BP 2 C is updated by performing SetRow and SetCol operations for any variables A, B and C being changed. P = D ∗ does not change since, as already observed, D is not affected by the change. Notice that new 1’s may appear in Q only in the i-th row and column due to this operation. Line 5: Set∗ is recursively called to propagate the changes of Q to E1 . We remark that E1 may change also outside of the i-th row and column due to this operation. Nevertheless, as we will see in Lemma 5, the fact that E1 is a closure implies that new 1’s appear in a very structured way. This makes it possible to propagate changes efficiently to any polynomial that, in turn, depends on E1 . Lines 6–9: polynomials F1 , G1 , H1 and R are updated by performing SetRow and SetCol operations for any variables E1 , B and C being changed. We recall

420

Algorithmica (2008) 51: 387–427

Fig. 5 Portions of Data Structure 4 affected during a Set∗ operation when: (a) 1 ≤ i ≤ n2 ; (b) n2 + 1 ≤ i ≤ n

that such operations take into account only the entries of E1 lying in the i-th row and in the i-th column, albeit other entries may be non-zero. Again, Lemma 5 and Lemma 6 will show that this is sufficient. Lines 10–13: H2 = R ∗ is not updated, but new 1’s that appear in R are lazily inserted in the data structure of H2 by calling LazySet∗ . Then LazySet operations are carried out on polynomials G2 , F2 , E2 to insert in the data structures that maintain them any new 1’s that appear in C, E1 and B. Lines 27–30. Recompute polynomials E, F , G and H from scratch. This is required as F1 , G1 and H2 may change everywhere and not only in a row and a column. Differently from the case of E1 , whose change is structured as it is a closure, we cannot exploit any particular structure of F1 , G1 and H2 that allows us to use SetRow and SetCol and we are forced to use Init. Note that, since E, F , G and H are all polynomials of degree 1, this is not a bottleneck in terms of running time.

Algorithmica (2008) 51: 387–427

421

Case 2: n2 + 1 ≤ i ≤ n. In this case an i-centered update of X may affect only the i-th row and the i-th column of D, the i-th row of C and the i-th column of B, while A is not affected at all by this kind of update (see Fig. 4). Operations performed by Set∗ are completely analogous to the case 1 ≤ i ≤ n2 , except for the fact that we need to rescale the index i in line 15 and we have also to perform a recursive call to update P in line 16.

Reset∗ Before describing our implementation of Reset∗ , we define a useful shortcut3 for performing simultaneous Reset operations on more than one variable in a polynomial P . procedure P.Reset(X1 , . . . , Xq ) 1. begin 2. P.Reset(X1 , X1 ) .. 3. . 4. P.Reset(Xq , Xq ) 5. end

Using this shortcut, we are now ready to define Reset∗ . We assume that, if M is a variable of a polynomial maintained in our data structure, M = Mold − Mcurr is the difference between the value Mold of M just before calling Reset∗ and the current value Mcurr of M. procedure Reset∗ (X) 1. begin 2. X ← X − X 3. P.Reset∗ (D) 4. Q.Reset(A, B, P , C) 5. E1 .Reset∗ (Q) 6. R.Reset(D, C, E1 , B) 7. H2 .Reset∗ (R) 8. F1 .Reset(E1 , B, P ) 9. { similarly for G1 , H1 , E2 , F2 , G2 , and then for E, F , G, H } 10. end

Reset∗ resets any entries of X as specified by X and runs through the closures and the polynomials in the data structure to propagate any changes of A, B, C, D to E, F , G, H . The propagation is done in accordance with a topological order τ of the graph of dependencies shown in Fig. 3 and is the same order followed by Init∗ , which has a similar structure. Actually, we could think of Reset∗ as a function that “undoes” any previous work performed by Init∗ and Set∗ on the data structure, leaving it as if the reset entries of X were never set to 1. 3 For the sake of simplicity, we use the same identifier Reset for both the shortcut and the native operation on polynomials, assuming to use the shortcut in defining Reset∗ .

422

Algorithmica (2008) 51: 387–427

Lookup∗ procedure Lookup∗ (x, y) 1. begin 2. return Y [x, y] 3. end

Lookup∗ simply returns the maintained value of Y [x, y]. 5.3 Analysis Now we discuss the correctness and the complexity of our implementation. Before providing the main claims, we give some preliminary definitions and lemmas that are useful for capturing algebraic properties of the changes that polynomials in our data structure undergo during a Set∗ operation. The next definition recalls a property of Boolean update matrices that is related to the operational concept of i-centered update. Definition 17 We say that a Boolean update matrix X is i-centered if X = IX,i + JX,i , i.e., all entries lying outside the i-th row and the i-th column are zero. If the variation X of some matrix X during an update operation is i-centered and X is a variable of a polynomial P that has to be efficiently reevaluated, then we can use P.SetRow and P.SetCol operations which are especially designed for doing so. But what happens if X changes by a X that is not i-centered? Can we still update efficiently the polynomial P without recomputing it from scratch via Init? This is the case of E1 and E1 while performing a Set∗ update with 1 ≤ i ≤ n2 . In the following we show that, under certain assumptions on X and X (which are satisfied by E1 and E1 ), we can still solve the problem efficiently. While the property of being i-centered is related to an update matrix by itself, the following two definitions are concerned with properties of an update matrix X with respect to the matrix X to which it is applied: Definition 18 If X is a Boolean matrix and X is a Boolean update matrix, we say that X is i-transitive with respect to X if IX,i = IX,i · X and JX,i = X · JX,i . Definition 19 If X is a Boolean matrix and X is a Boolean update matrix, we say that X is i-complete with respect to X if X = JX,i · IX,i + X · IX,i + JX,i · X. Intuitively, if X is i-transitive and i-complete with respect to X, then X encodes all the relevant information in its i-th row and column. As we will see, in this case, X can be updated efficiently by means of i-centered operations (SetRow and SetCol). Using the previous definitions we can show that the variation of X ∗ due to an i-centered update of X is i-transitive and i-complete with respect to X. Lemma 5 Let X be a Boolean matrix and let X be an i-centered update matrix. If we denote by X ∗ the matrix (X + X)∗ − X ∗ , then X ∗ is i-transitive and i-complete with respect to X ∗ .

Algorithmica (2008) 51: 387–427

423

Proof The following equalities prove the first condition of i-transitivity: IX∗ ,i · X ∗ = I(X+X)∗ −X∗ ,i · X ∗ = I(X+X)∗ ·X∗ −X∗ ·X∗ ,i = I(X+X)∗ −X∗ ,i = IX∗ ,i . The other conditions can be proved analogously. The hypothesis that X is icentered is necessary for the i-completeness.  The following lemma shows under what conditions for X and X it is possible to perform operations of the kind X ← X + X on a variable X of a polynomial by reducing such operations to i-centered updates even if X is not i-centered. Lemma 6 If X is a Boolean matrix such that X = X ∗ and X is an i-transitive and i-complete update matrix with respect to X, then X + X = (X + IX,i + JX,i )2 . Proof Since X = X ∗ it holds that X = X 2 and X = X + IX,i · JX,i . The proof 2 follows from Definition 18 and Definition 19 and from the facts that: IX,i ⊆ IX,i , 2  JX,i ⊆ JX,i and X = X + IX,i + JX,i . It follows that, under the hypotheses of Lemma 6, if we replace any occurrence of X in P with X 2 and we perform both P.SetRow(i, IX,i , X) and P.SetCol(i, JX,i , X), then new 1’s in P correctly appear. This is the reason why in Data Structure 4 we used E12 , H22 , and P 2 instead of E1 , H2 , and P , respectively. Before stating the main theorem of this section, which addresses the correctness of operations on our data structure, we discuss a general property of polynomials and closures over Boolean matrices that will be useful in proving the theorem. Definition 20 Let Bn be the set of n × n Boolean matrices and let P : Bn → Bn be a  is a relaxation of P if P (X) ⊆ P (X) function over Boolean matrices. We say that P for any X ∈ Bn .  Lemma 7 Let P and Q be polynomials or closures over Boolean matrices and let P  and Q be relaxations of P and Q, respectively. Then, for any X ∈ Bn : P (X)) ⊆ Q(P (X)). Q( = P (X) and Y = P (X). By definition, we have: Y  ⊆ Y and Q( Y ) ⊆ Proof Let Y ). By exploiting a monotonic behavior of polynomials and closures over Boolean Q(Y  ⊆ Y ⇒ Q(Y ) ⊆ Q(Y ). Thus: Q( Y ) ⊆ Q(Y ) ⊆ Q(Y ) ⇒ matrices, we have: Y Y ) ⊆ Q(Y ) ⇒ Q( P (X)) ⊆ Q(P (X)). Q(  Theorem 8 Operations Set∗ , LazySet∗ , Reset∗ and Init∗ maintain a matrix Y such that Y = X ∗ . Proof The proof is by induction on the size n of matrices in Data Structure 4. The base is trivial. We assume that the claim holds for instances of size n2 and we prove that it holds also for instances of size n.

424

Algorithmica (2008) 51: 387–427

• Init∗ : since Init∗ performs Init operations on each polynomial in Data Structure 4, then Y = H(X) = X ∗ . • Set∗ : we first prove that Y ⊆ X ∗ . Observe that Y is obtained as a result of a composition of functions that relax the correct intermediate values of polynomials and closures of Boolean matrices in our data structure allowing them to contain less 1’s. Indeed, by the properties of Lookup described in Sect. 3.1, we know that, if P is the correct value of a polynomial at any time, then P.Lookup() ⊆ P . Similarly, by inductive hypothesis, if K is a Kleene closure of an n2 × n2 Boolean matrix, then at any time K.Lookup∗ (x, y) = 1 ⇒ K[x, y] = 1. The claim then follows by Lemma 7, which states that the composition of relaxations computes values containing at most the 1’s contained in the values computed by the correct functions. To prove that X ∗ ⊆ Y , it suffices to prove that if H[x, y] flips from 0 to 1 due to operation Set∗ , then either X1∗ [x, y] flips from 0 to 1 (due to lines 4–8 when 1 ≤ i ≤ n2 ), or X2∗ [x, y] flips from 0 to 1 (due to lines 17–21 when n2 + 1 ≤ i ≤ n). Without loss of generality, assume that the Set∗ operation is performed with 1 ≤ i ≤ n2 (the proof is completely analogous if n2 + 1 ≤ i ≤ n). As shown in Fig. 4, sub-matrices A, B and C may undergo i-centered updates due to this operation and so their variation can be correctly propagated through SetRow and SetCol operations to polynomial Q (line 4) and to polynomials F1 , G1 and H1 (lines 6–8). As Q is also i-centered due to line 4, any variation of Q, that is assumed to be elsewhere correct from previous operations, can be propagated to closure E1 through a recursive call of Set∗ in line 5. By the inductive hypothesis, this propagation correctly reveals any new 1’s in E1 . We remark that E1 may contain less 1’s than E due to any previous LazySet operations done in line 23. Observe now that E1 occurs in polynomials F1 , G1 and H1 and that E1 is not necessarily i-centered. This would imply that we cannot propagate directly changes of E1 to these polynomials, as no efficient operation for doing so was defined in Sect. 3.1. However, by Lemma 5, E1 is i-transitive and i-complete with respect to E1 . Since E1 = E1∗ , by Lemma 6 performing both SetRow(i, IE1 ,i , E1 ) and SetCol(i, JE1 ,i , E1 ) operations on data structures F1 , G1 and H1 in lines 6–8 is sufficient to correctly reveal new 1’s in F1 , G1 and H1 . Again, note that F1 , G1 and H1 may contain less 1’s than F , G and H , respectively, due to any previous LazySet operations done in lines 23–26. We have then proved that lines 4–8 correctly propagate any i-centered update of X to X1∗ . To conclude the proof, we observe that E1 also occurs in polynomials E2 , F2 , G2 , R and indirectly affects H2 . Unfortunately, we cannot update H2 efficiently as R is neither i-centered, nor i-transitive/i-complete with respect to R. So in lines 9–13 we limit ourselves to update explicitly R and to log any changes of E1 by performing LazySet operations on polynomials G2 , F2 , and E2 and a LazySet∗ operation on H2 . This is sufficient to guarantee the correctness of subsequent Set∗ operations for n2 + 1 ≤ i ≤ n. • Reset∗ : this operation runs in judicious order through the different portions of the data structure and undoes the effects of previous Set∗ and Init∗ operations. Thus, any property satisfied by Y still holds after performing a Reset∗ operation.

Algorithmica (2008) 51: 387–427

425

 To conclude this section, we address the running time of operations and the space required to maintain an instance of our data structure. Theorem 9 Lookup∗ and Set∗ require O(n2 ) time. Init∗ requires O(n3 ) time, while Reset∗ requires O(1) time. All time bounds are amortized. The space required is O(n2 ). Proof Since all the polynomials in Data Structure 4 are of constant degree and involve a constant number of terms, the amortized cost of any SetRow, SetCol, and LazySet operation on them is quadratic in n2 (see Theorem 4). Let T (n) be the amortized time complexity of any Set∗ and LazySet∗ operation. Then: T (n) ≤ 3 T

n 2

+

c n2 4

for some suitably chosen constant c > 0. This implies that T (n) = O(n2 ). Init∗ recomputes recursively H from scratch using Init operations on polynomials. Let TInit∗ (n) be the amortized cost of Init∗ . Since the amortized cost of Init is O(n3 ) by Theorem 4, then TInit∗ (n) satisfies the relation: TInit∗ (n) ≤ 3 TInit∗

n 2

+

c n3 . 4

This implies that TInit∗ (n) = O(n3 ). As in Theorem 4, the cost of Reset∗ can be amortized against previous Init∗ and Set∗ operations. To conclude the proof, observe that if K(n) is the space used to maintain Data Structure 4, and M(n) is the space required to maintain a polynomial with the data structure of Sect. 3.1, then: n + 12 M(n). K(n) ≤ 3 K 2 Since M(n) = O(n2 ) by Theorem 4, then K(n) = O(n2 ).



Corollary 7 If we perform just one Init∗ operation in a sequence of (n) operations, or more generally one Init∗ operation every (n) Reset∗ operations, then the amortized cost of Init∗ is O(n2 ) per operation. Corollary 8 If we perform just one Init∗ operation in a sequence of (n2 ) operations, or more generally one Init∗ operation every (n2 ) Reset∗ operations, and we perform no Set∗ operations, then the amortized cost of Init∗ is O(n) per operation. In the traditional case where Op1 =Init∗ and Opi =Init∗ for any i > 1, i.e., Init∗ is performed just once at the beginning of the sequence of operations, previous corollaries state that the amortized time per operation is at most O(n2 ). In the

426

Algorithmica (2008) 51: 387–427

decremental case where only Reset∗ operations are performed, the amortized time per operation is at most O(n).

6 Conclusions In this paper we have presented new and improved algorithms for maintaining the transitive closure of a directed graph under edge insertions and edge deletions. As a main contribution, we have introduced a general framework for casting fully dynamic transitive closure into the problem of dynamically reevaluating polynomials over matrices when updates of variables are performed. Such technique has turned out to be very flexible and powerful, leading both to revisit the best known algorithm for fully dynamic transitive closure [14] from a completely different perspective, and to design new and faster algorithms for this problem. We remark that the algebraic approach to dynamic data structures for graph problems introduced in this paper has sparkled a wealthy of new results [17, 20, 21]. Designing efficient dynamic data structures for maintaining polynomials over Boolean matrices allowed us to devise the fairly complex deterministic algorithm described in Sect. 5, which supports updates in quadratic amortized time and queries with just one matrix lookup. Our algorithm improves the best bounds for fully dynamic transitive closure achieved in [14]. After our work, Roditty [19] reduced the Init bound from O(n3 ) to O(mn), and Sankowski [20] made our bounds worst case. Acknowledgements We are indebted to Garry Sagert and Mikkel Thorup for enlightening discussions, to Valerie King for insights on this work, and to the anonymous referees for useful comments and suggestions.

References 1. Aho, A.V., Hopcroft, J.E., Ullman, J.D.: The Design and Analysis of Computer Algorithms. AddisonWesley, Reading (1974) 2. Coppersmith, D., Winograd, S.: Matrix multiplication via arithmetic progressions. J. Symb. Comput. 9, 251–280 (1990) 3. Cormen, T.H., Leiserson, C.E., Rivest, R.L., Stein, C.: Introduction to Algorithms. MIT Press, Cambridge (2001) 4. Demetrescu, C., Italiano, G.F.: Fully dynamic transitive closure: breaking through the O(n2 ) barrier. In: Proc. of the 41st IEEE Annual Symposium on Foundations of Computer Science (FOCS’00), pp. 381–389 (2000) 5. Demetrescu, C., Italiano, G.F.: Trade-offs for fully dynamic reachability on dags: breaking through the O(n2 ) barrier. J. Assoc. Comput. Mach. (J. ACM) 52(2), 147–156 (2005) 6. Even, S., Shiloach, Y.: An on-line edge-deletion problem. J. ACM 28, 1–4 (1981) 7. Fischer, M.J., Meyer, A.R.: Boolean matrix multiplication and transitive closure. In: Conference Record 1971 Twelfth Annual Symposium on Switching and Automata Theory, East Lansing, Michigan, 13–15 October 1971, pp. 129–131. IEEE 8. Furman, M.E.: Application of a method of fast multiplication of matrices in the problem of finding the transitive closure of a graph. Sov. Math. Dokl. 11(5) (1970). English translation 9. Henzinger, M., King, V.: Fully dynamic biconnectivity and transitive closure. In: Proc. 36th IEEE Symposium on Foundations of Computer Science (FOCS’95), pp. 664–672 (1995) 10. Ibaraki, T., Katoh, N.: On-line computation of transitive closure for graphs. Inf. Process. Lett. 16, 95–97 (1983)

Algorithmica (2008) 51: 387–427

427

11. Italiano, G.F.: Amortized efficiency of a path retrieval data structure. Theor. Comput. Sci. 48(2–3), 273–281 (1986) 12. Italiano, G.F.: Finding paths and deleting edges in directed acyclic graphs. Inf. Process. Lett. 28, 5–11 (1988) 13. Khanna, S., Motwani, R., Wilson, R.H.: On certificates and lookahead in dynamic graph problems. Algorithmica 21(4), 377–394 (1998) 14. King, V.: Fully dynamic algorithms for maintaining all-pairs shortest paths and transitive closure in digraphs. In: Proc. 40th IEEE Symposium on Foundations of Computer Science (FOCS’99), pp. 81– 99 (1999) 15. King, V., Sagert, G.: A fully dynamic algorithm for maintaining the transitive closure. J. Comput. Syst. Sci. 65(1), 150–167 (2002) 16. La Poutré, J.A., van Leeuwen, J.: Maintenance of transitive closure and transitive reduction of graphs. In: Proc. Workshop on Graph-Theoretic Concepts in Computer Science. Lecture Notes in Computer Science, vol. 314, pp. 106–120. Springer, Berlin (1988) 17. Mucha, M., Sankowski, P.: Maximum matchings via Gaussian elimination. In: Proceedings of the 45th Annual IEEE Symposium on Foundations of Computer Science (FOCS’04), pp. 248–255 (2004) 18. Munro, I.: Efficient determination of the transitive closure of a directed graph. Inf. Process. Lett. 1(2), 56–58 (1971) 19. Roditty, L.: A faster and simpler fully dynamic transitive closure. In: Proceedings of the 14th ACMSIAM Symposium on Discrete Algorithms (SODA), Baltimore, MD, USA, pp. 404–412 (2003) 20. Sankowski, P.: Dynamic transitive closure via dynamic matrix inverse. In: Proceedings of the 45th Annual IEEE Symposium on Foundations of Computer Science (FOCS’04), Rome, Italy, pp. 509– 517 (2004) 21. Sankowski, P.: Subquadratic algorithm for dynamic shortest distances. In: Proceedings of the 11th Annual International Conference on Computing and Combinatorics (COCOON’05), pp. 461–470 (2005) 22. Yellin, D.M.: Speeding up dynamic transitive closure for bounded degree graphs. Acta Inf. 30, 369– 384 (1993)