THE UNIVERSITY OF MICHIGAN SYSTEMS ENGINEERING LABORATORY Department of Electrical Engineering College of Engineering SEL Technical Report No. 35 THE SOLUTION OF QUASI BIRTH AND DEATH PROCESSES ARISING FROM MULTIPLE ACCESS COMPUTER SYSTEMS by Victor L. Wallace March 1969 This research was supported by United States Air Force Contract AF 30(602)-3953.

This report was also a dissertation submitted in partial fulfillment of the requirements for the degree Doctor of Philosophy at The University of Michigan, 1969. ii

Acknowledgments I wish to express my appreciation to Professors Keki B. Irani, Kuei Chuang, Ralph L. Disney, Eugene L. Lawler, and Arch W. Naylor for their service on my doctoral committee, and for their advice and counsel. I particularly thank Professor Irani for his sustained interest and encouragement as chairman of the committee. For her expert job of typing the final manuscript, I also express my appreciation to Miss Sharon Bauerle. This work was supported by the Systems Engineering Laboratory through Air Force contract AF 30(602)-3953, granted by the Rome Air Development Center. This support was sincerely appreciated. Finally, I thank my wife, Mary, for the patience and understanding which made this possible. Victor Lew Wallace Ann Arbor, Michigan March 27, 1969 iii

Table of Contents List of Tables List of Figures Abstract Page vi vii ix Section 1 Introduction 2 The QBD Process 3 Q Matrices and Notational Conventions 4 Existence of a Quadratic Root 5 Boundary Leading QBD Processes 6 Determination of R* 7 Equilibrium Distributions 8 Recurrent Processes 9 Uniqueness of S* 10 Conclusions References 1 20 32 44 65 72 82 92 94 97 v

List of Tables Table Page 1.1 States of Figure 1. 1, when m=2, and m'=1 6 1. 2 The Transition Intensity Matrix of the Model of Figure 1.1 7 vi

List of Figures Figure Page 1.1 Illustrating a Model of a Multiple Access System 3 2.1 Illustrating the States of N = < N1, N2 > 12 5. 1 Graphs of Boundary Leading QBD Processes 45 5. 2 Illustrating Example 5. 1 49 5. 3 Illustrating Example 5. 2 50 5. 4 Illustrating Example 5. 3 51 5. 5 Illustrating Example 5. 4 52 5. 6 Graphs of U3 for Example of Figure 5. 1 54 5. 7 Graphs of Uk for a Particular N 59 7. 1 A Graph of a Process for which p(R) > 1 and -T is Convergent 74 8. 1 A Reducible QBD Process and One of its Subgraphs 87 vii

ABSTRACT Quasi birth and death processes (QBD-processes) are continuous time Markov chains whose transition intensity matrices can be represented in tri-diagonal block form, with the blocks along a diagonal being identical submatrices after the first. They are a general class of infinite state Markov chains which are useful as models of the "high traffic" behavior of large, multiple-access computer systems. They are particularly useful in analysis and design of those systems which are often described as "computer utilities". Classical queueing theory offers no general technique for determination of equilibrium distributions of QBD-processes, particularly not for QBD-processes of the degree of complexity likely to be encountered in computer applications. Furthermore, because the transition intensity matrix of a QBD-process is of infinite order, commonly used numerical techniques for determination of equilibrium probabilities are also not applicable. However the equilibrium equations can be shown to have a formal similarity to the equilibrium equations of simple birth and death processes with matrices and vectors in the former replacing scalars in the latter. R. V. Evans has proposed that this similarity be used to derive a numerical procedure which involves several separate finite operations, each of which does not tax the memory or computing capabilities of a moderately large digital computer. Fundamental to this procedure is the evaluation of matrix roots of a matrix "characteristic equation", and the use of those matrix roots to ix

generate "geometric" solutions to the equilibrium equations. The chief problem is to determine weak sufficient conditions for a QBDprocess to be amenable to this procedure. This dissertation rigorously develops the procedure, primarily through the use of modern analysis, matrix theory, and generalized difference equations. In the process, the notion of the "boundary leading" QBD-process is introduced. The boundary leading QBD-process is a very general QBD-process which includes all those which have a finite number of essential classes. In particular it includes all irreducible QBD-processes. It is shown that if a QBD-process is boundary leading then a suitable matrix root of the characteristic equation exists, as well as a solution of the equilibrium equations having the geometric form. This result is used to show that, for irreducible processes, the procedure can be applied, will determine an equilibrium distribution whenever one exists, and will provide an appropriate indication of failure when no equilibrium distribution exists. Other useful results regarding marginal equilibrium distributions, uniqueness of equilibrium distributions and roots of characteristic equations, and positive recurrence in QBD-processes are also developed. From this theory, one can conclude that the numerical solution of a very broad class of QBD-processes is feasible, and is well within the capability of a moderately large digital computer. Thus, these processes become available as models for the practice of computer system design and analysis and for other disciplines which use Markov chain models. x

1. Introduction Markovian models have found such considerable application [8, 9, 14, 15, 24] to the art of computer system design for two significant reasons. First, they have been applied because there is a sizable theory upon which to draw, and many powerful theoretical results which are available. Second, computer systems have come to a stage of evolution where they are predominantly "high-traffic"' systems and, as a consequence, the modeling and estimation of statistical performance has become a crucial aspect of their design. However, there are limitations to the models of this typelwhich are capable of analysis. The models which must be solved are often too complex for classical queueing theory, and often have too many states for modern numerical techniques [23]. The precision required over a wide range of parameter variations often rules out simulation methods as too costly. The study of multiple-access computer systems and the socalled "computer utilities" is particularly difficult in this regard. Furthermore, as these systems grow larger and more complex (through provision of more user circuits, introduction of multiprocessing, use of satellite computers and data-concentrators, etc.) this situation can be expected to become worse. Complexity of models, and largeness of state spaces will prevail. Yet these systems are the very systems for which analysis of Markov chain models can be 1

2 of greatest use to the system architect. Thus, new techniques are urgently needed. A typical example of such a model for a multiple-access computer system will serve to identify some of the problems. Suppose one were considering a multiple-access computer system with two thousand common-carrier trunk-lines leading to a bank of fast processors operating on a large multiprogrammed, paged main core memory. Suppose also that the memory is demand-paged with a large bulk-core "virtual" memory. Finally, suppose our purpose is to study the relationship between the expected number of queued unserviced requests for computation and the size of main core memory, the degree of multiprogramming (i. e. the number of service requests which have some main core memory committed to them), or the number of processors committed. A queueing model which could be used for this study is illustrated diagrammatically in Figure 1.1. This model anticipates that congestion in the channel to the bulk core will be a major problem [24]. (If it is not, of course, other models, similar in nature, would be examined.) It is assumed that connected users generate requests which "arrive" at a mean rate X per second, at times approximately governed by a Poisson process. Each request is queued in a "wait" queue (first in, first out), until it can be serviced. A maximum of m requests is permitted to be assigned main core memory at one time, and so requests are taken from the wait queue only if the

COMMAND ARRIVALS x L- nn i-KJ-^irnn -~ -> n3= n2^ 4j^~|3 r~42PI no nl+n2+n _ 3 n n p= <m 3 n2 n 3 P2 Figure 1.1 Illustrating a Model of a Multiple Access System

4 number so assigned is less than m (i. e. the "block" in the diagram is open). Once a request is given main core it will quickly demand that a number of pages be swapped ("burst paging") before serious execution can begin, and the bulk-core channel time for this to be accomplished has a mean 1/1 seconds. If all m' processors are busy, the request queues for processing after the burst. (This queue can assume at most the length m-m'.) Execution proceeds for 1/112 seconds (mean) on any one of the processors, until another page transfer is demanded. Again there is a queue for page transfers. The latter require 1/13 seconds (mean) of bulk core transfer time. If both paging and burst-paging demands exist simultaneously, the former is given preemptive priority (which is represented on the diagram by a blocked completion of the "burst paging"). The paging may have been the result of final output of a completed command, an executive-induced interruption of the process which suspends eligilibility for main core, or (most likely) a simple page demand to be followed by more execution. We assume each of these occurs with probability pl, p2, and p3, respectively, and result in exit, return to the wait queue, or return to the execution queue, respectively. This model will be a continuous-parameter Markov chain provided only that the time intervals described are independent, conditionally on state, and have one of the "nice" distributions: exponential, Erlang, hyperexponential, etc. For most purposes, the exponential distribution will be simplest and most expedient.

5 To make things as simple as possible for this illustration, we let the degree of multiprogramming m be equal to 2, the number of processors m' to be equal to 1, and we scale time to the mean page swap-time, so that L3=1. The states will consist of quadruples <n n 1, n2, n3 >, where nO is the number of commands in the wait queue, n1 is the number in the burst phase (0 or 1), n2 is the number in or waiting for execution, and n3 is the number in or waiting for paging. Those states which are possible are listed and consecutively numbered in Table 1. 1. Since nO can assume all nonnegative integer values, there are an infinite number of states, and states 10 to 15 repeat with n0=2, 3,.... The transition intensity matrix for this model is represented in Table 1. 2. It repeats the pattern within the dotted lines indefinitely. (The symbol x represents the diagonal entry, which for convenience is not evaluated. The dotted lines have been added to show the grouping of states for which n0=O, n0=1, n0=2, etc.) Note that this model offers small hope of successful closedform analytical solution. Furthermore, truncation to a finite matrix to approximate the process by one which could be solved numerically with standard methods is not possible, since nO can be expected to have a mean value of several hundred. For truncation to be a good approximation the range of no would need to be much larger. Even for the small m, m' used for illustration, the number of states would quickly become prohibitive for numerical work.

Table 1.1 States of Figure 1.1, when m=2, and m'=l n0 n1 n2 n3 0 0 0 2 1. 10. 2. 0 0 3. 0 0 4. 0 1 5. 0 1 6. 0 1 7. 0 0 8. 0 0 9. 0 0 1 1 2 0 11. 12. no n1 n2 n3 1 0 0 2 1 0 2 0 1 0 2 0 1 1 0 0 1 1 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 13. 14. 15. 16. 1g. 18. 19. 20. 21. 22. n n n2 n3 2 0 0 2 2 0 1 2 0 2 0 2 0 0 21 0 1 2 1 1 0 etc.

' 00000 0'00000 0 00 O0 X o o o o o o o0 0 0 X cC O o o o o o ol o o: 0o p OO0'b X 0 0 0 ~~~~I =I 0o 0o o o oi o o o o o o0 0 O O OO 00 0 0 0 0 X OO0 0000 0 0 0 0 0000 0 Y OD O OOO I'O O O O O O 0 X O O O O 0 C) - C OO WCD, 00000 0 X000 00 3, o o o oo ol x o o o o o o o to o O o o o o00 o 000 o0'o 0'00 0 I o o o o o o 00000 o o ooooooo oooUoo oooo0 z o oR oJ ol o o o o~ olw o o o o o! >o. ffI I o o oO 0 0 O o O 0 0 0 o oCo' 0 - o o o o" o 0o Yo00:o 0o, 000 ------------- -----------— ^ —-------------------— I —------------------ ------------------------- - o x oo o: o o o 00 Y 0 0 0 0'-'o'X o o o o 00 0 I L.'000000' 00 00t2: 0 0 o00 o o 0, o o o o o 0 o 0 o0 0 0 or. 0: 0 0 >' 0 o o' 01 J)I.. 0- Co 0 0'>*' 0''000001 0' I 1'! 0 o 0 o oCo o, o, Y o o i o o o O 0 W 0 00 o o ooo -: II| I 00 0' 0 O' t 0 0 > 0 0 0 >O 0 0 00 0 - - - - -- - - - - - -- - - — *- - -- - - -- - -- - - - -- - - - — r - - - - - - - - - -- - - - L

8 The need for large state spaces demonstrated in this example is a natural consequence of treating multiple-access systems in which the number of users simultaneously active is great. While the earliest multiple-access systems had about ten active users, it will not be surprising to see, in the near future, systems with several thousandz users simultaneously connected and generating requests at a rapid rate. This in turn implies that queues of waiting requests will be correspondingly long while performance is maintained at the same levels. For example, if one assumes an average system response time of one second, and an average request rate of one request every ten seconds for every user, then a system which has two thousand simultaneous users must be able to handle at least two hundred requests every second. It will also have an average of two hundred user requests awaiting service at any given time. Thus, it is often more reasonable, under such circumstances, to use a model which allows the number of waiting requests to assume any nonnegative integer value, as was done in the above example. Fortunately, it turns out that these problems frequently exhibit useful structural characteristics in common. In particular, the pattern of repetition in the matrix of Table 2. 2 should be particularly noted. The dotted blocks repeat diagonally, and the form of the matrix is block-wise tri-diagonal. It should also be observed that although the size of the blocks would increase considerably for larger values of m or m', the pattern described would be preserved. This

