Fast iterative solutions of Riccati and Lyapunov equations

• Nicholas Assimakis and Maria Adam
From the journal Open Mathematics

Abstract

In this article, new iterative algorithms for solving the discrete Riccati and Lyapunov equations are derived in the case where the transition matrix is diagonalizable with real eigenvalues. It is shown that the proposed iterative algorithms are faster than the classical ones, even for a small number of iterations.

MSC 2010: 15A24; 15B57; 93C55

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]:

(1) x ( k + 1 ) = Fx ( k ) + w ( k ) ,

(2) z ( k ) = Hx ( k ) + v ( k ) ,

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:

(3) P ( k + 1 | k ) = Q + FP ( k | k 1 ) F T FP ( k | k 1 ) H T [ HP ( k | k 1 ) H T + R ] 1 HP ( k | k 1 ) F T .

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:

(4) P = lim k P ( k | k 1 ) .

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

(5) P = Q + FP F T FP H T [ HP H T + R ] 1 HP F T .

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:

(6) P ( k + 1 | k ) = Q + F [ P 1 ( k | k 1 ) + B ] 1 F T ,

where

(7) B = H T R 1 H .

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

(8) P = Q + F [ P 1 + B ] 1 F T .

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

(9) P ( k + 1 | k ) = Q + FP ( k | k 1 ) F T ,

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

(10) P = Q + FP F T .

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

(11) F = WL W 1 ,

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

W 1 P ( k + 1 | k ) W T = W 1 Q W T + L W 1 P ( k | k 1 ) W T L L W 1 P ( k | k 1 ) W T W T H T

[ H WW 1 P ( k | k 1 ) W T W T H T + R ] 1

H WW 1 P ( k | k 1 ) W T L .

Setting in the latter equality

(12) P ( k | k 1 ) = W 1 P ( k | k 1 ) W T ,

(13) Q = W 1 Q W T ,

(14) H = HW

arises the modified Riccati equation.

(15) P ( k + 1 | k ) = Q + L P ( k | k 1 ) L L P ( k | k 1 ) H T [ H P ( k | k 1 ) H T + R ] 1 H P ( k | k 1 ) L .

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

(16) P = lim k P ( k + 1 | k ) = lim k P ( k | k 1 )

and satisfies the modified algebraic Riccati equation:

(17) P = Q + L P L L P H T [ H P H T + R ] 1 H P L .

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):

(18) P = W P W T .

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

W 1 P ( k + 1 | k ) W T = W 1 Q W T + L W 1 [ P 1 ( k | k 1 ) + B ] 1 W T L T W T W T = W 1 Q W T + L [ W T P 1 ( k | k 1 ) W + W T BW ] 1 L .

Setting in the latter equality

(19) P ( k | k 1 ) = W 1 P ( k | k 1 ) W T ,

Q = W 1 Q W T ,

(20) B = W T BW

Q arises the modified transformed Riccati equation:

(21) P ( k + 1 | k ) = Q + L [ P 1 ( k | k 1 ) + B ] 1 L .

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)

P = lim k P ( k |k 1 ) = lim k P ( k + 1 |k )

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

(22) P = Q + L [ P 1 + B ] 1 L .

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):

(23) P = W P W T .

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

P ( k + 1 | k ) W T = W 1 Q W T + L W 1 P ( k | k 1 ) W T L .

Setting in the latter equality

(24) P ( k | k 1 ) = W 1 P ( k | k 1 ) W T ,

Q = W 1 Q W T

arises

(25) P ( k + 1 | k ) = Q + L P ( k | k 1 ) L .

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)

P = lim k P ( k |k 1 ) = lim k P ( k + 1 |k )

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

(26) P = Q + L P L .

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):

(27) P = W P W T .

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

Table 1

Riccati and Lyapunov equations iterative solution algorithms

