Structured Doubling Algorithm for a Class of Large-Scale Discrete-Time Algebraic Riccati Equations with High-Ranked Constant Term

Abstract

Consider the computation of the solution for a class of discrete-time algebraic Riccati equations (DAREs) with the low-ranked coefficient matrix G and the high-ranked constant matrix H. A structured doubling algorithm is proposed for large-scale problems when A is of lowrank. Compared to the existing doubling algorithm of $O(2^{k}n)$ flops at the k-th iteration, the newly developed version merely needs $O(n)$ flops for preprocessing and $O({(k+1)}^3{m}^3)$ flopsfor iterations and is more proper for large-scale computations when $m\ll n$. The convergence and complexity of the algorithm are subsequently analyzed. Illustrative numerical experiments indicate that the presented algorithm, which consists of a dominant time-consuming preprocessing step and a trivially iterative step, is capable of computing the solution efficiently for large-scale DAREs.

1. Introduction

Consider a discrete-time control system

$$x_{k+1}=A{x}_k+B{u}_k,k=0,1,2,...,$$

where $A\in {\mathbb{C}}^{n\times n}$ and $B\in {\mathbb{C}}^{n\times l}$ with $l\le n$. Here, ${\mathbb{C}}^{n\times m}$ stands for sets of $n\times m$ complex matrices. The linear quadratic regulator (LQR) control minimizes the energy or the cost functional

$$J_c(x_k,u_k)\equiv \sum _{k=0}^{\infty }\left[x_k^{*}H{x}_k+u_k^{*}R{u}_k\right]$$

with the Hermitian constant term $H\in {\mathbb{C}}^{n\times n}$ being positive semi-definite [1]. Here, the symbol “*” is the conjugate transpose of a vector or a matrix.

The corresponding optimal control is

$$u_k=-F{x}_k$$

and the feedback gain matrix

$$F:={(R+B^{*}XB)}^{-1}(B^{*}XA)$$

can then be expressed in terms of the unique positive semi-definite stabilizing solution X of the discrete-time algebraic Riccati equation (DARE) [2]

$$\mathcal{D}(X)=-X+A^{*}X{(I+GX)}^{-1}A+H=0,$$

(1)

where $G=B{R}^{-1}{B}^{*}$ with $R\in {\mathbb{C}}^{l\times l}$, $H\in {\mathbb{C}}^{n\times n}$ is Hermitian and positive semi-definite. In many control problems, the matrix $A\in {\mathbb{C}}^{n\times n}$ is sparse in the sense that the matrix-vector product $Av$ and the inverse-vector product $A^{-1}v$ require $O(n)$ flops, respectively. The recent applications of the discrete-time control system can be found in [3] such as the wheeled robot and the airborne pursuer. There are also some applications (e.g., the singular Kalman filter) about the fractional Riccati equation, see [4,5] and the references therein.

The existence of the unique positive semi-definite solution X of DARE (1) has been well studied if $(A,B)$ is d-stabilizable and $(H,A)$ is observable, see [6,7] and their references for more details. The structure-preserving doubling algorithm (SDA) is one of the most efficient methods [7] to compute the unique positive semidefinite solution X via the following iteration