9 pattern is, in fact, common to a large quantity of useful models for multiple -access computing systems (and other applications). It is a characteristics feature of the infinite state Markov chains which we call quasi birth and death processes (QBD processes) and which are the primary subject of this dissertation. Disney [3] observed the occurrence of this pattern in some finite Markov queueing problems, and offered a method for their solution which appears to be most effective when the number of repetitions is not too large. The infinite (or QBD) process was first identified and treated by Evans [4], and his work represents the genesis of our approach. The chief purpose of this dissertation is to increase the arsenal of Markovian models available to system designers for the treatment of computer system (and other) design problems. A numerical technique has been proposed [4] for the solution of equilibrium distributions for QBD processes. These processes represent good models for many processes where one state variable assumes a very large number of possible values, as do computer utility problems when the "number of tasks awaiting completion of service" is a state variable. This variable can change by at most unity for any single state transition, and the probability intensity of each possible transition in state does not depend upon the number of tasks awaiting completion of service unless there is less than a certain number of them. Using structure imposed by these properties it will be seen that the equilibrium distributions of such a process can be solved by a procedure which

10 involves several separate, finite operations, each of which does not tax the capabilities of a moderately large digital computer. In this dissertation, the procedure will be rigorously derived, and sufficient conditions will be established for its application to a given QBD process to be valid. To this end, the notion of a "boundary leading" QBD process is introduced and shown to have very useful characteristics. (Almost every interesting QBD process is boundary leading. )

