#### 1. Introduction

When a mathematical model is not available for a dynamical system, it is reasonable to use a data-driven approach for analysis and control of the system. Thus, various data-driven solutions have been recently developed for several stability [1], [2], controllability [3], [4], observability [4], optimal control [5], [6], and so on.

Along this direction of research, the authors have recently developed a data-driven solution to Lyapunov equations in the form of \(PA+A^\top P= -Q\), where \(A\) and \(Q\) are given constant matrices and \(P\) is the unknown [7]. The Lyapunov equations are known to play an important role in control engineering, such as stability analysis and controllability analysis. The solution is based on the so-called *data-driven Lyapunov equation*, which is defined by not the matrix \(A\) but the state trajectories of the system \(\dot{x}(t)=Ax(t)+Bu(t)\). Using this framework, one can solve a Lyapunov equation even when a model is not available to the system.

On the other hand, the number of state trajectories to uniquely determine the solution is \(O(n^2)\) in [7], where \(n\) is the size of the matrix \(A\). In other words, the amount of data quadratically grows with \(n\). This prevents us from applying this method to a large-scale system, and thus it is preferable to reduce the amount of data.

This paper proposes a novel class of data-driven Lyapunov equations such that the number of state trajectories to uniquely determine the solution is \(O(n)\). In the previous method [7], one scalar equation, which corresponds to a part of the Lyapunov equation, is formulated by one state trajectory. On the other hand, in the proposed method, three scalar equations are formulated by any combination of two state trajectories. This idea drastically reduces the amount of data to construct a data-driven Lyapunov equation.

Finally, we note that our data-driven approach has an advantage that prior knowledge can be easily incorporated. In fact, our approach transforms a Lyapunov equation into an equation described by data, which allows us to take prior knowledge into account as additional equations. On the other hand, the model-based approach, which is based on a mathematical model, uniquely provides the solution to a Lyapunov equation from the model. Thus, it is hard to directly incorporate prior knowledge into the problem of solving a Lyapunov equation.

This paper is organized as follows. In Sect. 2, we introduce our problem and review the existing method and the amount of required data [7]. Our solutions are presented in Sects. 3 and 4, respectively. Finally, Sect. 5 concludes this paper.

**Notation:** (i) For a vector \(x \in \mathbb{R}^n\) and symmetric matrix \(P \in \mathbb{R}^{n \times n}\), let \(\mathrm{lvec}(x) \in \mathbb{R}^{1 \times \frac{1}{2}n(n+1)}\) and \(\mathrm{rvec}(P) \in \mathbb{R}^{\frac{1}{2}n(n+1)}\) be the pair of a row vector and column vector, such that \(x^\top P x = \mathrm{lvec}(x)\mathrm{rvec}(P)\). For example,

\[\begin{aligned} \mathrm{lvec}(x) &= \begin{bmatrix} x_1^2 & 2x_{1}x_{2} & 2x_{1}x_{3} & x_2^2 & 2x_{2}x_{3} & x_3^2 \end{bmatrix} , \\ \mathrm{rvec}(P) &= \begin{bmatrix} p_{11} & p_{12} & p_{13} & p_{22} & p_{23} & p_{33} \end{bmatrix}^\top \end{aligned}\] |

for

\[x = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} , P = \begin{bmatrix} p_{11} & p_{12} & p_{13} \\ p_{12} & p_{22} & p_{23} \\ p_{13} & p_{23} & p_{33} \\ \end{bmatrix} ,\] |

where \(x_i \in \mathbb{R}\) is the \(i\)-th component of \(x\), and \(p_{ij} \in \mathbb{R}\) is the \((i,j)\)-th component of \(P\).

(ii) Consider the equations \(f(x)=0\) and \(g(x) =0\), where \(x \in \mathbb{R}^n\). If their solution sets are equal, then the equations are said to be *equivalent*.

#### 2. Problem Formulation and Existing Result

Consider the linear system

\[\begin{equation*} \dot{x}(t) = Ax(t) , \tag{1} \end{equation*}\] |

where \(x(t) \in \mathbb{R}^n\) is the state and \(A \in \mathbb{R}^{n \times n}\) is a Hurwitz matrix.

For the system, a Lyapunov equation is given by

\[\begin{equation*} PA + A^\top P = -Q , \tag{2} \end{equation*}\] |