Equation Algorithm
Riccati equation Classical input: F , H , Q , R initialization: L = eigenvalues of F iteration P ( k + 1 |k ) = Q + FP ( k |k 1 ) F T FP ( k |k 1 ) H T [ HP ( k |k 1 ) H T + R ] 1 HP ( k |k 1 ) F T convergence: P = lim k P ( k |k 1 ) output: P
Proposed input: F , H , Q , R initialization: L = eigenvalues of F , W = eigenvectors of F , Q ̃ = W 1 Q W T , H ̃ = HW iteration P ̃ ( k + 1 | k ) = Q ̃ + L P ̃ ( k | k 1 ) L L P ̃ ( k | k 1 ) H ̃ T [ H ̃ P ̃ ( k | k 1 ) H ̃ T + R ] 1 H ̃ P ̃ ( k | k 1 ) L convergence: P ̃ = lim k P ̃ ( k |k 1 ) finalization: P = W P ̃ W T output: P
Transformed Riccati equation Classical input: F , H , Q , R initialization: L = eigenvalues of F , B = H T R 1 H iteration P ( k + 1 | k ) = Q + F [ P 1 ( k | k 1 ) + B ] 1 F T convergence: P = lim k P ( k |k 1 ) output: P
Proposed input: F , H , Q , R initialization: L = eigenvalues of F , W = eigenvectors of F , B = H T R 1 H , Q ̃ = W 1 Q W T , B ̃ = W T BW iteration P ̃ ( k + 1 | k ) = Q ̃ + L [ P ̃ 1 ( k | k 1 ) + B ̃ ] 1 L convergence: P ̃ = lim k P ̃ ( k |k 1 ) finalization: P = W P ̃ W T output: P
Lyapunov equation Classical input: F , Q initialization: L = eigenvalues of F iteration P ( k + 1 | k ) = Q + FP ( k |k 1 ) F T convergence: P = lim k P ( k |k 1 ) output: P
Proposed input: F , Q initialization: L = eigenvalues of F , W = eigenvectors of F , Q ̃ = W 1 Q W T iteration P ̃ ( k + 1 | k ) = Q ̃ + L P ̃ ( k | k 1 ) L convergence: P ̃ = lim k P ̃ ( k |k 1 ) finalization: P = W P ̃ W T output: P

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 2

Calculation burdens of matrix operations

Matrix operation Matrix dimensions Calculation burden
L = eigenvalues of F n × n CB L
W = eigenvectors of F n × n 1 6 ( 4 n 3 + 15 n 2 7 n )
C = A + B ( n × m ) + ( n × m ) nm
(1) S = A + B ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n
C = A B ( n × m ) ( m × l ) 2 nm l n l
(1) S = A B ( n × m ) ( m × n ) n 2 m + nm 1 2 n 2 1 2 n
(2) A D ( n × n ) ( n × n ) n 2
(2) D A ( n × n ) ( n × n ) n 2
(3) B = A 1 n × n 1 6 ( 16 n 3 3 n 2 n )

(1) S is a symmetric matrix.

(2) D is a diagonal matrix.

(3) For the general multi-dimensional case, where n 2 .

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.

Table 3

Calculation burdens of algorithms

Equation Algorithm Calculation burden
Riccati equation Classical CB CRE = CB L + 1 6 ( 16 m 3 3 m 2 m ) + 3 n 3 + 3 n 2 m + 3 n m 2 s
Proposed CB PRE = CB L + 1 6 ( 56 n 3 + 6 n 2 14 n ) + 2 n 2 m nm + 1 6 ( 16 m 3 3 m 2 m ) + 1 2 ( 5 n 2 + n ) + 3 n 2 m + 3 n m 2 s
Transformed Riccati equation Classical CB CTRE = CB L + 1 6 ( 16 m 3 3 m 2 m ) + 2 n m 2 + n 2 m 1 2 n 2 1 2 n + 1 6 ( 50 n 3 3 n 2 + n ) s
Proposed CB PTRE = CB L + 1 6 ( 16 m 3 3 m 2 m ) + 2 n m 2 + n 2 m + 1 6 ( 74 n 3 20 n ) + 1 6 ( 32 n 3 + 12 n 2 + 4 n ) s
Lyapunov equation Classical CB CLE = CB L + 3 n 3 s
Proposed CB PLE = CB L + 1 6 ( 56 n 3 + 6 n 2 14 n ) + 1 2 ( 5 n 2 + n ) s

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

(28) CB CRE CB PRE = n 6 { ( 18 n 2 15 n 3 ) s ( 56 n 2 + 6 n 14 ) ( 12 n 6 ) m } .

Note that in the case

s > 56 n 2 + 6 n 14 18 n 2 15 n 3 + 12 n 6 18 n 2 15 n 3 m = 3 + 2 n 2 + 51 n 5 18 n 2 15 n 3 + 12 n 6 18 n 2 15 n 3 m ,

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

CB CRE CB PRE > 0 .

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:

(29) CB CTRE CB PTRE = n 6 { ( 18 n 2 15 n 3 ) s ( 74 n 2 + 3 n 17 ) } .

Note that in the case

s > 74 n 2 + 3 n 17 18 n 2 15 n 3 = 4 + 2 n 2 + 63 n 5 18 n 2 15 n 3 ,

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

CB CTRE CB PTRE > 0 .

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:

