# Refinement Methods for State Estimation via Sylvester-Observer Equation

Authors: H. Saberi Najafi; A. H. Refahi Sheikhani
Journal: Advances in Numerical Analysis
 (2011)

---

## Abstract

We present new iterative methods based on refinement process for solving large sparse Sylvester-observer
equations applied in state estimation of a continuous-time system. These methods use projection methods to produce low-dimensional Sylvester-observer matrix equations that are solved by the direct methods. Moreover, the refinement process
described in this paper has the capability of improving the results obtained by any other methods. Some numerical results
will be reported to illustrate the efficiency of the proposed methods.


---

## Body

## 1. Introduction

Consider the continuous-time linear system(1.1)ẋ(t)=Ax(t)+Bu(t),y(t)=Cx(t),
where A∈Rn×n, B∈Rn×m, and C∈Rr×n.We well know that all the state-feedback problems,  such as the feedback stabilization, the LQR, and the state-feedbackH∞-control problems, require the state vector x(t) explicitly [1]. In most practical situations, the states are not fully accessible and the designer only knows the output and the input vectors. The unavailable states somehow need to be estimated accurately from the knowledge of the matrices A, B, and C, the output vector y(t), and the input vector u(t). There are two common procedures for state estimation: one via eigenvalue assignment (EVA) and the other via solution of the Sylvester-observer equation. The main step in state estimation via solution of the Sylvester-observer equation is solving the Sylvester-observer equation of the form (1.2)XA-FX=GC,
where X, A, F∈Rn×n, G∈Rn×r, and C∈Rr×n. Sylvester-observer equations (1.2) play vital roles in a number of applications such as control and communications theory [1], model reduction [2–4], numerical solution of matrix differential Riccati equations [5], and image processing [6].The analytical solution of the matrix equation (1.2) has been considered by many authors; see [7, 8]. Direct methods for solving the matrix equation (1.2) are attractive if the matrices are of small size. These methods are based on the Schur decomposition, by which the original equation is transformed into a form that is easy to be solved by a forward substitution. Moreover, the matrices of the large practical problem are very sparse. Since, the standard methods for solving the Sylvester equations destroy the sparsity of the problems, they are only applicable for the matrices of small size; see [9–11]. Iterative projection methods for solving large Sylvester-observer matrix equations have been developed during the past years; see [12–16]. These methods use the Arnoldi process to compute an orthonormal basis of certain Krylov subspace. In this paper, we extend the idea to propose a new projection method for solving (1.1) based on Weighted block Krylov subspace process. These methods are based on the reduction of the large sparse Sylvester-observer equation to the low-dimensional problem by orthogonal similarity that is solved by the direct methods. Moreover, refinement process presented in Sections 4 and 5 has the capability of improving the results obtained by any other methods.The remainder of the paper is organized as follows. In Section2, we describe some fundamental results about control theory. Then, we discuss how the states of a continuous-time system can be estimated in Section 3. In Sections 4 and 5, we introduce two new iterative methods for solving large sparse Sylvester-observer equation. These methods are based on the reduction of the large sparse Sylvester-observer equation to the low-dimension problem. Section 6 is devoted to some numerical tests. Some concluding remarks are given in Section 7.


## 2. Some Fundamental Results