$$\left\{\begin{array}{c}{A}_{k+1}=A_k{(I-G_k{H}_k)}^{-1}{A}_k,\hfill \\ H_{k+1}=H_k+A_k^{*}{H}_k{(I-G_k{H}_k)}^{-1}{A}_k,\hfill \\ G_{k+1}=G_k+A_k{(I-G_k{H}_k)}^{-1}{G}_k{A}_k^{*}\hfill \end{array}\right.$$

(2)

with $A_0=A$, $G_0=-G$, $H_0=H$. Regardless of the structure of coefficient matrices, the computational complexity of each iteration is about $O(n^3)$, obviously not fitting for large-scale problems. When the constant matrix H is low-ranked, the solution X is commonly numerically low-ranked and can be approximated by $H_k$ in terms of a series of decomposed matrix factors, making the SDA feasible for large-scale DAREs [8]. If only the feedback gain matrix F is required without outputting the solution X, an adaptive version of the SDA in [9] still works for large-scale problems even if H is high-ranked. In that case, the solution X is no longer numerically low-ranked but can be stored in a sequence of matrix-vector products [9]. In both situations, the computational complexity of the SDA at the kth iteration costs about $O(2^{k}n)$ flops (i.e., the exponential increase in k), resulting in the intolerable iteration time when k is large.

In this paper, we consider DAREs with A of the low-ranked structure (which may not be sparse)

$$A=C_{1}S{C}_2^{*}$$

(3)

with $C_1,C_2\in {\mathbb{C}}^{n\times m}$ and $S\in {\mathbb{C}}^{m\times m}$ ($m\ll n$). The motivation behind this is that the complexity of the SDA at the k-th iteration might be further reduced in this case and the DAREs, with the structure (3), have several applications in circuit-controlling areas, for example, the circuits system with $C_1$ and $C_2$ being the mesh inductance matrices, composed of the product of several mesh matrices (n is the number of meshes) and S being the resistance matrix [10]. To obtain the optimal feedback gain to control the circuit system, one is required to find the solution of the DARE (1).

The main contribution we made under the low-ranked structure (3) is that the computational complexity of the SDA at the k-th iteration can be further reduced to $O({(k+1)}^3{m}^3)$, far less than $O(2^{k}n)$ when $m\ll n$. As a result, the most time-consuming part of the SDA lies in the preprocessing step with a fixed computational complexity $O(n)$, and the other part for the iterations might be accordingly insignificant. Numerical experiments are implemented to validate the effectiveness of the presented algorithm, constituting a useful complement to the solver for computing the solution of DAREs.

The rest of the paper is organized as follows. In Section 2, we develop the structured SDA for DAREs with a low-ranked structure of A and construct its convergence. A detailed complexity analysis as well as the design of the termination criterion are established in Section 3. Section 4 is devoted to numerical experiments to indicate the efficiency of the proposed algorithm, and the conclusion is drawn in the last section.

Notation. Symbols ${\mathbb{R}}^{n\times n}$ and ${\mathbb{C}}^{n\times n}$ in this paper stand for sets of $n\times n$ real and complex matrices, respectively. $I_n$ is the $n\times n$ identity matrix. For a matrix $A\in {\mathbb{C}}^{n\times n}$, $\sigma (A)$ and $\rho (A)$ denote, respectively, the spectrum and spectral radius of A. A Hermitian matrix $A>0$ ($\ge 0$) when all its eigenvalues are positive (non-negative). Additionally, $M>N$ ($M\ge N$) if and only if $M-N>0$ ($\ge 0$).

We also need the concept of the numerically low-ranked matrix.

Definition 1.

([8]) A matrix A is said to be numerically low-ranked with respect to tolerance $ϵ>0$ if $\mathrm{rank}(A)\le c_{ϵ}$ for a constant $c_{ϵ}$ associated with ϵ but independent of the size of A.

2. Structured Doubling Algorithm

In this section, we describe the structured iteration scheme for DAREs with a high-ranked constant term and low-ranked A in (3). To avoid the inversion of large-scale matrices, the Sherman–Morrison–Woodbury formula (SMWF) [11,12] is first applied to the sparse-plus-low-ranked matrices to represent the corresponding structured matrices. Then, we aim at preserving the sparsity or the low-ranked structure of the iteration sequence rather than forming it explicitly. As a result, the SDA is capable of being implemented only with some small-scale matrices, referred to as kernels, and the complexity of the iteration can be ignored more easily than that of the preprocessing step for large-scale problems.

2.1. Iteration Scheme

Given the initial matrices $A_0=C_{1}S{C}_2$, $H_0=H$, $G_0=B{R}^{-1}{B}^{*}$, $S_0=S$, $T_0=0$, $G_0=R^{-1}$, and $B_0=B$, the SDA will be organized according to the following format:

$$\left\{\begin{array}{c}{A}_k=C_1{S}_k{C}_2^{*},\hfill \\ H_k=H+C_2{T}_k{C}_2^{*},\hfill \\ G_k=B_k{R}_k{B}_k^{*}\hfill \end{array}\right.$$

(4)

for $k\ge 0$, where $S_k,T_k\in {\mathbb{C}}^{m\times m}$, $B_k\in {\mathbb{C}}^{n\times (km+l)}$, and $R_k\in {\mathbb{C}}^{(km+l)\times (km+l)}$. One merit of the above scheme (4) is that the sizes of kernels $S_k$ and $T_k$ are always invariant (i.e., $m\times m$) during iterations. Although the column of $B_k$ and the size of $R_k$ increase linearly with respect to k, the enhanced scale is generally small due to the fast convergence of the SDA. Then, $G_k$ still hopefully maintains a low-ranked structure and could be derived and stored in an economic way.

Let ${\mathsf{\Sigma }}_k=R_k^{-1}-B_k^{*}{H}_k{B}_k$. By applying the Sherman–Morrison–Woodbury formula (SMWF) [11], we have

$${(I-G_k{H}_k)}^{-1}=I+B_k{\mathsf{\Sigma }}_k^{-1}{B}_k^{*}{H}_k,{(I-H_k{G}_k)}^{-1}=I+H_k{B}_k{\mathsf{\Sigma }}_k^{-1}{B}_k^{*}.$$

(5)

Insertion (5) into the SDA (2) with currently available $A_k$, $H_k$ and $G_k$ yield $A_{k+1}=C_1{S}_{k+1}{C}_2^{*}$, $H_{k+1}=H+C_2{T}_{k+1}{C}_2^{*}$, $G_{k+1}=B_{k+1}{R}_{k+1}{B}_{k+1}^{*}$ with

$$\left\{\begin{array}{c}{S}_{k+1}=S_k(C_2^{*}{C}_1+{\mathsf{\Phi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Psi }}_k)S_k,\hfill \\ T_{k+1}=T_k+S_k^{*}(C_1^{*}{H}_k{C}_1+{\mathsf{\Psi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Psi }}_k)S_k,\hfill \\ B_{k+1}=[C_1,B_k],\hfill \\ R_{k+1}=\left[\begin{array}{cc}{\tilde{R}}_k& \\ & R_k\end{array}\right]\hfill \end{array}\right.$$

(6)

and

$${\mathsf{\Phi }}_k=C_2^{*}{B}_k,{\mathsf{\Psi }}_k=C_1^{*}{H}_k{B}_k,{\tilde{R}}_k=S_k{\mathsf{\Phi }}_k{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Phi }}_k^{*}{S}_k^{*}.$$

The main computational task of (6) is the update of $H_k{B}_k$, $B_k^{*}{H}_k{B}_k$ in ${\mathsf{\Psi }}_k$, ${\mathsf{\Phi }}_k$, ${\mathsf{\Sigma }}_k$ and the solutions of two linear system associated with ${\mathsf{\Sigma }}_k$. Regardless of the concrete structure, the complexity of such calculations is $O(2^{k}n)$ [8,9]. A deeper observation made here will show that such computations can be further down to the complexity of $O({(k+1)}^3{m}^3)$, far less than that of the preprocessing for large-scale problems with $m\ll n$. In fact, by setting $B_0=B$, it follows from (6) that

$$B_k=[\stackrel{\stackrel{\mathrm{km}}{⏞}}{C_1,C_1,...,C_1},\stackrel{\stackrel{l}{⏞}}{B}]n$$

and thus

$${\mathsf{\Phi }}_k=[\stackrel{\stackrel{\mathrm{km}}{⏞}}{C_2^{*}{C}_1,C_2^{*}{C}_1,...,C_2^{*}{C}_1},\stackrel{\stackrel{l}{⏞}}{C_2^{*}B}]m.$$

(7)

Analogously, we have

$$\begin{array}{cc}{H}_k{B}_k\hfill & =H{B}_k+C_2{T}_k{C}_2^{*}{B}_k\hfill \\ & =[\stackrel{\stackrel{\mathrm{km}}{⏞}}{H{C}_1+C_2{T}_k{C}_2^{*}{C}_1...H{C}_1+C_2{T}_k{C}_2^{*}{C}_1}\stackrel{\stackrel{l}{⏞}}{HB+C_2{T}_k{C}_2^{*}B}]n\hfill \end{array}$$

and

$$\begin{array}{c}{\mathsf{\Psi }}_k=[\stackrel{\stackrel{\mathrm{km}}{⏞}}{C_1^{*}H{C}_1+C_1^{*}{C}_2{T}_k{C}_2^{*}{C}_1...C_1^{*}H{C}_1+C_1^{*}{C}_2{T}_k{C}_2^{*}{C}_1}\stackrel{\stackrel{l}{⏞}}{C_1^{*}HB+C_1^{*}{C}_2{T}_k{C}_2^{*}B}]m.\hfill \end{array}$$

Furthermore, as

$$\begin{array}{cccc}& B_k^{*}{H}_k{B}_k=B_k^{*}H{B}_k+B_k^{*}{C}_2{T}_k{C}_2^{*}{B}_k\hfill & & \\ =& \stackrel{km}{\overbrace{[\begin{array}{ccc}{C}_1^{*}H{C}_1+C_1^{*}{C}_2{T}_k{C}_2^{*}{C}_1& \dots & C_1^{*}H{C}_1+C_1^{*}{C}_2{T}_k{C}_2^{*}{C}_1\\ ⋮& \ddots & ⋮\\ C_1^{*}H{C}_1+C_1^{*}{C}_2{T}_k{C}_2^{*}{C}_1& \dots & C_1^{*}H{C}_1+C_1^{*}{C}_2{T}_k{C}_2^{*}{C}_1\\ B^{*}H{C}_1+B^{*}{C}_2{T}_k{C}_2^{*}{C}_1& \dots & B^{*}H{C}_1+B^{*}{C}_2{T}_k{C}_2^{*}{C}_1\end{array}}}\hfill & \stackrel{l}{\overbrace{\begin{array}{c}{C}_1^{*}HB+C_1^{*}{C}_2{T}_k{C}_2^{*}B\\ ⋮\\ C_1^{*}HB+C_1^{*}{C}_2{T}_k{C}_2^{*}B\\ B^{*}HB+B^{*}{C}_2{T}_k{C}_2^{*}B\end{array}]}}& \begin{array}{l}\left.\begin{array}{l}\\ \\ \end{array}\right\}km\\ \left.\right\}l\end{array},\end{array}$$