(30) CB CLE CB PLE = n 6 { ( 18 n 2 15 n 3 ) s ( 56 n 2 + 6 n 14 ) } .

Note that in the case

s > 56 n 2 + 6 n 14 18 n 2 15 n 3 = 3 + 2 n 2 + 51 n 5 18 n 2 15 n 3 ,

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

CB CLE CB PLE > 0

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.

Acknowledgements

1. Funding information: The authors state that no funding was involved.

2. Author contributions: The authors applied the SDC approach for the sequence of authors. NA: conceptualization and performed the simulation. MA: performed the formal analysis and prepared the manuscript.

3. Conflict of interest: The authors state no conflict of interest.

4. Data availability statement: All data generated or analyzed during this study are included in this published article.

Appendix

The calculation burdens of the classical and the proposed algorithms for solving the Riccati and Lyapunov equations, for the general multi-dimensional case (where n 2 , m 2 ), are analytically presented.

Riccati equation – classical iterative algorithm

 Initialization L = eigenvalues of F CB L Iteration ( s iterations) M 1 = HP ( k | k − 1 ) ( m × n ) ⋅ ( n × n ) 2 n 2 m − nm M 2 = M 1 H T ( m × n ) ⋅ ( n × m ) n m 2 + nm − 1 2 m 2 − 1 2 m M 3 = M 2 + R ( m × m ) + ( m × m ) 1 2 m 2 + 1 2 m M 4 = M 3 − 1 ( m × m ) 1 6 ( 16 m 3 − 3 m 2 − m ) M 5 = M 4 M 1 ( m × m ) ⋅ ( m × n ) 2 n m 2 − nm M 6 = M 1 T M 5 ( n × m ) ⋅ ( m × n ) n 2 m + nm − 1 2 n 2 − 1 2 n M 7 = P ( k |k − 1 ) − M 6 ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n M 8 = F M 7 ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 M 9 = M 8 F T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n P ( k + 1 |k ) = Q + M 9 ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n CB CRE = CB L + 1 6 ( 16 m 3 − 3 m 2 − m ) + 3 n 3 + 3 n 2 m + 3 n m 2 s

Riccati equation – proposed iterative algorithm

 Initialization L = eigenvalues of F CB L W = eigenvectors of F 1 6 ( 4 n 3 + 15 n 2 − 7 n ) W − 1 ( n × n ) 1 6 ( 16 n 3 − 3 n 2 − n ) W − 1 Q ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 Q ∼ = W − 1 Q W − T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n H ∼ = HW ( m × n ) ⋅ ( n × n ) 2 n 2 m − nm Iteration ( s iterations) M 1 = H ∼ P ∼ ( k | k − 1 ) ( m × n ) ⋅ ( n × n ) 2 n 2 m − nm M 2 = M 1 H ∼ T ( m × n ) ⋅ ( n × m ) n m 2 + nm − 1 2 m 2 − 1 2 m M 3 = M 2 + R ( m × m ) + ( m × m ) 1 2 m 2 + 1 2 m M 4 = M 3 − 1 ( m × m ) 1 6 ( 16 m 3 − 3 m 2 − m ) M 5 = M 4 M 1 ( m × m ) ⋅ ( m × n ) 2 n m 2 − nm M 6 = M 1 T M 5 ( n × m ) ⋅ ( m × n ) n 2 m + nm − 1 2 n 2 − 1 2 n M 7 = P ∼ ( k |k − 1 ) − M 6 ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n M 8 = L M 7 ( n × n ) ⋅ ( n × n ) n 2 M 9 = M 8 L ( n × n ) ⋅ ( n × n ) n 2 P ∼ ( k + 1 |k ) = Q ∼ + M 9 ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n Finalization W P ∼ ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 P = W P ∼ W T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n CB PRE = CB L + 1 6 ( 56 n 3 + 6 n 2 − 14 n ) + 2 n 2 m − nm + 1 6 ( 16 m 3 − 3 m 2 − m ) + 1 2 ( 5 n 2 + n ) + 3 n 2 m + 3 n m 2 s