The two basic concepts in control theory are controllability and observability of a control system. In this section, we will state some well-known facts about controllability and observability for convenient use later in the paper. For an excellent account of controllability and observability and related results, see [1].Definition 2.1.
The system (1.1) is said to be controllable if, starting from any initial state x(0), the system can be driven to any final state xf in some finite time tf, choosing the input vector u(t), 0≤t≤tf, appropriately.Observability is the dual concept of controllability.Definition 2.2.
The system (1.1) is said to be observable if there exists t1>0 such that the initial state x(0) can be uniquely determined from the knowledge of u(t) and y(t) for all0≤t≤t1.Remark 2.3.
Since matrixC does not have any role in the definition of controllability, the controllability of the system (1.1) (see [1]), is often referred to as the controllability of the pair (A,B). Similarly, since matrix B does not have any role in the definition of observability, the observability of the system (1.1) is often referred to as the observability of the pair (A,C).Some well-known criteria of controllability and observability are now stated in the next two theorems. The proofs of Theorems2.4 and 2.5 can be found in [1]. In the following, A is n×n, B is n×m(m≤n), and C is r×n(r≤n).Theorem 2.4.
The pair(A,B) is controllable if and only if any one of the following equivalent conditions holds.
(1)
The controllability matrix(2.1)CM=(B,AB,A2B,…,An-1B),
has rank n.(2)
Rank(A-λI,B)=n for every eigenvalue λ of A.(3)
Let(λ,x) be an eigenpair of A, that is, Ay=λy; then Cy≠0.Proof.
See [1].Theorem 2.5.
The pair(A,C) is observable if and only if any one of the following equivalent conditions holds
(1)
The observability matrix(2.2)OM=(CCACA2⋮CAn-1)
has rank n.(2)
The matrix(λI-AC) has rank n for every eigenvalue λ of A.(3)
Let(λ,y) be an eigenpair of A, that is, Ay=λy; then Cy≠0.Proof.
See [1].Definition 2.6.
A matrixA is called a stable matrix if all of the eigenvalues of A have negative real parts.


## 3. State Estimation

In this section, we discuss how the states of a continuous-time linear system (1.1) can be estimated. The discussions here apply equally to the discrete-time systems possibly with some minor changes. So we concentrate on the continuous-time case only. We describe two common procedures for state estimation: one via eigenvalue assignment (EVA) and the other,via solution of the Sylvester-observer equation.


### 3.1. State Estimation via Eigenvalue Assignment

Now consider the linear time-invariant system (1.1). Let x̂(t) be an estimate of the state vector x(t). Obviously, we would like to construct the vector x̂(t) in such a way that the error e(t)=x(t)-x̂(t) approaches zero as fast as possible, for all initial states x(0) and for every input u(t). The following theorem shows that the problem of state estimation can be solved by finding a matrix K such that the matrix A-KC has a suitable desired spectrum. See [1].Theorem 3.1.
If(A,C) is observable, then the states x(t) of the system (1.1) can be estimated by
(3.1)x̂̇(t)=(A-KC)x̂(t)+Ky(t)+Bu(t),
where K is constructed such that A-KC is a stable matrix. The error e(t)=x(t)-x̂(t) is governed by
(3.2)ė(t)=(A-KC)e(t),
and e(t)→0 as t→∞.Proof.
See [1].


### 3.2. State Estimation via Sylvester Equation

There is another approach for state estimation. KnowingA, B, C, u(t), and y(t), let us construct the system (3.3)ż(t)=Fz(t)+Gy(t)+Pu(t),
where F is n×n, G is n×r, and P is n×m, in such a way that for some constant n×n nonsingular matrix X the error vector e(t)=z(t)-Xx(t)→0 for all x(0), z(0), and for every input u(t). The vector z(t) will then be an estimate of Xx(t). The system (3.3) is then said to be the state observer for the system (1.1). D. Luenberger originated the idea and is hence referred to in control theory as the Luenberger observer; see [1].We now show that the system (3.3) will be a state observer if the matrices X, F, G, and P satisfy certain requirements.Theorem 3.2.
The system (3.3) is a state observer of the system (1.1), that is, z(t) is an estimate of Xx(t) in the sense that the error e(t)=z(t)-Xx(t)→0 as t→∞ for any initial conditions x(0), z(0), and u(t) if
(1)XA-FX=GC,
(2)P=XB,
(3)F is stable.Proof.
See [1].Definition 3.3.
The matrix equation(3.4)XA-FX=GC,
where X, A, F∈Rn×n, G∈Rn×r, and C∈Rr×n, is called the Sylvester-observer equation.Theorem3.2 suggests the following method for the observer design; see [1].Algorithm 3.4 (full-order observer design via Sylvester-observer equation).
One has the following.
Inputs
The system matricesA, B, and C of order n×n, n×m, and r×n, respectively.Output
An estimatex̂(t) of the state vector x(t).Assumptions.
(A,C) is observable.Step 1.
Find a nonsingular solution of the Sylvester-observer equation (1.2) by choosing F as a stable matrix and choosing G in such a way that the resulting solution X is nonsingular.Step 2.
ComputeP=XB.Step 3.
Construct the observerz(t) by solving the system of differential equations
(3.5)ẑ(t)=Fz(t)+Gy(t)+Pu(t),z(0)=z0.Step 4.
Find an estimatex̂(t) of x(t): x̂(t)=X-1z(t).


