RADI: a low-rank ADI-type algorithm for large scale algebraic Riccati equations

This paper introduces a new algorithm for solving large-scale continuous-time algebraic Riccati equations (CARE). The advantage of the new algorithm is in its immediate and efficient low-rank formulation, which is a generalization of the Cholesky-factored variant of the Lyapunov ADI method. We discuss important implementation aspects of the algorithm, such as reducing the use of complex arithmetic and shift selection strategies. We show that there is a very tight relation between the new algorithm and three other algorithms for CARE previously known in the literature—all of these seemingly different methods in fact produce exactly the same iterates when used with the same parameters: they are algorithmically different descriptions of the same approximation sequence to the Riccati solution.


Introduction
The continuous-time algebraic Riccati equation, where Q = C * C, G = B B * , A ∈ R n×n , B ∈ R n×m , C ∈ R p×n , appears frequently in various aspects of control theory, such as linear-quadratic optimal regulator problems, H 2 and H ∞ controller design and balancing-related model reduction. While the equation may have many solutions, for such applications one is interested in finding a so-called stabilizing solution: the unique positive semidefinite solution X ∈ C n×n such that the matrix A − G X is stable (i.e. all of its eigenvalues belong to C − , the left half of the complex plane). If the pair (A, G) is stabilizable (i.e. rank[A − λI, G] = n, for all λ in the closed right half plane), and the pair (A, Q) is detectable (i.e. (A * , Q * ) is stabilizable), then such a solution exists [11,19]. These conditions are fulfilled generically, and we assume they hold throughout the paper.
There are several algorithms for finding the numerical solution of (1). In the case when n is small enough, one can compute the eigen-or Schur decomposition of the associated Hamiltonian matrix and use an explicit formula for X , see [11,20]. However, if the dimensions of the involved matrices prohibit the computation of a full eigenvalue or Schur decomposition, specialized large-scale algorithms have to be constructed. In such scenarios, Riccati equations arising in applications have additional properties: A is sparse, and p, m n, thus making the matrices Q and G of very-low rank compared to n. In practice, this often implies that the sought-after solution X will have a low numerical rank [3], and allows for construction of iterative methods that approximate X with a series of matrices stored in low-rank factored form. Most of these methods are engineered as generalized versions of algorithms for solving a large-scale Lyapunov equation [10,33], which is a special case of (1) with G = 0.
The alternating directions implicit (ADI) method [35] is a well established iterative approach for computing solutions of Lyapunov and other linear matrix equations. There exists an array of ADI methods [7,8,21,22,27,30,35,36], covering both the ordinary and the generalized case. All of these methods have simple statements and efficient implementations [27]. One particular advantage of ADI methods is that they are very well suited for large-scale problems: the default formulation which works with full-size dense matrices can be transformed into a series of iteratively built approximations to the solution. Such approximations are represented in factored form, each factor having a very small rank compared to the dimensions of the input matrices. This makes ADI methods very suitable for large-scale applications.
Likewise, A ≥ (≤)B refers to A − B ≥ (≤)0. If not stated otherwise, · is the Euclidean vector or subordinate matrix norm, and κ(·) is the associated condition number.

