from Transport Theory
Peter Chang-Yi Weng
∗Hung-Yuan Fan
†Eric King-wah Chu
‡Abstract
We consider the solution of the large-scale nonsymmetric algebraic Riccati equation XCX − XD − AX + B = 0 from transport theory (Juang 1995), with M ≡ [D, −C; −B, A] ∈ R2n×2n being a nonsingular M-matrix. In addition, A, D are rank-1 updates of diagonal matrices, with the products A−1u, A−>u, D−1v and D−>v computable in O(n) complexity, for some vectors u and v, and B, C are rank 1. The structure-preserving doubling algorithm by Guo, Lin and Xu (2006) is adapted, with the appropriate applications of the Sherman-Morrison-Woodbury formula and the sparse-plus-low-rank representations of various iter-ates. The resulting large-scale doubling algorithm has an O(n) computational complexity and memory requirement per iteration and converges essentially quadratically, as illustrated by the numerical examples.
Keywords. doubling algorithm, Krylov subspace, M-matrix, nonsymmetric algebraic Riccati equation
AMS subject classifications. 15A24, 65F50
1
Introduction
Consider the nonsymmetric algebraic Riccati equation (NARE)
R(X) ≡ XCX − XD − AX + B = 0, (1) with A = ∆ − eq>, B = ee>, C = qq>, D = Γ − qe> being n × n, ∆ = diag{d1, · · · , dn},
Γ = diag{γ1, · · · , γn}, e = [1, · · · , 1]>, q = [q1, · · · , qn]> and di, γi, qi > 0 (i = 1, · · · , n). From
[19], for the solvability of (1), we assume that the matrix
M = D −C −B A ∈ R2n×2n
is a nonsingular M-matrix. We are interested in the minimal nonnegative solution X [19]. Moti-vated by the low-ranked cases from the applications in transport theory [19, 23, 24], we further
∗School of Mathematical Sciences, Building 28, Monash University 3800, Australia; [email protected],
phone: +61399054465, FAX: +61399054403
†(Corresponding author) Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan;
‡School of Mathematical Sciences, Building 28, Monash University 3800, Australia; [email protected]
assume that n is large, thus A and D are sparse-like (with the products A−1u, A−>u, D−1v and D−>v computable in O(n) complexity, for some vectors u and v). We adapt the structure-preserving doubling algorithm (SDA) [17] to the NARE in (1), resulting in a large-scale doubling algorithm (SDA ls) with an O(n) computational complexity and memory requirement per itera-tion. Note that the orthodox SDA in [17] has a computational complexity of O(n3) per iteration. Also, the Newton’s algorithm and the iterative methods in [4, 16] are of O(n2) in computational
complexity per iteration, and appear to be incapable of benefitting from the low-rank structures in (1) and the numerically low-ranked X (implied constructively by the SDA; see the discussion at the end of Section 2.1). The fact that X = T ◦ (uv>), the Hadamard product of the Cauchy matrix T =(di+ γj)−1 and the rank-1 matrix uv>, is clever and appropriate in previous work
but is inhibiting for problems of large size n here.
More generally, algebraic Riccati equations arise in many important applications, including the total least squares problems with or without symmetric constraints [9], the spectral factorizations of rational matrix functions [10], the linear and nonlinear optimal controls [3], the contractive rational matrix functions [21], the structured complex stability radius [18], the transport theory [19, 23, 24], the Wiener-Hopf factorization of Markov chains [36], and the optimal solutions of linear differential systems [22]. Symmetric algebraic Riccati equations have been the topic of extensive research, and the theory, applications and numerical solution of these equations are the subject of the monographs [22, 32]. The minimal positive solution to the NARE in (1), for medium size problems without the sparseness and low-rank assumptions, has been studied recently by many authors, employing functional iterations, Newton’s method and the structure-preserving algorithm; for example, see [1, 4, 12, 13, 14, 15, 16, 17, 19, 23, 24, 29, 30, 33] and references therein. There have also been generalizations of (1) in [23, 24], for one-dimensional multi-state and two-dimensional problems.
2
Large-Scale Doubling Algorithm
2.1
Main Ideas
Borrowing from [25, 26, 27], the main ideas behind the algorithm are:
(1) The appropriate application of the Sherman-Morrison-Woodbury formula in order avoid the inversion of large or unstructured matrices.
(2) The use of sparse-plus-low-rank matrices.
(3) The computation of matrix operators recursively, to preserve the corresponding sparsity or low-rank structures, instead of forming them explicitly.
(4) The compression and truncation of fast growing yet diminishing low-ranked components in the iterates, to trade off negligible amount of accuracy for better efficiency.
(5) The careful organization of convergence control in the algorithm, so as to preserve the O(n) computational complexity and memory requirement per iteration.
Recall the low-rank decompositions
The SDA [17] has the form E0 = Eτ, F0= Fτ, G0= Gτ, H0= Hτ; Ek+1 = Ek(I − GkHk)−1Ek, Fk+1= Fk(I − HkGk)−1Fk, Gk+1 = Gk+ Ek(I − GkHk)−1GkFk, Hk+1= Hk+ Fk(I − HkGk)−1HkEk; (2) with Eτ= I − 2τ Vτ−1, Gτ = 2τ D−1τ CWτ−1, Fτ = I − 2τ Wτ−1, Hτ= 2τ Wτ−1BDτ−1, (3)
where the parameter τ ≥ maxi,j{aii, djj} and
Aτ = A + τ I, Dτ = D + τ I, Wτ = Aτ− BDτ−1C, Vτ= Dτ− CA−1τ B. (4)
From (2), we can summarize our algorithm in the following sentence. It is evident that (I − GkHk)−1 and (I − HkGk)−1 are low-rank updates of the identity, and Ek, Fk Gk and
Hk are low-rank updates of Ek−12 , F 2
k−1, Gk−1 and Hk−1, respectively. This, together with the
compression and truncation of the Krylov subspaces described in Section 2.3, form the basis of the SDA for the large-scale NARE in (1).
Similar to the approach in [25, 26, 27], we apply the Sherman-Morrison-Woodbury formula (SMWF):
( eA ± U eCV )−1= eA−1∓ eA−1U ( eC−1± V eA−1U )−1V eA−1,
to various inverses of matrices in sparse-plus-low-rank (splr) form, required in the SDA. With A, D being sparse-like, as Aτ and Dτ can be “inverted” (or the corresponding linear equations
solved) efficient in O(n) complexity. Consequently, we obtain the rank-1 updates of the inverses of Aτ and Dτ, respectively,
Wτ−1 = A−1τ + cW · A−1τ e · q >A−1
τ
Vτ−1 = D−1τ + cV · Dτ−1q · e>Dτ−1, (5)
with the constants
cW ≡ δτ 1 − δτατ , cV ≡ ατ 1 − ατδτ ; δτ ≡ e>Dτ−1q, ατ ≡ q>A−1τ e.
With ∆τ ≡ ∆ + τ I and Γτ ≡ Γ + τ I, note that
A−1τ = ∆−1τ + cA· ∆−1τ e · q >∆−1 τ , cA≡ (1 − eδτ)−1, eδτ ≡ q>∆−1τ e; D−1τ = Γ−1τ + cD· Γ−1τ q · e>Γ−1τ , cD≡ (1 −γeτ) −1, e γτ≡ q>Γ−1τ e. (6)
Consequently, for arbitrary vectors u and v, the products A−1τ u, A−>τ u, Dτ−1v and D−>τ v, or
V−1
τ u, Vτ−>u, Wτ−1v and Wτ−>v, can all be computed in O(n) complexity. Indeed, from (6), we
see that
A−1τ e = cA(∆−1τ e), q>A−1τ = cA(q>∆−1τ ),
D−1τ q = cD(Γ−1τ q), e >D−1
Furthermore, it follows from (5)–(7) that ατ= q>A−1τ e = cA· (q>∆−1τ e) = cAδ˜τ= cA(1 − 1 cA ) = cA− 1, δτ= e>Dτ−1q = cD· (e>Γ−1τ q) = cD˜γτ = cD(1 − 1 cD ) = cD− 1,
respectively. Hence, the scalar coefficints cV and cW in (5) become
cV = ατ 1 − ατδτ = cA− 1 1 − (cD− 1)(cA− 1) = 1 − cA cAcD− cA− cD , cW = δτ 1 − δτατ = cD− 1 1 − (cA− 1)(cD− 1) = 1 − cD cAcD− cA− cD , (8)
respectively. Consequently, it follows from (7) and (8) that
Wτ−1= ∆−1τ + cτ(∆−1τ e)(q >∆−1 τ ), Vτ−1= Γ−1τ + cτ(Γ−1τ q)(e >Γ−1 τ ), (9) with cτ = cAcD cA+ cD− cAcD .
Substituting (9) into (3), we obtain
Eτ = (I − 2τ Γ−1τ ) − 2τ cτ(Γ−1τ q)(e>Γ−1τ ), Fτ = (I − 2τ ∆−1τ ) − 2τ cτ(∆−1τ e)(q>∆−1τ ),
Gτ = 2τ cτ(Γ−1τ q)(q>∆−1τ ), Hτ= 2τ cτ(∆−1τ e)(e>Γ−1τ ). (10)
Note that the matrices Eτ and Fτare in the form of a diagonal matrix plus a rank-1 update, and
Gτ, Hτ are rank-1. All these four matrices can be computed in O(n) complexity, dependent on
only the four vectors ∆−1τ e, ∆−1τ q, Γ−1τ e and Γ−1τ q, as well as two inner products eδτ andeγτ in (6). From [17], as k → ∞, we have the convergence result
Ek, Fk→ 0, Hk→ X, Gk → Y,
with Y being the solution to the adjoint NARE. Together with the formulae derived above, it is clear that Hkand Gk are numerically low-ranked, with diminishing components keep being added
for increasing values of k. This exponential growth of small components in Hk and Gk can be
controlled as in [25, 27], by compressing and truncating the low-ranked components. Importantly, the computation complexity and memory requirement of the resulting algorithm will be O(n) per iteration.
2.2
Algorithm
For k = 0, 1, · · · , we organize the SDA so that the iterates have the recursive forms
Gk = C1kTkC2k>, Hk= B1kRkB2k> (11)
with Eik ∈ Rn×rEk and Fik ∈ Rn×rFk (i = 1, 2), and the kernels Rk ∈ Rmk×mk and Tk ∈
Rlk×lk. We shall see later that r
Ek= mk and rFk= lk from (19a)–(22), after the truncation and
compression of Bik and Cik (i = 1, 2) in Section 2.3. (Note that Hk = B1kRkB2k>, contrasting
with Hk = CkTk−1Ck> in [25, 26, 27], due to the inconsistent use of symbols for different algebraic
Riccati equations and the use of inverses to represent kernels.)
We should compute products like Eku, E>ku, Fkv, Fk>v, for some vectors u and v, by applying
(12) recursively. Without actually forming Ek and Fk, we avoid any possible deterioration of
their sparsity or other structures as k grows, and preserve the O(n) computational complexity of the algorithm. As a trade-off, we need to store all the Bik, Rk, Cik and Tk for all previous k, as
we shall see below.
Applying the SMWF again, we obtain the low-ranked updates of I:
(I − GkHk)−1 = I + GkB1kRk Imk− B > 2kGkB1kRk −1 B2k> (13a) = I + C1k Ilk− TkC > 2kHkC1k −1 TkC2k>Hk, (13b) and (I − HkGk)−1 = I + B1k Imk− RkB > 2kGkB1k −1 RkB>2kGk (14a) = I + HkC1kTk Ilk− C > 2kHkC1kTk −1 C2k>. (14b)
From (2), (13a)–(14b), we can choose the matrices in (12) recursively as
B1,k+1= [B1k, FkB1k], B2,k+1= [B2k, Ek>B2k], (15) Rk+1 = Rk⊕Rk+ (Imk− RkB > 2kGkB1k)−1RkB2k>GkB1kRk (16a) = Rk⊕Rk+ RkB2k>C1kTk(Ilk− C > 2kHkC1kTk)−1C2k>B1kRk , (16b) C1,k+1 = [C1k, EkC1k], C2,k+1= [C2k, Fk>C2k], (17) Tk+1 = Tk⊕Tk+ TkC2k>B1kRk(Imk− B > 2kGkB1kRk)−1B2k>C1kTk (18a) = Tk⊕Tk+ (Ilk− TkC > 2kHkC1k)−1TkC2k>HkC1kTk , (18b) E1,k+1 = EkGkB1kRk Imk− B > 2kGkB1kRk −1 (19a) = EkC1k Ilk− TkC > 2kHkC1k −1 TkC2k>B1kRk, (19b) E2,k+1 = Ek>B2k, (20) and F1,k+1 = FkB1k Imk− RkB > 2kGkB1k −1 RkB2k>C1kTk (21a) = FkHkC1kTk Ilk− C > 2kHkC1kTk −1 , (21b) F2,k+1= Fk>C2k. (22)
Different but equivalent organizations for Eik and Fik (i = 1, 2) are obviously possible.
For the initial values as in (2), we have E0 = Eτ, F0 = Fτ, and G0 = Gγ = C10T0C20> and
H0= Hγ = B10R0B20> with
C10 = Γ−1τ q, C20= ∆−1τ q, T0= −2τ cτI;
B10 = ∆−1τ e, B20= Γ−1τ e, R0= −2τ cτI; (23)
from (10).
Ultimately from (3)–(5), we see that the SDA is dominated by the computation of products like Eku, Ek>u, Fkv, Fk>v, for some vectors u and v. By applying (12) recursively, products like
Eτu, Eτ>u, Fτv, Fτ>v have to be computed. This can be realized using (3)–(5) in O(n) complexity
and memory requirement, because of our assumptions on A, B, C, D.
Similar to the SDA ls in [25, 26, 27] for CAREs, DAREs and large-scale Stein and Lyapunov equations, in our SDA ls for NAREs, the convergence of the doubling iteration competes with the exponential growth in Bik and Cik. The success of the SDA ls depends on the trade-off between
the accuracy of the approximate solution and its CPU-time and memory requirements, controlled by the compression and truncation of Bik and Cikin Section 2.3.
The SDA ls for NAREs, realizing the iteration in (2) with the help of (3), (5), (10), (11), (12), and the rank mkformulae in (15)–(22) (involving small inverses of size mk, ignoring the size
lk alternatives), with initial values in (23), the compression and truncation in (24) and (25) in
Section 2.3, and the convergence control in Section 3.1. We summarize the SDA ls in Algorithm 1 below.
Algorithm 1 (SDA ls)
Input: A, B, C, D ∈ Rn×n; shift τ , positive tolerances τB, τC and , and m
max, lmax;
Output: Bi∈ Rn×m, R∈ Rm×m, Ci∈ Rn×l and T∈ Rl×l (i = 1, 2),
with B1RB2> and C1TC2> approximating, respectively, the solutions
X and Y to the large-scale NARE in (1) and its adjoint; Compute ∆−1τ e, Γ−>τ e, ∆−1τ q, Γτ−>q, e>(Γ−1τ q), q>(∆−1τ e), for
Bi0, Ci0 (i = 1, 2), R0, T0, and E0, F0 (as in (2), (10), (23));
Set k = 0,re0= 2; E0= Eτ, F0= Fτ, G0= Gτ, H0= Hτ;
Do until convergence:
If the relative residual ˜rk = |rk/(hk+mek+ h)| < ,
Set Bi= Bik, R= Rk, Ci= Cikand T= Tk (i = 1, 2); Exit End If Compute B1,k+1= [B1k, FkB1k], B2,k+1= [B2k, Ek>B2k], Rk+1= Rk⊕Rk+ (Imk− RkB > 2kGkB1k)−1RkB>2kGkB1kRk; C1,k+1= [C1k, EkC1k], C2,k+1 = [C2k, Fk>C2k], Tk+1= Tk⊕Tk+ TkC2k>B1kRk(Imk− B > 2kGkB1kRk)−1B>2kC1kTk;
with Ek = E2k−1+ E1kE2k>, Fk= Fk−12 + F1kF2k> (as in (12), (19a)–(22));
Compress Bi,k+1and Ci,k+1(i = 1, 2), using the tolerances τB and τC,
and modify Rk+1and Tk+1, as in (24) and (25);
Compute k ← k + 1, rk, hk andmek, as in (27) in Section 3.1; End Do
2.3
Truncation and Compression of B
ikand C
ikIn Gk= C1kTkC2k> and Hk = B1kRkB2k>, let n B ikand n
C
ikbe the widths of Bik and Cik (i = 1, 2),
respectively. Apply the QR decomposition with column pivoting to obtain (for i = 1, 2)
Bik = QBikM B ik+ eQ B ikMfikB, Cik = QCikM C ik+ eQ C ikMf C ik, with QB ik∈ Rn×r B ik and QC ik∈ Rn×r C ik(i = 1, 2) being unitary, MB ik∈ Rr B ik×n B ik and MC ik∈ Rr C ik×n C ik
being upper triangular,
k fMikBk ≤ τB, k fMC ikk ≤ τ C (i = 1, 2) and rBik = rank Bik≤ nBik≤ mmax n, rCik = rank Cik≤ nCik≤ lmax n.
Here τB and τC are some small positive tolerances controlling the compression and truncation
process, and mmax and lmax are some preset upper bounds.
Truncating Cik and Bik by discarding the small components undere, we have
Hk = QB1kM B 1kRk(M2kB) > (QB 2k) >+ O(τB), Gk = QC1kM1kCTk(M2kC)> (QC2k)>+ O(τC). (24)
By the replacements or re-organization of symbols: (for i = 1, 2)
Bik ← QBik, Rk← M1kBRk(M2kB) >;
Cik ← QCik, Tk← M1kCTk(M2kC)>; (25)
we sacrifice a hopefully small bit of accuracy, in return for the control of the growth in Bik and
Cik, as well as the corresponding demand on CPU-time and computer memory. We can also
restrict Rk and Tk to be square and nonsingular, applying further compression and truncation if
necessary. At the end, we relabel the widths of Bik and Cikas mk and lk (i = 1, 2), respectively.
In the SDA ls, we also prevent any excessive growth of mk and lk by setting the respective upper
bounds mmax and lmax.
2.4
SDA and Krylov Subspaces
There is an interesting relationship between the SDA ls and Krylov subspaces. Define the Krylov subspaces
Kk(A, V ) ≡
span{V } (k = 0), span{V, AV, A2V, · · · , A2k−1V } (k > 0). From (3), (5)–(5), (15) and (19a)–(22), we see that
B10 = Wτ−1e ∈ K0(A−1τ , A−1τ e), B20= D−>τ e ∈ K0(Dτ−>, D−>τ e); B11 = [B10, F0B10] ∈ K1(A−1τ , A −1 τ e), B21= [B20, E0>B20] ∈ K1(Dτ−>, D −> τ e).
(We have abused notations, with V ∈ Kk(A, B) meaning span{V } ∈ Kk(A, B).) Similarly, it is
easy to show that
B1k ∈ Kk(A−1τ , A−1τ e), B2k ∈ Kk(D−>τ , Dτ−>e);
C1k ∈ Kk(Dτ−1, Dτ−1q), C2k∈ Kk(A−>τ , A−>τ q). (26)
Note that, for k = 0, 1, · · · ,
Kk(A−1τ , A −1 τ e) = Kk(∆−1τ , ∆ −1 τ e), Kk(D−>τ , D −> τ e) = Kk(Γ−>τ , Γ −> τ e), Kk(D−1τ , D−1τ q) = Kk(Γ−1τ , Γ−1τ q), Kk(A−>τ , A−>τ q) = Kk(∆−>τ , ∆−>τ q).
In summary, the general SDA is closely related to approximating the solutions X and Y using Krylov subspaces, with additional components vanishing quadratically. However, for problems of moderate size n, Bk and Ck become full-ranked after a few iterations.
2.5
Error of SDA ls
We limit the rank of the approximation to X, trading off the accuracy in X with the efficiency of the SDA ls. Assume that the compression and truncation in (24) and (25) creates errors of O(τ ), τ = max{τB, τC}, in Gk and Hk. It is easy to see from (2) that errors of the same magnitude
will propagate through to Ek+1, Fk+1, Gk+1and Hk+1.
Let δEk, δFk, δGk and δHk be the errors in Ek, Fk, Gk and Hk, respectively, created from
round-off errors or the compression and truncation of Bk and Ck as described in Section 2.3.
From (2), with
N1k ≡ I − GkHk, ∆1k ≡ GkδHk+ δGkHk+ δGkδHk;
N2k ≡ I − HkGk, ∆2k ≡ HkδGk+ δHkGk+ δHkδGk;
and ignoring higher order terms, we have
δEk+1 ≈ δEkN1k−1Ek+ EkN1k−1∆1kN1k−1Ek+ EkN1k−1δEk,
δFk+1 ≈ δFkN2k−1Fk+ FkN2k−1∆2kN2k−1Fk+ FkN2k−1δFk,
δGk+1 ≈ δGk+ δEkN1k−1GkFk+ EkN1k−1∆1kN1k−1GkFk+ EkN1k−1δGkFk+ EkN1k−1GkδFk,
δHk+1 ≈ δHk+ δFkN2k−1HkEk+ FkN2k−1∆2kN2k−1HkEk+ FkN2k−1δHkEk+ FkN2k−1HkδEk.
Let δk ≡ max{kδFkk, kδEkk, kδGkk, kδHkk}, ck ≡ max{kN1k−1k, kN −1 2k k} · max{kFkk, kEkk} and αk≡ max{kGkk, kHkk}, we have k∆1kk, k∆2kk ≤ 2αkδk+ δk2 and δk+1≤ max2ck(1 + αkck), 1 + ck(2αk+ 2ckα2k+ kEkk) · δk+ O(δ2k).
Recall that Fk, Ek→ 0 and Gk, Hk converge to Y, X respectively. Eventually, ck, kEkk converge
to zero, αk converges to a constant and the O(τ ) errors δFk, δEk, δGk and δHk then pass into
δGk+1, δHk+1, δEk+1and δFk+1, creating errors of the same order. Note that ignoring the higher
that Ek, Fk→ 0 will contribute towards diminishing the errors. From our numerical experience,
the trade-off between the ranks of Gk and Hk and the accuracy of the approximate solutions
to X and Y is the key to the success of our computation. If these ranks grow out of control, unnecessary and insignificant small additions to the iterates overwhelm the computation in terms of flop counts and memory requirement. Limiting the ranks will obviously reduce the accuracy of the approximate solution.
An error analysis for the truncated SDA ls can also be written down, applying the theory for Galerkin methods as in [25, 26, 27], different to those in [20]. This difference flows from the fact than SDA ls generates automatically the Krylov subspaces in (26), from which Bik and Cik
(i = 1, 2), thus the approximate solutions Gk and Hk, are generated respectively. In [20], the
Krylov subspaces have to be generated explicitly by some Arnoldi processes, resulting in different convergence control and residuals.
3
Computational Issues
3.1
Residual and Convergence Control
For the convergence control in Algorithm 1, we should compute residuals and differences of iterates carefully in O(n) complexity.
Consider the difference of successive iterates:
δGk≡ C1kTkC2k> − C1,k−1Tk−1C2,k−1> = eC1,kTekCe2,k>
with
e
Ci,k≡ [Cik, Ci,k−1] (i = 1, 2), Tek ≡ Tk⊕ (−Tk−1) .
All the calculations in this subsection involve the norms of similar low-ranked matrices as in δGk. To compute kδGkk efficiently, orthogonalize eCi,k (i = 1, 2) and modify eTk accordingly,
compressing as in (24) and (25) without the truncation, we then have
kδGkk = k eC1,kTekCe2,k> k = k eTkk
in 2- or F-norm. Note that eTk is very small in size relative to δGk.
Similarly, we have
δHk ≡ B1kRkB2k> − B1,k−1Rk−1B>2,k−1= eB1,kRekBe2,k> , with
e
Bi,k≡ [Bik, Bi,k−1] (i = 1, 2), Rek ≡ Rk⊕ (−Rk−1) .
For the residual rk≡ kR(Hk)k of the NARE, the corresponding relative residual equals
e rk ≡ rk hk+mek+ h with hk≡ kAHk+ HkDk, mek ≡ kHkCHkk, h ≡ kBk. We then have rk = kHkCHk− AHk− HkD + Bk = k bB1kRbkBb2k>k,
with ˇ B1k ≡ [AB1k, B1k] , Bˇ2k ≡D>B2k, B2k , Bb1k ≡Bˇ1k, B1 , Bb2k≡Bˇ2k, B2 ; ˘ Rk ≡ RkB2k>CB1kRk, Rˇk≡ 0 Rk Rk 0 , Rbk ≡ 0 −Rk −Rk R˘k ⊕ R.
After orthogonalizing eBi,k, bBik, Bik, Bi, ˘Bik(= Bik) and ˇBikand transforming respectively eRk,
b
Rk, Rk, R, ˘Rk and ˇRk, we have efficient formulae in 2- or F-norm:
kδHkk = k eRkk, rk= k bRkk, kHkk = kRkk,
h = kRk, mek= k ˘Rkk, hk= k ˇRkk. (27)
Using similar techniques, we may replace kA>Hk+ HkDk with kA>Hkk + kHkDk if one desires.
3.2
Operation and Memory Counts
We assume that cτmn flops are required in the products M R or M>R (with M = A, D), or the
solution of M Z = R or M>Z = R (with M = Aτ, Dτ), and R ∈ Rn×m. A start up cost of 10n
flops is made up of the following:
(1) set up Aτ = A + τ I and Dτ = D + τ I, requiring 2n flops; and
(2) set up Bi0 and Ci0 as in (23) with the help of (5) and (10), requiring 8n flops.
The operation and memory counts of Algorithm 1 (SDA ls) for the kth iteration are summa-rized in Table 1 below. In the third column, the number of variables is recorded. Only O(n) operations or memory requirement are included. Note that most of the work is done in the com-putation of Bi,k+1and Ci,k+1(i = 1, 2), for which FkB1k, Ek>B2k, EkC1k and Fk>C2kin (15) and
(17) have to be calculated recursively, as Ekand Fkin (12) are not available explicitly. In Table 1,
we use the notation Nk ≡P k
j=1(lj+ mj). The operation count for the QR decomposition of an
n × r matrix is 4nr2 flops [11, p. 250].
Computation Flops Memory
Bi,k+1, Ci,k+1 2k+2cτ(lk+ mk) + 2(6lk+ 1)Nk n 2Nk+1n
Rk+1, Tk+1 8lkmkn O(l2k+ m2k)
Ei,k+1, Fi,k+1 4lkmkn −
Compress Bi,k+1, Ci,k+1 8(l2k+ m 2 k)n − Modify Rk+1, Tk+1 O(lk3+ m 3 k) − ˜ rk+1 8(2mk+ 1)2+ 40m2k n − Total 2k+2cτ(lk+ mk) + 2(6lk+ 1)Nk+ 2Nk+1n 4(lk2+ m2k+ 3lkmk) + 8(2mk+ 1)2+ 40m2k n
Table 1: Operation and Memory Counts for the kth Iteration in Algorithm 1 (SDA ls)
With lk and mk controlled by the compression and truncation in Section 2.3, the operation
examples in Section 4, the flop count near the end of Algorithm 1 dominates, with the work involved in one iteration approximately doubled that of the previous one. This corresponds to the 2k+2 factor in the total flop count. However, the last iteration is virtually free, as there is no need to prepare for the next iteration.
4
Numerical Examples
The numerical examples were constructed as in [19] with (di, γi, qi) randomly chosen while
sat-isfying the solvability conditions. Only non-critical cases were tested, to illustrate the feasibility of the SDA ls. The critical cases will behave similarly, although with much slower convergence. All the examples are the same except for the values of n = 10000, 50000, 100000, 500000, 1000000. The numerical results were computed using MATLAB [31] Version R2010b, on a HP SL390s workstation with a 3.60GHz Intel Xeon X5687 4-core processor and 48GB RAM, with machine accuracy eps = 2.22 × 10−16.
We have including tables containing kδHkk/kHkk, erk, mk, δtk (the CPU-time for the kth iteration) and tk (≡ P
k
i=1δti, the total CPU-time up to the kth iteration), for various k. As
discussed before, the operation counts are dominated by that for the last but one iteration, as each iteration doubles the previous one in work. The last iteration is virtually free because there is no need to prepare the updated Krylov bases for the next iteration.
Comparing with CAREs of similar sizes, at a particular iteration, the CPU-times required for NAREs roughly doubled. This is due to the lack of symmetry in NAREs, reflected by the 2k+2 factor in Table 1, instead of the 2k+1 factor in [26].
Example 1
We tried various parameters in τ1, τ2, mmaxand lmax, looking for a good balance between accuracy
and efficiency. We may choose large τi, and restrict mmax, lmax, to find out the relative numerical
ranks of X. Equally, we may choose small mmax, lmax, and relax τi with small values, to find
out the corresponding achievable accuracy. For the numerical results we present here, we try to achieve accuracies as high as possible, within a few thousand seconds. For higher accuracies, smaller τi and larger mmax, lmax are required, demanding more computing resources. In the
following examples, the choice τi = 10−12 or 10−15 makes very little different in the numerical
results.
k kδHkk/kHkk erk mk δtk tk
1 1.7e−03 1.6e−01 2 1.50e−01 1.50e−01
2 3.5e−01 6.0e−02 4 1.90e−01 3.40e−01
3 1.6e−01 1.1e−02 8 2.90e−01 6.30e−01
4 3.4e−02 5.6e−04 16 4.30e−01 1.06e+00
5 1.9e−03 2.2e−06 32 1.34e+00 2.40e+00
6 7.7e−06 5.3e−11 59 4.32e+00 6.72e+00
7 1.9e−10 2.5e−15 104
For n = 10, 000, a near machine accuracy of 10−15 was achieved within 6.72 seconds and 7 iterations, as summarized in Table 2. The CPU time tk roughly double that for the previous
iteration, as predicted in Table 1.
Example 2
k kδHkk/kHkk erk mk δtk tk
1 5.8e−01 2.1e−01 2 3.50e−01 3.50e−01
2 4.2e−01 9.0e−02 4 2.70e−01 6.20e−01
3 2.2e−01 2.2e−02 8 6.50e−01 1.27e+00
4 6.3e−02 1.8e−03 16 2.10e+00 3.37e+00
5 5.8e−03 1.8e−05 29 6.43e+00 9.80e+00
6 5.9e−05 2.5e−09 53 2.41e+01 3.39e+01
7 8.5e−09 5.9e−15 88
Table 3: Example 2 (n = 50, 000; γ = 196.25; τ1, τ2= 10−15; mk, lk ≤ 500)
For n = 50, 000, a near machine accuracy of 10−15 was achieved within 33.9 seconds and 7 iterations, as summarized in Table 3. The CPU time tk roughly double that for the previous
iteration for the first few iterations, as predicted in Table 1. For k ≥ 4, tk ≈ 3tk−1, possibly
because of the swapping of virtual memory between the RAM and the hard disk, because of the large amount of data involved. This observation holds true also for the later larger examples.
Example 3
k kδHkk/kHkk erk mk δtk tk
1 2.4e+00 2.0e−01 2 5.30e−01 5.30e−01
2 4.0e−01 7.8e−02 4 4.70e−01 1.00e+00
3 1.9e−01 1.7e−02 8 1.37e+00 2.37e+00
4 5.1e−02 1.2e−03 16 4.51e+00 6.88e+00
5 3.9e−03 8.2e−06 32 1.32e+01 2.01e+01
6 2.8e−05 5.6e−10 59 5.49e+01 7.50e+01
7 2.0e−09 1.2e−14 104 1.73e+02 2.48e+02
8 2.7e−14 9.2e−15 165
Table 4: Example 3 (n = 100, 000; γ = 136.05; τ1, τ2= 10−15; mk, lk≤ 500)
For n = 100, 000, a near machine accuracy of 10−15 was achieved within 248 seconds and 8 iterations, as summarized in Table 4. The CPU time tk roughly double that for the previous
iteration for the first few iterations. For k ≥ 4, tk ≈ 3tk−1 because of the large amount of data
k kδHkk/kHkk erk mk δtk tk
1 7.9e+01 3.0e−01 2 1.33e+00 1.33e+00
2 5.5e−01 1.6e−01 4 1.63e+00 2.96e+00
3 3.4e−01 5.5e−02 8 7.92e+00 1.09e+01
4 1.4e−01 8.7e−03 16 2.22e+01 3.30e+01
5 2.5e−02 2.8e−04 32 6.98e+01 1.03e+02
6 8.4e−04 3.3e−07 59 2.42e+02 3.45e+02
7 1.0e−06 5.8e−13 104 7.32e+02 1.08e+03
8 1.8e−12 3.0e−14 165
Table 5: Example 4 (n = 500, 000; γ = 59.58; τ1, τ2= 10−15; mk, lk ≤ 500)
Example 4
For n = 500, 000, a near machine accuracy of 10−14 was achieved within 1,080 seconds and 8 iterations, as summarized in Table 5. The CPU time tk roughly double that for the previous
iteration for the first few itera tions. For k ≥ 3, tk ≈ 3tk−1 because of the large amount of data
involved.
Example 5
k kδHkk/kHkk erk mk δtk tk
1 1.6e+02 3.0e−01 2 2.16e+00 2.16e+00
2 5.5e−01 1.6e−01 4 2.67e+00 4.83e+00
3 3.4e−01 5.5e−02 8 1.40e+01 1.88e+01
4 1.4e−01 8.7e−03 16 3.76e+01 5.64e+01
5 2.5e−02 2.8e−04 32 1.17e+02 1.73e+02
6 8.4e−04 3.3e−07 59 4.15e+02 5.88e+02
7 1.0e−06 6.0e−13 104 1.32e+03 1.91e+03
8 1.8e−12 2.3e−13 165
Table 6: Example 5 (n = 1, 000, 000; γ = 59.59; τ1, τ2= 10−15; mk, lk ≤ 300)
For n = 1, 000, 000, an accuracy of 10−13 was achieved within 1,910 seconds and 8 iterations, as summarized in Table 6. The CPU time tk roughly double that for the previous iteration for
the first few iterations. For k ≥ 3, tk ≈ 3tk−1because of the large amount of data involved. We
accept the slightly lower accuracy in Examples 4 and 5 in order to avoid excessive amount of computation. One or two more iterations will produce better answers but approximately two or four times more computing time will be required.
5
Conclusions
We have proposed a structure-preserving doubling algorithm, the SDA ls, for the large-scale nonsymmetric algebraic Riccati equation in (1), with A and D being large and sparse(-like), and
B and C being low-ranked. The trick is to apply the Sherman-Morrison-Woodbury formula when appropriate and not to form Ekand Fk(the iterates for E and F ) explicitly. For well-behaved
non-critical NAREs, low-ranked approximations to the solutions X and Y can be obtained efficiently. The convergence of the SDA ls is quadratic, ignoring the compression and truncation of Bk and
Ck, as shown in [17, 28]. The computational complexity and memory requirement are both O(n)
per iteration, provided that the growth of Bk and Ck is controlled.
Similar to the methods in [25, 26, 27] for other large-scale algebraic Riccati and linear matrix equations, our technique can be applied when A and D are large and sparse(-like), or are products (inverses) of such matrices. The feasibility of the SDA ls depends on whether A−1u, A−>u, D−1v
and D−>v can be formed efficiently, for arbitrary vectors u and v.
The proposed SDA ls solves NAREs efficiently. For instance in Section 4, Example 5, of dimension n = 1, 000, 000 with one trillion variables in the solution X, was solved using MATLAB on a HP SL390s workstation within 1, 910 seconds to a near-machine accuracy of O(10−13).
The adaptation of the SDA to other large-scale NAREs is under investigation and will be reported elsewhere.
Acknowledgement
This research work is partially supported by the National Science Council and the National Centre for Theoretical Sciences, as well as the TaiDa Institute of Mathematical Sciences and the Center of Advanced Study for Theoretical Sciences (at the National Taiwan University), Taiwan. The first author has been supported by a Monash Graduate Scholarship and a Monash International Postgraduate Research Scholarship. Part of the work was completed when the third author visited the National Centre for Theoretical Sciences at Tainan and the National Chiao Tung University, and we would like to thank these institutions.
References
[1] z.-z. bai, y.-h. gao and l.-z. lu. Fast iterative schemes for nonsymmetric algebraic Riccati equations arising from transport theory, SIAM J. Sci. Comput., 30 (2008) 804–818.
[2] a. berman and r.j. plemmons. Nonnegative Matrices in the Mathematical Sciences, SIAM, Philadelphia, 1994.
[3] d.s. bernstein and w.m. haddad. LQG control with an H∞performance bound: a Riccati
equation approach, IEEE Trans. Automat. Control, 34 (1989) 293–305.
[4] d.a. bini, b. iammazzo and f. poloni. A fast Newtons method for a nonsymmetric alge-braic Riccati equation, SIAM J. Matrix Anal. Appl., 30 (2008) 276–290.
[5] c.-y. chiang, e.k.-w. chu, c.-h. guo, t.-m. huang, w.-w. lin and s.-f. xu. Conver-gence analysis of the doubling algorithm for several nonlinear matrix equations in the critical case, SIAM J. Matrx Anal. Applic., 31 (2009) 227–247.
[6] e.k.-w. chu, h.-y. fan and w.-w. lin. A structure-preserving doubling algorithm for continuous-time algebraic Riccati equations, Linear Alg. Appl., 396 (2005) 55–80.
[7] e.k.-w. chu, h.-y. fan, w.-w. lin and c.-s. wang. A structure-preserving doubling al-gorithm for periodic discrete-time algebraic Riccati equations, Int. J. Control, 77 (2004) 767–788.
[8] e.k.-w. chu, t.m. huang, w.-w. lin and c.-t. wu. Palindromic eigenvalue problems: a brief survey, Taiwanese J. Math., (to appear).
[9] b. de moor and j. david. Total linear least squares and the algebraic Riccati equations, Systems Control Lett., 18 (1992) 329–337.
[10] i. gohberg and m.a. kaashoek. An inverse spectral problem for rational matrix functions and minimal divisibility, Integral Equations Operator Theory, 10 (1987) 437–465.
[11] g.h. golub and c.f. van loan. Matrix Computations, 2nd Ed., Johns Hopkins University Press, Baltimore, MD, 1989.
[12] c.-h. guo. Nonsymmetric algebraic Riccati equations and Wiener-Hopf factorization for M-matrices, SIAM J. Matrix Anal. Applic., 23 (2001) 225–242.
[13] c.-h. guo. A new class of nonsymmetric algebraic Riccati equations, Lin. Alg. Applic., 426 (2007) 636–649.
[14] c.-h. guo and n.j. higham. Iterative solution of a nonsymmetric algebraic Riccati equation, SIAM J. Matrix Anal. Appl., 29 (2007) 396–412.
[15] c.-h. guo and p. lancaster. Analysis and modification of Newton’s method for algebraic Riccati equations, Math. Comp., 67 (1998) 1089–1105.
[16] c.-h. guo and w.-w. lin. Convergence rates of some iterative methods for nonsymmetric algebraic Riccati equations arising in transport theory, Lin. Alg. Appl., 432 (2010) 283–291.
[17] x.-x. guo, w.-w. lin and s.-f. xu. A structure-preserving doubling algorithm for nonsym-metric algebraic Riccati equations, Numer. Math., 103 (2006) 393–412.
[18] d. hinrichsen, b. kelb and a. linnemann. An algorithm for the computation of the structured complex stability radius, Automatica, 25 (1989) 771–775.
[19] j. juang. Existence of algebraic matrix Riccati equations arising in transport theory, Lin. Alg. Applic., 230 (1995) 89–100.
[20] k. jbilou. Low rank approximate solutions to large Sylvester matrix equations, Appl. Math. Comp., 177 (2006) 365–376.
[21] p. lancaster and l. rodman. Solutions of the continuous and discrete-time algebraic Riccati equations: a review, in The Riccati Equations, S. Bittanti, A.J. Lamb and J.C. Willems eds., Springer-Verlag, Berlin, 1991.
[22] p. lancaster and l. rodman. Algebraic Riccati Equations, Clarendon Press, Oxford, 1995. [23] t. li, e.k.-w. chu, j. juang and w.-w. lin. Solution of a nonsymmetric algebraic Riccati equation from a two-dimensional transport model, Lin. Alg. Applic., 434 (2011) 201–214.
[24] t. li, e.k.-w. chu, j. juang and w.-w. lin. Solution of a nonsymmetric algebraic Riccati equation from a one-dimensional multi-state transport model, IMA J. Numer. Anal., 31 (2011) 1453–1467.
[25] t. li, e.k.-w. chu and w.-w. lin. Solving large-scale discrete-time algebraic Riccati equa-tions by doubling, Technical Report, NCTS, National Chiao Tung University, Hsinchu, Tai-wan, 2011.
[26] t. li, e.k.-w. chu, w.-w. lin and c.-y. weng. Solving large-scale continuous-time alge-braic Riccati equations by doubling, Technical Report, NCTS, National Chiao Tung Univer-sity, Hsinchu, Taiwan, 2011.
[27] t. li, c.-y. weng, e.k.-w. chu and w.-w. lin. Solving large-scale Stein and Lyapunov equations by doubling, Technical Report, NCTS, National Chiao Tung University, Hsinchu, Taiwan, 2011.
[28] w.-w. lin and s.-f. xu. Convergence analysis of structure-preserving doubling algorithms for Riccati-type matrix equations. SIAM J. Matrix Anal. Appl., 28 (2008) 26–39.
[29] y. lin, l. bao and y. wei. A modified Newton method for solving non-symmetric algebraic Riccati equations arising in transport theory. IMA J. Numer. Anal., 29 (2008) 215–224.
[30] l.-z. lu. Newton iterations for a non-symmetric algebraic Riccati equation, Numer. Lin. Alg. Applic., 12 (2005) 191–200.
[31] mathworks. MATLAB User’s Guide, 2010.
[32] v. mehrmann. The Autonomous Linear Quadratic Control Problem, Lecture Notes in Con-trol and Information Sciences 163, Springer-Verlag, 1991.
[33] v. mehrmann and h. xu. Explicit solutions for a Riccati equations from transport theory, SIAM J. Matrix Anal., 30 (2008) 1339–1357.
[34] r.a. smith. Matrix equation XA + BX = C, SIAM J. Appl. Math., 16 (1968) 198–201. [35] n. truhar and r.-c. li. On ADI method for Sylvester equations, Technical Report 2008-02,
Department of Mathematics, University of Texas at Arlington, 2008.
[36] d. williams. A potential-theoretical note on the quadratic Wiener-Hopf equation for Q-matrices, in Seminar on Probability XVI, Lecture Notes in Mathematics 920, pp. 91–94, Springer-Verlag, Berlin, 1982.