(8)

, the update of the matrix

$${\Sigma }_k=R_k^{-1}-B_k^{*}{H}_k{B}_k=\left[\begin{array}{cc}{\tilde{R}}_{k-1}^{-1}& \\ & R_{k-1}^{-1}\end{array}\right]-B_k^{*}{H}_k{B}_k$$

will be of size $(km+l)\times (km+l)$. Now, suppose that matrices $C_1^{*}H{C}_1$, $C_1^{*}HB$, $B^{*}HB$, $C_2^{*}{C}_1$, and $C_2^{*}B$ are available in the preprocessing step, then ${\mathsf{\Phi }}_k$ in (7) does not require additional computations. Additionally, ${\mathsf{\Psi }}_k$ and ${\mathsf{\Sigma }}_k$ can be obtained via updating several small scale matrix multiplications of the size $m\times m$, i.e., ${(C_2^{*}{C}_1)}^{*}{T}_k(C_2^{*}{C}_1)$, ${(C_2^{*}{C}_1)}^{*}{T}_k(C_2^{*}B)$, and ${(C_2^{*}B)}^{*}{T}_k(C_2^{*}B)$, and replicating them $km+l$ times (here ${\tilde{R}}_{k-1}^{-1}$ and $R_{k-1}^{-1}$ are assumed to be available in the last iteration for computing ${\mathsf{\Sigma }}_k$). Consequently, the left computation lies in solving two linear systems ${\mathsf{\Sigma }}_{k}U={\mathsf{\Phi }}_k$ and ${\mathsf{\Sigma }}_{k}V={\mathsf{\Psi }}_k$ of size $(km+l)\times (km+l)$. We summarize the whole process in Algorithm 1 as below; the concrete complexity analysis in the next section shows that the iteration only costs about $O({(k+1)}^3{m}^3)$ flops.