A new low-rank factored iteration
We start this section by stating various methods for solving Lyapunov and Riccati equations, which will be used throughout the paper. First, consider the Lyapunov equation Here we assume that n is much larger than p. When A is a stable matrix, the solution X lya is positive semidefinite. The Lyapunov ADI algorithm [36] generates a sequence of approximations (X lya k ) k to X lya defined by Lyapunov ADI (4) We will assume that the initial iterate X lya 0 is the zero matrix, although it may be set arbitrarily. The complex numbers σ k ∈ C − are called shifts, and the performance of ADI algorithms depends strongly on the appropriate selection of these parameters [30]; this is further discussed in the context of the RADI method in Sect. 4.5. Since each iteration matrix X lya k is of order n, formulation (4) is unsuitable for large values of n, due to the amount of memory needed for storing X lya k ∈ C n×n and to the computational time needed for solving n linear systems in each half-step. The equivalent low-rank algorithm [4,6] generates the same sequence, but represents X lya k in factored form X lya k = Z lya k (Z lya k ) * with Z lya k ∈ C n× pk : low-rank Lyapunov ADI (5) Initially, the matrix Z lya 0 is empty. This formulation is far more efficient for large values of n, since the right-hand side of the linear system in each step involves only the tall-and-skinny matrix R lya k−1 ∈ C n× p . We now turn our attention to the Riccati equation (1). Once again we assume that p, m n, and seek to approximate the stabilizing solution X . Wong and Balakrishnan [37,38] suggest the following quadratic ADI iteration (abbreviated as qADI): quadratic ADI (6) Again, the initial approximation is usually set to X adi 0 = 0. Note that by inserting G = 0 in the quadratic iteration we obtain the Lyapunov ADI algorithm (4). As mentioned in the introduction, we will develop a low-rank variant of this algorithm such that inserting G = 0 will reduce it precisely to (5). To prepare the terrain, we need to introduce two more methods for solving the Riccati equation.
In the small scale setting, the Riccati equation (1) is usually solved by computing the stable invariant subspace of the associated 2n × 2n Hamiltonian matrix (2). To be more precise, let 1 ≤ k ≤ n, and where P k , Q k ∈ C n×k and the matrix Λ k ∈ C k×k is stable. For k = n, the stabilizing solution of (1) is given by X = −Q k P −1 k . In the large scale setting, it is computationally too expensive to compute the entire n-dimensional stable invariant subspace of H . Thus an alternative approach was suggested in [1]: for k n, one can compute only a k-dimensional, stable, invariant subspace and use an approximation given by the formula Clearly, X inv n = X . This approach was further studied and justified in [3], where it was shown that if p = 1 and the shifts used for the qADI iteration coincide with the Hamiltonian eigenvalues of the matrix Λ k . In fact, properties of the approximate solution X inv k given in [3] have lead us to the definition of the low-rank variant of the qADI iteration that is described in this paper.
Here X cay 0 is some initial approximation and σ k ∈ C − are any chosen shifts (Here we have adapted the notation to fit the one of this paper). In Sect. 3 we will show that this method is also equivalent to the qADI and the new low-rank RADI iterations.

Derivation of the algorithm
The common way of converting ADI iterations into their low-rank variants is to perform a procedure similar to the one originally done by Li and White [22] for the Lyapunov ADI method. A crucial assumption for this procedure to succeed is that the matrices participating in the linear systems in each of the half-steps mutually commute for all k. For the Lyapunov ADI method this obviously holds true for the matrices A * + σ k+1 I . However, in the case of the quadratic ADI iteration (6), the matrices A * + σ k+1 I − X adi k G do not commute in general, for all k, and neither do A * + σ k+1 I − X adi k+1/2 G. Thus we take a different approach in constructing the low-rank version. A common way to measure the quality of the matrix Ξ ∈ C n×n as an approximation to the Riccati solution is to compute the norm of its residual matrix The idea for our method is to repetitively update the approximation Ξ by forming a so-called residual equation, until its solution converges to zero. The background is given in the following simple result.

Theorem 1
Let Ξ ∈ C n×n be an approximation to a solution of (1).

