# Markov Parameter Identification via Chebyshev Approximation

Jiayun Li  Yilin Mo  Department of Automation and BNRist, Tsinghua University,
Beijing, P.R.China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper proposes an identification algorithm for Single Input Single Output
(SISO) Linear Time-Invariant (LTI) systems. In the noise-free setting, where
the first $T$ Markov parameters can be precisely estimated, all Markov
parameters can be inferred by the linear combination of the known $T$ Markov
parameters, of which the coefficients are obtained by solving the uniform
polynomial approximation problem, and the upper bound of the asymptotic
identification bias is provided. For the finite-time identification scenario,
we cast the system identification problem with noisy Markov parameters into a
regularized uniform approximation problem. Numerical results demonstrate that
the proposed algorithm outperforms the conventional Ho-Kalman Algorithm for
the finite-time identification scenario while the asymptotic bias remains
negligible.

###### keywords:

Stochastic system identification, Identification for control, Linear systems,
Time-invariant systems.

††thanks: This work is supported by National Natural Science Foundation of
China under grant No. 62192752.

##  1 Introduction

Linear Time-Invariant (LTI) systems are an important class of models in many
technical fields, e.g., industry (Schroeck et al., 2001), automobile
(Kritayakirana and Gerdes, 2012), robotics (Şen and Kalyoncu, 2019) and so on.
Although conventional control methods for the LTI systems, which may not need
accurate system model, are proved to be effective (e.g. the PID controller),
for tasks requiring high agility and performance, precise models become a
necessity. Therefore, the system identification problem has caught great
attention. While classical identification research mainly focuses on the
identifiability of the system, as well as the asymptotic performance of
certain identification methods, e.g., the Ordinary Least Square (OLS) based
method, there has recently been growing attention on the finite-time
performance of the identification algorithms, which can be further divided
into two categories, depending on whether the state of the system is directly
obtained.

When the state of the system is precisely observed, the identification of
system parameters using the classic Least Square (LS) method is proved to be
near optimal, where several recent works are devoted to deriving concentration
analysis on the identification error of the finite-time LS algorithm
(Faradonbeh et al., 2017; Shirani Faradonbeh et al., 2018; Rantzer, 2018;
Simchowitz et al., 2018; Sarkar and Rakhlin, 2019). Besides the LS methods,
Wagenmaker and Jamieson (2020) propose an active learning method for system
identification, and a theoretical upper bound on the identification error of
system parameters is provided.

When the state of the system cannot be directly measured, the system
identification problem becomes “more challenging” (Oymak and Ozay, 2019). A
great line of work aims to address the identification of a finite number of
Markov parameters from the limited sample trajectories and analyze the
performance of specific identification algorithms. Simchowitz et al. (2019)
propose to identify finite Markov parameters using the pre-filtered least
square method, and an upper bound on the identification error is further
derived. Furthermore, Zheng and Li (2021) leverage the LS method to identify
finite Markov parameters from the sample trajectories of both stable and
unstable systems, where the performance of the LS method is further provided.
In addition, another line of work approximates the transfer function of the
system by recovering Finite Impulse Response (FIR) (Ljung and Yuan, 1985;
Ljung and Wahlberg, 1992; Goldenshluger, 1998) or truncated Infinite Impulse
Response (IIR) (Yin et al., 2021; Iannelli et al., 2021) from the sampled
data. The performance of the FIR methods are guaranteed by Tu et al. (2017)
when only a limited amount of data is available. However, leveraging a finite
number of Markov parameters to approximate either the input-output
relationship of the system or the transfer function suffers from truncation
error (Tu et al., 2017).

On the other hand, several methods are proposed to address the above problem
by identifying system parameters of the state-space model using finite Markov
parameters. Oymak and Ozay (2019) propose the finite-time Ho-Kalman algorithm,
which recovers a balanced realization of the original system by conducting SVD
decomposition to the Hankel matrix formed by finite Markov parameters, and
they further provide a theoretical upper bound on the identification error of
system parameters. Besides the LS-based methods, subspace methods are another
common approach to recover system parameters from Markov parameters, and
Tsiamis and Pappas (2019) derive an upper bound on the identification error of
subspace methods. There is also a line of work that focuses on the
identification performance of the classic MOESP algorithm (Ikeda and Oku,
2015; Chiuso and Picci, 2004), which is a well-known subspace method.

However, our initial work (Li et al., 2022) provides theoretical analysis on
the sample complexity of the Ho-Kalman algorithm with respect to system
dimension, and prove that the algorithm is ill-conditioned for high-
dimensional systems. Besides, the ill-conditionedness of subspace methods is
also observed both in theory and practice, e.g., see Chiuso and Picci (2004);
Hachicha et al. (2014). As a result, although the above algorithms are proved
to be asymptotically unbiased, the finite-time performance may be poor,
especially for identification problems of high-dimensional systems. Numerical
results further demonstrate the finite-time performance of the Ho-Kalman
algorithm in Section 5.

Besides the specific identification algorithms, our initial works also prove
that the ill-conditionedness is actually rooted in the identification problem
itself, and that the identification problem is ill-conditioned for high-
dimensional systems using any unbiased identification algorithms (Sun et al.,
2022; Sun and Mo, 2022; Li et al., 2022).