where \(Q \in \mathbb{R}^{n \times n}\) is a given symmetric matrix and \(P \in \mathbb{R}^{n \times n}\) is the unknown. There exists a unique solution \(P \in \mathbb{R}^{n \times n}\) to (2) if \(A\) is Hurwitz [8].

If the matrix \(A\) is known to us, the solution can be easily obtained; otherwise, it is reasonable to use a data-driven solution to the equation. The problem of finding the solution \(P\) from the behavioral data of (1) is formulated as follows.

**Problem 1:** Consider the system in (1) and the Lyapunov equation in (2). Assume that \(A\) is unknown but Hurwitz. Suppose that the data of state trajectories \(x(t,x_{0i}),\; t \in [0,T]\;(i=1,2,\ldots,q)\) are given, where \(T\) is a positive real number and \(x(t,x_{0i})\) is the state \(x(t)\) of the system (1) for the initial state \(x(0) = x_{0i}\). Find the solution of (2). \(\tag*{◻}\)

For Problem 1, the following solution [7] has been presented. Let us briefly review the existing result. Consider the Lyapunov equation in (2). Let us multiply the both sides of (2) by \(x(t,x_{0i})^\top\) from the left and by \(x(t,x_{0i})\) from the right and integrate them on the interval \([0,T]\). As the result, we obtain the following equation.

\[\begin{multline*} x^\top(T,x_{0i}) P x(T,x_{0i}) - x_{0i}^\top P x_{0i} = -\int_0^T \!\!\! x^\top(t,x_{0i}) Q x(t,x_{0i})dt\\ (i=1,2,\ldots,q) \tag{3} \end{multline*}\] |

This is called the *data-driven Lyapunov equation*.

By using the functions \(\mathrm{lvec}\) and \(\mathrm{rvec}\) introduced at the end of Sect. 1, (3) is equivalently transformed into

\[\begin{multline*} [\mathrm{lvec}(x(T,x_{0i})) - \mathrm{lvec}(x_{0i})]\mathrm{rvec}(P) \\ = -\int_0^T x^\top(t,x_{0i}) Q x(t,x_{0i})dt\\ (i=1,2,\ldots,q) , \tag{4} \end{multline*}\] |

which is represented as

\[\begin{equation*} D\,\mathrm{rvec}(P) = d \tag{5} \end{equation*}\] |

for

\[\begin{align} D &:= \begin{bmatrix} \mathrm{lvec}(x(T,x_{01})) - \mathrm{lvec}(x_{01})\\ \mathrm{lvec}(x(T,x_{02})) - \mathrm{lvec}(x_{02})\\ \vdots\\ \mathrm{lvec}(x(T,x_{0q})) - \mathrm{lvec}(x_{0q})\\ \end{bmatrix} \in \mathbb{R}^{q \times \frac{1}{2}n(n+1)} , \tag{6} \\ d &:= \begin{bmatrix} -\int_0^T x^\top(t,x_{01}) Q x(t,x_{01})dt \\ -\int_0^T x^\top(t,x_{02}) Q x(t,x_{02})dt\\ \vdots\\ -\int_0^T x^\top(t,x_{0q}) Q x(t,x_{0q})dt \end{bmatrix} \in \mathbb{R}^{q}. \tag{7} \end{align}\] |

Equation (5) is in the standard form of a linear equation with respect to the unknown vector \(\mathrm{rvec}(P) \in \mathrm{R}^{\frac{1}{2}n(n+1)}\). Thus, if

\[\begin{equation*} \mathrm{rank}(D) = \frac{1}{2}n(n+1) , \tag{8} \end{equation*}\] |

then (5) has a unique solution. Based on this idea, the solution to Problem 1 in [7] is given as follows.

**Lemma 1:** Consider the equation in (3). If (8) holds, then (3) has a unique solution, and it is equal to the solution to the Lyapunov equation in (2). \(\tag*{◻}\)

From this result, it turns out that the solution to the Lyapunov equation in (2) is obtained by solving a linear equation in (3) defined by the data of state trajectories. Since \(D \in \mathrm{R}^{q \times \frac{1}{2}n(n+1)}\) and (8) imply \(q \geq (1/2)n(n+1)\), the data-driven method [7] requires at least \((1/2)n(n+1)\) state trajectories, i.e., it requires data of \(O(n^2)\).