### 3.3. Characterization of Nonsingular Solutions of the Sylvester Equation

In this section, we describe some necessary conditions for a unique solution of the Sylvester equation to have such properties. The following theorem was proved by Bhattacharyya and de Souza. The proof here has been taken from [1].Theorem 3.5.
LetA, F, G, and C, respectively, be of order n×n, n×n, n×r, and r×n. Let X be a unique solution of the Sylvester-observer equation (1.2). Then, necessary conditions for X to be nonsingular are that (A,C) is observable and (F,G) is controllable.Proof.
See [1].Corollary 3.6.
IfG is n×1 and C is 1×n, then necessary and sufficient conditions for the unique solution X of (1.2) to be nonsingular are that (F,G) is controllable and (A,C) is observable.Proof.
See [1].Remark 3.7.
According to Theorem3.5 and Corollary 3.6, the controllability of (F,G) and observability of (A,C) guarantee the existence of a nonsingular solution of the Sylvester-observer equation (1.2) in Step 1 of Algorithm 3.4. Moreover, there are other choices for F and G in Step 1 of Algorithm 3.4 provided (F,G) is controllable. Also, we can use Theorems 2.4 and 2.5 for analyzing the controllability of (F,G) and the observability of (A,C).Example 3.8.
In this example we show how the Sylvester-observer equation (1.2) can be applied in state estimation of a continuous time system (1.1). In this sense, at first we use the MATLAB function ode23 for directly solving (1.1) with u(t) as the unit step function and x(0)=(6,0)T. Then, we apply Algorithm 3.4 for computing the estimate x̂(t). Also, the differential equation (3.12) was solved with z(0)=0 and MATLAB function ode23. The comparison of the state x(t) and estimate x̂(t) for this example is shown in Figure 1. The solid line corresponds to the exact state and the dotted line corresponds to the estimate state. Let
(3.6)A=(1111),B=(10),C=(10).
According to criteria 1 of Theorem 2.5, (A,C) is observable. Thus, we can use Algorithm 3.4 for state estimation of (1.1).The (a) first and (b) second variables of the statex(t) and estimate x̂(t) for Example 3.8.
(a)(b)Step 5.
Choose(3.7)G=(13),F=diag⁡(-1,-3).
According to criteria 1 of Theorem  2.4, (F,G) is controllable. Thus, by Corollary 3.6, the nonsingular  solution X of XA-FX=GC is
(3.8)X=(0.6667-0.33330.8000-0.2000),
(computed by MATLAB function lyap).Step 6.
One has(3.9)P=XB=(0.66670.8000).Step 7.
An estimatex̂(t) of x(t) is
(3.10)x̂(t)=X-1z(t)=(-1.52.5-65)(z1(t)z2(t))=(-1.5z1+2.5z2-6z1+5z2),
where
(3.11)z(t)=(z1(t)z2(t)),
is given by
(3.12)ẑ(t)=(-100-3)z(t)+(13)y(t)+(0.66670.8000)u(t),z(0)=z0.Remark 3.9.
According to Algorithm3.4 the most important step is Step 1.  In the case that n is small, there are many reliable algorithm s for solving the Sylvester-observer equation and the states of a continuous-time system can be estimated. However, for large and sparse systems solving the Sylvester-observer equation by the available methods can be costly. In Sections 4 and 5, we introduce two iterative refinement methods for solving large sparse Sylvester-observer equations.


## 3.1. State Estimation via Eigenvalue Assignment

Now consider the linear time-invariant system (1.1). Let x̂(t) be an estimate of the state vector x(t). Obviously, we would like to construct the vector x̂(t) in such a way that the error e(t)=x(t)-x̂(t) approaches zero as fast as possible, for all initial states x(0) and for every input u(t). The following theorem shows that the problem of state estimation can be solved by finding a matrix K such that the matrix A-KC has a suitable desired spectrum. See [1].Theorem 3.1.
If(A,C) is observable, then the states x(t) of the system (1.1) can be estimated by
(3.1)x̂̇(t)=(A-KC)x̂(t)+Ky(t)+Bu(t),
where K is constructed such that A-KC is a stable matrix. The error e(t)=x(t)-x̂(t) is governed by
(3.2)ė(t)=(A-KC)e(t),
and e(t)→0 as t→∞.Proof.
See [1].