Remark 1.

The output matrices $B_{ϵ}$ and $R_{ϵ}$ are numerically low-ranked with respect to the tolerance ϵ. $\hat{T}$ is the matrix from the convergence of $T_k$ given in the next subsection.

Remark 2.

The QR decomposition of $C_2$ is for the derivation of the relative residual and also could be implemented in the preprocessing step. The computational complexity of the preprocessing part is about $O(n)$ flops, taking the dominant CPU time compared with the iteration part.

Remark 3.

The computations of the iteration part and of the relative residual in the DARE cost about $O({(k+1)}^3{m}^3)$ and $O(m^3)$ flops, respectively, much less than $O(n)$ of the preprocessing part when $m\ll n$. Hence, the main computation of Algorithm 1 concentrates on the preprocessing part.

2.2. Convergence

To establish the convergence of Algorithm 1, we first review some results for iteration format (2).

Algorithm 1. Structured SDA for DAREs.
Input:	$C_1$, $C_2$, S, B, $R^{-1}=R^{-*}$, H and tolerances ${\tau }_g$ and $ϵ$, and $m_{\max }$;
Output:	$B_{ϵ}\in {\mathbb{C}}^{n\times m_{ϵ}}$, $R_{ϵ}\in {\mathbb{C}}^{m_{ϵ}\times m_{ϵ}}$, $T_{ϵ}\in {\mathbb{C}}^{m_{ϵ}\times m_{ϵ}}$, and normalized relative resi-
	dual ${\tilde{r}}_{ϵ}$;
Preprocess:	Compute $C_1^{*}H{C}_1$, $C_1^{*}HB$, $B^{*}HB$, $C_2^{*}{C}_1$, $C_2^{*}B$ and the economic QR de-
	composition of $C_2$.
Iteration:	Set $T_0=0$, $S_0=S$, $R_0=-R^{-1}$, $B_0=B$, $H_0=H$, ${\mathsf{\Sigma }}_0=-(R+B^{*}HB)$,
	${\mathsf{\Phi }}_0=C_2^{*}B$, ${\mathsf{\Psi }}_0=C_1^{*}HB$ and $k=0$;
	For $k\ge 1$, do until convergence:
	Compute the relative residual ${\tilde{r}}_k$ as in (11).
	If ${\tilde{r}}_k\le ϵ$, set $B_{ϵ}=B_k$, $R_{ϵ}=R_k$, $T_{ϵ}=T_k$ and ${\tilde{r}}_{ϵ}={\tilde{r}}_k$; Exit;
	End If
	Compute
	$S_{k+1}=S_k(C_2^{*}{C}_1+{\mathsf{\Phi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Psi }}_k)S_k$;
	$T_{k+1}=T_k+S_k^{*}(C_1^{*}{H}_k{C}_1+{\mathsf{\Psi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Psi }}_k)S_k$;
	$R_{k+1}=\left[\begin{array}{cc}{S}_k{\mathsf{\Phi }}_k{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Phi }}_k^{*}{S}_k^{*}& \\ & R_k\end{array}\right]$;
	Obtain $B_{k+1}^{*}{H}_{k+1}{B}_{k+1}$ in (8) with preprocessed matrices.
	${\mathsf{\Sigma }}_{k+1}^{-1}={(I-R_{k+1}{B}_{k+1}^{*}{H}_{k+1}{B}_{k+1})}^{-1}{R}_{k+1}$,
	${\mathsf{\Phi }}_{k+1}=[C_2^{*}{C}_1,{\mathsf{\Phi }}_k]$,
	${\mathsf{\Psi }}_{k+1}=[C_1^{*}H{C}_1+C_1^{*}{C}_2{T}_{k+1}{C}_2^{*}{C}_1,...,C_1^{*}H{C}_1+C_1^{*}{C}_2{T}_{k+1}{C}_2^{*}{C}_1$,
	$C_1^{*}HB+C_1^{*}{C}_2{T}_{k+1}{C}_2^{*}B];$
	Set $k\leftarrow k+1$.
	End Do

