# Finite-State-Space Truncations for Infinite Quasi-Birth-Death Processes

Authors: Hendrik Baumann
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

For dealing numerically with the infinite-state-space Markov chains, a truncation of the state space is inevitable, that is, an approximation by a finite-state-space Markov chain has to be performed. In this paper, we consider level-dependent quasi-birth-death processes, and we focus on the computation of stationary expectations. In previous literature, efficient methods for computing approximations to these characteristics have been suggested and established. These methods rely on truncating the process at some levelN, and for N⟶∞, convergence of the approximation to the desired characteristic is guaranteed. This paper’s main goal is to quantify the speed of convergence. Under the assumption of an f-modulated drift condition, we derive terms for a lower bound and an upper bound on stationary expectations which converge quickly to the same value and which can be efficiently computed.


---

## Body

## 1. Introduction

In many applications, discrete-time or continuous-time Markov chains are used to model real-life problems. For these models, characteristics of interest can be determined, enabling us to understand relationships between parameters and characteristics, solve optimization problems, and so on. When Markov models are used, the characterization of the dynamics by its probability transition matrixP or generator matrix Q, respectively, is quite simple due to all theoretical problems being solved. However, in applications, the determination of interesting characteristics of the process is very important. In many models, there is no chance of finding explicit analytical representations of these characteristics, and hence, we have to use numerical calculations or simulation methods. In most practical situations, interesting characteristics refer to stationary probabilities or stationary expectations which can be derived from the invariant distribution. In queueing theory, elementary examples are given by moments of the stationary number of customers in a queueing system or queueing network, or by the stationary probability that the number of customers exceeds some threshold. Similarly, such characteristics are interesting in Markovian population models, epidemic models, etc. Formally, stationary expectations are given by πf=∑xϵEπxfx, where E is the state space of the Markov chain, π is the invariant distribution, and f is some cost or reward function.The invariant distribution for an irreducible and positive recurrent Markov chain is characterized by the unique probability vector which solves a certain homogeneous system of linear equations. The number of states coincides with the number of equations and the number of variables. Unfortunately, in realistic models, we have a very large number of states or even infinitely many states. In these situations, numerically computing stationary probabilities or stationary expectations requires a truncation of the system of equations (which means a truncation of the state space), and this truncation results in inevitable errors.Hence, numerical calculations only compute some approximationA to the desired characteristic πf, and from both a mathematical and a practical point of views, bounds on the truncation error πf−A are interesting. In this paper, we will derive lower and upper bounds on πf which is equivalent to computing an approximation to πf and bounds on πf−A. Since all values are computed numerically, our results and methods provide a posteriori error estimates. Competing methods for finding a posteriori error estimates were given in [1–3], and an approach for finding a priori estimates can be found in [4, 5]. In Section 8, we will briefly compare these methods to the method suggested in this paper.Although our theoretical results will be quite general, we will only discuss the efficient implementation of the computation of the bounds for level-dependent quasi-birth-death processes (LDQBDs). These processes are characterized by the block-tridiagonal structure of the transition probability matrix (discrete-time) or generator matrix (continuous time). This structure is very typical for a large class of queueing models, population and epidemic models, etc.The rest of the paper is organized as follows: After introducing some notation in Section2, we will present our theoretical results on lower and upper bounds on stationary expectations πf in Section 3. In Section 4, we will derive a method for computing these bounds efficiently for LDQBDs. Afterwards, we will apply our method to the elementary M/M/1 queue (Section 5), a variant of the M/PH/1 queue (Section 6), and a popular retrial queueing model (Section 7). Whereas we have an explicit analytical representation of the invariant distribution and of stationary characteristics in the first example, numerical computations are the only chance to find stationary expectations in the latter two examples. Finally, we give a comparison of the suggested method to competing ones (Section 8).


## 2. Preliminaries

### 2.1. Notations and an Auxiliary Result

In order to introduce our goals formally, we review some well-known facts on the limit behaviour of Markov chains. In this context, and in the whole paper, we will use the following notations:(i)
I is a finite or infinite identity matrix. The dimension will become clear from the context(ii)
1 denotes a finite or infinite (column) vector with all entries being 1. Alternatively, 1 refers to a function with constant value 1. The meaning will be clear from the context(iii)
1A denotes the indicator function of set A, that is, 1Ax=1 if x∈A, and 1Ax=0 otherwise. With slight abuse of notation, we use 1B for Boolean expressions B. If B is true, 1B is 1, and otherwise, 1B is 0(iv)
id is the identity function(v)
We use the notationℕ=1,2,3,⋯ and ℕ0=0,1,2,3,⋯.Furthermore, we will use a probabilistic proof of invertibility frequently throughout the paper: IfP is finite and transient, I−P is invertible. We give some details: A finite substochastic matrix P=piji,j=1m is referred to as transient if limn→∞Pn⟶0 (component-wise). This is true if and only if P contains no stochastic submatrix piji,j∈F for some F⊂1,⋯,m. Hence, P is transient if and only if, for all i∈1,⋯,m, there are i0,i1,⋯,ik∈1,⋯,m with i0=i,pir−1,ir>0 for r=1,⋯,k, and ∑j=1mpik,j<1.Due toP being finite, Pn⟶0 entails λ<1 for all eigenvalues of P, and hence, ∑n=0∞Pn converges, and by a standard argument, we obtain
(1)I−P∑n=0∞Pn=limN→∞I−P∑n=0NPnlimN→∞I−PN+1=I,which means that I−P is invertible.


### 2.2. Results on the Limit Behaviour of the Markov Chains

We briefly summarize some results on the limit behaviour of Markov chains since these results are the reason why it is important to be able to compute (approximations to) stationary expectationsπf.We consider Markov chainsXnnϵℕ0 in discrete time or Markov chains Ytt≥0 in continuous time with countable state space E. Throughout the paper, we assume that Xnn∈ℕ0 is irreducible and recurrent with transition probability matrix P=pxyx,y∈E or that Ytt≥0 is nonexplosive, irreducible and recurrent with generator matrix Q=qxyx,y∈E. Then, there is an invariant measure ψ, that is, a positive vector ψ with ψP=ψ or ψQ=0, respectively. In both cases, ψ is unique up to constant multiples. For f:E⟶ℝ, we introduce the notation ψf≔∑xϵEψxfx. If ψf and ψg converge absolutely with ψg≠0 for functions f,g:E⟶ℝ, a fundamental result on the limit behaviour of Markov chains is
(2)limn→∞∑k=0n−1fXk∑k=0n−1gXk=ψfψgorlimn→∞∫0tfYsds∫0tgYsds=ψfψg.in the sense of almost sure convergence (see [6], pp. 85-86, 203-209). If we have the stronger assumption of positive recurrence, the sum of the entries of ψ is finite. Hence, by multiplying with an appropriate constant, we obtain the uniquely determined invariant distribution π which still is an invariant measure, but additionally satisfies ∑xϵEπx=1. Then, we can set ψ=π and g=1 constant, and obtain
(3)limn→∞1n∑k=0n−1fXk=πforlimt→∞1t∫0tfYsds=πf,almost surely if πf converges absolutely. If f measures costs or rewards that are associated with the time spent in state x, πf is the long-run average of costs or rewards per time unit. In particular, for A⊂E, we can specify
(4)fx=1Ax=1,x∈A,0,x∉A,and π1A is the long-run average of the proportion of time spent in set A. For A=x, π1x=πx is the long-run average of the proportion of time spent in state x. Note that there is another interpretation of the entries of the invariant distribution π: For irreducible and positive recurrent discrete-time Markov chains which are additionally aperiodic, Xn converges to some random variable X∞ in distribution, and for nonexplosive, irreducible, and positive recurrent continuous-time Markov chains, Yt converges to some random variable Y∞ in distribution. Then, π is the distribution of X∞ or Y∞ (see any textbook on Markov chains, e.g., [6, 7]), that is, ℙX∞=x=πx and ℙY∞=x=πx. Hence, if πf converges absolutely for some function f:E⟶ℝ, we have πf=EfX∞ or πf=EfY∞, respectively.Due to these results on the limit behaviour, in many applications, all characteristics of interest can be deduced fromψfℓ for functions f0,⋯,fL, where we will set f0=1 in case of positive recurrence. Furthermore, by splitting functions into positive and negative part, we can assume fℓ≥0 for all ℓ=0,⋯,L without restriction. We remark that ψf is vector valued, ψf=ψf0,⋯,ψfL.


### 2.3. Notations for Block-Structured Markov Chains

