# Iterative Algorithms for Symmetric Positive Semidefinite Solutions of the Lyapunov Matrix Equations

Authors: Min Sun; Jing Liu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

It is well-known that the stability of a first-order autonomous system can be determined by testing the symmetric positive definite solutions of associated Lyapunov matrix equations. However, the research on the constrained solutions of the Lyapunov matrix equations is quite few. In this paper, we present three iterative algorithms for symmetric positive semidefinite solutions of the Lyapunov matrix equations. The first and second iterative algorithms are based on the relaxed proximal point algorithm (RPPA) and the Peaceman–Rachford splitting method (PRSM), respectively, and their global convergence can be ensured by corresponding results in the literature. The third iterative algorithm is based on the famous alternating direction method of multipliers (ADMM), and its convergence is subsequently discussed in detail. Finally, numerical simulation results illustrate the effectiveness of the proposed iterative algorithms.


---

## Body

## 1. Introduction

Linear matrix equations are often encountered in various engineering disciplines, such as control theory [1, 2], system theory [3, 4] and stability analysis [5, 6]. For example, the stability of the first-order autonomous system x˙t=Axt is determined by whether the associated continuous-time Lyapunov equation XA+A⊤X=−M exists a positive definite solution or not, and the stability of the discrete-time system xk+1=Axk can be determined by solving the associated discrete-time Lyapunov matrix equation X−A⊤XA=M, where M is a given positive definite matrix with approximate size [7].In this paper, we consider the following constrained (continuous-time) Lyapunov matrix equations (termed as CLMEs): finding a matrixX∗∈S+m such that(1)AX∗+X∗A⊤=C,where S+m denotes the symmetric positive semidefinite (SPSD) matrix cone, A∈ℝm×m and C∈ℝm×m are two given constant matrices, and X∗∈ℝm×m is an unknown matrix to be determined. To solve problem (1), many iterative methods have been developed in the literature. For example, based on the hierarchical identification principle in parameter estimation, Ding et al. [8] designed a gradient-based iterative algorithm and a least-squares iterative algorithm for problem (1) without X∗∈S+m. Niu et al. [9] and Xie and Ma [10] extended the method in [8] by designing some accelerated gradient-based iterative algorithms for problem (1) without X∗∈S+m, which often have more efficient numerical performance by choosing suitable parameters. Based on the hierarchical identification principle and a new fixed point equation, Sun et al. [7] designed two least-squares iterative algorithms for problem (1) without X∗∈S+m. Conjugate gradient method is an influential numerical method for nonlinear optimization programming, which has also been applied to solve problem (1) without X∗∈S+m. For example, Huang and Ma [11] proposed a conjugate gradient method to solve problem (1) whether it is consistent or not. An appealing characteristic of the method proposed by Huang and Ma [11] is that: (i) it has finite termination property in the absence of round-off errors; (ii) it can obtain least Frobenius norm solution of problem (1) when it adopts some special kind of initial matrix. Furthermore, Ke and Ma [12] extended the alternating direction method of multipliers (ADMM), which is a famous numerical method in separable optimization programming, to solve problem (1) with the nonnegative constraint X∗≥0, which is obviously motivated by the design principle proposed by Xu et al. [13].Though the above iterative methods often have high efficient numerical performance, they only solve problem (1) without the constraint X∗∈S+m. However, we often need to solve problem (1) with the constraint X∗∈S+m in practice. In this paper, we are going to deal with this issue and propose three iterative algorithms, which are based on the famous relaxed proximal point algorithm (RPPA), the augmented Lagrangian method (ALM) and the alternating direction method of multipliers (ADMM), respectively. By transforming CLMEs into some optimization programming, we can determine the existence of solutions to CLMEs by judging that the optimal value of the transformed optimization programming equals to zero or not.Before ending this section, let us introduce the notations used in this paper. For any matrixN∈ℝm×n, we use N⊤ and tr(N) to denote its transpose and trace, respectively. The symbol Im stands for an identity matrix of order m. Let A,B=trB⊤A=∑ijaijbij, where the matrices A and B has the same size. The Frobenius norm N is defined as N=trN⊤N. The symbol A⊗B is defined as A⊗B=aijB, which represents the Kronecker product of matrices A and B. For a matrix A∈ℝm×n, the vectorization operator vec(A) is defined by vec(A) = a1⊤,a2⊤,…,an⊤⊤, where ak is the k–th column of the matrix A. According to the property of the Kronecker product, for any matrices M,N and X with approximate size, it holds that(2)vecMXN=N⊤⊗MvecX.The Kronecker product satisfies the following properties [14]:Property 1.
LetA∈ℝm×n and B∈ℝp×q. Then,(3)A⊗B=A⊗IpIn⊗B=Im⊗BA⊗Iq.Property 2.
LetA∈ℝm×n,C=∈ℝn×p,B∈ℝq×r, and D∈ℝr×s. Then,(4)A⊗BC⊗D=AC⊗BD.
LetPS+m⋅ denote the projection operator of a matrix onto the convex set S+m, i.e.,(5)PS+mA=arg minB∈S+mA−B2,which satisfies the following property [15]:(6)PS+mA−A,B−PS+mA≥0,∀B∈S+m.
Furthermore,PS+mA has a closed-form expression, i.e.,(7)PS+mA=UΛ+U⊤,U,Λ=eigA,Λ+=max0,Λ,and the computational load of the above procedure is 10m3 flops.
The remainder of this paper is organized as follows. In Sections2 and 3, we present the first and second iterative algorithm for problem (1) and analyze their global convergence. Section 4 discusses the third iterative algorithm for problem (1) and presents its global convergence. Section 5 presents some numerical results to illustrate the efficiency of the proposed algorithms. Finally, Section 6 concludes the entire paper and discusses future research opportunities.


