Accurate Numerical Solution for Shifted M-Matrix Algebraic Riccati Equations

Abstract

An algebraic Riccati equation (are) is called a shifted M-matrix algebraic Riccati equation (mare) if it can be turned into an mare after its matrix variable is partially shifted by a diagonal matrix. Such an are can arise from computing the invariant density of a Markov modulated Brownian motion. Sufficient and necessary conditions for an are to be a shifted mare are obtained. Based on the conditions, a highly accurate implementation of the alternating directional doubling algorithm (adda) is established to compute the extremal solution of a shifted mare, as well as a quantity needed for computing the invariant density in the application, with high entrywise relative accuracy. Numerical examples are presented to demonstrate the theory and algorithms.

Appendix

1.1 Markov Modulated Brownian Motion

Consider a general Markov modulated Brownian motion (MMBM) $\{(F(t),J(t)):\ t\ge 0\}$, where F(t) is the level of fluid in a container at time t and J(t) is the state at time t of a continuous time Markov chain modulating the level process. Let Q be the infinitesimal generator of the Markov chain which has a finite state space partitioned as

$$\begin{aligned} \mathbb {S}=\mathbb {S}_b\cup \mathbb {S}_u\cup \mathbb {S}_d\cup \mathbb {S}_0. \end{aligned}$$

(A.1)

With respect to this state space partition, the dynamics of the fluid level can be described as follows:

• 1. when the chain is in a state $i\in \mathbb {S}_b$, the level evolves according to a Brownian motion process with drift $r_i$ and diffusion coefficient $\frac{1}{2}v_i>0$, where $r_i$ can be positive, negative, or 0;

• 2. when the chain is in a state $i\in \mathbb {S}_u$, the level increases at rate $|r_i|$ with $r_i>0$;

• 3. when the chain is in a state $i\in \mathbb {S}_d$, the level decrease at rate $|r_i|$ with $r_i<0$, provided the container is nonempty;

• 4. when the chain is in a state $i\in \mathbb {S}_0$, the level remains unchanged.

When the Markov chain falls within the subset of states $\mathbb {S}_u\cup \mathbb {S}_d\cup \mathbb {S}_0$, the fluid level changes very much like a Markov modulated fluid flow model (MMFF).

In what follows, we adopt the convention that $v_i=0$ for $i\in \mathbb {S}_u\cup \mathbb {S}_d\cup \mathbb {S}_0$, and, without loss of generality, we may assume that $\mathbb {S}_0=\emptyset$ because it can be censored out [1]. Therefore, instead of (A.1), we will have

$$\begin{aligned} \mathbb {S}=\mathbb {S}_b\cup \mathbb {S}_u\cup \mathbb {S}_d. \end{aligned}$$

(A.2)

Let $N=|\mathbb {S}|$, the cardinality of $\mathbb {S}$, i.e., the number of states in $\mathbb {S}$, and, for future references,

$$\begin{aligned} n_b=|\mathbb {S}_b|,\,\, n_u=|\mathbb {S}_u|,\,\, n_d=|\mathbb {S}_d|,\,\, N=|\mathbb {S}|=n_b+n_u+n_d, \end{aligned}$$

(A.3a)

$$\begin{aligned} m=n_b+n_d,\,\, n=n_b+n_u. \end{aligned}$$

(A.3b)

Besides Q, let

$$\begin{aligned} {\hat{R}}={{\,\mathrm{diag}\,}}(r_i)_{i\in \mathbb {S}},\quad V={{\,\mathrm{diag}\,}}\left( {v_i}/{2}\right) _{i\in \mathbb {S}}. \end{aligned}$$

We assume the infinitesimal generator Q is irreducible and denote by $\pmb {\pi }\in \mathbb {R}^N$ its stationary distribution, i.e. $\pmb {\pi }> 0$, $\pmb {\pi }^{{{\,\mathrm{T}\,}}} Q=0$ and $\pmb {\pi }^{{{\,\mathrm{T}\,}}}\pmb {1}=1$. In the model, $\nu =\pmb {\pi }^{{{\,\mathrm{T}\,}}}{\hat{R}}\mathbf {1}<0$ which means the fluid level is positive recurrent and has a steady-state distribution.

We are interested in computing its invariant density vector $\pmb {p}: \mathbb {R}_{0+}\rightarrow \mathbb {R}_{0+}^{1\times N}$. It satisfies the following second-order ordinary differential equation