Proof
(a) This is a straightforward computation which follows by insertingX = X − Ξ and the formula for the residual of Ξ into (10), see also [2,25]. (b) The first part follows as in (a). If Ξ ≥ 0 andX is a stabilizing solution to (10), then X = Ξ +X ≥ 0 and A − G X =Ã − GX is stable, which makes X the stabilizing solution to (1). Our algorithm will have the following form: 1. Let Ξ = 0. 2. Form the residual equation (10) for the approximation Ξ .
To complete the derivation, we need to specify Step 3 in a way thatX 1 ≥ 0 and R(Ξ +X 1 ) ≥ 0. With these two conditions imposed, Theorem 1 ensures that the residual equation in Step 2 always has a unique stabilizing solution and that the approximation Ξ is kept positive semidefinite and monotonically increasing towards the stabilizing solution of (1). The matrixX 1 fulfilling these conditions can be obtained by computing a 1-dimensional invariant subspace for the Hamiltonian matrix associated with the residual equation, and plugging it into formula (8).
More precisely, assume that R(Ξ ) =C * C ≥ 0, G = B B * , and that r, q ∈ C n satisfy where λ ∈ C − is such that −λ is not an eigenvalue ofÃ. Equivalently, From the second equation we get q = (Ã * + λI ) −1C * (Cr). Let Multiply the first equation by q * from the left, and the transpose of the second by r from the right; then add the terms to obtain Expression (8) has the formX 1 = −q(q * r ) −1 q * . When p = 1 andCr = 0, the terms containingCr cancel out, and we get The derivation above is valid when λ is an eigenvalue of the Hamiltonian matrix, similarly as in [3,Theorem 7] which also studies residual Riccati equations related to invariant subspaces of Hamiltonian matrices. Nevertheless, the expression (11) is welldefined even when λ is not an eigenvalue of the Hamiltonian matrix, and for p > 1 as well. The Hamiltonian argument here serves only as a motivation for introducing (11), and the following proposition shows that the desired properties of the updated matrix Ξ +X 1 still hold for any λ in the left half-plane which is not an eigenvalue of −Ã, and for all p. Proposition 1 Let Ξ ≥ 0 be such that R(Ξ ) =C * C ≥ 0, and let λ ∈ C − not be an eigenvalue of −Ã. The following holds true for the update matrixX 1 as defined in (11):
We use these expressions to obtain: We are now ready to state the new RADI algorithm. Starting with the initial approximation X 0 = 0 and the residual R(X 0 ) = C * C, we continue by selecting a shift σ k ∈ C − and computing the approximation (11), adapted to approximate the solution of the residual equation with Ξ = X k−1 , i.e. withÃ = A − B B * X k−1 and C * = R k−1 . Proposition 1 provides a very efficient update formula for the low-rank factor R k of the residual. The whole procedure reduces to the following: (12) Note that any positive semi-definite X 0 can be used as an initial approximation, as long as its residual is positive semi-definite as well, and its low-rank Cholesky fac-torization can be computed. From the derivation of the RADI algorithm we have that When p = 1 and all the shifts are chosen as eigenvalues of the Hamiltonian matrix associated with the initial Riccati equation (1), the update described in Proposition 1 reduces to [3,Theorem 5]. Thus in that case, the RADI algorithm reduces to the invariant subspace approach (8).
Furthermore, iteration (12) clearly reduces to the low-rank Lyapunov ADI method (5) when B = 0; in that case Y k = I . The relation to the original qADI iteration (6) is not clear unless p = 1 and the shifts are chosen as eigenvalues of H , in which case both of these methods coincide with the invariant subspace approach. We discuss this further in the following section.

Equivalences with other Riccati methods
In this section we prove that all Riccati solvers introduced in Sect. 2 in fact compute exactly the same iterations, which we will refer to as the Riccati ADI iterations in the remaining text. This result is collected in Theorem 2; we begin with a simple technical lemma that provides different representations of the residual factor.
Then R Proof From the definition of the RADI iteration (12) it is obvious that R k .

Theorem 2 If the initial approximation in all algorithms is zero, and the same shifts are used, then for all k,
If rank C = 1 and the shifts are equal to distinct eigenvalues of H , then for all k, Proof We first use induction to show that X k = X adi k , for all k. where First, note that (14) can be rewritten as Subtracting this from the expression for the Riccati residual, Equation (13) can be reorganized as Replace the second term on the right-hand side with the right-hand side of (14). Thus, it remains to prove or after some rearranging, and by using (15), The right-hand side of (16) is, thus, equal to which turns out to be precisely the same as the left-hand side of (16) once we use the identity This completes the proof of X k = X adi k . Next, we use induction once again to show X and suggestively introduce X and This is the same relation as the one defining X adi k−1/2 , and thus X cay k−1/2 = X adi k−1/2 . Next, equating the leftmost and the rightmost matrix in (17), it follows that Multiply (21) from the right by M −1 k to obtain and multiply (20) from the left by X cay k−1/2 and from the right by Adding (22) and (23) yields which is the same as the defining equation for X adi k . Thus X cay k = X adi k , so both the Cayley subspace iteration and the qADI iteration generate the same sequences.
In the case of rank C = 1 and shifts equal to the eigenvalues of H , the equality X inv k = X adi k is already shown in [3]. Equality among the iterates generated by the other methods is a special case of what we have proved above.
It is interesting to observe that [23] also provides a low-rank variant of the Cayley subspace iteration algorithm: there, formulas for updating the factors of X ∈ C n× pk and Y cay k ∈ C pk× pk , are given. The contribution of [24] was to show that the same formulas can be derived from a control-theory point of view. The main difference in comparison to our RADI variant of the low-rank Riccati ADI iterations is that, in order to compute Z cay k , one uses This way, the need for using the Sherman-Morrison-Woodbury formula is avoided. However, as a consequence, the matrix Y cay k looses the block-diagonal structure, and its update formula becomes much more involved. Also, it is very difficult to derive a version of the algorithm that would use real arithmetic. Another disadvantage is the computation of the residual: along with Z cay k and Y cay k , one needs to maintain a QRfactorization of the matrix [C * A * Z cay k Z cay k ], which adds significant computational complexity to the algorithm.
Each of these different statements of the same Riccati ADI algorithm may contribute when studying theoretical properties of the iteration. For example, directly from our definition (12) of the RADI iteration it is obvious that Also, the fact that the residual matrix R(X k ) is low-rank and its explicit factorization follows naturally from our approach. On the other hand, approaching the iteration from the control theory point of view as in [24] is more suitable for proving that the non-Blaschke condition for the shifts, is sufficient for achieving the convergence when A is stable, i.e. lim k→∞ X k = X.
We conclude this section by noting a relation between the Riccati ADI iteration and the rational Krylov subspace method [34]. It is easy to see that the RADI iteration also uses the rational Krylov subspaces as the basis for approximation. This fact also follows from the low-rank formulation for X cay k as given in [23], so we only state it here without proof.