This result motivates us to avoid identifying state-space parameters, and
instead derive a biased Markov parameter identification algorithm. We
circumvent ill-conditionedness by minimizing the mean squared identification
error of Markov parameters using Chebyshev approximation. Taking notice of the
fact that infinite Markov parameters contain full information about the input-
output relationship of the original system, we choose to recover all Markov
parameters using the proposed algorithm. Furthermore, an upper bound on the
asymptotic identification bias is provided, which indicates that the bias is
negligible compared to the impact of noise in the original data.

The main contribution of this paper is as follows:

  * •

We formulate the Markov parameter identification as a uniform polynomial
approximation problem, assuming the first $T$ Markov parameters can be
precisely observed.

  * •

When only noisy Markov parameters from finite sample trajectories are
available, we further cast the identification problem into a regularized
polynomial approximation problem.

  * •

An upper bound on the asymptotic identification bias of the proposed algorithm
is provided, which is negligible compared to the impact of noise in the
original data.

The rest of the paper is organized as follows: In Section 2, we formulate the
identification problem in the SISO case. In addition, Section 3 illustrates
our algorithm in the noise-free setting, and a theoretical upper bound on the
identification bias is provided. Section 4 analyzes the mean squared
identification error in the noisy scenario, by which we formulate the finite-
time identification problem as regularized uniform approximation. Numerical
results are depicted in Section 5, which verifies the finite-time performance
of our algorithm, and Section 6 concludes the paper.

##  2 Problem Formulation

Consider a Single Input Single Output (SISO) discrete-time LTI system with the
following state-space realization:

$\displaystyle x_{t+1}$ $\displaystyle=Ax_{t}+Bu_{t}+w_{t},$ (1)
$\displaystyle y_{t}$ $\displaystyle=Cx_{t}+v_{t},$

where $x_{t}\in\mathbb{R}^{n},u_{t}\in\mathbb{R},y_{t}\in\mathbb{R}$ are the
system state, input and output at time $t$ respectively, and $n$ is the
dimension of the system. $\\{w_{t}\\}$ and $\\{v_{t}\\}$ are independent and
identically distributed (i.i.d.) Gaussian noise with zero mean and covariance
$Q$, $R$, (with $Q,R\geq 0$), respectively. $A,B$ and $C$ are system
parameters with proper dimensions.

The overarching goal of this paper is to identify all Markov parameters
$\\{H_{t}\\}_{t=1}^{\infty}$, which contain full information about the input-
output relationship of the original system, from the first $T$ exact Markov
parameters $\\{H_{t}\\}_{t=1}^{T}$ in Section 3 or from the estimated Markov
parameters $\\{\tilde{H}_{t}\\}_{t=1}^{T}$ in Section 4.

We further make the following assumptions:

###### Assumption 1

System matrix $A$ is real diagonalizable, and all eigenvalues of $A$ lie on
the real interval $[-\rho,\rho]$, where $\rho>0$ is known.

Since similarity transformation preserves the input-output relationship of the
original system, for convenience of notation, we assume that the system takes
the diagonal canonical form, i.e.,
$A=\text{diag}\\{\lambda_{1},\lambda_{2},\cdots,\lambda_{n}\\}$, where
$\lambda_{i},i=1,\cdots,n$ are the eigenvalues of system matrix $A$. $B$ is an
all one column vector, i.e.,
$B=\begin{bmatrix}1&1&\cdots&1\end{bmatrix}^{\top}$.
$C=\begin{bmatrix}c_{1}&c_{2}&\cdots&c_{n}\end{bmatrix}$, and we further
define the modified system energy as

$C_{m}\triangleq\|C\|_{1}.$ (2)

Note that though the roots $\lambda_{i}$ are denoted distinctly, the
eigenvalues with different indices can have the same value.

###### Remark 1

The identification of state-space parameters $A,B$ and $C$ and the estimation
of all Markov parameters are two independent schemes that are able to recover
full input-output relationship of the original system. Previous identification
algorithms mainly consider the estimation of state-space parameters, while we
choose to directly identify all Markov parameters of the LTI system.

###### Remark 2

For simplicity and convenience of notation, we only discuss the identification
of SISO systems in this paper. However, it is worth noting that the proposed
algorithm is directly applicable to Multiple-Input Multiple-Output (MIMO)
systems without significant modification.

###### Remark 3

The identification of systems with complex poles is left for future research.

##  3 Identification algorithm in the noise-free case

This section considers the identification problem in the noise-free case,
where we make the following assumption:

###### Assumption 2

The first $T$ Markov parameters $\\{H_{t}\\}_{t=1}^{T}$ are perfectly known.

The proposed algorithm aims to identify $H_{k}$ through the linear combination
of the first $T$ Markov parameters $\\{H_{t}\\}_{t=1}^{T}$, where $k>T$:

$\hat{H}_{k}=\sum_{t=1}^{T}\alpha_{t-1}H_{t},$ (3)

such that the mean squared identification error is minimized:

$\epsilon\triangleq|H_{k}-\hat{H}_{k}|^{2}=\left|H_{k}-\sum_{t=1}^{T}\alpha_{t-1}H_{t}\right|^{2}.$
(4)

According to Cayley-Hamilton Theorem, for each $k$, $H_{k}$ can be expressed
as a linear combination of the first $n$ Markov parameters. Thus, $\epsilon=0$
can be reached by choosing system-specific coefficients $\alpha_{t}$ when
system parameters are known, which is not satisfied by the identification
problem. On the other hand, noticing that the $k$-th Markov parameter of
system (1) has the following relationship with state-space parameters:

$H_{k}=CA^{k-1}B=\sum_{j=1}^{n}c_{j}\lambda_{j}^{k-1},$ (5)