**Remark:** In Problem 1, \(A\) is assumed to be Hurwitz. Thus, the problem is not solved to determine whether the system in (1) is stable or not. One application of Problem 1 is to analyze the observability of (1), since the solution of (2) for \(Q=C^{\top}C\) is equal to the observability Gramian of (1) [8]. Another application is to obtain a Lyapunov function of a stable system, which can be used, for example, to construct an event-triggered controller [9], [10].

#### 3. Small-Data Solution to Problem 1

##### 3.1 Main Results

In the data-driven Lyapunov equation in (3), each scalar equation is constructed by one state trajectory. However, as will be shown later, three scalar equations can be formulated by any combination of two distinct state trajectories. This fact gives the following result.

**Theorem 1:** (i) Consider the following equation:

\[\begin{multline*} x^\top(T,x_{0i}) P x(T,x_{0j}) - x_{0i}^\top P x_{0j} = -\int_0^T \!\!\!\!\! x^\top(t,x_{0i}) Q x(t,x_{0j})dt\\ (i=1,2,\ldots,q,\;j=i,i+1,\ldots,q) \tag{9} \end{multline*}\] |

which is a modified version of (3). If

\[\begin{equation*} \mathrm{rank}(X_0) = n , \tag{10} \end{equation*}\] |

then (9) has a unique solution, and it is equal to the solution to (2), where \(X_0 := \begin{bmatrix} x_{01} & x_{02} & \cdots & x_{0q} \end{bmatrix}\) is the matrix defined by the initial states of the data \(x(t,x_{0i}),\;t \in [0,T]\;(i=1,2,\ldots,q)\).

(ii) The minimum number of state trajectories to determine the unique solution to (9) is \(n\), i.e., \(O(n)\). \(\tag*{◻}\)

This means that the solution to the Lyapunov equation can be obtained with much less data than the conventional method [7] (Lemma 1). This fact is summarized in Table 1.

**Proof of Theorem 1:** (i) It is the consequence of the following facts.

- Equation (2) has a unique solution if \(A\) is Hurwitz [8].
- Equation (2) is equivalent to
\[\begin{equation*} e^{A^\top T} P e^{AT} - P = -\int_0^T e^{A^\top t} Q e^{At}dt . \tag{11} \end{equation*}\] - If (10) holds, then (11) is equivalent to
\[\begin{equation*} X_0^\top e^{A^\top T} P e^{AT} X_0 - X_0^\top P X_0 = -\int_0^T X_0^\top e^{A^\top t} Q e^{At} X_0dt . \tag{12} \end{equation*}\] - Equation (12) is equivalent to (9).

Let us prove Facts (II)-(IV).

(II) First, we multiply the both sides of (2) by \(e^{A^\top t}\) from the left and by \(e^{A t}\) from the right and integrate them on the interval \([0,T]\). As the result, we obtain (11). Thus, if \(P\) is a solution to (2), then it is also a solution to (11). Meanwhile, the converse holds if both (2) and (11) are unique solutions. Since (11) is a discrete-time Lyapunov equation with the Schur matrix \(e^{AT}\) (because \(A\) is Hurwitz), it has a unique solution. These facts and (I) imply (II).

(III) Under (10), the matrix \(X_0\) has full row rank and \(X_0^\top\) has full column rank, which, together with Lemma 2 in the appendix, proves (III).

(IV) It is trivial from the facts that (9) is the element-wise expression of (12) and \(x(t,x_{0i})=e^{At}x_{0i}\).

(ii) Equation (10) implies \(q \geq n\). This proves (ii). \(\tag*{◻}\)

##### 3.2 Examples

###### 3.2.1 Case with Noise-Free data

Consider the linear system

\[\begin{equation*} \dot{x}(t)= \begin{bmatrix} 0 & 1\\ -2 & -3\\ \end{bmatrix} x(t) \tag{13} \end{equation*}\] |

and the Lyapunov equation in (2) for \(Q := I\), which has the unique solution

\[\begin{equation*} P^*= \begin{bmatrix} 1.250 & 0.250\\ 0.250 & 0.250\\ \end{bmatrix}. \tag{14} \end{equation*}\] |