$$\begin{aligned} \pmb {p}''(x)V-\pmb {p}'(x){\hat{R}}+\pmb {p}(x)Q=0, \end{aligned}$$

(A.4)

with suitable boundary conditions. It has been proved [17, 20] that $\pmb {p}(x)$ takes the matrix-exponential form

$$\begin{aligned} \pmb {p}(x)=\pmb {c}e^{Kx}\begin{bmatrix} I_{n}&\Gamma \end{bmatrix}, \end{aligned}$$

(A.5)

where $\pmb {c}\in \mathbb {R}^{1\times n}$ can be determined by the boundary conditions, $K\in \mathbb {R}^{n\times n}$ and $\Gamma \in \mathbb {R}^{n\times n_d}$. The pair $(K,\ \Gamma )$ has a probabilistic meaning: $\Gamma \ge 0$ is the matrix recording the first-return probabilities of the time-reversed process, and K is the sub-generator matrix for the downward-record process with the property that $-K$ is a Z-matrix and $K\pmb {1}\le 0$, which implies that $-K$ is an M-matrix.

Substituting (A.5) into (A.4) gives

$$\begin{aligned} K^2\begin{bmatrix} I&\Gamma \end{bmatrix}V-K\begin{bmatrix} I&\Gamma \end{bmatrix}{\hat{R}}+\begin{bmatrix} I&\Gamma \end{bmatrix}Q=0, \end{aligned}$$

(A.6)

where unknowns are K and $\Gamma$. The computation of $\pmb {p}$ is thus reduced to solving the matrix Eq. (A.6). Nguyen and Poloni [24] proposed to first transform (A.6) into a matrix equation of the form $X^2\widehat{A}-X\widehat{B}+\widehat{C}=0$ with $\widehat{A},\ \widehat{C}\in \mathbb {R}_{0+}^{N\times N}$, $\widehat{B}\in \mathbb {R}^{N\times N}$, and $\widehat{B}-\widehat{A}-\widehat{C}$ being a regular M-matrix and then solve it by the cyclic reduction (CR) that can produce a solution with high entrywise relative accuracy. Later K and $\Gamma$ are recovered from the computed solution X. It is an elegant and effective approach. The only drawback is perhaps the transformed equation has a larger scale than (A.6).

Remark 8.1

When $\mathbb {S}=\mathbb {S}_b$, V is nonsingular and (A.6) reduces to $K^2V-K{\hat{R}}+Q=0$, or equivalently, $K^2-K{\hat{R}}V^{-1}+QV^{-1}=0$. Guo [11] investigated the quadratic matrix equation of the form $X^2+XE+F=0$ with a diagonal matrix E and an M-matrix F and showed such an equation has a unique M-matrix solution.

Inspired by the idea in Ahn and Ramaswami [1], in what follows, we seek to transform (A.6) into an $n\times m$are which turns out to be a shifted mare and then propose to solve the shifted mare by the highly accurate ADDA detailed in Sects. 4 and 5.

Consistently with (A.2), we partition Q, V and ${\hat{R}}$ as

$$\begin{aligned} Q=\begin{bmatrix} Q_{bb} &{}\quad Q_{bu} &{}\quad Q_{bd}\\ Q_{ub} &{}\quad Q_{uu} &{} \quad Q_{ud}\\ Q_{db} &{} \quad Q_{du} &{}\quad Q_{dd} \end{bmatrix},\quad V=\begin{bmatrix} V_b &{}\quad &{}\quad \\ &{}\quad 0&{}\quad \\ &{}\quad &{}\quad 0 \end{bmatrix},\quad \widehat{R}=\begin{bmatrix} \widehat{R}_b &{}\quad &{}\quad \\ &{}\quad \widehat{R}_u &{}\quad \\ &{} \quad &{}\quad \widehat{R}_d \end{bmatrix}, \end{aligned}$$

(A.7)

where $V_b={{\,\mathrm{diag}\,}}(v_i)_{i\in \mathbb {S}_b}$ with all $v_i>0$, $\widehat{R}_b={{\,\mathrm{diag}\,}}(r_i)_{i\in \mathbb {S}_b}$, $\widehat{R}_u={{\,\mathrm{diag}\,}}(r_i)_{i\in \mathbb {S}_u}$ with all $r_i>0$, and $\widehat{R}_d={{\,\mathrm{diag}\,}}(r_i)_{i\in \mathbb {S}_d}$ with all $r_i<0$. Similarly, corresponding to $\mathbb {S}_b\cup \mathbb {S}_u$, we partition K and $\Gamma$ as

