Fast iterative solutions of Riccati and Lyapunov equations

1 Introduction

The discrete time Riccati equation arises by implementing the discrete time Kalman filter, which is the most famous estimation algorithm associated with time invariant models of the form [1]:

for k ≥ 0 (discrete time).

In these models, x ( k ) is the n dimensional state and z ( k ) is the m dimensional measurement. The model describes the relation between two successive states through the n × n transition matrix F and the relation between the state and the measurement through the m × n output matrix H . In addition, w ( k ) and v ( k ) are the state and measurement noises, which are zero mean Gaussian processes with known covariances Q and R , respectively. It is worth to note that Q and R are symmetric and nonnegative definite of dimensions n × n and m × m , respectively. So, the model parameters F , H , Q and R are time invariant (constant) and known.

Kalman filter produces iteratively the state prediction with the associated prediction error covariance as well as the state estimation with the associated estimation error covariance. The equations of the discrete time Kalman filter result in the discrete time Riccati equation:

Riccati equation relates two successive values of the prediction error covariance P ( k |k − 1 ) , which is symmetric and nonnegative definite of dimensions n × n .

It is well known [1] that if the system is asymptotically stable, i.e., all eigenvalues of the transition matrix F lie inside the unit circle, then for any nonnegative symmetric initial condition P 0 = P ( 0 | − 1 ) , there exists a unique steady state or limiting solution P of the Riccati equation:

This limiting solution satisfies the algebraic Riccati equation (or steady-state Riccati equation):

It is worth to note that inverse of the matrix HP ( k |k − 1 ) H T + R in (3) exists when R is positive definite, which means that no measurement is accurate, or when the initial condition P 0 is positive definite, which seems not to be a problem, since we are able to choose a convenient initial condition.

The nonsingularity of Q and R ensures the nonsingularity of P ( k |k − 1 ) , which then it is positive definite. In this (reasonable and general) case, by using the matrix inversion lemma,^[1] the Riccati equation (3) can be transformed to the transformed Riccati equation:

where

The limiting solution P of the transformed Riccati equation (6) satisfies the transformed algebraic Riccati equation:

In the infinite measurement noise case, the discrete time Lyapunov equation is derived:

with a limiting solution P , that satisfies the algebraic Lyapunov equation:

The importance of the Riccati and Lyapunov equations is undoubtable: the Riccati equation plays a very important role in various problems of stochastic filtering, statistics, ladder networks, and dynamic programming [1,2], and has many applications in the process of obtaining optimal control and determining system stability [3,4]; Lyapunov equations play a very important role in the stability theory of discrete systems [1,5].

Significant bibliography exists concerning iterative as well as noniterative solutions of the Riccati and Lyapunov equations [1,2,6,7,8,9,10,11,12,13]. The existence of the unique solution of the Riccati equation and of the Lyapunov equation requires the knowledge of the eigenvalues of the transition matrix. The classical iterative solution of the Riccati equation implements the iterative Riccati equation (3) or the iterative transformed Riccati equation (6) till the steady-state solution is reached. Similarly, the classical iterative solution of the Lyapunov equation implements the iterative Lyapunov equation (7) till the steady-state solution is reached.

The aim of this article is to study the optimal algorithm for solving the discrete Riccati and Lyapunov equations, in the case where the transition matrix is diagonalizable with real eigenvalues. The required condition appears in various application areas, see the eigenvalues of the transition matrix in the radar tracking system [14], in the eye movement prediction model [15], and in the mobile position tracking model [16,17]; the state and parameter estimation has been derived using a transition matrix, which is a diagonal itself [18]. The novelty of this article concerns the exploitation of the knowledge of the eigenvalues and eigenvectors of the transition matrix, which leads to the development of the new algorithms for iterative solving of the Riccati and Lyapunov equations, which are superior to classical algorithms, since the proposed algorithms are in general faster than the classical ones.

This article is organized as follows: In Section 2, new iterative solutions of the Riccati and Lyapunov equations are derived by the eigenvalues and eigenvectors of the transition matrix. In Section 3, the proposed algorithms are compared to the classical ones with respect to their computational burdens. Section 4 summarizes the conclusions.

2 New iterative solution algorithms for Riccati equation and Lyapunov equation

Recall that the model parameters F , H , Q and R are time invariant (constant) and known. The Riccati equation and the Lyapunov equation have an unique symmetric nonnegative definite solution P if and only if | λ i | < 1 for all i = 1 , … , n , where λ i , i = 1 … n are the eigenvalues of F .

