Articulo 11.dvi CUBO A Mathematical Journal Vol.12, No¯ 02, (169–187). June 2010 A new solution algorithm for skip-free processes to the left Claus Bauer Dolby Laboratories, San Francisco, 94103, USA email: cb@dolby.com ABSTRACT This paper proposes a new solution algorithm for steady state models describing skip-free processes to the left where each level has one phase. The computational complexity of the algorithm is independent of the number of levels of the system. If the skip parameter of the skip-free process is significantly smaller than the number of levels of the system, our algorithm numerically outperforms existing algorithms for skip-free processes. The proposed algorithm is based on a novel method for applying generalized Fibonacci series to the solution of steady state models. RESUMEN Este art́ıculo propone un nuevo algoritmo solución para modelos estado-steady describi- endo procesos libres-salto para la izquierda donde todo nivel tiene una fase. La compleji- dad computacional del algoritmo es independiente del número de niveles del sistema. Si el parámetro de salto de los procesos libre-salto es significativamente pequeño respecto del número de niveles del sistema, nuestro algoritmo numérico supera algoritmos existentes para procesos libre-salto. El algoritmo propuesto se basa en un método reciente para aplicar series de Fibonacci generalizados para la solución de modelos-steady. Key words and phrases: Skip-free processes, Markovian environment, stationary distribution AMS 2000 Subj. Class.: 60J10, 60J99 170 Claus Bauer CUBO 12, 2 (2010) 1 Introduction Skip-free processes have been widely researched as a a way to study packet based communication systems. In particular, skip-free processes to the left have been applied to model the dynamics of finite buffers of switches in data networks that experience the arrival of several data packets at a time while they are only being able to forward one packet per time slot. Using a discrete time model and common queueing theoretic terminology, a skip-free process to the left defines a state model where the (queue) occupancy can move down by one level in a time slot, but might move up by several levels in a time slot. The dynamics of these finite queues are commonly modeled as M/G/1/K queues where K denotes the size of the queue. Similarly, buffers that experience only one packet arrival at a time, but can forward more than one packet per time slot can be modeled as skip-free processes to the right. Algorithms to solve steady state models for skip-free processes were first investigated in [14], [15] by Neuts. In these papers, skip-free processes are considered as a subset of a more general class of steady state models and the application of matrix-geometric methods to solve this class of steady state models is investigated. A steady state analysis that takes explicitly into account the special structure of skip-free processes is proposed in [13, chap. 13]. Combining methods from [5] and [6], skip-free processes are modeled as a special case of Quasi-Birth-and-Death Processes (QBDs). Following the terminology introduced in [13], we define a QBD process as a skip-free process where the system can not move more than one level both downwards or upwards. Several other variations of skip-free processes have been researched in [1], [3], [4], [20], [22], [23]. In this paper, we investigate homogeneous finite skip-free processes to the left, i.e., we assume a finite number of levels and we assume that all levels have the same number of phases. In particular, we assume that each level has one phase. We first give a result of primarily theoretical interest by presenting a new way to derive a closed-form solution solution of a steady state model for skip-free processes to the left. In a second step, we show that the structure of this closed-form solution can be characterized by specific linear recurrent equations. We then apply the theory of general Fibonacci sequences [19] to these linear recurrent equations,. This allows us to derive the main contribution of this paper which is a new solution algorithm for skip-free processes to the left. This solution algorithm has a complexity that is independent of the number of levels of the system. Previous solution algorithms for skip-free models have complexities that also depend on the number of levels of the system. Prior to our work, solution algorithms for steady state models that have a complexity independent of the number of levels of the system were only known for specific classes of QBDs [13, chap. 10.4], [7], [18]. These classes of QBDs do not contain the QBDs used to model skip-free processes in [13, chap. 13]. Finally, we perform numerical experiments to compare the numerical complexity of our algorithm with the numerical complexity of previously known algorithms. Our experiments show that the complexity of our method is lower than the complexity of previous algorithm if the skip parameter is small compared to the overall number of levels of the system. The skip parameter is defined as the maximum number of levels that the queue occupancy can increase in a time slot. We note - without presenting any details in this paper - that our methods can also be applied to derive corresponding results for skip-free processes to the right. In the next section, we provide the exact problem definition. In sec. 3 - 5, we prove and analyze our Theorem 3.1 which gives a closed-form solution for skip-free processes to the left. In sec. 6, we show that the obtained solution can be described by a set of linear recurrent equations. Based CUBO 12, 2 (2010) A new solution algorithm for skip-free processes to the left 171 on this observation, we apply the theory of generalized Fibonacci sequences to obtain a new solution algorithm for skip-free processes to the left. In sec. 7 and 8, we compare the computational complexity of the solution algorithm with previous techniques. We summarize our results in sec. 9. 2 Problem formulation We consider a discrete time Markov chain Xt, t ∈ N on the two-dimensional state space {(n, i) : 0 ≤ n ≤ N, 1 ≤ i ≤ p}, which we partition as ⋃ n≥1 l(n), where l(n) = {(n, 1), (n, 2), ..., (n, p)} for 0 ≤ n ≤ N. The first coordinate n is called the level, and the second coordinate j is called the phase of the state (n, j). In this paper, we only consider the case p = 1. We suppose that the discrete Matrix chain is described by a transition matrix T = (Ti,j )0≤i,j≤N , where Ti,j is the transition probability from the level i to level j. We assume T to be of the following form: j 0 → N T = 0 i ↓ N                         E B2 B3 .. Bm B0 B1 B2 .. Bm−1 Bm B0 B1 .. .. Bm−1 Bm B0 .. .. .. .. Bm .. .. .. .. .. .. B0 Bm−1 Cm B0 Bm−2 Cm−1 .. .. .. .. .. .. .. .. B0 B1 C2 B0 C1                         (2.1) Here the expressions Bj , 0 ≤ j ≤ m, E = B0 + B1, Ci = m ∑ j=i Bj , 1 ≤ i ≤ m are real numbers in the open interval (0, 1). We see that the structure of the matrix T depends essentially on the skip parameter m as the parameter defines which levels can be reached from a current level in a one-step transition. Any level i can be reached from all stages i − j, −1 ≤ j ≤ m − 1 if m − 1 ≤ i ≤ N − 1. By the definition of a transition matrix, we have N ∑ j=1 Ti,j = 1, ∀i, 0 ≤ i ≤ N, 172 Claus Bauer CUBO 12, 2 (2010) which implies that m ∑ l=0 Bl = 1. (2.2) Defining the steady state vector as π = (π0, .., πN ), we can write the steady state equation π = πT as π1B0 = π0(1 − E), (2.3) πkB0 = πk−1(1 − B1) − k ∑ l=2 πk−lBl, 2 ≤ k ≤ m − 1, (2.4) πkB0 = πk−1(1 − B1) − m ∑ l=2 πk−lBl, m ≤ k ≤ N − 1, (2.5) πN (C1 − 1) = − m−1 ∑ l=1 πN−lCl+1. (2.6) We now define new variables C = B0, (2.7) Am−l = 1 − l ∑ n=0 Bn, 1 ≤ l ≤ m − 1. (2.8) Using eqn. (2.2), we see that A1 = Bm = Cm. Using the definitions (2.7) and (2.8), we can rewrite the eqn. (2.3) - (2.5) as follows: π1C = π ∗ 0 Am−1, π ∗ 0 = π0(1 − E)A −1 m−1 (2.9) πkC = πk−1(Am−1 + C) + k ∑ l=2 πk−l(Am−l − Am−(l−1)), 2 ≤ k ≤ m − 1, (2.10) πkC = πk−1(Am−1 + C) + m−1 ∑ l=2 πk−l(Am−l − Am−(l−1)) − πk−mA1, (2.11) m ≤ k ≤ N − 1, πN (C1 − 1) = − m−1 ∑ l=1 πN−lCl+1. (2.12) 3 Closed form solution of the steady state model In this section, we give a closed form solution of the steady state model defined via the transition matrix T in (2.1) or equivalently in (2.9) - (2.12). We show the following Theorem: Theorem 3.1. π∗0 = ( 1 + N−1 ∑ k=1 ∑ bl k E(bm−1, .., b1) + Y (C0 − 1) −1 )−1 , (3.1) πk = π ∗ 0 ∑ bl k E(bm−1, ..., b1), 1 ≤ k ≤ N − 1, (3.2) πN = π ∗ 0 Y (C1 − 1) −1, (3.3) CUBO 12, 2 (2010) A new solution algorithm for skip-free processes to the left 173 where E(bm−1, ..., b1) = C −d ( d cm−1 cm−2 .... c1 ) m−1 ∏ l=1 Acll , (3.4) cl = bl − 2bl−1 + bl−2, 3 ≤ l ≤ m − 1, (3.5) c2 = b2 − 2b1, (3.6) c1 = b1, (3.7) d = bm−1 − bm−2. (3.8) The summation ∑ bl k runs over all bl ≥ 0, 1 ≤ l ≤ m − 1 such that cl ≥ 0 for 1 ≤ l ≤ m − 1, In this summation, the variable bm−1 only takes the fixed value bm−1 = k. Further, Y = − m−1 ∑ l=1 Cl+1 ∑ bl N−l E(bm−1, ..., b1). (3.9) Remarks: We state some remarks for later use in this paper: 1. By definition, m−1 ∑ l=1 cl = d. (3.10) 2. For any fixed bl+1 over which is summed in the summation ∑ bl k , the number of bl over which is summed in ∑ bl k is limited by bl ≤ l l + 1 bl+1. (3.11) This follows from the summation condition c2 ≥ 0 and the eqn. (3.6) for l = 1. For l ≥ 2, we apply the induction principle. Assuming that eqn. (3.11) holds for l, then by the summation condition cl+2 ≥ 0 2bl+1 ≤ bl + bl+2 ≤ l l + 1 bl+1 + bl+2, which implies eqn. (3.11) for l + 1. 3. The eqn. (3.11) implies that d ≥ 1. (3.12) 4 Proof of Theorem 3.1 We first recall a well-known identity for multinomial coefficients. For any set of strictly positive integers a1, .., ak with n = k ∑ j=1 aj , there is ( n a1 a2 ... ak ) = k ∑ j=1 ( n − 1 a1 .. (aj − 1) ... ak ) . (4.1) 174 Claus Bauer CUBO 12, 2 (2010) For any integer x, 1 ≤ x ≤ m − 1 and a given set of variables b1, ..., bm−1 as defined in (3.4), we introduce the additional variables bxl and c x l defined as follows: bxl = bl − max(0, x − (m − 1 − l)), l ≥ 1.. (4.2) cxl = b x l − 2b x l−1 + b x l−2, 3 ≤ l ≤ m − 1, (4.3) cx2 = b x 2 − 2b x 1 , (4.4) cx1 = b x 1 , (4.5) dx = bxm−1 − b x m−2. (4.6) We now state two lemmas which we will need for the proof of Theorem 3.1. Lemma 4.1. For 1 ≤ x ≤ m − 1, cxl = cl − { 1 if l = m − x, 0 else . } (4.7) dx = d − 1. (4.8) Proof of Lemma 4.1: We prove the eqn. (4.7) by considering different ranges of the variables l, m and x : Range 1: For l ≥ m − x + 1, l ≥ 3, cxl = cl − x + m − 1 − l + 2x − 2m + 2 + 2l − 2 − x + m − 1 − l + 2 = cl. Range 2: For l = m − x, l ≥ 2, cxl = cl − x + m − 1 − l + 2x − 2m + 2 + 2l − 2 = cl + x − m + l − 1 = cl − 1. Range 3: For l = m − x = 1 cxl = cl − 1. Range 4: For l ≤ m − x − 1, l ≥ 1, cxl = cl. We now prove the eqn. (4.8): dx = bm−1 − x − (bm−2 − max(0, (x − 1)) = d − 1. For later usage, we note that the eqn. (3.10), (4.7), and (4.8) imply that for 1 ≤ x ≤ m − 1, there is m−1 ∑ l=1 cxl = d x. (4.9) Lemma 4.2. E(bm−1, .., b1) = m−1 ∑ x=1, cm−x 6=0 Ex(bxm−1, ..., b x 1 )Am−xC −1. CUBO 12, 2 (2010) A new solution algorithm for skip-free processes to the left 175 Proof of Lemma 4.2: We first note that in the definition (3.4), if for any x there is cm−x = 0, then C−d ( d cm−1 cm−2 cm−x+1 cm−x cm−x−1.... c1 ) m−1 ∏ l=1 Acll = C−d ( d cm−1 cm−2 cm−x+1 cm−x−1.... c1 ) m−1 ∏ l=1, l 6=m−x A cl l , i.e., we can neglect the contribution of the quantity cm−x. Now, Lemma 4.2 follows from the definition (3.4) and the eqn. (3.12), (4.1), (4.7), and (4.8). Proof of Theorem 3.1: First, we prove eqn. (3.2) by induction over k. For k = 1, we see bm−1 = cm−1 = d = 1, whereas bl = cl = 0 for l ≤ m − 2. Thus, we obtain from eqn. (3.2) that π1 = π ∗ 0 Am−1C −1 as required by eqn. (2.9). For the inductive step, we now assume that eqn. (3.2) holds for πi, 1 ≤ i ≤ k − 1, and prove it for i = k < N. In view of the eqn. (2.10) and (2.11), we see that in order to show that the eqn. (3.2) holds for πk, it is sufficient to prove the following two equations: πk = min(m−1,k) ∑ x=1 πk−xAm−xC −1, (4.10) πk−1 = min(m−1,k−1) ∑ x=1 πk−1−xAm−xC −1. (4.11) For the proof of the eqn. (4.10), we first consider the case m − 1 < k. Replacing the probabilities πl, k − m + 1 ≤ l ≤ k, in eqn. (4.10) with the right hand side of eqn. (3.2), we obtain ∑ bl k E(bm−1, ..., b1) = m−1 ∑ x=1 Am−xC −1 ∑ b̃l k−x E(b̃m−1, ..., b̃1). (4.12) Applying Lemma 4.2 to the left hand side of eqn. (4.12), we see that in order to prove eqn. (4.10) we have to show the following identity: m−1 ∑ x=1, cm−x 6=0 Am−xC −1 ∑ bl k E(bxm−1, ..., b x 1 ) = m−1 ∑ x=1 Am−xC −1 ∑ b̃l k−x E(b̃m−1, ..., b̃1). (4.13) For the proof of eqn. (4.13), it is sufficient to prove the following two claims: 1. For each integer x, 1 ≤ x ≤ m − 1, and for each set bxm−1, ..., b x 1 derived via the relation (4.2) from a set bm−1, ..., b1 over which is summed in the sum ∑ bl k , and for which cm−x 6= 0, there exists a set b̃m−1, ..., b̃1 over which is summed in the sum ∑ b̃l k−x for which holds cxl = c̃l for all l, 1 ≤ l ≤ m − 1, and dx = d̃. 2. For each integer x, 1 ≤ x ≤ m − 1, and for each set b̃m−1, ..., b̃1 over which is summed in the sum ∑ b̃l k−x , there exists a set bxm−1, ..., b x 1 derived via the relation (4.2) from a set bm−1, ..., b1 176 Claus Bauer CUBO 12, 2 (2010) over which is summed in the sum ∑ bl k , for which cm−x 6= 0, and c̃l = c x l for all l, 1 ≤ l ≤ m − 1, and d̃ = dx. In order to show claim 2, we set eb̃l = b̃l + max(0, x − (m − 1 − l)). Further, we define e c̃ l and e d̃ l via e b̃ l in the same way cl and d are defined via bl as in (4.3) - (4.6). By definition, we see that (e b̃ l ) x - defined as in (4.2) - equals b̃l which implies that c̃l = (e c̃ l ) x and d̃l = (e d̃ l ) x. Now, in order to prove claim 2, it remains to show that the set eb̃m−1, ..., e b̃ 1 belongs to the summation ∑ bl k , i.e., eb̃m−1 = k, e c̃ m−l ≥ 0 if l 6= x and ec̃m−x > 0. The first claim is obvious by the definition of e b̃ m−1 = b̃m−1 + x = k − x + x = k. We note that reversing the proof of the relation (4.7), we can show that ec̃l = c̃l + { 1 l = m − x 0 else. } (4.14) As c̃l ≥ 0, the eqn. (4.14) implies that e c̃ l ≥ 0, 1 ≤ l ≤ m − 1, l 6= m − x, and e c̃ m−x > 0, q.e.d. In order to show claim 1, we have to show that the set bxm−1, ..., b x 1 belongs to the summation ∑ bl k−x . For this purpose, we have to show that bxm−1 = k − x and c x l ≥ 0, 1 ≤ x ≤ m − 1. The first inequality follows from eqn. (4.2). The second relation follows from eqn. (4.7) and the fact that on the left-hand side of eqn. (4.13) we only sum over cm−x 6= 0. In the case m − 1 ≥ k, we argue as above and have to show that k ∑ x=1, cm−x 6=0 Am−xC −1 ∑ bl k E(bxm−1, ..., b x 1 ) = k ∑ x=1 Am−xC −1 ∑ b̃l k−x E(b̃m−1, ..., b̃1). (4.15) Eqn. (4.15) is shown in the same way as eqn. (4.13). Eqn. (4.11) is shown in the same way as eqn. (4.10). Eqn. (3.3) follows from (2.12) and (3.2). Eqn. (3.1) follows from (3.2), (3.3), and the eqn. N ∑ i=0 πi = 1. 5 Complexity of the calculation of πk using Theorem 3.1 The inequality (3.11) gives an upper bound for the number of bl over which is summed in the sum- mation ∑ bl k . Using bm−1 = k, eqn. (3.11), and the well-known relation, N ∑ k=1 km = 1 m + 1 N m+1 + O(N m), we obtain ⌊ (m−2)k m−1 ⌋ ∑ bm−2=1 ⌊ (m−3)bm−2 m−2 ⌋ ∑ bm−3=1 ..... ⌊ b2 2 ⌋ ∑ b1=1 1 = 1 (m − 1)! m−1 ∏ l=2 ( l − 1 l )l km−1 + O(km−2). CUBO 12, 2 (2010) A new solution algorithm for skip-free processes to the left 177 Thus, the complexity of the calculation of the steady state probability of π∗0 using the relation (3.1) is ∼ 1 (m−1)! m−1 ∏ l=2 ( l−1 l )l N−1 ∑ k=1 km−1 ∼ 1 m! m−1 ∏ l=2 ( l−1 l )l N m. Here, we assume that the values of the multino- mial coefficients have been pre-computed as they do not depend on the actual values of the transition matrix T. Once π∗0 is calculated, the calculation of πk for any 2 ≤ k ≤ N − 1 requires O(1) steps if we assume that the values of the sums ∑ bl k already determined for the calculation of π∗0 have been stored. The calculation of πN requires O(m) steps. In summary, we see that the overall computational complexity of Theorem 3.1 is ∼ 1 m! m−1 ∏ l=2 ( l−1 l )l N m. For m > 3 and large values of N, this complexity is too high for any prac- tical applications. In the next section, we will show how the complexity of Theorem 3.1 can be significantly reduced. 6 A reduction of the complexity of Theorem 3.1 In this section, we apply the theory of linear recurrence equations to Theorem 3.1. This will allows to find a computationally very efficient way of computing the sum ∑ bl k E(bm−1, .., b1). First, we recall some facts from the theory of linear recurrent equations: Lemma 6.1. We consider a sequence of real numbers xn, n ≥ 1 defined via a set of real numbers a1, .., an−1 and a set of strictly positive real numbers b1, .., bn−1 as follows: xi = ai, 1 ≤ i ≤ n − 1, (6.1) xi = n−1 ∑ k=1 bkxi−k, i ≥ n. (6.2) Then, we have: xn = w ∑ u=1 fu−1 ∑ v=0 cu,vn vrnu . (6.3) The numbers ru, fu, and w are defined via the characteristic polynomial yn−1 + n−1 ∑ k=1 bky n−1−k. (6.4) ru are the roots, fu are the multiplicities of the roots ru, and w is the number of different roots of the characteristic polynomial (6.4). Obviously, w ∑ u=1 fu = n − 1. The numbers cu,v are the solutions of the (n − 1)− equations ai = w ∑ u=1 fu−1 ∑ v=0 cu,vi vriu, 1 ≤ i ≤ n − 1. (6.5) Proof: The proof can be found in [10]. We now apply Lemma 6.1 to the eqn. (4.12). Setting xk = ∑ bl k E(bm−1, .., b1) in (6.1) and (6.2), we obtain from Lemma 6.1: 178 Claus Bauer CUBO 12, 2 (2010) Lemma 6.2. ∑ bl k E(bm−1, .., b1) = w ∑ u=1 gu−1 ∑ a=0 tu,ak ahku, Here, the numbers gu denote the respective multiplicities of the w different roots hu of the characteristic polynomial xm−1 − m−1 ∑ u=1 Am−uC −1xm−1−u. (6.6) The numbers tu,h are the solutions of the (m − 1)- linear equations ∑ bl k E(bm−1, .., b1) = m−1 ∑ u gu−1 ∑ h=0 tu,hv hhvu, 1 ≤ v ≤ m − 1. We now use Lemma 6.2 to simplify Theorem 3.1 as follows: Theorem 6.1. π∗0 = ( 1 + N−1 ∑ k=1 Fk + Y ∗(1 − C0) )−1 , (6.7) πk = π ∗ 0 Fk, 0 ≤ k ≤ N − 1, (6.8) πN = π ∗ 0 Y ∗(C1 − 1) −1, (6.9) where Fk = w ∑ u=1 gu−1 ∑ a=0 tu,ak ahku. (6.10) The numbers w, hu, gu, tu,h are defined as in Lemma 6.2, and Y ∗ = − m−1 ∑ l=1 Cl+1FN−l. Proof: The eqn. (6.7) - (6.10) follow directly from Theorem 3.1 and Lemma 6.2. 7 The complexity of Theorem 6.1 7.1 Preliminaries In this section, we investigate the computational complexity of Theorem 6.1. We see from the formu- lation of Theorem 6.1 that the complexity of the calculation of a steady state probability πk mainly derives from the calculation of the sum N−1 ∑ k=1 Fk, the solution of the (m − 1) eqn. (6.5), and the cal- culation of the roots of the characteristic polynomial (6.6). It is well-known that the solution of the system of equations (6.5) using Gaussian elimination requires O(m3) operations. We consider the complexity of the other two operations separately in the next two sections. CUBO 12, 2 (2010) A new solution algorithm for skip-free processes to the left 179 7.2 The calculation of the sum N−1 ∑ k=1 Fk We see from eqn. (6.10) that the sum N−1 ∑ k=1 Fk can be rewritten as follows: N−1 ∑ k=1 Fk = w ∑ u=1 gu−1 ∑ a=0 tu,a N−1 ∑ k=1 kahku := w ∑ u=1 gu−1 ∑ a=0 tu,aHN−1,a,hu , (7.1) where HN,q,x := N ∑ k=1 kqxk. For the calculation of HN,q,x we prove the following lemma: Lemma 7.1. HN,q,1 = q+1 ∑ k=1 (−1)q−k+1BEq−k+1 p + 1 ( p + 1 k ) N k, (7.2) HN,0,x = x xN − 1 x − 1 , (7.3) HN,q,x = x (N + 1)qxN + (N + 1)q − 2 x − 1 − x x − 1 q−1 ∑ b=0 HN,b,x ( q b ) , x 6= 1, q > 1. (7.4) In eqn. (7.2), BEk denotes the k − th Bernoulli number [2]. Proof: The relation (7.2) is explained in [2] and the relation (7.3) is a known formula for geometric sums. For the proof of eqn. (7.4), we will use the technique of partial summation defined in [17] as follows: For any complex series ak, bk, M < n ≤ N, there is ∑ M 4 no general exact solutions exists. Thus, the solution of polynomials of higher degree requires the application of approx- imative methods. An overview of these methods is given in [8]. Among these methods, Laguerre’s method has been widely researched [16]. Whereas its convergent behavior is well understood, the speed of its convergence is not well described in the literature. Therefore, we performed numerical experiments in order to understand the computational com- plexity of the calculation of the roots of the characteristic polynomial (6.6). For each degree m, we evaluated 10000 polynomials of degree m, and for each polynomial we used a random generator to generate the coefficients of the polynomial. For each polynomial, we measured the execution time s(m) - measured in milliseconds - needed to determine all roots of the polynomial on a Pentium IV, 2.4 Ghz PC. The graph in fig. 1 shows the logarithm - to the base 10 - of s(m) as a function of the logarithm - to the base 10 - of the degree m of the polynomial. An least mean square analysis of the simulation data shows that log(s(m)) can be described as a linear function of log m as follows: log(s(m)) = 1.99 log m + c for some constant c. Thus, the execution time s(m) can be approximately described as a second order polynomial of the polynomial degree m. 2 2.5 3 3.5 4 4.5 5 5.5 6 0 1 2 3 4 5 6 7 Logarithm of degree m of polynomial Lo ga rit hm o f e xe cu tio n tim e in m s Figure 1: Speed of convergence of the calculation of the roots of the characteristic polynomial 7.4 Combined complexity of Theorem 6.1 The combinations of the results of sec. (7.1) - (7.3) shows that the overall computational complexity required to calculate a specific steady state probability πk using Theorem 6.1 is O(m 3). 8 Comparison of Theorem 6.1 with previous work In this section, we compare the computational complexity of Theorem 6.1 with known methods to solve steady state models as defined via the transition matrix (2.1). First, we present a well-known recursive way of solving the system described by the matrix (2.1). 182 Claus Bauer CUBO 12, 2 (2010) Then, we follow an idea from [6] and model the Markov process defined via (2.1) as a QBD process with u := ⌈ N +1 m−1 ⌉ levels. We discuss two known methods to solve QBD models due to Gaver et. al [5], [13, chap. 13] and Ye and Li [21], respectively. Finally, we compare the computational complexity of these approaches with the complexity of Theorem 6.1. 8.1 Recursive calculation The equations (2.3) - (2.6) allow a recursive calculation of the steady state probabilities πk, 0 ≤ k ≤ N. We see from these equations that each πk can be expressed as a linear function of all πl, 0 ≤ l ≤ k − 1. Thus, using a recursive argument each πl can be expressed as the product of π0 and a real number kl. As N ∑ l=0 kl = π −1 0 , we can can calculate all steady state probabilities πk. These calculations require O(N 2) calculations. We note that a similar approach can be used to derive a closed - form solution for πk. In particular, we consider the unconstrained random walk corresponding to (2.1) without a reflecting barrier on the right, i.e., with an infinite number of levels. In order to define this random walk formally, we generalize the transition probabilities Bn defined in (2.1) to Dn = { Bn, if 0 ≤ n ≤ m, 0, n ∈ Z\[0, m]. } . (8.6) We define the transition probabilities from state i to state j of the unconstrained random walk as pij = { D0 + D1, if i = j = 0, Dj−i+1, else. } . We denote by {mj : j ≥ 0} the invariant measure of this Markov chain and set m0 = 1. Then, mj = ∑ i≥0 pij mi. Using the definition (8.6), we see m1 = 1 − D0 − D1 D0 , D0mj+1 = mj − j ∑ i=0 miDj−i+1 (j ≥ 1). (8.7) The first N columns of the transition matrix defined by (8.6) are identical to the first N columns of the transition matrix T defined in (2.1). Thus, there exists a constant K such that πj = Kmj if j < N. There is πN = K m−1 ∑ i=1 Cm−i+1mN−m+i + πN C1, and πN C1 = πN (1 − D0) = πN − D0 + D0K N−1 ∑ j=0 mj , (8.8) which implies K = D0 D0 N−1 ∑ i=0 mi + m−1 ∑ i=1 Cm−i+1mN−m+i . CUBO 12, 2 (2010) A new solution algorithm for skip-free processes to the left 183 Let D(s) := ∑ k≥0 skDk be the generating function of the series Dk. Using standard manipulations with generating functions yields M (s) := ∑ k≥0 mks k = D0 1 − s B(s) − s , (8.9) valid in the disc |s| < q, where q is the least positive solution of D(s) = s. As Dn = 0 for n > m, B(s) is a polynomial and therefore M (s) is meromorphic function with a finite number of poles. Using a partial fraction expansion of M (s) and picking the coefficients of sk on the right hand side of (8.9) gives a representation of the mk. The calculation of the coefficients uses the recursive relations (8.7) and relation (8.8). Thus, it has a complexity of O(N 2). 8.2 Modeling of the Markov process as a QBD process with u := ⌈ N +1 m−1 ⌉ levels First, we describe how the matrix T can be redefined as a block matrix that describes a QBD process. For the sake of simplicity, we first assume that N + 1 = u(m − 1) for some integer u. Then, we can rewrite T as defined in (2.1) as a u × u block matrix as follows: T =                    P1,1 A0 0 0 .. .. .. .. 0 A2 A1 A0 .. .. .. .. .. 0 0 A2 A1 .. .. .. .. .. .. 0 .. .. .. .. .. .. .. .. 0 .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. A1 A0 0 .. .. .. .. .. .. A2 A1 A0 0 0 .. .. .. .. .. A2 Pu,u                    , (8.10) A0 =                    Bm+1 0 0 .. .. .. .. .. .. Bm Bm+1 0 .. .. .. .. .. .. Bm−1 Bm Bm+1 .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. B5 .. .. .. .. .. Bm+1 0 0 B4 B5 .. .. .. .. Bm Bm+1 0 B3 B4 .. .. .. .. .. Bm Bm+1                    , 184 Claus Bauer CUBO 12, 2 (2010) A1 =                    B2 B3 B4 .. .. .. .. .. Bm B1 B2 B3 .. .. .. .. .. Bm−1 0 B1 B2 .. .. .. .. .. Bm−2 0 0 B1 B2 .. .. .. .. Bm−3 .. .. .. B1 B2 .. .. .. .. .. .. .. .. B1 B2 .. .. .. .. .. .. .. .. B1 B2 .. .. .. .. .. .. .. 0 B1 B2 B3 0 0 0 0 0 0 0 B1 B2                    , (8.11) A2 =                    .. .. .. .. .. .. .. 0 B1 .. .. .. .. .. .. .. 0 0 .. .. .. .. .. .. .. .. 0 .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. ..                    , P1,1 =                    E B3 B4 .. .. .. .. .. Bm B1 B2 .. .. .. .. .. .. Bm−1 0 B1 B2 .. .. .. .. .. Bm−2 0 0 B1 .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. B1 B2                    , (8.12) Pu,u =                    B2 B3 B4 .. .. .. .. Bm−1 Cm−1 B1 B2 B3 .. .. .. .. Bm−2 Cm−2 0 B1 B2 .. .. .. .. Bm−3 Cm−3 0 0 B1 B2 .. .. .. Bm−4 Cm−4 .. .. .. B1 B2 .. .. .. .. .. .. .. .. B1 B2 .. .. .. .. .. .. .. .. B1 B2 .. .. .. .. .. .. .. 0 B1 B2 C2 0 0 0 0 0 0 0 B1 C1                    . CUBO 12, 2 (2010) A new solution algorithm for skip-free processes to the left 185 Table 1: Comparison of complexity of different solution algorithms Method Complexity Recursive calculation O(N 2) Linear level reduction O(N m2) Method of folding O(m3log N m ) Theorem 6.1 O(m3) The matrices A0, A1, A2, P1,1 and Pu,u are (m − 1) × (m − 1) matrices. In the case that N + 1 is not an integer multiple of m − 1, i.e, N + 1 = u(m − 1) − c, 0 < c < m − 1, the presentation (8.10) changes as follow. The matrix Pu,u is replaced by the matrix (m − 1 − c) × m matrix P ∗ n,n which consists of the m− 1 − c right columns of Pn,n. Similarly, the matrix A0 in right-most column of (8.10) is replaced by the (m − c) × m matrix A∗0 which consists of the m − c right columns of A2. Gaver et. al [5], [13, chap. 13] propose the linear level reduction method to solve QBD systems. Applied to (8.10), the complexity of this solution algorithm is O(um3) = O(N m2). Ye and Li [21] propose the method of folding. The application of this algorithm to solve (8.10) requires O(m3log2 u) = O(m3log2 N m ) computational steps. 8.3 Comparison of Complexities We see from Table 1 that the complexity of Theorem 6.1 is lower than the complexity of the recursive calculations if m = o(N 2/3). Its complexity is lower than the complexity of the linear reduction based method if m = o(N ). Theorem 6.1 also outperforms the method of folding if N is sufficiently larger than m. However, the performance gain is only of the order (log N m )−1, whereas Theorem 6.1 achieves a performance gain of the order m/N compared to the linear level reduction method. We note that in contrast to the recursive calculation, the linear level reduction approach and the folding method, the number of computations required to calculate a specific steady state probability πk using Theorem 6.1 does not depend on the number of phases N + 1. In summary, we see that Theorem 6.1 outperforms the other approaches presented in this section if N is significantly larger than m. Thus, the applicability of Theorem 6.1 to the efficient solution of Markov models for skip-free processes depends on the specific choice of the parameters N and m describing the skip-free process. 9 Conclusions This paper proposes a new solution algorithm for skip-free processes. The complexity of the algorithm is independent of the number of levels of the system. In consequence, it numerically outperforms previous algorithms if the skip parameter is small compared to the number of system levels. It is based on a novel method for deriving a closed-form solution for skip-free processes and analyzing this solution using Fibonacci sequences. 186 Claus Bauer CUBO 12, 2 (2010) Received: July 2008. Revised: April 2009. References [1] A. Adhikari, Skip free processes. Ph.D thesis, Berkeley, 1986. [2] C. B. Boyer, Pascal’s Formula for the Sums of Powers of the Integers. Scripta Math. 9, 237-244, 1943. [3] P.J. Brockwell, The extinction time of a birth, death and catastrophe process and of a related diffusion model. Advances Applied Probability, Vol. 17, 1985, 17 - 42. [4] M.F. Chen, Single birth processes, Chinese Ann. Math. Ser. A, 20:77-82, 1999. [5] D.P. Gaver, P.A. Jacobs and G. Latouche, Finite birth-and-death models in randomly changing environments. Advances in Applied Probability, 16:715-731, 1984. [6] W.K Grassmann and D.A. Stanford, Matrix analytic methods. In: Computational Probabil- ity, WK Grassmann (Ed.), Kluwer, 2000, pp 153-202. [7] L. Guen and A.M Makowski, Matrix-geometric solution for finite capacity queues with phase- type distributions. Performance’87, edited by Courtois, P.J.; Latouche, G.; Amsterdam, 1987, p. 269 - 282. [8] E. Hansen, M. Patrick and J. Rusnack, Some modifications of Laguere’s method. BIT, 17(1977), 409 -417. [9] T. W. Hungerford, Algebra. Springer Publishing House, 1974. [10] W.G. Kelley and A.C. Peterson, Difference equations, An introduction with applications, Academic Press, Inc. 1991. [11] P.A. Jacobs, D.P. Gaver and G. Latouche, Finite markov chain models skip free in one direction. Naval Research Logistics Quarterly, 31, 1984, pp. 571-588. [12] E. Kuntz, Algebra., Vieweg Verlag, Braunschweig, 1991. [13] G. Latouche and V. Ramaswami, Introduction to matrix analytic methods in stochastic modeling. Society for Industrial and Applied Mathematics, 1999. [14] M.F. Neuts, Matrix-geometric solutions in stochastic models. An Algorithmic Approach. The John Hopkins University Press, Baltimore, MD, 1981. [15] M.F. Neuts, Structured stochastic matrices of M/G/1 type and their applications. Marcel Dekker, New York, 1989. [16] Y.V. Pan, Solving a polynomial equation : Some history and recent progress. Siam Rev., Vol. 39, No. 2, pp. 187 -220, June 1997. [17] K. Prachar, Primzahlverteilung., Berlin, Heidelberg, New York, Springer Verlag, 1978. CUBO 12, 2 (2010) A new solution algorithm for skip-free processes to the left 187 [18] W.J. Stewart, On the use of numerical methods for ATM model. Modeling and performance evaluation of ATM technology, edited by Perros, H.; Pujolle, G.; Takahashi, p. 375 - 396. [19] D.A. Wolfram, Solving generalized Fibonacci recurrences. The Fibonacci Quarterly 36.2. May 1998, 129-45. [20] S.J. Yan and M.F. Chen, Multidimensional Q-processes. Chin. Ann. Math. Ser. A,7:90-110, 1986 [21] J. Ye and S.Q. Li, Folding algorithm: A computational method for finite QBD processes with level dependent transitions. IEEE Transactions on Communications., 42:625-639, 1994. [22] J.K. Zhang, Generalized birth-death processes. Acta Mathematica Sinica, Vol. 46, 1984, 241 - 259. [23] Y.H. Zhang, Strong ergodicity for single-birth processes. Journal of Applied Probability, Vol. 38, 2001, 207 - 277.