## 2. The First Iterative Algorithm and Its Convergence

In this section, motivated by the relaxed proximal point algorithm (RPPA) in [16], we shall present an iterative algorithm for the constrained Lyapunov matrix equations and discuss its convergence.Firstly, let us transform the CLMEs into the following constrained matrix optimization (termed as CMO):(8)minX∈S+m12AX+XA⊤−C2.Remark 1.
Obviously, the objective function of CMO (8) has lower bound zero, therefore, CMO (8) must have optimal solution, and from its optimal value, we can determine the existence of solutions of CLMEs. In fact, if the optimal value of CMO (8) equals to zero, CLMEs has at least one solution; otherwise, CLMEs does not have solution.
Since CMO (8) can be viewed as the matrix form of the traditionary convex optimization with vector variables, we apply the relaxed proximal point algorithm (termed as RPPA) proposed by Gol’shtein and Tret’yakov in to solve CMO (8), and get the following RPPA.Algorithm 1.
The RPPA for problem (8)Step 1. Input a parameterγ∈0,2, a symmetric positive semidefinite matrix G∈ℝm2×m2, a tolerance ε>0. Initialize X0∈S+m. Set k≔0.Step 2. Generate a temporal iterateX˜k by(9)X˜k=arg minX∈S+m12Im⊗A+A⊗ImvecX−vecC2+12vecX−vecXkG2,in whichxG2 denotes x⊤Gx for some vector x∈ℝm2.Step 3. If(10)Xk−X˜k<ε,then, stop and return an approximate solutionXk+1 of (8).Step 4. Compute the new iterate by(11)Xk+1=Xk−γXk−X˜k.Setk≔k+1, and goto Step 2.In the following, we set the matrixG in RPPA as(12)G=τIm2−Δ⊤Δ,Δ=Im⊗A+A⊗Im,τ≥Δ⊤Δ.Now, let us elaborate on the optimal subproblem in (9).For the givenXk, by some simple manipulations, the X subproblem in (9) is amount to solving the following problem:(13)X˜k=argminX∈S+mτ2vecX−Im2−Δ⊤ΔτvecXk−Δ⊤vecCτ2=argminX∈S+mX−Xk−A⊤AXk+A⊤XkA⊤+AXkA+XkA⊤Aτ+AC+CA⊤τ2=PS+mXk−A⊤AXk+A⊤XkA⊤+AXkA+XkA⊤Aτ+AC+CA⊤τ,where the second equality follows from the two properties in Section1, i.e.,(14)Δ⊤ΔvecXk=Im⊗A⊤A+A⊗A⊤+A⊤⊗A+A⊤A⊗ImvecXk=A⊤AXk+A⊤XkA⊤+AXkA+XkA⊤A,Δ⊤vecC=Im⊗A+A⊗ImvecC=A⊤C+CA.
Then,X˜k can be computed by the formula (7).
From [16], we have the following theorem, which means that the iterative sequence Xk generated by RPPA is Fejér monotone with respective to the solution set. The proof of its global convergence follows immediately from this property and the results in [17].Theorem 1.
Let the sequencesXk be generated by the RPPA. Then, for any k≥0, one has(15)Xk+1−X∗≤Xk−X∗−γ2−γXk−X˜k,where X∗ is a solution of CMO.