Proposition 2 For a matrix M, (block-)vector v and a tuple
denote the rational Krylov subspace generated by M and the initial vector v. Then the columns of X k belong to From the proposition we conclude the following: if U k contains a basis for the rational Krylov subspace K (A * , C * , − → σ k ), then both the approximation X kry k of the Riccati solution obtained by the rational Krylov subspace method and the approximation X adi k obtained by a Riccati ADI iteration satisfy The columns of both X kry k and X adi k belong to the same subspace, and the only difference between the methods is the choice of the linear combination of columns of U k , i.e. the choice of the small matrix Y k . The rational Krylov subspace method [34] generates its Y kry k by solving the projected Riccati equation, while the Riccati ADI methods do it via direct formulas such as the one in (12).

Implementation aspects of the RADI algorithm
There are several issues with the iteration (12), stated as is, that should be addressed when designing an efficient computational routine: how to decide when the iterates X k have converged, how to solve linear systems with matrices A * − X k−1 G + σ k I , and how to minimize the usage of complex arithmetic. In this section we also discuss the various shift selection strategies.

Computing the residual and the stopping criterion
Tracking the progress of the algorithm and deciding when the iterates have converged is very simple, and can be computed cheaply thanks to the expression R(X k ) = R k R * k = R * k R k . This is an advantage compared to the Cayley subspace iteration, where computing R(X k ) is more expensive because a low-rank factorization along the lines of Proposition 1 (b) is currently not known. The RADI iteration is stopped once the residual norm has decreased sufficiently relative to the initial residual norm CC * of the approximation X 0 = 0. Algorithm 1: The RADI iteration using complex arithmetic.
Obtain the next shift σ ; 4 if first pass through the loop then

Solving linear systems in RADI
During the iteration, one has to evaluate the expression Here the matrix A is assumed to be sparse, while There are different options on how to solve this linear system; if one wants to use a direct sparse solver, the initial expression can be adapted by using the Sherman-Morrison-Woodbury (SMW) formula [14]. We introduce the approximate feedback matrix K k := X k B and update it during the RADI iteration: k also appears in the update of the residual factor, and that V * k B appears in the computation ofỸ k , so both have to be computed only once. The initial expression is rewritten as Thus, in each RADI step one needs to solve a linear system with the coefficient matrix A * + σ k I and p + m right hand sides. A very similar technique is used in the low-rank Newton ADI solver for the Riccati equation [7,9,15,29]. In the equivalent Cayley subspace iteration [23,24], linear systems defined by A * + σ k I and only p right hand sides have to be solved, which makes their solution less expensive than their counterparts in RADI.
The RADI algorithm, implementing the techniques described above, is listed in Algorithm 1. Note that, if only the feedback matrix K is of interest, e.g. if the CARE arises from an optimal control problem, there is no need to store the whole low-rank factors Z , Y since Algorithm 1 requires only the latest blocks to continue. This is again similar to the low-rank Newton ADI solver [7], and not possible in the current version of the Cayley subspace iteration.