Then (A.6) can be expressed as

$$\begin{aligned}&\begin{bmatrix} K_{11} &{}\quad K_{12}\\ K_{21} &{}\quad K_{22} \end{bmatrix}^2 \begin{bmatrix} I &{}\quad 0 &{}\quad \Gamma _{1}\\ 0 &{}\quad I &{}\quad \Gamma _{2} \end{bmatrix} \begin{bmatrix} V_b &{}\quad &{}\quad \\ &{}\quad 0&{}\quad \\ &{}\quad &{}\quad 0 \end{bmatrix} -\begin{bmatrix} K_{11} &{}\quad K_{12}\\ K_{21} &{}\quad K_{22} \end{bmatrix}\begin{bmatrix} I &{}\quad 0 &{}\quad \Gamma _{1}\\ 0 &{} \quad I &{}\quad \Gamma _{2} \end{bmatrix}\begin{bmatrix} \widehat{R}_b &{}\quad &{}\quad \\ &{}\quad \widehat{R}_u &{}\quad \\ &{}\quad &{}\quad \widehat{R}_d \end{bmatrix}\nonumber \\&\quad +\begin{bmatrix} I &{} \quad 0 &{}\quad \Gamma _{1}\\ 0 &{}\quad I &{}\quad \Gamma _{2} \end{bmatrix}\begin{bmatrix} Q_{bb} &{}\quad Q_{bu} &{}\quad Q_{bd}\\ Q_{ub} &{}\quad Q_{uu} &{}\quad Q_{ud}\\ Q_{db} &{}\quad Q_{du} &{}\quad Q_{dd} \end{bmatrix}=0, \end{aligned}$$

(A.8)

where the unknowns are all $K_{ij}$ and $\Gamma _i$. To reduce the number of unknowns, we claim that $\begin{bmatrix} K_{12} \\ K_{22} \end{bmatrix}$ can be obtained by some affine transformation on $\begin{bmatrix} \Gamma _{1}\\ \Gamma _{2} \end{bmatrix}$. In fact, equating the second block columns on both sides of (A.8), we arrive at

$$\begin{aligned} -\begin{bmatrix} K_{11} &{}\quad K_{12}\\ K_{21} &{}\quad K_{22} \end{bmatrix}\begin{bmatrix} I &{}\quad 0 &{}\quad \Gamma _{1}\\ 0 &{}\quad I &{}\quad \Gamma _{2} \end{bmatrix}\begin{bmatrix} 0\\ \widehat{R}_u\\ 0 \end{bmatrix}+\begin{bmatrix} I &{}\quad 0 &{}\quad \Gamma _{1}\\ 0 &{}\quad I &{}\quad \Gamma _{2} \end{bmatrix}\begin{bmatrix} Q_{bu}\\ Q_{uu}\\ Q_{du} \end{bmatrix}=0, \end{aligned}$$

or equivalently,

$$\begin{aligned} -\begin{bmatrix} K_{12}\\ K_{22} \end{bmatrix} \widehat{R}_u+\begin{bmatrix} Q_{bu}\\ Q_{uu} \end{bmatrix}+\begin{bmatrix} \Gamma _{1}\\ \Gamma _{2} \end{bmatrix} Q_{du}=0 \quad \Rightarrow \quad \begin{bmatrix} K_{12}\\ K_{22} \end{bmatrix}=\begin{bmatrix} Q_{bu}\\ Q_{uu} \end{bmatrix}\widehat{R}_u^{-1}+\begin{bmatrix} \Gamma _{1}\\ \Gamma _{2} \end{bmatrix}Q_{du}\widehat{R}_u^{-1}. \end{aligned}$$

(A.9)

Thus $\begin{bmatrix} K_{12}\\ K_{22} \end{bmatrix}$ can be eliminated from (A.8), leaving $\begin{bmatrix} K_{11}\\ K_{21}\end{bmatrix}$ and $\begin{bmatrix}\Gamma _1\\ \Gamma _2\end{bmatrix}$ as the ones to be determined. Next, we group the unknowns into one⁶

$$\begin{aligned} X:= \begin{bmatrix} K_{11} &{} \quad \Gamma _{1}\\ K_{21} &{}\quad \Gamma _{2} \end{bmatrix}\in \mathbb {R}^{n\times m} \end{aligned}$$

(A.10)

and seek to establish an equation in X, based on (A.8). To this end, we extract the first and the third block columns in (A.8) to yield