Then we address Problem 1 for the dataset shown in Fig. 1, where the state trajectories are generated by \(x_{01} = [1 0]^\top\) and \(x_{02} = [0 1]^\top\).

Let us apply Theorem 1 to the problem. In this case, we have \(q=n=2\) and (10) holds for \(x_{01}\) and \(x_{02}\). Then the data-driven Lyapunov equation in (9) is given by

\[\begin{align} &\nonumber \begin{bmatrix} 0.600 & -0.465 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.600 \\ -0.465 \end{bmatrix}\\ & - \begin{bmatrix} 1.000 & 0.000 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.000 \\ 0.000 \end{bmatrix} = -0.885 , \tag{15} \\ &\nonumber \begin{bmatrix} 0.233 & -0.097 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.600 \\ -0.465 \end{bmatrix}\\ & - \begin{bmatrix} 0.000 & 1.000 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.000 \\ 0.000 \end{bmatrix} = -0.106 , \tag{16} \\ &\nonumber \begin{bmatrix} 0.233 & -0.097 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.233 \\ -0.097 \end{bmatrix}\\ & - \begin{bmatrix} 0.000 & 1.000 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.000 \\ 1.000 \end{bmatrix} = -0.191 , \tag{17} \end{align}\] |

which provides

\[\begin{equation*} P= \begin{bmatrix} 1.250 & 0.250\\ 0.250 & 0.250\\ \end{bmatrix} . \tag{18} \end{equation*}\] |

This agrees with the solution to the original Lyapunov equation, i.e., \(P^*\) in (14).

###### 3.2.2 Case with Noisy Data

We consider the case where noisy data is given as shown in Fig. 2. The data are the state trajectories in Fig. 1 corrupted by Gaussian noise with a mean of 0 and a variance of 0.05.

The data-driven Lyapunov equation in (9) is given by

\[\begin{align} &\nonumber \begin{bmatrix} 0.551 & -0.444 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.551 \\ -0.444 \end{bmatrix}\\ & - \begin{bmatrix} 0.943 & 0.036 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.943 \\ 0.036 \end{bmatrix} = -0.878 , \tag{19} \\ &\nonumber \begin{bmatrix} 0.227 & -0.098 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.551 \\ -0.444 \end{bmatrix}\\ & - \begin{bmatrix} -0.049 & 0.960 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.943 \\ 0.036 \end{bmatrix} = -0.110 , \tag{20} \\ &\nonumber \begin{bmatrix} 0.227 & -0.098 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.227 \\ -0.098 \end{bmatrix}\\ & - \begin{bmatrix} -0.049 & 0.960 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} -0.049 \\ 0.960 \end{bmatrix} = -0.199 . \tag{21} \end{align}\] |

which provides

\[\begin{equation*} P= \begin{bmatrix} 1.300 & 0.317\\ 0.317 & 0.306\\ \end{bmatrix} \tag{22} \end{equation*}\] |

It turns out that this approximates the true value \(P^*\) in (14).

On the other hand, the existing method [7] requires more data, i.e., three state trajectories, because three scalar equations are constructed by three state trajectories, as shown in (3). For example, for the three state trajectories in Figs. 2 and 3, the data-driven Lyapunov equation in (3) is given by (19), (21), and

\[\begin{align} &\nonumber \begin{bmatrix} 1.036 & -0.683 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.036 \\ -0.683 \end{bmatrix}\\ &- \begin{bmatrix} 0.969 & 1.977 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.969 \\ 1.977 \end{bmatrix} = -2.066, \tag{23} \end{align}\] |

which provides

\[\begin{equation*} P = \begin{bmatrix} 1.384 & 0.229\\ 0.229 & 0.305 \end{bmatrix}. \tag{24} \end{equation*}\] |

Although the solution is derived by more data, the accuracy is similar to the proposed method.

Finally, we show that the proposed method provides better performance with the data as shown in Figs. 2 and 3. In this case, the data-driven Lyapunov equation is given by (19), (20), (21), (23),

\[\begin{align} &\nonumber \begin{bmatrix} 0.551 & -0.444 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.036 \\ -0.683 \end{bmatrix}\\ & - \begin{bmatrix} 0.943 & 0.036 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.969 \\ 1.977 \end{bmatrix} = -1.085 , \tag{25} \\ &\nonumber \begin{bmatrix} 0.227 & -0.098 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.036 \\ -0.683 \end{bmatrix}\\ & - \begin{bmatrix} -0.049 & 0.960 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.969 \\ 1.977 \end{bmatrix} = -0.494 . \tag{26} \end{align}\] |