Transformed Riccati equation – classical iterative algorithm

 Initialization L = eigenvalues of F CB L R − 1 ( m × m ) 1 6 ( 16 m 3 − 3 m 2 − m ) H T R − 1 ( n × m ) ⋅ ( m × m ) 2 n m 2 − nm B = H T R − 1 H ( n × m ) ⋅ ( m × n ) n 2 m + nm − 1 2 n 2 − 1 2 n Iteration ( s iterations) P − 1 ( k | k − 1 ) ( n × n ) 1 6 ( 16 n 3 − 3 n 2 − n ) P − 1 ( k | k − 1 ) + B ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n [ P − 1 ( k | k − 1 ) + B ] − 1 ( n × n ) 1 6 ( 16 n 3 − 3 n 2 − n ) F [ P − 1 ( k | k − 1 ) + B ] − 1 ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 F [ P − 1 ( k | k − 1 ) + B ] − 1 F T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n P ( k + 1 | k ) = Q + F [ P − 1 ( k | k − 1 ) + B ] − 1 F T ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n CB CTRE = CB L + 1 6 ( 16 m 3 − 3 m 2 − m ) + 2 n m 2 + n 2 m − 1 2 n 2 − 1 2 n + 1 6 ( 50 n 3 − 3 n 2 + n ) s

Transformed Riccati equation – proposed iterative algorithm

 Initialization L = eigenvalues of F CB L W = eigenvectors of F 1 6 ( 4 n 3 + 15 n 2 − 7 n ) W − 1 ( n × n ) 1 6 ( 16 n 3 − 3 n 2 − n ) W − 1 Q ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 Q ∼ = W − 1 Q W − T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n R − 1 ( m × m ) 1 6 ( 16 m 3 − 3 m 2 − m ) H T R − 1 ( n × m ) ⋅ ( m × m ) 2 n m 2 − nm B = H T R − 1 H ( n × m ) ⋅ ( m × n ) n 2 m + nm − 1 2 n 2 − 1 2 n W T B ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 B ∼ = W T BW ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n Iteration ( s iterations) P ∼ − 1 ( k | k − 1 ) ( n × n ) 1 6 ( 16 n 3 − 3 n 2 − n ) P ∼ − 1 ( k | k − 1 ) + B ∼ ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n [ P ∼ − 1 ( k | k − 1 ) + B ∼ ] − 1 ( n × n ) 1 6 ( 16 n 3 − 3 n 2 − n ) L [ P ∼ − 1 ( k | k − 1 ) + B ∼ ] − 1 ( n × n ) ⋅ ( n × n ) n 2 L [ P ∼ − 1 ( k | k − 1 ) + B ∼ ] − 1 L ( n × n ) ⋅ ( n × n ) n 2 P ∼ ( k + 1 | k ) = Q ∼ + L [ P ∼ − 1 ( k | k − 1 ) + B ∼ ] − 1 L ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n Finalization W P ∼ ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 P = W P ∼ W T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n CB PTRE = CB L + 1 6 ( 16 m 3 − 3 m 2 − m ) + 2 n m 2 + n 2 m + 1 6 ( 74 n 3 − 20 n ) + 1 6 ( 32 n 3 + 12 n 2 + 4 n ) s

Lyapunov equation – classical iterative algorithm

 Initialization L = eigenvalues of F CB L iteration ( s iterations) FP ( k | k − 1 ) ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 FP ( k | k − 1 ) F T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n P ( k + 1 | k ) = Q + FP ( k | k − 1 ) F T ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n CB CLE = CB L + 3 n 3 s

Lyapunov equation – proposed iterative algorithm

 Initialization L = eigenvalues of F CB L W = eigenvectors of F 1 6 ( 4 n 3 + 15 n 2 − 7 n ) W − 1 ( n × n ) 1 6 ( 16 n 3 − 3 n 2 − n ) W − 1 Q ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 Q ∼ = W − 1 Q W − T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n Iteration ( s iterations) L P ∼ ( k + 1 | k ) ( n × n ) ⋅ ( n × n ) n 2 L P ∼ ( k + 1 | k ) L ( n × n ) ⋅ ( n × n ) n 2 P ∼ ( k + 1 | k ) = Q ∼ + L P ∼ ( k | k − 1 ) L ( n × n ) + ( n × n ) 1 2 n 2 + 1 2 n Finalization W P ∼ ( n × n ) ⋅ ( n × n ) 2 n 3 − n 2 P = W P ∼ W T ( n × n ) ⋅ ( n × n ) n 3 + 1 2 n 2 − 1 2 n CB PLE = CB L + 1 6 ( 56 n 3 + 6 n 2 − 14 n ) + 1 2 ( 5 n 2 + n ) s

References

[1] B. D. O. Anderson and J. B. Moore, Optimal Filtering, Dover Publications, New York, 2005.Search in Google Scholar

[2] N. Assimakis and M. Adam, Fast doubling algorithm for the solution of the riccati equation using cyclic reduction method, 2020 International Conference on Mathematics and Computers in Science and Engineering (MACISE), 2020, pp. 1–5, https://doi.org/10.1109/MACISE49704.2020.00007.10.1109/MACISE49704.2020.00007Search in Google Scholar