Theorem 1.

([13]) Assume that X and Y are the Hermitian and positive semi-definite solutions of the DARE (1) and its dual equation

$${\mathcal{D}}_d(Y)=-Y+AY{(I+HY)}^{-1}{A}^{*}+G=0,$$

(9)

respectively. Let $P:={(I+GX)}^{-1}A$ and $Q:={(I+HY)}^{-1}{A}^{*}$. Then, the matrix sequences $\left\{A_k\right\}$, $\left\{G_k\right\}$ and $\left\{H_k\right\}$ generated by the SDA (2) satisfy

$$\begin{array}{c}(1)A_k=(I+G_{k}X)P^{2^k};\hfill \\ (2)H\le H_k\le H_{k+1}\le X,X-H_k={\left(P^{*}\right)}^{2^k}(X+X{G}_{k}X)P^{2^k};\hfill \\ (3)G\le G_k\le G_{k+1}\le Y,Y-G_k={\left(Q^{*}\right)}^{2^k}(Y+Y{H}_{k}Y)Q^{2^k}.\hfill \end{array}$$

(10)

It follows from (10) that

$$\begin{array}{c}\parallel A_k\parallel \le (1+\parallel X\parallel ·\parallel Y\parallel )\parallel P^{2^k}\parallel ,\hfill \\ \parallel H_k-X\parallel \le \parallel X\parallel (1+\parallel X\parallel ·\parallel Y\parallel )\parallel P^{2^k}{\parallel }^2,\hfill \\ \parallel G_k-Y\parallel \le \parallel Y\parallel (1+\parallel X\parallel ·\parallel Y\parallel )\parallel Q^{2^k}{\parallel }^2.\hfill \end{array}$$

indicating that sequences $\left\{A_k\right\}$, $\left\{H_k\right\}$ and $\left\{G_k\right\}$ converge quadratically to zero, X, and Y, respectively, if $\rho (P)<1$ and $\rho (Q)<1$. By noting the decomposition $A_k=C_1{S}_k{C}_2^{*}$, the sequence $\left\{S_k\right\}$ must converge to zero. On the other hand, the decomposition $H_k=H+C_2{T}_k{C}_2^{*}$ implies that the sequence $\left\{T_k\right\}$ converges to some matrix $\hat{T}\in {\mathbb{C}}^{m\times m}$ such that the solution of the DARE $X=H+C_2\hat{T}{C}_2^{*}$. At last, the decomposition

$$\begin{array}{ccc}\hfill G_k& =& B_k{R}_k{B}_k^{*}\hfill \\ & =& [C_1,...,C_1,B]·\left[\begin{array}{cccc}{S}_k{\mathsf{\Phi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Phi }}_k{S}_k^{*}& & & \\ & \ddots & & \\ & & S_1{\mathsf{\Phi }}_1^{*}{\mathsf{\Sigma }}_1^{-1}{\mathsf{\Phi }}_1{S}_1^{*}& \\ & & & S_0{\mathsf{\Phi }}_0^{*}{\mathsf{\Sigma }}_0^{-1}{\mathsf{\Phi }}_0{S}_0^{*}\end{array}\right]·\left[\begin{array}{c}{C}_1^{*}\\ ⋮\\ C_1^{*}\\ B^{*}\end{array}\right]\hfill \\ & =& B{S}_0{\mathsf{\Phi }}_0^{*}{\mathsf{\Sigma }}_0^{-1}{\mathsf{\Phi }}_0{S}_0^{*}{B}^{*}+\sum _{i=1}^k{C}_1{S}_k{\mathsf{\Phi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Phi }}_k{S}_k^{*}{C}_1^{*}\hfill \end{array}$$

indicates that the solution Y of the dual DARE has a numerically low-ranked decomposition $Y\approx B_{ϵ}{R}_{ϵ}{B}_{ϵ}^{*}$ with respect to a sufficient small tolerance $ϵ>0$. So, we have the following corollary.

Corollary 1.

Suppose that X and Y are the Hermitian and positive semi-definite solutions of the DARE (1) and its dual form (9), respectively. Then, for Algorithm 1, the sequence $\left\{S_k\right\}$ converges to zero matrix quadratically, and $\left\{T_k\right\}$ converges to some matrix $\hat{T}$ with $X=H+C_2\hat{T}{C}_2^{*}$. Moreover, for sufficiently large k, the matrix $R_k$ is numerically low-ranked with respect to tolerance ϵ. That is, the solution Y of the dual Equation (9) has the low-ranked approximation $Y\approx B_{ϵ}{R}_{ϵ}{B}_{ϵ}^{*}$, where matrices $B_{ϵ}$ and $R_{ϵ}$ associate with ϵ but independently of the size of Y.

3. Computational Issues

3.1. Residual and Stop Criterion

Recalling the low-ranked structures of G and A, the residual of the DARE is