we instead formulate the minimization of $\epsilon$ in (4) as an approximation
problem below:

###### Problem 1

$\boldsymbol{\alpha}^{*}=\text{argmin}_{\\{\alpha_{t}\\}_{t=0}^{T-1}}\left\|\lambda^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}\lambda^{t}\right\|_{\infty,\rho},$
(6)

in which $\|\cdot\|_{\infty,\rho}$ denotes the uniform norm of functions,
i.e., $\|f\|_{\infty,\rho}=\sup\\{|f(\lambda)|:\lambda\in[-\rho,\rho]\\}$,
$\boldsymbol{\alpha}^{*}\triangleq\begin{bmatrix}\alpha^{*}_{0}&\alpha^{*}_{1}&\cdots&\alpha^{*}_{T-1}\end{bmatrix}^{\top}$
is the optimal solution.

Therefore, $\boldsymbol{\alpha}^{*}$ becomes a universal solution without the
dependence on system parameters $A,B,C$ and $n$, and is only associated with
$T,\rho$ and $k$.

Furthermore, the estimation of the $k$-th Markov parameter can be calculated
as:

$\displaystyle\hat{H}_{k}=\sum_{t=1}^{T}\alpha_{t-1}^{*}H_{t}.$ (7)

###  3.1 Performance Analysis

In this subsection, we analyze the performance of the above algorithm, where
the proof of the following theorem is provided in Appendix A.

###### Theorem 1

The mean squared identification error of the $k$-th Markov parameter of system
(1) using the first $T$ Markov parameters $\\{H_{t}\\}_{t=1}^{T}$ with the
algorithm in (6) has the following upper bound, $\forall k>T$:

$\displaystyle\left|H_{k}-\hat{H}_{k}\right|^{2}\leq
C_{m}^{2}\rho^{2k-2}\min\left(4\exp\left(-\frac{2(T-1)^{2}}{k-1}\right),\frac{1}{4}\right).$
(8)

where $\hat{H}_{k}=\sum_{t=1}^{T}\alpha_{t-1}^{*}H_{t}$ is the estimation of
$H_{k}$, the coefficients $\boldsymbol{\alpha}^{*}$ are obtained by (6) and
$C_{m}$ is the modified system energy defined in (2).

###### Remark 3.1

When the system is strictly stable and $\rho<1$, the finite $\mathcal{L}_{2}$
identification error of all Markov parameters can be obtained by putting all
estimation of Markov parameters to be $0$, i.e., $\hat{H}_{k}=0,\forall k>T$
(which shall be referred to as the truncation method in the following):

$\displaystyle\sqrt{\sum_{k=1}^{\infty}|H_{k}-\hat{H}_{k}|^{2}}$
$\displaystyle=\sqrt{\sum_{k=T+1}^{\infty}\left|\sum_{j=1}^{n}c_{j}\lambda_{j}^{k-1}\right|^{2}}$
(9)
$\displaystyle\leq\sqrt{\sum_{k=T}^{\infty}C_{m}^{2}\rho^{2k}}=\frac{C_{m}\rho^{T}}{\sqrt{1-\rho^{2}}}.$

Furthermore, Theorem 1 reveals that the upper bound on the identification
error of each Markov parameter is better than the truncation method:

$\displaystyle|H_{k}-\hat{H}_{k}|^{2}$ $\displaystyle\leq
C_{m}^{2}\rho^{2k-2}\min\left(4\exp\left(-\frac{2(T-1)^{2}}{k-1}\right),\frac{1}{4}\right)$
(10) $\displaystyle\leq C_{m}^{2}\rho^{2k-2},\forall k>T,$

which indicates that the $\mathcal{L}_{2}$ identification error of the
proposed algorithm is also finite and is strictly better than
$\frac{C_{m}\rho^{T}}{\sqrt{1-\rho^{2}}}$. Moreover, the exponential term
$\exp\left(-\frac{2(T-1)^{2}}{k-1}\right)$ in (8) guarantees that the upper
bound on the identification error of the first few unknown Markov parameters,
which dominates the $\mathcal{L}_{2}$ error of the truncation method, is
drastically reduced by our algorithm.

On the other hand, when the system is unstable, since the value of Markov
parameter $H_{k}$ explodes to infinity with the increase of the index $k$,
only a finite number of Markov parameters are worth identifying, where the
identification error is also bounded by Theorem 1.

###### Remark 3.2

It is worth noting that the identification bias can be efficiently reduced by
increasing the number of known Markov parameters $T$, which is further
illustrated considering the following three cases:

  * •

For a fixed index $k$, the identification error of $H_{k}$ attenuates
superlinearly w.r.t. the number of known Markov parameters $T$:

$|H_{k}-\hat{H}_{k}|^{2}\leq
4C_{m}^{2}\rho^{2k-2}\exp\left(-\frac{2(T-1)^{2}}{k-1}\right).$ (11)

  * •

When the system is strictly stable and $\rho<1$, the identification bias has
the following supremum among all Markov parameters:

$\sup_{k}|H_{k}-\hat{H}_{k}|^{2}\leq 4C_{m}^{2}\exp\left(-4(T-1)\sqrt{\log
1/\rho}\right),$ (12)

which attenuates exponentially w.r.t. $T$.

  * •

When the system is strictly stable and $\rho<1$, in order to obtain
identification accuracy $\delta^{2}$, by (12), the number of known Markov
parameters required satisfies

