Componentwise accurate fluid queue computations using doubling algorithms

Abstract

Markov-modulated fluid queues are popular stochastic processes frequently used for modelling real-life applications. An important performance measure to evaluate in these applications is their steady-state behaviour, which is determined by the stationary density. Computing it requires solving a (nonsymmetric) M-matrix algebraic Riccati equation, and indeed computing the stationary density is the most important application of this class of equations. Xue et al. (Numer Math 120:671–700, 2012) provided a componentwise first-order perturbation analysis of this equation, proving that the solution can be computed to high relative accuracy even in the smallest entries, and suggested several algorithms for computing it. An important step in all proposed algorithms is using so-called triplet representations, which are special representations for M-matrices that allow for a high-accuracy variant of Gaussian elimination, the GTH-like algorithm. However, triplet representations for all the M-matrices needed in the algorithm were not found explicitly. This can lead to an accuracy loss that prevents the algorithms from converging in the componentwise sense. In this paper, we focus on the structured doubling algorithm, the most efficient among the proposed methods in Xue et al., and build upon their results, providing (i) explicit and cancellation-free expressions for the needed triplet representations, allowing the algorithm to be performed in a really cancellation-free fashion; (ii) an algorithm to evaluate the final part of the computation to obtain the stationary density; and (iii) a componentwise error analysis for the resulting algorithm, the first explicit one for this class of algorithms. We also present numerical results to illustrate the accuracy advantage of our method over standard (normwise-accurate) algorithms.

Appendices

Appendix: Doubling and censoring

In this section, we derive some results connecting the doubling map $\mathcal {F}$ to the concept of censoring (or stochastic complementation) in Markov chains, which we shall need in the next section when we prove the numerical accuracy of Algorithm 2, our componentwise accurate SDA. In particular, Theorem 6.1 and Lemma 6.2 are standard results on censoring that we shall need.

Theorem 6.1

[21, Theorems 2.1 and 2.2] Consider $\mathcal {S} = \{1, \ldots , m\} = \mathcal {S}_1 \cup \mathcal {S}_2$, where $\mathcal {S}_1 \cap \mathcal {S}_2 = \emptyset$, and let $P$ be an $m \times m$ stochastic matrix partitioned as

$$\begin{aligned} P=\begin{bmatrix} P_{11}&\quad \! P_{12} \\ P_{21}&\quad \! P_{22} \end{bmatrix}\!, \end{aligned}$$

(29)

where $P_{\ell m}$ contains the transition probabilities from $\mathcal {S}_{\ell }$ to $\mathcal {S}_m$. Suppose that $I-P_{11}$ is invertible. Then,

• (i) $P':=P_{22}+P_{21}(I-P_{11})^{-1}P_{12}$ is stochastic.

• (ii) If a row vector $\varvec{\xi }=\begin{bmatrix} \varvec{\xi }_1&\quad \! \varvec{\xi }_2 \end{bmatrix}\ge \varvec{0}$ is such that $\varvec{\xi } P=\varvec{\xi }$, then $\varvec{\xi }_2 P'=\varvec{\xi }_2$.

The matrix $P'$ is said to be the matrix obtained by censoring the set $\mathcal {S}_1$ or the submatrix $P_{11}$. Censoring can alternatively be described in fully linear algebraic terms: the matrix $I-P'$ is the Schur complement of the submatrix $I-P_{11}$ in $I-P$.

Lemma 6.2

(Quotient property, [10]) Let $P$ be an $m \times m$ stochastic matrix, and $\mathcal {T}$ and $\mathcal {U}$ be two disjoint subsets of $\mathcal {S} = \{1,2,\dots ,m\}$, with $\mathcal {S} \backslash (\mathcal {T} \cup \mathcal {U}) \ne \emptyset$. Define $P'$ to be the matrix obtained from $P$ by censoring $\mathcal {T}$, and $P''$ the matrix obtained from $P'$ by censoring $\mathcal {U}$.

Then, $P''$ is also the matrix obtained from $P$ by censoring $\mathcal {T}\cup \mathcal {U}$.

It was proved in [9] that the doubling map $\mathcal {F}$ (13) corresponds to a step in another, similar iteration, Cyclic Reduction, on the matrices

$$\begin{aligned} A_+ :=\begin{bmatrix} E&\quad \! 0\\ 0&\quad \! 0 \end{bmatrix}\!, \quad A_= :=\begin{bmatrix} 0&\quad \! G\\ H&\quad \! 0 \end{bmatrix}\!, \quad A_-:=\begin{bmatrix} 0&\quad \! 0\\ 0&\quad \! F \end{bmatrix}\!. \end{aligned}$$

(30)

Moreover, it was also established [5, 8] that Cyclic Reduction corresponds to Schur complementation on a suitable tridiagonal block-Toeplitz matrix; thus, it should not be surprising that there is a direct interpretation of doubling as Schur complementation/censoring, which is our next result. This result is slightly different from the existing ones, since we are interested with block-circulant matrices rather than block-tridiagonal ones here; nevertheless, the proof is similar (see, for instance, [8, Section 7.3]).

Let $Z_m$ be the $m \times m$ circulant generator matrix, defined by