## 3.2. State Estimation via Sylvester Equation

There is another approach for state estimation. KnowingA, B, C, u(t), and y(t), let us construct the system (3.3)ż(t)=Fz(t)+Gy(t)+Pu(t),
where F is n×n, G is n×r, and P is n×m, in such a way that for some constant n×n nonsingular matrix X the error vector e(t)=z(t)-Xx(t)→0 for all x(0), z(0), and for every input u(t). The vector z(t) will then be an estimate of Xx(t). The system (3.3) is then said to be the state observer for the system (1.1). D. Luenberger originated the idea and is hence referred to in control theory as the Luenberger observer; see [1].We now show that the system (3.3) will be a state observer if the matrices X, F, G, and P satisfy certain requirements.Theorem 3.2.
The system (3.3) is a state observer of the system (1.1), that is, z(t) is an estimate of Xx(t) in the sense that the error e(t)=z(t)-Xx(t)→0 as t→∞ for any initial conditions x(0), z(0), and u(t) if
(1)XA-FX=GC,
(2)P=XB,
(3)F is stable.Proof.
See [1].Definition 3.3.
The matrix equation(3.4)XA-FX=GC,
where X, A, F∈Rn×n, G∈Rn×r, and C∈Rr×n, is called the Sylvester-observer equation.Theorem3.2 suggests the following method for the observer design; see [1].Algorithm 3.4 (full-order observer design via Sylvester-observer equation).
One has the following.
Inputs
The system matricesA, B, and C of order n×n, n×m, and r×n, respectively.Output
An estimatex̂(t) of the state vector x(t).Assumptions.
(A,C) is observable.Step 1.
Find a nonsingular solution of the Sylvester-observer equation (1.2) by choosing F as a stable matrix and choosing G in such a way that the resulting solution X is nonsingular.Step 2.
ComputeP=XB.Step 3.
Construct the observerz(t) by solving the system of differential equations
(3.5)ẑ(t)=Fz(t)+Gy(t)+Pu(t),z(0)=z0.Step 4.
Find an estimatex̂(t) of x(t): x̂(t)=X-1z(t).


## 3.3. Characterization of Nonsingular Solutions of the Sylvester Equation

In this section, we describe some necessary conditions for a unique solution of the Sylvester equation to have such properties. The following theorem was proved by Bhattacharyya and de Souza. The proof here has been taken from [1].Theorem 3.5.
LetA, F, G, and C, respectively, be of order n×n, n×n, n×r, and r×n. Let X be a unique solution of the Sylvester-observer equation (1.2). Then, necessary conditions for X to be nonsingular are that (A,C) is observable and (F,G) is controllable.Proof.
See [1].Corollary 3.6.
IfG is n×1 and C is 1×n, then necessary and sufficient conditions for the unique solution X of (1.2) to be nonsingular are that (F,G) is controllable and (A,C) is observable.Proof.
See [1].Remark 3.7.
According to Theorem3.5 and Corollary 3.6, the controllability of (F,G) and observability of (A,C) guarantee the existence of a nonsingular solution of the Sylvester-observer equation (1.2) in Step 1 of Algorithm 3.4. Moreover, there are other choices for F and G in Step 1 of Algorithm 3.4 provided (F,G) is controllable. Also, we can use Theorems 2.4 and 2.5 for analyzing the controllability of (F,G) and the observability of (A,C).Example 3.8.
In this example we show how the Sylvester-observer equation (1.2) can be applied in state estimation of a continuous time system (1.1). In this sense, at first we use the MATLAB function ode23 for directly solving (1.1) with u(t) as the unit step function and x(0)=(6,0)T. Then, we apply Algorithm 3.4 for computing the estimate x̂(t). Also, the differential equation (3.12) was solved with z(0)=0 and MATLAB function ode23. The comparison of the state x(t) and estimate x̂(t) for this example is shown in Figure 1. The solid line corresponds to the exact state and the dotted line corresponds to the estimate state. Let
(3.6)A=(1111),B=(10),C=(10).
According to criteria 1 of Theorem 2.5, (A,C) is observable. Thus, we can use Algorithm 3.4 for state estimation of (1.1).The (a) first and (b) second variables of the statex(t) and estimate x̂(t) for Example 3.8.
(a)(b)Step 5.
Choose(3.7)G=(13),F=diag⁡(-1,-3).
According to criteria 1 of Theorem  2.4, (F,G) is controllable. Thus, by Corollary 3.6, the nonsingular  solution X of XA-FX=GC is
(3.8)X=(0.6667-0.33330.8000-0.2000),
(computed by MATLAB function lyap).Step 6.
One has(3.9)P=XB=(0.66670.8000).Step 7.
An estimatex̂(t) of x(t) is
(3.10)x̂(t)=X-1z(t)=(-1.52.5-65)(z1(t)z2(t))=(-1.5z1+2.5z2-6z1+5z2),
where
(3.11)z(t)=(z1(t)z2(t)),
is given by
(3.12)ẑ(t)=(-100-3)z(t)+(13)y(t)+(0.66670.8000)u(t),z(0)=z0.Remark 3.9.
According to Algorithm3.4 the most important step is Step 1.  In the case that n is small, there are many reliable algorithm s for solving the Sylvester-observer equation and the states of a continuous-time system can be estimated. However, for large and sparse systems solving the Sylvester-observer equation by the available methods can be costly. In Sections 4 and 5, we introduce two iterative refinement methods for solving large sparse Sylvester-observer equations.