Throughout this paper, let the state spaceE=E0∪E1∪E2∪⋯ be partitioned into finite and pairwise disjoint sets Ei which will be referred to as levels. We introduce some notation:
(i)
We fixK∈ℕ0. K will be chosen such that the f-modulated drift condition (see below) holds on EN+1∪EN+2∪⋯. Furthermore, throughout this paper, x0∈EK is some arbitrary, but fixed state within level EK(ii)
ψ=ψxx∈E denotes the invariant measure with ψx0=1. In case of positive recurrence, π=πxx∈E is the invariant distribution(iii)
Fori,j∈ℕ0, Pij or Qij, respectively, shall denote the transition probabilities or transition rates, respectively, for transitions from states ∈Ei to states ∈Ej.(iv)
Fori∈ℕ0, we write ψi=ψxx∈Ei,πi=πxx∈Ei,fi=fℓxx∈Ei0≤ℓ≤L and gi=gℓxx∈Ei0≤ℓ≤L for functions f,g:E⟶ℝL+1.Later on, multidimensional functionsf will be useful for computing approximations to several stationary expectations at the same time. With this notation, we have ψf=∑i=0∞ψifi, and we are interested in computing approximations to ψi and use them for approximating ψf, or we want to find an approximation to ψf directly by other means. As in other approaches for finding error bounds on the approximation of stationary characteristics, we assume that an f-modulated drift condition is satisfied, that is,
(5)∑j=K+1∞Pijgj−gi≤−fi,i>K,(6)∑j=K+1∞Qijgj≤−fi,i>K.In scalar notation, (6) reads as ∑y∈EK+1∪EK+2∪⋯qxygy≤−fx for xϵEK+1∪EK+2∪⋯. f-modulated drift conditions are a popular tool to prove positive recurrence forfx=1 and convergence of ψf or πf. In our examples in Sections 5–7, we will use this tool, too. For details, we refer to Theorem 1 and Theorem 2 in [8]. Note that the idea of making use of f-modulated drift conditions for this purpose is much older (see, e.g., [9–13], but some of these classical references formally require fx=1, whereas the results in [8] only require fx≥0.Finally, throughout the paper, we will require that(7)Pij=0fori>K>jorQij=0fori>K>j.We will comment this structural requirement in Section8.


## 2.1. Notations and an Auxiliary Result

In order to introduce our goals formally, we review some well-known facts on the limit behaviour of Markov chains. In this context, and in the whole paper, we will use the following notations:(i)
I is a finite or infinite identity matrix. The dimension will become clear from the context(ii)
1 denotes a finite or infinite (column) vector with all entries being 1. Alternatively, 1 refers to a function with constant value 1. The meaning will be clear from the context(iii)
1A denotes the indicator function of set A, that is, 1Ax=1 if x∈A, and 1Ax=0 otherwise. With slight abuse of notation, we use 1B for Boolean expressions B. If B is true, 1B is 1, and otherwise, 1B is 0(iv)
id is the identity function(v)
We use the notationℕ=1,2,3,⋯ and ℕ0=0,1,2,3,⋯.Furthermore, we will use a probabilistic proof of invertibility frequently throughout the paper: IfP is finite and transient, I−P is invertible. We give some details: A finite substochastic matrix P=piji,j=1m is referred to as transient if limn→∞Pn⟶0 (component-wise). This is true if and only if P contains no stochastic submatrix piji,j∈F for some F⊂1,⋯,m. Hence, P is transient if and only if, for all i∈1,⋯,m, there are i0,i1,⋯,ik∈1,⋯,m with i0=i,pir−1,ir>0 for r=1,⋯,k, and ∑j=1mpik,j<1.Due toP being finite, Pn⟶0 entails λ<1 for all eigenvalues of P, and hence, ∑n=0∞Pn converges, and by a standard argument, we obtain
(1)I−P∑n=0∞Pn=limN→∞I−P∑n=0NPnlimN→∞I−PN+1=I,which means that I−P is invertible.


## 2.2. Results on the Limit Behaviour of the Markov Chains

We briefly summarize some results on the limit behaviour of Markov chains since these results are the reason why it is important to be able to compute (approximations to) stationary expectationsπf.We consider Markov chainsXnnϵℕ0 in discrete time or Markov chains Ytt≥0 in continuous time with countable state space E. Throughout the paper, we assume that Xnn∈ℕ0 is irreducible and recurrent with transition probability matrix P=pxyx,y∈E or that Ytt≥0 is nonexplosive, irreducible and recurrent with generator matrix Q=qxyx,y∈E. Then, there is an invariant measure ψ, that is, a positive vector ψ with ψP=ψ or ψQ=0, respectively. In both cases, ψ is unique up to constant multiples. For f:E⟶ℝ, we introduce the notation ψf≔∑xϵEψxfx. If ψf and ψg converge absolutely with ψg≠0 for functions f,g:E⟶ℝ, a fundamental result on the limit behaviour of Markov chains is
(2)limn→∞∑k=0n−1fXk∑k=0n−1gXk=ψfψgorlimn→∞∫0tfYsds∫0tgYsds=ψfψg.in the sense of almost sure convergence (see [6], pp. 85-86, 203-209). If we have the stronger assumption of positive recurrence, the sum of the entries of ψ is finite. Hence, by multiplying with an appropriate constant, we obtain the uniquely determined invariant distribution π which still is an invariant measure, but additionally satisfies ∑xϵEπx=1. Then, we can set ψ=π and g=1 constant, and obtain
(3)limn→∞1n∑k=0n−1fXk=πforlimt→∞1t∫0tfYsds=πf,almost surely if πf converges absolutely. If f measures costs or rewards that are associated with the time spent in state x, πf is the long-run average of costs or rewards per time unit. In particular, for A⊂E, we can specify
(4)fx=1Ax=1,x∈A,0,x∉A,and π1A is the long-run average of the proportion of time spent in set A. For A=x, π1x=πx is the long-run average of the proportion of time spent in state x. Note that there is another interpretation of the entries of the invariant distribution π: For irreducible and positive recurrent discrete-time Markov chains which are additionally aperiodic, Xn converges to some random variable X∞ in distribution, and for nonexplosive, irreducible, and positive recurrent continuous-time Markov chains, Yt converges to some random variable Y∞ in distribution. Then, π is the distribution of X∞ or Y∞ (see any textbook on Markov chains, e.g., [6, 7]), that is, ℙX∞=x=πx and ℙY∞=x=πx. Hence, if πf converges absolutely for some function f:E⟶ℝ, we have πf=EfX∞ or πf=EfY∞, respectively.Due to these results on the limit behaviour, in many applications, all characteristics of interest can be deduced fromψfℓ for functions f0,⋯,fL, where we will set f0=1 in case of positive recurrence. Furthermore, by splitting functions into positive and negative part, we can assume fℓ≥0 for all ℓ=0,⋯,L without restriction. We remark that ψf is vector valued, ψf=ψf0,⋯,ψfL.


## 2.3. Notations for Block-Structured Markov Chains

Throughout this paper, let the state spaceE=E0∪E1∪E2∪⋯ be partitioned into finite and pairwise disjoint sets Ei which will be referred to as levels. We introduce some notation:
(i)
We fixK∈ℕ0. K will be chosen such that the f-modulated drift condition (see below) holds on EN+1∪EN+2∪⋯. Furthermore, throughout this paper, x0∈EK is some arbitrary, but fixed state within level EK(ii)
ψ=ψxx∈E denotes the invariant measure with ψx0=1. In case of positive recurrence, π=πxx∈E is the invariant distribution(iii)
Fori,j∈ℕ0, Pij or Qij, respectively, shall denote the transition probabilities or transition rates, respectively, for transitions from states ∈Ei to states ∈Ej.(iv)
Fori∈ℕ0, we write ψi=ψxx∈Ei,πi=πxx∈Ei,fi=fℓxx∈Ei0≤ℓ≤L and gi=gℓxx∈Ei0≤ℓ≤L for functions f,g:E⟶ℝL+1.Later on, multidimensional functionsf will be useful for computing approximations to several stationary expectations at the same time. With this notation, we have ψf=∑i=0∞ψifi, and we are interested in computing approximations to ψi and use them for approximating ψf, or we want to find an approximation to ψf directly by other means. As in other approaches for finding error bounds on the approximation of stationary characteristics, we assume that an f-modulated drift condition is satisfied, that is,
(5)∑j=K+1∞Pijgj−gi≤−fi,i>K,(6)∑j=K+1∞Qijgj≤−fi,i>K.In scalar notation, (6) reads as ∑y∈EK+1∪EK+2∪⋯qxygy≤−fx for xϵEK+1∪EK+2∪⋯. f-modulated drift conditions are a popular tool to prove positive recurrence forfx=1 and convergence of ψf or πf. In our examples in Sections 5–7, we will use this tool, too. For details, we refer to Theorem 1 and Theorem 2 in [8]. Note that the idea of making use of f-modulated drift conditions for this purpose is much older (see, e.g., [9–13], but some of these classical references formally require fx=1, whereas the results in [8] only require fx≥0.Finally, throughout the paper, we will require that(7)Pij=0fori>K>jorQij=0fori>K>j.We will comment this structural requirement in Section8.


## 3. Upper and Lower Bounds onψf: Theoretical Results

Our goal is to find a lower bound and an upper bound onψf such that
(i)
Both bounds converge toψf as the truncation level N tends to ∞(ii)
Both bounds can be computed numerically in an efficient manner.In this section, we focus on the “mathematical” criteria, that is, we derive exact terms for the lower and upper bounds in which both converge toψf. The latter goal of computational efficiency will be considered in Section 4 for the special case of QBDs.First, we focus on the discrete-time setting. We use the notation for block-structured Markov chains as introduced above, and we begin with a representation of invariant measures and their approximations. For this purpose, we introduce the hitting timesσN=infn∈ℕ:Xn∈EK∪EN+1∪EN+2∪⋯, τ=infn∈ℕ:Xn∈EK, and τ0=infn∈ℕ:Xn=x0.Lemma 1.
LetXnn∈ℕ0 be an irreducible recurrent discrete-time Markov chain with block-structured transition probability matrix P as in the preliminaries, let K,N∈ℕ0 with K≤N and x0∈EK. Let
(8)SKN=I,SjN=∑i=0NSiNPij,j∈ℕ0∖K,and TN=∑i=0NSiNPi,K.
(a)
SiN is well defined, and we have
(9)SiN=E∑n=1σN1Xn=yX0=xx∈EKy∈Ei,i≠K,TN=E∑n=1σN1Xn=yX0=xx,y∈EK.
TN and all SiN increase monotonically in N.
(b)
Si∞=limn→∞SiN exists for all i∈ℕ0 with SK∞=I and
(10)Si∞=E∑n=1τ1Xn=yX0=xx∈EKy∈Ei,i≠K.
T∞=limn→∞TN exists with
(11)T∞=E∑n=1τ1Xn=yX0=xx,y∈EK.
For the limits, we haveSi∞=∑i=0∞Si∞Pij for allj≠K and T∞=∑i=0∞Si∞PiK.
(c)
T∞ is stochastic and irreducible. Let ψK=ψxx∈EK be the invariant measure for T∞ subject to ψx0=1, and let ψi=ψKSi∞ for i≠K. Then ψ=ψ0,ψ1,⋯ is the unique invariant measure for P subject to ψx0=1.(d)
We haveψf=ψKF where F=∑i=0∞Sifi.Proof.
We refer to the representations ofSiN, TN, Si, and T as “probabilistic interpretations.”
(a)
The definition ofSiN implies
(12)SiNi=0i≠KNI=Piji,j=0i≠KN=PKjj=0j≠KN.
SincePiji,j=0i≠KN is finite and transient, I=Piji,j=0i≠KN is invertible, and hence, SiN is well defined for i=0,⋯,K−1,K+1,⋯,N. Afterwards, the definition of SiN for i>N is explicit. We omit the proof of the probabilistic interpretation since such considerations are standard in the context of finding probabilistic interpretations for invariant measures (see, e.g., [7, 14]), and since these considerations are similar to those in the proof of Lemma 3.2b) or Lemma 4 below. Note that the probabilistic interpretation directly implies the monotonicity since σN increases monotonically.
(b)
More precisely,σN converges to τ almost surely. By monotone convergence, the probabilistic interpretation of Si∞ and T∞ follows. Note that these expectations are finite due to recurrence. Si∞=∑i=0∞Si∞Pij for j≠K and T∞=∑i=0∞Si∞PiK can be proved by means of the probabilistic interpretations. Again, we omit all further details.(c)
The probabilistic interpretation ofT∞ can be rewritten as
(13)T∞=ℙτ<∞,Xτ=yX0=xx,y∈EK.
Due to recurrence, we haveτ<∞ almost surely, and stochasticity and irreducibility follow easily. Since T∞ is finite, it admits an invariant measure ψK which becomes unique by requiring ψx0=1. Multiplying Si∞=∑i=0∞Si∞Pij by ψK immediately leads to ψj=∑i=0∞ψiPij for i≠K, and analogously, we find ψK=ψKT∞=∑i=0∞ψiPiK. Hence, ψ=ψ0,ψ1,⋯ is invariant for P. Note that we can find the probabilistic interpretation
(14)ψy=E∑n=1τ01Xn=yX0=x0,by standard considerations (again, we omit details). This is the standard representation of invariant measures with ψx0=1 (see, e.g., [7, 14]).
(d)
The statement is a direct consequence ofψf=∑i=0∞ψifi and the representation of ψi found in (c).
BothψK and F depend on Pij for all i,j∈ℕ0. Our goal is to find bounds both on ψK and F which can be computed from S0N,⋯,SNN, since these matrices only depend on the finite matrix PN=Piji,j=0N. We begin with bounds on ψK.Lemma 2.
Let all requirements of Lemma1 hold and choose A=AN=axyx,y∈EK such that AI−TN is an invertible diagonal matrix.
(a)
Such a matrixA exists, and for two such matrices A,A~, we have A~=DA with an invertible diagonal matrix D.(b)
For any such choice ofA, we have
(15)φ¯x≔ax0,xax0,x0≤ψx≤axxax0,x0≕φ¯x,x∈EK.(c)
Letφ¯xN=φ¯x and φ¯xN=φ¯x be defined as in (b). Then φ¯xN increases monotonically in N,φ¯xN decreases monotonically in N, and both converge to ψx for all x∈EK.Proof.
We adapt ideas for proving classical bounds on quotients of entries of the invariant measure (see, e.g., [15]).
(a)
The probabilistic interpretation ofTN implies that TN is transient. Hence, I−TN is invertible and we could choose A=I−TN−1. If both AI−TN and A~I−TN are invertible diagonal matrices, so is A~I−TNAI−TN−1=A~·A−1.(b)
First, we note thatφ¯x and φ¯x do not depend on the choice of A since a left-hand multiplication by some diagonal matrix has no impact on the quotients under consideration. Now, we will directly construct a matrix A with axx=1 for all x∈EK. With TN=txyx,y∈EK and T∞=txy∞x,y∈EK, this works as follows: We set
(16)λxy1=txy,λxyn=∑z≠xλxzn−1tzy,n≥2,μxy1=txy∞,μxyn=∑z≠xλxzn−1tzy∞,n≥2,axy=∑n=1∞λxyn,gxy=∑n=1∞μxyn,y≠x.Additionally, we setaxx=gxx=1. Then for all, x∈EK,axyy∈EK is the standard construction (see, e.g., [7, 14]) of the minimal subinvariant measure for the substochastic matrix T subject to axx=1. Indeed,
(17)∑z∈EKaxztzy=txy+∑z≠x∑n=1∞λxzntzy=λxy1+∑n=2∞∑z≠xλxzn−1tzy=∑n=1∞λxyn=axy,for y≠x. Hence, AI−TN is a diagonal matrix. Analogously, the rows of G=gxyx,y∈EK are subinvariant measures for T∞. This latter matrix is finite and stochastic, and hence, every subinvariant measure is invariant. In particular, every row of G is a constant multiple of ψK, that is,
(18)ψyψx=xy,for all x,y∈EK. Due to txy≤txy∞, we have λxyn≤μxyn by a trivial induction, and thus, axy≤gxy. From axx=1, we obtain
(19)axyaxx≤ψyψx,for all x,y∈EK. By setting x=x0, we obtain φ¯y≤ψy (due to ψx0=1), and by setting y=x0, we obtain 1/φ¯x≤1/ψx.(c)
We havelimn→∞TN=T∞, implying limn→∞λxyn=μxyn, and by monotone convergence, we obtain
(20)limn→∞axyaxxgxy=ψyψx.
By settingy=x0 or x=x0, respectively, the statement follows.
Lemma2 gives a lower and an upper bound on all entries of ψK. It remains to find bounds on F. A lower bound is a direct consequence of Lemma 1.Lemma 3.
Let all requirements of Lemma1 hold, let f≥0, and let
(21)F¯N=∑i=0NSiNfi.(a)
We haveF¯N≤F.(b)
F¯N increases monotonically in N with limit F.As pointed out above, for finding an upper bound, we have to use some information onpxy with x∈EN+1∪EN+2∪⋯ or y∈EN+1∪EN+2∪⋯. This is done by using the f-modulated drift condition (5). Note that the function f is given, and the function g solving the drift condition can often be found by “educationally guessing” even if there is no chance to find explicit analytical terms for ψ or ψf.Lemma 4.
Let all requirements of Lemma1 hold, let the f-modulated drift condition (5) hold, and let the structural requirement (7) be fulfilled. Furthermore, define f¯i=fi+∑j=N+1∞pijgjfori≤NandF¯N=∑i=0NSif¯i.
(a)
ForN≥K, we have F¯N≥F.(b)
Additionally, letψg<∞. Then F¯N converges to F.Proof.
(a) First, remember(22)F=∑i=0∞Si∞fi=fK+∑i=0∞i≠KSi∞fi=fℓx+∑y∈E∖EKΨxyfℓyx∈EK0<ℓ≤L,where for y∉<EK(23)Ψxy=E∑n=1τ1Xn=yX0=x=E∑n=0τ−11Xn=yX0=x,implying
(24)Ψxyfℓy=E∑n=0τ−11Xn=yfℓyX0=x.
Fory∈EK, we have
(25)E∑n=0τ−11Xn=yfℓyX0=x=δxyfℓx,and together, we obtain
(26)F=∑yϵEE∑n=0τ−11Xn=yfℓyX0=xxϵEK0≤ℓ<L=E∑n=0τ−1fℓXnX0=xxϵEK0≤ℓ<L.
Analogously, we find(27)F¯N=E∑n=0σNf¯ℓXnX0=xxϵEK0≤ℓ<L.
Hence, it remains to showE∑n=0τ−1fℓXnX0=x≤E∑n=0σN−1f¯ℓXnX0=x for all
x∈EK. For this purpose, we first write
(28)E∑n=0τ−1fℓXnX0=x=E∑n=0σN−1fℓXnX0=x+E1σN<τ∑n=σNτ−1fℓXnX0=x=E∑n=0σN−1fℓXnX0=x+E∑k=0σN−11Xk+1∈EN+1∪EN+1∪⋯∑n=k+1τ−1fℓXn∣X0=x.
Letg∞ℓx=E∑n=0τ−1fℓXn∣X0=x. By means of the tower property of conditional expectations, we find
(29)E∑k=0σN−11Xk+1∈EN+1∪EN+2∪⋯∑n=k+1τ−1fℓXn∣X0=x=E∑k=0σN−11Xk+1∈EN+1∪EN+2∪⋯E∑n=k+1τ−1fℓXn∣Xk+1∣X0=x=E∑k=0σN−11Xk+1∈EN+1∪EN+2∪⋯g∞ℓXk+1∣X0=x=E∑k=0σN−1E1Xk+1∈EN+1∪EN+2∪⋯g∞ℓXk+1∣Xk∣X0=x=E∑k=0σN−1∑y∈EN+1∪EN+2∪⋯pXk,yg∞ℓy∣X0=x.
Together, we obtain(30)E∑n=0τ−1fℓXn∣X0=x=E∑n=0σN−1fℓXn+∑y∈EN+1∪EN+2∪⋯pXk,yg∞ℓyX0=x,for the entries of F. Below, we will show g∞ℓx≤gℓx for x∈EK+1∪EK+2∪⋯, and due to f¯x=fx+∑y∈EN+1∪EN+2∪⋯pxygy and N≥K, the desired statement follows.
For provingg∞ℓx≤gℓx, we define gMℓx=E∑n=0maxτ−1,MfℓXnX0=x for x∈EK+1∪EK+2∪⋯. Then g0ℓx=fx≤gℓx for x∈EK+1∪EK+2∪⋯, and by induction, we obtain
(31)gMℓx=fx+E1X1∈EK+1∪EK+2∪⋯∑n=1maxτ−1,MfℓXnX0=x+E1X1∈E0∪⋯∪EK+1∑n=1maxτ−1,MfℓXnX0=x=fℓx+E1X1∈EK+1∪EK+2∪⋯∑n=1maxτ−1,MfℓXnX1X0=x+0=fℓx+E1X1∈EK+1∪EK+2∪⋯gM−1ℓx1X0=x≤fℓx=fℓx+∑y∈EK+1∪EK+2∪⋯pxygM−1ℓy≤gℓx,for all x∈EK+1∪EK+2∪⋯ and all M∈ℕ0, since transitions from EK+1∪EK+2∪⋯toE0∪⋯∪EK−1 have probability 0. Due to recurrence, τ is finite almost surely, implying maxτ−1,M⟶τ almost surely, and by monotone convergence, gMℓx⟶g∞ℓx, implying g∞ℓx≤gℓx for all x∈EK+1∪EK+2∪⋯ as desired.
(b) SinceF¯N converges to F, it suffices to prove that F¯N−F¯N converges to 0. For this purpose, we write
(32)F¯N−F¯N=∑i=0NSiNfi¯−fi=∑i=0NSiN∑j=N+1∞Pijgj=∑j=N+1∞∑i=0NSiNPijgj=∑j=N+1∞SiNgj≤∑j=N+1∞Sj∞gj.
Again, letSj∞=Ψxyy∈Ejx∈EK. For x=x0, we have Ψxz≤ψz (due to the respective probabilistic interpretation, see proof of Lemma 1 for the interpretation of ψz), and since we assume that ψg converges, we have convergence to 0 for N⟶∞. Since the choice of x0∈EK is arbitrary, this is true for any x∈EK.
By summarizing all previous results, we obtain our main result for the discrete-time setting where we will omit the dependency ofSiF¯,F¯,⋯ on N. Furthermore, we will phrase the result directly for both discrete-time and continuous-time setting. All previous results may be transferred into the context of continuous-time Markov chains easily by means of the embedded jump chain. Since similar considerations have been used in the literature frequently (see, e.g., [4, 16, 17]), we omit further details.Theorem 5.
LetX=XnnϵN0 be an irreducible and recurrent discrete-time Markov chain, or let Y=Ytt≥0 be a nonexplosive, irreducible, and recurrent continuous-time Markov chain, respectively. Let the transition probability matrix P=Piji,j=0∞ or Q=Qiji,j=0∞, respectively, be block-structured as introduced in the preliminaries. Let the drift condition (5) or (6), respectively, and the structural requirement (7) hold. Furthermore, let
(i)
Si∈ℝdK×di be defined by SK=I and Si=∑j=0NSjPji or 0=∑j=0NSjQjiresp.,i=0,⋯,K−1,K+1,⋯N(33)T=∑j=0NSjPjKorT=−∑j=0NSjQjK,resp.(ii)
A=axyx,y∈EK be chosen such that AI−T or AT, respectively, is a diagonal matrix, let φ¯x=ax0,x/ax0,x0 and φ¯x=axx/ax,x0 for x∈EK, and let φ¯=φ¯xx∈EK andφ¯=φ¯xx∈EK as row vectors(iii)
∑i=K+1∞Pijgj+fi≤gi or ∑i=K+1∞Qijgj+fi≤0, resp., for i≥K+1(iv)
f¯i≔fi+∑j=N+1∞Pijgj or f¯i≔fi+∑j=N+1∞Qijgj, resp., be component-wise finite for all i=0,⋯N(v)
and(34)F¯=∑i=0NSifiandF¯=∑i=0NSif¯i.(a)
Si is uniquely defined for i=0,⋯N, and I−T or T, resp., is invertible, in particular, there is a matrix A with the desired properties(b)
We haveφ¯F¯≤ψf≤φ¯F¯(c)
Ifψg<∞, we have
(35)limN→∞φ¯F¯=ψF=limN→∞φ¯F¯.Remark 6.
Theorem5 is phrased for Markov chains with an infinite state space, and the goal is to find approximations to ψf which rely on the transition probabilities for transitions within a finite subset of the state space. Since even finite state spaces can be so large that an exact computation of ψf cannot be performed, we might want to apply Theorem 5 to Markov chains with a finite state space in order to reduce the state space to a smaller one.
Indeed, forE=E0∪⋯∪EM, K≤N<M, ∑i=K+1MPijgj+fi≤gi for i=K+1,⋯,M, and f¯i=fi+∑j=K+1MPijgj, (a) and (b) in Theorem 5 remain true. Note that formally, we cannot allow N=M since SK=Hxyx,y∈Ga would become stochastic, and I−T would not be invertible anymore. However, N=M means that we have no reduction of the state space, and hence, this setting makes no sense. Part (c) becomes obsolete.


## 4. Efficiently Computing the Bounds for Quasi-Birth-Death Process (Discrete-Time Setting)

We turn now to developing methods for computing both bounds onψf efficiently (and simultaneously). As pointed out in the introduction, we focus on quasi-birth-death processes (QBDs). QBDs are characterized by Pij=0 or Qij=0, resp., for i,j∈N0 with ∣i−j∣≥2, that is, each jump changes the level at most by 1.Level-independent QBDs were analysed by Neuts [18] by means of matrix-geometric methods. First algorithmic approaches for level-dependent QBDs (LDQBDs) are due to Bright and Taylor [17] and Hanschke [19]. The approach in [17] generalizes Neuts’ probabilistic interpretations of the matrices which arise in the matrix-geometric method, whereas the approach in [19] is motivated by the relationship between second-order vector-matrix difference equations and matrix-valued continued fractions. Remarkably, both methods are equivalent (up to suggestions regarding some initializations). More details and comparisons can be found in [20, 21]. All these methods intend to find (approximations to) an invariant measure or the invariant distribution. In [22], an algorithm was suggested which allows to compute stationary expectations directly. In the continuous-time setting, this method reads as follows:
(i)
Choose largeN, set R=0 and F=fN(ii)
Fori=N−1,N−2,⋯,0, replace R by Qi,i+1−Qi+1,i+1−RQi+2,i+1−1 and F by fi+RF(iii)
Determineψ0 as an (approximative) solution to x−Q00−RQ10=0(iv)
Returnψ0F.Since the memory requirement does not depend directly onN (with di=∣Ei∣, it depends on max{di:0≤i≤N}), the truncation level N can be chosen very large. Despite the possibility to choose large N, results on the truncation error still remained desirable from both a mathematical and a practical points of view.It is not difficult to prove that (with the notation of the present paper) the method from [22] computes the lower bound on ψf for K=0 (K has no impact on the lower bound). The goal of this section is to generalize this method in such a way that any value K∈0,⋯,N is allowed, and that the upper bound will be computed, too. For means of conciseness, we focus on the continous-time setting.Note that the QBD property directly impliesQij=0 for i>K>j, that is, (7) holds. The drift condition (6) is met if
(36)Qi,i−1gi−1+Qiigi+Qi,i+1gi+1+fi≤0,i≥K+1,and the definition of f¯i simplifies to
(37)f¯i=fii<N,fN+QN,N+1gN+1,i=N.Most importantly, the computation of the matricesSi according to the system
(38)0=∑i=0NSiQij,j≠K,becomes much easier (the matrix-analytic methods in, e.g., [17, 19–21] use this fact). Many of the following identities are consequences of the probabilistic interpretation of the involved matrices, and we omit these proofs since the considerations are similar to those in [17] anyway.
(i)
LetRK,⋯,RN−1 be defined by RN−1=QN−1,N−QNN−1 and
(39)Ri−1=Qi−1,i−Qii−RiQi+1,i−1K<i<N.Then, the inverses in (39) exist, and we have
(40)Si=RKRK+1⋯Ri−1,i=K+1,⋯N(ii)
Similarly, forB1,⋯,BK defined by B1=Q10−Q00−1 and
(41)Bi+1=Qi+1,iI−Qii−BiQi−1,i−1,0<i<K,we have
(42)Si=BKBK−1⋯Bi+1,i=0,⋯,K−1(iii)
The structure ofQ allows to write T=−QKK−BKQK−1,K−RKQK+1,K.If the level sizesEi are relatively small, the computational effort induced by (39) and (41) is acceptable. Then, it seems natural to compute all Ri and Bi, then all Si and finally F¯ and F¯. However, in [22], it was pointed out that for K=0, computing F¯ can be performed in a much more efficient way (in particular with respect to memory requirement) by using a Horner-type scheme. Here, we generalize this procedure slightly with respect to two issues: We consider an arbitrary K, and we consider the simultaneous computation of F¯ and F¯. Introduce
(43)Zn=fn+∑m=n+1N∏r=nm−1Rrfm,n=K,K+1,⋯N,Z¯n=f¯n+∑m=n+1N∏r=nm−1Rrf¯m,n=1,⋯K.Wn=∑m=0n−1∏r=0n−1−mBn−rfm,n=1,⋯,K.With this notation, we have(44)F¯=WK+ZK,F¯=WK+Z¯K,ZN=fN,Zn=fn+RnZn+1,n=K,K+1,⋯,N−1,Z¯N=f¯N,Z¯n=f¯n+RnZ¯n,n=K,K+1,⋯,N−1,W1=B1f0,Wn=BnWn−1+fn−1,n=2,⋯,K−1.The recurrence schemes forRn, Zn, and Z¯n all start at n=N, and are used in the sequel for computing the other values for n=N−1,N−2,⋯,K. Similarly, the computation of Bn and Wn starts at n=1, and then the values for n=2,3,⋯,K can be computed. For finding φ¯,φ¯,F¯,F¯, we only need the values for n=K. Hence, for n>K, Rn is only used for finding Zn, Z¯n, and Rn−1. Similarly, the other matrices are only used in a single step. Hence, there is no need to store these matrices. Instead, we suggest to use the recurrence scheme for all these matrices as an “update” procedure. In total, we suggest the following method.
(i)
ChooseN, set R=0, Z=fN, and Z¯=fN+QN,N+1gN+1.(ii)
Fori=N−1,N−2,⋯,K, update
(45)−R=Qi,i+1,−Qi+1,i+1−RQi+2,i+1−1,−Z=fi+RZ,andZ¯=fi+rZ¯.(iii)
SetB=Q10−Q00−1, and W=Bf0.(iv)
Fori=2,3,⋯,K, update
(46)−B=Qi,i−1−Qi−1,i−1−BQi−2,i−1−1,–andW=Bfi−1+W.(v)
SetT=−QKK−BQK−1,K−RQK+1,K and determine A=axyx,y∈EK such that AT is a diagonal matrix, define φ¯ by φ¯x=ax0,x/ax0,x0, and define φ¯ by φ¯x=ax,,x/ax,x0.(vi)
Returnφ¯W+Z as a lower bound and φ¯W+Z¯ as an upper bound.Usually the matricesQij and fi can be generated when they are needed. Up to four of these matrices are needed at the same time, and we need memory for saving R,Z,B,W. In total, a small number of finite matrices have to be stored at the same time. In particular, if ∣di∣=Ei is bounded by d, the memory requirement is ≤10d2+5dL+1, and this bound does not depend on K or N. Note that the “price” for this low memory requirement is that we only calculate the values ψfℓ for cost/reward functions fℓ which have to be specified before the computation procedure begins. Since we do not store the values ψx or the matrices Si, adding any “new interesting” function requires a complete restart of the method.A discussion of all numerical details of the algorithm is beyond the scope of this paper. For two specific aspects (avoiding ill-conditioned problems when computingϕ and avoiding instabilities in the update step for B), we refer to the Appendix.


## 5. Application to theM/M/1 Queue

We first consider an example where we have explicit terms forψx,Hxy,⋯. Of course, we would never use numerical methods for finding bounds on ψf in such a situation, but the following considerations illustrate how the method works, and they show how sharp the upper bound can be. Numerical examples for situations in which we do not have an explicit analytical representation for ψx will follow in Section 6 and Section 7.Consider theM/M/1 queue, that is, we have E=N0 and
(47)Q=−λλμ−λ+μλμ−λ+μλ⋱⋱⋱.λ>0 is referred to as arrival rate, and μ>0 is called service rate. Q is bounded, and therefore, nonexplosive. For any λ,μ>0, Q is irreducible, and we will assume positive recurrence which is equivalent to ρ≔λ/μ<1.Here, we haveEi=i, and we will first consider choices of g which allow choosing K=0.The invariant measureψ with ψ0=1 is given by ψi=ρi, and the 1×1-matrices Si shall solve 0=∑j=0NSjQji for i>0 which reads as
(48)0=Si−1λ−Siλ+μ+Si+1μ,i=1,⋯,N−1,0=SN−1λ−SNλ+μ.This system of equations can be solved explicitly (e.g., by using standard results concerning linear difference equations with constant coefficients), and we find(49)Si=ρi−ρN+11−ρN+1,i=0,⋯,N.Due to∣E0∣=1, we do not have to find the upper and lower bounds on ψxx∈E0=1, and instead, we focus on bounds on F.First, considerfi=1. As pointed out above, we have positive recurrence, and therefore, ψf should be finite. From the explicit term, we obtain ψf=1/1−ρ<∞. To show the finiteness of ψf by means of the drift criteria, set gi=1/μ−λ. Then, for all i≥1, we have
(50)qi,i−1gi−1+qiigi+qi,i+1gi+1+fi=μi−1−λ+μi+λi+1μ−λ+1=λ−μμ−λ+1=0.By Theorem 2 in [8] (or older references, see remark in the preliminaries), we obtain positive recurrence. Furthermore, the requirements of Theorem 5 are met. Therefore, we have F≤ψf≤F¯, where
(51)F¯=∑i=0NSifi=∑i=0Nρi−ρN+11−ρN+1=11−ρ−N+1ρN+11−ρN+1,F¯=F¯+SNqN,N+1gN+1=F¯+ρN−ρN+11−ρN+1⋅λ⋅N+1μ−λ=11−ρ−N+1ρN+11−ρN+1+ρN1−ρ1−ρN+1⋅ρN+11−ρ=11−ρ.Hence, the lower bound converges toψf=1/1−ρ, and the upper bound coincides with ψf which is the best possible result. It is clear that this can only happen if we have ∑qijgj+fi=gi (instead of mere ≤ for all i. Consider now ψf where f=1A, that is, fx=1 for x∈A and fx=0 otherwise. Obviously, we can still use the same function g, and with the same computation as before, we obtain
(52)F¯=∑i=0Nρi−ρN+11−ρN+11Ai,F¯=∑i=0Nρi−ρN+11−ρN+11Ai+ρN+1N+11−ρN+1.LetA=0,⋯,M, and let N>M. If we are interested in computing the bounds on ψ1A/ψ1, we combine the previous bounds, and the numerical algorithm would return
(53)∑i=0Nρi−ρN+1/1−ρN+11Ai1/1−ρ≤ψ1Aψ1≤∑i=0Nρi−ρN+1/1−ρN+11Ai+ρN+1N+1/1−ρN+11/1−ρ−N+1ρN+1/1−ρN+1.For any choice ofA, both bounds will converge to 1−ρ∑i∈Aρi which is the stationary probability for set A. As an easy choice, let A=j for some j∈ℕ0. Then, we know that the stationary probability is given by πj=1−ρρj, and for N≥j, we obtain the bounds
(54)1−ρ1−ρN+1ρj−ρN+1≤πj≤1−ρ⋅ρj−ρN+1+ρN+1N+11−ρN+1−1−ρN+1ρN+1.Next, letf=id, that is, fi=i for all i∈ℕ0. Note that ψid/ψ1 is the expected stationary number of customers in the system which is one of the most important performance measures in queueing theory. We set gi=i2/2μ−λ+iμ+λ/2μ−λ2, and then we have
(55)qi,i−1gi−1−qiigi+qi,i+1gi+1+fi=μi−12−λ+μi2+λi+122μ−λ+μ+λμi−1−λ+μi+λi+12μ−λ2+i=−2μ−λi+μ+λ2μ−λ+−μ+λμ−λ2μ−λ2+i=0.The exact value isψf=ρ/1−ρ2, the computed lower bound is given by
(56)F¯=f0+∑i=1NSifi=∑i=1Nρi−ρN+11−ρN+1⋅i=ρ1−ρ21−ρN+1NρN+1−N+1ρN+1−ρN+1NN+121−ρN+1=ρ1−ρ2+ρ1−ρ21−ρN+1NρN+1−N+1ρN+ρN+1−ρN+1NN+121−ρN+1=ψf−ρN+1N+11−ρ1−ρ21−ρN+1−ρN+1NN+121−ρN+1=ψf−ρN+1N+121−ρ1−ρN+12+N1−ρ,and the upper bound is given by
(57)F¯=F¯+SN⋅λ⋅gN+1=F¯+ρN−ρN+11−ρN+1ρN+1221−ρ+ρN+11+ρ21−ρ2=F¯+ρN+11−ρN+121−ρ21−ρN+1N+11−ρ+1+ρ=ψf.Hence again, we find the best possible upper boundF¯=ψf. Appropriately combining the bounds on ψ1 and ψid leads to the upper and lower bounds on the expected stationary number πf of customers in the system.At a first glance, the choicegi=i2/2μ−λ+iμ+λ/2μ−λ2 might be “difficult to guess,” but note that we have used a function g of the form gi=ci for dealing with fi=1, and it is quite natural to choose a function g of the form gi=c1i2+c2i for dealing with fi=i. Afterwards, it is not difficult to determine c1,c2 such that the f-modulated drift condition is met. Nevertheless, in other situations, we will not be able to find such an optimal function g. Therefore, we demonstrate next that the method still works if we use some function g which is far from being optimal.Let us directly considerf=1,1A,id, that is, L=2,f0i=1,f1i=1Ai, and f2i=i. Then, ψf1/ψf0 is the stationary probability for set A, and ψf2/ψf0 is the expected stationary number of customers in the system. Let us choose gℓi=ai for i∈ℕ0 and ℓ=0,1,2. Then
(58)qi,i−1gℓi−1+qiigℓi+qi,i+1gℓi+1+fℓi=μai−1−λ+μai+λai+1+fℓi=aiλa−μa−1+fℓi≤0,for ℓ=0,1,2 and all i>K for some sufficiently large K if we choose a∈1,μ/λ=1,1/ρ. Due to K>0, the values Si will change. Now, we have
(59)0=−S0λ+S1μ,0=Si−1λ−Siλ+μ+Si+1μ,i=1,⋯,K−1,K+1,⋯,N,0=SN−1λ−SNλ+μ.The solution to this system is given by(60)Si=ρi−ki=0,⋯,Kρi−k−ρN+1−K1−ρN+1−k,i=K,⋯,N.We omit an explicit representation ofF¯ and F¯ (which is strictly larger than ψf here), but we remark that
(61)F¯−F¯=SNqN,N+1gN+1=ρN−K1−ρ1−ρN+1−K⋅λ⋅aN+1,and due to ρ<ρa<1, this difference tends to 0 as N⟶∞. Note that the speed of convergence of F¯−F¯ to 0 is still exponentially fast, but slower than for the optimal choices of gℓ.


## 6. Application to a Variant of theM/PH/1 Queue

Next, we consider an application of our method to a queueing model where we do not know the exact invariant distribution and have no chance but to use numerical methods. We consider a variant of theM/PH/1 queue where arriving customers decide whether to join the queue or not depending on the number of customers they find in the queue. Precisely, we assume that
(i)
Customers arrive according to the Poisson process with intensityλ(ii)
An arriving customer joins the queue with probabilityαn if there are n other customers in the system upon the arrival(iii)
The service time is phase-type distributed with parametersβ∈ℝ1×d and B∈ℝd×d, that is, the cumulative distribution function of a service time is given by t↦1−βeBt1 and its expectation by 1/μ≔β−B−11(iv)
There is one server.This queueing system can be modelled as the Markov chainY=Ytt≥0, where Yt=Nt,Ut is two-dimensional: Nt is the number of customers in the system at time t, and for Nt>0,Ut∈1,⋯,d is the current service phase. By setting the service phase to 1 for Nt=0, we obtain levels Ei with E0=0,1 and Ei=0,1,⋯,0,d for i≥1. The generator matrix is given by
(62)Q=−λα0λα0β−B1B−λα1Iλα1I−B1βB−λα3Iλα2I−B1βB−λα3Iλα3I⋱⋱⋱.Forαn=1, we obtain a level-independent QBD (the classical M/PH/1 queue) for which an explicit analytical representation of the solution exists (see, e.g., [23]). For nonconstant αn, we obtain an LDQBD, and we have to use numerical methods. In what follows, we assume that
(63)λ.limn→∞supαn<µ,which is equivalent to the existence of some λ~<μ with λαn≤λ~ for all n≥K0 with some K0∈N0. Without restriction, we assume K0≥2. Condition (63) guarantees that we have positive recurrence and that the stationary number of customers in the system has finite expectation as we will prove below. Hence, under condition (63), it makes sense to compute
(i)
the stationary probabilityπf1 that an arriving customer joins the queue where fi1=1·αi for all i(ii)
the stationary probabilityπf2 that at least M customers are in the system where fi2=1 for i≥M and fi2=0andi<M(iii)
the expectationπf3 of the stationary number of customers in the system where fi3=1·i for all i.In order to use our algorithmic method, we setfi0=1 additionally. For applying our method and for proving ψfℓ<∞ for ℓ=0,1,2,3 under condition (63), we check some f-modulated drift conditions. Set hi=Qi,i−1gi−1+Qii+Qi,i+1gi+1 for i≥1, and then we have hi=−B1βgi−1+Bgi+λαigi+1−gi≤−B1βgi−1+Bgi+λ~gi+1−gi for i≥K0 and all choices of g for which gi+1−gi≥0 (component-wise). For ℓ=0,1,2, we have fiℓ≤1, and we set giℓ=1⋅i+−B−11μ⋅1/μ−λ~ for i>1. Due to β1=1, we obtain
(64)hiℓ+fiℓ≤−B1β1⋅i−1−B1β−B−11⋅μ+B1⋅i−1⋅μ+1⋅λ~⋅1μ−λ~+fi=−B1⋅i−1−B1+B1⋅i−1μ−λ~⋅1μ−λ~+1=0,for i≥K0. From the special case fi0=1, we obtain positive recurrence from Theorem 2 in [8], and by the same result, we have ψfℓ<∞ for ℓ=1,2. For ℓ=3, we have fi3=1⋅i, and we set
(65)gi3=1⋅i2+−B−11⋅2μi⋅1μ−λ~,for i≥1 (note that −B−1 is a nonnegative matrix, and hence, we have gi+13≥gi3, implying
(66)hi3⋅μ−λ~≤−B1β1i+12−B1β−B−11⋅2μi+1+B1i2−1⋅2μi+1λ~i+12−i2+−B−11⋅2λ~μ=−B1i2−2i+1−B12i−2+B1i2−1⋅2μi+1⋅2λ~i+1λ~+−B−11⋅2λ~μ=−1⋅2μ−λ~i−B1+1λ~+−B−11⋅2λ~μ,for i≥K0. Hence,
(67)hi3+fi3≤−B1+1λ~+−B−11⋅2λ~μ.1μ−λ~−1⋅i.ChooseK≥K0 such that
(68)−B1+1λ~+−B−11⋅2λ~μ⋅1μ−λ~≤1⋅K+1.Then, we havehi3+fi3≤0 for i>K. The existence of such a value K proves ∑i=0∞πi1·i<∞, that is, the stationary number of customers has finite expectation.For applying our algorithmic method, we have to specifyK. In our numerical examples, we set λ=5, μ=1, αn=1/10+9/10n+1. This allows to choose λ~=3/4, and it is easy to derive that λαn≤λ~ for n≥K0=17. Furthermore, we consider the special case where the service times are Erlang-2-distributed, that is, β=1,0 and
(69)B=−2μ2μ0−2μ,implying
(70)−B1=02μ,−B−1=12μ12μ012μ,and
(71)−B−11⋅μ=112.Then, (68) is guaranteed if 2μ+3λ~/μ−λ~−K+1≤0⇔K≥2μ+3λ~/μ−λ~−1=μ+4λ~/μ−λ~. For μ=1 and λ~=3/4, this results in K≥8. Since we have to choose K≥K0, we set K=17. Note that giℓ simplifies to
(72)giℓ=1⋅i+−B−11μ⋅1μ−λ~=i+1i+124forℓ=0,1,2,gi3=1⋅i2+−B−11⋅2μi⋅1μ−λ~=i2+2i2+1⋅4.In Table1, some numerically computed results are listed. The parameters are chosen as specified above, and additionally, we set M=30 for πf2. The results clearly illustrate the convergence of both bounds to the same value. For obtaining precise results for πf2, the truncation level N has to be chosen a bit larger than for obtaining quite precise bounds on πf1 and πf3. This is not surprising since due to the fact that πf2 is the stationary probability for at least M=30 customers, the truncation level should be chosen significantly larger than 30. However, the results demonstrate that even for the computation of πf2, the truncation level 55 yields very precise results.Table 1
Numerical results for the considered variant of theM/PH/1 queue.
Nπf1πf2 ≥πf3Lower boundUpper boundLower boundUpper boundLower boundUpper bound200.1679710.32970000.1119187.8018712.3849250.1967640.2119700.0105579.312679.79899300.1996980.2009151.05844×10−58.76397×10−49.491819.53472350.1999270.2000026.09687×10−51.15438×10−49.507739.51063400.199940.1999446.8538×10−57.13126×10−59.50889.50895450.1999410.1999416.90867×10−56.92067×10−59.508859.50886500.1999410.1999416.91157×10−56.91203×10−59.508869.50886550.1999410.1999416.91170×10−56.91172×10−59.508869.50886


## 7. Application to a Retrial Queueing Model

As a final example, we consider theM/M/c/d−1 queue with retrials which can be seen as some kind of “prototype” for LDQBDs.
(i)
Customers arrive according to the Poisson process with intensityλ(ii)
The service times are independent and identically exponentially distributed with parameterμ(iii)
There arec servers(iv)
The system capacity isd−1≥c, that is, d−1−c customers can wait at the same time(v)
Customers who cannot enter the system due to lack of waiting capacity enter the “orbit” of retrials(vi)
Each customer in the orbit will retry to enter the system after a time which is exponentially distributed with parameterν and independent of all other random variables(vii)
Retrying customers who cannot enter the system due to lack of waiting capacity stay in the orbit of retrying customers.LetYt=Ot,Nt, where Ot is the number of customers in the orbit of retrials, and Nt is the number of customers in the queue (including service). Then, Y=Ytt≥0 is a continuous-time Markov chain with state space E∪i=0∞Ei, where Ei=i,0,⋯,i,d−1. Obviously, a transition with a positive rate will change the level at most by 1, and hence, we have a quasi-birth-death process with the following state transitions from state (i,u):
(i)
Arrivals occur with rateλ. For u<d−1, the arriving customer enters the queueing system, and the new state is i,u+1. For u=d−1, the arriving customer enters the orbit of retrials, and the new state is (i+1,d−1)(ii)
Foru>0, service completions occur with rate maxu,c⋅µ, and the new state is i,u−1(iii)
Fori>0 and u<d−1, successful retrials occur with rate i⋅ν, and the new state is i−1,u+1.In matrix notation, we have(73)Qi,i+1=0⋯00⋮⋱⋮⋮0⋯000⋯0λ,for i≥0,
(74)Qi,i−1=0vi0vi⋱⋱0vi0,for i≥1, and
(75)Qii=−λ−viλμ−λ−μ−viλ2μ−λ−2μ−viλ⋱⋱⋱cμ−λ−cμ−viλ⋱⋱⋱cμ−λ−cμ−viλcμ−λ−cμ,for i≥1.Note that retrial queues have been discussed intensively in the literature. We refer to the textbook [24] which gives an overview on retrial queueing models and computational methods for determining invariant distributions, stationary expectations, etc. Computational methods are very important in this context since even for d−1=c, there are only explicit analytical representations of the invariant distribution if c≤2 (see [24]).Of course, there are many interesting characteristics such as the mean number of customers in the orbit and the mean number of customers in the queueing system. Our purpose is to demonstrate that our method provides reliable results for interesting characteristics, but we do not want to discuss the retrial queueing model in details. Therefore, we focus on a single characteristic, the stationary probability for arriving customers to find the queuing system full and being forced to join the orbit. Ifπ=πi,ui,u∈E denotes the invariant distribution (if there is one), this probability is given by
(76)pb≔∑i=0∞πu,d−1.The invariant distribution exists if and only if we have positive recurrence. Below, we will restate a proof for the fact that this is the case ifλ<cμ. Note that λ<cμ is also a necessary requirement for stability/positive recurrence; we refer to [24] for more details.We demonstrate how to prove thatλ<cμ is sufficient for positive recurrence by means of drift criteria. The function g which we use for this purpose can also be used for finding the upper bounds on ψ1 and ψ1A for any A⊂E. In particular, by setting A=i,d−1:i∈ℕ0, we obtain bounds on pb≔ψ1A/ψ1.In any case (both for proving recurrence ifλ<cμ and for applying our truncation method), we want to find a function g such that we have hi+fi≔Qi,i−1gi−1+Qiigi+Qi,i+1gi+1+fi≤0 for i>K with some sufficiently large K, where 0≤fi≤1. We use the approach gi,u=γi+δu. For i≥1 and u<d−1, the entries of hi compute as
(77)hi,u=νigi−1,u+1−gi,u+λgi,u+1−gi,u+maxc,uµgi,u−1−gi,u=νiδ−γ+λδ−maxc,uµδ≤νiδ−γ+λδ=hi,0.If we chooseγ>δ, we surely have hi,u+fi,u≤0 for i>K with some sufficiently large K (due to fi,u≤1). Additionally, we have to take
(78)hi,d−1=λgi+1,d−1−gi,d−1+cµgi,d−2−gi,d−1=λγ−cµδ,into account. Hence, we have to guarantee λγ−cμδ≤−1. Due to λ<cμ, such a choice with γ>δ is possible. For example, we can choose δ=2/cμ−λ>0 and γ=1+cμ/λ/cμ−λ>2/cμ−λ, since then we have
(79)λγ=cμδ=1cμ−λλ+cμ−2λ=−1.For this choice, we havehi,d−1+1≤0 for all i≥1.As pointed out above, we have to guarantee that(80)hi,u+1≤hi,0+1=−viγ−δ+λδ+1=−viλ⋅λ1+cμ/λ−2cμ−λ+2λ+cμ−λcμ−λ=−viλ⋅i+cμ+λcμ−λ≤0,for i≥K+1. Hence, we choose
(81)K+1=λv⋅cμ+λcμ−λ.In total, we set(82)fi=10⋮⋮1011,gi=γiγiγi+δγi+δ⋮⋮γi+δd−1γi+δd−1,where γ=1/λ⋅cμ+λ/cμ−λ and δ=2/cμ−λ, and we choose K as above. Then, we can apply our algorithm which returns bounds on ψf0 and ψf1 where ψ is the invariant measure with ψK,0=1. By dividing lower and upper bounds appropriately, we find bounds on
(83)pb=ψf1ψf0.In Tables2 and 3, some numerically computed values for lower and upper bound are listed for different choices of the arrival rate λ and fixed values μ,ν,c,d. In both cases, we see that lower and upper bounds converge to the same limit (which coincides with the value computed in Table 3.8 in the textbook [24] by other means). In Table 2, we have low traffic due to λ/cμ=0.2. Then, the blocking probability is quite low, and—more important for the evaluation of our method—we obtain precise results for low truncation levels. In Table 3, we have more traffic due to λ/cμ=0.8. Hence, the blocking probability is much larger, and we need higher truncation levels to obtain precise results.Table 2
Numerical results forλ=1.0, μ=1.0, ν=0.5, c=d−1=5, and K=4.
Truncation levelNLower bound onpbUpper bound onpb51.85061×10−35.62599×10−362.79784×10−33.71132×10−373.12382×10−33.32205×10−383.20060×10−33.24191×10−393.21688×10−33.22540×10−3103.22025×10−33.22200×10−3113.22094×10−33.22130×10−3123.22109×10−33.22116×10−3133.22112×10−33.22113×10−3143.22112×10−33.22112×10−3Table 3
Numerical results forλ=4.0, μ=1.0, ν=0.5, c=d−1=5, and K=73.
Truncation levelNLower bound onpbUpper bound onpb752.67599×10−57218.46808.89475×10−52171.49907.39895×10−4261.0431005.79399×10−333.33531104.23642×10−24.559141200.2035450.9489031300.3849330.5017611400.4320690.4470221500.4385780.4403871600.4393740.4395901700.4394690.4394941800.4394800.4394831900.4394810.4394822000.4394820.439482The fact that we have to choose higher truncation levelsN in the situation of Table 3 is anything but surprising since due to more traffic, the Markov chain will assume states within higher levels much more often. However, it is important to remark that it is not clear how large N should be chosen. The suggested method can be iteratively applied with increased truncation levels until a prescribed accuracy is achieved. Unfortunately, when increasing N, we have to restart most of the calculations (B and W do not change, but R, Z, and Z¯ do). However, for each such calculation, the memory requirement is the same since it does not depend on N. In particular, in all situations considered in Tables 2 and 3, the memory requirement coincides.As pointed out above,pb is an important performance measure, but there are more interesting characteristics for this model. Note that with the same choice for g, we can find bounds on the stationary probability for any finite set of (finite or infinite) sets. In particular, if pu denotes the stationary probabilities for u customers in the queueing system (waiting or in service), we could determine bounds on pu for u=0,⋯,d−1 simultaneously with little additional effort. Other performance measures (e.g., moments of the numbers of retrying customers) will require other choices for g.We conclude this discussion of the basic retrial queueing model by remarking that there are many extensions of the retrial queueing model which improve the applicability. For example, we could consider Markovian arrival processes, phase-type distributed service times, impatient customers who leave the queue and enter the orbit, impatient customers who leave the orbit, queueing networks where declined external arrivals join the orbit of retrying customers, etc. For all such models, the suggested algorithm allows to compute important performance measures precisely and efficiently. Note that such retrial queues provide good models for many problems in telecommunications and computer networks, e.g., in the design of wireless networks. However, a deep discussion of a realistic (and hence complex) a model justifies a separate publication and is far beyond the scope of this methodological paper.


## 8. Evaluation of the Method

### 8.1. Thef-Modulated Drift Condition

The major requirement for finding the bounds in Theorem5 and for the resulting algorithmic procedure is the f-modulated drift condition (5) or (6). We give some comments on this condition.
(i)
As pointed out above,f-modulated drift conditions are very popular for proving convergence/finiteness of ψf or πf. For f≥0, we remark that the existence of an f-modulated drift condition is even equivalent to the convergence of ψf. Hence, from a purely theoretical point of view, the requirement that an f-modulated drift condition is satisfied is no restriction(ii)
Most often, finding a functiong which satisfies the f-modulated drift condition is the easiest way (or even the only way) to prove positive recurrence or to prove ψf<∞. If convergence of ψf is not guaranteed, it makes no sense to compute approximations to ψf (which are always finite) by numerical means. Hence, also from a practical point of view, finding a function g which satisfies the appropriate f-modulated drift condition does not require any additional effort in many situations(iii)
Every functiong which can be used for proving ψf<∞ can be used in Theorem 5. Note that this is not true for all truncation bounds which rely on f-modulated drift conditions, see Section 8.6.


### 8.2. The Structural Requirement (7)

We give some comments on condition (7). Remember that K is determined by the drift condition. If (7) is not met, we can always restructure the levels by setting E~0=E0∪⋯∪EK and E~i=EK+i for i≥1, and set K~=0. Then (7) is fulfilled. If we have the QBD structure for the levels Ei, this is preserved for the levels E~i. Consequently, we have d~0=d0+⋯+dK. This is unfortunate since dK (which becomes d~0 after the restructuring) should be quite small due to the operations involving matrices of dimension dK×dK.In total, from a purely mathematical point of view,K=0 can be guaranteed without restriction, and then the structural requirement (7) can be omitted. From a practical point of view, it is advantageous to keep d0,d1,⋯, and in particular, dK small, and in this context, we benefit from allowing K>0. Then, the price is the requirement (7).


### 8.3. Asymptotic Properties of the Bounds

In order to compare our bounds to competing methods, we consider the asymptotic behaviour of the difference between upper and lower bounds.In order to keep things simple and concise, we focus on the approximation ofπf1=ψf1/ψf0 (where f0=1) for positive recurrent continuous-time Markov chains, and we assume dK=1 (which is a clearly restrictive condition), implying φ¯=φ¯=ψK=1. The bounds on πf1 are given by F¯1/F¯0 and F¯1/F¯0, and for N⟶∞, the difference behaves asymptotically like
(84)F¯1F¯0−F¯1F¯0=F¯1−F¯1F¯0+F¯1F¯0⋅F¯0−F¯0F¯0=∑i=0NSi∑j=i+1∞Pijgj1+F¯1F¯0⋅∑i=0NSi+Pijgj0F¯0∼∑i=0Nπi∑j=i+1∞Pijgj1πfℓ⋅∑i=0Nπi∑j=i+1∞Pijgj0,where we have used that Si/F¯0=ψi/ψKF¯0 converges to πi. In particular, if gj1=gj0 (e.g., if f1 is an indicator function), we have the asymptotic expansion
(85)F¯1F¯0−F¯1F¯0∼1+πf1∑i=0Nπi∑j=1+1∞Pijgj0.These considerations clearly motivate to choosegjℓ as small as possible.


### 8.4. Augmented Matrices versus Nonaugmented Matrices

Theorem5 and the resulting method for computing bounds on ψf or πf rely on considering the northwest-corner truncation PN=Piji,j=0N which is substochastic. From PN, we obtain an approximation ψN to ψ, where ψN is uniquely determined by ψx0N=1 and by satisfying ψNPN=ψN up to the scalar equation corresponding to state x0.There are many authors who prefer to “repair” the row sums ofPN, that is, they deal with an augmented stochastic matrix P~N≥PN. Due to finiteness and stochasticity, it is possible to determine an invariant distribution π~N of P~N directly. Under some mild constraints, π~N converges to π (see, e.g., [3]). The fact that P~N admits an invariant distribution (in contrast to PN) can be interpreted as an advantage of methods relying on augmented matrices. There are more theoretical issues why considering augmented matrices is popular, and additionally, the augmented matrices often allow meaningful probabilistic interpretations. A more detailed list of advantages can be found in [3] (for the continuous-time setting).Despite these undeniable advantages, the approach relying on nonaugmented matrices was chosen in this paper. The simple reason is thatψN converges (pointwise) monotonically to ψ which makes lower bounds on ψf trivial. As we will see in the next lines, this fact is more important than all the advantages of augmentation procedures.


### 8.5. Comparison to Other A Posteriori Bounds on Truncation Errors

The most recent results for truncation bounds on stationary expectations can be found in [1–3] where the authors consider continuous-time Markov chains. They consider northwest-corner truncations Q=Qiji,j=0N and their augmented versions Q~N (with row sums 0). With our notation and the invariant distribution π~N of Q~N, Theorem 2.1 in [2] reads as
(86)∑x∈Eπx−π~xNfx≤1+πfinffx⋅π~NQ~−Q⋅g+bβϕ¯β.The authors require that(i)
the Markov chain is irreducible and positive recurrent(ii)
f≥1(iii)
Q~=Q~N=q~xyx,y∈E is constructed such that q~xyx,y∈E0∪⋯∪EN arises from truncation and augmentation and the other entries are arbitrary such that Q~ is conservative(iv)
e.g., all other entries are 0(v)
g satisfies the f-modulated drift condition ∑y∈Eqxygy≤−fxforx∈EK+1∪EK+2∪⋯(87)b=max∑y∈Eqxygy+fx:x∈E0∪⋯EK<∞(vi)
β>0 is arbitrary(vii)
ϕ¯β=supj∈Emini∈E0∪⋯∪EKϕijβ, where Φβ=ϕijβi,j∈E=∫0∞βe−βtPtdt with the transition function t↦Pt arising from Q.The right-hand side in (86) has some similarities to the asymptotic expansion (85). Asymptotically, the entries of π~NQ~−Qg behave like ∑i=0Nπi∑j=i+1∞Pijgj. Hence, (86) mainly differs by the factor inffx and the summand b/βϕ¯β1. Furthermore, note that the right-hand side of (86) provides a bound on the difference between πf and its approximation. Hence, “upperbound−lowerbound” is larger by factor 2.Finally, it is important to remark thatΦβ cannot be determined by means of numerical computations. Therefore, it is suggested in [2] to use I−1/β⋅QN−1 that increases monotonically in N with limit Φβ. Hence, the entries of I−1/β⋅QN−1 can be used to find computable bounds on ϕ¯β, and thus, on ∑x∈Eπx−π~xNfx.ReplacingΦβ in this manner has two effects:
(i)
There is no evidence that the bound in (86) is sharp. Even if it was, the bound which results from this replacement cannot be sharp anymore(ii)
Evaluating the computable bound which arises from (86) requires computing at least some rows of I−1/β⋅QN−1 which can be costly in terms of computational effort.Together, we see that the algorithm suggested in this paper will provide slightly better bounds which are significantly easier to compute at the same time. Finally,f≥1 in [1–3] is clearly restrictive. In particular, this requirement prevents us from using indicator functions in order to obtain stationary probabilities.


### 8.6. Comparison to A Priori Bounds on Truncation Errors

Theorem5 provides a posteriori truncation bounds (similar to (86) from [2]), since both the approximation to πf and the bounds are computed numerically (simultaneously).A priori bounds do not rely on numerical computations. Hence, before any computation scheme starts, a priori bounds would allow to determine the truncation levelN in such a way that the (relative) error of the approximation is smaller than a prescribed bound.An approach in this direction was introduced in [5] and generalized in [4] where a method was developed which determines a level N such that
(88)∑i=0Nψifi∑i=0∞ψifi≥1−∈,for some prescribed ∈>0 and some scalar function f>0. The method works as follows:
(i)
Choosesome functiong and determine (continuous-time setting) hx=∑y∈Eqxygyforallx∈E.(ii)
SetC0=x:hx>0 and c=maxx∈C0hx/fx.(iii)
Setγ=c⋅1−∈/∈ and C=x∈E:hx>−γfx.(iv)
IfC is finite, we have ∑x∈Cψxfx/∑x∈Eψxfx≥1−∈. In particular, we can find N such that E0∪⋯∪EN⊃C and obtain (88).(v)
IfC is not finite, we have to change our choice of g.At a first glance, this method provides a priori error bounds, but there are some problems. In the original papers [4, 5], it was already pointed out that it is not easy to find appropriate functions g, and that very often, the bounds are quite pessimistic. In order to illustrate these effects, we consider the M/M/1 queue from Section 5 and set f=1. Whereas Theorem 5 can be applied with gi=b⋅i for an appropriate constant b, this choice will fail here: We obtain h0=bλ and hi=bλ−μ for i≥1, that is, hi is constant for i≥1. In particular, if ∈>0 is very small (and hence, γ is large), we have C=i∈E:hi>−γfi=ℕ0. Note that the choice of b is unimportant since it cancels out.The next intuitive choice isgi=i2, resulting in h0=λ and hi=2λ−µi+λ+μ. If λ<μ<2λ, we have c=h1=3λ−μ, γ=3λ−μ⋅1−∈/∈, and
(89)C=i:2λ−μi+λ+μ>−3λ−μ⋅1−∈∈=0,1,⋯,3λ−μ2μ−λ⋅1∈.Forλ/μ=2/3, we obtain C=0,1,⋯,N where N=3/2⋅1/∈. Hence, the method guarantees that ∑i=0Nπi≥1−∈ whereas the exact value is ∑i=0Nπi=1−2/3N. For ∈⟶0 and N=N∈, the latter term converges to 1 much faster than 1−∈. For example, for ∈=0.05, we have N=30, and hence, 1−2/3N≈1−5.215⋅10−6.Summarizing this brief numerical example, the most intuitive choice forg does not work, and the next choice provides very pessimistic estimates. It might be possible to find functions g which provide sharp estimates, but it does not seem likely that such functions can be “guessed.”There is an even more serious problem when interpreting the results from [4, 5] as a priori estimates for the truncation error: The numerator in (88) refers to the exact invariant measure ψ which cannot be computed by numerical means (up to some trivial exceptions like the M/M/1 queue). For honest a priori bounds, we should be interested in guaranteeing
(90)∑i=0NψiNfi∑i=0∞ψifi≥1−∈,with a numerically calculated approximation ψN to ψ. Unfortunately, the methods in [4, 5] cannot be modified in this direction. One might argue that this slight modification is negligible in view of the pessimistic bounds which are obtained. However, if we guessed an optimal function g by accident, we would not have a bound on ∑i=0NψiNfi/∑i=0∞ψifi anymore. In total, the results from [4, 5] may be exploited to find theoretical estimates on stationary probabilities or expectations. With respect to numerical means, they can only be used for rough preliminary considerations (how large should the truncation level N be chosen) before applying the method suggested in this paper. Even that might be hard due to the increased difficulty to find appropriate functions g (in comparison to finding appropriate functions g for applying Theorem 5).


## 8.1. Thef-Modulated Drift Condition

The major requirement for finding the bounds in Theorem5 and for the resulting algorithmic procedure is the f-modulated drift condition (5) or (6). We give some comments on this condition.
(i)
As pointed out above,f-modulated drift conditions are very popular for proving convergence/finiteness of ψf or πf. For f≥0, we remark that the existence of an f-modulated drift condition is even equivalent to the convergence of ψf. Hence, from a purely theoretical point of view, the requirement that an f-modulated drift condition is satisfied is no restriction(ii)
Most often, finding a functiong which satisfies the f-modulated drift condition is the easiest way (or even the only way) to prove positive recurrence or to prove ψf<∞. If convergence of ψf is not guaranteed, it makes no sense to compute approximations to ψf (which are always finite) by numerical means. Hence, also from a practical point of view, finding a function g which satisfies the appropriate f-modulated drift condition does not require any additional effort in many situations(iii)
Every functiong which can be used for proving ψf<∞ can be used in Theorem 5. Note that this is not true for all truncation bounds which rely on f-modulated drift conditions, see Section 8.6.


## 8.2. The Structural Requirement (7)

We give some comments on condition (7). Remember that K is determined by the drift condition. If (7) is not met, we can always restructure the levels by setting E~0=E0∪⋯∪EK and E~i=EK+i for i≥1, and set K~=0. Then (7) is fulfilled. If we have the QBD structure for the levels Ei, this is preserved for the levels E~i. Consequently, we have d~0=d0+⋯+dK. This is unfortunate since dK (which becomes d~0 after the restructuring) should be quite small due to the operations involving matrices of dimension dK×dK.In total, from a purely mathematical point of view,K=0 can be guaranteed without restriction, and then the structural requirement (7) can be omitted. From a practical point of view, it is advantageous to keep d0,d1,⋯, and in particular, dK small, and in this context, we benefit from allowing K>0. Then, the price is the requirement (7).


## 8.3. Asymptotic Properties of the Bounds

In order to compare our bounds to competing methods, we consider the asymptotic behaviour of the difference between upper and lower bounds.In order to keep things simple and concise, we focus on the approximation ofπf1=ψf1/ψf0 (where f0=1) for positive recurrent continuous-time Markov chains, and we assume dK=1 (which is a clearly restrictive condition), implying φ¯=φ¯=ψK=1. The bounds on πf1 are given by F¯1/F¯0 and F¯1/F¯0, and for N⟶∞, the difference behaves asymptotically like
(84)F¯1F¯0−F¯1F¯0=F¯1−F¯1F¯0+F¯1F¯0⋅F¯0−F¯0F¯0=∑i=0NSi∑j=i+1∞Pijgj1+F¯1F¯0⋅∑i=0NSi+Pijgj0F¯0∼∑i=0Nπi∑j=i+1∞Pijgj1πfℓ⋅∑i=0Nπi∑j=i+1∞Pijgj0,where we have used that Si/F¯0=ψi/ψKF¯0 converges to πi. In particular, if gj1=gj0 (e.g., if f1 is an indicator function), we have the asymptotic expansion
(85)F¯1F¯0−F¯1F¯0∼1+πf1∑i=0Nπi∑j=1+1∞Pijgj0.These considerations clearly motivate to choosegjℓ as small as possible.


## 8.4. Augmented Matrices versus Nonaugmented Matrices

Theorem5 and the resulting method for computing bounds on ψf or πf rely on considering the northwest-corner truncation PN=Piji,j=0N which is substochastic. From PN, we obtain an approximation ψN to ψ, where ψN is uniquely determined by ψx0N=1 and by satisfying ψNPN=ψN up to the scalar equation corresponding to state x0.There are many authors who prefer to “repair” the row sums ofPN, that is, they deal with an augmented stochastic matrix P~N≥PN. Due to finiteness and stochasticity, it is possible to determine an invariant distribution π~N of P~N directly. Under some mild constraints, π~N converges to π (see, e.g., [3]). The fact that P~N admits an invariant distribution (in contrast to PN) can be interpreted as an advantage of methods relying on augmented matrices. There are more theoretical issues why considering augmented matrices is popular, and additionally, the augmented matrices often allow meaningful probabilistic interpretations. A more detailed list of advantages can be found in [3] (for the continuous-time setting).Despite these undeniable advantages, the approach relying on nonaugmented matrices was chosen in this paper. The simple reason is thatψN converges (pointwise) monotonically to ψ which makes lower bounds on ψf trivial. As we will see in the next lines, this fact is more important than all the advantages of augmentation procedures.


## 8.5. Comparison to Other A Posteriori Bounds on Truncation Errors

The most recent results for truncation bounds on stationary expectations can be found in [1–3] where the authors consider continuous-time Markov chains. They consider northwest-corner truncations Q=Qiji,j=0N and their augmented versions Q~N (with row sums 0). With our notation and the invariant distribution π~N of Q~N, Theorem 2.1 in [2] reads as
(86)∑x∈Eπx−π~xNfx≤1+πfinffx⋅π~NQ~−Q⋅g+bβϕ¯β.The authors require that(i)
the Markov chain is irreducible and positive recurrent(ii)
f≥1(iii)
Q~=Q~N=q~xyx,y∈E is constructed such that q~xyx,y∈E0∪⋯∪EN arises from truncation and augmentation and the other entries are arbitrary such that Q~ is conservative(iv)
e.g., all other entries are 0(v)
g satisfies the f-modulated drift condition ∑y∈Eqxygy≤−fxforx∈EK+1∪EK+2∪⋯(87)b=max∑y∈Eqxygy+fx:x∈E0∪⋯EK<∞(vi)
β>0 is arbitrary(vii)
ϕ¯β=supj∈Emini∈E0∪⋯∪EKϕijβ, where Φβ=ϕijβi,j∈E=∫0∞βe−βtPtdt with the transition function t↦Pt arising from Q.The right-hand side in (86) has some similarities to the asymptotic expansion (85). Asymptotically, the entries of π~NQ~−Qg behave like ∑i=0Nπi∑j=i+1∞Pijgj. Hence, (86) mainly differs by the factor inffx and the summand b/βϕ¯β1. Furthermore, note that the right-hand side of (86) provides a bound on the difference between πf and its approximation. Hence, “upperbound−lowerbound” is larger by factor 2.Finally, it is important to remark thatΦβ cannot be determined by means of numerical computations. Therefore, it is suggested in [2] to use I−1/β⋅QN−1 that increases monotonically in N with limit Φβ. Hence, the entries of I−1/β⋅QN−1 can be used to find computable bounds on ϕ¯β, and thus, on ∑x∈Eπx−π~xNfx.ReplacingΦβ in this manner has two effects:
(i)
There is no evidence that the bound in (86) is sharp. Even if it was, the bound which results from this replacement cannot be sharp anymore(ii)
Evaluating the computable bound which arises from (86) requires computing at least some rows of I−1/β⋅QN−1 which can be costly in terms of computational effort.Together, we see that the algorithm suggested in this paper will provide slightly better bounds which are significantly easier to compute at the same time. Finally,f≥1 in [1–3] is clearly restrictive. In particular, this requirement prevents us from using indicator functions in order to obtain stationary probabilities.


## 8.6. Comparison to A Priori Bounds on Truncation Errors

Theorem5 provides a posteriori truncation bounds (similar to (86) from [2]), since both the approximation to πf and the bounds are computed numerically (simultaneously).A priori bounds do not rely on numerical computations. Hence, before any computation scheme starts, a priori bounds would allow to determine the truncation levelN in such a way that the (relative) error of the approximation is smaller than a prescribed bound.An approach in this direction was introduced in [5] and generalized in [4] where a method was developed which determines a level N such that
(88)∑i=0Nψifi∑i=0∞ψifi≥1−∈,for some prescribed ∈>0 and some scalar function f>0. The method works as follows:
(i)
Choosesome functiong and determine (continuous-time setting) hx=∑y∈Eqxygyforallx∈E.(ii)
SetC0=x:hx>0 and c=maxx∈C0hx/fx.(iii)
Setγ=c⋅1−∈/∈ and C=x∈E:hx>−γfx.(iv)
IfC is finite, we have ∑x∈Cψxfx/∑x∈Eψxfx≥1−∈. In particular, we can find N such that E0∪⋯∪EN⊃C and obtain (88).(v)
IfC is not finite, we have to change our choice of g.At a first glance, this method provides a priori error bounds, but there are some problems. In the original papers [4, 5], it was already pointed out that it is not easy to find appropriate functions g, and that very often, the bounds are quite pessimistic. In order to illustrate these effects, we consider the M/M/1 queue from Section 5 and set f=1. Whereas Theorem 5 can be applied with gi=b⋅i for an appropriate constant b, this choice will fail here: We obtain h0=bλ and hi=bλ−μ for i≥1, that is, hi is constant for i≥1. In particular, if ∈>0 is very small (and hence, γ is large), we have C=i∈E:hi>−γfi=ℕ0. Note that the choice of b is unimportant since it cancels out.The next intuitive choice isgi=i2, resulting in h0=λ and hi=2λ−µi+λ+μ. If λ<μ<2λ, we have c=h1=3λ−μ, γ=3λ−μ⋅1−∈/∈, and
(89)C=i:2λ−μi+λ+μ>−3λ−μ⋅1−∈∈=0,1,⋯,3λ−μ2μ−λ⋅1∈.Forλ/μ=2/3, we obtain C=0,1,⋯,N where N=3/2⋅1/∈. Hence, the method guarantees that ∑i=0Nπi≥1−∈ whereas the exact value is ∑i=0Nπi=1−2/3N. For ∈⟶0 and N=N∈, the latter term converges to 1 much faster than 1−∈. For example, for ∈=0.05, we have N=30, and hence, 1−2/3N≈1−5.215⋅10−6.Summarizing this brief numerical example, the most intuitive choice forg does not work, and the next choice provides very pessimistic estimates. It might be possible to find functions g which provide sharp estimates, but it does not seem likely that such functions can be “guessed.”There is an even more serious problem when interpreting the results from [4, 5] as a priori estimates for the truncation error: The numerator in (88) refers to the exact invariant measure ψ which cannot be computed by numerical means (up to some trivial exceptions like the M/M/1 queue). For honest a priori bounds, we should be interested in guaranteeing
(90)∑i=0NψiNfi∑i=0∞ψifi≥1−∈,with a numerically calculated approximation ψN to ψ. Unfortunately, the methods in [4, 5] cannot be modified in this direction. One might argue that this slight modification is negligible in view of the pessimistic bounds which are obtained. However, if we guessed an optimal function g by accident, we would not have a bound on ∑i=0NψiNfi/∑i=0∞ψifi anymore. In total, the results from [4, 5] may be exploited to find theoretical estimates on stationary probabilities or expectations. With respect to numerical means, they can only be used for rough preliminary considerations (how large should the truncation level N be chosen) before applying the method suggested in this paper. Even that might be hard due to the increased difficulty to find appropriate functions g (in comparison to finding appropriate functions g for applying Theorem 5).


## 9. Conclusion and Further Research

In this paper, we have developed an exact and efficient method for numerically computing lower and upper bounds onψf where ψ is an invariant measure for an irreducible recurrent discrete-time or continuous-time QBD and f is a multidimensional nonnegative cost/reward function. By combining the lower and upper bounds appropriately, we obtain lower and upper bounds on πf in the case of positive recurrence, where π is the unique invariant distribution and f is an arbitrary cost/reward function. The term πf can be interpreted as the stationary expectation of the Markov chain measured by f, or as the long-run average of cost/rewards measured by f. If f=1A is chosen to be an indicator function, we obtain stationary probabilities. The computation of πf2 in Section 6 illustrates that even very small probabilities can be computed with high relative precision.The algorithm was developed as an extension of previous matrix-analytic methods which also dealt algorithmically with level-dependent QBDs, but for which no error bounds were known up to now. Due to the fact that we compute lower and upper error bounds, we obtain automatically an approximation toπf and (a posteriori) error bounds. In our examples, we have seen that the difference between lower and upper bounds converges to 0 very fast, that is, we obtain very precise results. For optimal choices of gℓ, we have observed that the upper bound coincides with the true value of ψf.Obvious directions of future publications concern detailed applications of our algorithm to specific Markov models, e.g., we could consider enhanced retrial queuing models with Markovian arrival processes, phase-type distributed service times, impatient customers, etc., but there is a huge variety of other practical processes which are modelled by infinite LDQBDs. If we focus on the analysis of a specific model, we should spend effort into finding very good functionsg which may result in obtaining very precise results for small truncation levels.In this context, it may be worth to find not only an upper truncation level but also a lower one. Then, the numerical computation would focus on a very small region of the state space, and the method would become even more efficient.Note that the main results of this paper are not restricted to QBDs. Theorem5 only requires that no transitions from Ei to Ej occur with positive probability/rate if i>K>j. Hence, these results also apply if P or Q, respectively, has upper block-Hessenberg structure (also referred to as the block-M/G/1 structure). We have focused on QBDs here since this specialized structure allows a very efficient computation of lower and upper bounds. Future research could deal with developing efficient algorithms for more general structures.We remark that some of the key results also transfer to nondiscrete state spaces: LetE=E0∪E1∪⋯, and let the levels Ei be petite in the sense of [10] (this property is satisfied by compact sets for many Markov chains). Then Pij becomes a (sub-Markovian) kernel, and Si does so, too. The multiplication SjPji is defined straightforward by SjPjix,dz=∫ySjx,dyPjiy,dz. We still have ψf=ψKF, where ψK solves ψKlimN→∞I−T=0, and
(91)F=limN→∞fK+∑i=0Ni≠KSifi,where Sifiℓx=∫ySis,dyfiℓy. Under the structural assumptions (no transitions from Ei to Ej for i>K>j), the results on the bounds on F remain unchanged. However, for finding S0,⋯,SN by means of numerical calculations (and for finding a lower and an upper bound on ψK), it seems that some sort of discretization is inevitable.We conclude with a remark concerning the relationship to nonprobabilistic topics: As mentioned above, forK=0, the basic algorithm for computing approximations to the subvectors ψi of the invariant measure ψ can be deduced from general considerations concerning (matrix-valued) continued fractions. Hence, in some sense, our computation of bounds on ψf is related to the speed-of-convergence statements for matrix-valued continued fractions. Note that the probabilistic interpretation of continued fractions (arising in nonprobabilistic contexts) and their convergents has led to new convergence criteria for continued fractions and their generalizations (see [25, 26]), and with further research, considerations similar to those in this paper might provide speed-of-convergence results.



---

*Source: 2678374-2020-08-07.xml*