[3] D. S. Bernstein, Matrix Mathematics: Theory Facts and Formulas with Application to Linear Systems Theory, Princeton University Press, Princeton, NJ, 2005.Search in Google Scholar

[4] C. Kojima, K. Takaba, O. Kaneko, and P. Rapisarda, A characterization of solutions of the discrete-time algebraic riccati equation based on quadratic difference forms, Linear Algebra Appl. 416 (2006), 1060–1082.10.1016/j.laa.2005.11.027Search in Google Scholar

[5] A. Nakhmani, Modern Control: State-Space Analysis and Design Methods, McGraw Hill, New York, 2020.Search in Google Scholar

[6] N. D. Assimakis, D. G. Lainiotis, S. K. Katsikas, and F. L. Sanida, A survey of recursive algorithms for the solution of the discrete time Riccati equation, Nonlinear Anal. (1997), no. 4, 2409–2420.10.1016/S0362-546X(97)00062-XSearch in Google Scholar

[7] D. G. Lainiotis, N. D. Assimakis, and S. K. Katsikas, A new computationally effective algorithm for solving the discrete Riccati equation, J. Math. Anal. Appl. 186 (1994), no. 3, 868–895.10.1006/jmaa.1994.1338Search in Google Scholar

[8] D. G. Lainiotis, N. D. Assimakis, and S. K. Katsikas, Fast and numerically robust recursive algorithms for solving the discrete time Riccati equation: The case of nonsingular plant noise covariance matrix, Neural, Parallel Sci. Comput. 3 (1995), 565–584.Search in Google Scholar

[9] N. Assimakis, S. Roulis, and D. Lainiotis, Recursive solutions of the discrete time riccati equation, Neural, Parallel Sci. Comput. 11 (2003), 343–350.Search in Google Scholar

[10] V. Dragan, The linear quadratic optimization problem for a class of singularly stochastic systems, Int. J. Innov. Comput. Inf. Control 1 (2005), 53–63.Search in Google Scholar

[11] N. Komaroff, Iterative matrix bounds and computational solutions to the discrete algebraic Riccati equation, IEEE Trans. Autom. Control 39 (1994), 1676–1678.10.1109/9.310049Search in Google Scholar

[12] J. Zhang and J. Liu, New upper and lower bounds, the iteration algorithm for the solution of the discrete algebraic Riccati equation, Adv. Differ. Equ. 313 (2015), 1–17, https://doi.org/10.1186/s13662-015-0649-6.10.1186/s13662-015-0649-6Search in Google Scholar

[13] N. Assimakis and M. Adam, Closed form solutions of Lyapunov equations using the vech and veck operators, WSEAS Trans. Math. 20 (2021), 276–282, https://doi.org/10.37394/23206.2021.20.28.10.37394/23206.2021.20.28Search in Google Scholar

[14] M. S. Grewal and A. P. Andrews, Kalman Filtering: Theory and Practice Using MATLAB, John Wiley and Sons, Inc, New York, 2001.10.1002/0471266388Search in Google Scholar

[15] T. Grindinger, Eye movement analysis & prediction with the Kalman filter, (MS thesis), Clemson University, 2006.Search in Google Scholar

[16] J. P. Dubois, J. S. Daba, M. Nader, and C. El Ferkh, GSM position tracking using a Kalman filter, Int. J. Electron. Commun. Eng. 6 (2012), no. 8, 867–876, https://publications.waset.org/vol/68.Search in Google Scholar

[17] N. Assimakis and M. Adam, Mobile position tracking in three dimensions using Kalman and Lainiotis filters, Open Math. J. 8 (2015), 1–6, http://doi.org/10.2174/1874117701508010001.10.2174/1874117701508010001Search in Google Scholar

[18] R. Zanetti and C. D’Souza, Recursive implementations of the Schmidt-Kalman ‘consider’ filter, J. of Astronaut. Sci. 60 (2013), no. 3–4, 672–685, http://doi.org/10.1007/s40295-015-0068-7.10.1007/s40295-015-0068-7Search in Google Scholar

[19] R. W. Farebrother, Linear Least Squares Computations, STATISTICS: Textbooks and Monographs, Marcel Dekker, New York, 1988.Search in Google Scholar

[20] N. Assimakis and M. Adam, Discrete time Kalman and Lainiotis filters comparison, Internat. J. Math. Anal. 1 (2007), no. 13, 635–659.Search in Google Scholar