## 4. Block Refinement Method

As we already mentioned, so far many numerical methods have been developed by different authors; see [1, 13, 17].  For example, the Hessenberg-Schur method is now widely used as an effective computational method for the Sylvester-observer equation. But numerical stability of this method has not been investigated. As the iterative methods are very efficient for the solution of computational problems, we thought it will be good idea to create an iterative method for solving the Sylvester-observer equation XA-FX=GC where A, F, X∈Rn×n, G∈Rn×r, and C∈Rr×n. In this section we propose to show that the obtained approximate solution of the Sylvester-observer equation by any method can be improved, in other words the accuracy can be increased. If this idea is applicable, then we have found an iterative method for solving of the Sylvester-observer equation.Theorem 4.1.
LetX0 be the approximate solution obtained by an arbitrary method for the matrix equation (1.2), and let R(X0)=GC-(X0A-FX0) be the corresponding residual. If m steps of the block Arnoldi process for matrices AT and F have been run (m≪n) and Ym∈Rm×m is the solution of the low-dimensional Sylvester equation
(4.1)YmÂmT-F̂mYm=WmTR(X0)Vm,
then
(4.2)R(X1)=R(X0+WmYmVmT)=0.Proof.
LetVm=[v1,…,vm] and Wm=[w1,…,wm] be two orthogonal bases constructed by the block Arnoldi process for the matrices AT and F, respectively. Thus, we have
(4.3)VmTVm=Im,WmTWm=Im.
Also, the square block Hessenberg matrices Âm and F̂m (m=rℓ, where r and ℓ are the dimensions of blocks) whose nonzero entries are the scalars âij and f̂ij, constructed by the Block Arnoldi process, can be expressed as
(4.4)Âm=VmTATVm,F̂m=WmTFWm.
If we set(4.5)X1=X0+WmYmVmT,
then the corresponding residual R(X1)=GC-(X1A-FX1) satisfies
(4.6)R(X1)=GC-((X0+WmYmVmT)A-F(X0+WmYmVmT)),=R(X0)-WmYmVmTA+FWmYmVmT,=R(X0)-WmYmÂmTVmT+WmF̂mYmVmT,=R(X0)-Wm(YmÂmT-F̂mYm)VmT.
Since Ym is the solution of (4.1) and by using (4.3), we have
(4.7)R(X1)=R(X0)-WmWmTR0VmVmT=0.According to Theorem4.1, we can develop an iterative method for solving the Sylvester-observer equation when the matrices A, F, G, and C are large and sparse. For achieving this idea, if we choose m<n, then instead of solving XA-FX=GC we can solve (4.1). In other words, in this method, first we transform the initial Sylvester-obsever equation to another Sylvester equation with less dimensions and then in each iteration step solve this matrix equation and extend the obtained solution to the solution of the initial equation by (4.5). The algorithm  is as follows.Algorithm 4.2 (block refinement method).
(1) Start: choose an initial approximate solutionX0, and a tolerance ε.
(2) Select two numbersr and ℓ for dimensions of block and set m=rℓ(m<n).
(3) ComputeR(X0)=GC-(X0A-FX0).
(4) Construct the orthonormal basesVm and Wm∈Rn×m by the block Arnoldi process, such that
(4.8)Âm=VmTATVm,F̂m=WmTFWm.
(5) Solve the low-dimensional Sylvester-observer equationYmÂmT-F̂mYm=WmTR(X0)Vm with the Hessenberg-Schur method.
(6) SetX1=X0+WmYmVmT.
(7) ComputeR(X1)=GC-(X1A-FX1).
(8) If∥R(X1)∥F/∥R(X0)∥F≤ϵ stop, otherwise set X0=X1, R(X0)=R(X1) and go to step (3).Remark 4.3.
By choosingm<n, Algorithm 4.2 reduces the original large sparse Sylvester-observer equation (1.2) to a low-dimensional Sylvester-observer equation (4.1). In step (5), we solve this low-dimensional matrix equation by any direct method such as the Hessenberg-Schur method. Then, in step (6) by using relation (4.5), we extend the obtained solution to the solution of the original matrix equation. Also, according to Theorem 4.1, Algorithm 4.2 is the convergence for any initial matrix X0.