The least-squares solution of these equations is

\[\begin{equation*} P= \begin{bmatrix} 1.289 & 0.252\\ 0.252 & 0.278\\ \end{bmatrix} , \tag{27} \end{equation*}\] |

which is closer to \(P^*\) than (24).

###### 3.2.3 Application to Input Node Design of a Large-Scale System

Consider the model of a 100-order gene regulatory network system as shown in Fig. 4, which is developed in [11]. For the system, we address the problem of finding a gene to which the control input is applied so that the trace of the controllability Gramian is maximized.

As shown in [7], the solution is the index of the maximum diagonal element of the solution \(P\) of (2) for \(Q = I\). Thus, let us solve the corresponding data-driven Lyapunov equation. Our result provides gene 34 as the solution with the data of 100 state trajectories. On the other hand, the existing method [7] requires 5050 state trajectories to obtain the solution. This suggests that our method is more practical than the existing method.

#### 4. Extension to the Case with Discrete-Time Data

##### 4.1 Main Results

Suppose that the state trajectory \(x(t,x_{0})\) is sampled for the period \(h\) and the discrete data \(\bar{x}_k := x(kh,x_{0})\) is given. Then we consider the discrete-time Lyapunov equation

\[\begin{equation*} \bar{A}^\top P \bar{A} - P = -Q , \tag{28} \end{equation*}\] |

where \(\bar{A}:= e^{Ah}\), \(Q \in \mathbb{R}^{n \times n}\) is a given symmetric matrix, and \(P \in \mathbb{R}^{n \times n}\) is the unknown. There exists a unique solution \(P \in \mathbb{R}^{n \times n}\) to (28) if \(A\) is Hurwitz, i.e., \(\bar{A}\) is Schur [8].

Now, we consider the following problem.

**Problem 2:** Consider the system in (1) and the discrete-time Lyapunov equation in (28). Assume that \(\bar{A}\) is unknown and Schur, but the sampled data \(\bar{x}_k := x(kh,x_{0})\;(k=0,1,\dots,s)\) of the state trajectories of (1) are given, where \(h \in (0,\infty)\) is the sampling period. Find the solution of (28). \(\tag*{◻}\)

To derive a solution, we represent the data \(\bar{x}_k\) \((k=0,1,\dots,s-1)\) in the following matrix form:

\[\begin{equation*} \bar{X}_0 := \begin{bmatrix} \bar{x}_0 & \bar{x}_1 & \cdots & \bar{x}_{s-1} \end{bmatrix}. \tag{29} \end{equation*}\] |

Then a solution to Problem 2 is obtained as follows.

**Theorem 2:** (i) Consider the following equation:

\[\begin{multline*} \bar{x}_k^\top P \bar{x}_l - \bar{x}_{k-1}^\top P \bar{x}_{l-1} = -\bar{x}_{k-1}^\top Q \bar{x}_{l-1}\\ (k=1,2,\ldots,s,\;l=k,k+1,\ldots,s). \tag{30} \end{multline*}\] |

If

\[\begin{equation*} \mathrm{rank}(\bar{X}_0) = n, \tag{31} \end{equation*}\] |

then (30) has a unique solution, and it is equal to the solution to (28).

(ii) The minimum number of the sampled data to determine the unique solution to (30) is \(n+1\), i.e., \(O(n)\). \(\tag*{◻}\)

**Proof of Theorem 2:** (i) It is the consequence of the following facts.

- Equation (28) has a unique solution if \(A\) is Hurwitz [8].
- If (31) holds, then (28) is equivalent to
\[\begin{equation*} \bar{X}_0^\top \bar{A}^\top P \bar{A} \bar{X}_0 - \bar{X}_0^\top P \bar{X}_0 = -\bar{X}_0^\top Q \bar{X}_0 . \tag{32} \end{equation*}\] - Equation (32) is equivalent to (30).

Let us prove Facts (II) and (III).

(II) Under (31), the matrix \(\bar{X}_0\) has full row rank and \(\bar{X}_0^\top\) has full column rank, which, together Lemma 2 in the appendix, proves (II).