Reducing the use of complex arithmetic
To increase the efficiency of Algorithm 1, we reduce the use of complex arithmetic. We do so by taking shift σ k+1 = σ k immediately after the shift σ k ∈ C \ R has been used, and by merging these two consecutive RADI steps into a single one. This entire procedure will have only one operation involving complex matrices: a linear solve with the matrix A * − X k−1 G + σ k I to compute V k . There are two key observations to be made here. First, by modifying the iteration slightly, one can ensure that the matrices K k+1 , R k+1 , and Y k+1 are real and can be computed by using real arithmetic only, as shown in the upcoming technical proposition. Second, there is no need to compute V k+1 at all to proceed with the iteration: the next matrix V k+2 will once again be computed by using the residual R k+1 , the same way as in Algorithm 1.
The residual factor R k+1 and the matrix K k+1 can be computed as Proof The basic idea of the proof is similar to [5, Theorem 1]; however, there is a major complication involving the matrix Y , which in the Lyapunov case is simply equal to the identity. Due to the technical complexity of the proof, we only display key intermediate results.
To simplify notation, we use indices 0, 1, 2 instead of k − 1, k, k + 1, respectively. We start by taking the imaginary part of the defining relation for R We use this and the Sherman-Morrison-Woodbury formula to compute V 2 : where we used Im (σ ) A −1 0 V 1 = − Im (V 1 ) to obtain the last line. Next, We first compute T −1 by using the SMW formula once again: , and, applying the congruence transformation with T −1 to By using (24), it is easy to show Inserting this into the formula forŶ 2 , all terms containing inverses ofỸ 1 and S cancel out. By rearranging the terms that do appear in the formula, we get the expression from the claim of the proposition. Deriving the formulae for R 2 and K 2 is straightforward: where the last line is due to the structure of the matrix T . Since the matrixŶ 2 is real, so are K 2 and R 2 .

RADI iteration for the generalized Riccati equation
Before we state the final implementation, we shall briefly mention the adaptation of the RADI algorithm for handling generalized Riccati equations Multiplying (25) by E − * from the left and by E −1 from the right leads to The generalized RADI algorithm is then easily derived by running ordinary RADI iterations for the Riccati equation (26), and deploying standard rearrangements to Algorithm 2: The RADI iteration, with reduced use of complex arithmetic.
Input: matrices A, E ∈ R n×n , B ∈ R n×m , C ∈ R p×n . Output: approximation X ≈ ZY −1 Z * for the solution of the generalized equation The matrices Z and Y are real.
Obtain the next shift σ ; 4 if first pass through the loop then lower the computational expense (such as solving systems with the matrix A * + σ E * instead of (AE −1 ) * + σ I ). Algorithm 2 shows the final implementation, taking into account Proposition 3 and the handling of generalized Eq. (25).

Shift selection
The problem of choosing the shifts in order to accelerate the convergence of the Riccati ADI iteration is very similar to the one for the Lyapunov ADI method. Thus we apply and discuss the techniques presented in [3,6,27,30] in the context of the Riccati equation, and compare them in several numerical experiments. It appears natural to employ the heuristic Penzl shifts [27]. There, a small number of approximate eigenvalues of A are generated. From this set the values which lead to the smallest magnitude of the rational function associated to the ADI iteration are selected in a heuristical manner. Simoncini and Lin [23] have shown that the convergence of the Riccati ADI iteration is related to a rational function built from the stable eigenvalues of H . This suggests to carry out the Penzl approach, but to use approximate eigenvalues for the Hamiltonian matrix H instead of A. Note that, due to the low rank of Q and G, we can expect that most of the eigenvalues of A are close to the eigenvalues of H , see the discussion in [3]. Thus in many cases the Penzl shifts generated by A should suffice as well. Penzl shifts require significant preprocessing computation: in order to approximate the eigenvalues of M = A or M = H , one has to build Krylov subspaces with matrices M and M −1 . All the shifts are computed in this preprocessing stage, and then simply cycled during the RADI iteration. Here, we will mainly focus on alternative approaches generating each shift just before it is used. This way we hope to compute a shift which will better adapt to the current stage of the algorithm.