$$\begin{aligned} (Z_m)_{ij} = {\left\{ \begin{array}{ll} 1 &{} j-i \equiv 1 \mod m,\\ 0 &{} \text {otherwise}, \end{array}\right. } \end{aligned}$$

and denote by $\otimes$ the Kronecker product.

Theorem 6.3

Let $P$ be a $n\times n$ stochastic matrix, and define $E$, $F$, $G$, and $H$ as in (11), and $A_+$, $A_=$, and $A_-$ as in (30). Then, for $k \ge 0$ the matrix $\mathcal {F}^k(P)$ is obtained by censoring the top $(2^k-1)n\times (2^k-1)n$ block of the $2^kn\times 2^kn$ matrix

$$\begin{aligned} P^{(2^k)}&:= Z_{2^k}^{-1}\otimes A_- + I_{2^k} \otimes A_= + Z_{2^k} \otimes A_+ \\&= \begin{bmatrix} A_=&\quad \! A_+&\quad \!&\quad \!&\quad \! A_-\\ A_-&\quad \! A_=&\quad \! A_+&\quad \!&\quad \! \\&\quad \! A_-&\quad \! A_=&\quad \! \ddots&\quad \!&\quad \!\\&\quad \!&\quad \!\ddots&\quad \! \ddots&\quad \! A_+\\ A_+&\quad \!&\quad \!&\quad \!A_-&\quad \! A_= \end{bmatrix}\!. \nonumber \end{aligned}$$

(31)

Here and in the following, for clarity we omit blocks of zeros from some large block matrices such as the one in (31).

Proof

We prove the theorem by induction. For $k=1$,

Then, using the identity

$$\begin{aligned} \left( I-\begin{bmatrix} 0&\quad \! G\\ H&\quad \! 0 \end{bmatrix}\right) ^{-1} = \begin{bmatrix} I&\quad \! G\\ H&\quad \! I \end{bmatrix}\begin{bmatrix} (I-GH)^{-1}&\quad \! 0 \\ 0&\quad \! (I-HG)^{-1} \end{bmatrix} \end{aligned}$$

we can verify the result directly.

Now, we assume that the result holds for a certain $k$ and prove it for $k+1$. Let $\widehat{P}:=\mathcal {F}(P)$. We first prove that the matrix obtained by censoring out the components corresponding to the blocks with odd index numbers from $P^{(2^{k+1})}$ is exactly $\widehat{P}^{(2^k)}$, that is, the matrix with the same structure (31) built starting from $\widehat{P}$. The result then follows from Lemma 6.2.

We first reorder the blocks in $P^{(2^{k+1})}$ to put the odd-numbered ones on the top, obtaining

Next, we censor the top part, $I \otimes A_{=}$, to obtain

$$\begin{aligned}&I \otimes A_= + (I \otimes A_- + Z \otimes A_+)(I\otimes (I - A_=))^{-1}(Z^{-1} \otimes A_- + I\otimes A_+)\\&\quad =Z^{-1} \otimes A_-(I-A_=)^{-1}A_- + I \otimes (A_= + A_-(I-A_=)^{-1}A_+ \\&\quad \quad +\,A_+(I-A_=)^{-1}A_-) + Z \otimes A_+(I-A_=)^{-1}A_+\\&\quad = Z^{-1} \otimes \begin{bmatrix} 0&\quad \! 0\\ 0&\quad \! \widehat{F} \end{bmatrix} + I \otimes \begin{bmatrix} 0&\quad \! \widehat{G}\\ \widehat{H}&\quad \! 0 \end{bmatrix} + Z \otimes \begin{bmatrix} \widehat{E}&\quad \! 0\\ 0&\quad \! 0 \end{bmatrix}\!, \end{aligned}$$

where $\widehat{E},\widehat{F},\widehat{G},$ and $\widehat{H}$ are the blocks of $\widehat{P}$ as in (12). The last expression is precisely $\widehat{P}^{(2^k)}$. By the inductive hypothesis, the matrix obtained by censoring the top $(2^k-1)n \times (2^k-1)n$ of this matrix is $\mathcal {F}^k(\widehat{P})=\mathcal {F}^{k+1}(P)$, which completes the proof. $\square$

Remark 6.4

One can obtain a different, probabilistic proof for this result relying on a physical interpretation of doubling algorithms in terms of the underlying fluid queue. We are currently working on an explanation for these algorithms from a probabilistic point of view, giving some insight that complements the linear-algebraic point of view presented in this paper.

Similarly, we can express the initial values of the doubling algorithms as the censoring of a suitable stochastic matrix.

Theorem 6.5

Let $(T,C)$ be the transition and rate matrices of a fluid queue, and let $\gamma \le \left( \max _{i} Q_{ii}\right) ^{-1}$. Then, $P_0$ as in (15) is obtained by censoring the first $n$ states out of

$$\begin{aligned} S=\begin{bmatrix} I-\gamma Q&\quad \! \gamma R\\ I&\quad \! 0 \end{bmatrix}\!. \end{aligned}$$

(32)

Proof

Clearly, $0+I(I-(I-\gamma Q))^{-1}\gamma R=Q^{-1}R$, so it suffices to prove that $S$ is stochastic. For the nonnegativity, we need only to verify that the diagonal entries of $R$ are nonnegative, which holds because $\alpha \le \alpha _{\mathrm {opt}},\beta \le \beta _{\mathrm {opt}}$, and that the diagonal entries of $I-\gamma Q$ are nonnegative, which holds because of the choice of $\gamma$. Moreover, substituting the definitions of $Q$ and $R$ from (15), and using $T\mathbf {1}=\varvec{0}$, one sees that $S\mathbf {1}=\mathbf {1}$. $\square$

Corollary 6.6

Let $\varvec{\xi }\ge 0$ be such that $\varvec{\xi }T=\varvec{0}$. Then, for $k \ge 0$, $P_k$ is stochastic and $\varvec{\xi }|C|P_k= \varvec{\xi }|C|$.

Proof

By their definitions and the fact that $\varvec{\xi }T=\varvec{0}$,

$$\begin{aligned} \varvec{\xi } Q =\varvec{\xi }|C|=\varvec{\xi } R, \end{aligned}$$

and consequently that $\begin{bmatrix} \varvec{\xi }&\quad \! \gamma \varvec{\xi }|C| \end{bmatrix}S=\begin{bmatrix} \varvec{\xi }&\quad \! \gamma \varvec{\xi }|C| \end{bmatrix}$. Hence, by applying Theorem 6.1 to $P_0$ as a censoring of $S$, the corollary holds for $k=0$.

For $k > 0$, one can verify directly that $(\mathbf {1}^\top \otimes \varvec{\xi })P^{(2^k)}=(\mathbf {1}^\top \otimes \varvec{\xi })$, so the result follows again from applying Theorem 6.1 to $P_k$ as a censoring of $P^{(2^k)}$. $\square$

Appendix: Componentwise stability of SDA

In this section, we prove the results on componentwise stability of SDA, thus giving the proof of Theorem 4.4. We begin with a perturbation bound that tells us how the iterates are affected by a small change in the initial values.

1.1 Componentwise perturbation bounds

We focus on first-order results, assuming $\varepsilon$ to be a small parameter and hence ignoring all terms containing $\varepsilon ^2$ and higher powers of $\varepsilon$. For brevity, we use the notation $A \mathrel {\dot{\le }}B$ to denote $A \le [1+O(\varepsilon ^2)] B$. Moreover, we focus at first on the case $n_+\approx n_-$; hence, we replace $n_+$ and $n_-$ liberally with $n_{\max }:=\max (n_+,n_-)$ in our bound. The results of a more accurate analysis that keeps track of the differences between $n_+$ and $n_-$ are reported later in Sect. 7.3.

We state here two useful lemmas. The first is a small variant of [1, Lemma 2.2], while the second is a simple consequence of the triangle inequality.

Lemma 7.1

Let $({{\mathrm{offdiag}}}(A),\mathbf {1}, \varvec{w})$ and $({{\mathrm{offdiag}}}(\tilde{A}),\mathbf {1},$ $\tilde{\varvec{w}})$ be triplet representations for two $m\times m$ M-matrices $A$ and $\tilde{A}$ such that

$$\begin{aligned} |{{\mathrm{offdiag}}}(A)-{{\mathrm{offdiag}}}(\tilde{A})|\le \varepsilon |{{\mathrm{offdiag}}}(A)| \quad \text{ and } \quad |\varvec{w}-\tilde{\varvec{w}}| \le \varepsilon \varvec{w}. \end{aligned}$$

Then,

$$\begin{aligned} |A^{-1}-\tilde{A}^{-1}| \mathrel {\dot{\le }}(2m-1)\varepsilon A^{-1}\!. \end{aligned}$$

(33)

Lemma 7.2

Let $A,B,\tilde{A},$ and $\tilde{B}$ be nonnegative matrices of suitable sizes such that $|A-\tilde{A}|\mathrel {\dot{\le }}a\varepsilon A$ and $|B-\tilde{B}|\mathrel {\dot{\le }}b\varepsilon B$. Then,

• (i) $|A+B-(\tilde{A}+\tilde{B})| \mathrel {\dot{\le }}\max (a,b) \varepsilon (A+B)$.

• (ii) $|AB-\tilde{A}\tilde{B}| \mathrel {\dot{\le }}(a+b)\varepsilon AB$.

The following lemma provides a perturbation bound for censoring.

Lemma 7.3

Consider two stochastic matrices

$$\begin{aligned} P=\begin{bmatrix} P_{11}&\quad \! P_{12}\\ P_{21}&\quad \! P_{22} \end{bmatrix}\quad \text { and }\quad \tilde{P}=\begin{bmatrix} \tilde{P}_{11}&\quad \! \tilde{P}_{12} \\ \tilde{P}_{21}&\quad \! \tilde{P}_{22} \end{bmatrix}\!, \end{aligned}$$

with $P_{11},\tilde{P}_{11}\in \mathbb {R}^{m\times m}$, $|{{\mathrm{offdiag}}}(P-\tilde{P})|\le \varepsilon {{\mathrm{offdiag}}}(P),$ and suppose that $I-P_{11}$ is invertible. Then,

$$\begin{aligned} |\tilde{P}_{21}(I-\tilde{P}_{11})^{-1}\tilde{P}_{12}-P_{21}(I-P_{11})^{-1}P_{12}| \mathrel {\dot{\le }}(2m+1)\varepsilon P_{21}(I-P_{11})^{-1}P_{12}. \end{aligned}$$

If, in addition, $|\mathrm {diag}(P_{22})-\mathrm {diag}(\tilde{P}_{22})|\le (2m+1)\varepsilon P_{22}$, then

$$\begin{aligned}&|\tilde{P}_{22}+\tilde{P}_{21}(I-\tilde{P}_{11})^{-1}\tilde{P}_{12}-P_{22}-P_{21}(I-P_{11})^{-1}P_{12}|\\&\quad \mathrel {\dot{\le }}(2m+1)\varepsilon (P_{22}+P_{21}(I-P_{11})^{-1}P_{12}). \end{aligned}$$

Proof

We have

$$\begin{aligned} |(I-P_{11})\mathbf {1}-(I-\tilde{P}_{11})\mathbf {1}|=|P_{12}\mathbf {1}-\tilde{P}_{12}\mathbf {1}|\le \varepsilon P_{12}\mathbf {1}. \end{aligned}$$

Thus, $I-P_{11}$ and $I-\tilde{P}_{11}$ satisfy the hypotheses of Lemma 7.1 and hence $|(I-\tilde{P}_{11})^{-1}-(I-P_{11})^{-1}|\le (2m-1)\varepsilon (I-P_{11})^{-1}$. The results now follow by propagating errors with the help of Lemma 7.2. $\square$

We start from a componentwise perturbation bound for $P_0$.

Lemma 7.4

Let $(T,C)$ and $(\tilde{T},\tilde{C})$ be two pairs of the transition matrix and the rate matrix of two fluid queues such that

$$\begin{aligned} |{{\mathrm{offdiag}}}(T)-{{\mathrm{offdiag}}}(\tilde{T})|\le \varepsilon {{\mathrm{offdiag}}}(T) \quad \text{ and } \quad |C-\tilde{C}|\le \varepsilon |C|. \end{aligned}$$

Consider $P_0=Q^{-1}R$ and its analogous $\tilde{P}_0=\tilde{Q}^{-1}\tilde{R}$ computed using the same formulas (15) but starting from $\tilde{T}$ and $\tilde{C}$, and with the same values of $\alpha$ and $\beta$, chosen such that $\alpha \le \alpha _{\mathrm {opt}},\beta \le \beta _{\mathrm {opt}}$. Then,

$$\begin{aligned} |P_0-\tilde{P}_0| \mathrel {\dot{\le }}(2n+1) \varepsilon P_0. \end{aligned}$$

Proof

We have

$$\begin{aligned} |\tilde{R}_{ii}-R_{ii}|&= |\tilde{C}_+ + \beta \tilde{T}_{ii} - C_+ - \beta T_{ii}| \\&\le |\tilde{C}_+ -C_+| + \beta |\tilde{T_{ii}} - T_{ii}| \\&\le \varepsilon C_+ + \beta \sum _{j\ne i} |\tilde{T}_{ij}-T_{ij}| \\&\le \varepsilon C_+ + \beta \varepsilon \sum _{j\ne i} T_{ij} \\&= \varepsilon R_{ii}. \end{aligned}$$

Hence, for $S$ from (32) and its analog $\tilde{S}$ (where we take the same value of $\gamma$ for both), we have the bound $|{{\mathrm{offdiag}}}(\tilde{S})-{{\mathrm{offdiag}}}(S)|\le \varepsilon {{\mathrm{offdiag}}}(S)$. By Lemma 7.3, we obtain the required result. $\square$

Next, we derive a perturbation bound for the result of $k$ steps of SDA.

Lemma 7.5

Let $P$ and $\tilde{P}$ be two $n\times n$ stochastic matrices such that $|\tilde{P}-P|\le \varepsilon P$, and $\mathcal {F}$ be the doubling map (13). Then,

$$\begin{aligned} |\mathcal {F}^k(\tilde{P})-\mathcal {F}^k(P)|\mathrel {\dot{\le }}n 2^k \varepsilon \mathcal {F}^k(P). \end{aligned}$$

Proof

Consider the matrices $P^{(2^k)}$ and its analog $\tilde{P}^{(2^k)}$ defined in the same way but starting from $\tilde{P}$. They are stochastic matrices, and they satisfy the hypotheses of Lemma 7.3. We apply the second part of this lemma, which completes the proof. $\square$

By applying the first part of Lemma 7.3 rather than the second, we get a corresponding bound for

$$\begin{aligned} J_k:=G_{k}-G_{k-1}=E_{k-1}(I-G_{k-1}H_{k-1})^{-1}G_{k-1}F_{k-1} > 0. \end{aligned}$$

Corollary 7.6

With the hypotheses of Lemma 7.5, we have

$$\begin{aligned} |\tilde{J}_{k}- J_{k}| \mathrel {\dot{\le }}n 2^k\varepsilon J_k. \end{aligned}$$

The above bound gives us immediately $|\tilde{G}_k-G_k|\mathrel {\dot{\le }}n2^k\varepsilon G_k$ for the doubling iterates starting from two nearby sets of matrices. This bound, however, is quite weak: in particular, it would allow for $\lim _{k\rightarrow \infty } |\tilde{G}_k-G_k|=\infty$, while one would expect this quantity to converge, since $G_k\rightarrow \Psi$ and $\tilde{G}_k\rightarrow \tilde{\Psi }$.

Here, we overestimate the error: essentially, we compute a new term of the sum $G_k = G_0 + J_1 + J_2 + \dots + J_k$ at each step, and at the same time make the following estimate:

$$\begin{aligned} |\tilde{G}_k-G_k|&\le |\tilde{G}_0-G_0| + |\tilde{J}_1-J_1| + |\tilde{J}_1-J_1| + |\tilde{J}_2-J_2| + \dots + |\tilde{J}_k-J_k| \\&\mathrel {\dot{\le }} \varepsilon G_0 + n\varepsilon J_1 + 2n \varepsilon J_2 + \dots + 2^kn \varepsilon J_k \\&\le 2^kn(G_0+J_1+J_2+\dots + J_k) \end{aligned}$$

Indeed, the last line is how errors in sums are bounded using Lemma 7.2. However, the factors $J_k$ decay exponentially like $O(\delta ^{2^k})$, and this compensates for the growing powers of two. The following theorem formalizes this intuition to improve this bound.

Theorem 7.7

With the hypotheses of Lemma 7.5, we have

$$\begin{aligned} |\tilde{G}_k-G_k| \mathrel {\dot{\le }}\left( 1 + \frac{4\widehat{K}_0n}{1-\delta }\right) \varepsilon \Psi . \end{aligned}$$

for a suitable constant $\widehat{K}_0$.

Proof

We know from Theorem 3.1 that $0\le \Psi -G_k\le K_0 \delta ^{2^k}$, and that $\Psi -G_k\le \Psi$. Hence we can write, for a suitable constant $\widehat{K}_0$ (different in general from $K_0$), $\Psi -G_k \le \widehat{K}_0 \delta ^{2^k}\Psi$.

$$\begin{aligned} |\tilde{G}_k-G_k|&\le |\tilde{G}_0-G_0|+ \sum _{h=1}^k |\tilde{J}_{h} - J_{h}| \\&\mathrel {\dot{\le }} \varepsilon G_0 + n \varepsilon \sum _{h=1}^k 2^h J_h \\&\le \varepsilon G_0 + n \varepsilon \sum _{h=1}^k 2^h (\Psi -G_{h-1}) \\&\le \varepsilon G_0 + n \varepsilon \sum _{h=1}^k 2^h \widehat{K}_0 \delta ^{2^{h-1}}\Psi \\&\le \varepsilon G_0 + 2\widehat{K}_0n \varepsilon \sum _{h=1}^k 2^{h-1} \delta ^{2^{h-1}}\Psi \\&\le \varepsilon G_0 + 4\widehat{K}_0n \varepsilon \sum _{j=1}^{2^k} \delta ^{j} \Psi \le \varepsilon \Psi + \frac{4\widehat{K}_0n}{1-\delta } \varepsilon \Psi \end{aligned}$$

$\square$

Clearly, an analogous result for $H_k$ holds.

1.2 Numerical stability

We now establish a more complex result that keeps track also of the numerical errors in an implementation of the algorithm in floating-point arithmetic. We denote by ${\mathrm {u}}$ the machine precision, and ignore second-order terms in ${\mathrm {u}}$ by hiding them again in the syntax $A \mathrel {\dot{\le }}B$.

We shall rely on the two following lemmas. The first one is a version of Lemma 7.2 that keeps track of floating point arithmetic errors as well.

Lemma 7.8

Let $A$ and $B$ be nonnegative matrices of suitable sizes, and let $\tilde{A}$ and $\tilde{B}$ be matrices of floating point numbers such that $|A-\tilde{A}|\mathrel {\dot{\le }}a {\mathrm {u}} A$, $|B-\tilde{B}|\mathrel {\dot{\le }}b {\mathrm {u}} B$.

• 1. Let $\tilde{S}$ be the result of the matrix addition $\tilde{A}+\tilde{B}$ performed in floating point arithmetic. Then,

  $$\begin{aligned} |\tilde{S}-(A+B)|\mathrel {\dot{\le }}(\max (a,b)+1){\mathrm {u}}(A+B). \end{aligned}$$

• 2. Let $\tilde{P}$ be the result of the matrix product $\tilde{A}\tilde{B}$ performed in floating point arithmetic. Then,

  $$\begin{aligned} |\tilde{P}-(AB)|\mathrel {\dot{\le }}(a+b+m){\mathrm {u}}(AB), \end{aligned}$$

  where $m$ is the number of columns of $A$ (and rows of $B$).

The second one is a slight variation of Theorem 4.1 to convert it to the same form as the bounds in Lemma 7.8.

Lemma 7.9

Let $A, \tilde{A} \in \mathbb {R}^{m\times m}$ be two M-matrices with triplet representations $({{\mathrm{offdiag}}}(A),\mathbf {1},\varvec{w})$ and $({{\mathrm{offdiag}}}(\tilde{A}),\mathbf {1},\tilde{\varvec{w}})$, and $B, \tilde{B} \in \mathbb {R}^{m\times k}$ be nonnegative, such that

$$\begin{aligned} |{{\mathrm{offdiag}}}(A)-{{\mathrm{offdiag}}}(\tilde{A})|&\mathrel {\dot{\le }} a{\mathrm {u}}|{{\mathrm{offdiag}}}(A)|, \\ |B-\tilde{B}|&\mathrel {\dot{\le }} b{\mathrm {u}}B, \\ |\varvec{w}-\tilde{\varvec{w}}|&\mathrel {\dot{\le }} a{\mathrm {u}}\varvec{w}, \end{aligned}$$

and all of ${{\mathrm{offdiag}}}(\tilde{A}), \tilde{B},$ and $\tilde{\varvec{w}}$ are floating-point numbers.

Let $\tilde{X}$ be the result of the operation $\tilde{A}^{-1}\tilde{B}$ performed in floating-point arithmetic with the GTH-like algorithm. Then,

$$\begin{aligned} |\tilde{X}-A^{-1}B| \mathrel {\dot{\le }}((2m-1)a+b+\psi (m)){\mathrm {u}} A^{-1}B. \end{aligned}$$

Proof

By Lemma 7.1, $|A^{-1}-\tilde{A}^{-1}|\mathrel {\dot{\le }}(2m-1){\mathrm {u}}A^{-1}$. Using this result and Lemma 7.2, we get

$$\begin{aligned} |\tilde{A}^{-1}\tilde{B}-A^{-1}B| \mathrel {\dot{\le }}((2m-1)a+b){\mathrm {u}} A^{-1}B. \end{aligned}$$

By Theorem 4.1,

$$\begin{aligned} |\tilde{X}-\tilde{A}^{-1}\tilde{B}| \mathrel {\dot{\le }}\psi (m){\mathrm {u}}\tilde{A}^{-1}\tilde{B}\mathrel {\dot{\le }}\psi (m){\mathrm {u}}A^{-1}B. \end{aligned}$$

The result follows using the triangle inequality. $\square$

Our first result towards an accuracy proof concerns the initial values.

Proposition 7.10

Let Algorithm 2 be performed in floating-point arithmetic with $\alpha \le 2\alpha _\mathrm {opt}$ and $\beta \le 2\beta _\mathrm {opt}$. Then, the computed approximation $\tilde{P}_0$ to $P_0$ satisfies the bound

$$\begin{aligned} |\tilde{P}_0-P_0| \mathrel {\dot{\le }}K_1{\mathrm {u}}P_0, \end{aligned}$$

with $K_1:=\psi (n)+5n+2$ and $\psi (n)$ defined as in Theorem 4.1.

Proof

We assume without loss of generality that $\alpha$ and $\beta$ are exact machine numbers (otherwise, we simply replace them with their machine representation). The machine representations of $T$ and $c_i$ satisfy $|\tilde{T}-T|\le \mathrm {u}T$ and $|\tilde{c}_i-c_i|\le \mathrm {u}c_i$.

Let $\tilde{Q}$ and $\tilde{R}$ be the computed versions of $Q$ and $R$ in floating-point arithmetic. They are computed by taking the machine representation of the corresponding entry of $T$ and multiplying it by $\alpha$ or $\beta$, hence

$$\begin{aligned} |{{\mathrm{offdiag}}}(\tilde{Q})-{{\mathrm{offdiag}}}(Q)| \mathrel {\dot{\le }}2\mathrm {u}{{\mathrm{offdiag}}}(Q) \end{aligned}$$

and

$$\begin{aligned} |{{\mathrm{offdiag}}}(\tilde{R})-{{\mathrm{offdiag}}}(R)| \mathrel {\dot{\le }}2\mathrm {u}{{\mathrm{offdiag}}}(Q). \end{aligned}$$

The quantities $d_i:=\beta \sum _{j\ne i} T_{ij}$, for $i\in \mathcal {S}_+$, are computed (with $n-1$ additions and one multiplication) up to an error $|\tilde{d}_i-d_i|\mathrel {\dot{\le }}(n+1) \mathrm {u}d_i$, as we get by applying repeatedly Lemma 7.8. Moreover, due to our choice of $\beta$, $c_i \ge 2d_i$, which implies $d_i\le c_i-d_i$ and $c_i\le 2(c_i-d_i)$. Hence we have

$$\begin{aligned} |\mathrm{fl }(\tilde{c}_i-\tilde{d}_i) - (c_i-d_i)|&= |(\tilde{c}_i-\tilde{d}_i)(1+\varepsilon )-(c_i-d_i)|\\&\mathrel {\dot{\le }} |\tilde{c}_i-c_i|+|\tilde{d}_i-d_i| + \mathrm {u}|\tilde{c}_i-\tilde{d}_i|\\&\mathrel {\dot{\le }} \mathrm {u}c_i + (n+1)\mathrm {u}d_i + \mathrm {u}(c_i-d_i)\\&\le (2+n+1+1)\mathrm {u}(c_i-d_i)=(n+4)\mathrm {u}(c_i-d_i). \end{aligned}$$

A similar bound applies for $i\in \mathcal {S}_-$. Putting everything together, we have $|\tilde{R}-R|\le (n+4)R$, $|{{\mathrm{offdiag}}}(\tilde{Q})-{{\mathrm{offdiag}}}(Q)|\mathrel {\dot{\le }}2\mathrm {u}{{\mathrm{offdiag}}}(Q)$, and $|\tilde{C}\mathbf {1}-C\mathbf {1}|\le \mathrm {u}C\mathbf {1}$. Hence, by Lemma 7.9

$$\begin{aligned} |\tilde{P}_0-P_0|\le (2(2n-1)+(n+4)+\psi (n))\mathrm {u}P_0. \end{aligned}$$

$\square$

We are now ready to assess the numerical error a single SDA step.

Lemma 7.11

Let $\tilde{\mathcal {F}}$ be the map that computes approximately the result of one iteration of SDA implementing Equations (12) in machine arithmetic. For any $n\times n$ stochastic matrix $P$, we have

$$\begin{aligned} |\tilde{\mathcal {F}}(P)-\mathcal {F}(P)| \mathrel {\dot{\le }}K_2\mathrm {u}\mathcal {F}(P), \end{aligned}$$

with $K_2:=\psi (n_{\max }) + 4n_{\max }^2+2n_{\max }-1$.

Proof

Let $E,F,G,H$ be the blocks of $P$, $\widehat{E},\widehat{F},\widehat{G},\widehat{H}$ those of $\mathcal {F}(P)$, and $\tilde{E},\tilde{F},\tilde{G},\tilde{H}$ those of $\tilde{\mathcal {F}}(P)$. For simplicity, with a slight abuse of notation, for any matrix $X\ge 0$ appearing in the algorithm we denote here by $c(X)$ its computed value in machine arithmetic along the SDA step computation. Applying Lemma 7.8 repeatedly, we have $|c(GH)-GH| \le n_-{\mathrm {u}} \tilde{G}\tilde{H}$ and

$$\begin{aligned} |c(F\mathbf {1}) - F\mathbf {1}|&\mathrel {\dot{\le }} n_-{\mathrm {u}} F\mathbf {1},\\ |c(GF\mathbf {1}) - GF\mathbf {1}|&\mathrel {\dot{\le }} 2n_- {\mathrm {u}} GF\mathbf {1},\\ |c(E\mathbf {1}) - E\mathbf {1}|&\mathrel {\dot{\le }} n_+{\mathrm {u}} E\mathbf {1},\\ |c(GF\mathbf {1}+E\mathbf {1}) - GF\mathbf {1}- E\mathbf {1}|&\mathrel {\dot{\le }} (2n_{\max }+1){\mathrm {u}} (GF\mathbf {1}- E\mathbf {1}). \end{aligned}$$

Hence

$$\begin{aligned}&|c((I-GH)^{-1}G)-(I-GH)^{-1}G| \\&\quad \mathrel {\dot{\le }}((2n_{\max }-1)(2n_{\max }+1) + \psi (n_{\max })) \mathrm {u}(I-GH)^{-1}G, \end{aligned}$$

and carrying on the error propagation in the same fashion leads to

$$\begin{aligned} |\tilde{G}-\widehat{G}|&= |c(G+E(I-GH)^{-1}GF)- G - E(I-GH)^{-1}GF| \\&\mathrel {\dot{\le }}(\psi (n_{\max }) + 4n_{\max }^2+2n_{\max }-1) \mathrm {u}\widehat{G}. \end{aligned}$$

The error computations for the other three iterates $\tilde{E}$, $\tilde{F}$, $\tilde{H}$ are analogous. $\square$

We are now ready to tackle several doubling steps together.

Theorem 7.12

Let $\tilde{P}_k=\tilde{\mathcal {F}}^k(\tilde{P}_0)$ be the approximation to $P_k=\mathcal {F}^k(P_0)$ computed by Algorithm 2 performed in floating-point arithmetic with machine precision ${\mathrm {u}}$, using values of the parameters $\alpha ,\beta$ such that $\alpha \le 2\alpha _\mathrm {opt}$ and $\beta \le 2\beta _\mathrm {opt}$. Then,

$$\begin{aligned} |\tilde{P}_k-P_k| \le n 2^k (K_1+K_2)\mathrm {u}P_{k}. \end{aligned}$$

Proof

We prove the result by induction on $k$.

The following manipulation is nothing but a formal version of the statement that when considering first-order error bounds, we can simply add up the local errors in every step of the algorithm. We start from the telescopic sum

$$\begin{aligned} |\tilde{\mathcal {F}}^k(\tilde{P}_0)-\mathcal {F}^k(\tilde{P}_0)|&\le \sum _{h=1}^k |\mathcal {F}^{h-1}\tilde{\mathcal {F}}^{k-h+1}(\tilde{P}_0)-\mathcal {F}^{h}\tilde{\mathcal {F}}^{k-h}(\tilde{P}_0)| \nonumber \\&= \sum _{h=1}^k |\mathcal {F}^{h-1}\tilde{\mathcal {F}}(\tilde{P}_{k-h})-\mathcal {F}^{h-1}\mathcal {F}(\tilde{P}_{k-h})|. \end{aligned}$$

(34)

By Lemma 7.11,

$$\begin{aligned} |\tilde{\mathcal {F}}(\tilde{P}_{k-h}) - \mathcal {F}(\tilde{P}_{k-h})| \mathrel {\dot{\le }}K_2\mathrm {u}\mathcal {F}(\tilde{P}_{k-h}); \end{aligned}$$

hence we can apply Lemma 7.5 with $\varepsilon =K_2\mathrm {u}$, to obtain for each $h \ge 1$

$$\begin{aligned} |\mathcal {F}^{h-1}\tilde{\mathcal {F}}(\tilde{P}_{k-h})-\mathcal {F}^{h-1}\mathcal {F}(\tilde{P}_{k-h})|&\mathrel {\dot{\le }} n2^{h-1} K_2 \mathrm {u}\mathcal {F}^h(\tilde{P}_{k-h})\\&\mathrel {\dot{\le }} n2^{h-1} K_2 \mathrm {u}\mathcal {F}^k(P_0). \end{aligned}$$

In the last line, we have substituted $P_{k-h}$ for $\tilde{P}_{k-h}$, which is possible because the two expressions differ only for a term of magnitude $O(\mathrm {u})$ by inductive hypothesis.

Insert this inequality into (34) to obtain

$$\begin{aligned} |\tilde{\mathcal {F}}^k(\tilde{P}_0)-\mathcal {F}^k(\tilde{P}_0)|&\mathrel {\dot{\le }} \sum _{h=1}^k n2^{h-1} K_2\mathrm {u}\mathcal {F}^k(P_0)\nonumber \\&< n 2^k K_2\mathrm {u}\mathcal {F}^k(P_0). \end{aligned}$$

(35)

Moreover, using again Lemma 7.5 we get

$$\begin{aligned} |\mathcal {F}^k(\tilde{P}_0)-\mathcal {F}^k(P_0)| \le n2^k K_1\mathrm {u}\mathcal {F}^k(P_0). \end{aligned}$$

(36)

Combining (35) and (36) we get the required result. $\square$

Again, this error bound is exponential in $k$, which is quite unwelcome. We can improve it using the same trick as in Theorem 7.7. Since computing $J_k$ is just removing the final addition from the algorithm, the previous bound can be easily modified to yield $|\tilde{J}_k-J_k| \le n2^k(K_1+K_2)\mathrm {u}J_k$. Now we proceed as in the proof of Theorem 7.7 to prove Theorem 4.4.

Proof (of Theorem 4.4) The matrix $\tilde{G}_k$ is the result of the addition (in machine arithmetic) of the matrices $\tilde{G}_0, \tilde{J}_1, \tilde{J}_2, \dots , \tilde{J}_k$, performed in this order at different steps along the algorithm. For this addition, we have

$$\begin{aligned} |\tilde{G}_k-G_k|&\mathrel {\dot{\le }} |(\tilde{G}_0 + \sum _{h=1}^k \tilde{J}_h)(1+k\mathrm {u}) - (G_0 + \sum _{h=1}^k J_h)| \\&\mathrel {\dot{\le }} k\mathrm {u}G_k + |\tilde{G}_0-G_0|+\sum _{h=1}^k |\tilde{J}_h-J_h| \\&\mathrel {\dot{\le }} k\mathrm {u}G_k + K_1\mathrm {u}G_0+\sum _{h=1}^k n2^h (K_1+K_2)\mathrm {u}J_h \\&\le k\mathrm {u}G_k + K_1\mathrm {u}G_0+(K_1+K_2)n\sum _{h=1}^k 2^h \mathrm {u}\tilde{K}_0\delta ^{2^h}\Psi \\&\le k\mathrm {u}\Psi + K_1\mathrm {u}\Psi + \frac{2(K_1+K_2)n}{1-\delta }\tilde{K}_0\Psi . \end{aligned}$$

$\square$

1.3 The unbalanced case

In the previous bounds we have replaced $n_+$ and $n_-$ with $n_{\max }$ for simplicity; hence, these bounds are overestimated whenever $n_-\ne n_+$. It turns out that there is a different strategy that we can adopt in this case, which leads to a slightly smaller value of $K_2$. First of all, we recall from [29] that the following alternative formulas equivalent to (12) can also be used for the doubling step

$$\begin{aligned} E_{k + 1}&:= E_k(I+G_k(I - H_k G_k)^{-1}H_k)E_k, \end{aligned}$$

(37a)

$$\begin{aligned} F_{k + 1}&:= F_k(I+H_k(I - G_k H_k)^{-1}G_k)F_k, \end{aligned}$$

(37b)

$$\begin{aligned} G_{k + 1}&:= G_k + E_k G_k(I - H_kG_k)^{-1}F_k, \end{aligned}$$

(37c)

$$\begin{aligned} H_{k + 1}&:= H_k + F_kH_k(I - G_kH_k)^{-1}E_k. \end{aligned}$$

(37d)

If $n_+=n_-$, the formulas (12) are the ones with the lower computational cost; however, if $n_+ < n_-$, one can invert the $n_+\times n_+$ matrix $I-G_k H_k$ rather than the larger $n_-\times n_-$ matrix $I-H_kG_k$, that is, use (12a), (37b), (12c), (37d), and vice versa if $n_- < n_+$.

By using the appropriate formulas, one can derive a slightly tighter version of Lemma 7.11. We report only the final result here, since the proof is basically unchanged.

Theorem 7.13

If we choose for each iterate the formula with the smallest matrix to invert among (12) and (37), then in Lemma 7.11 (as well as in Theorems 7.12 and 4.4) we may replace the factor $K_2$ with

$$\begin{aligned} \tilde{K}_2&:= \psi (n_{\min })+(2n_{\max }+1)(2n_{\min }-1) +2n_{\max }+n_{\min } \\&= \psi (n_{\min })+4n_{\max }n_{\min } + 3n_{\min }. \end{aligned}$$