Highly accurate doubling algorithms for M-matrix algebraic Riccati equations

Abstract

The doubling algorithms are very efficient iterative methods for computing the unique minimal nonnegative solution to an M-matrix algebraic Riccati equation (MARE). They are globally and quadratically convergent, except for MARE in the critical case where convergence is linear with the linear rate 1 / 2. However, the initialization phase and the doubling iteration kernel of any doubling algorithm involve inverting nonsingular M-matrices. In particular, for MARE in the critical case, the M-matrices in the doubling iteration kernel, although nonsingular, move towards singular M-matrices at convergence. These inversions are causes of concerns on entrywise relative accuracy of the eventually computed minimal nonnegative solution. Fortunately, a nonsingular M-matrix can be inverted by the GTH-like algorithm of Alfa et al. (Math Comp 71:217–236, 2002) to almost full entrywise relative accuracy, provided a triplet representation of the matrix is known. Recently, Nguyen and Poloni (Numer Math 130(4):763–792, 2015) discovered a way to construct triplet representations in a cancellation-free manner for all involved M-matrices in the doubling iteration kernel, for a special class of MAREs arising from Markov-modulated fluid queues. In this paper, we extend Nguyen’s and Poloni’s work to all MAREs by also devising a way to construct the triplet representations cancellation-free. Our construction, however, is not a straightforward extension of theirs. It is made possible by an introduction of novel recursively computable auxiliary nonnegative vectors. As the second contribution, we propose an entrywise relative residual for an approximate solution. The residual has an appealing feature of being able to reveal the entrywise relative accuracies of all entries, large and small, of the approximation. This is in marked contrast to the usual legacy normalized residual which reflects relative accuracies of large entries well but not so much those of very tiny entries. Numerical examples are presented to demonstrate and confirm our claims.

Appendix: Sparsity of \(E_k,F_k,X_k,Y_k\)

For convenience, we introduce a partial ordering on nonnegative matrices with respect to their entrywise nonzero patterns. For matrices $P\ge 0$, $Q\ge 0$ of the same size, we say that Q majorizes P in the entrywise nonzero pattern, written as $P \,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Q$, if $Q_{(i,j)}=0$ implies $P_{(i,j)}=0$, and write $P\,{\mathop {=}\limits ^{\scriptstyle 0}}\,Q$ if $P \,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Q$ and $Q \,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,P$. Evidently, $0\le P\le Q$ implies $P \,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Q$, but not the other way around.

Lemma A.1

• 1. Given $0\le P_i\le Q_i$ for $i=1,2$, all of the same size, we have

  $$\begin{aligned} P_1+P_2\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Q_1+Q_2, \quad P_1P_2\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Q_1Q_2. \end{aligned}$$

• 2. Given $0\le P\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Q$, we have $P+Q\,{\mathop {=}\limits ^{\scriptstyle 0}}\,Q$.

Lemma A.2

Let P, Q be nonsingular M-matrices, and split P and Q as

$$\begin{aligned} P&=D_P-N_P, \quad D_P={{\mathrm{diag}}}(P), \end{aligned}$$

(A.1a)

$$\begin{aligned} Q&=D_Q-N_Q, \quad D_Q={{\mathrm{diag}}}(Q). \end{aligned}$$

(A.1b)