$\displaystyle
T\sim\mathcal{O}\left(\frac{\log(2C_{m})-\log\delta}{\sqrt{\log(1/\rho)}}\right),$
(13)

and when $\rho$ is close to $1$,

$T\sim\mathcal{O}\left(\frac{\log(2C_{m})-\log\delta}{\sqrt{1-\rho}}\right).$
(14)

The above result indicates that the number of known Markov parameters $T$
required for a specific system increases at an efficient speed of
$T\sim\mathcal{O}\left(\log\left(\frac{1}{\delta}\right)\right).$

##  4 Identification in the noisy scenario

In this section, we address the identification problem where only a finite
amount of data is available, and the estimate of Markov parameters
$\\{\tilde{H}_{t}\\}_{t=1}^{T}$ is contaminated by noise. We make the
following assumption on the obtained data:

###### Assumption 3

The sample trajectories are collected episodically from the system, where each
episode is reset to steady state and is independent from each other.

###  4.1 Mean Squared Identification Error

We first illustrate one method to recursively identify $\tilde{H}_{t}$, where
$t\leq T$, using the sample trajectories obtained from the episodic
experiments.

In each episode of the experiment, we inject the unit impulse function to the
system, i.e.,

$u_{t}=\left\\{\begin{aligned} 1,\quad&t=0,\\\ 0,\quad&t\geq
1,\end{aligned}\right.$ (15)

and we denote the $t$-th output in the $\ell$-th episode as
$y_{t}^{(\ell)},t=1,\cdots,T$. The estimation of the $t$-th Markov parameter
can be updated recursively by:

$\tilde{H}_{t}^{(\ell)}=\tilde{H}_{t}^{(\ell-1)}+\frac{1}{\ell}\left(y_{t}^{(\ell)}-\tilde{H}_{t}^{(\ell-1)}\right),\tilde{H}_{t}^{(0)}=0.$
(16)

The following theorem quantifies the mean squared identification error of
$H_{k}$ using the linear combination of the first $T$ Markov parameters
$\\{\tilde{H}_{t}^{(N)}\\}_{t=1}^{T}$ estimated from $N$ sample trajectories,
where the proof is reported in Appendix B.

###### Theorem 2

Let $\hat{H}_{k}=\sum_{t=1}^{T}\alpha_{t-1}\tilde{H}_{t}^{(N)}$ be the
approximation of the $k$-th Markov parameter $H_{k}(k>T)$. Then, the mean
squared approximation error of the $k$-th Markov parameter has the following
upper bound:

$\mathbb{E}\left(H_{k}-\hat{H}_{k}\right)^{2}\leq
C_{m}^{2}\left\|\lambda^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}\lambda^{t}\right\|_{\infty,\rho}^{2}+\frac{\Sigma}{N}\|\boldsymbol{\alpha}\|_{1}^{2},$
(17)

where $\Sigma$ is the variance of the estimation of Markov parameters:

$\Sigma=\mathbb{E}\left(H_{t}-y_{t}^{(\ell)}\right)^{2}=CPC^{\top}+R,\text{
where }P=APA^{\top}+Q,$ (18)

$\forall\ell=1,\cdots,N,t=1,\cdots,T$, $A,C,Q$ and $R$ are system parameters
defined in Section 2.

###### Remark 4.1

The assumption that each experiment is reset to steady state is only for
simplicity of the analysis in Theorem 2. Generally, the proposed algorithm is
able to work without this assumption.

###### Remark 4.2

Although $\Sigma$ cannot be calculated by (18) since the system parameters are
unknown, it can be estimated from data using the following equation:

$\hat{\Sigma}=\frac{1}{TN-1}\sum_{t=1}^{T}\sum_{\ell=1}^{N}\left(y_{t}^{(\ell)}-\bar{y}_{t}\right)^{2},$
(19)

where $\bar{y}_{t}=\frac{1}{N}\sum_{\ell=1}^{N}y_{t}^{(\ell)}$.

###### Remark 4.3

Note that the proposed estimation scheme regarding the first $T$ Markov
parameters can be replaced by other Markov parameter estimators.

###  4.2 Identification Algorithm

We now propose our regularized identification algorithm in the noisy scenario,
which aims to minimize the mean squared identification error of the Markov
parameter $H_{k}$:

###### Problem 4.4

$\boldsymbol{\alpha}^{*}=\text{argmin}_{\\{\alpha_{t}\\}_{t=0}^{T-1}}\left\|\lambda^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}\lambda^{t}\right\|_{\infty,\rho}^{2}+\gamma\|\boldsymbol{\alpha}\|_{1}^{2},$
(20)

in which $\gamma$ is a parameter indicating the relative importance between
the above two terms, and shall be chosen as
$\gamma=\frac{\hat{\Sigma}}{C_{m}^{2}N}$ according to (17).

Suppose $N$ sample trajectories are leveraged to estimate the first $T$ Markov
parameters $\\{H_{t}\\}_{t=1}^{T}$, the inference of $H_{k}$ can be calculated
as:

$\hat{H}_{k}=\left\\{\begin{aligned} &\tilde{H}_{k}^{(N)},&k\leq T;\\\
&\sum_{t=1}^{T}\alpha_{t-1}^{*}\tilde{H}_{t}^{(N)},&k>T.\end{aligned}\right.$
(21)

Please refer to Algorithm 1 for a detailed illustration of the proposed
algorithm.

Algorithm 1 Complete procedure of our regularized identification algorithm
proposed in Section 4.2.

The index of Markov parameter to be identified $k\geq 1$, maximal steps of
each episode $T$, total iteration epochs $N$, modified system energy $C_{m}$

Identified $k$-th Markov parameter $\hat{H}_{k}$

Initialization $\tilde{H}_{1}\leftarrow 0,\cdots,\tilde{H}_{T}\leftarrow 0$,
$\ell\leftarrow 1$

while $\ell\leq N$ do$\triangleright$ Identification of the first $T$ Markov
parameters

$\tilde{H}_{t}\leftarrow\tilde{H}_{t}+\frac{1}{\ell}\left(y_{t}^{(\ell)}-\tilde{H}_{t}\right),t=1,\cdots,T$

$\ell\leftarrow\ell+1$

end while

if $k>T$ then$\triangleright$ Identification of $H_{k}$ via polynomial
approximation

$\bar{y}_{t}\leftarrow\frac{1}{N}\sum_{\ell=1}^{N}y_{t}^{(\ell)},t=1,\cdots,T$

$\hat{\Sigma}\leftarrow\frac{1}{TN-1}\sum_{t=1}^{T}\sum_{\ell=1}^{N}\left(y_{t}^{(\ell)}-\bar{y}_{t}\right)^{2}$

$\gamma\leftarrow\frac{\hat{\Sigma}}{C_{m}^{2}N}$

$\boldsymbol{\alpha}^{*}\leftarrow\text{argmin}_{\\{\alpha_{t}\\}_{t=0}^{T-1}}\left\|\lambda^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}\lambda^{t}\right\|_{\infty,\rho}^{2}+\gamma\|\boldsymbol{\alpha}\|_{1}^{2}$

$\hat{H}_{k}\leftarrow\sum_{t=1}^{T}\alpha_{t-1}^{*}\tilde{H}_{t}$

else

$\hat{H}_{k}\leftarrow\tilde{H}_{k}$

end if

##  5 Numerical Results

This section compares the finite-time performance of the unbiased Ho-Kalman
algorithm (Oymak and Ozay, 2019), the truncation method introduced in Remark
3.1 and the proposed algorithm using numerical experiments.

We consider the following $6$-dimensional SISO system with $\rho=0.95$:

$\displaystyle A$
$\displaystyle=\text{diag}(0.94,0.75,-0.75,-0.69,0.46,0.42),$ (22)
$\displaystyle B$
$\displaystyle=\begin{bmatrix}1&1&1&1&1&1\end{bmatrix}^{\top},C=\begin{bmatrix}1&1&1&1&1&1\end{bmatrix}.$

We iteratively estimate the target Markov parameter after each episode, and
let the weighting coefficient $\gamma$ in (20) attenuate at a rate of
$\mathcal{O}\left(\frac{1}{N}\right)$, where $N$ is the number of epochs, to
illustrate the finite-time performance of our algorithm with the growth of the
amount of data sampled. The identification result of $H_{13}$ and $H_{22}$
using $\\{\tilde{H}_{i}\\}_{i=1}^{12}$ by our algorithm, the Ho-Kalman
algorithm, and the truncation method are compared in Figure 1.

Figure 1: Identification error of $H_{13}$ and $H_{22}$ using
$\\{\tilde{H}_{i}\\}_{i=1}^{12}$ by the Ho-Kalman algorithm, the proposed
algorithm and the truncation method respectively, which is depicted in a log-
log plot. The yellow, blue, and gray lines depict the identification error of
the Ho-Kalman algorithm, the proposed algorithm, and the truncation method,
respectively. Furthermore, the orange line denotes the asymptotic
identification bias of the proposed algorithm when $N$ tends to infinity, and
the weighting parameter in (20) tends to $0$.

As can be observed from Figure 1, although the Ho-Kalman algorithm is proved
to be asymptotically unbiased, the finite-time identification result
oscillates to at least $10$ times larger than the true value of $H_{13}$, and
shows little decay during the first $1000$ episodes. On the contrary, our
algorithm smoothly identifies the Markov parameter and quickly converges to
the asymptotic value. Our algorithm also outperforms the naive truncation
method, which reveals the effectiveness of the proposed optimization problem.

Similar results are obtained from the identification process of $H_{22}$
plotted in Figure 1, where our algorithm significantly outperforms the Ho-
Kalman algorithm. In contrast, the bias of our algorithm remains negligible
compared to the impact of noise in the original data. Note that the asymptotic
bias can be further reduced by slightly increasing the number of known Markov
parameters $T$, where the impact of $T$ on the identification bias is
illustrated in Section 3.1.

Finally, the identification error of $\\{H_{i}\\}_{i=13}^{50}$ using the first
$12$ Markov parameters estimated from $1000$ sample trajectories of system
(22) in a single experiment are shown in Figure 2. The result reveals that the
identification error of the Ho-Kalman algorithm explodes exponentially w.r.t.
the index of Markov parameters due to the fact that some poles of the
identified $A$ matrix are unstable. On the other hand, the identification
result of our algorithm outperforms both the Ho-Kalman algorithm and the naive
truncation method, while the bias of our algorithm remains imperceptible.

Figure 2: Identification error of the Ho-Kalman algorithm, the proposed
algorithm and the truncation method w.r.t. the index of Markov parameters
using the first $12$ Markov parameters from $1000$ sample trajectories of
system (22) in a single experiment, which is depicted in a log plot. The
yellow, blue, and gray lines depict the identification error of the Ho-Kalman
algorithm, the proposed algorithm, and the truncation method, respectively.
The orange line denotes the asymptotic identification error of the proposed
algorithm when $N$ tends to infinity, and the weighting parameter in (20)
tends to $0$.