$$\begin{aligned}&\begin{bmatrix} K_{11} &{}\quad K_{12}\\ K_{21} &{}\quad K_{22} \end{bmatrix}^2\begin{bmatrix} I &{} \quad 0 &{}\quad \Gamma _{1}\\ 0 &{}\quad I &{}\quad \Gamma _{2} \end{bmatrix}\begin{bmatrix} V_b &{}\quad 0 \\ 0 &{}\quad 0\\ 0&{}\quad 0 \end{bmatrix} -\begin{bmatrix} K_{11} &{} \quad K_{12}\\ K_{21} &{}\quad K_{22} \end{bmatrix}\begin{bmatrix} I &{} \quad 0 &{}\quad \Gamma _{1}\\ 0 &{}\quad I &{} \quad \Gamma _{2} \end{bmatrix}\begin{bmatrix} \widehat{R}_b &{}\quad 0 \\ 0 &{}\quad 0\\ 0&{}\quad \widehat{R}_d \end{bmatrix}\nonumber \\&\quad +\begin{bmatrix} I &{} \quad 0 &{}\quad \Gamma _{1}\\ 0 &{}\quad I &{}\quad \Gamma _{2} \end{bmatrix}\begin{bmatrix} Q_{bb} &{}\quad Q_{bd}\\ Q_{ub} &{}\quad Q_{ud}\\ Q_{db} &{} \quad Q_{dd} \end{bmatrix}=0. \end{aligned}$$

(A.11)

Using (A.9), we get

$$\begin{aligned} \begin{bmatrix} K_{11} &{}\quad K_{12}\\ K_{21} &{}\quad K_{22} \end{bmatrix}=\begin{bmatrix} 0_{n_b\times n_b} &{}\quad Q_{bu}\widehat{R}_u^{-1}\\ 0 &{}\quad Q_{uu}\widehat{R}_u^{-1} \end{bmatrix}+\begin{bmatrix} K_{11} &{}\quad \Gamma _{1}\\ K_{21} &{}\quad \Gamma _{2} \end{bmatrix}\begin{bmatrix} I &{}\quad \\ &{}\quad Q_{du}\widehat{R}_u^{-1} \end{bmatrix}, \end{aligned}$$

(A.12)

plugging which into (A.11), we finally obtain an are in the form of (1.1) with

$$\begin{aligned} A&=\begin{bmatrix} 0_{n_b\times n_b} &{}\quad -\,Q_{bu}\widehat{R}_u^{-1}\\ 0 &{}\quad -\,Q_{uu}\widehat{R}_u^{-1} \end{bmatrix}, \quad B=\begin{bmatrix} \widehat{R}_bV_b^{-1}&{}\quad \\ -\,Q_{db}V_b^{-1}&{}\quad -\,Q_{dd}(-\widehat{R}_d)^{-1} \end{bmatrix}, \end{aligned}$$

(A.13a)

$$\begin{aligned} C&=\begin{bmatrix} Q_{bb}V_b^{-1} &{}\quad Q_{bd}(-\widehat{R}_d)^{-1}\\ Q_{ub}V_b^{-1} &{}\quad Q_{ud}(-\widehat{R}_d)^{-1} \end{bmatrix}, \quad D=\begin{bmatrix} I &{}\quad \\ &{}\quad Q_{du}\widehat{R}_u^{-1} \end{bmatrix}. \end{aligned}$$

(A.13b)

We now take a close look at X in (A.10). The matrices $K_{21}$, $\Gamma _{1}$ and $\Gamma _2$ are nonnegative, while $K_{11}$ has nonnegative off-diagonal entries and negative diagonal entries. If we add a sufficiently large diagonal matrix, say $\varLambda$, to $K_{11}$, then the resulting matrix

$$\begin{aligned} X_{\varOmega }:= \begin{bmatrix} K_{11}+\varLambda &{}\quad \Gamma _{1}\\ K_{21} &{}\quad \Gamma _{2} \end{bmatrix}=\begin{bmatrix} K_{11} &{} \quad \Gamma _{1}\\ K_{21} &{}\quad \Gamma _{2} \end{bmatrix}+\begin{bmatrix} \varLambda &{}\quad 0\\ 0 &{}\quad 0 \end{bmatrix}=:X+\varOmega \end{aligned}$$

(A.14)

can be made nonnegative. Plug $X=X_{\varOmega }-\varOmega$ into (1.1) to yield an are in $X_{\varOmega }$ in the form of (1.5) with $p=n_b$. The associated matrix of are (1.5)

(A.15)