In the following, we consider that the n × n transition matrix F has real eigenvalues with the algebraic multiplicity of every eigenvalue identified by its geometric; thus, the associated eigenvectors are linearly independent vectors.

Consider the diagonalization formula of the transition matrix F given by

where L is the n × n diagonal matrix containing the eigenvalues of matrix F and the n × n matrix; W consists of columns of the associated linearly independent eigenvectors of F.

It is worth noting that equation (11) is essential for the derivation of the new iterative algorithms for solving the discrete Riccati and Lyapunov equations, since they use the eigenvalues and eigenvectors of the transition matrix instead of the transition matrix itself, reducing the computations due to the similarity of the matrices F and L .

Concerning the Riccati equation, by multiplying the left-hand side of equation (3) by W − 1 and the right-hand side by W − T and using (11) and the property of the diagonal matrix L, L = L T , as well as the equalities W T W − T = W − T W T = I , where I denotes the identity matrix, we derive

Setting in the latter equality

arises the modified Riccati equation.

It is important to note that the Riccati equation (3) with parameters F , H , Q and R has been modified to equation (15), which is a Riccati equation with well-defined parameters L , H ∼ , Q ∼ and R .

Then, the limiting solution of (15) is

and satisfies the modified algebraic Riccati equation:

It is clear that the solution of the Riccati equation (3) can be computed by the solution of the modified Riccati Equation (15) due to (16) and (12):

Concerning the transformed Riccati equation, by multiplying the left-hand side of equation (6) by W − 1 and the right-hand side by W − T , using (11) and the property of the diagonal matrix L, L = L T , as well as the equalities W T W − T = W − T W T = I , we derive

Setting in the latter equality

Q arises the modified transformed Riccati equation:

It is important to note that the transformed Riccati equation (6) with parameters F , H , Q and R has been modified to Equation (21), which is a transformed Riccati equation with well-defined parameters L , H ∼ , Q ∼ and R .

Then, using (16)

arises the limiting solution of (21), which satisfies the algebraic Riccati equation:

It is clear that the solution of the transformed Riccati equation (6) can be computed by the solution of the Riccati equation (22) due to (16) and (19):

Concerning the Lyapunov equation, by multiplying the left-hand side of (9) by W − 1 and the right-hand side by W − T , using (11) and the property of the diagonal matrix L, L = L T , as well as the equalities W T W − T = W − T W T = I , we derive

Setting in the latter equality

arises

It is important to note that the Lyapunov equation (9) with parameters F and Q has been modified to equation (25), which is a Lyapunov equation with well-defined parameters L and Q ∼ .

Then, using (16)

arises the limiting solution of (25), which satisfies the modified algebraic Lyapunov equation:

It is clear that the solution of the Lyapunov equation (9) can be computed by the solution of the modified Lyapunov equation (25) due to (16) and (24):

Table 1 summarizes the classical and the derived iterative algorithms for solving the Riccati and Lyapunov equations.

3 Comparison of classical and proposed algorithms

It is established that the proposed iterative algorithms for solving the Riccati and Lyapunov equations have been derived by the classical iterative algorithms. Thus, the classical and the derived algorithms are equivalent with respect to their behavior since they calculate theoretically the same solutions executing the same number of iterations s .

We are going to compare the classical and the derived algorithms with respect to their calculation burdens. Scalar operations are involved in matrix manipulation operations, which are needed for the implementation of the algorithms. Table 2 summarizes the calculation burdens of needed matrix operations. The computation of the eigenvalues is achieved by solving the characteristic equation using the Newton-Raphson method, Laguerre’s method, Regula-Falsi method or Bernoulli’s method, because each method requires different computations for solving the characteristic equation. In the following, we consider that all the algorithms apply the same method; thus the computation of the eigenvalues has the same calculation burden, which is denoted as CB L . The eigenvectors computation is achieved using Gauss elimination and backward substitutions. In fact, Gauss elimination requires 1 6 ( 4 n 3 + 9 n 2 − 7 n ) scalar operations; analytically, 1 2 ( n 2 + n ) divisions, 1 6 ( 2 n 3 + 3 n 2 − 5 n ) multiplications, and 1 6 ( 2 n 3 + 3 n 2 − 5 n ) subtractions; the details are given in [19]. Backward substitutions require n 2 scalar operations: analytically: ∑ i = 1 n − 1 i = ( n − 1 ) n 2 multiplications and ∑ i = 1 n − 1 i = ( n − 1 ) n 2 and n divisions. Thus, the eigenvalues computation requires 1 6 ( 4 n 3 + 15 n 2 − 7 n ) scalar operations. The details for the matrix operations are given in [20].