$$\begin{array}{ccc}& & -H_k+A^{*}{H}_k{(I+G{H}_k)}^{-1}A+H\hfill \\ & =& C_2(-T_k+S^{*}{C}_1^{*}{H}_k{(I+B{R}^{-1}{B}^{*}{H}_k)}^{-1}{C}_{1}S)C_2^{*}\hfill \\ & =& C_2(-T_k+S^{*}({\mathsf{\Pi }}_k-{\mathsf{\Xi }}_k{\mathsf{\Theta }}_k^{-1}{\mathsf{\Xi }}_k^{*})S)C_2^{*}\hfill \end{array}$$

with

$$\begin{array}{ccc}\hfill {\mathsf{\Pi }}_k& =& C_1^{*}{H}_k{C}_1=C_1^{*}H{C}_1+{(C_2^{*}{C}_1)}^{*}·T_k·C_2^{*}{C}_1,\hfill \\ \hfill {\mathsf{\Xi }}_k& =& C_1^{*}{H}_{k}B=C_1^{*}HB+{(C_2^{*}{C}_1)}^{*}·T_k·C_2^{*}B,\hfill \\ \hfill {\mathsf{\Theta }}_k& =& R+B^{*}{H}_{k}B=R+B^{*}HB+{(C_2^{*}B)}^{*}·T_k·C_2^{*}B.\hfill \end{array}$$

Let $C_2=Q_{{}_{C_2}}{R}_{{}_{C_2}}$ ($Q_{{}_{C_2}}\in {\mathbb{C}}^{n\times m}$, $R_{{}_{C_2}}\in {\mathbb{C}}^{m\times m}$) be the economic QR decomposition of $C_2$, derived from the preprocessing step. The matrix norm of the residual is

$$r_k=\parallel R_{{}_{C_2}}(-T_k+S^{*}({\mathsf{\Pi }}_k-{\mathsf{\Xi }}_k{\mathsf{\Theta }}_k^{-1}{\mathsf{\Xi }}_k^{*})S)R_{{}_{C_2}}^{*}\parallel$$

and Algorithm 1 can be terminated by the normalized relative residual

$$\mathrm{NRRes}=r_k/(t_k+s_k+m_k):={\tilde{r}}_k<ϵ$$

(11)

with

$$r_k=\parallel R_{{}_{C_2}}{T}_k{R}_{{}_{C_2}}^{*}\parallel ,s_k=\parallel R_{{}_{C_2}}{S}_k^{*}{\mathsf{\Pi }}_k{S}_k{R}_{{}_{C_2}}^{*}\parallel ,m_k=\parallel R_{{}_{C_2}}{S}_k^{*}{\mathsf{\Xi }}_k{\mathsf{\Theta }}_k^{-1}{\mathsf{\Xi }}_k^{*}{S}_k{R}_{{}_{C_2}}^{*}\parallel .$$

Note that the calculation of NRRes only associates with several matrix operations with the small-scale $m\times m$, requiring $O(m^3)$ flops and far less than $O(n)$ when $m\ll n$.

3.2. Complexity Analysis

The main flops of Algorithm 1 come from the preprocessing step of forming matrices $C_1^{*}H{C}_1$, $C_1^{*}HB$, $B^{*}HB$, $C_2^{*}{C}_1$, $C_2^{*}B$ and QR decomposing $C_2=Q_{{}_{C_2}}{R}_{{}_{C_2}}$ with the Householder transformation in [14,15].

Table 1. Complexity and memory of the preprocessing step in Algorithm 1.

Items	Flops	Memory
$H{C}_1$, $HB$	$2m(m+l)n$	$(m+l)n$
$C_1^{*}H{C}_1$, $C_1^{*}HB$, $B^{*}HB$	$2(m^2+ml+l^2)n$	$m^2+ml+l^2$
$C_2^{*}{C}_1$, $C_2^{*}B$	$2m(m+l)n$	$m^2+ml$
$C_2=Q_{{}_{C_2}}{R}_{{}_{C_2}}$	$(10m^2+4ml+4l^2)n$	$m^2$
Total	$\begin{array}{c}(16m^2+10ml+6l^2)n\end{array}$	$\begin{array}{c}3m^2+2ml+l^2+(m+l)n\end{array}$

lists the details, where only the matrix $R_{{}_{C_2}}\in {\mathbb{C}}^{m\times m}$ stored as $Q_{{}_{C_2}}\in {\mathbb{C}}^{n\times m}$ is orthornormal satisfying $Q_{{}_{C_2}}^{*}{Q}_{{}_{C_2}}=I_m$.

It is seen from Table 1 that the computation and the storage are both of $O(n)$ flops when $m,l\ll n$. We subsequently analyze the complexity of the iteration part. Assume that the LU decomposition is employed for solving the linear system $MZ=N$ with $M,Z,N\in {\mathbb{C}}^{(km+l)\times (km+l)}$. The flops and memory of the kth iteration are summarized in

Table 2. Complexity and memory at kth iteration in Algorithm 1.