Residual Hamiltonian shifts
One such approach is motivated by Theorem 1 and the discussion about shift selection in [3]. The Hamiltonian matrixH is associated to the residual equation (10), where Ξ = X k is the approximation after k steps of the RADI iteration. If (λ, r q ) is a stable eigenpair ofH , and σ k+1 = λ is used as the shift, then In order to converge as fast as possible to X , it is better to choose such an eigenvalue λ for which the update is largest, i.e. the one that maximizes q(q * r ) −1 q * . Note that as the RADI iteration progresses and the residual matrix R(Ξ ) =C * C converges to zero, the structure of eigenvectors ofH that belong to its stable eigenvalues is such that q becomes smaller and smaller. Thus, one can further simplify the shift optimality condition, and use the eigenvalue λ such that the corresponding q has the largest norm-this is also in line with the discussion in [3].
However, in practice it is computationally very expensive to determine such an eigenvalue, since the matrixH is of order 2n. We can approximate its eigenpairs through projection onto some subspace. If U is an orthonormal basis of the chosen subspace, then is the projected residual Riccati equation with the associated Hamiltonian matrix Approximate eigenpairs ofH are (λ, Ur Uq ), where (λ, rq ) are eigenpairs ofH proj . Thus, a reasonable choice for the next shift is suchλ, for which Uq = q is the largest.
We still have to define the subspace span{U }. One option is to use V k (or, equivalently, the last p columns of the matrix Z k ), which works very well in practice unless p = 1. When p = 1, all the generated shifts are real, which can make the convergence slow in some cases. Then it is better to choose the last columns of the matrix Z k ; usually already = 2 or = 5 or a small multiple of p will suffice. An ultimate option is to use the entire Z k , which we denote as = ∞. This is obviously more computationally demanding, but it provides fast convergence in all cases we tested.

Residual minimizing shifts
The two successive residual factors are connected via the formula Our goal is to choose the shifts so that the residual drops to zero as quickly as possible. Locally, once X k is computed, this goal is achieved by choosing σ k+1 ∈ C − so that R k+1 is minimized. This concept was proposed in [6] for the low-rank Lyapunov and Sylvester ADI methods and refined later in [18]. In complete analogy, we define a rational function f in the variable σ by and wish to find is also a function of σ . Since f involves large matrices, we once again project the entire equation to a chosen subspace U , and solve the optimization problem defined by the matrices of the projected problem. The optimization problem is solved numerically. Efficient optimization solvers use the gradient of the function f ; after a laborious computation one can obtain an explicit formula for the case p = 1: Here For p > 1, a similar formula can be derived, but one should note that the function f is not necessarily differentiable at every point σ , see, e.g., [26]. Thus, a numerically more reliable heuristic [18] is to artificially reduce the problem once again to the case p = 1. This can be done in the following way: let v denote the right singular vector corresponding to the largest singular value of the matrix R k . Then R k ∈ C n× p in (28) is replaced by the vector R k v ∈ C n . Since numerical optimization algorithms usually require a starting point for the optimization, the two shift generating approaches may be combined: the residual Hamiltonian shift can be used as the starting point in the optimization for the second approach. However, from our numerical experience we conclude that the additional computational effort invested in the post-optimization of the residual Hamiltonian shifts often does not contribute to the convergence. The main difficulty is the choice of an adequate subspace U such that the projected objective function approximates (28) well enough. This issue requires futher investigation. The rationale is given in the following example.
Example 1 Consider the Riccati equation given in Example 5.2 of [32]: the matrix A is obtained by the centered finite difference discretization of the differential equation on a unit cube with 22 nodes in each direction. Thus A is of order n = 10648; the matrices B ∈ R n×m and C ∈ R p×n are generated at random, and in this example we set m = p = 1 and B = C * .
Suppose that 13 RADI iterations have already been computed, and that we need to compute a shift to be used in the 14th iteration. Let the matrix U contain an orthonormal basis for Z 13 . Figure 1a shows a region of the complex plane; stars are at locations of the stable eigenvalues of the projected Hamiltonian matrix (27). The one eigenvalue chosen as the residual Hamiltonian shift σ ham is shown as 'x'. The residual minimizing shift σ opt is shown as 'o'. Each point σ of the complex plane is colored according to the , where R proj is the residual for the projected Riccati equation. The ratio ρ proj (σ ham ) ≈ 0.54297 is not far from the optimal ratio ρ proj (σ opt ) ≈ 0.53926.
On the other hand, Fig. 1b shows the complex plane colored according to ratios for the original system of order 10648, ρ(σ ) = R 14 (σ ) / R 13 . Neither of the values ρ(σ ham ) ≈ 0.71510 and ρ(σ opt ) ≈ 0.71981 is optimal, but they both offer a reasonable reduction of the residual norm in the next step. In this case, σ ham turns out even to give a slightly better residual reduction for the original equation than σ opt , making the extra effort in running the numerical optimization algorithm futile.