## 5. Weighted Block Refinement Method

In this section we discuss a new iterative method based upon a modified block refinement method. The new process uses instead of the Euclidean scalar product another one, denoted by(·,·)D where D is a chosen diagonal matrix. The idea of changing the inner product is to accelerate the convergence of the components of the residual which are far away from zero. To achieve this, an appropriate weight is associated to each term of the inner product. A natural choice of these weights is the entries of the first residual. The following method is based on reduction of A and F to the Hessenberg matrix with the use of weighted block Arnoldi process. Before giving a complete description of the new algorithm, let us define the D-scalar product.Ifu and v are two vectors of Rn, their D-scalar product is (5.1)(u,v)D=vTDu=∑i=1ndiuivi,
where D=diag(d1,d2,…,dn) is a diagonal matrix.This inner product is well defined if and only if the matrixD is positive definite, that is, di>0, for all i∈{1,…,n}.In this case, we can define theD-norm ∥·∥D associated with this inner product by (5.2)‖u‖D=(u,u)D=uTDu=∑i=1ndiui2∀u∈Rn.Theorem 5.1.
LetX0 be the approximate solution obtained by an arbitrary method for the matrix equation (1.2), and let R(X0)=GC-(X0A-FX0) be the corresponding residual. If m(m≪n) steps of the weighted block Arnoldi process by the diagonal matrices D and D̂, respectively, for matrices A and F have been run and Ym∈Rm×m is the solution of the low-dimensional Sylvester equation
(5.3)YmÂm-F̂mTYm=VmTR(X0)Um,
then
(5.4)VmTR(X1)Um=VmTR(X0+D̂VmYmUmTD)Um=0.Proof.
By using the weighted Arnoldi process, we generate the basesUm=[u1,…,um] and Vm=[v1,…,vm] that are, respectively, D-orthonormal and D̂-orthonormal; thus it holds that
(5.5)UmTDUm=Im,VmTD̂Vm=Im,
where Um,Vm∈Rn×m and D,D̂∈Rn×n are two diagonal matrices.
Moreover, the square Hessenberg matricesÂm and F̂m whose nonzero entries are the scalars âij and f̂ij, constructed by the weighted Arnoldi process, can be expressed in the form
(5.6)Âm=UmTDAUm,F̂m=VmTD̂FVm.
Now, we set
(5.7)X1=X0+E0,
where E0=D̂VmYmUmTD and Ym∈Rm×m is the solution of the Sylvester-obsever equation (5.3). Thus, the new residual matrix becomes
(5.8)R(X1)=GC-(X1A-FX1),=GC-((X0+E0)A-F(X0+E0)),=R(X0)-(E0A-FE0),=R(X0)-(D̂VmYmUmTDA-FD̂VmYmUmTD).
Multiplying the above relation on the left by VmT and on the right by Um, we have
(5.9)VmTR(X1)Um=VmTR(X0)Um-(VmTD̂VmYmUmTDAUm-VmTFD̂VmYmUmTDUm).
Now, by using (5.3), (5.5), and (5.6) we get
(5.10)VmTR(X1)Um=VmTR(X0)Um-(YmÂm-F̂mTYm)=0.In order to getYm∈Rm×m, we need to solve the low-dimensional Sylvester equation (5.3). According to the results, we can develop an iterative method for solving of the Sylvester-observer equation. The algorithm  is as follows.Algorithm 5.2 (weighted block refinement (WBR) method).
(1) Start: choose an initial solutionX0, new dimension m lesser than n and a tolerance ϵ.
(2) ComputeR(X0)=GC-(X0A-FX0).
(3) Construct theD-orthonormal basis Um∈Rn×m and D̂-orthonormal basis Vm∈Rn×m by the weighted Arnoldi process, such that
(5.11)Âm=UmTDAUm,F̂m=VmTD̂FVm.
(4) Solve the low-dimensional Sylvester-observer equationYmÂm-F̂mTYm=VmTR(X0)Um with the Hessenberg-Schur method.
(5) SetX1=X0+D̂VmYmUmTD.
(6) Compute residual matrixR(X1)=GC-(X1A-FX1).
(7) If∥R(X1)∥F/∥R(X0)∥F≤ϵ stop, otherwise set X0=X1, R(X0)=R(X1) and go to step (2).Remark 5.3.
By choosingm<n, Algorithm 5.2 reduces the original large sparse Sylvester-observer equation (1.2) to a low-dimensional Sylvester-observer equation (5.3). In step (4), we solve this low-dimensional matrix equation by any direct method such as the Hessenberg-Schur method. Also, according to Theorem 5.1, Algorithm 5.2 is the convergence for any initial matrix X0.