## 3. The Second Iterative Algorithm and Its Convergence

In this section, based on the Peaceman–Rachford splitting methods in [18, 19], we shall present another iterative algorithm for the constrained Lyapunov matrix equations and discuss its convergence.By introducing an auxiliary matrix variabley∈ℝm, CMO (8) can be further written as the following separable convex optimization (termed as SCO):(16)min12AX+XA⊤−C2s.t.X−Y=0,X∈ℝm×m,Y∈S+m.Since SCO (16) can be viewed as the matrix form of the traditionary convex optimization with vector variables, we apply the Peaceman–Rachford splitting method in [18, 19] to solve SCO (15), and get the following Peaceman–Rachford splitting method (termed as PRSM).Algorithm 2.
The PRSM for problem (16)Step 1. Input two parametersα∈0,1,β>0, a tolerance ε>0. Initialize X0∈ℝm×m, Y0∈S+m, Λ0∈ℝm×m. Set k≔0.Step 2. Generate the new iterateWk+1=Xk+1,Yk+1,Λk+1 by(17)Xk+1=argminX∈ℝm×m12AX+XA⊤−C2−Λk⊤X+β2X−Yk2,Λk+1/2=Λk−αβXk+1−Yk,Yk+1=argminY∈S+mΛk+1/2⊤Y+β2Xk+1−Y2,Λk+1=Λk+1/2−αβXk+1−Yk+1.Step 3. If(18)maxXk−Xk+1,Yk−Yk+1,Λk−Λk+1<ε,then stop and return an approximate solutionXk+1,Yk+1 of (16); else set k≔k+1, and goto Step 2.Now, let us elaborate on the two optimization subproblems in (17).For the givenYk, Λk, the X subproblem in (17) is amount to solving the following normal equation:(19)A⊤R+RA+βX−Yk−Λkβ=0,where(20)R=AX+XA⊤−C.Thus, we have(21)vecXk+1=Im⊗A⊤A+A⊗A⊤+A⊤⊗A+AA⊤⊗Im+βIm2−1vecA⊤C+CA+βYk+Λk.For the givenXk+1, Λk+1/2, the Y subproblem in (17) is amount to solving the following optimization problem:(22)Yk+1=argminY∈S+m12Y−Xk+1−Λk+1/2/β+Xk+1−Λk+1/2/β⊤22=PS+mXk+1−Λk+1/2/β+Xk+1−Λk+1/2/β⊤2,which can be computed by the formula given in (7).Remark 2.
Though the computation ofXk+1 involves an inverse of a matrix. However, since this term is invariant in each iteration, we only need to compute it once before all iterations.
From [18, 19], we have the following theorem, which means that the iterative sequence Wk generated by PRSM is Fejér monotone with respective to the solution set. Obviously, SCO is a convex optimization problem. Thus, its global convergence follows immediately from this property and the results in [17].Theorem 2.
LetWk be the sequence generated by the PRSM. Then, for any k≥0, one has(23)Wk+1−W∗H2≤Wk−W∗H2−Wk−W˜kN2,where H and N are two positive definite matrices, W∗=X∗,Y∗,Λ∗ with X∗,Y∗ being a solution of SCO, Λ∗ being the corresponding Lagrangian multiplier, and W˜k is an auxiliary sequence defined as(24)W˜k=X˜kY˜kΛ˜k=Xk+1Yk+1Λk−βXk+1−Yk.Remark 3.
Following the above procedure, some famous iterative algorithms for convex programming, such as the proximal point algorithm (PPA), the relaxed PPA [16], can also be used to solve SCO.


## 4. The Third Iterative Algorithm and Its Convergence

