An Improvement of the Newton Method for Solving Symmetric Algebraic Riccati Equations

In this work, we analyze a strategy for solving symmetric algebraic Riccati equation based on the use of efficient high-order iterative scheme. This iterative scheme is more efficient than Newton’s method. Then, we propose two iterative two-stage predictor–corrector schemes using an iterative scheme with good accessibility as the predictor iteration and a high-order iterative scheme as the corrector iteration. The iterative schemes constructed turn out to be competitive compared to the commonly used Newton’s method. The efficiency of these methods is illustrated by a numerical example.


Introduction
The study of algebraic Riccati equation (ARE) is motivated by the important role played in many applications from different areas, such as optimal filter design and control theory [20], queueing models [23], numerical solutions of the transport theory [14], problems with or without symmetric constraints, etc. The algebraic Riccati equation is given by ( [18]): where X ∈ C m×n is the unknown, and the coefficients are A ∈ C n×n , B ∈ C m×m , D ∈ C n×m and C ∈ C m×n . Equation (1) is known as the nonsymmetric algebraic Riccati equation (NARE), which is distinguished from the symmetric one: where A * denotes the transpose conjugate of the matrix A. The symmetric term refers to a matrix equation R(X) = 0 such that R(X) * = R(X * ). In particular, Eq. (1) turns to a symmetric equation if A = B * , D = D * and C = C * . Besides, if m = n, then Eq. (2) is known as a continuous-time algebraic Riccati equation (CARE). In this case, the solution X of interest is a Hermitian matrix. A typical problem where a (CARE) appears involved is the linear-quadratic optimal control problems.
In the literature can be found different techniques for solving these equations. Many papers dedicated to solve ARE equations are based on algebraic techniques and theory of matrices, see for instance [19] where the Schur method was introduced. And, other techniques like the use of the hierarchical matrices for solving large scale Riccati equations, see [10]. In fact, in the large and sparse case, it is a common approach to use Newton's method to find a solution of (2), see for instance [7,8,12,16,21,24] and references given therein.
For solving (CARE), at every Newton step, a Lyapunov equation ( [18]) must be solved: where F and W are given matrices.
In this work, we analyze a strategy for solving symmetric algebraic Riccati equation based on the use of efficient high-order iterative scheme [3,6]. We propose two iterative two-stage predictor-corrector schemes [17] using an iterative scheme with good accessibility as the predictor iteration and a high-order iterative scheme as the corrector iteration which accelerates the convergence. An advantage of this strategy is that at every step, the matrix of the Lyapunov equations which must be solved, is the same. Thus, the use of this iterative scheme reduce the number of iterations and hence the total amount of work.
Throughout the work, we suppose the pair (A, D) is stabilizable, that is a (feedback) matrix K ∈ R n×n exists such that all eigenvalues of A − DK are in the open left half-plane. Note that the matrix K can be chosen to be symmetric if the matrix D is symmetric (see [18]). Under these conditions, the existence of Hermitian solutions X of (2) can be characterized using spectral properties of the matrix −A D C A * , see [18]. For this, we denote H the set of Hermitian matrices in R n×n . For any matrix, norm H is a Banach space, and R is a mapping from H into itself. From now on, we denote by · a general norm in the set of Hermitian matrices. Then, the first Frèchet derivative of R at a matrix X ∈ H is a linear map R (X) : H → H, given by Also, the second derivative at X ∈ H, is a bilinear map R (X) : H × H → H, given by One of the most well-known iterative schemes for solving Eq. (2) is the Newton method [1,15]: Thus, taking into account (3), to approximate a solution of equation (2) via the Newton method is equivalent to solve the Lyapunov equation: Our main aim is to approximate a solution of the equation (2), improving the results obtained by the classic Newton's method. To do this, first, we think of an iterative scheme that improves the quadratic speed of convergence that Newton's method has. But, obviously, this is not enough, we must consider an iterative scheme with reduced operational cost thus, it is more efficient than Newton's method. Thus, we consider the M5 method [4] that has local fifth order of convergence and reduced operational cost. In fact, we prove that it is a more efficient iterative scheme than Newton's method. However, when choosing an iterative scheme, in addition to its speed of convergence and its operational cost, there is another very important aspect, which is its accessibility. The measure of this accessibility is given by the set of starting points that make the iterative scheme convergent. Therefore, a commonly used procedure is to obtain a local convergence result for the iterative scheme. This result gives us the well-known convergence ball [2] that indicates a domain of starting points that make the iterative scheme convergent. Thus, comparing the convergence balls of the iterative schemes we compare their accessibility. Then we check that the method M5 has a reduced accessibility and therefore we modify this scheme. Thus, to solve this problem, we consider two-stage predictor-corrector iterative schemes. Thus, first, it is applied to an iterative scheme which has a wide accessibility region, and second, it is applied in the M5 method with better computational efficiency than Newton's method. In this way, we obtain two interesting improvements of the M5 method considering the Newton method and the modified Newton method [15] as predictor schemes.
The paper is organized as follows-First, in Sect. 2, we study the M5 method, proving that it has local fifth order of convergence for solving Eq. (2). Next, in Sect. 3, from a local convergence result given for the iterative scheme M5, we study its domain of parameters and we prove that its accessibility is poor. Then, we justify a first improvement of the M5 method through an iterative two-stage predictor-corrector scheme. In Sect. 4, we obtain a second improvement of the M5 method by means of the modified Newton method from a reduction of the operational cost. To finish, in Sect. 5, we use the predictor-corrector iterative schemes considered to approximate a solution of a particular symmetric ARE of the benchmark collection. We show that these iterative schemes considered are competitive with respect to Newton's