Items	Flops	Memory
${\mathsf{\Sigma }}_k^{-1}{\mathsf{\Phi }}_k^{*}$, ${\mathsf{\Sigma }}_k^{-1}{\mathsf{\Psi }}_k^{*}$	$16{(km+l)}^{2}m$	$2(km+l)m$
${\mathsf{\Phi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Psi }}_k^{*}$, ${\mathsf{\Phi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Phi }}_k^{*}$, ${\mathsf{\Psi }}_k^{*}{\mathsf{\Sigma }}_k^{-1}{\mathsf{\Psi }}_k^{*}$	$6m^2(km+l)$	$3m^2$
$S_{k+1}$	$4m^3$	$m^2$
$T_{k+1}$	$8m^3$	$m^2$
$R_{k+1}$	$4m^3$	$m^2$
$B_{k+1}^{*}{H}_{k+1}{B}_{k+1}$	$2m(2m^2+2ml+l^2)$	${((k+1)m+l)}^2$
${\mathsf{\Sigma }}_{k+1}$	$2{((k+1)m)}^3$	${((k+1)m)}^2$
${\mathsf{\Phi }}_{k+1}$	—	$((k+1)m+l)m$
${\mathsf{\Psi }}_{k+1}$	$4m^2(m+l)$	$((k+1)m+l)m$
Total	$\begin{array}{c}(24+2{(k+1)}^3)m^3\\ +2ml(4m+l)\\ +2m(km+l)\\ (8(km+l)+3m)\end{array}$	$\begin{array}{c}(4{(k+1)}^2+2k+6)m^2\\ +(2k+6)ml\\ +l^2\end{array}$

below.

Table 2 shows that the complexity of the kth iteration in Algorithm 1 is about $O({(k+1)}^3{m}^3)$, far less than $O(n)$ of the preprocessing step when $m\ll n$. Thus, the dominantly calculating cost of Algorithm 1 locates at the preprocessing step; however, it is still far less than the exponentially increasing complexity $O(2^{k}n)$ [8,9] when k grows large.

4. Numerical Experiments

In this section, we will show the effectiveness of Algorithm 1 to calculate the solution X of the large-scale DARE (1). The code was programmed by Matlab 2014a [16], and all computations were implemented in a ThinkPad notebook with 2.4 GHz Intel i5-6200 CPU and 8G memory. The stop criterion is the NRRes in (11) with a proper tolerance $ϵ$. To show the location of the dominant computations in Algorithm 1, we record the ratio of iteration time and total time in the percentage

$${\mathrm{R}}_t=\frac{\mathrm{TIME}-\mathrm{I}}{\mathrm{TIME}-\mathrm{P}+\mathrm{TIME}-\mathrm{I}}\times 100\%,$$

(12)

where “TIME-P” represents the pre-processing time elapsed for forming matrices associated with n, and “TIME-I” stands for the costed CPU time for iterations.

Example 1.

The first example is devised to measure the actual error between the true solution X and the approximated solution $H_k$ computed from Algorithm 1. Let $S=1$, $C_1=\mathbf{1}/\parallel \mathbf{1}\parallel \in {\mathbb{R}}^{n\times 1}$ and $C_2\in {\mathbb{R}}^{n\times 1}$ be a vector such that $C_1^{*}{C}_2=0$ and $C_i^{*}{C}_i=1$ ($i=1,2$), where $\mathbf{1}$ is a vector with all elements 1. Set $B^{*}=[0,0,...,0,1]\in {\mathbb{R}}^{1\times n}$, $R=1$ and $H=I_n$. Then, the solution of the DARE is

$$X=I_n-U{U}^{*}$$

with

$$U=w{C}_2$$

and

$$w^2=\frac{C_{2_n}^2-2+\sqrt{{(C_{2_n}^2-2)}^2+4C_{2_n}^2(2-C_{1_n}^2)}}{2C_{2_n}^2}$$

being the root of the equation $(1-w^2)(2+w^2*C_{2_n}^2)-C_{1_n}^2=0$. Here, $C_{1_n}$ and $C_{2_n}$ represent the n-th element of $C_1$ and $C_2$, respectively. The coefficient matrices are $A=C_{1}S{C}_2^{*}$ and $G=B{R}^{-1}{B}^{*}$. The principle of selecting the above vectors and matrices is for the convenient construction of the true solution of the DARE. Then, we can evaluate the error between the computed approximated solution and the true solution.

We consider the medium scales with $n=1000$, 3000, and 5000 to test the accuracy of Algorithm 1, which is terminated when the NRRes is less the prescribed $ϵ=1.0\times {10}^{-13}$. Numerical experiments show that Algorithm 1 always takes three iterations to obtain the approximate solution for all tested dimensions n. The obtained results on NRRes and Errors are listed in

Table 3. Residual and actual errors in Example 1.

n	1000	3000	5000
${\mathrm{NRRes}}_1$	$9.99\times {10}^{-1}$	$9.99\times {10}^{-1}$	$9.99\times {10}^{-1}$
${\mathrm{NRRes}}_2$	$2.49\times {10}^{-7}$	$2.77\times {10}^{-7}$	$9.99\times {10}^{-8}$
${\mathrm{NRRes}}_3$	$6.24\times {10}^{-17}$	$3.62\times {10}^{-17}$	$4.48\times {10}^{-18}$
$\parallel H_1-X\parallel$	$9.99\times {10}^{-1}$	$9.99\times {10}^{-1}$	$9.99\times {10}^{-1}$
$\parallel H_2-X\parallel$	$2.26\times {10}^{-7}$	$3.33\times {10}^{-7}$	$2.75\times {10}^{-7}$
$\parallel H_3-X\parallel$	$1.24\times {10}^{-14}$	$1.25\times {10}^{-14}$	$1.24\times {10}^{-14}$
$R_t$	$37.5\%$	$12.6\%$	$6.8\%$