##  6 Conclusion

This paper proposes a Markov parameter identification algorithm for SISO LTI
systems from the first $T$ Markov parameters, which is suitable for both the
noise-free and the noisy scenarios. We further provide the upper bound on the
asymptotic identification bias when the first $T$ Markov parameters are
precisely estimated. In addition, the mean squared identification error in the
noisy scenario is derived, by which we formulate the finite-time
identification problem as a regularized uniform polynomial approximation
problem. The identification algorithm for systems with complex roots is left
for future research.

## References

  * Chiuso and Picci (2004) Chiuso, A. and Picci, G. (2004).  On the ill-conditioning of subspace identification with inputs.  _Automatica_ , 40(4), 575–589.  10.1016/j.automatica.2003.11.009.  URL https://linkinghub.elsevier.com/retrieve/pii/
S0005109803003674.

  * Faradonbeh et al. (2017) Faradonbeh, M.K.S., Tewari, A., and Michailidis, G. (2017).  Finite time analysis of optimal adaptive policies for linear-quadratic systems.  _arXiv preprint arXiv:1711.07230_.
  * Goldenshluger (1998) Goldenshluger, A. (1998).  Nonparametric estimation of transfer functions: rates of convergence and adaptation.  _IEEE Transactions on Information Theory_ , 44(2), 644–658.  10.1109/18.661510.  Conference Name: IEEE Transactions on Information Theory.
  * Hachicha et al. (2014) Hachicha, S., Kharrat, M., and Chaari, A. (2014).  N4sid and moesp algorithms to highlight the ill-conditioning into subspace identification.  _International Journal of Automation and Computing_ , 11(1), 30–38.
  * Iannelli et al. (2021) Iannelli, A., Yin, M., and Smith, R.S. (2021).  Experiment design for impulse response identification with signal matrix models.  URL http://arxiv.org/abs/2012.08126.  ArXiv:2012.08126 [cs, eess].
  * Ikeda and Oku (2015) Ikeda, K. and Oku, H. (2015).  Estimation error analysis of system matrices in some subspace identification methods.  In _2015 10th Asian Control Conference (ASCC)_ , 1–6.  10.1109/ASCC.2015.7244557.
  * Kritayakirana and Gerdes (2012) Kritayakirana, K. and Gerdes, J.C. (2012).  Using the centre of percussion to design a steering controller for an autonomous race car.  _Vehicle System Dynamics_ , 50(sup1), 33–51.  10.1080/00423114.2012.672842.  URL https://doi.org/10.1080/00423114.2012.672842.
  * Li et al. (2022) Li, J., Sun, S., and Mo, Y. (2022).  Fundamental limit on siso system identification.  In _2022 IEEE 61st Conference on Decision and Control (CDC)_ , 856–861.  10.1109/CDC51059.2022.9993203.
  * Ljung and Wahlberg (1992) Ljung, L. and Wahlberg, B. (1992).  Asymptotic Properties of the Least-Squares Method for Estimating Transfer Functions and Disturbance Spectra.  _Advances in Applied Probability_ , 24(2), 412–440.  10.2307/1427698.  URL https://www.jstor.org/stable/1427698.  Publisher: Applied Probability Trust.
  * Ljung and Yuan (1985) Ljung, L. and Yuan, Z.D. (1985).  Asymptotic properties of black-box identification of transfer functions.  _IEEE Transactions on Automatic Control_ , 30(6), 514–530.  10.1109/TAC.1985.1103995.  Conference Name: IEEE Transactions on Automatic Control.
  * Oymak and Ozay (2019) Oymak, S. and Ozay, N. (2019).  Non-asymptotic identification of lti systems from a single trajectory.  In _2019 American control conference (ACC)_ , 5655–5661. IEEE.
  * Rantzer (2018) Rantzer, A. (2018).  Concentration Bounds for Single Parameter Adaptive Control.  In _2018 Annual American Control Conference (ACC)_ , 1862–1866.  10.23919/ACC.2018.8431891.  ISSN: 2378-5861.
  * Saibaba (2021) Saibaba, A.K. (2021).  Approximating monomials using chebyshev polynomials.  _arXiv preprint arXiv:2101.06818_.
  * Sarkar and Rakhlin (2019) Sarkar, T. and Rakhlin, A. (2019).  Near optimal finite time identification of arbitrary linear dynamical systems.  _arXiv:1812.01251 [cs]_.  URL http://arxiv.org/abs/1812.01251.  ArXiv: 1812.01251.
  * Schroeck et al. (2001) Schroeck, S.J., Messner, W.C., and McNab, R.J. (2001).  On compensator design for linear time-invariant dual-input single-output systems.  _IEEE/ASME Transactions on mechatronics_ , 6(1), 50–57.
  * Shirani Faradonbeh et al. (2018) Shirani Faradonbeh, M.K., Tewari, A., and Michailidis, G. (2018).  Finite time identification in unstable linear systems.  _Automatica_ , 96, 342–353.  10.1016/j.automatica.2018.07.008.  URL https://linkinghub.elsevier.com/retrieve/pii/
S0005109818303546.

  * Simchowitz et al. (2019) Simchowitz, M., Boczar, R., and Recht, B. (2019).  Learning Linear Dynamical Systems with Semi-Parametric Least Squares.  In _Proceedings of the Thirty-Second Conference on Learning Theory_ , 2714–2802. PMLR.  URL https://proceedings.mlr.press/v99/simchowitz