If $N_P\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,N_Q$, then $P^{-1}\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Q^{-1}$. In particular, if $N_P\,{\mathop {=}\limits ^{\scriptstyle 0}}\,N_Q$, then $P^{-1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,Q^{-1}$.

Proof

Since P, Q are nonsingular M-matrices, $P^{-1}$ and $Q^{-1}$ admit the following series respresentations

$$\begin{aligned} P^{-1}=D_P^{-1}\sum _{k=0}^{\infty } \left( N_PD_P^{-1}\right) ^k,\quad Q^{-1}=D_Q^{-1}\sum _{k=0}^{\infty } \left( N_QD_Q^{-1}\right) ^k. \end{aligned}$$

(A.2)

If $N_P\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,N_Q$, then $(N_PD_P^{-1})^k\,{\mathop {=}\limits ^{\scriptstyle 0}}\,N_P^k \,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,N_Q^k \,{\mathop {=}\limits ^{\scriptstyle 0}}\,(N_QD_Q^{-1})^k$ for all $k\ge 0$ and thus $P^{-1}\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Q^{-1}$. In the case when $N_P\,{\mathop {=}\limits ^{\scriptstyle 0}}\,N_Q$, we also have $N_Q\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,N_P$ and thus $Q^{-1}\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,P^{-1}$ at the same time. $\square$

Lemma A.3

For a nonsingular M-matrix P, $P^{-1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,P^{-k}$ for $k\ge 1$, and $(\alpha I-P^{-1})^{-1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,P^{-1}$ for $\alpha >\rho (P^{-1})$.

Proof

In this proof and that of Lemma A.4 later, the series $\sum _{i=0}^{\infty } N_P^i$ is used only symbolically for its entrywise nonzero pattern since the series itself may not even converge. Again we have (A.2) which yields that $P^{-1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,\sum _{i=0}^{\infty } N_P^i$, and that

$$\begin{aligned} P^{-k}=(P^{-1})^k \,{\mathop {=}\limits ^{\scriptstyle 0}}\,\left( \sum _{i=0}^{\infty } N_P^i\right) ^k \,{\mathop {=}\limits ^{\scriptstyle 0}}\,\sum _{i=0}^{\infty } N_P^i \,{\mathop {=}\limits ^{\scriptstyle 0}}\,P^{-1}. \end{aligned}$$

For $\alpha >\rho (P^{-1})$, we have $(\alpha I-P^{-1})^{-1} =\alpha ^{-1}\sum _{i=0}^{\infty }(\alpha P)^{-i} \,{\mathop {=}\limits ^{\scriptstyle 0}}\,P^{-1},$ as expected. $\square$

Lemma A.4

Let P be a nonsingular M-matrix and $Q\ge 0$, and split P, Q as

$$\begin{aligned} P&=D_P-N_P, \quad D_P={{\mathrm{diag}}}(P), \end{aligned}$$

(A.3a)

$$\begin{aligned} Q&=D_Q+\widetilde{N}_Q, \quad D_Q={{\mathrm{diag}}}(Q). \end{aligned}$$

(A.3b)

If all of the diagonal entries of $D_Q$ are positive and $\widetilde{N}_Q\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,N_P$, then

$$\begin{aligned} P^{-1}Q\,{\mathop {=}\limits ^{\scriptstyle 0}}\,P^{-1}. \end{aligned}$$

Proof

Since $P^{-1}\ge 0$ and $Q\ge 0$, $P^{-1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,P^{-1}D_Q\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,P^{-1}D_Q+P^{-1}\widetilde{N}_Q=P^{-1}Q$. On the other hand,

$$\begin{aligned} P^{-1}Q\,{\mathop {=}\limits ^{\scriptstyle 0}}\,P^{-1}+P^{-1}\widetilde{N}_Q\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,P^{-1}+P^{-1}N_P \,{\mathop {=}\limits ^{\scriptstyle 0}}\,\sum _{i=0}^{\infty } N_P^i \,{\mathop {=}\limits ^{\scriptstyle 0}}\,P^{-1}, \end{aligned}$$

as was to be shown. $\square$

Theorem A.1

Let $E_0,F_0,X_0,Y_0$ be as in (3.6b), and let $E_k,F_k,X_k,Y_k$ be produced by the doubling iteration (3.1). Then for $k\ge 0$,

$$\begin{aligned} \begin{bmatrix} E_{k+1}&\quad Y_{k+1}\\ X_{k+1}&\quad F_{k+1} \end{bmatrix} \,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,\begin{bmatrix} E_k&\quad Y_k \\ X_k&\quad F_k \end{bmatrix}, \end{aligned}$$

(A.4)

$X_{k+1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,X_k\,{\mathop {=}\limits ^{\scriptstyle 0}}\,\Phi$, and $Y_{k+1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,Y_k\,{\mathop {=}\limits ^{\scriptstyle 0}}\,\Psi$.

Proof

It is clear that $X_k\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,X_{k+1}$ and $Y_k\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,Y_{k+1}$ because of (3.1c) and (3.1d). So if (A.4) is proven true, then we will immediately have $X_{k+1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,X_k\,{\mathop {=}\limits ^{\scriptstyle 0}}\,\Phi$ and $Y_{k+1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,Y_k\,{\mathop {=}\limits ^{\scriptstyle 0}}\,\Psi$.

It remains to prove (A.4). Because of (3.6), it suffices to prove (A.4) with a hat over every symbol. Set $Q=\begin{bmatrix} N_B&D\\ C&N_A \end{bmatrix}.$ We have

$$\begin{aligned} \begin{bmatrix} B+{\hat{\alpha }} I_m&\quad -D\\ -C&\quad A+{\hat{\beta }} I_n \end{bmatrix}&=\begin{bmatrix} D_B+{\hat{\alpha }} I_m&\quad 0 \\ 0&\quad D_A+{\hat{\beta }} I_n \end{bmatrix} -Q, \\ \begin{bmatrix} {\hat{\beta }} I_m-B&\quad D\\ C&\quad {\hat{\alpha }} I_n-A \end{bmatrix}&=\begin{bmatrix} {\hat{\beta }} I_m-D_B&\quad 0 \\ 0&\quad {\hat{\alpha }} I_n-D_A \end{bmatrix}+Q. \end{aligned}$$

Noting (3.5) and using Lemmas A.2 and A.4, we see

$$\begin{aligned} \begin{bmatrix} \widehat{E}_0&\quad Y_0\\ X_0&\quad \widehat{F}_0 \end{bmatrix}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,(\gamma _1 I-Q)^{-1} \end{aligned}$$

for $\gamma _1>\rho (Q)$. By Lemmas A.2 and A.3,

$$\begin{aligned} \begin{bmatrix} I_m&\quad -Y_0\\ -X_0&\quad I_n \end{bmatrix}^{-1} \,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,\left( \begin{bmatrix}\gamma _2 I_m&\quad 0 \\ 0&\quad \gamma _2 I_n \end{bmatrix}-\begin{bmatrix} \widehat{E}_0&\quad Y_0\\ X_0&\quad \widehat{F}_0 \end{bmatrix} \right) ^{-1}\,{\mathop {=}\limits ^{\scriptstyle 0}}\,(\gamma _1 I-Q)^{-1} \end{aligned}$$

for $\gamma _2$ large enough. Using (3.20) with $k=0$ and noting

$$\begin{aligned} \begin{bmatrix} 0&\quad Y_0\\ X_0&\quad 0 \end{bmatrix} \,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,(\gamma _1I-Q)^{-1},\quad \begin{bmatrix} \widehat{E}_0&\quad 0 \\ 0&\quad \widehat{F}_0 \end{bmatrix}\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,(\gamma _1I-Q)^{-1}, \end{aligned}$$

we have

$$\begin{aligned} \begin{bmatrix} \widehat{E}_1&\quad Y_1 \\ X_1&\quad \widehat{F}_1 \end{bmatrix}&\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,(\gamma _1I-Q)^{-1}+(\gamma _1I-Q)^{-1}(\gamma _1I-Q)^{-1}(\gamma _1I-Q)^{-1}\\&\,{\mathop {=}\limits ^{\scriptstyle 0}}\,(\gamma _1I-Q)^{-1}+(\gamma _1I-Q)^{-1}\\&\,{\mathop {=}\limits ^{\scriptstyle 0}}\,(\gamma _1I-Q)^{-1} \\&\,{\mathop {=}\limits ^{\scriptstyle 0}}\,\begin{bmatrix} \widehat{E}_0&\quad Y_0\\ X_0&\quad \widehat{F}_0 \end{bmatrix}. \end{aligned}$$

This proves (A.4) for $k=0$.

Now suppose (A.4) holds for $k=\ell$. We can show that it also holds for $k=\ell +1$ by using the following majorizations in the entrywise nonzero pattern

$$\begin{aligned} \begin{bmatrix} \widehat{E}_{\ell +1}&\quad 0 \\ 0&\quad \widehat{F}_{\ell +1} \end{bmatrix}\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,\begin{bmatrix} \widehat{E}_{\ell }&\quad 0 \\ 0&\quad \hat{F}_{\ell } \end{bmatrix},\quad \begin{bmatrix} 0&\quad Y_{\ell +1} \\ X_{\ell +1}&\quad 0 \end{bmatrix}\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,\begin{bmatrix} 0&\quad Y_{\ell } \\ X_{\ell }&\quad 0 \end{bmatrix} \\ \begin{bmatrix} I_m&\quad -Y_{\ell +1} \\ -X_{\ell +1}&\quad I_n \end{bmatrix}^{-1}\,{\mathop {\preceq }\limits ^{\scriptstyle 0}}\,\begin{bmatrix} I_m&\quad -Y_{\ell } \\ -X_{\ell }&\quad I_n \end{bmatrix}^{-1}. \end{aligned}$$

By induction, (A.4) holds for all k. $\square$