is a nonsingular M-matrix (for sufficiently large $\varLambda$), making (1.5) an mare, and, as a result, are (1.1) with (A.13) is a shifted mare. The interesting solution $X_{\varOmega }$ corresponding to the interesting one of (1.1) for the underlying application turns out to be the unique minimal nonnegative solution of (1.5). This is because $-K$ is an M-matrix and, by (A.9) and (A.10),

$$\begin{aligned} -K=(A-XD)=A_{\varOmega }-X_{\varOmega } D, \end{aligned}$$

(A.16)

following the solution properties of an mare as summarized in Proposition 2.2.

Remark 8.2

A few comments are in order about are (1.1) with (A.13).

• (a) It can be verified that A, B, C, and D given by (A.13) admit the forms of those in (3.2) with the properties specified in Theorem 3.2.

• (b) Assumption 4.1(a), i.e., that the diagonal matrices $A_{11}$, $B_{11}$ and $D_{11}$ are accurately known with high entrywise relative accuracy, is satisfied by those in (A.13). In particular, $A_{11}=0$ and $D_{11}=I$ here.

• (c) To realize Assumption 4.1(b), we notice

  $$\begin{aligned} W_0=-\begin{bmatrix} Q_{dd}&{}\quad Q_{db} &{}\quad Q_{du}\\ Q_{bd}&{}\quad Q_{bb} &{}\quad Q_{bu}\\ Q_{ud}&{}\quad Q_{ub} &{}\quad Q_{uu} \end{bmatrix}\begin{bmatrix} -\,{\hat{R}}_d^{-1} &{}\quad &{}\quad \\ &{}\quad V_b^{-1} &{}\quad \\ &{}\quad &{} \quad {\hat{R}}_u^{-1} \end{bmatrix}. \end{aligned}$$

  (A.17)

  Consistently with (A.2), we partition the stationary distribution $\pmb {\pi }$ of the infinitesimal generator Q as $\pmb {\pi }=\begin{bmatrix} \pmb {\pi }_b&\pmb {\pi }_u&\pmb {\pi }_d \end{bmatrix}\in \mathbb {R}^{1\times N}$ which can be computed with high entrywise relative accuracy by the GTH algorithm [25, Theorem 2]. We have $\begin{bmatrix} \pmb {\pi }_d&\pmb {\pi }_b&\pmb {\pi }_u \end{bmatrix}W_0=0$, a very natural thing to have. This leads to a left triplet representation of $W_0$, as needed for computing $\varPhi _{\varOmega }$ and $K=-(A-\varPhi D)$ entrywise accurately by the algorithms in Sects. 4 and 5.

What we have described so far is for a type of shifted mares as arising from computing the invariant density vector $\pmb {p}$ of MMBM as determined by (A.4). It is about a steady-state analysis. In the time-dependent analysis of MMBM, Ahn and Ramaswami [1] established ares of the first passage time distribution. These ares provide ample examples of shifted mares, too. They are in the form of (1.1) with

$$\begin{aligned} A&=\begin{bmatrix} -\,\widehat{R}_bV_b^{-1} &{}\quad -\,\sqrt{2}V_b^{-\frac{1}{2}}Q_{bu}\\ 0 &{}\quad -\,\widehat{R}_u^{-1}(Q_{uu}-sI) \end{bmatrix}, \quad B=\begin{bmatrix} 0&{}\quad 0 \\ \widehat{R}_d^{-1}Q_{db}&{}\quad \widehat{R}_d^{-1}(Q_{dd}-sI) \end{bmatrix}, \end{aligned}$$

(A.18a)

$$\begin{aligned} C&=\begin{bmatrix} \sqrt{2}V_b^{-\frac{1}{2}}(Q_{bb}-sI) &{}\quad \sqrt{2}V_b^{-\frac{1}{2}}Q_{bd}\\ \widehat{R}_u^{-1}Q_{ub}&{}\quad \widehat{R}_u^{-1}Q_{ud} \end{bmatrix}, \quad D=\begin{bmatrix} (2V_b)^{-\frac{1}{2}} &{}\quad \\ &{}\quad (-\widehat{R}_d)^{-1}Q_{du} \end{bmatrix}, \end{aligned}$$

(A.18b)

where $s\in \mathbb {C}_{0+}$ in general but we will consider $s\ge 0$ only. Now we shift the corresponding are (1.1) as in (1.4) – (1.7) with $\varLambda ={\hat{R}}_b(2V_b)^{-1/2}+\sqrt{{\hat{R}}_b^2(2V _b)^{-1}+2\varLambda _{Q_{bb}}+2sI}$ to get