In this section, motivated by the alternating direction methods of multiplier (ADMM) in [12, 13], we shall present the third iterative algorithm for the constrained Lyapunov matrix equations and discuss its convergence.Firstly, CMO (8) can be written as the following equivalent form:(25)min12AX+YA⊤−C2s.t.X−Y=0,Y−Z=0,X∈ℝm×m,Y∈ℝm×m,Z∈S+m.The augmented Lagrangian function of (25) is(26)ℒX,Y,Z,Λ,Π=12AX+YA⊤−C2−Λ,X−Y−Π,Y−Z+α2X−Y2+β2Y−Z2,where the matrices Λ,Π∈ℝm×m are Lagrangian multipliers and the constants α,β>0 are penalty parameters. Then, we can derive the second iterative algorithm (termed as ADMM) for CMO by alternating minimizing one variable in ℒX,Y,Z,Λ,Π when other variables is fixed.Algorithm 3.
The ADMM for problem (25)Step 1. Input three parametersα, β>0, γ∈0,1+5/2, a tolerance ε>0. Initialize X0, Y0, Λ0, Π0∈ℝm×m, and Z0∈S+m. Set k≔0.Step 2. Compute the new iterateWk+1=Xk+1,Yk+1,Zk+1,Λk+1,Πk+1 by(27)Xk+1=argminX∈ℝm×mℒX,Yk,Zk,Λk,Πk,Yk+1=argminY∈ℝm×mℒXk+1,Y,Zk,Λk,Πk,Zk+1=argminZ∈S+mℒXk+1,Yk+1,Z,Λk,Πk,Λk+1=Λk−γαXk+1−Yk+1,Πk+1=Πk−γβYk+1−Zk+1.Step 3. If(28)maxXk−Xk+1,Yk−Yk+1,Zk−Zk+1,Λk−Λk+1,Πk−Πk+1<ε,then stop and return an approximate solutionXk+1,Yk+1,Zk+1 of (25); else, set k≔k+1, and goto Step 2.Now, let us elaborate on the three optimization subproblems in (27).For the givenYk, Zk, Λk, and Πk, the X subproblem in (27) is amount to solving the following normal equation:(29)A⊤AX+YkA⊤−C−Λk+αX−Yk=0,from which, we have(30)Xk+1=αIm+A⊤A−1A⊤C−A⊤YkA⊤+Λk+αYk.For the givenXk+1, Zk, Λk, and Πk, the Y subproblem in (27) is amount to solving the following normal equation:(31)AXk+1+YA⊤−CA+Λk−Πk+αY−Xk+1+βY−Zk=0,from which, it holds that(32)Yk+1=CA−AXk+1A−Λk+Πk+αXk+1+βZkα+βIm+A⊤A−1.For the givenXk+1, Yk+1, Λk, and Πk, the Z subproblem in (27) is amount to solving the following optimization problem:(33)Zk+1=arg minZ∈S+m12Z−Yk+1−Πk/β+Yk+1−Πk/β⊤22=PS+mYk+1−Πk/β+Yk+1−Πk/β⊤2,
which can be computed by the formula given in (7).Remark 4.
It is obvious that the two matricesαIm+A⊤A and α+βIm+A⊤A are symmetric positive definite for any real constants α, β>0, which are invariant during the iteration, and we only need to compute their inverses once before all iterations.
The KKT conditions for problem (25) are:(34)A⊤AX+YA⊤−C−Λ=0,AX+YA⊤−CA+Λ−Π=0,Π−Γ=0,X−Y=0,Y−Z=0,Γ,Z=0,Z∈S+m,Γ∈S+m,where Λ, Π, and Γ are three Lagrangian multipliers corresponding to the constraints X−Y=0, Y−Z=0 and Z∈S+m in (25), respectively. Eliminating the variable Γ, the above KKT system can be written as:(35)A⊤AX+YA⊤−C−Λ=0,(36)AX+YA⊤−CA+Λ−Π=0,(37)X−Y=0,(38)Y−Z=0,(39)Π,Z=0,Z∈S+m,Π∈S+m.
If the pointX,Y,Z,Λ,Π satisfies the above conditions, then X,Y,Z is called KKT point, and Λ,Π is the corresponding multiplier.
To simplify notation, we setW=X,Y,Z,Λ,Π. Now, let us analyze the convergence of ADMM. Due to the objective function of (25) is not separable, we only get the following weak convergence result.Theorem 3.
LetWk be the sequence generated by ADMM which satisfies the condition(40)limk⟶∞Wk+1−Wk=0.Then, any accumulation point ofWk satisfies the KKT conditions for problem (25).Proof.
From (11)–(13), it holds that(41)αIm+A⊤AXk+1−Xk=A⊤C−AXk−YkA⊤+Λk+αYk−Xk,Yk+1−Ykα+βIm+A⊤A=C−AXk+1−YkA⊤A−Λk+Πk+αXk+1−Yk+βZk−Yk,Zk+1−Zk=PS+mYk+1−Πk/β+Yk+1−Πk/β⊤2−Zk,Λk+1−Λk=−γαXk+1−Yk+1,Πk+1−Πk=−γβYk+1−Zk+1.
Fromlimk⟶∞Wk+1−Wk=0, we can find that both the sides in the above five equations tend entirely to zero, i.e.,(42)A⊤C−AXk−YkA⊤+Λk⟶0,(43)C−AXk+1−YkA⊤A−Λk+Πk⟶0,(44)PS+mYk+1−Πk/β+Yk+1−Πk/β⊤/2−Zk⟶0,(45)Xk+1−Yk+1⟶0,(46)Yk+1−Zk+1⟶0.
For any limitW˜=X˜,Y˜,Z˜,Λ˜,Π˜ of the sequence Wk, from (42), (43), (45) and (46), it is obvious that the first four equations (14)–(17) in the KKT conditions are satisfied.
In the following, we only need to prove that (39) also hold for W˜. From (35) and (36), it holds that(47)A⊤AX˜+Y˜A⊤−C−Λ˜=0,AX˜+Y˜A⊤−CA+Λ˜−Π˜=0.
This together withX˜=Y˜ implies that(48)Π˜=AX˜A+X˜A⊤A−CA+A⊤AX˜+A⊤X˜A⊤−A⊤C.
This and the symmetry ofX˜ and C imply that the matrix Π˜ is also symmetric. From (44) and (46), it holds that(49)Z˜=PS+mZ˜−Π˜β.
From this and (6), we have(50)Π˜,B−Z˜≥0,∀B∈S+m.
SettingB=0 and B=2Z˜ in the above inequality, we get(51)Π˜,Z˜≤0,Π˜,Z˜≥0.
Thus,Π˜,Z˜=0. Let λi i=1,2,…,m with λ1≥λ2≥⋯≥λm be the m eigenvalues of the matrix Π˜. Suppose the matrix Π˜ is not positive semidefinite, then λm<0. Let P∈ℝm×m be an orthogonal matrix, such that Π˜=Pdiagλ1,λ2,…,λmP⊤. Setting B=Z˜+Pdiag0,0,…,tP⊤∈S+m with t>0 in (50), we have(52)Π˜,Pdiag0,0,…,tP⊤≥0.
Thus,(53)P⊤Π˜P,diag0,0,…,t≥0,
That is,(54)diagλ1,λ2,…,λm,diag0,0,…,t≥0.
Then,tλm≥0, which contradicts with t>0 and λm<0. Thus, the matrix Π˜∈S+m. This completes the proof.Remark 5.
Based on the proximal ADMM [20], we can design an inverse-free ADMM for the studied problem by attaching some proximal terms to the subproblems with respective to X and Y in ADMM. Furthermore, based on the parallel or partially parallel splitting method or their proximal versions [21, 22], we can design some inverse-free (partially) parallel splitting methods for the studied problem [23, 24].