## 6. Numerical Experiments

In this section, we present some numerical examples to illustrate the effectiveness of our new iterative methods for solving large sparse Sylvester-obsever equation. In Examples6.1 and 6.2, we apply Algorithms 2 and 3 for solving matrix equation (1.2). In Example 6.3, we compare the Hessenberg-Schur method described in [1] with our new algorithms for solving large sparse Sylvester-obsever equation. In order to show the efficiency of our algorithms, we choose the matrices A and C arbitrary in these three examples. But in Example 6.4, we use four matrices from MATLAB matrix collection with the large estimation of condition numbers.The initial approximate solution isX0=0n×n. The error is monitored by means of the test(6.1)‖R(X1)‖F‖R(X0)‖F≤ϵ,
with the value of ϵ depending on the examples. The time is given in seconds for all examples. All numerical tests are performed in MATLAB software on a PC with 2.20 GHz with main memory 2 GB.Example 6.1.
For the first test we use two arbitrary matricesA and C. We choose the matrices F and G completely satisfying the controllability requirement of the pair (G,C). Now We apply block refinement method for solving Sylvester-observer equation XA-FX=GC with n=200. Also, we take ϵ=10-6. In Table 1, we report the results for different values of m. In Table 1, the results show that by increasing the values of m and l, the number of iterations decreases. The last column of Table 1 also shows the decreasing of time consumption. Note that the fourth and fifth columns of this table are the errors of the orthogonalization method. The desired accuracy has been chosen as 10-6, but the model works well with any choice of 10-t:
(6.2)A=(101.2.42.82.3.80……01.8101.2.42.82.3.80⋱01.61.8⋱⋱⋱⋱⋱⋱0⋮1.641.6⋱⋱⋱⋱⋱⋱.801.31.64⋱⋱⋱⋱⋱⋱2.3.81.611.3⋱⋱⋱⋱⋱⋱.82.301.61⋱⋱⋱⋱⋱⋱.42.8⋮0⋱⋱⋱⋱⋱⋱1.2.42⋮⋮01.611.31.641.61.8101.200…01.611.31.641.61.810)n×nC=(1,0,0,…,0).Table 1
Implementation of block refinement method to solve the Sylvester equation with different values ofm.
mrℓ∥VmTATVm-Âm∥∥WmTFWm-F̂m∥IterationTime10254.11E-0143.71E-01648212.03202106.32E-0145.97E-01539110.08302151.16E-0137.34E-0153198.26402203.65E-0138.89E-0142256.03502252.01E-0129.39E-0141714.47602304.73E-0121.44E-013672.46Example 6.2.
Consider Example6.1 again. We apply the weighted Block refinement method for solving XA-FX=GC and take ϵ=10-6. In Table 2, we report the results for different values of m.Table 2
Implementation of weighted block refinement method to solve the Sylvester-observer equation with different values ofm.
mrℓ∥UmTDAUm-Âm∥∥VmTD̂FVm-F̂m∥IterationTime10259.34E-0162.27E-0153117.34202101.97E-0152.98E-0151745.89302153.17E-0153.12E-0151064.48402205.22E-0156.05E-015673.37502257.74E-0159.59E-015402.19602308.82E-0151.11E-014110.984Example 6.3.
According to the results in Tables1 and 2, we see that the weighted block refinement method in comparison with block refinement method works better. Now, consider that A and C are the same matrices used in Example 6.1. We apply our two iterative methods with 2 iterations and the Hessenberg-Schur method to solve the Sylvester-observer equation when the dimensions of the matrices are large. Results are shown in Table 3.Table 3
Implementation of new Iterative methods with 2 iterations and the Hessenberg-Schur method for solving the Sylvester equation.
nHessenberg-Schur methodBlock refinementWeighted block refinementCond (F)ErrorTimeErrorTimeErrorTime4003.87E-074.844.12E-093.581.78E-0151.813.17E+0048006.12E-00261.893.55E-00448.247.32E-01332.112.13E+007120011.01139.540.04396.254.13E-01271.431.74E+010160026.29208.2411.01134.321.22E-09112.551.27E+15200071.98347.1439.91242.115.67E-05194.787.12E+021Example 6.4.
In this example we show that the convergence of our proposed algorithms independent of the matrices structure. In this sense, we use four matrices from MATLAB collection for the matrixA. The first matrix is a sparse, random finite element matrix with the condition number 1.84E+03. The second matrix is a symmetric, positive semidefinite (SPD) Toeplitz matrix that is composed of the sum of 800 rank 2 SPD Toeplitz matrices with the condition number 6.18E+04. The third matrix is a row diagonally dominant matrix with the condition number 1.35E+010. The last matrix is a sparse singular, row diagonally dominant matrix resulting from discrediting the Neumann problem with the usual five-point operator on a regular mesh. The estimated condition number is 5.56E+017. For all of these examples, the matrix C is C=sprand(1,n,d), where is a random, sparse matrix with approximately d·n·n uniformly distributed nonzero entries with d=0.5. We choose the matrices F and G completely satisfying the controllability and observability requirement of the pairs (F,G) and (G,C). We apply the Hessenberg-Schur method and weighted block refinement method (Algorithm 5.2) with 3 iterations for solving the Sylvester-observer equation XA-FX=GC. The results are shown in Table 4.
It is also obvious from Table4 that the performance of the weighted block refinement method is much better than that of the Hessenberg-Schur method, specifically for the ill-conditioned matrices.Table 4
Effectiveness of the Hessenberg-Schur method and WBR algorithm  with 3 iterations for randomly generated matrices.
MatrixHessenberg-SchurWBR methodErrorCPU-timeErrorCPU-timeA = gallery (“wathen”, 736, 736)2.19E-0666.241.01E-01649.95A = gallery (“toeppd”, 800, 800, r and (1,800)6.01E-0572.194.98E-01654.81A = gallery (“dorr”, 1200, 0.01)7.93E-01137.1136.12E-014101.02A = gallery (“neumann”, 1156)3.24E-01122.719.21E-01594.29


## 7. Comments and Conclusion

In this paper, we propose two new iterative algorithm s for solving the large sparse Sylvester-obsever matrix equations. The existing projection methods use the Arnoldi process, but the methods described in this paper are based on the weighted block Arnoldi process. Moreover, the refinement process presented in Sections4 and 5 has the capability of improving the results obtained by an arbitrary method. The numerical examples show the efficiency of the proposed schemes.



---

*Source: 184314-2011-09-12.xml*