A High-Order Iterative Scheme: M5 Method
In order to accelerate the speed of convergence, we consider the efficient fifth-order iterative scheme M5: Given an initial guess X 0 ∈ H, given that the maps R (X n ) are all invertible. Observe that, from the aforesaid, the iteration (5) is equivalent to solve a single Lyapunov equation, L F (X) = W , with different independent terms: where L F (X) := XF + F * X with F and W given matrices and W is Hermitian if X k is so. In this section, we prove that method (5), under certain conditions at X, is convergent to the solution X of (2) with at least local fifth order of convergence.
In what follows, we denote E k = X k − X the error in the k-th iteration, the equation E k+1 = ME p k , where M is a p-linear operator, M : H×· · · (p · · ·× H → H, is called the error equation and p is the local order of convergence. Notice that E p k is E k × · · · (p · · · × E k . Next result proves the local order of convergence of the M5 method, using the Taylor expansions, and obtain the error equation.  (2) and X a solution of (2). Suppose that A − D X is stable and R is nonsingular at X. Then, the sequence {X k }, given by the M5 method, converges to X with local fifth order of convergence. Moreover, the error equation satisfies Proof. We expand R and R in the Taylor series around the solution X. We denoting by and . Next, taking into account that and taking norms, it follows One of the most important aspects to consider when choosing an iterative scheme to solve a nonlinear equation is its efficiency, see [9,11]. The parameters that are taken into account to analyze the efficiency of an iterative scheme are the order of convergence and the operational cost.
A widely used efficiency index is the computational efficiency (CE), see [25]. This efficiency index is defined by CE = ord 1/op , where ord is the order of convergence and op is the amount of operations per iteration associated with the method.

Computational Efficiency
In this section, we analyze the computational efficiency of method (6) and compare it with that of Newton method.
To analyze the computational efficiency, we need to know the order of operations that these methods perform per step.
On the other hand, Newton's method requires solving a Lyapunov equation at each step. Applying the M5 method involves a Lyapunov equation with three different independent terms.
One of the most known and effective algorithms to the sequential solution of Lyapunov equation with different independent terms is the Bartels-Stewart [5]. The main idea of this algorithm is to reduce the matrix system to Schur or upper Hessenberg form. In this way, it transforms the Lyapunov equation into a triangular system which can be solved efficiently by forward or backward substitutions. Thus, the algorithm leaves the matrix system XF + F * X = W , transformed into the form

A First Improvement for M5 Method
Another interesting aspect when selecting an iterative scheme is its accessibility, that is, the set of starting matrices that make the iterative scheme a convergent process. This fact is measured from a local convergence result that shows us the convergence ball of an iterative scheme.

Local Convergence
Usually, the local convergence results for iterative schemes require conditions on the operator R and a solution X of Eq. (2). Notice that a local result provides the convergence ball which we denote by B( X, S). From the value S, the convergence ball gives information about the accessibility of the solution X, since the convergence of the iterative scheme is ensured from any starting point belonging to the ball B( X, S).
We started our study obtaining some previous results.
Proof. Taking a Hermitian matrix X, to obtain the existence of the lineal operator R ( X) −1 , we have to solve a Lyapunov equation in this way It is well known [18], that this equation has solution if A − D X is stable, that is, all their eigenvalues have real negative part. In this case, the solution is known: On the other hand, as D = D * , C = C * and the pair ( Hence, from now on, we consider that [R ( To prove that {X n }, given by (5), converges to X, we will previously prove that for each element of the sequence, it is verified that there exists [R(X n )] −1 and X n ∈ B( X, S).

Lemma 3. Under conditions of previous Lemma
Proof. Firstly, we have Second, from by the Perturbation Lemma [21], then we get the thesis.
From the algorithm (5) and the Taylor expansions, it is easy to check the following technical result. Lemma 4. If X n , Y n , Z n ∈ H, then the following equalities are satisfied for n 0 :

Lemma 5.
Under conditions of Lemma 2, assuming that X n ∈ B( X, S) are well defined, the following items are satisfied for n ≥ 1: and the auxiliary real functions: Ψ 0 (t) = 3t, Proof. From algorithm (5) and the Taylor expansions, it is easy to check Now, taking norms, we obtain which proves (i n ).
Proof. Note that if β D S < 0.14, then Ψ 2 (Δ(S)) < 1, and therefore { X k − X } is a strictly decreasing sequence of positive real numbers. Therefore, the thesis is obtained.

Accessibility
Notice that, the local convergence result, Theorem 6, is based on demanding conditions to β D S. This result provides the so-called domain of parameters [1] corresponding to the conditions required that guarantee the convergence of sequence {X n } to a solution X of symmetric ARE, (2). In this case, the condition is β D S < 0.14. Therefore, the domain of parameters measures the accessibility of the M5 method, taking x = β D and y = S. That is, the location of starting approximations, X 0 , from which the M5 method converges to a solution of symmetric ARE, (2). Thus, the accessibility of the M5 method from Theorem 6 is given by its domain of parameters, that we can see in Fig. 3. It is clear that to compare the accessibility of two iterative schemes, we must consider two local convergence results in the same conditions for the operator considered. An interesting local result for Newton's method is given in [13] Theorem 7. Let X be a Hermitian solution of (2) and A−D X is stable.Then, the sequence {X n } generated by the Newton method (4) converges to X from every matrix X 0 ∈ B( X, S), with β D S < 1 3 . Thus, if we consider the previous result for Newton method (4), in the same conditions that Theorem 6, its domain of parameters is given by As we can see in Fig. 4, there is an important difference between the domains of parameters associated with these two methods. And therefore, we can affirm that the accessibility of Newton's method is greater than that of the M5 method.

A Predictor-Corrector Iterative Scheme
As we have obtained previously, the iterative scheme M5 has better computational efficiency than Newton's method. However, in our analysis of accessibility, it has been shown that the iterative scheme M5 has a reduced accessibility region, contrary to what happens to Newton's method. Next we build an iterative scheme that uses the best features of both the M5 and Newton iterative schemes. Thus, we consider a two-stage predictor-corrector iterative scheme. Firstly it is applied Newton's method, which has a wide accessibility region, and second, it is applied the M5 method with better computational efficiency than Newton's method.
Given an initial guess X 0 ∈ H, whose expression, according to the practical application by solving Lyapunov equations, is given as follows: Given an initial guess X 0 ∈ H, Obviously, the study of the local convergence of this iterative scheme requires the determination of the existence of the value N 0 hence we can ensure that the iteration W 0 is in the accessibility region of the iterative scheme M5.
Let X 0 ∈ D N , D N given in (10), that is, X 0 ∈ B( X, S) with β||D||S < 1/3. Being X, a Hermitian solution of the equation of (2), where A − D X is stable.
First, from previous Lemma 3, for n 1 if X n ∈ B( X, S), with β||D||S < 1/3, then there exists the operator [R (X n )] −1 . Secondly, from Taylor series we have: Thirdly, from the algorithm of Newton's iterative scheme (4), we have: Thus, the following decomposition is verified: Hence, taking norms, if X n ∈ B( X, S), we obtain Therefore, under the previous conditions, for n = 0, we have As Δ(S) < 1, then X 1 ∈ B( X, S). To continue, by applying an inductive procedure, it follows that X n ∈∈ B( X, S) for all n 1. Moreover, Now, as Δ(S) < 1, there will be a value N 0 ∈ N such that Then W 0 = X N0+1 ∈ D M 5 , given in (9), and then the sequence {W n } converges to X.
In this way, we have obtained the following local convergence result for the iterative scheme given in (11).

Theorem 8. Let X be a Hermitian solution of (2) and A − D X is stable.
Then, the sequence {W n } generated by the iterative scheme (11) converges to X from every matrix X 0 ∈ B( X, S), with β D S < 1 3 . As we can see, from this result, the iterative scheme (11) has the same accessibility as Newton's method, maintaining, except for the first N 0 iterations, the same computational efficiency as the iterative scheme M5.

A Second Improvement for M5 Method
As we have already seen in the case of the iterative scheme M5 (5), highorder iterative schemes are known to have a small region of accessibility associated with them. Therefore, locating starting points for them is a difficult problem to solve. For this reason, in this new improvement, we maintain the idea of improving the accessibility of the M5 method in the first stage of the predictor-corrector iterative scheme. But, in this case, we also intend to reduce the operational cost. The first improvement of the M5 method that we have established in the previous section, we have considered Newton's method as a predictor method. This iterative scheme improved accessibility but for its application it requires the resolution of a Lyapunov equation in each of the N 0 iterations that we must carry out. Well, we propose a second improvement of the M5 method in which we are going to consider the wellknown modified Newton method [15]. This method also improves accessibility and only needs to solve a Lyapunov equation, with different independent terms, in the N 0 iterations that we must carry out, which significantly reduces the operational cost of the predictor-corrector iterative scheme. Thus, we consider the modified Newton method as predictor iterative scheme: To obtain the (N 0 + 1)-th step of the modified Newton method (13), is equivalent to solve the Lyapunov equation, [18]: where P D : H → H, is the operator defined by P D (X) := XDX. Therefore, the modified Newton method given in (13) has a low operational cost, the cost of solving a Lyapunov equation (14) with N 0 different independent terms. Hence, we show that it has a good accessibility region, although smaller than that of Newton's method. But, it has a problem, its linear convergence, whose influence will depend on the number of iterations to perform in the first stage of the predictor-corrector method.
Next, we study the accessibility region of the method.
Theorem 9. Let X be a Hermitian solution of (2) and A − D X is stable. If X 0 ∈ B( X, S) with β D S < 0.2, then the sequence the sequence {X k }, given by the modified Newton method (13), converges to X.
Proof. Notice that, from Lemma 3, for n 1 if X n ∈ B( X, S), with β||D||S < 0.2, then there exists the operator [R (X n )] −1 . Now, from the algorithm of the modified Newton method (13), we have: Next, taking into account the Taylor series, we have: Thus, the following decomposition is verified: Thus, taking norms, if X n ∈ B( X, S), we obtain Hence, if we consider n = 0, we obtain that and as 3Δ(S) < 1, since β||D||S < 0.2, then X 1 ∈ B( X, S) with X 1 − X < X 0 − X . To continue, by applying an inductive procedure, it follows that X n+1 ∈ B( X, S) and X n+1 − X < X n − X for all n 1. Therefore, { X n+1 − X } is a strictly decreasing sequence of positive real numbers and then the result is proved.
From the previous Theorem, we obtain that the domain of parameters for the modified Newton method is given by Therefore, the domain of parameters for the modified Newton method is smaller than the domain of parameters for Newton's method given in (10), but it is greater than that of the M5 method. In this case, Therefore, the modified Newton method is an interesting iterative scheme to consider as a predictor scheme. It has better accessibility than the M5 method and also has a lower operating cost than the Newton method.

A Predictor-Corrector Iterative Scheme
From the previous comments, to obtain an efficient iterative scheme with a suitable region of accessibility, we consider the predictor-corrector iterative scheme given by Given an initial guess X 0 ∈ H, whose expression, according to the practical application by solving Lyapunov equations, is given as follows: Given an initial guess X 0 ∈ H, Next, we study the local convergence of this iterative scheme (15). If X 0 ∈ B( X, S) satisfies condition β||D||S < 0.2, and W 0 = X N0 ∈ B( X, S) satisfies the condition β||D|| S < 0.14, in this situation, the sequence {W n }, given by the M5 method, converges to X a solution of symmetric ARE (2). To continue, we show how to do this procedure.
Theorem 10. Let X be a Hermitian solution of (2) and A − D X is stable. Then, the sequence {W n } generated by the iterative scheme (15) converges to X from every matrix X 0 ∈ B( X, S), with β D S < 0.2.

Numerical Experiments
We show a numerical experiment related to the control for a tubular ammonia reactor, which can be found in the benchmark collection for CARE [22]. We apply the Newton method and the constructed predictor-corrector methods, to approximate a solution of equation (2) related with a ninth-order discrete state-space model of a tubular ammonia reactor, where the matrices A and B are the following: Numerical experiment shows the high precision and accuracy of the constructed predictor-corrector methods compared to the well-known Newton's method. The Newton method is the iterative process commonly used to approximate solutions of nonlinear equations, in particular Eq. (2). In tables 1 and 2, is showed the number of iterations, errors X * − X k F , with stopping criteria X * − X k F < 10 −40 , and residuals R(X k ) F , of the iterative processes considered.
First, we consider the starting matrix This starting matrix is found in the accessibility domain of the Newton Method Thus, performing a step with the Newton method, we already obtain a matrix that is found in the accessibility domain of the M5 method. Therefore, we obtain N 0 = 1 for predictor-corrector method (12). In Table 1, we can verify that the predictor-corrector method (12) obtains greater precision in a smaller number of iterations than Newton's method. In addition, we see that it uses a smaller number of operations. Therefore, we can say that the predictor-corrector method (12) improves Newton's method.
Second, we consider X 0 = 3I 9 /2, where I 9 is the identity matrix of size 9. This matrix is found in the accessibility domain of the modified Newton method. And, performing two step with the modified Newton method, we already obtain a matrix that is found in the accessibility domain of the M5 method. Therefore, we obtain N 0 = 2 for the predictor-corrector method (16). In Table 2, we can observe that with fewer iterations, the predictorcorrector method (16)  number of operations. Therefore, we can say that the predictor-corrector method (16) improves Newton's method. Finally, we can state that the constructed predictor-corrector methods, (12) and (16), improve Newton's method for approximating a solution of the equation (2).
Author contributions This article is a joint research between both authors.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. The research has been partially supported by the project MTM2018-095896-B-C21 of the Spanish Ministry of Science.

Conflict of interest
The authors declare no competing of interest.

Ethical approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.