## 5. Numerical Results

In this section, we evaluate the performance of the three iterative algorithms designed in this paper for solving the Lyapunov matrix equations, and report some preliminary numerical results. The codes of our algorithms were written entirely in Matlab R2014a, and were executed on a ThinkPad notebook with Pentium(R) Dual-Core CPU T4400@2.2 GHz, 2 GB of memory.The parameters in the three tested algorithms are set as follows:For RPPA:τ=Δ⊤Δ, γ=1.8For PRSM:α=0.9, β=2For ADMM:α=2, β=2, γ=1.618Problem 1.
Consider the following constrained Lyapunov matrix equations:(55)AX∗+X∗A⊤=−C,X∗∈S+m,where A and C are set as the following three cases:(56)A=1221,X∗=2111,C=−AX∗−X∗A⊤=−8−8−8−6,A=1221,X∗=2011,C=−AX∗−X∗A⊤=−6−6−8−4,A=1221,X∗=−1111,C=−AX∗−X∗A⊤=−2−2−2−6.
All initial matrices are set zero matrix. Due to the iterateYk generated by PRSM satisfies the constraint Y∈S+m and(57)limk⟶∞Xk−Yk=0,
and RPPA, ADMM also have similar property, then we adopt the following stopping criterion:(58)Errork=AXk+XkA⊤+C<10−6,
or the number of iterationsk>1000.
The numerical results of the three tested algorithms for the first case of Problem 1 are plotted in Figure1, from which we find that they successfully solve this problem within 66, 35, and 37 iterations, respectively, and the final iterates of the three tested algorithms are(59)XRPPA66=1.9999995881055641.0000001285789551.0000001285789551.000000128578955,XPRSM35=1.9999995614815410.9999995614815410.9999995614815400.999999561481540,XADMM37=2.0000460663423171.0000460696498651.0000460627317251.000046066039320.
The corresponding absolute errors are(60)XRPPA66−X∗≅9.3941×10−7,XPRSM35−X∗≅8.7704×10−7,XADMM37−X∗≅6.0869×10−7.
Though the matrix(61)X∗=2011.
Satisfies the equationsAX+XA⊤=−C, it is not symmetric. Then, it is not a solution of CMO. The numerical results of the three tested algorithms for the second case of Problem 1 are plotted in Figure 2, from which we find that the Error(k)s generated by the three algorithms do not converge to zero, and the final iterate is about(62)X1000=20.50.51.
Similarly, the matrix(63)X∗=−1111,satisfies the equations AX+XA⊤=−C, it is not positive definite. Then, it is not a solution of CMO. The numerical results of the three tested algorithms for the third case of Problem 1 are plotted in Figure 3, from which we find that the Error(k)s generated by the three algorithms do not converge to zero, and the final iterate is about(64)X1000=0.156530.410450.410451.0763.Figure 1
Numerical results of RPPA, PRSM, and ADMM for the first case of Problem 1.Figure 2
Numerical results of RPPA, PRSM and ADMM for the second case of Problem 1.Figure 3
Numerical results of RPPA, PRSM, and ADMM for the third case of Problem 1.Problem 2.
Consider the following constrained Lyapunov matrix equations:(65)AX∗+X∗A⊤=−C,X∗∈S+m,where A and C are set as(66)A=1321,X∗=2111,C=−AX∗−X∗A⊤=−10−9−10−6.
Three curves in Figure4 indicates that for CLMEs with nonsymmetric matrix A, RPPA and PRSM fall into the premature convergence, while ADMM successfully solve this problem.Figure 4
Numerical results of RPPA, PRSM, and ADMM for Problem 2.Problem 3.
Consider the following constrained Lyapunov matrix equations:(67)AX∗+X∗A⊤=−C,X∗∈S+m,where A and X∗ are set by(68)A=randn,A=A+A′,D=diagrandn,1,U=orthrandn,n,X∗=U′∗D∗U.
Then,C is set as C=−AX∗+X∗A⊤. Furthermore, we set n=20.
The numerical results of the three tested algorithms are plotted in Figure5, from which we find that the performance of ADMM is a litter better than that of PRSM, which is obviously better than that of RPPA.Figure 5
Numerical results of RPPA, PRSM, and ADMM for Problem 3.


## 6. Conclusions

In this paper, three iterative algorithms have been proposed to solve constrained Lyapunov matrix equations, and their convergence is discussed in detail. Furthermore, preliminary numerical results have been included, which verify the efficiency of the proposed algorithms.Among various challenges for large-scale constrained Lyapunov matrix equations, high performance is a must, not an option. Therefore, an object of our future investigations is to design more efficient iterative algorithms and propose some valid criteria, from which we can directly determine whether the Lyapunov matrix equations has a symmetric positive semidefinite solution or not.



---

*Source: 6968402-2020-02-10.xml*