(III) It is trivial from the facts that (30) is the element-wise expression of (32) and \(\bar{x}_k =\bar{A}\bar{x}_{k-1}\).

(ii) Equation (31) implies \(s \geq n\). This proves (ii). \(\tag*{◻}\)

##### 4.2 Example

Consider the linear system in (13) and the discrete-time Lyapunov equation in (28) for \(Q := I\) and \(h=0.1\). The solution is given by

\[\begin{equation*} P^*= \begin{bmatrix} 13.000 & 2.508 \\ 2.508 & 3.050 \end{bmatrix} . \tag{33} \end{equation*}\] |

Then we consider Problem 2 for the dataset composed of

\[\begin{equation*} \bar{x}_{0} = \begin{bmatrix} 1.000 \\ 1.000 \end{bmatrix},\, \bar{x}_{1} = \begin{bmatrix} 1.077 \\ 0.560 \end{bmatrix},\, \bar{x}_{2} = \begin{bmatrix} 1.116 \\ 0.225 \end{bmatrix}. \tag{34} \end{equation*}\] |

Let us apply Theorem 2 to the problem. In this case, (31) holds for \(\bar{x}_{0}\) and \(\bar{x}_{1}\). Then the data-driven discrete-time Lyapunov equation in (30) is given by

\[\begin{align} &\nonumber \begin{bmatrix} 1.077 & 0.560 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.077 \\ 0.560 \end{bmatrix}\\ & - \begin{bmatrix} 1.000 & 1.000 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.000 \\ 1.000 \end{bmatrix} = -2.000 , \tag{35} \\ &\nonumber \begin{bmatrix} 1.116 & 0.225 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.077 \\ 0.560 \end{bmatrix}\\ & - \begin{bmatrix} 1.077 & 0.560 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.000 \\ 1.000 \end{bmatrix} = -1.637 , \tag{36} \\ &\nonumber \begin{bmatrix} 1.116 & 0.225 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.116 \\ -0.225 \end{bmatrix}\\ & - \begin{bmatrix} 1.077 & 0.560 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.077 \\ 0.560 \end{bmatrix} = -1.474 , \tag{37} \end{align}\] |

which provides

\[\begin{equation*} P= \begin{bmatrix} 13.000 & 2.508 \\ 2.508 & 3.050 \end{bmatrix} . \tag{38} \end{equation*}\] |

This agrees with the solution to the original discrete-time Lyapunov equation, i.e., \(P^*\) in (33).

#### 5. Conclusions

In this paper, we have developed a new class of data-driven Lyapunov equations, which are defined by data of amount \(O(n)\). Since the amount of data is much less than that of the existing result [7], i.e., \(O(n^2)\), this result broadens the applicability of the framework of data-driven Lyapunov equations.

On the other hand, it is still open whether \(O(n)\) is the minimum or not to determine the unique solution to a Lyapunov equation in a data-driven manner. In the future, we plan to answer to this question.

#### Acknowledgments

This research was supported by JST Moonshot R&D #JPMJMS2021.

#### References

[1] U.S. Park and M. Ikeda, “Stability analysis and control design of lti discrete-time systems by the direct use of time series data,” Automatica, vol.45, no.5, pp.1265-1271, 2009.

CrossRef

[2] F. Zhang and D. Söffker, “A data-driven online stability monitoring method for unknown discrete-time nonlinear systems,” 50th IEEE Conference on Decision and Control and European Control Conference, pp.7356-7361, 2011.

CrossRef

[3] V.K. Mishra, I. Markovsky, and B. Grossmann, “Data-driven tests for controllability,” IEEE Control Syst. Lett., vol.5, no.2, pp.517-522, 2021.

CrossRef

[4] Z. Wang and D. Liu, “Data-based controllability and observability analysis of linear discrete-time systems,” IEEE Trans. Neural Netw., vol.22, no.12, pp.2388-2392, 2011.

CrossRef

[5] D. Vrabie, O. Pastravanu, M. Abu-Khalaf, and F. Lewis, “Adaptive optimal control for continuous-time linear systems based on policy iteration,” Automatica, vol.45, no.2, pp.477-484, 2009.

CrossRef

[6] Y. Jiang and Z.P. Jiang, “Computational adaptive optimal control for continuous-time linear systems with completely unknown dynamics,” Automatica, vol.48, no.10, pp.2699-2704, 2012.