.

It is seen from the table that Algorithm 1 is efficient to calculate the solution of the DARE. In fact, for different dimensions, the actual error between $H_k$ and the solution X is less than the prescribed accuracy after three iterations, and the derived relative residual is down to a lower level about ${10}^{-17}$ to ${10}^{-18}$. Especially, the value of $R_t$ gradually decreases with the rising scale of n, indicating that the CPU time for iterations takes only a small part of the whole for large-scale problems.

Example 2.

Randomly generate matrices $C_1,C_2,B\in {\mathbb{R}}^{n\times m}$ and define

$$C_1:=\frac{C_1}{\sqrt{2}{\parallel C_1\parallel }^{1/2}},C_2:=\frac{C_2}{\sqrt{2}{\parallel C_2\parallel }^{1/2}},B:=\frac{B}{{\parallel B\parallel }^{1/2}}.$$

Set $R=S=I_m$ and consider the DARE (1) with $A=C_1{C}_2^{*}$, $G=B{B}^{*}$, and $H=I-B{B}^{*}-\frac{1}{2}{C}_2{C}_2^{*}+C_2{C}_1^{*}B{B}^{*}{C}_1{C}_2^{*}$. It is not difficult to see the solution of the DARE is $X=I-B{B}^{*}$. Similarly, the principle of selecting the above matrices is for the convenience of evaluating the error.

We take $n=5000,6000,7000$ to test the error between the true solution and the computed solution. The obtained results together with the NRRes are plotted in Figure 1, Figure 2 and Figure 3. Still, ${\mathrm{R}}_{\mathrm{t}}$ represents the ratio of the iteration time and the total time.

Figure 1, Figure 2 and Figure 3 show that as the number of iterations increases, the NRRes and errors decrease exponentially and Algorithm 1 terminates at the 6th iteration. In all experiments, the preprocessing time for three cases varied from 0.1 to 0.2 s, while the iterative time only took from 0.0032 to 0.0035 s, costing a small part of the whole CPU time. More experiments also indicated that the ratio $R_t$ became smaller as the scale of the problem increased.

Example 3.

This example comes from a proper modification of the circuits from the magneto-quasistatic Maxwell equations ([17,18]). The matrix $S\in {\mathbb{R}}^{632\times 632}$ represents the DC resistance matrix of each current filament (see Figure 4) and $C_1$ as well as $C_2\in {\mathbb{R}}^{n\times 632}$ associated with the mesh matrices. Let $R=1$, $B^{\top }=[1,0,...,0]\in {\mathbb{R}}^{1\times n}$ and $H=I_n$. We randomly generate the matrix $U\in {\mathbb{R}}^{n\times 632}$ and define

$$U=\frac{U}{{\parallel U\parallel }^{1/2}},C_1=U,C_2=U.$$

The tolerance ϵ is taken as ${10}^{-14}$, and the dimensions are $n=i\times {10}^5$ ($i=1,2,...,6$).

For all cases in our experiments, Algorithm 1 was observed attaining the relative residual level below $1.01\times {10}^{-16}$ at the 4-th iteration. The elapsed CPU time and the ratio ${\mathrm{R}}_t$ are plotted in Figure 5, where “${\mathrm{T}}_{pre}$” and “${\mathrm{T}}_{it}$” record the CPU time for the preprocessing and for the iteration, respectively. One can see from the figure that as the scale n rises, the preprocessing time becomes more dominant (about 112 s at $n=600,000$) but the iteration time remains almost unchanged (about 3.5 s for all n). The gradually reduced ratio ${\mathrm{R}}_t$ also illustrates that the main computations of Algorithm 1 when solving the large-scale problems lie in the preprocessing step of $O(n)$ flops, much less than the exponentially increasing one of $O(2^{k}n)$ in [8,9].

5. Conclusions

We have proposed an efficient algorithm to solve the large-scale DAREs with low-ranked matrices A and G and a high-rank matrix H. Compared with the SDA of the complexity $O(2^{k}n)$ in [8,9], the newly developed algorithm only requires preprocessing step of $O(n)$ flops and iteration step of $O({(k+1)}^3{m}^3)$ flops. For large-scale problems with $m\ll n$, the main computations of the whole algorithm lie in the preprocessing step with several matrix multiplications and an economic QR decomposition, while the elapsed CPU time for the iteration part is trivial. Some numerical experiments validate the effectiveness of the proposed algorithm. For future work, we may investigate the possibility of the SDA for solving large-scale DAREs with the structure of sparse-plus-low-rank in A, where the possible difficulty might be understanding the concrete structure of the iterative matrix and knowing how to compute and store it efficiently.