Algorithmica (2010) 56: 180–197 DOI 10.1007/s00453-008-9166-2

Fast Dynamic Transitive Closure with Lookahead Piotr Sankowski · Marcin Mucha

Received: 27 January 2006 / Accepted: 4 January 2008 / Published online: 25 January 2008 © Springer Science+Business Media, LLC 2008

Abstract In this paper we consider the problem of dynamic transitive closure with lookahead. We present a randomized one-sided error algorithm with updates and queries in O(nω(1,1,)− ) time given a lookahead of n operations, where ω(1, 1, ) is the exponent of multiplication of n × n matrix by n × n matrix. For  ≤ 0.294 we obtain an algorithm with queries and updates in O(n2− ) time, whereas for  = 1 the time is O(nω−1 ). This is essentially optimal as it implies an O(nω ) algorithm for boolean matrix multiplication. We also consider the offline transitive closure in ω planar graphs. For this problem, we show an algorithm that requires O(n 2 ) time to 1 process n 2 operations. We also show a modification of these algorithms that gives faster amortized queries. Finally, we give faster algorithms for restricted type of updates, so called element updates. All of the presented algorithms are randomized with one-sided error. All our algorithms are based on dynamic algorithms with lookahead for matrix inverse, which are of independent interest. Keywords Transitive closure · Dynamic algorithms · Lookahead · Algebraic algorithms · Matrix inverse · Planar graphs

1 Introduction In this paper we consider the problem of dynamic transitive closure with lookahead. A dynamic graph algorithm maintains a given property of a graph subject to dynamic P. Sankowski () · M. Mucha Institute of Informatics, Warsaw University, ul. Banacha 2, 02-097 Warsaw, Poland e-mail: [email protected] M. Mucha e-mail: [email protected]

Algorithmica (2010) 56: 180–197

181