CrossRef

[7] I. Banno, S. Azuma, R. Ariizumi, T. Asai, and J. Imura, “Data-driven estimation and maximization of controllability Gramians,” 2021 60th IEEE Conference on Decision and Control, pp.5053-5058, 2021.

CrossRef

[8] B.N. Datta, Numerical Methods for Linear Control Systems: Design and Analysis, Academic Press, 2004.

CrossRef

[9] I. Banno, S. Azuma, R. Ariizumi, and T. Asai, “Data-driven sparse event-triggered control of unknown systems,” 2021 American Control Conference, pp.3383-3388, 2021.

CrossRef

[10] W.P.M.H. Heemels, K.H. Johansson, and P. Tabuada, “An introduction to event-triggered and self-triggered control,” IEEE Conference on Decision and Control, pp.3270-3285, 2012.

CrossRef

[11] X. Shen, M. Morishita, J. Imura, M. Oku, and K. Aihara, “Low-sample-size data-driven re-stabilization of gene network systems,” 10th IFAC Symposium on Robust Control Design ROCOND 2022, vol.55, no.25, pp.241-246, 2022.

CrossRef

[12] R. Piziak and P.L. Odell, “Full rank factorization of matrices,” Mathematics Magazine, vol.72, no.3, pp.193-201, 1999.

CrossRef

#### Appendix: On Equivalence of Two Equations

For a matrix-valued function \(f: \mathbb{R}^{n \times n} \to \mathbb{R}^{n \times n}\) and matrices \(A \in \mathbb{R}^{m \times n}\), \(B \in \mathbb{R}^{n \times l}\), we consider the two equations

\[\begin{align} &f(X) = 0_{n \times n} , \tag{A$\cdot $1} \\ &Af(X)B = 0_{m \times l} , \tag{A$\cdot $2} \end{align}\] |

where \(0_{n \times n}\) is an \(n \times n\) zero matrix and \(0_{m \times l}\) is similarly defined. Let \(S \subset \mathbb{R}^{n \times n}\) be the set of the solutions to (A\(\cdot\)1) and \(\tilde{S} \subset \mathbb{R}^{n \times n}\) be that to (A\(\cdot\)2), i.e.,

\[\begin{align} S &:= \{ X \in \mathbb{R}^{n \times n} \mid f(X) = 0_{n \times n} \} , \tag{A$\cdot $3} \\ \tilde{S} &:= \{ X \in \mathbb{R}^{n \times n} \mid Af(X)B = 0_{m \times l} \} . \tag{A$\cdot $4} \end{align}\] |

**Lemma 2:** Consider the sets in (A\(\cdot\)3) and (A\(\cdot\)4). If \(A\) has full column rank and \(B\) has full row rank, then \(S=\tilde{S}\). \(\tag*{◻}\)

**Proof:** The relation \(S \subset \tilde{S}\) is trivial from \(X \in S \Rightarrow f(X) = 0_{n \times n} \Rightarrow Af(X)B = 0_{m \times l} \Rightarrow X \in \tilde{S}\).

Next, we prove \(S \supset \tilde{S}\). If \(X \in \tilde{S}\), then

\[\begin{equation*} Af(X)B = 0_{m \times l} . \tag{A$\cdot $5} \end{equation*}\] |

If \(A\) has full column rank, then \(A^\top A\) is non-singular [12]. Similarly, if \(B\) has full row rank, then \(B B^\top\) is non-singular. Therefore, multiplying the both sides of (A\(\cdot\)5) by \((A^\top A)^{-1}A^\top\) from the left and by \(B^\top(B B^\top)^{-1}\) from the right, we get

\[\begin{equation*} (A^\top A)^{-1}A^\top Af(X)B B^\top(B B^\top)^{-1} = 0_{n \times n} , \tag{A$\cdot $6} \end{equation*}\] |

which is equal to

\[\begin{equation*} f(X) = 0_{n \times n} , \tag{A$\cdot $7} \end{equation*}\] |

i.e., \(X \in S\). This proves \(S \supset \tilde{S}\).

Since \(S \subset \tilde{S}\) and \(S \supset \tilde{S}\), we obtain \(S=\tilde{S}\). \(\tag*{◻}\)