19a.html.  ISSN: 2640-3498.

  * Simchowitz et al. (2018) Simchowitz, M., Mania, H., Tu, S., Jordan, M.I., and Recht, B. (2018).  Learning Without Mixing: Towards A Sharp Analysis of Linear System Identification.  _arXiv:1802.08334 [cs, math, stat]_.  URL http://arxiv.org/abs/1802.08334.  ArXiv: 1802.08334.
  * Sun and Mo (2022) Sun, S. and Mo, Y. (2022).  Fundamental identification limit on diagonal canonical form for siso systems.  In _2022 IEEE 17th International Conference on Control & Automation (ICCA)_, 728–733.  10.1109/ICCA54724.2022.9831907.
  * Sun et al. (2022) Sun, S., Mo, Y., and You, K. (2022).  Fundamental identification limit of single-input and single-output linear time-invariant systems.  In _2022 13th Asian Control Conference (ASCC)_ , 2157–2162.  10.23919/ASCC56756.2022.9828137.
  * Tsiamis and Pappas (2019) Tsiamis, A. and Pappas, G.J. (2019).  Finite Sample Analysis of Stochastic System Identification.  In _2019 IEEE 58th Conference on Decision and Control (CDC)_ , 3648–3654. IEEE, Nice, France.  10.1109/CDC40024.2019.9029499.  URL https://ieeexplore.ieee.org/document/9029499/.
  * Tu et al. (2017) Tu, S., Boczar, R., Packard, A., and Recht, B. (2017).  Non-Asymptotic Analysis of Robust Control from Coarse-Grained Identification.  _arXiv:1707.04791 [cs, math]_.  URL http://arxiv.org/abs/1707.04791.  ArXiv: 1707.04791.
  * Wagenmaker and Jamieson (2020) Wagenmaker, A. and Jamieson, K. (2020).  Active Learning for Identification of Linear Dynamical Systems.  In _Proceedings of Thirty Third Conference on Learning Theory_ , 3487–3582. PMLR.  URL https://proceedings.mlr.press/v125/wagenmaker
20a.html.  ISSN: 2640-3498.

  * Yin et al. (2021) Yin, M., Iannelli, A., and Smith, R.S. (2021).  Maximum Likelihood Estimation in Data-Driven Modeling and Control.  _IEEE Transactions on Automatic Control_ , 1–1.  10.1109/TAC.2021.3137788.  Conference Name: IEEE Transactions on Automatic Control.
  * Zheng and Li (2021) Zheng, Y. and Li, N. (2021).  Non-Asymptotic Identification of Linear Dynamical Systems Using Multiple Trajectories.  _IEEE Control Systems Letters_ , 5(5), 1693–1698.  10.1109/LCSYS.2020.3042924.  URL https://ieeexplore.ieee.org/document/9284539/.
  * Şen and Kalyoncu (2019) Şen, M.A. and Kalyoncu, M. (2019).  Grey wolf optimizer based tuning of a hybrid lqr-pid controller for foot trajectory control of a quadruped robot.  _Gazi University Journal of Science_ , 32(2), 674 – 684.

##  Appendix A Proof of theorem 1

Using the specific expression of $H_{k}$ in (5), the bias of the estimation
$\hat{H}_{k}$ has the following relationship with the polynomial approximation
error:

$\displaystyle\left|H_{k}-\sum_{t=1}^{T}\alpha_{t-1}^{*}H_{t}\right|$
$\displaystyle=\left|\sum_{j=1}^{n}c_{j}\lambda_{j}^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}^{*}\sum_{j=1}^{n}c_{j}\lambda_{j}^{t}\right|$
(23)
$\displaystyle=\left|\sum_{j=1}^{n}c_{j}\left(\lambda_{j}^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}^{*}\lambda_{j}^{t}\right)\right|$
$\displaystyle\leq
C_{m}\left\|x^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}^{*}x^{t}\right\|_{\infty,\rho},$

where $\boldsymbol{\alpha}^{*}$ is obtained by (6).

  * •

Proof of $|H_{k}-\hat{H}_{k}|^{2}\leq
4C_{m}^{2}\rho^{2k-2}\exp\left(-\frac{2(T-1)^{2}}{k-1}\right)$

According to Saibaba (2021), the Chebyshev approximation problem has the
following upper bound on the approximation error, where $\mu\in[-1,1]$:

$\left\|\mu^{k-1}-\sum_{t=0}^{T-1}\beta_{t}^{*}\mu^{t}\right\|_{\infty,1}\leq
2\exp\left(-\frac{(T-1)^{2}}{k-1}\right),$ (24)

and the optimal solution $\beta_{t}^{*}$ is unique. Since $\mu\in[-1,1]$, take
$x=\rho\mu$, then $x\in[-\rho,\rho]$, and the following inequality holds by
putting $\mu=\rho^{-1}x$ into (24):

$\displaystyle\left\|\rho^{1-k}x^{k-1}-\sum_{t=0}^{T-1}\rho^{-t}\beta^{*}_{t}x^{t}\right\|_{\infty,\rho}\leq
2\exp\left(-\frac{(T-1)^{2}}{k-1}\right).$ (25)

Therefore,

$\displaystyle\left\|x^{k-1}-\sum_{t=0}^{T-1}\rho^{k-t-1}\beta^{*}_{t}x^{t}\right\|_{\infty,\rho}\leq
2\rho^{k-1}\exp\left(-\frac{(T-1)^{2}}{k-1}\right).$ (26)