Table 3 summarizes the calculation burdens of the classical and the proposed algorithms, for the general multidimensional case, where n ≥ 2 and m ≥ 2 . The details are given in the Appendix. In this table, the calculation burden CB L for the computation of the eigenvalues is the same for all algorithms, since we consider that the computation is achieved by solving the characteristic equation using the same method.

Concerning the Riccati equation solution algorithms, from Table 3, we can write:

Note that in the case

since 18 n 2 − 15 n − 3 > 0 , for every n ≥ 2 , equation (28) yields

The aforementioned notation leads to the conclusion that the proposed algorithm is faster than the classical algorithm for a small number of iterations, especially when the state and measurement dimensions are large.

Figure 1 depicts the faster algorithm with respect to the state dimension n , the measurement dimension m , and the number of iterations s the algorithms execute. The proposed algorithm is faster than the classical algorithm in the region under the curves and is slower than the classical algorithm in the region above the curves, confirming the aforementioned note.

Figure 1

Comparison of the Riccati equation solution algorithms.

Concerning the transformed Riccati equation solution algorithms, from Table 3, we obtain:

Note that in the case

since 18 n 2 − 15 n − 3 > 0 , for every n ≥ 2 , equation (29) yields

The aforementioned notation leads to the conclusion that the proposed algorithm is faster than the classical algorithm for a small number of iterations. In fact, the proposed algorithm is always faster than the classical algorithm, when s > 7 . In addition, the proposed algorithm outperforms the classical one, when the number of iterations is greater than 4 as the state dimension increases.

Figure 2 depicts the faster algorithm with respect to the state dimension n and the number of iterations s the algorithms execute. The proposed algorithm is faster than the classical algorithm in the region under the curve and is slower than the classical algorithm in the region above the curve, confirming the aforementioned note.

Figure 2

Comparison of transformed Riccati equation solution algorithms.

Concerning the Lyapunov equation solution algorithms, from Table 3, we obtain:

Note that in the case

since 18 n 2 − 15 n − 3 > 0 , for every n ≥ 2 , equation (30) yields

The aforementioned notation leads to the conclusion that the proposed algorithm is faster than the classical algorithm for a small number of iterations. In fact, the proposed algorithm is always faster than the classical algorithm, when s > 5 . In addition, the proposed algorithm outperforms the classical one, when the number of iterations is greater than 3, as the state dimension increases.

Figure 3 depicts the faster algorithm with respect to the state dimension n and the number of iterations s the algorithms execute. The proposed algorithm is faster than the classical algorithm in the region under the curve and is slower than the classical algorithm in the region above the curve, confirming the aforementioned note.

Figure 3

Comparison of Lyapunov equation solution algorithms.

4 Conclusions

Riccati and Lyapunov equations are very important in many fields of science. The existence of their solutions involves the knowledge of the eigenvalues of the transition matrix. We have developed new fast iterative algorithms for solving the Riccati and Lyapunov equations. The proposed algorithms take advantage of the knowledge of the eigenvalues and eigenvectors of the transition matrix. The proposed algorithms hold in the case where the transition matrix is diagonalizable in the form = WL W − 1 , where the diagonal matrix L contains the real eigenvalues of F and the matrix W consists of columns containing the linearly independent eigenvectors of F. It is shown that the proposed algorithms may be faster than the classical ones depending on the state dimension, the measurement dimension, and the number of iterations executed, even for a small number of iterations. The proposed algorithms are faster than the classical ones for such number of iterations that decreases as the state dimension increases. This is rational because the operations that involve the eigenvalue matrix L are of the order of O ( n 2 ) , while operations that involve the eigenvalue matrix F are of the order of O ( n 3 ) .

More specifically,

1. Concerning the Riccati equation, the proposed solution is faster than the classical one for a small number of iterations, especially when the state and measurement dimensions are large.

2. Concerning the transformed Riccati equation, the proposed solution is always faster than the classical one when the number of iterations is greater than 7; the proposed algorithm outperforms the classical one when the number of iterations is greater than 4, as the state dimension increases.

3. Concerning the Lyapunov equation, the proposed solution is always faster than the classical one when the number of iterations is greater than 5; the proposed algorithm outperforms the classical one when the number of iterations is greater than 3, as the state dimension increases.

It is evident that the proposed iterative algorithms outperform the classical ones even for a small number of iterations.