$$\begin{aligned} A_{\varOmega }&=\begin{bmatrix} -\,\widehat{R}_b(2V_b)^{-1}+(2V_b)^{-\frac{1}{2}}\sqrt{(2V_b)^{-1}\widehat{R}_b^2+2\varLambda _{Q_{b,b}}+2sI} &{}\quad -\,\sqrt{2}V_b^{-\frac{1}{2}}Q_{bu}\\ 0 &{}\quad -\,\widehat{R}_u^{-1}(Q_{uu}-sI) \end{bmatrix},\end{aligned}$$

(A.19a)

$$\begin{aligned} B_{\varOmega }&=\begin{bmatrix} (2V_b)^{-1}\widehat{R}_b+(2V_b)^{-\frac{1}{2}}\sqrt{(2V_b)^{-1}\widehat{R}_b^2+2\varLambda _{Q_{b,b}}+2sI} &{}\quad 0 \\ \widehat{R}_d^{-1}Q_{db}&{}\quad \widehat{R}_d^{-1}(Q_{dd}-sI) \end{bmatrix},\ \end{aligned}$$

(A.19b)

$$\begin{aligned} C_{\varOmega }&=\begin{bmatrix} \sqrt{2}V_b^{-\frac{1}{2}}(Q_{bb}+\varLambda _{Q_{b,b}}) &{} \quad \sqrt{2}V_b^{-\frac{1}{2}}Q_{bd}\\ \widehat{R}_u^{-1}Q_{ub}&{} \quad \widehat{R}_u^{-1}Q_{ud} \end{bmatrix}, \end{aligned}$$

(A.19c)

where $\varLambda _{Q_{bb}}={{\,\mathrm{diag}\,}}(-[Q_{bb}]_{(i,i)})_{i\in \mathbb {S}_b}$. The resulting are (1.5) is the same as the one in [1, Theorem 5.1] which provably is an mare. Ahn and Ramaswami [1] suggested to solve the mare, as a nonlinear equation, by Newton’s method.

Remark 8.3

The first two comments in Remark 8.2 still applies to are (1.1) with (A.18). To realize Assumption 4.1(b), we now have

$$\begin{aligned} W_0=&\begin{bmatrix} -\,\widehat{R}_d^{-1} &{}\quad &{}\quad \\ &{}\quad \sqrt{2} V_b^{-1/2} &{}\quad \\ &{}\quad &{}\quad \widehat{R}_u^{-1}\end{bmatrix} \begin{bmatrix} sI-Q_{dd} &{}\quad -\,Q_{db}&{}\quad -\,Q_{du}\\ -\,Q_{bd} &{}\quad sI-Q_{bb} &{}\quad -\,Q_{bu}\\ -\,Q_{ud} &{}\quad -\,Q_{ub} &{}\quad sI-Q_{uu} \end{bmatrix} \nonumber \\&\times \begin{bmatrix} I_{n_d} &{}\quad &{}\quad \\ &{}\quad (2V_b)^{-1/2} &{}\quad \\ &{}\quad &{}\quad I_{n_u} \end{bmatrix}. \end{aligned}$$

(A.20)

Again let $\pmb {\pi }$ be as in Remark 8.2(c). Then $\pmb {u}^{{{\,\mathrm{T}\,}}} W_0=\pmb {v}^{{{\,\mathrm{T}\,}}}$, where

$$\begin{aligned} \pmb {u}^{{{\,\mathrm{T}\,}}}=\begin{bmatrix} \pmb {\pi }_d&\pmb {\pi }_b&\pmb {\pi }_u \end{bmatrix}{{\,\mathrm{diag}\,}}\Big (-\widehat{R}_d,\frac{1}{\sqrt{2}}V_b^{1/2}, \widehat{R}_u\Big ),\quad \pmb {v}^{{{\,\mathrm{T}\,}}}=s\begin{bmatrix} \pmb {\pi }_d,&\pmb {\pi }_b(2V_b)^{-1/2},&\pmb {\pi }_u \end{bmatrix}. \end{aligned}$$

(A.21)

We commented that $\pmb {\pi }$ can be computed with high entrywise relative accuracy. Since $\widehat{R}_d$, $V_b$, and $\widehat{R}_u$ are diagonal, $\pmb {u}$ and $\pmb {v}$ as in (A.21) can be computed with high entrywise relative accuracy, too.