Thus, the upper bound can be proved by combining (26) with (23).

  * •

Proof of $|H_{k}-\hat{H}_{k}|^{2}\leq\frac{1}{4}C_{m}^{2}\rho^{2k-2}$

The bound can be derived by considering the following two cases. When $k$ is
odd, it can be verified that the following inequality holds:

$\left|x^{k-1}-\frac{1}{2}\rho^{k-1}\right|\leq\frac{1}{2}\rho^{k-1},\forall
x\in[0,\rho].$ (27)

Since LHS is an even function, the polynomial approximation error can be upper
bounded by:

$\left\|x^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}^{*}x^{t}\right\|_{\infty,\rho}\leq\left\|x^{k-1}-\frac{1}{2}\rho^{k-1}\right\|_{\infty,\rho}\leq\frac{1}{2}\rho^{k-1}.$
(28)

Similarly, when $k$ is even, the following inequality holds:

$\left|x^{k-1}-\frac{1}{2}\rho^{k-2}x\right|\leq\frac{1}{2}\rho^{k-1},\forall
x\in[0,\rho].$ (29)

Thus,

$\left\|x^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}^{*}x^{t}\right\|_{\infty,\rho}\leq\left\|x^{k-1}-\frac{1}{2}\rho^{k-2}x\right\|_{\infty,\rho}\leq\frac{1}{2}\rho^{k-1}.$
(30)

Therefore, the second upper bound in (8) can be derived by combining (23),
(28) and (30).

##  Appendix B Proof of Theorem 2

The mean squared approximation error could be expanded as:

$\displaystyle\mathbb{E}\left(H_{k}-\hat{H}_{k}\right)^{2}=\mathbb{E}\left(H_{k}-\sum_{t=1}^{T}\alpha_{t-1}\tilde{H}_{t}^{(N)}\right)^{2}$
(31)
$\displaystyle=\mathbb{E}\left(H_{k}-\sum_{t=1}^{T}\alpha_{t-1}H_{t}+\sum_{t=1}^{T}\alpha_{t-1}H_{t}-\sum_{t=1}^{T}\alpha_{t-1}\tilde{H}_{t}^{(N)}\right)^{2}$
$\displaystyle=\left(H_{k}-\sum_{t=1}^{T}\alpha_{t-1}H_{t}\right)^{2}+\mathbb{E}\left[\sum_{t=1}^{T}\alpha_{t-1}\left(H_{t}-\tilde{H}_{t}^{(N)}\right)\right]^{2},$

which is the bias-variance decomposition.

Therein, the bias can be written as follows according to (23):

$\left(H_{k}-\sum_{t=1}^{T}\alpha_{t-1}H_{t}\right)^{2}\leq
C_{m}^{2}\left\|x^{k-1}-\sum_{t=0}^{T-1}\alpha_{t}x^{t}\right\|_{\infty,\rho}^{2}.$
(32)

On the other hand, for random variables $X$ and $Y$, the following
inequalities hold:

$\text{Var}(X+Y)=\text{Var}(X)+\text{Var}(Y)+2\text{Cov}(X,Y)$
$\displaystyle\leq\text{Var}(X)+\text{Var}(Y)+2\sqrt{\text{Var}(X)\text{Var}(Y)}$
(33)
$\displaystyle\leq\left(\sqrt{\text{Var}(X)}+\sqrt{\text{Var}(Y)}\right)^{2}.$

Thus, the variance term satisfies:

$\displaystyle\mathbb{E}\left[\sum_{t=1}^{T}\alpha_{t-1}\left(H_{t}-\tilde{H}_{t}^{(N)}\right)\right]^{2}=\text{Var}\left[\sum_{t=1}^{T}\alpha_{t-1}\tilde{H}_{t}^{({N})}\right]$
(34)
$\displaystyle\leq\left(\sum_{t=1}^{T}\sqrt{\text{Var}\left(\alpha_{t-1}\tilde{H}_{t}^{(N)}\right)}\right)^{2}.$

By expanding $\tilde{H}_{t}^{(N)}$ as the mean of the $t$-th system output
from $N$ independent experiments, we can obtain the following result:

$\displaystyle\text{Var}\left(\tilde{H}_{t}^{(N)}\right)$
$\displaystyle=\text{Var}\left(\frac{1}{N}\sum_{\ell=1}^{N}y_{t}^{(\ell)}\right)$
(35)
$\displaystyle=\frac{1}{N^{2}}\sum_{\ell=1}^{N}\text{Var}\left(y_{t}^{(\ell)}\right)=\frac{\Sigma}{N},$

where $\Sigma$ denotes the variance of system output $y_{t}$ in each episode,
which is the same for each $t$ due to the fact that the input in the
experiment is determined. Furthermore, since each episode is reset to steady
state, $\Sigma$ is equal to:

$\Sigma=CPC^{\top}+R,\text{where}\ P=APA^{\top}+Q.$ (36)

Therefore, equation (34) can be further derived as:

$\left(\sum_{t=1}^{T}\sqrt{\alpha_{t-1}^{2}\text{Var}\left(\tilde{H}_{t}^{(N)}\right)}\right)^{2}\leq\frac{\Sigma}{N}\left(\sum_{t=0}^{T-1}|\alpha_{t}|\right)^{2}=\frac{\Sigma}{N}\|\boldsymbol{\alpha}\|_{1}^{2}.$
(37)

Thus, the theorem can be proved by combining (32) with the result above.