Numerical experiments
In this section we show a number of numerical examples, with several objectives in mind. First, our goal is to compare different low-rank implementations of the Riccati ADI algorithm mentioned in this paper: the low-rank qADI proposed in [37,38], the Cayley transformed subspace iteration [23,24], and the complex and real variants of the RADI iteration (12). Second, we compare performance of the RADI approach against other methods for solving large-scale Riccati equations, namely the rational Krylov subspace method (RKSM) [34], the extended block Arnoldi (EBA) method [16,32], and the Newton-ADI algorithm [7,9,15].
Finally, we discuss various shift strategies for the RADI iteration described in the previous section.
The numerical experiments are run on a desktop computer with a four-core Intel Core i5-4690K processor and 16GB RAM. All algorithms and testing routines are implemented and executed in MATLAB R2014a, running on Microsoft Windows 8.1.
Example 2 Consider again the Riccati benchmark CUBE from Example 1. We use three versions of this example: the previous setting with n = 10648, m = p = 1 and m = p = 10, and later on a finer discretization with n = 74088 and m = 10, p = 1. Table 1 collects timings in seconds for four different low-rank implementations of the Riccati ADI algorithm. The table shows only the time needed to run 80 iteration steps; time spent for computation of the shifts used by all four variants is not included (in this case, 20 precomputed Penzl shifts were used). All four variants compute exactly the same iterates, as we have proved in Theorem 2. Clearly, the real variant of iteration (12), implemented as in Algorithm 2, outperforms all the others. 1 Thus we use this implementation in the remaining numerical experiments. The RADI algorithms mostly obtain the advantage over the Cayley subspace iteration because of the cheap computation of the residual norm. In the latter algorithm, costly orthogonalization procedures are required for this task, and after some point these compensate the computational gains from the easier linear systems (cf. Sect. 4.2). Also, the times for the algorithm of Wong and Balakrishnan shown in the table do not include the (very costly) computation of the residuals at all, so their actual execution times are even higher.
For the RKSM, we have implemented the algorithm so that the use of complex arithmetic is minimized by merging two consecutive steps with complex conjugate shifts [28]. We use the adaptive shift strategy, implemented as described in [13]. EBA and RKSM require the solution of a projected CARE which can become expensive if p > 1. Hence, this small scale solution is carried out only periodically in every 5th or every 10th step-in the results, we display the variant that was faster. For all methods, the threshold for declaring convergence is reached once the relative residual is less than tol = 10 −11 . A summary of the results for all the different methods and strategies is shown in Table 3. The column "final subspace dimension" displays the number of columns of the matrix Z , where X ≈ Z Z * is the final computed approximation. Dividing this number by p (for EBA, by 2 p), we obtain the number of iterations used in a particular method. Just for the sake of completeness, we have also included a variant of the Newton-ADI algorithm [7] with Galerkin projection [9]. Without the Galerkin projection, the Newton-ADI algorithm could not compete with the other methods. The recent developments from [15], which make the Newton-ADI algorithm more competitive, are beyond the scope of this study.
It is interesting to analyze the timing breakdown for RADI and RKSM methods. These timings are listed in Table 2 for the CUBE example with m = p = 10 where a significant amount of time is spent for tasks other than solving linear systems.
As increases, the cost of computing shifts in RADI increases as well-the projection subspace gets larger, and more effort is needed to orthogonalize its basis and compute the eigenvalue decomposition of the projected Hamiltonian matrix. This effort is, in the CUBE benchmark, awarded by a decrease in the number of iterations. However, there is a trade-off here: the extra computation does outweigh the saving in the number of iterations for sufficiently large . Convergence history for CUBE is plotted in Fig. 2; to reduce the clutter, only the selected few methods are shown.  The fact that in each step RADI solves linear systems with p + m right hand side vectors, compared to only p vectors in RKSM, may become noticable when m is larger than p. This effect is shown in Table 3 for CUBE with m = 10 and p = 1. Unlike these two methods, EBA can precompute the LU factorization of A, and win by a large margin in this test case.   Example 3 Next, we run the Riccati solvers for the well-known benchmark example CHIP. All coefficient matrices for the Riccati equation are taken as they are found in the Oberwolfach Model Reduction Benchmark Collection [17]. Here we solve the generalized Riccati equation (25). The cost of precomputing shifts is very high in case of CHIP. One fact not shown in the table is that all algorithms which compute shifts dynamically have already solved the Riccati equation before "RADI-Penzl" has even started.