2. The QBD-Process Let mo and m be positive integers, with m > m. (2.1) Let N = {N(t), t > 0} be a continuous parameter Markov chain with stationary transition probability functions continuous at t= 0, and with state space ~l, where = {<0, 1>, <0, 2>,..., <0, mo>, <1,1>, <1, 2>,..., <m, 1>, <2, 1>,..., <2,m>,..,...}. (2. 2) The process N is called a QBD process if its transition intensities {Ujk; j, k E d' }, satisfy the properties: (a) If j = <jl, j2> and k = <k1, k2> are two states for which Ijl-k1l > 2, then ujk = 0. (b) Ifj = <jl, j2> and k = <k1, k2> are two states for which Ij1-k1 < 1, if(jl+kl) > 3, and ifj'= <jl-l, j and k' = <k1-1, k2>, then ujk = Ujk,. A continuous-parameter Markov chain will also be called a QBD process if a mere one-one mapping of its states results in a QBD process. However, no loss in generality will result from consideration of only state spaces of the form (2. 2), represented by the set n1. The state space'9 is more graphically described as a twodimensional integer space whose first dimension is countably infinite, whose second dimension is finite, and which has the configuration of 11

12 Figure 2.1. States in the set = {<0, 1>, <0, 2>,... <0, m>} will be called boundary states and the set / will be called the boundary of N. mo 0 2 ( 0 0 0 0 3 0 0 0 0 2 0 0 0 0 m o o o o.. 0 1 2 3 4 *' n, Figure 2.1 Illustrating the States of N = < N1, N2 > If the lexicographic ordering of the states used in equation (2. 2) is preserved, then the conditions (a) and (b) are equivalent to a specification that the transition intensity matrix U = {ujk; j, k E 1 } is an infinite matrix which can nevertheless be displayed in the partitioned form E A A C B A 0 U = (2.3) 0

13 where the submatrix E is an mo x mo matrix; where A, B, and C are mx m matrices; and where A and C are m x m and mx m 0 0 matrices, respectively, and are further partitioned into the forms A A A = (2.4) and C = [C 0] (2.5) Many QBD processes whose transition intensity matrices do not have the appearance of Equation (2. 3) can nevertheless be identified as QBD processes by simply regrouping the states and reidentifying them. Two examples are suggestive of the power of the model. Example 2.1: The transition intensity matrix E'F' G'B A C B A C B A where E is an mo' x m0' matrix, where A, B, and C are m x m, where FT and GT are mo' x m and m x m' respectively, and where mo' is not necessarily greater than m, can be converted to the form of Equation (2. 3) by reidentifying states so that

14 B G' E = F' EJ Then E is m x m, wherem = m + m > m. Example 2. 2: The transition intensity matrix E' F' G' H' A' D' C' B' A' D' C' B' A' can be converted to the form of Equation (2. 3) by reidentifying its states so that' F' E = B G' H 0 0 A A' O B' A' B C' B and

15 DI C' C = 0 DI provided the degrees of E' and H' together are greater than twice the degree of B'. (If not, the result is the given form in example 2. 1, and can be further altered.) Of course, the transition matrices of QBD processes need not be so easily converted to the standard form. Practically speaking, though, the application will usually suggest the appropriate state identification whereby every possible jump in state is at most a unit step in the first coordinate. When mo = m = 1, the QBD process is readily recognized to be a simple birth and death process. The "QBD" in the name "QBD process" is an acronym for the somewhat barbaric, though apt descriptor "quasi birth and death", in recognition of this fact. (The prefix "quasi" is frequently used in matrix theory in describing matrix forms created by replacing scalar elements by submatrices, as in "quasi diagonal" and "quasi triangular" matrices.) In order to calculate an equilibrium distribution of a QBD process N, reasoning reminiscent of that commonly followed in solving birth and death processes can be followed. This reasoning and the calculation technique that results is essentially the same, in outline, as that proposed by Evans [4, 5, 6]. An equilibrium distribution is a sequence f = {rk' k e O } of

16 nonnegative real numbers which satisfies the equation r 2L = 0 and whose components have unit sum. If we let 7- = [ %,, 1'.2.. ] be a partitioning of that sequence for which 00 is an mo-component row vector, and 1,'2',..., are m-component row-vectors, then the equation ir U = 0 is equivalent to the equations %fE + 01C = 0 (2.6a) oA + 01B + 2C = 0 (2.6b) ~n-1A + nB + n+lC = (2. 6c) n=2, 3,... Now suppose that there is a solution [ 0t' 1',2''. ] to Equations (2. 6) of the form A 1 = %oR (2.7a) n ='1Rn-1 n=2, 3,..., (2.7b) A where R is an m x m matrix, where R is an mo x m matrix partitioned to the form A R R=, (2.8) and where 0o is some row-vector whose components are not all zero. Then it develops that, by substitution of Equations (2. 7) in Equation (2. 6), the matrix R and the vector f must satisfy the equations k%(E + RC) = 0 (2. 9a) (A + B RC) = (2. 9b) P(A + RB + R RC) = 0 (2.9b)

17 A n-2 2 o RR (A + RB + R C) = 0. (2.9c) Furthermore, since a solution of Equations (2. 6) where %0, O1, 2',..., have only nonnegative components is sought, (o and R must guarantee that nonnegativity as well. If we can show that such (o and R exist, and yield a solution 7T which is an equilibrium distribution, then the assumption of Equations (2. 7) will be vindicated. The Equations (2. 9) suggest a sufficient condition for sucn a solution to be found. This condition is that a matrix R exists having the properties: (a) All entries of R are nonnegative ^ A (b) E + RC is singular, and has a nonnegative, nonzero vector (0 in its null-space (so that Equation (2. 9a) is satisfied). 2 (c) All entries of A + RB + R C are zero. (This latter property is reminiscent of the classical characteristic roots of a difference equation, while the first two are reminiscent of the application boundary conditions to difference equation solutions.) It will be our purpose to show sufficient conditions that a matrix R does exist having the properties (a), (b) and (c). Remarkably, those conditions are very general, and include almost all interesting QBD processes. Under those conditions, an equilibrium distribution can be calculated by the four steps (1) Determine a matrix R having properties (a), (b) and (c). (2) Determine a vector 0o

18 (3) Normalize o, if possible, so that 7T will have unit sum. (4) Determine as many of the n as desired using Equations (2. 7a) and (2. 7b). It will be shown also that such a matrix R can be determined by an iterative process in every interesting case, and that.o can also be calculated by known techniques. Thus, the problem of determining an equilibrium distribution of a QBD process (which has an infinite state space) will be shown to reduce to several routine finite problems, each of which does not tax the memory or computing capabilities of a digital computer of moderate-to-large size. This dissertation represents a significant extension of the very original work of Evans. We extend his work in several ways. In the first place, he considered only QBD processes for which A was a scalar multiple of the identity matrix (A = XI). Secondly, he required that the process N be "irreducible with minimum path of order 2m". Thirdly, he considered only positive recurrent processes. In what follows, the matrix A is arbitrary (consistent, however, with U being a transition intensity matrix), the QBD process N is assumed only to be "boundary leading" (an assumption weaker than "irreducible"), and existence of a matrix R having the properties (a), (b), (c) is proved without any reference to "recurrence". (We do, however, still require positive recurrence,or something very close to it,for finding R.)

19 In addition, results are given here providing (1) necessary and sufficient conditions on R for an irreducible QBD process to be positive recurrent, (2) sufficient conditions for R and ir to be unique, with r equal to the limiting distribution. These relaxations of Evans' assumptions considerably increase the depth of the results. The presentation, relying as it does almost exclusively on algebra and modern analysis rather than probability theory, also offers a sharper, more rigorous definition for computational purposes.

3. Q-matrices and Notational Conventions The proofs of this work are primarily algebraic, and are concerned in large measure with a type of matrix known as a Q-matrix. In this chapter, a convenient notation for some frequently used matrix relations will be established, the term "Q-matrix" will be defined, and several of its more important properties described. The machinery thus provided is necessary to the work that follows. For any two matrices X = {jk, j e J, k e K} and Y = {Yjk, j e J, ke K} of equal dimensions (perhas infinite) the relations X = Y, X > Y, X > Y will be defined to signify X=Y: xjk = Yjk for every j J and k K, X > Y: xjk > Yjk for every j e J and k K, X > Y: xjk > Yjk for every j J and ke K. The special relation X > Y will signify that the relation X > Y holds, but that the relation X = Y does not. In other words, when X > Y, every element xjk is greater than or equal to Yjk, but at least one element exhibits inequality. The symbol Q will be used for a matrix whose elements are all zero. Similarly, the symbols 1 and 0 will be used for row or column vectors whose components are all 1 or 0, respectively. The square matrix I, of course, is the identity matrix. Where it will not be ambiguous, the context will determine dimensions, and whether row- or column-vectors are intended. 20

21 A matrix or vector X will be called nonnegative (non-positive) if every element is nonnegative (non-positive). In other words X is called nonnegative if X > Q. A matrix will be called positive (negative) if X > 0 (X < 0). The row-sums of a matrix X are the components of the vector XI, and are found by adding all the elements in each row of the matrix. A square matrix X = {Xjk; j, k e K} is said to be stochastic if it is nonnegative and its row-sums are all unity (i.e. XI = 1). A square matrix Q = {qjk; j k e K} is called a Q-matrix when its non-diagonal elements are not negative, and its row sums are not positive. In other words, the matrix Q is a Q-matrix when qjk > 0 whenever j i k, and Q1 < 0. The transition intensity matrix of any continuous parameter Markov chain is always a Q-matrix, as are the matrices U, B, and E of the QBD process N. By the same token, the matrices A and C are nonnegative. A matrix X = {Xjk; j, k e K} is said to be reducible if its index set K can be divided into two disjoint nonempty subsets K1 and K2 such that xjk = for all j e K1 and ke K2. The matrix X is said to be irreducible otherwise. We say, further, that index j leads to index k (under X) whenever the indices K cannot be divided into two disjoing nonempty subsets K1, K2 containing j and k, respectively, and such that xjk, = 0 for all j' K1 and k' e K2. Further, we say that j leads directly to k (under X) if qjk 0. Two indices j, k are said to communicate, on the other hand, if j leads to k and k leads to j. By

22 fiat, an index of a matrix will be said to communicate with itself. Therefore, the relation "i communicates with j under X't is symmetric, transitive, and reflexive, and partitions the indices into disjoint subsets called communicating index classes (or, simply, index classes) such that two indices belong in the same class if and only if they communicate. An index which communicates with every index it leads to is called essential; otherwise inessential. An index class which contains an essential index must contain only essential indices, and therefore the class is called an essential index class. (Some authors [17] prefer the term "minimal closed".) An irreducible matrix must have but a single index class, and that class must be essential. An index class K of a Q-matrix Q is called a conservative index class if qjk = KK for every j e K. The matrix Q is called a conservative Q-matrix if every index class of Q is conservative, or in other words if Q1 = 0. For every finite conservative Q-matrix Q, the equation xQ = 0. is true for at least some row-vector x > 0. It should be noted that a continuous parameter Markov chain is irreducible (reducible) if its transition intensity matrix is irreducible (reducible). Furthermore, the indices of the transition

23 intensity matrix are t'states" of the chain, and thus all terms for indices of a Q-matrix apply also to states of a chain. Note, however, that this terminology has been developed so that it can be used on Qmatrices regardless of whether or not they can be interpreted as transition intensity matrices. The remainder of this chapter is concerned only with finite matrices. The singularity or non-singularity of a finite Q-matrix is, of course, an important property. Since every row of a conservative Q-matrix has zero sum, its columns must be linearly dependent vectors, and Q clearly must be singular. A much more useful result is the following theorem: Theorem 3.1 A Q-matrix is nonsingular if and only if it has no essential index class which is conservative. The proof of this theorem will make use of two lemmas, the first of which is very minor. Lemma 3.1 Let M be a reducible matrix, partitioned in the form M 0 M = (3.1) M21 M22 where Ml and M22 are square suomatrices. Then M is nonsingular iff Ml1 and M22 are both nonsingular.

24 Proof: This lemma is only a minor alteration of the known result [11, page 431] that the determinant of a quasi triangular matrix (as Equation (3. 1)) is equal to the product of the determinants of its diagonal submatrices. Obviously M has nonzero determinant only if both M11 and M22 do also, and conversely if M11 and M22 have nonzero determinant then M does also. Thus, the lemma is true. The second lemma is due to Collatz [ 2 ] and Taussky [20], among others. It will not be proved here. Lemma 3.2 An irreducible nonconservative Q-matrix is nonsingular, and its inverse is (strictly) negative. With that we can proceed to the proof of Theorem 3.1. Proof (Theorem 3. 1): Suppose that a Q-matrix Q has an essential index class K1 which is conservative. If K1 contains all indices of Q, then Q is a conservative matrix, hence singular. If, on the other hand, K1 does not contain all indices of Q, then we can permute its rows and columns so that it can be partitioned to the form Q11 0 Q = (3.2) Q21 Q22 where the row sums of Q11 are all 0. Qll is clearly a conservative Q-matrix, and singular. By Lemma 3.1, then, Q is singular. This proves (modus ponens) that Q is nonsingular only if it has no

25 essential index class which is conservative. To prove the converse, let Q be n x n and have no essential index class which is conservative. We proceed by induction on the number of index classes. Thus, if it has only one index class then Q is irreducible and the index class must be essential. Therefore, the matrix must be nonconservative and,by Lemma 3.2, nonsingular. Now, suppose Q has i+1 (1 < i < n-1) index classes, and it is true that a Q-matrix is nonsingular if it has i index classes, none of which is both essential and conservative. At least one class must be essential, and therefore nonconservative. Then the rows and columns of Q can again be permuted so that it can be partitioned to the form of Equation (3.2) where Qll is irreducible and nonconservative (its indices corresponding to an essential index class of Q), and where Q21> Q. Every index class of Q22 must correspond to an index class of Q, and therefore Q22 must have i index classes. It can have no essential index class which is conservative, for if it did the corresponding rows of Q21 would necessarily be zero, and the corresponding indices of Q would be both essential and conservative, contradicting the main hypothesis that Q has no such classes. Consequently, Q22 is nonsingular. Since Q11 and Q22 are nonsingular, by Lemma 3. 1 Q is nonsingular. Therefore, by induction, if the Q-matrix Q has any number of index classes and has no essential index class which is conservative, it is nonsingular. This completes the proof of Theorem 3.1.

26 We will have a need for more detailed properties of the inverse of a Q-matrix which do not appear to be common currency, mostly because they do not assume irreducibility of the matrix. But first, we quote a well-known lemma which is the basis for the "method of partitioning" (or "escalator method") for the numerical inversion of matrices. We refer to [13, p. 78] for proof. Lemma 3.3 Let M be a partitioned nonsingular matrix M11 M12 M = (3.3) M21 M22 Let Ml1 and M22 be nonsingular, and let the matrix W = M be partitioned with the same partitioning as M, W11 W 12 W = (3.4) W21 W22 Then W11 and W22 are nonsingular and w (1-M - (3.5) Wll = (Mll - M12M22M21) (5) -1 W12 = - M1l 12W22 (3.6)'-1 W21 = (M22 M21MlW1 (3. 7) 22 = (M22 - M21M l1 (3 8)

27 The next two lemma follow almost directly. Lemma 3. 4 Let Q be a nonsingular Q-matrix. Then the inverse of Q is non-positive Q- 0 (3.9) and its diagonal elements are strictly negative. Proof: This proof also follows an induction on the number of index classes in Q. Because of its similarity to the proof of the converse part of Theorem 3. 1, this proof will only be sketched. For one index class, the result follows from Lemma 3.2. Assuming that nonsingular Q-matrix with i index classes is non-positive with negative diagonals, we observe that a (i+1) index class nonsingular Q-matrix can always be partitioned (after permutation) as Q11 Q12 Q Q21 Q22 where Q11 is irreducible, and Q22 has i classes. Then we can proceed to show that the matrices Ql' (Q11 - Q12Q22Q21)' and (Q22 - Q21QllQ12) all are nonconservative Q-matrices and have non-positive inverses. This can be done by observing row-sums, under the constraint that h(Q11 + Q12) and (Q22 + Q21) mus have nonpositive row-sums and that at least one of those sums must be nonzero, and by observing that the number of index classes in a

28 Q-matrix (i.e. Q1l and Q22' in this case) cannot be increased by adding a nonnegative matrix. Application of Equations (3. 5)-(3. 8) of Lemma 3.3 then demonstrates the non-positivity. The diagonal property is the result of the fact that (Q11 - Q12Q22Q21) 1 and are strictly negative (because they are irreducible), and that (Q22 - Q21QlQ12) and Q22 have negative diagonals (by the induction hypothesis). This completes the sketch of the induction. Lemma 3.5 Let M be a nonsingular matrix. Then index j leads to index k under M if and only if index j leads to index k under M. Proof: Suppose j does not lead to k under M. Then (by definition of "leads to") M can be partitioned into the form M11 ~ M = M21 M22 where j is in the upper partition of indices and k in the lower. By Lemma 3. 3, the inverse of M can be shown to be correspondingly partitioned to W-1 SWl o WClearl21 W Clearly j cannot lead to k under Q. The converse is proven by

29 identical reasoning, replacing M by W and vice versa. This completes the proof (modus ponens). With these results, we have finally: Theorem 3.2 Let Q be a nonsingular Q-matrix. Then the entries Wjk of its inverse W = Q1 are nonzero if and only if jk j leads to k under Q. Proof: By Lemma 3.5, wjk C 0 implies that j leads to k under Q, so that only the converse needs to be proved. We therefore suppose that wjk = 0 and examine Q. Since WQ = I, we can write jk wjiqi + Wjtqq = 0 id for every 2 [ j. (3.10) Let I. be the set Ij = {i: wi 0}and note that since w.. < 0 (by Lemma 3. 4), it follows that j e I.. Note further that wji < 0 (by Lemma 3. 4) for all i and j, that ql < 0 (Q is nonsingular!) for every 2, and that qi. > 0 for all i ] i. It follows that the summation over i j=j in Equation (3. 10) must be zero for every f 4 I., and that qi = 0 for every 2 4 I. and every i e I.. (Observe that 2 4 I. implies 2 ] j.) Thus no index in j j I. leads to one not in I.. But, by hypothesis wjk = 0 so that k I, and, as observed previously, j e I.. Thus, the theorem is completed.

30 Putting together the equivalences from Lemma 3. 5 and Theorem 3. 2, we get the obvious corollary: Corollary 1 If a Q-matrix Q is nonsingular, then index j leads to index k under Q1 if and only if j leads directly to k under Q. (Equivalently, j leads to k under W = Q if and only if the entry wjk is nonzero.) We conclude this section with a simple theorem about sums of Q-matrices or nonnegative matrices. Theorem 3.3 If X is either a nonnegative matrix or a Q-matrix, and if Y is either a nonnegative matrix or a Q-matrix having the same index set as X, then every essential index class of X + Y contains at least one essential index class of X and at least one essential index class of Y. Proof: Each index class of X and each index class of Y must be contained in an index class of X+Y, since the addition of positive, nondiagonal entries to a matrix cannot zero-out non-diagonal entries already present in the original and since, therefore, states which communicate under X or Y communicate under X+Y. Let K be an essential index class of X+Y. Then no indices of K lead to indices not in K under X (nor, separately, under Y either) for otherwise the same

31 would be true for X+Y. Then K must contain an essential index of X (and, separately, Y), and therefore one of the index classes of X (or Y) contained in K must be essential. This proves the theorem.

4. Existence of a Quadratic Root A central question in this theory is whether or not a matrix R exists having the properties (a), (b), and (c) described in Chapter 2. These properties can be rewritten, using the notation of Chapter 3, as the properties: (a) R > 0, A A (b) (E + RC) is lpsingular, and a row-vector To > 0 exists for which 7o(E + RC) = 0, (c) A + RB + R2C = 0. If a matrix R having those properties does not exist, then our construction scheme cannot be applied to finding the equilibrium distribution XT of a QBD process N, and the analogy between the solution of QBD processes and the solution of simple birth and death processes breaks down. In this chapter and the next we are concerned with some sufficient conditions that such a matrix exists. Put in other terms, the question is whether there exists a root of the quadratic equation (c) which satisfies properties (a) and (b). In the scalar (birth and death process) case, with A=X, B = -(X++L), C=g, and E = -X, the quadratic becomes the familiar one A - R(X+/) + R2g = 0, and has the two roots R = Xk/ and R=1. Only the root X/g satisfies condition (b), so that there is exactly one 1x1 matrix R satisfying all 32

33 three conditions. The question is not So simple when m > 1. Unlike a scalar quadratic, a matrix quadratic equation A + XB + X C = 0 may have more than two roots. In fact, it has potentially as many roots as there are combinations of 2m objects taken m at a time, where m is the degree of the coefficient matrices. (See, for example, the treatment of matrix polynomials in [11, pp. 227-239].) This multiplicity is not of direct concern, since only a restricted subset of the roots (those satisfying condition (b)) is of interest. However, it does offer an indication that a study of the roots of A + XB + X C = Q must proceed by means different from those which apply to scalar quadratics. The conditions (a) and (b) will be strengthened somewhat in this section, for convenience. (Since the existence of roots satisfying stronger conditions implies their existence under the weaker conditions, no loss of significance results. Our proof will be adequately general.) Let, be the set of nonnegative matrices R' R' > 0 (4.1) for which the condition R'C1 = Al (4. 2) is true. Then one can say: Lemma 4.1 If R e ~, then conditions (a) and (b) are true.

34 Proof: Condition (a) is true by definition. To prove condition (b) is true, we first observe that the transition intensity matrix U is an infinite, conservative Q-matrix. Therefore, by Equation 2. 3, the following properties hold for the matrices A, B, C, and E: (A + B + C) = 0, (4. 3) (E +A)l = 0. (4. 4) Equation (4. 4) implies that -E = (4. 5) and Equations (2. 5 and (2. 8) imply that ^ ^ RC 0 RC RC Q 1= -. (4.6) Since R is in 6, then RC1 = Al, and ^ A Al RC1 = - El. (4.7) 0 Thus AA (E + RC)l = 0. (4. 8) But E is a Q-matrix, and R and C are nonnegative. It follows, AA then, that E + RC has nonnegative off-diagonal entries, and zero row A A sum. Thus E + RC is an mo x mo conservative Q-matrix, and a nonnegative row-vector To > 0 exists for which

35 of(E + RC) = 0. This proves the lemma. The existence of a quadratic root R e - is explored through a study of a fixed point of a certain mapping. Let a be the set of nonnegative m x m matrices S which satisfy the restrictions S1 = Al (4. 9) and (B + S) is nonsingular. (4.10) (This set may be empty, but we will show later that in all interesting cases it is not.) The matrix B+S is clearly a Q-matrix for every S e,, since (1) the addition of a nonnegative matrix cannot make the off-diagonal elements of the Q-matrix B nonnegative, and (2) the row sums of B+S are non-positive by equations (4. 3), (4. 9) and (4. 10), and by the construction - Cl = (B+A)1 = (B+S)1 (4.11) The matrix -(B+S) C is also stochastic. First, a multiplication of Equation (4. 11) by (B+S) shows -(B+S) Cl = 1. Second, Lemma 3. 4 shows that (B+S)-1 < Q. Thus -(B+S) C is nonnegative, and has unit row-sums. These results will be expressed, for future reference, as the lemma:

36 Lemma 4. 2 If S is a nonnegative matrix S > 0 for which Sl = Al, then the matrix (B+S) is a Q-matrix. If, in addition, (B+S) is nonsingular (and therefore S is in j ), then -(B+S) 1C is a stochastic matrix. Now, let T be a mapping on 4 defined by T(S) = - A(B+S) C (4.12) Then we have a third lemma which will be found useful in what follows: Lemma 4. 3 For each S e, T(S) is nonnegative, and T(S)_ = Al. Proof: Since B+S is a nonsingular Q-matrix for all Se c, then by Lemma 3. 4, -(B+S)1 is nonnegative. Therefore T(S) = -A(B+S)- C is a product of nonnegative matrices, and is nonnegative. Multiplying the left and right of Equation (4. 11) by A(B+S) gives - A(B+S) C1 = Al, showing that T(S)1 = Al, and proving the lemma. The connection between fixed points of T and quadratic roots in O/ can now be shown. A fixed point S* of T is a matrix S* e for which S* T(S*) BS-1C. (4.13) S* = T(S*) = -A(B+S*) C. (4. 13)

37 Theorem 4. 1 If T has a fixed point S* e o, then R* = -A(B+S*)-1 2 is a root of the quadratic A + XB + X C = and is in (6. Proof: Replacing X by -A(B+S*)1 in the quadratic, and applying Equation (4.13) once, gives A + R*B + (R*) C = A - A(B+S*) B + A(B+S*)-A(B+S*) C = A - A(B+S*) B - A(B+S*)-S* = A[I- (B+S*)-1(B+S*)] = A[I- I] = Q, (4.14) proving that R* is a root of the quadratic. To prove that R* is in 6, we must prove that it is nonnegative, and that R*C1 = Al. The matrix S* is in /. Consequently, by Lemma 4. 2 and Lemma 3. 4, (B+S*) is a Q-matrix and its inverse is non-positive. Thus R* = - A(B+S*)1 is nonnegative. Replacing S by S* in Equation (4. 11), multiplying the right side of Equation (4. 11) by -A(B+S*)1 and multiplying the left side of (4. 11) by R* demonstrates that -R*C1 = -Al (4.15) and completes the proof. That such a fixed point exists will be proved using the wellknown theorem of Brouwer [18, pp. 42-45] which will not be proved here.

38 Theorem 4 2 (Brouwer's Fixed-Point Theorem) Let f be a continuous mapping of a non-empty, convex, compact subset A' of a finite-dimensional normed space 7/ into itself. Then f has a fixed point in It will be shown that a particular subset of / plays the role of' in this theorem, and that T plays the role of ~. When that can be done, we will have shown that a fixed point S* exists in the subset of <, (consequently, also in /), and that therefore a root R* = - A(B+S*) of the quadratic also exists. The set Ai of m x m matrices is a normed space (cf. [Fa, p. 58]), under the norm m IIMII max Z Imijl, (4.16) i=l,...,m j=l where M = (mij) is a matrix in. There is a class of subsets of 97t which are compact and convex. Let M be an m x m matrix. Let MM be the set of m x m matrices S which satisfy the restrictions S1 = Al (4.17) S > M > 0 (4.18) Theorem 4.3 demonstrates the compactness and convexity of cAM. Theorem 4. 3 For every nonnegative m x m matrix M, e'M is a compact, convex subset of )2.

39 Proof: Since CM is finite dimensional, in order to show compactness it is sufficient to show that it is closed and bounded. That it is convex is readily shown, since for each S1, S2 e <M, and each a e [0, 1], the matrix S = aS1 + (l-a)S2 has the properties SI = aS1 + (l-a)S1 = aAl + (l-a)Al = Al and S = aS1 + (l-a)S2 > aM + (l-a)M = M. That 2M is bounded is similarly easily shown, since every member of 6M is nonnegative and has the same row sums, and therefore every member of <M has the same norm. I. e. S1, S2 e M implies Ilsll1 = IIS2l = IIAl. o That AM is closed is proved by showing that if So = (s)ij is a limit point of.AM then it is in 4. First, it is noted that, for every s e J-M, S1 = Al. Consequently, any limit point So of <M clearly has S 1 = Al. Second, if S is a limit point but S > M, then there is at least one entry of S -M which is negative. Say, j - mi*j* = - 6 < 0. Then II (S-M)-(So-M) I > 6 for every S e AM' since (S-M) is nonnegative, and the absolute sum of the i*-th row is at least 6. Consequently, IIS-S ol > 6 for everySe jE M' and there is some e -neighborhood (e. g. e =6/2) of So which does not contain any member of ~M. This contradicts that So is a limit point, proves that

40 So > M and that S e GM' and completes the proof of the theorem. Having shown that aM qualifies as the subset for the Brouwer theorem for arbitrary M > Q, the next requirement is to examine the continuity of T. Theorem 4. 4 Let M be a nonnegative mx m matrix. If rM C / then T is a continuous mapping on 4M. Proof: By the definition of A, (B+S) is nonsingular for every S e CtM Let S be any member of aM n and consider any S in eM such that lS-Soll < 6 = e/IIA(B+So) 11. Then T(S) - T(So) = -A(B+S) C + A(B+So)1C = -A(B+So) [(!+So) (B+S) -I] C -A(B+So)[(B+So)-(B+S)](B+S) C =-A(B+So) (So-S)(B+S) C and IIT(S)-T(So)II < IIA(B+So)-11 IIS-Sojl i I3+S)-IlcJ < EII(B+S) -CI However, by Lemma 4.2, -(B+S) -C is a stochastic matrix. It follows that II (B+S)1C II = 1 and that IIT(S)-T(So)II < IIA(B+So)-l Is IS-Sll < which proves continuity at every SO e M concluding the theorem Several statements are useful alternatives to the antecedent of Theorem 4. 4.

41 Theorem 4. 5 Let M be a nonnegative m x m matrix. Then the following are equivalent: (a) (B+S) is nonsingular for every S E M. (b) eIMCC I When M1 < Al, these are equivalent to: (c) B+M has no essential index class whose corresponding entries to C1 are identically zero. Proof: The equivalence of (a) and (b) follows trivially from the definitions of I and CM. The equivalence of (c) is the important result. Since M is nonnegative, and M1 < Al, then (B+M)1 < -C1 and B+M is a Q-matrix. For every S e aC M, B+S> B+M, and both B+S and B+M are Q-matrices. Therefore, by Theorem 3. 3, every essential index class of B+S contains at least one essential index class of B+M. Thus condition (c) implies that for every S e eM' B+S has no essential index class whose corresponding entries to C1 are identically zero. But (B+S)l = -C1. Consequently, for every S e ciM' B+S has no such index class which is conservative. By Theorem 3. 1, B+S is nonsingular and thus (a) is implied by (c). Conversely, let D be a matrix formed by adding the entries of (A1-M1) to the corresponding diagonal entries of M. Then D1 = Al and D > M, so that D e /-M. By theorem 3.1 and Lemma 4.2, conditions (a) and (c) imply that B+D has no conservative index class, and hence no class for which corresponding entries of C1 are zero.

42 But the addition of a diagonal matrix does not change the communication properties of a matrix, hence B+M has the same properties as B+D. This proves that (a) and (b) imply (c), and completes the theorem. We are now prepared to routinely prove the principal result of this section: Theorem 4. 6 If a nonnegative matrix M exists such that (a) AM is not empty, (b) AM C (c) Tmaps M into M, then a root R of the quadratic A +XB + X2C = exists in. Proof: The conditions of the Brouwer theorem (Theorem 4. 2) are met for the normed space Ai, the nonempty, convex, compact subset <caM' and the continuous mapping T of CM into itself. Therefore there exists a fixed point S* of T, where S* e oM. Since eM C, then S* is also in, and by Theorem 4.1 the root R exists in A. This proves the theorem. The condition of Theorem 4. 6 is easily shown to be met by QBD-process for which B is irreducible, and for others where B has other readily identifiable properties. These classes of processes are quite useful, and are special cases of an even broader class which will be shown in Chapter 5 to also meet the conditions of

43 Theorem 4. 6. However, since the proof of the more general case is quite complicated, it will be instructive and somewhat satisfying to prove the special case now, as a simple corollary of Theorem 4.6. Corollary 1 If B has no essential index class whose corresponding entries to C1 are all zero, then a root R e R of the quadratic equation A + XB + X2C = Q exists. Proof: Choose M = Q. The set o contains the nonnegative matrix A, and is therefore not empty. Furthermore, it is clear that M1 < Al and thus, by Theorem 4. 5, po C J. Finally, by Lemma 4.2, T(S) > 0 and T(S)1 = Al for every S e co. Consequently, T maps 0 into 0 and the antecedent of Theorem 4.6 is satisfied, proving the corollary. A special case of QBD processes satisfying the antecedent of this corollary occurs when B is irreducible and C is not a zero matrix (Q). Such problems occur frequently in applications.

5. Boundary Leading QBD Processes The sufficient conditions for the existence of the quadratic root R which were defined by Theorem 4. 6 are satisfied by a very general type of QBD process, the "boundary leading" process. Definition A QBD process N is boundary leading if every nonboundary state leads to a boundary state. This condition is a very weak one, and encompasses all irreducible QBD-processes as well as a considerable variety of reducible ones. Transition graphs of a number of simple boundary leading QBD processes are shown in Figure 5. 1, demonstrating this variety. These graphs represent every state by a vertex, and every nonzero, non-diagonal transition intensity ujk (j [ k) by a directed branch from state j to state k. The states are arranged as in Figure 2. 1. Two way paths in Figure 5. 1 are shown by heavy lines, and boundary states by square vertices. These graphs are useful in visualizing the communication properties of particular QBD-processes, since a state j "leads to" a state k if and only if there is a directed path in the graph which can be followed from j to k. Note that only the graphs (e) and (f) describe an irreducible process. The rest are reducible. It should be particularly noted that the definition of boundary leading processes makes no mention of paths between boundary 44

45 (a) (d) (b) (e) (c) (f) ~(c) Figure 5. ) Figure 5. 1 Graphs of Boundary Leading QBD-processes

46 states. Thus, branches between square vertices in Figure 5.1 can be removed or added at will without altering the boundary leading character of the examples. This exercise offers still other interesting examples. The most interesting aspect of all the examples, however, is the way that the horizontal repetitiveness, characteristic of QBD-processes, forces certain communication properties. The most pervasive of these properties can be stated in the obvious "lemma". Lemma 5.1 Let j and k be two states, not both boundary states, of a QBD-process. Then j = < il, i2 > leads to k = < k, k2 > if and only if < j+l+, j2 > leads to < k-l+, k2 >, where f is a positive integer. (We omit proof.) The assumption that N is boundary leading is so general that it automatically includes every QBD process which has a finite number of communicating classes, as seen by the theorem: Theorem 5. 1 If a QBD process N has a finite number of communicating classes of states, then it is boundary leading. Proof Let nl {<nl, n2>: n2=1, 2,..., m} for each n=1, 2, Let 2K ={<0, n>: n2=1,2, 2,... m}. Suppose N is not boundary leading. Then there is a set of nonboundary states which do not lead to a boundary state. Let < al a2> be a member of that

47 set which has least value al > 0. It cannot lead to any state in T a, since those are either boundary states or lead to boundary states. The state <a1+1, a2> can not lead to any state in 7a, by Lemma 5. 1. Similarly states <a1+2, a cannot lead to states in I'C1, etc. Each of these states, therefore, cannot communiia+1' cate with one another, and must be in distinct communicating classes. This proves the existence of an infinite number of such classes, contradicting the antecedent and proving the theorem. The converse is not true, however, since Figure 5. l(d) describes a boundary leading process which (degenerately) has an infinite number of communicating classes.(Note that each state in the lowest row is a communicating class by itself.) A weaker form of the converse is somewhat interesting, and easily proved. Theorem 5. 2 If N is boundary leading, then N has a finite number of essential communicating classes. Furthermore, each such class includes at least one boundary state. Proof: By definition, every essential state communicates with every state it leads to. Since every state of a boundary leading process leads to a boundary state, every essential state communicates with a boundary state and every essential class contains a boundary state. No state may be in more than one class, so there must be no more essential classes than there are boundary states, a finite number. This proves the theorem.

48 This theorem shows that when a boundary leading process has an infinite number of communicating classes, the classes must necessarily be inessential (except for a finite number). A useful corollary shows that every boundary leading QBD process has some essential states. Corollary 1 If N is boundary leading, then N has at least one essential communicating class of states. Proof: It is sufficient, by the theorem, to show that at least one boundary state is essential. Every boundary state is either essential, or leads to another state which does not lead to it. If the latter, it must lead to another boundary state, because if the state it leads to is a non boundary state, that state must lead to a boundary state. Since the boundary states are finite in number, all cannot lead to another without communicating. Thus one must be essential, and the corollary is proved. The boundary leading QBD process has been singled out for study here because it represents one of the most general processes for which a quadratic root matrix R exists in 0. That the root in 6&adoes not exist in general is illustrated by a pair of examples. The first has an infinite number of essential communicating classes of states, and the second has an infinite number of inessential classes. (Recall that non-boundary leading processes must have an infinite number of communicating classes )

49 Example 5.1: Let A = B = C = E = 0 1 Lo -4_ 2 I_ 1 -2_j Then the graph of N becomes that of Figure 5. 2. Figure 5.2 Illustrating Example 5. 1 An equilibrium solution 7r can be expressed as = < f0' ~1, 0~, 2' 0,... > where 0', 1'2,'... are arbitrary positive constants. A bruteforce calculation of roots to A + XB + X C (found by expanding the matrix product and solving the resulting four equations in Aij) shows that every root is of the form R = 1 2 4 where ~ is an arbitrary constant. There RC1 = Al, and therefore no root R e -. Al = [1], B+S is necessarily singular for is therefore empty. is no root for which Furthermore, si nce every S E e, and j

50 Example 5. 2: Let 1 1 -2 0 -2 Q A =, B, C = E = -0 1_ 0 -1 1 -2_ so that N has the graph of Figure 5. 3. In this case the quadratic Figure 5. 3 Illustrating Example 5. 2 reduces to A + XB = 0, so that R = -AB is the only root:.5 1 R = 0 1_ Since C = 0, RCl is obviously not equal to Al, so that there are no roots in i. In this case, however, ~ is not empty, since it contains at least K0 2 S = These examples notwithstanding, boundary leadingness is not a necessary condition for the existence of a quadratic root in 1. This is attested by the further examples: Example 5. 3: Let 1 0~ -1 0 0 0 -2 1 A - B - C = E =, 1 0 -2_ _0 1 h 1 -2_ so that N has the graph of Figure 5. 4. Using the procedure outlined

51 Figure 5. 4 Illustrating Example 5. 3 by Gantmacher [11, pp. 227-230], all of the roots of A + XB + X2C are found to be: 01 0, [1 1!.1 2] R=:J LD U4l 2fl Of these, the second (and only the second) is in A, and for that one A^^ -2 2 E +RC = _ 1 -1 Thus, ro = < 1, 2>, and an equilibrium solution is = <1,2,3,3,6,6,12,2,12,... > We have also, the example Example 5. 4: Let 10 -1 0 0 -2 1 A = B = C= E = 1 - 0 - with the graph of N given in Figure 5. 5. In this case, using Gantmacher's procedure, we have the roots P1 0 5 1.5 R =,, 1~ 1- -~ t5 LO - 5J

52 Figure 5. 5 llustrating Example 5. 4 In this case, the last is in,. These latter examples appear to indicate that perhaps boundary leading should be extended to permit a "path" via infinity, so as to include these cases. However, the proof which follows uses only the definition originally given. We will show that for boundary leading QBD processes there exists a matrix V such that (a) V is not empty, (b) V C, and(c) T( V) CC V. This proof unfortunately requires considerable machinery for the construction of V. Let n be a positive integer, and let Un be the mn x mn matrix B A C B A Un =... (5.1) C B. This matrix is a Q-matrix. To construct V, we consider the inverse W(n) (assuming it exists) of Un. Let W(n) be partitioned so that its diagonal submatrices are m x m and m(n-l) x m(n-l)

53 Wl(n) W3(n)Un =W(n) = (5.2) W2(n) W4(n) We contend that the matrix B - AWl(n)C is a Q-matrix, and that for some n=a the nonnegative mxm matrix B - AWl(a)C has no essential index class whose corresponding entries to C1 are identically zero. (By Theorem 4. 4 this property is sufficient to show that B+S is nonsingular for all S e J.) We let -AW1(a)C V = -AWl(a)C, (5.3) and contend that <V has the desired properties. These contentions will be proved in a sequence of six lemmas. But first, observe that the matrix Un represents the transitions between the states in a subset ra of the state set of N, o, n where a is a positive integer and t(c - = {a,a+l,...,a+n-1} x {l,2,...,m}. Ca, n To aid visualization, the transitions of Un can be interpreted by a graph which is a "slice" of width n taken from the middle of the graph of N. (For example, the graphs of Figure 5. 6 interpret U3 for each corresponding graph of N in Figure 5.1.) With this background, we are ready to prove the first lemma:

54 (a) (b) (c) ( c ) (d) (e) Figure 5. 6 Graphs of U3 for Examples of Figure 5.1

55 Lemma 5. 2 If N is boundary leading, then Un is nonsingular for all n > 1. In particular, B is nonsingular. Proof: Suppose U is singular. Since it is a Q-matrix, by Theorem 3. 1 it has an essential index class which is conservative. This would imply that N has an essential class of states entirely within?an',since no transitions from states in 9 n are possible a, n a, n without producing row sums of U which are greater than the corresponding ones of Un. l does not contain any boundary states. nn o U(a,n Therefore there is an essential class which does not include a boundary state, contradicting Theorem 5.2. This proves the main consequent of the lemma. The remainder is proved by merely noting that U1 = B. Thus, the lemma is proved. A recursive relation for the matrices W1(n), n=l, 2,..., can be found by making use of the formula for inverting partitioned matrices. Lemma 5. 3 If N is boundary leading, then Wl(n) can be determined from the recursion Wl(n+l) = (B - AW(n)C), where W1(1) = B1 and n=1, 2,.... -1 Proof: That W1(1) = B1 is a trivial result of the definitions. When Un+I is partitioned so that its diagonal submatrices are m x m and mn x mn, the resulting matrix is written

56 A B A n Un+ (5.4) Cn Un A where An is the m x mn matrix n A A = [A... 0], (5. 5) ~~~~n and C is the mn x m matrix C 0 A Cn =. (5.6) Q But both B and Un are nonsingular, and the matrix Wl(n+1) is the upper diagonal submatrix of the inverse of Un+l. Therefore, by Lemma 3. 3, Wl(n+l) is given by the formula -1 -1A W(n+l) = (B-AnUn C) (5. 7) nnnn A A1 A for n=l, 2,... (the non singularity of B - AnU C being guarA A anteed by the lemma). In view of the zero elements of An and C nnn only the m x m upper diagonal submatrix of Un need be evaluated, and Equation (5. 7) can be simplified to W1(n+l) = (B - AW1(n)C)1 (5. 8) (Nonsingularity is, again, guaranteed.) This proves the lemma. The matrix Un is a nonsingular Q-matrix. Therefore its inverse is non-positive, and so also is its submatrix Wl(n). Furthermore, B - AW1(n)C is a Q-matrix, by the lemma:

57 Lemma 5. 4 If N is boundary leading, the following inequalities hold for all n=l, 2,... (a) -AWi(n)C > 0 (b) -AW1(n)C1 < Al, and therefore B - AWl(n)C is a Q-matrix. Proof: The non negativity of -AW1(n)C (inequality (a)) follows immediately from the non positivity of W1(n), and the non negativity of A and C. Recalling (Equation 4.3) that (A + B + C) = 0, it follows that - Cl - Al = B -1 and, multiplying on the left by AB -AB-1C1-AB-1A1 = Al. Since B is non positive, -AB Al is nonnegative, and therefore -AB Cl < Al. (5.9) Since W1(l) = B, this proves (b) for n=l. If - AW1(n)Cl < Al, then (B - AW1(n)C)l < (B + A)1, (B - AW1(n)C)l < - C1 -1 and, multiplying on the left by A(B - AW1(n)C), which is nonpositive Al > - A(B - AW(n)C) Cl = -AW(n+l)Cl.

58 This completes the induction for (b). The matrix B - AW1(n)C is the sum of a nonnegative matrix and a Q-matrix, so that its off-diagonal elements are nonnegative. Its row sums are non positive, so it is clearly a Q-matrix. This completes the proof. The next step is to show that for some a the matrix Wl(a) has no essential index class for which the corresponding entries of the vector Cl are zero. (This is needed in order to show that B+S is nonsingular for all S e cV' which in turn will be used to show that aV -C-. ) Our objective will be pursued by closer examination of the communication properties of Un. The matrix Wl(n) has nonzero entries wij(n) if and only if i leads to j under Un, where i and j are members of the first m indices of Un. This is a consequence of Theorem 3.2, Equations (5. 1) and (5. 2), and the fact (Lemma 5. 2) that Un is a nonsingular Q-matrix. Now, what we are interested in is whether or not all those first m indices of U lead to other indices whose entries in C1 are not zero. n To make this clearer, observe Figure 5. 7, which shows graphs of U1, U2, U3, and U4 for a particular N. The index?"4" is the only one of the first m (m=4) whose entry in C1 is not zero, since it has a leftward leading branch leaving it. In this graph it is seen that every other index in {1, 2, 3, 4} of U3 leads to index 4 (see the darkened paths), but that index 1 does not in U2, and indices 1 and 2 do not in U1. Thus, W1(3), has nonzero entries such that every index leads to an index whose entry in C1 is not zero.

2 3 4 U, U2 U3 Figure 5. 7 Graphs of Uk for a Particular N. U4

60 Let JC be the set of indices (in 1, 2,..., m) which have nonzero entries to C1. (I. e., in the above example J = {4}.) Let I(n) be the set of indices in 1, 2,..., m which do not lead (under U n) to any index in JC. Then we can generalize what we have just shown in example by the following lemma. Lemma 5. 5 If N is boundary leading, then there exists an integer a < m such that I(a+ 1) is empty. Proof: First, we note that if an index is a member of I(n+l) it must have been a member of I(n), since the increase of n to n+1 merely appears to add branches to the right of the graph to those already present in Un and can only add to the indices which lead to JC. Symbolically I(n+l) C I(n). Second, we note that if I(n+l) = I(n), then I(n+&) = I(n) for every positive ~, since: the increase of n to n+1 can also be regarded as an addition of a column of vectors and their branches to the left of the graph Un, and if I(n+l) = I(n), then no new paths are added. Consequently, the addition of a new column to the left of that can do no more than the last, because of the repetitive character of the graphs, etc. (Strictly we are here appealing to Lemma 5. 3). Third, we observe that if I(n+C) = I(n) for all positive l, then the indices in I(n) do not lead to an index in JC for any value of n. This implies that the states {< a, i >, ie I(n)} do not lead to any states which lead "leftward" from <a, and hence implies that the process N is not boundary leading. Consequently, I(n+l) C I(n) unless

61 I(n) = j. Finally, we observe that since I(1) could contain at most m members, and each successive set has less members, then the first set I(a+l) which is empty must occur before n=m, and thus a < n. This proves the lemma. The next lemma completes the construction of V. Lemma 5. 6 If N is boundary leading and a is the integer referred to in Lemma 5. 5, then B - AW1(a)C has no essential class whose corresponding entries to C1 are identically zero. Proof: The index i leads to index j under W (a+1) if and only if i leads to j under U + If it does, then, by the previous lemma, a+l' every index of W1(c+1) leads to one in JC. But by Equation (5. 8) Wl(a+l) = (B - AWl(a)C) 1 Again using Theorem 3. 2, the matrix B - AWl(a)C also has no index which does not lead to one in JC. Every index in an essential class must lead to one in JC' and the class must contain at least one. If an index in an essential class leads to another index, that latter index must be in the class. Therefore every essential class contains at least one index in JC' proving the lemma. One last lemma is required before the main theorem. It shows that {Wl(n), n=1, 2,... }, is a monotone sequence.

Lemma 5.7 If N is boundary leading, - Wl(n+l) > - Wl(n), for all n=1, 2,.... Proof: The lemma follows by induction. First, noting from Equation (5. 2) that W1(l) B1, we have W1(2) - W(l1) = Wi(2)[(W(1))1 - (W1(2)) ]W1(l) = W1(2)[B - (B-AB 1C)]W1(1) -1 = W(2)(AB C)W1(1) < 0 Second, if Wl(n) - Wl(n-1) < Q, then Wl(n+l) -wl(n) = Wl(n+l)[(B-AWl(n-l)C)-(B-AWl(n)C)]W1(n) = Wl(n+l)A(W1(n) - W1(n-1))C W1(n) _ 0 This proves the lemma. Theorem 5. 3 If N is boundary leading, a is the integer referred to in Lemma 5. 5,and V = - AW1(a)C, then (a) Iv <, (b) V C:, and (c) T( ) ) v. Proof: / is the set of all matrices S for which S1 = Al and S > V. By Lemma 5. 4, V1 = - AW(a)Cl < Al.

63 Therefore, there is a nonnegative matrix (e. g. the diagonal matrix whose diagonal elements are the entries of Al - VI) which, when added to V, produces a matrix in V. Thus (a) is true. By Lemma 5. 6, B+V has no essential index class for which the corresponding entries of C1 are identically zero. If S e ciV' then (B+S) > (B+V), and (Theorem 3. 3) every essential class of (B+S) contains at least one essential class of (B+V). But since (B+S)1 = C1, then no essential class of (B+S) is conservative. With Theorem 4. 5, this proves (b). To prove (c), we first observe that (Lemma 4. 3) for every S e 7, and therefore for every S E CV, T(S)l = Al. Second we observe (Lemma 5. 7) that - W1(c+l) > - W1(a), and that therefore - AW1(a+l)C > - AW1(a)C = V. Finally, we observe that for all S e < S - V> 0 and T(s) + AW1((+l)C = - A(B+S) C - A(B+V) 1C = - A(B+S)-[I - (B+S)-A(B+V)- ] = - A(B+S) 1[(B+V) - (B+S)](B+V)-1C = A(B+S) -(S-V)(B+V)- C > 0 and that, therefore T(S) - AW1(a+l)C > V. Thus, for all S e V T(s) e V. That completes the proof.

64 Theorem 5. 4 If N is boundary leading, then (a) There exists a fixed point S* of T, in V, and therefore also in ci C V' (b) There exists a quadratic root R e (, satisfying Equation (4. 8). Proof: This theorem follows immediately from Theorems 4. 1, 4.6 and 5.3.

6. Determination of R* It is, of course, not very useful to know that a root matrix R exists in L unless there is the means to calculate it. One such means, proposed by Evans, is to repeatedly apply the mapping T to some initial iterate So, hoping that the sequence T(S ), T (S ),..., will converge to a fixed point S*. Once S* has been determined, a root matrix R* is found by use of Theorem 4. 1 and the formula R* = - A(B+S*). (6. 1) The question naturally arises whether the iteration will converge to a fixed point S* of T, and for what initial iterates. Clearly, for T(So) to converge, T must not map any iterate T (S ) into a matrix S for which (B+S) is singular, for otherwise T +1(S ) = T(S) would be undefined and the sequence S, T(So), T (S ), T (S ),..., is without meaning. To avoid this difficulty, the initial iterate must be chosen from a subset, c a for which the mapping T has the property T( j ). The strongest convergence theorem we have succeeded in proving follows. In it we adopt the notation p(M) for the spectral radius of a matrix M. Theorem 6. 1 Let be a subset A C ~ for which T(g ). Let S* be a fixed point of T for which p(A(B+S*) 1) 1. Then 65

66 lim T(S) = S* n- oo for every So e. The following well known lemma is required for the proof of this theorem. The proof of the lemma will be found in [21, p. 13]. Lemma 6. 1 If M is an m x m complex matrix, then lim Mn = 0 n-to if and only if p(M) < 1. With this lemma, we proceed to the proof of Theorem 6. 1. Proof (Theorem 6. 1): The difference T(S) - S* is given by T(S) - S* = T(S) - T(S*) = - A(B+S)1C + A(B+S*)-1C - A[(B+S)1 - (B+S*) 1]C = - A(B+S*) [(B+S*) - (B+S)](B+S)- C, for every SE 6. Consequently, T(S) - S* [- A(B+S*) 1](S-S*)[- (B+S)-C] (6.2) and, by repeated application of Equation (6. 2) we have, for every SoE A, S 1 Tn S* = in n-1ink — T(S) - S* [-A(B+S*) ]n(So-S*) I [-(B4Tk (S)) C]. (6.3) k=0 However, since -(B+S)- C is a stochastic matrix for every S e A, -(B+T (S ))- C is a stochastic matrix for every k=0, 1, 2,..., and 0

67 II- (B+Tk(So))1CI = 1. Therefore the norm of the product on the right is less than or equal to unity for every n. Consequently, from Equation (6. 3) we get IITn(So) - s*lI < II[- A(B+S*)-l]nll IS0 - s*11. (6.4) But if p(A(B+S*) 1) < 1, then by Lemma 6. 1 lim[- A(B+S*)-]n 0, n-oo and thus lim IITn(S) -S*ll 0. (6.5) n -co For a finite dimensional normed space, convergence in norm is equivalent to convergence in components, so that the theorem is proved. Using Theorem 5. 3, we get the corollary Corollary 1 Let N be boundary leading. Then if a fixed Toint S* of T exists for which p(A(B+S*) ) < 1, we have lim T (S) = S* n -oo for every So E ~V. The theorem as given appears to be excessively restrictive. Experimentation with the iteration, with the purpose of discovering useful examples and counterexamples, produced remarkably consistent results. The only failures of convergence either occurred in non-boundary-leading processes or were the result of singularity

68 of B+Sn where Sn was some iterate. The latter were presumably the result of choice of initial iterate which was not in V' or some other suitable X. In some cases this was checked and -1 proved to be true. Specifically, the restriction "p(A(B+S*) ) < 1T in the theorem appears to be unnecessary. Convergence to rather large root matrices were rather frequently obtained, as a result of deliberate challenges to the program. However, the failure to find a stronger theorem, removing the spectral radius assumption, has not lead to serious difficulties in this theory. It will be seen (Chapter 8) that in almost all interesting cases, p(A(B+S*) ) < 1 holds. Of that, more later. There is one theoretical problem resulting from Theorem 6. 1 which is more troublesome. It is very awkward to work with the set V. A separate calculation of V is required and, with large m, it is rarely plain when o iterations have been completed and V has been obtained. One escape from this difficulty is to take one's chances with singularities in B+T n(S ) and simply restart with a perturbed T n(S ) when it occurs. This will seem more attractive with the observation that the interruption will only occur early in the iteration process. This will follow from the readily proved fact that: Lemma 6.2 If N is boundary leading, then Tn(S) > - AW1(n)C (6.6)

69 for every S e d and for every integer n for which Tn(S) is defined. Proof: The proof will be by induction on n. From Lemma 5. 1 we have Wl(n+l) =(B - AWi(n)C)-1 and Wl(1) = B Thus, -1 -1 T(S) + AWi(1)C = -A(B+S) C+ AB 1C = A[B- - (B+S)-1C =AB-1[(B+S) - B](B+S)-C = A[-B- S[-(B+S) 1]C Since every factor is a nonnegative matrix, we have T(S) +AW1(1)C > 0 and therefore T(S) > -AWi(1)C. For general n, we have, similarly, Tn+i(S) + AW (N+1)C = - A(B+Tn(S)) -1 AW(n)C)1 + A(B(- AWC(n)C) = A(B -AW1(n)C) [T (S) +AW(n)C](B+T (S))- C = A[- Wl(n+l)][T (S) + AWi(n)C][- (B+T (S))- ]. But [- W1(n+l)] is nonnegative, as is [- (B+Tn(S)) 1], and [T (S) + AW1(n)C] is nonnegative by the induction hypothesis. Therefore T (S) + AWl(n+l)C >

70 and the induction is complete, proving that Tn(S) > - AW1(n)C for all integers n for which Tn(S) is defined (i. e. all n such that Tk(S) e A fork=l, 2,...,n-1.) The set V is the set of matrices S' in < for which S' > V = - AW (a)C, where a is the integer defined in Lemma = 1 5. 5. Lemma 6. 2 shows that if N is boundary leading and if Tk(S ) is in J for every k=, 2,..., a, then Ta(So) is in V because T (S) V =- AWi(a)C. By Theorem 5.3 T( V) C V' so that Tk(So) must also be defined for k=u+l, (+2,... In short, if N is boundary leading, and the iteration "survives" a iterations without leaving the set c, then the iteration is "safe" from leaving J thereafter. Usually a is a very small integer, since it is at most equal to the number of essential classes of B. An alternative escape from the difficulty of finding V is to find another set which can play the role of j in Theorem 6. 1. It need only satisfy A D z ) T(4 ). Unlike:V' i' does not have to be convex, compact, or closed. For example, if N is boundary leading, the set {Tn(A), n=l, 2,...} can be shown to qualify, so that the iteration will remain in < if Note that this proof can also be used to prove that - AW1(n)C = Tn(0), if the domain of T is extended to the set {S: S > 0 and (B+S) is nonsingular}.

71 A is chosen as the initial iterate. Of course, convergence depends on whether the spectral radius property holds.

7. Equilibrium Distributions The construction of an equilibrium solution offers little theoretical difficulty. We have the theorem: Theorem 7. 1 If there exists a root R c /{ of the matrtix quadratic A +XB+X C, then (a) There exists a sequence 7T = {7k' k e 9 }, which satisfies 7T > 0 and 7TU = 0. (b) The sequence iT is of the partitioned form - = [0', 1',.. ] where 00 satisfies 00 > 0 and A A 0O(E + RC) = 0, (7. 1) and where A 01 ='R (7.2) On = n- R = 01(R) (7. 3) n=2, 3, 4,... (c) If p(R) < 1, then vT is convergent. The proof of this theorem requires the following lemma, due to Ostrowski. (See [21, p. 82] for proof.) Lemma 7. 1 Let M be an arbitrary complex matrix. Then the 2 series I + M + M +... converges if and only if p(M) < 1. If p(M) < 1, then (I-M) is nonsingular and the series converges to (I-M)1. 72

73 Having this lemma, we proceed with the proof. Proof (Theorem 7. 1): In the early discussion of Chapter 4, it was observed that if a root R of A + XB + X C = 0 exists in, then Equation (7. 1) has a solution <0 > 0. It was observed in Chapter 2 that the partitioned sequence iT = [, 1,... ] as defined by Equations (7. 1), (7. 2) and (7. 3) satisfies 7TU = 0, and that 7T> 0. Thus, conclusions (a) and (b) were already proven. Clearly, iT converges if 01 + 2 +. converges. Using Equation (7. 3) we have <1 + 02 +''' = 01 + 01R + <1R + * (7. 4) and, provided the series I + R + R2 +... converges, we can factor out the vector 01 so that 01 + 2 + * = 01(I + R+R +. (7. 5) Since by hypothesis p(R) < 1, then by Lemma 7. 1 the series I + R + R +... does converge. Therefore, fr is convergent, (c) is true, and the proof is complete. The converse of (c), i. e. that if vT is convergent then p(R) < 1, is not generally true for boundary leading QBD processes. Example 7. 1- The process whose transition graph is shown in Figure 7. 1 has a root 2 0 R = 0.5 so that p(R) = 2, and yet it has an equilibrium solution

74 IT = {0, 1,0,. 50,.25,...} which is convergent. 2 2 ~2 2 5.5.5.5 5 Figure 7. 1 A Graph of a Process for which p(R) > 1 and T is Convergent A A The matrix (E+RC) is a conservative Q-matrix for every R E I% (c. f. Chapter 4), and therefore Equation (7. 1) could be reA A garded as defining a finite Markov chain with E + RC as its transition intensity matrix. Consequently, any technique satisfactory for determining equilibrium solutions for finite Markov chains can be applied to determine a vector 00 satisfying (7. 1). Once 00 is determined, Equations (7. 2) and (7. 3) can be applied routinely to obtain a solution T, or as much of that solution as one wishes to expand. Theorem 7. 1 shows, further, that if p(R) < 1, then that solution is necessarily convergent. One technique for determining 00 from Equation (7. 1) follows that of the "Recursive Queue Analyzer" [23], and develops the iterations 1 = n[ (E + RC)+ I] (7. 6) where n=, 1, 2,, is a row vector with m entries, is a where fi, n=O, 1, 2,... is a row vector with m~ entries% A is a

75 positive scalar slightly greater than the absolute value of the larA A gest diagonal entry to E + RC, and g0 is an otherwise arbitrary rowvector for which 0 > 0. Proof that the iteration (7. 6) converges to an equilibrium solution 00 is found in [23]. (If the matrix A A E + RC is reducible, the 00 found is not unique and will depend on the choice of %.O Nevertheless, any sequence v constructed from any 00 is an equilibrium solution, so that the solution T is similarly not unique.) If p(R) < 1, then an equilibrium distribution can be found by normalization, through the use of the following theorem. Theorem 7. 2 If, under the conditions of Theorem 7. 1, p(R) < 1, then = [,0' 01'.. ] is an equilibrium distribution, where ~0 = %0/s (7.7) O = OR (7.8) - n-il(7.9) On = n-1 R = O1(R) (79) and where s is a scalar: (I-R) 1! s = jo (7.10) Proof: We show that S = Z k, (7.11) k E 9 so that it normalizes fT = IT/s to a distribution. The sum (7. 11) can

76 also be written k =k o + 01 +... = 01 + (01 +02 + ) (7.12) since it is (absolutely) convergent. By Equations (7. 2) and (7. 5) <0 o1 + ( + 2 +' -)1 = 0 1 + 0R (I + R + R2 +...)1! (7.13) The right hand side is equivalent to the partitioned form 1+ 22 2 1(R + R2 +.)1 (I+R+ + + ) 00 + 0 =0 1 o0 1 Finally, by Lemma 3. 5, we write (I+R+R2+...)1 (I- R) 1 ^0 = o o proving the theorem. It is possible to determine the equilibrium distribution more directly, while at the same time determining "marginal" equilibrium distributions on the second (finite) component of the process N = <N1, N2>. Let y be a row vector whose i-th entry is the sum of the i-th entries of all n' n=0, 1, 2,.... = + [ 00] + [0] + [~1 0] + 2 ~ + (7.14) We shall refer to y as the marginal equilibrium distribution on N2, since the sum of its entries is unity. Using Equations (7. 8) and (7. 9), and again using the absolute convergence of Equation (7. 4), it is easily verified that

77 y = 0o (I-R) o 0 0 I (7. 15) and that (I-R) is nonsingular. =0 = 7 Substitution of Equation (7. 16) partitioned so that E = where E1l is m x m and E22 i V(I-R)(E11 Equation (7. 15) implies that I-R 0 (7.16) 0 I into Equation (7. 1), and letting E be E11 E12 (7. 17) E21 E22 s(mo-m) x (m0-m), gives the equation - RC) E12 =0 (7.18) y (I-R)E21 E22It can be solved under the restriction y1 = 1, (7. 19) whereupon %0 is automatically found from Equation (7. 16). Usually, (7. 18) and (7. 19) are more easily solved than (7. 1), (7. 7) and (7. 10) since at least one matrix inversion is eliminated. Sometimes y is all that is desired. We express these results as the Theorem 7. 3 which has just been proved.

78 Theorem 7. 3 If, under the conditions of Theorem 7. 1, the spectral radius p(R) < 1, then for every equilibrium distribution 7T there is a marginal equilibrium distribution y which is a solution of (I-R)(E11-RC) E12 y = Q (7.18) (I-R)E2 E22 and yl = 1. (7.19) Furthermore, (I-R) 0 o 0 = Y It often happens that the equilibrium mean of N1 conditioned on N2 is desired. One can express this mean as an m-component rowvector y1 whose k-th entry is the equilibrium mean of N1 conditional on N2=k (Observe that for k > m, N1 can have only value zero): 71 = + 22 + 33 +... (I2R+3R2...). = 01(I + 2R + 3R +...). The calculation of y1 is facilitated by the theorem

79 Theorem 7. 4 Let M be an arbitrary complex matrix. Then the series I + 2M + 3M2 +... converges if and only if -2 p(M) = r < 1. The series converges to (I-M). Proof: This theorem is analogous to Lemma 7. 1. We first assume that p(M) < 1. We have the identity (I+2M+3M2+... + nMn-1(I-M)2 - I = [n(I-M) + I]Mn By multiplying both sides by [(I-M) -, (which we write (I-M) -2), we have (I+2M+M2... +nM ) -(I-M)2 = [n(I-M) +I]Mn(I-M)2 Thus l (I+.2M +3M2... + +nM 1) -(I-M) 211 < II [n(I-M) + I]MnI I (I-M)211, (7.17) for all r > 1. But I [n(I-M)+I]Mn I I = [(n+l)rn] and therefore the left hand side of (7. 17) converges to zero. Thus the series -2 converges (it is finite dimensional) to (I-M). Conversely, if the series converges, let / be an eigenvalue, corresponding to the eigenvector x. Then (I + 2M + 3M2 +... )x = (1 + 2 +32 +...)x, and the convergence of the series implies convergence of 2 3 ui + 2 + 3/3j +... for any eigenvalue of M. It is well known that this implies I tu < 1 for all eigenvalues of M, and thus that p(M) < 1, completing the proof.

80 As a consequence of this theorem we have a1 = E1(I-R)2 (7.18) and, from Equations (7. 2) and (7.16), we get R(I-R)1 Y1 = Y (7.19) 0 (recalling that R(I-R) = (I-R)R). The equilibrium second moment of N1 conditional on N2, which is correspondingly:'2 = \1 + 402 + 903 + *= p1(I+4R+9R2 +...) can similarly be calculated by the use of the theorem Theorem 7.5 Let M be an arbitrary complex matrix. Then the series I + 4M + 9M + 16M +... converges if and only if p(M) = r < 1. Then the series converges to (I+M)(I-M) 3 The proof of this theorem follows the same manipulations as previously, starting the with identity 3 2 2Mn-1 3 2 (I-M)3(I+ 4M + 9M2 +... + n2M ) - (I+M) [n3M - (n2+2n-1)I]Mn It need not be repeated here. Using this theorem, we have'2 = 01(I+R)(I-R)3 (7.20) and

81 R(I+R)(I-R)72 = 7 (7.21) (again recalling the commutativity of products of matrix polynomials).

8. Recurrent Processes A state j of a continuous-parameter Markov chain Z is said to be recurrent if o0 f p j(t)dt= oo, (8.1) 0 Where, by definition, Pjk(t) = pr[Z(t) = k Z(0) = j]. (8.2) The state j is said to be non-recurrent otherwise. A recurrent state is said to be positive recurrent if, in addition, lim p. (t) > 0, (8.3) t-oo JJc and null recurrent if lim p j(t) = 0. (8.4) t- oo J These properties are shared by all states ina communicating class. Consequently, we can say that a communicating class of states is non-recurrent, positive recurrent, or null recurrent if the states in in the class are non-recurrent, positive recurrent, or null recurrent, respectively. Every recurrent class is essential, and every finite essential class is positive recurrent. However, essential classes containing an infinite number of states are not necessarily positive recurrent, or even null recurrent. When a QBD process is used as a model of a physical system, the presence of an essential state (and therefore an essential class) 82

83 which is not positive recurrent signifies that the system is "unstable" in the sense that, with nonzero probability, the state of the system increases without bound as time elapses. This is a condition generally to be avoided in design of systems, and its mere occurrence generally makes further calculation unnecessary. However, the presence of an essential class which is not positive recurrent is often difficult to detect by casual observation or physical argument, and therefore an important capability of a computational technique is to provide such detection automatically. Concomitantly, it is not usually important that a computational technique provide solutions when such a class is detected. A key result of this section will be to show that a necessary and sufficient condition for an irreducible QBD process to be positive recurrent is that there be a fixed point S* of T for which p(A(B+S*) ) < 1. Since the latter was already shown (Chapter 6) to be sufficient for the convergence of the iteration process T to a unique fixed point S*, it follows that this procedure for calculating S* will work for every irreducible, positive recurrent QBD process, and that failure to converge, or convergence to a limit S* for which p(A(B+S*) ) > 1 is conclusive evidence that the process is not positive recurrent. Hence, we will have a test for positive recurrence, and also assurance that the computational technique works in every interesting case. For boundary leading QBD processes similar, but weaker, results will be demonstrated.

84 The following theorem, due to R. G. Miller [16] provides the key to (1) the relationship between positive recurrence of a continuous-parameter Markov chain and existence of convergent equilibrium solutions, (2) the uniqueness of equilibrium solutions, and (3) the relationship between the limiting distribution and the equilibrium solutions. It is an adaptation to continuous-parameter Markov chains of a very powerful theorem of F. G. Foster (see [10], and [17, p. 2561) for discrete-parameter Markov chains. The proof will be omitted here. Theorem 8.1 (R. G. Miller) If a continuous-parameter Markov chain Z is uniquely determined by its transition intensity matrix Q, and if Q is irreducible and conservative, then Z is positive recurrent if and only if -tQ = 0 has a convergent solution 7i > 0. The solution is unique (up to a multiplicative constant) and positive, and is equal to the limiting distribution (within a multiplicative constant). Specializing this theorem to QBD processes, we obtain the following corollary. Corollary 1 If N is an irreducible QBD process, then it is positive recurrent if and only if TrU = 0 has a convergent solution r > 0. The solution is unique (up to a multiplicative constant) and positive, and is equal to the limiting

85 distribution (within a multiplicative constant). Proof: Every continuous parameter Markov chain whose transition intensity matrix is bounded (i. e. whose entries are bounded) and conservative is uniquely determined by that intensity matrix. (See Chung [ 1, p. 237]). Since no element of U, the intensity matrix of N, can be greater in absolute value than the diagonal elements of E or B, and since there are finitely many of those, then the intensity matrix of N is bounded. By definition, every QBD process is conservative. Hence N is uniquely determined by its intensity matrix, and the antecedent of Theorem 8. 1 is satisfied, proving the corollary. Corollary 2 If N is an irreducible QBD process, then it is positive recurrent if and only if there exists a fixed point S* of T for which p(A(B+S*) ) < 1. Proof: If there exists such a fixed point S*, then by Theorem 4. 1, R* = - A(B+S*)-1 is a root of the quadratic A + XB + X2C. But then p(R*) < 1, and by Theorem 7. 1 a convergent sequence rT > 0 can be constructed which satisfies 7U = 0. By Corollary 1, N is positive recurrent. To prove the converse, we observe by Corollary 1 that if N is positive recurrent then there exists a convergent solution i > 0 to 7TU = 0, which is unique and positive. However, by Theorem 5. 4, there exists a root R e 1 of A + XB + X2C (recall that irreducible QBD processes are boundary leading). Theorem

86 7. 1 then tells us that iT = [ O, 1' 02'... ] is a solution to 7rU = 0, and 7r > 0. It must, therefore, be the unique, positive, convergent solution of Corollary 1. Moreover, since it is convergent oo oo - o00 n E_ On = nE >R 1 must converge. Since 1 is positive, 01R n=l n=l n=l can converge only if lim Rn =0 n - oo This implies that p(R) < 1 (cf. Lemma 6. 1), completing the proof. The preceding two corollaries give the desired assurances for irreducible processes. It is clear that if a QBD process is positive recurrent, the conditions for convergence of the iteration on T are satisfied, so that a solution can always be constructed. It is also clear that convergence of T and construction of a convergent 7T (or determination that IRI I < l)is sufficient assurance that the solution is the desired limiting distribution, and that the process is positive recurrent. Conversely, failure of T to converge (an occurrence which has not been observed to date), or divergence of the constructed 7T, is sufficient evidence that the irreducible QBD process is not positive recurrent, and that its limiting distribution is therefore pathological. These results extend to boundary leading QBD processes, with appropriate changes in terminology. First, for a model to be useful in estimating limiting probabilities, it is only necessary for all the essential states to be positive recurrent. Inessential states

87 are necessarily non-recurrent, but do not (at least in boundary leading processes) admit the possibility of tinstability" or "continued growth" referred to earlier. Inessential states always have a limiting probability of zero, and the equilibrium solutions reflect this property, as will be demonstrated. Secondly, when there is more than one essential class, the limiting distribution does not exist. Rather the limiting probabilities depend on the initial distribution. Thus, the equilibrium solution should not be expected to be unique. In general, an essential communicating class of states of a reducible QBD process does not separately define an irreducible QBD process, as is the case for Markov chains in general. Figure 8. la illustrates the transition graph of a boundary leading process whose subgraphs using any particular essential class of states is not a QBD process. Figure 8. lb shows a subgraph, and it is seen (a) 1 2 1 2 2 3 2 3 (b) Figure 8. 1 A Reducible QBD Process and One of its Subgraphs

88 that the transition intensities do not repeat properly. Therefore, the boundary leading QBD processes must be treated as a distinctive entity, rather than relying on the usual arguments for Markov chains which justify study of only irreducible processes by tsplitting up" the original process. The analogue of Corollary 1 can be expressed as follows. Theorem 8. 2 If N is boundary leading, then every essential state of N is positive recurrent if and only if U = 0 has a convergent solution = {7j: j e 9 } such that the components ur. are positive for every essential j and zero for every inessential j. Proof: Since the matrix U = {ujk: j, k e 2 } is conservative, then equivalently L Ujk = 0 je 2Z (8.11) k 7 Jk Then let 0 be the number of essential communicating classes in N and consider an arbitrary essential communicating class I, a=l, 2,..., 0. Since it is essential, then ujk = 0 whenever j e I jk - a and k k I. Therefore U jk O je Ia (8.12) a a and the matrix Q = {Ujk: j,k e I } is conservative. It is also aL jk a

89 irreducible (by definition of "communicating class"), bounded, and uniquely defines a Markov chain Z. By Theorem 8. 1 that Markov chain is positive recurrent if and only if there is a convergent vector v > 0 for which vQ =0. Now, suppose that irU = 0 has a convergent solution 7r = {Tr: j e 5f } for which Tj > 0 for every essential j and 7r = 0 for every inessential j. Then E 7TU, 0, kE u, (8.13) and since for every j a I either 7T. = 0 (inessential j) or ujk = 0 a jk for all k E I (essential j), then a 7Tjujk = O, ke Ia. (8.14) j 61 j jk This demonstrates that v = {7j: j Ia} is a positive solution to vQa = 0. The sequence fr is absolutely convergent, so that v is necessarily convergent. Thus the Markov chain Za is positive recurrent and v > 0. This, in turn, demonstrates that every state a in I is positive recurrent and a. > 0 for every essential j. (Recall a j that positive recurrence of j E I is defined by a property of pr[N(t) = j I N(O) = j], and that since every sample function of N for which N(O) = j must remain forever in I, then pr[N(t) = j | N(O) = j] = pr[Za(t) = j I Z (0) = j]. As a consequence, positive recurrence of Z implies positive recurrence of each a

90 j e I.) Every essential state is. in some essential class Ia, so that this proves the necessity of positive recurrence. Conversely, suppose every essential state is positive recurrent. Then every Z is also positive recurrent, and there exists a convergent solution v > 0 to v Q = 0. Write that solution as v =-{j: j 6 Ic} Then jUj = 0 ke I (8.15) jeI Now let I be the set of inessential states of N. Since uj = 0 for every j. I and k e I, 7Tjuj=0, ke. (8.16) jE j The sequence {7j: j e I} = 0 is surely a solution to Equation (8. 16). By Equations (8.15) and (8.16) we have demonstrated that 7T = {7tj: j E, } is a solution to rU=O for which ur. > 0 for every J J essential j, and 7i. = 0 for every inessential j. This completes the proof. It should be noted that the restriction that "ij. > 0 for every essential state" could be weakened to "Tf. > 0 for at least one state in each essential class" without affecting the proof. In that case, "t. > 0 for every essential state" becomes an additional conclusion J of the theorem, as it is in Corollary 1 of Theorem 8. 1. Every component solution v is unique within a constant. If the above a

91 qualification were not applied, it would be possible to demonstrate a convergent solution for a process which had a null-recurrent class I by the simple expedient of setting j. = 0 for every j E I. Example 7. 1 illustrates why the boundary leading analogue of Corollary 2 of Theorem 8.1 must be significantly weaker. The existence of a fixed point S* for which p(A(B+S*) 1) < 1 is clearly not necessary for the essential states of N to be positive recurrent. It appears, however, to be sufficient. The difficulty stems from the existence of inessential states of N which might have been essential with a different E-matrix. (Observe that properties of T, S* and R* do not depend at all upon E.) It appears that if every state of a boundary leading N is positive recurrent, then it is necessary that p(A(B+S*) ) < 1. But if every state is positive recurrent then every state is essential, so that we have a very different kind of a process. Weaker conditions that imply that p(A(B+S*) ) < 1 can be found. For example, it appears to hold if one assumes simply that "every infinite communicating class (is) positive recurrent. " (Proof of these suppositions requires development of relationships between communication properties of indices of S* and of boundary states of N, and will not be given.)

9. Uniqueness of S* The question of whether or not there is more than one root of 2 A + RB + R C in, is not central to this theory. The central issue is whether the equilibrium solution rT is unique, and whether it is given by the procedure of Chapter 7. If such a solution is unique, the existence of more than one R* which generates that solution is not disturbing. If the solution is not unique, it is similarly not too disturbing if the multiplicity comes from multiple solutions to the boundary equation (Equation 7. 1), or from multiplicity of R*. Nevertheless, a better intuitive feel for the mathematics of QBD-processes is a likely by-product of considering questions of uniqueness. Also, the answers to these questions may be useful in future development of QBD theory. Theorem 6.1 implies that if T has a fixed point S* e J for which p(A(B+S*) ) < 1, then there is no fixed point of ~, other than S*, in any set A for which T(<f ) C X C J. (Otherwise T (So) would not converge to S* for every SE. ) This result can be extended to the theorem: Theorem 9.1 If there exists a fixed point S* of T for which p(A(B+S*) ) < 1, then there is exactly one fixed point of T 92

93 Proof: Suppose that S' is a fixed point of T. Then T (S') = S' and the set 4 = {S'} has the property T ( )C, C g. Furthermore lim Tn(S') = S' n-coo However, by Theorem 6.1 lim Tn(S') = S* n —oo Therefore S' = S*, and the theorem is proved. Note that this theorem does not require any property of N other than that a fixed point of T exist with the given property. Even that property may be unnecessary, but a proof of a stronger theorem is not available. (An example of a QBD process with multiple fixed points of T has not been found, either.) It does not follow from the uniqueness of S* that there is, necessarily, a unique root of the quadratic A + XB + X C in the set X. The matrix R* = - A(B+S*)1 is, of course, a root. Furthermore, it can be readily shown that if R is a root and RC is in <I, then RC is a fixed point of T. (Thus S* = R*C.) However, since C is often singular, it is conceivable that R = R* is not the only root matrix for which RC = S*. Furthermore, if R is a root in 6i~, it does not necessarily follow that RC is in A, since B+RC may be singular. Thus, even if there is only one root with RC e J, there may still be more than one root in 6i.

10. Conclusions It has been shown that, under common and predictable circumstances, the QBD-processes lend themselves to expeditious numerical solution, and are therefore a useful class of stochastic models for computer system analysis (or other queueing application). The equilibrium distributions and other related measures are determined by a procedure which involves operations with only finite matrices, and which is analogous to procedures used in solution of the difference equation of simple birth and death processes. This conclusion has been demonstrated by a comprehensive algebraic theory whose two main theoretical conclusions can be summarized as: Conclusion 1. If a QBD process is boundary leading, then the matrix 2 quadratic A + XB +X C is, indeed, analogous to the characteristic equation for the birth and death process, having a root matrix R which generates the characteristic solution "f = 01 R n=2,3,... n 1 The m-vector ~1 satisfies a boundary condition %o(E + RC) = 0 A 01 = 0R and the equilibrium solution is given by 94

95 I: = [00 0' 0 * * *] The root matrix is in the known set &. If r is convergent, the equilibrium distribution is merely a scalar multiple of tr. Conclusion 2o If the QBD process is irreducible and positive recurrent, then the root R can be found by a well-defined iteration process S = T(Sn) R = C(B + S ) - A 00 and the resulting solution t7 will be convergent. Conversely, if the process is irreducible and a root R is found for which the resulting 7r is convergent, then the process is positive recurrent. Much more has been said, of course, but these results are the most significant and general. The theory presented is broad enough so that there are numerous specific questions about QBD processes which can be answered by intelligent use of the theorems, corollaries, and lemmas explicitly supplied. Thus, this theory provides a good base on which to build even greater usefulness. The most important area for further work is in the process for determining the root matrix. Firstly, it was observed that the interaction on T has been found, empirically, to converge even in nonrecurrent irreducible, and in reducible, but boundary leading, QBD

96 processes. Both of these cases do not have the advantage of valid proofs of convergence (except, in the latter case, when the spectral radius property p(A(B+S) ) < 1 holds). Stronger proofs should be available. Secondly, it does not appear likely that the iteration T, requiring, as it does, an inversion on every iteration, is the most efficient procedure possible for finding the root R. A case in point is the iteration, Rn+ - DB (A + CRn) n+1 n n which has been studied extensively. (Dn is a diagonal matrix chosen to keep Rn1 in t. ) A series of 1800 cases with diverse choices of A, B, and C have been run, and this iteration has converged in every boundary leading case. No proof that this should be true has been found. Nevertheless, the inverse of B can be calculated once, and no matrix inversion is required thereafter. Thus, this iteration is probably more efficient than the iteration of T, which requires inversion on every iteration. In the last analysis, when one considers the variety of birth and death processes which have been studied, and the results concerning them which are available, one realizes the potential for future work that is here. How much of birth and death theory can be generalized by the act of replacing scalar transition-intensities by submatrices? This work has probably only scratched the surface in providing the answer. However, it should provide a firm base from which to generalize and develop a larger theory.

References [1] Chung, K. L., Markov Chains with Stationary Transition Probabilities, (Hamburg: Springer-Verlag, 1960). [2] Collatz, L., "Aufgaben monotoner Art, " Arch. Math., Vol. III, 1952, pp. 366-376. [3] Disney, R. L., "Some Multichannel Queueing Problems with Ordered Entry, " J. Industrial Engineering, (Technical Notes) Vol. XIII, No. 1, January-February, 1962. [4] Evans, R. V., "The Structure of Some Two-Dimensional Queueing Systems, " Research Report 85, Management Sciences Research Project, University of California, Los Angeles, July 1963. [5] Evans, R. V., "Numerical Methods for Queues with Complex Service Rates, " Research Report 82, Management Sciences Research Project, University of California, Los Angeles. [6] Evans, R. V., "Geometric Distribution in Some Two-Dimensional Queueing Systems, " Operations Research, Vol. 15, No. 5, September-October 1967, pp. 830-846. [7] Fadeeva, V. N., Computational Methods of Linear Algebra (Transl. by C. D. Benster), (New York: Dover, 1959). [8] Fife, D. W. and Rosenberg, R. S., "Queueing in a Memory Shared Computer, " Proc. of the 19th National Conf. of the ACM Philadelphia, 1964. [9] Foley, J. D., "A Markovian Model of The University of Michigan Executive System," Comm. of the ACM, Vol. 10, No. 9, September 1967. [10] Foster, F. G., "On the Stochastic Matrices Associated with Certain Queueing Processes, " Ann. Math. Stat., Vol. 24, 1953, pp. 355-360. [11] Gantmacher, F. R., The Theory of Matrices, Vol. 1, (New York: Chelsea, 1959). [12] Gantmacher, F. R., Applications of the Theory of Matrices (Transl. by J. L. Brenner), (London: Intersciences, 1959). 97

98 [13] Householder, A. S., Principles of Numerical Analysis, (New York: McGraw-Hill, 1953) - [14] Kleinrock, L., "Time Shared Systems: A Theoretical Treatment," J. ACM, Vol. 14, No. 2, April 1967, pp. 242-261. [15] Krishnamoorthi, B. and Wood, R. C., "Time Shared Computer Operations with Both Interarrival and Service Times Exponential, " J. ACM, Vol. 13, No. 3, July 1966, pp. 317-338. [16] Miller, R. G., "Stationary Equations in Continuous Time Markov Chains, " Technical Report No. 80, Applied Math. and Stat. Labs, Stanford University, Stanford, California. [17] Parzen, E., Stochastic Processes, (San Francisco: HoldenDay, 1962). [18] Saaty, T. L., and J. Bram, Nonlinear Mathematics, (New York: McGraw-Hill, 1964). [19] Smith, J. L., "An Analysis of Time-Sharing Computer Systems Using Markov Models, " Proc. AFIPS 1966 Spring Joint Computer Conference, Vol. 28, pp. 97-104. [20] Taussky, 0., "A Recurring Theorem on Determinants, Amer. Math. Monthly, Vol. 56, pp. 672-676. [21] Varga, R. S., Matrix Interactive Analysis, (Englewood Cliffs, New York: Prentice-Hall, 1962). [22] Vulikh, B. Z., Introduction to Functional Analysis, (Transl. by I. N. Sneddon), (Oxford: Pergamon Press, 1963). [23] Wallace, V. L. and Rosenberg, R. S., "RQA-1, The Recursive Queue Analyzer, " Technical Report No. 2, Systems Engineering Laboratory, The University of Michigan, Ann Arbor (1966). [24] Wallace, V. L. and Mason, D. L., "The Optimal Degree of Multiprogramming in Page-on-Demand Computer Systems," Comm. of the ACM, Vol. 12, No. 5, May 1969.

UNIVERSITY OF MICHIGAN IIII 9II0 IIIII IE 10II 3 9015 03627 7450