changes, such as edge insertions and deletions. Dynamic algorithms with lookahead are provided with some kind of information about future operations. A dynamic algorithm for the transitive closure problem supports the following operations on a directed graph G = (V , E): insert(u, v): inserts an edge from u to v, delete(u, v): deletes an edge from u to v, query(u, v): answers “yes” if v is reachable from u, and “no” otherwise. We also consider a more general transitive closure problem with vertex operations. An algorithm for this problem supports the following operations: insert(Ev ): inserts a set Ev of outgoing edges of v, delete(Ev ): deletes a set Ev of incoming edges of v, query(v): returns the set of vertices reachable from v. The dynamic transitive closure problem is one of the most fundamental problems in graph algorithms and has numerous application. In many of these applications, a natural lookahead is available (see [8, 19]), which allows more efficient solutions. In this paper we study the relationship between the size of the lookahead and the complexity of the transitive closure problem. 1.1 Previous Work The problem of fully dynamic transitive closure was first solved in the case of planar graphs. The algorithm of Subramanian [18] is based on the separator decomposition ˜ 23 ) time. The first fully dynamic transitive and supports queries and updates in O(n closure algorithm in general graphs is due to Henzinger and King [6]. They used the data structure of Even and Shiloach [5] for decrementally maintaining a breadth-first search tree together with fast matrix multiplication for stitching together short paths. 0.58 ) amor˜ Their algorithm is Monte Carlo with one-sided error and requires O(nm tized time per update, and (n/ log n) time per query. Khanna, Motwani and Wilson [8] gave a deterministic algorithm with update time (n2.18 ) when a lookahead of (n0.18 ) updates is given. The next improvement was made by King and Sagert [10], who developed an algorithm for general directed graphs, supporting queries in constant time and updates in O(n2.26 ) time. They also showed an algorithm for acyclic graphs with O(n2 ) update time. These bounds were improved by King [9], who presented a deterministic algorithm with O(n2 log n) amortized update time and O(1) query time. In [4] Demetrescu and Italiano further improved these bounds by presenting an algorithm for general digraphs with O(n2 ) amortized update time and O(1) query time. They introduced a data structure for dynamically maintaining polynomials over matrices and used it for computing the transitive closure. A similar result was later presented in [12]. In [4] the authors also gave the first algorithm with subquadratic update time for directed acyclic graphs. This algorithm can answer queries in O(n ) time and perform updates in O(nω(1,,1)− + n1+ ) time, for any  ∈ [0, 1]. The current best bounds on ω(1, , 1) [7] imply an O(n0.575 ) query time and O(n1.575 ) update time. This subquadratic algorithm is randomized with one-sided error. More recently Roditty and Zwick [13] presented an algorithm with

182

Algorithmica (2010) 56: 180–197

√ √ O(m n) update time and O( n) query time and another [14] with O(m + n log n) update time and O(n) query time. The first worst-case O(n2 ) algorithm was shown recently by Sankowski [15]. He also gives two subquadratic algorithms for this problem: one with O(nω(1,,1)− + n1+ ) update time and O(n ) query time, and another with O(nω(1,,1)− + n2 ) update time and O(n2 ) query time. In particular, using current best bounds on ω(1, , 1) this gives an algorithm with O(n1.575 ) update time and O(n0.575 ) query time and another with O(n1.495 ) time for updates and queries. 1.2 Our Results In this paper we present new algorithms for the dynamic transitive closure problem with lookahead. We use the reduction of this problem to the problem of dynamic matrix inverse, introduced by Sankowski [15]. An algorithm for this problem maintains a dynamic matrix A subject to the following operations: update(i, j, x): changes Ai,j to x, query(i, j ): returns (A−1 )i,j . We also consider a dynamic matrix inverse problem with column operations, where the following operations are supported: update(i, vi ): changes (A)i — the i-th column of A to vi , query(i): returns (A−1 )i . We give two algorithms for the dynamic matrix inverse problem with column operations. The first algorithm supports updates and queries in O(nω(1,1,)− ) time, given a lookahead of n operations. The second algorithm requires a lookahead of a block of operations containing n updates and gives a trade-off between amortized update ˆ ˜ ω(1,1,)− ) time and query complexities. In particular it can support updates in O(n ω(1,1,)−1 ˆ ˜ and queries in O(n ). Here, the ωˆ is the upper bound on ω given by Huang and Pan [7] and not the actual value of ω (see Sect. 3). We also show that we can achieve better bounds if we allow only element updates. The query and update complexity is then O(nω(1,1,μ)−μ ), where μ satisfies ω(1, 1, μ) − μ = ω(1, μ, ) − . Our results are summarized on Fig. 1. A matrix is planar if its non-zero structure corresponds to a planar graph. We give an algorithm for planar case of dynamic matrix inverse problem which requires ω−1 1 O(n 2 ) amortized time to process n 2 operations. These algorithms for the dynamic matrix inverse problem yield randomized algorithms for the transitive closure problem with the same update and query complexities. Column operations translate to vertex operations, and element operations translate to edge operations. The only known result on dynamic transitive closure algorithms with lookahead is the one by Khanna, Motwani and Wilson [8], which has already been superseded by dynamic algorithms without lookahead. Thus, we can only compare the time bounds of our algorithms with algorithms without lookahead, which is obviously a bit unfair. All known algorithms with subquadratic query and update times work only for edge operations (see [4, 15]). We achieve subquadratic times even for vertex operations, provided that we are given a lookahead of (n ) operations, for any  > 0.

Algorithmica (2010) 56: 180–197

183

Fig. 1 The dependence of the time complexity nx on the size of lookahead available to the algorithm n

The best known algorithms for vertex operations are the algorithm of Demetrescu and Italiano [4] and the algorithm of Sankowski [15]. Both support O(n2 ) updates and O(n) queries, which is a special case (no lookahead) of one of our algorithms. Moreover, in case of lookahead of n operations our algorithms are essentially optimal with respect to the O(nω ) static time algorithm. We are not aware of any result of this kind for dynamic algorithms without lookahead for the transitive closure problem. By optimality we mean that using our algorithms the static transitive closure problem can be solved within the best known time bounds. In order to compute the static transitive closure we start with an empty graph, perform n vertex updates to get a given graph G, and then use n vertex queries to retrieve its transitive closure, which gives O(nnω−1 ) = O(nω ) time. For edge operations our algorithms improve over the algorithms of Sankowski [15] for a lookahead of (n ) operations, for any  > 0. As for our algorithm for planar graphs, it answers updates and queries in ω−1 O(n 2 ) = O(n0.69 ) time, so it is a little bit slower than the algorithm of Subra˜ 23 ) time, and without a lookahead. Still, we manian [18] which does the same in O(n believe that there is room for improvement in our approach. Also, even very small improvement of ω would lead to a faster algorithm.

2 Dynamic Transitive Closure Our algorithms for transitive closure use the reduction presented in [15]. The main idea is to construct a certain adjacency matrix, so that the inverse of this matrix encodes with high probability the transitive closure. Then, we can directly apply any algorithm for the dynamic matrix inverse for computing the transitive closure. This reduction is stated in the following theorem. Theorem 1 If there exists an algorithm for dynamic matrix inverse in non-singular case with update cost of O(nα ) arithmetic operations and query cost of O(nβ ) arithmetic operations, then there exists an algorithm for dynamic transitive closure prob-

184

Algorithmica (2010) 56: 180–197

lem with updates in O(nα ) time and queries in O(nβ ) time. The algorithm is randomized with one-sided error. This reduction also works for algorithms with lookahead. Importantly, in order to get the transitive closure algorithms we only need to device algorithms for dynamically computing the inverse matrix only in the simpler non-singular case. In the case of the algorithms with no lookahead the singular case can be easily reduced to non-singular case [16]. However, in the case of lookahead such an reduction in not known. Moreover, one should note that by applying the above theorem we always get an randomized algorithm even if the base algorithm for the inverse matrix was deterministic. In the case of planar graphs we need a stronger theorem. In Sect. 7 we give an offline algorithm for planar matrices. It can be used to solve an offline transitive closure problem in planar graphs. The algorithm requires certain assumptions about the matrix it works on. In Sect. 8 we show that these assumptions hold in our case.

3 Fast Multiplication of Rectangular Matrices In the construction of our algorithms we use algorithms for fast rectangular matrix multiplication. These algorithms were developed by Huang and Pan [7] and are based on the previous results of Coppersmith and Winograd [3] and Coppersmith [2]. Let us denote by ω(a, b, c) the exponent of the multiplication of an na × nb matrix by an nb ×nc matrix. For a = b = c = 1 we get the exponent of square matrix multiplication ω = ω(1, 1, 1). The best bound on ω is currently ω < 2.376 [3]. Coppersmith [2] showed that it is possible to multiply a n × n0.294 matrix by a n0.294 × n matrix in ˜ 2 ) arithmetic operations. O(n Let α = sup{a : ω(1, a, 1) = 2 + o(1)} ≤ 0.294. By combining bounds on α and ω Huang and Pan [7] showed the following bound on ω(a, b, c). Lemma 1 (Huang and Pan [7]) Let 1 ≥ b ≥ c, then  1+b ω(1, b, c) ≤ ω(1, ˆ b, c) = β + γ b + δc where β = 1, γ =

1+α−αω 1−α

≤ 0.844, δ =

ω−2 1−α

if 0 ≤ c ≤ αb, if αb ≤ c ≤ 1,

≤ 0.533.

Note that we introduced ωˆ to denote the Huang-Pan bound on ω. This will be useful in Sect. 5. We also use the fact, shown by Bunch and Hopcroft [1], that matrix inverse can be computed in matrix multiplication time.

4 Dynamic Matrix Inverse In this section we recall an algorithm given by Sankowski [15], supporting column updates in O(n2 ) algebraic operations. As usual, we assume that the matrix remains

Algorithmica (2010) 56: 180–197

185

non-singular during updates. Let us denote by (A)i the i-th column of A. The matrix A after the update of the i-th column to the vector v is given by A = A + (v − (A)i )eiT , where ei is the i-th versor. Thus, we can use Algorithm 1 for updating the inverse. Algorithm 1 Matrix inverse update procedure UPDATE-INVERSE(v, i): b := A−1 (v − (A)i ) B −1 := (I + beiT ) A−1 := B −1 A−1 A := A + (v − (A)i )eiT Let us verify that the procedure UPDATE-INVERSE works correctly. We have: B −1 A−1 = (I + beiT )−1 A−1 = (A(I + beiT ))−1 = (A + AbeiT )−1 = (A + AA−1 (v − (A)i )eiT )−1 = (A + (v − (A)i )eiT )−1 = A−1 . Indeed, the updated inverse is computed correctly. The computation of the vector b in the first step of the procedure takes O(n2 ) arithmetic operations. In the second step we compute the inverse of matrix B, which is given by B −1 = (I + beiT )−1 . The matrix B −1 has only one non-zero column outside the diagonal, and its computation can be done in O(n) time (see [15] for details). Next, we have to perform the multiplication A−1 = B −1 A−1 . Since B −1 has only O(n) non-zero entries, this can be done using O(n2 ) arithmetic operations. Theorem 2 The dynamic matrix inverse problem with non-singular column updates can be solved with the following costs: • initialization O(nω ) arithmetic operations, • update O(n2 ) arithmetic operations (worst-case), • query O(1) arithmetic operations (worst-case). Remark 3 Note that using the same algorithm we can answer the following additional types of queries: determinant We only have to recompute the determinant after each update using the formula det(A ) = det(B) det(A).

(1)

adjoint Here, we query the columns (or entries) of the adjoint matrix adj(A) of a A. The adjoint matrix is adj(A)i,j = (−1)i+j det Aj,i , where Aj,i is A with j -th row and i-th column removed. If the matrix is non-singular, then we can use the following identity adj(A)ij = det(A)(A−1 )ij .

(2)

186

Algorithmica (2010) 56: 180–197

linear system of equations In this problem, we query the solution vector x of a linear system of equations Ax = b, where A and b are subject to changes. After each update we can recompute the solution vector using the formula x = A−1 b.

(3)

For a more detailed treatment of these problems, see [15]. All the algorithms presented in the remainder of this paper also support these additional types of queries.

5 Matrix Inverse with Lookahead In this section we show a dynamic algorithm with update and query time O(nω(1,1,)− ), given a lookahead of n operations. In particular, for a lookahead of n operations, this gives updates and queries in O(nω−1 ) time. In the next section we show how to modify this algorithm to achieve faster queries. 5.1 Update-Only Algorithm We begin with presenting an algorithm computing the vectors b introduced in Algorithm 1 for a sequence of N = n updates in time O(nω(1,1,) ). We assume that N ≤ n as larger values of N do not give speed improvement. Let (i0 , v0 ), . . . , (iN −1 , vN −1 ) be a sequence of N updates and consider the k-th update (ik , vk ). The corresponding vector bk can be expressed as: −1 −1 bk = Bk−1 Bk−2 · · · B0−1 A−1 (vk − (A)ik ).

This follows from formula on b in Algorithm 1 by plugging in the value of A just before the k-th update. Let us define −1 −1 Bp = Bq−1 Bq−2 · · · Bp−1 q

and bk = B0 A−1 (vk − (A)ik ), (q)

q

so that (k)

bk = bk

and bk = A−1 (vk − (A)ik ). (0)

(0)

Now, the vectors bk can be computed for all k = 0, . . . , N − 1 during the initial(q) ization phase of the algorithm and then the bk for q > 0 can be computed once we q (q) know the corresponding B0 . Note that for each k, we only compute bk for selected values of q. (q) The key to computing the bk vectors faster is to: (q)

q

• compute bk -s using Bp -s for large q − p, rather than using Br -s, (q) • compute many bk -s using a single matrix multiplication. Algorithm 2 shows how this can be done. Before analyzing the complexity of the UPDATE-ONLY procedure let us prove the following lemma.

Algorithmica (2010) 56: 180–197

187

Algorithm 2 Update only algorithm UPDATE-ONLY-INTERVAL(p, q): if q − p = 1 then return let m :=  p+q 2  UPDATE-ONLY-INTERVAL(p, m)     (p) (p) (m) (m) bm | . . . |bq−1 := Bpm bm | . . . |bq−1 T −1 m+1 := (I + b Bm m ei m ) UPDATE-ONLY-INTERVAL(m, q) q q Bp := Bm Bpm (m)

UPDATE-ONLY:     (0) (0) b0 | . . . |bN −1 := A−1 (v0 − (A)i0 )| . . . |(vN −1 − (A)iN−1 ) UPDATE-ONLY-INTERVAL(0, N ) A−1 := B0N A−1 q

Lemma 2 For all 0 ≤ p < q ≤ N the matrix Bp − I has at most q − p non-zero columns. m+1 = Proof We prove this fact by the induction on q − p. For q − p = 1 we have Bm (m) T −1 (I + bm eim ) and the induction hypothesis holds. When q − p > 1, we have q

q

Bp = Bm Bpm , for some m, such that p < m < q and q

q

q

Bp − I = (Bm − I )(Bpm − I ) + (Bm − I ) + (Bpm − I ). q

Notice that the product (Bm − I )(Bpm − I ) has the same non-zero columns as q (Bpm − I ). Thus, the number of non-zero columns in (Bp − I ) is not greater than q the number of non-zero columns in both (Bm − I ) and (Bpm − I ). By the induction hypothesis this is at most (q − m) + (m − p) = q − p.  Note that by padding the sequence of updates with trivial updates, we can guarantee that N = 2r for some r ∈ N , so that q − p is always a power of 2. We now analyze the complexity of the internal UPDATE-ONLY-INTERVAL procedure. (m) (m) Let us first consider the matrices used to compute (bm | . . . |bq−1 ). The matrix (p)

(p)

(bm | . . . |bq−1 ) is an n × (q − m) matrix and Bpm is an n × n matrix, but it only has m − p non-zero columns outside the diagonal, so in fact we need a n × (m − p) by q (m − p) × (q − m) multiplication. Similarly, to compute Bp we need a n × (q − m) by (q − m) × (m − p) multiplication.

188

Algorithmica (2010) 56: 180–197

Let T (x) denote the time required to perform UPDATE-ONLY-INTERVAL for a subsequence of updates of length x. We have T (2k ) = 2T (2k−1 ) + 2(n/2k−1 )(2k−1 )ω = 2T (2k−1 ) + 2n2(k−1)(ω−1) . The first equality follows from the fact that multiplication of a n × 2k−1 matrix by a 2k−1 × 2k−1 matrix can be done in time (n/2k−1 )(2k−1 )ω by splitting the first matrix into n/2k−1 square matrices and performing n/2k−1 fast matrix multiplications. By solving the above recurrence, we get T (N) = O(nN ω−1 ) = O(nω(1,,) ). Let us now analyze the UPDATE-ONLY procedure. The UPDATE-ONLYINTERVAL call takes T (N ) time and both matrix multiplications take O(nω(1,1,) ) time. These matrix multiplications are thus the most time consuming part of the algorithm and its time complexity is O(nω(1,1,) ). The cost per update is O(nω(1,1,)− ). 5.2 Answering Queries We will now show how the UPDATE-ONLY algorithm can be modified to support queries. Notice that the i-th column of A−1 after k-th update is equal to B0k A−1 ei . This is exactly the value of bk we would get if vk − (A)ik = ei . Thus, the UPDATEONLY algorithm can be used to answer queries in the following way: • if the k-th operation is a query to the i-th column, then it is substituted with an update (ik , vk ), where ik is arbitrary and vk is such that vk − (A)ik = ei , • these new updates do not modify A so they should be ignored in the computations q of Bp matrices. We get the following: Theorem 4 Given a lookahead of n operations, The dynamic matrix inverse problem with non-singular column updates can be solved with the following costs: • initialization—O(nω ) arithmetic operations, • update—O(nω(1,1,)− ) arithmetic operations (amortized), • query—O(nω(1,1,)− ) arithmetic operations (amortized). The above theorem suggests that queries and updates are equally difficult, but this is not true, as the following theorem shows. Theorem 5 Given a lookahead of a block of operations containing n updates, the dynamic matrix inverse problem with non-singular column updates can be solved with the following costs: • initialization—O(nω ) arithmetic operations, ˆ ˜ + nβ+γ +δ−1 ) arithmetic operations (amortized), ˜ ω(1,1,)−1 ) = O(n • query—O(n ω(1,1,)− ˆ ˜ β+γ +(δ−1) ) arithmetic operations (amortized), ˜ ) = O(n • update—O(n

Algorithmica (2010) 56: 180–197

189

where β, γ , δ are the constants introduced in Lemma 1. Proof To solve the dynamic matrix inverse problem with operation costs as stated in the theorem, we use the FAST-QUERIES algorithm, which is based on the UPDATEONLY algorithm. This algorithm works on a sequence of updates (i0 , v0 ), . . . , (iN −1 , vN −1 ) interspersed with queries, where N = n . Let Qi be the set of queries between the i-th and the i + 1-th update, and let C(Q) p denote the set of columns queried in a set Q of queries. Let ak be the k-th column p of A−1 after p-th update. Recall that computing the ak -s, required to answer the p queries, is equivalent to computing the bk -s for special dummy updates. Also note p that we only need to update a single copy of each ak , even if it is queried multiple times. The FAST-QUERIES algorithm is based on this simple observation, but careful analysis is required to prove the time bounds given in the theorem statement. Algorithm 3 An algorithm with fast amortized queries FAST-QUERIES-INTERVAL(p, q): if q − p = 1 then return let m :=  p+q 2  FAST-QUERIES-INTERVAL(p, m)     (p) (p) (m) (m) bm | . . . |bq−1 := Bpm bm | . . . |bq−1 m+1 := (I + b(m) eT )−1 Bm m im let C = {c0 , . . . , cr } = C(Qm ∪ · · · ∪ Qq−1 )     (p) (m) (m) (m) ac0 | . . . |acr := Bpm ac0 | . . . |acr

answer queries in Qm FAST-QUERIES-INTERVAL(m, q) q q Bp := Bm Bpm FAST-QUERIES:     (0) (0) b0 | . . . |bN −1 := A−1 (v0 − (A)i0 )| . . . |(vN −1 − (A)iN−1 ) for c ∈ C(Q0 ∪ · · · ∪ QN −1 ) ac := (A−1 )c FAST-QUERIES-INTERVAL(0, N ) A−1 := B0N × A−1 (0)

We have already computed the time needed to perform the update-related parts of this algorithm. Also, initialization of the ac(0) -s and answering the queries takes time linear in the number of queries. Thus, we only need to compute the time used to (p) compute the values ac . By a block of operations we will mean a sequence of updates

190

Algorithmica (2010) 56: 180–197

and queries processed in a single call to FAST-QUERIES-INTERVAL procedure, i.e. updates numbered p, . . . , q − 1 and queries from C(Qp ∪ · · · ∪ Qq−1 ). (p)

Lemma 3 The total cost of computing the ac -s for a block of na updates and nb queries can be distributed among these operations so that queries have cost O(n + nβ+γ +δ−1 ) and updates have cost O(nβ+γ +(δ−1) ). (p)

Proof of Lemma 3 Computing the ac -s requires a multiplication of a n × na matrix by a na × nmin(b,1) matrix (in fact, we use the lower-half updates and upper-half queries, but this only changes the time complexity by a constant factor). Let b = min(b, 1). We consider three cases: (i) a ≤ αb : (see Sect. 3 for the definition of α)  In this case our matrix multiplication takes time O(n1+b ) (see Lemma 1), which is O(nb+1 ). (ii) αb < a ≤ b :  In this case our matrix multiplication takes time O(nβ+γ b +δa ). This would give  +δa β+(γ −1)b β+γ +δ−1 O(n ) cost per query. If this is O(n ), then we are done. So let us assume that β + (γ − 1)b + δa > β + γ + δ − 1, which is (1 − γ )b − δa < 1 − γ − δ, and, using the fact that

1−γ δ

= α, we get αb − a < α − .

(4)

If we distribute the cost of the matrix multiplication among the updates we get  O(nβ+γ b −(1−δ)a ) cost per update. Now: β + γ b − (1 − δ)a = β + γ b − α(1 − δ)b + α(1 − δ)b − (1 − δ)a = β + (γ − α(1 − δ))b + (1 − δ)(αb − a) ≤ β + γ − α(1 − δ) + (1 − δ)(α − ) = β + γ + (δ − 1), which is the claimed update cost. (iii) b < a:  In this case our matrix multiplication takes time O(nβ+γ a+δb ) = O(nβ+(γ +δ)a ). β+(γ +δ−1)a β+γ +(δ−1) ). This gives O(n ) cost per update, which is O(n  Since every operation is contained in a logarithmic number of blocks, the overall amortized costs of both operations are at most a logarithmic factor higher than the ones given by the above lemma. 

Algorithmica (2010) 56: 180–197

191

Now, let us see what happens if we want linear time queries. Theorem 6 Given a lookahead of a block of operations containing n updates, the dynamic matrix inverse problem with non-singular column updates can be solved with the following costs: • initialization—O(nω ) arithmetic operations, ˜ • query—O(n) arithmetic operations (amortized),  )−  ˆ ˜ ω(1,1, ˜ β+γ +(δ−1)  ) arithmetic operations (amortized), • update—O(n ) = O(n where   = min(, α). Proof The proof is similar to that of Theorem 5. It is enough to prove the following lemma, similar to Lemma 3. (p)

Lemma 4 The total cost of computing the ac -s for a block of na updates and nb queries can be distributed among these operations so that queries have cost O(n)  and updates have cost O(nβ+γ +(δ−1) ). (p)

Proof of Lemma 4 Computing the ac -s requires a multiplication of a n × na matrix  by a na × nb matrix. As in the proof of Lemma 3 we consider three cases: (i) a ≤ αb :  In this case our matrix multiplication takes time O(n1+b ) (see Lemma 1), which is O(nb+1 ). (ii) αb < a ≤ b :  In this case our matrix multiplication takes time O(nβ+γ b +δa ). This gives  −(1−δ)a β+γ b O(n ) cost per update. Now: β + γ b − (1 − δ)a = β + γ b − α(1 − δ)b + α(1 − δ)b − (1 − δ)a = β + (γ − α(1 − δ))b + (1 − δ)(αb − a) ≤ β + γ − α(1 − δ) ≤ β + γ + (δ − 1)  , which is the claimed update cost. (iii) b < a:  In this case our matrix multiplication takes time O(nβ+γ a+δb ) = O(nβ+(γ +δ)a ). This gives O(nβ+(γ +δ−1)a ) cost per update, which is O(nβ+γ +(δ−1) ) =  O(nβ+γ +(δ−1) ).   6 Matrix Inverse with Lookahead for Element Updates In [4], Demetrescu and Italiano introduced an idea of keeping a matrix in a certain lazy form. This idea was later used by Sankowski [15] in his algorithm for dynamic

192

Algorithmica (2010) 56: 180–197

matrix inverse. It can also be used together with the algorithms presented in the previous section to obtain fast algorithms for the case of element updates. Theorem 7 Given a lookahead of n operations, the dynamic matrix inverse problem with non-singular element updates and column queries can be solved with the following costs: • initialization—O(nω ) arithmetic operations, • query—O(nω(1,1,μ)−μ ) arithmetic operations (amortized), • update—O(nω(1,1,μ)−μ ) arithmetic operations (amortized), where μ satisfies the equality ω(1, 1, μ) − μ = ω(1, μ, ) − . Proof We keep the inverse matrix A−1 in the following lazy form A−1 = N −1 M −1 . The lazy form is recomputed after nμ operations using the following formula M −1 := N −1 M −1 , N := I. We assume that μ ≥ , i.e. the recomputations take place after at least one lookahead block has been processed. The matrix N is maintained using Algorithm 2. Queries and updates to N have to be translated in the following way, in order to take the matrix M into account: update(i, j, x) is translated to update(j, M −1 ei (x − Ai,j )). query(j ) is translated to a query for the linear combination of columns given by M −1 ej . Notice that querying a linear combination of columns of N −1 is equivalent to solving a system of linear equations N x = M −1 ej , so we can still use the Algorithm 2 (see Remark 3). The answers to the queries are obviously correct. Let us verify that the translated updates are correct as well. The matrix A after the update is given by A = A + (x − Ai,j )ei ejT . We have N

−1

M −1 = (N + (M −1 ei (x − Ai,j ))ejT )−1 M −1 = (M(N + (M −1 ei (x − Ai,j ))ejT ))−1 = (MN + MM −1 ei (x − Ai,j )ejT )−1 = (MN + ei (x − Ai,j )ejT )−1 = (A + ei (x − Ai,j )ejT )−1 = A

−1

.

We need the following observation. Lemma 5 The matrix N −1 − I has at most nμ non-zero columns.

Algorithmica (2010) 56: 180–197

193

Proof Note that the matrix N − I after nμ updates has at most nμ non-zero columns. We have N −1 − I = N −1 (I − N ) = −N −1 (N − I ). The product N −1 (N − I ) has at most as many non-zero columns as N − I .



This special form of the matrix N −1 − I allows us to improve the bound on the cost of both multiplications in the UPDATE-ONLY procedure. These matrix multiplications now require only O(nω(1,μ,) ) arithmetic operations. The cost of processing nμ operations has the following three components: • The cost of recomputing the lazy form of A−1 . Computing the product N −1 M −1 requires O(nω(1,1,μ) ) arithmetic operations. • The cost of UPDATE-ONLY-INTERVAL procedure is T (n ) = O(nω(1,,) ) times nμ− . • The cost of matrix multiplications in UPDATE-ONLY procedure is O(nω(1,μ,) ) times nμ− . The amortized cost of one operation is   nω(1,1,μ) nω(1,,)+μ− nω(1,μ,)+μ− O + + = O(nω(1,1,μ)−μ + nω(1,μ,)− ). nμ nμ nμ  This theorem is based on Algorithm 2. We can also use the algorithm with faster queries and this leads to a similar trade-off between query and update time.

7 Offline Planar Matrix Inverse In this section we consider the case of planar matrices, i.e., matrices with non-zero structure corresponding to planar graphs. Gaussian elimination (and thus an LU factorization) of nonsingular planar matrices can be done in O(nω/2 ) using the nested dissection algorithm [11]. This algorithm is based on the separator decomposition of the graph underlying the matrix M. The idea is to reorder the row and columns before the elimination is performed, so as to minimize the running time of the elimination. In particular the rows and columns corresponding to the top-level separator of the decomposition are the last ones. Nested dissection only works if we can guarantee that all diagonal matrix entries are non-zero during the elimination. Otherwise row or column pivoting is necessary and special properties of the separator based row and column orderings are lost. In the remainder of this paper, when we say that we can “guarantee elimination with no pivoting”, it means that diagonal elements are non-zero during the elimination. √ Let us consider a sequence of n update operations on a planar matrix M. We can include all the rows and columns used in these operations in the top-level separator. Let us denote by S the set of all vertices of the graph corresponding to this enlarged separator, and by V the set of all vertices. We will also identify vertices with corresponding rows and columns. Note that the enlarged separator S still has

194

Algorithmica (2010) 56: 180–197

√ size O( n). We now perform the nested dissection algorithm on M, but stop before eliminating the columns and rows corresponding to S. This way, we get the following factorization of M M = LXU. Here, L is lower triangular, U is√upper triangular and X only has non-zero entries on the diagonal and in the last O( n) rows and columns. For any matrix Z, let ZX,Y denote Z restricted to rows from X and columns from Y . We have



IV −S,V −S 0 MV −S,V −S MV −S,S =L U. MS,V −S MS,S 0 XS,S The elements of MS,S are not used in the elimination, so they give additive contribution into XS,S . In other words we have



I MV −S,V −S MV −S,S 0 = L V −S,V −S U, (5) MS,V −S MS,S + 0 XS,S + where is an |S| × |S| matrix representing any additive change of MS,S . We now take a closer look at the inverse of M. M −1 =



(U −1 )V −S,V −S 0

(U −1 )V −S,S (U −1 )S,S



I 0

0

(X−1 )S,S



(L−1 )V −S,V −S (L−1 )S,V −S

0

(L−1 )S,S

.

We get the following equation (M −1 )S,S = (U −1 )S,S (X −1 )S,S (L−1 )S,S .

(6)

While the updates of the part of M corresponding to S can be easily represented√by the updates of X using (5), queries are a bit tricky. Consider a sequence of s = O( n) queries given by pairs (ik , jk ). The k-th query should return the value in ik -th row and jk -th column in (M −1 )S,S . We use the algorithm of Theorem 4 to maintain the inverse of the matrix XS,S . We have to translate the answers to the queries and the queries themselves so as to take the (U −1 )S,S and (L−1 )S,S matrices in (6) into account: • The translated queries are no longer queries of a single matrix entry. They query a linear combination of columns. Coefficients of such a combination can be represented by a vector vk . These vectors are obtained from the updates (ik , jk ) using the formula   (v0 | . . . |vs ) = (L−1 )S,S ej0 | . . . |ejs . Notice that querying a linear combination vk of columns of the inverse (M −1 )S,S is equivalent to solving a system of linear equations (M −1 )S,S x = vk , so we can still use the algorithm of Theorem 4 (see Remark 3). • The answers bk to the queries vk are translated back using the formula T  (a0 , . . . , as ) = ei0 | . . . |eis (U −1 )S,S (b0 | . . . |bs ) , where (a0 , . . . , as ) is a vector of final answers to the queries (ik , jk ).

Algorithmica (2010) 56: 180–197

195

√ All matrices in the above multiplications have size O( n) and so they require O(nω/2 ) arithmetic operations. Using Theorem 4 with  = 1 to update the matrix XS,S we obtain the algorithm for planar matrices with operation costs given by the following theorem. 1

Theorem 8 Given a lookahead of n 2 operations the planar case of the dynamic matrix inverse problem with element updates can be solved with the following costs: • initialization—O(nω ) arithmetic operations, • update—O(n • query—O(n

ω−1 2

ω−1 2

) arithmetic operations (amortized),

) arithmetic operations (amortized),

provided that we can guarantee Gaussian elimination with no pivoting when it is 1 performed (every n 2 operations).

8 Dynamic Transitive Closure in Planar Graphs In Theorem 8 we assumed that we can guarantee Gaussian elimination with no pivoting. Here we show that the reduction of the dynamic transitive closure problem to the dynamic matrix inverse problem works even if we require this assumption. We first need to recall some key concepts behind Theorem 1. The basic idea is the following correspondence between the transitive closure and the adjoint matrix. For any matrix A, let A(x) denote A with a unique variable in place of each non-zero entry. Theorem 9 Let A be the adjacency matrix of a graph G. Then, there exists a path in G from i to j , iff adj(I − A(x))ij is not identically zero (as a polynomial). This is true for finite fields as well. The polynomial adj(I − A(x))ij may have exponentially many terms, so we cannot use symbolic computations to check if it is non-zero. The following lemma, due Zippel [20] and Schwartz [17], gives a general technique of resolving such problems using randomization. Lemma 6 If p(x1 , . . . , xm ) is a non-zero polynomial of degree d with coefficients in a field and S is a subset of the field, then the probability that p evaluates to 0 on a random element (s1 , s2 , . . . , sm ) ∈ S m is at most d/|S|. In order to compute the transitive closure we substitute random numbers from {1, . . . , p} for the variables in A(x), where p is a prime number of order (n1+c ), ˜ over a finite for c > 0. Let A˜ be the resulting matrix. Next, we compute adj(I − A) 1 ˜ field Zp . With probability higher than 1 − nc the element adj(I − A)ij gives a correct answer to a question if j is reachable from i. To prove Theorem 1 we have to maintain ˜ dynamically. adj(I − A)

196

Algorithmica (2010) 56: 180–197

The matrix I − A˜ is non-singular with high probability by Lemma 6. However, as we mentioned at the beginning of this section, we need to guarantee Gaussian elimination with no pivoting. The following lemma shows that it is enough randomize the identity matrix I . Lemma 7 For any column and row order in the matrix A, Gaussian elimination of the matrix I˜ − A˜ fails with small probability. Proof Consider the i-th step of the Gaussian elimination. The eliminated element is 0 with probability p1 , because the random value of I˜i,i gives additive contribution to this element’s value. All random values in I˜ − A˜ are independent, so the probability 1 ).  that the elimination fails on any diagonal element is 1 − (1 − p1 )n = O( nc−1 The non-zero structure of adj(I (x) − A(x)) is obviously identical with the nonzero structure of adj((I − A)(x)). We have shown that Theorem 1 works even if the dynamic matrix inverse algorithms require the guarantee of Gaussian elimination with no pivoting. In particular, Theorem 8 gives an algorithm for the dynamic transitive closure problem with the same time bounds. 9 Concluding Remarks and Open Problems We have presented new algorithms for the problem of matrix inverse with lookahead. Our results show that a subquadratic update and query time in the case of column updates is possible when a lookahead of (n ) size is available, for any  > 0. Our technique can be used also in case of the planar matrices. The algorithms for the matrix inverse can be used to devise new algorithms for the transitive closure problem yielding the most efficient algorithms in the case when any lookahead is available. For a lookahead of n operations the algorithm is essentially optimal. Any improvement in this case would lead to a faster algorithm for boolean matrix multiplication. However, there are still some open problems. In the planar case our algorithm is 1 restricted to a lookahead of n 2 operations. It is not clear if using similar methods can give an algorithm for an arbitrary lookahead. The most intriguing question, however, is whether it is possible to drop the nonsingularity assumption without loss of efficiency. In this case, we would obviously have to work on the adjoint matrix instead of the inverse matrix. This would have several applications, e.g. to a dynamic algorithm for perfect matching testing. Acknowledgements The authors would like to thank their favorite supervisor Krzysztof Diks for his unwavering support. The research was supported by the Polish Ministry of Science, grant KBN-1P03A01830.

References 1. Bunch, J., Hopcroft, J.: Triangular factorization and inversion by fast matrix multiplication. Math. Comput. 28, 231–236 (1974)

Algorithmica (2010) 56: 180–197

197

2. Coppersmith, D.: Rectangular matrix multiplication revisited. J. Complex. 13(1), 42–49 (1997) 3. Coppersmith, D., Winograd, S.: Matrix multiplication via arithmetic progressions. J. Symb. Comput. 9(3), 251–280 (1990) 4. Demetrescu, C., Italiano, G.F.: Fully dynamic transitive closure: Breaking through the O(n2 ) barrier. In: Proceedings of the 40th IEEE Symposium on Foundations of Computer Science, pp. 381–389 (2000) 5. Even, S., Shiloach, Y.: An on-line edge-deletion problem. J. ACM 28(1), 1–4 (1981) 6. Henzinger, M.R., King, V.: Fully dynamic biconnectivity and transitive closure. In: Proceedings of the 35th IEEE Symposium on Foundations of Computer Science, pp. 664–672 (1995) 7. Huang, X., Pan, V.Y.: Fast rectangular matrix multiplication and applications. J. Complex. 14(2), 257–299 (1998) 8. Khanna, S., Motwani, R., Wilson, R.H.: On certificates and lookahead in dynamic graph problems. Algorithmica 21(4), 377–394 (1998) 9. King, V.: Fully dynamic algorithms for maintaining all-pairs shortest paths and transitive closure in digraphs. In: Proceedings of the 40th IEEE Symposium on Foundations of Computer Science, pp. 81–91 (1999) 10. King, V., Sagert, G.: A fully dynamic algorithm for maintaining the transitive closure. In: Proceedings of the Thirty-First Annual ACM Symposium on Theory of Computing, pp. 492–498. ACM, New York (1999) 11. Lipton, R.J., Rose, D.J., Tarjan, R.E.: Generalized nested dissection. SIAM J. Numer. Anal. 16, 346– 358 (1979) 12. Roditty, L.: A faster and simpler fully dynamic transitive closure. In: Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 404–412. Society for Industrial and Applied Mathematics, Philadelphia (2003) 13. Roditty, L., Zwick, U.: Improved dynamic reachability algorithms for directed graphs. In: Proceedings of the 43rd IEEE Symposium on Foundations of Computer Science, p. 679 (2002) 14. Roditty, L., Zwick, U.: A fully dynamic reachability algorithm for directed graphs with an almost linear update time. In: Proceeding of the 36th Annual ACM Symposium on Theory of Computing. ACM, New York (2004) 15. Sankowski, P.: Dynamic transitive closure via dynamic matrix inverse. In: 45th Annual IEEE Symposium on Foundations of Computer Science (2004) 16. Sankowski, P.: Faster dynamic matchings and vertex connectivity. In: SODA ’07: Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 118–126. Society for Industrial and Applied Mathematics, Philadelphia (2007) 17. Schwartz, J.T.: Fast probabilistic algorithms for verification of polynomial identities. J. Algorithms 10, 701–717 (1980) 18. Subramanian, S.: A fully dynamic data structure for reachability in planar digraphs. In: Proceedings of the First Annual European Symposium on Algorithms, pp. 372–383. Springer, New York (1993) 19. Yannakakis, M.: Graph-theoretic methods in database theory. In: Proceedings of the Ninth ACM SIGACT-SIGMOD-SIGART Symposium on Principles of Database Systems, pp. 230–242. ACM, New York (1990) 20. Zippel, R.E.: Probabilistic algorithms for sparse polynomials. In: Proceedings of EUROSAM 79. Lecture Notes in Computer Science, vol. 72, pp. 216–226. Springer, New York (1979)