Example 4
We use the IFISS 3.2. finite-element package [31] to generate the coefficient matrices for a generalized Riccati equation. We choose the provided example T-CD 2 which represents a finite element discretization of a two-dimensional convection diffusion equation on a square domain. The leading dimension is n = 66049, with E symmetric positive definite, and A non-symmetric. The matrix B consists of m = 5 randomly generated columns, and C = [C 1 , 0] with random C 1 ∈ R 5×5 ( p = 5).
In this example, the RADI iteration with Penzl shifts converges very slowly. The RADI iteration with dynamically generated shifts are quite fast, and the final subspace dimension is smallest among all methods. On the other hand, the version with residual minimizing shifts does not converge for = 6 p, ∞: it quickly reaches the relative residual of about 10 −7 , and then gets stuck by continually using shifts very close to zero. Figure 3 shows the convergence history for some of the used shift strategies.

Example 5
The example RAIL is a larger version of the steel profile cooling model from the Oberwolfach Model Reduction Benchmark Collection [17]. A finer finite element discretization was used for the heat equation resulting in a generalized CARE n = 317377, m = 7, p = 6, and E and A symmetric positive and negative definite, respectively. Once again, there is a trade-off between (questionably) better shifts with larger and faster computation with lower .

Example 6
The final example LUNG from the UF Sparse Matrix Collection [12] models temperature and water vapor transport in the human lung. It provides matrices with leading dimension n = 109460, E = I , A nonsymmetric, and B, C are generated as random matrices with m = p = 10. This example shows the importance of proper shift generation: precomputed shifts are completely useless, while dynamically generated ones show different rates of success. The projection based methods (RKSM, EBA) encountered problems at the numerical solution of the projected ARE. Either the complete algorithm broke down or convergence speed was reduced. Similar issues were encountered at the Galerkin acceleration stage in the Newton-ADI method.
Let us summarize the findings from these and a number of other numerical examples we used to test the algorithms. Clearly, using dynamically generated shifts for the RADI iteration has many benefits compared to the precomputed Penzl shifts. Not only that the number of iterations and running time are reduced, but the convergence is more reliable. Further, there frequently exists a small value of for which one or both of the dynamical shift strategies converge in a number of iterations comparable to runs with = ∞, and in far less time. However, an a-priori method of determining a sufficiently small with such properties is still to be found, and a topic of our future research. RADI appears to be quite competitive with other state of the art algorithms for solving large scale CAREs. It frequently generates solutions using the lowest dimensional subspace. Since RADI iterations do not require any orthogonalization nor solving projected CAREs, the algorithm may outperform RKSM and EBA in problems where the final subspace dimension is high. On the other hand, the later methods may have an advantage when the running time is dominated by solving linear systems. It seems that for now, there is no single algorithm of choice that would consistently and reliably run fastest.

Conclusion
In this paper, we have presented a new low-rank RADI algorithm for computing solutions of large scale Riccati equations. We have shown that this algorithm produces exactly the same iterates as three previously known methods (for which we suggest the common name "Riccati ADI methods"), but it does so in a computationally far more efficient way. As with other Riccati solvers, the performance is heavily dependent on the choice of shift parameters. We have suggested several strategies on how this may be done; some of them show very promising results, making the RADI algorithm competitive with the fastest large scale Riccati solvers.