# Identification of Linear Systems with Multiplicative Noise from Multiple
Trajectory Data

Yu Xing  Benjamin Gravell  Xingkang He  Karl Henrik Johansson  Tyler Summers
Division of Decision and Control Systems, School of Electrical Engineering and
Computer Science,
KTH Royal Institute of Technology, and Digital Futures, Stockholm, Sweden.
Department of Mechanical Engineering, The University of Texas at Dallas,
Richardson, TX, USA.  Department of Electrical Engineering, University of
Notre Dame, South Bend, IN, USA.

###### Abstract

The paper studies identification of linear systems with multiplicative noise
from multiple-trajectory data. An algorithm based on the least-squares method
and multiple-trajectory data is proposed for joint estimation of the nominal
system matrices and the covariance matrix of the multiplicative noise. The
algorithm does not need prior knowledge of the noise or stability of the
system, but requires only independent inputs with pre-designed first and
second moments and relatively small trajectory length. The study of
identifiability of the noise covariance matrix shows that there exists an
equivalent class of matrices that generate the same second-moment dynamic of
system states. It is demonstrated how to obtain the equivalent class based on
estimates of the noise covariance. Asymptotic consistency of the algorithm is
verified under sufficiently exciting inputs and system controllability
conditions. Non-asymptotic performance of the algorithm is also analyzed under
the assumption that the system is bounded. The analysis provides high-
probability bounds vanishing as the number of trajectories grows to infinity.
The results are illustrated by numerical simulations.

###### keywords:

linear system identification, multiplicative noise, multiple trajectories,
non-asymptotic results

AND OR NOT TO COMMENT BODY IF ELSE ELSIF FOR WHILE

††thanks: Y. Xing and B. Gravell contributed equally to this paper.
_Email addresses:_<EMAIL_ADDRESS>(Y. Xing<EMAIL_ADDRESS>(B.
Gravell<EMAIL_ADDRESS>(X. He<EMAIL_ADDRESS>(K. H. Johansson),
<EMAIL_ADDRESS>(T. Summers).

, , , ,

##  1 Introduction

The study of stochastic systems with multiplicative noise (i.e., system states
and inputs multiplied by noise) has a long history in control theory (Wonham,,
1967), and is re-emerging in the context of complex networked systems and
learning-based control. In contrast to the additive-noise setting, the
multiplicative-noise modeling framework has the ability to capture the
coupling between noise and system states. This situation occurs in modern
control systems as diverse as robotics with distance-dependent sensor errors
(Du Toit and Burdick,, 2011), networked systems with noisy communication
channels (Antsaklis and Baillieul,, 2007; Hespanha et al.,, 2007), modern
power networks with high penetration of intermittent renewables (Guo and
Summers,, 2019), turbulent fluid flow (Lumley,, 2007), and neuronal brain
networks (Breakspear,, 2017). Linear systems with multiplicative noise are
particularly attractive as a stochastic modeling framework because they remain
simple enough to admit closed-form expressions for stabilization (Boyd et
al.,, 1994) and optimal control (Wonham,, 1967; Kleinman,, 1969; Gravell et
al.,, 2021).

It is important to study identification of linear systems with multiplicative
noise, because, when solving problems such as control design of
multiplicative-noise linear quadratic regulator (LQR), system parameters
including the nominal system matrices and the noise covariance matrix,
especially the latter, generally need be known (Gravell et al.,, 2021). In
contrast, for the design problem of additive-noise LQR, the covariance matrix
of additive noise needs not be known (Dean et al.,, 2019). Moreover, the
identification problem requires further investigation; for instance, it is
unclear how to formally quantify identifiability issues resulting from
coupling between system states and multiplicative noise, and how to design
identification algorithms to efficiently tackle the influence of
multiplicative noise.

Another issue that must be addressed is how to perform system identification
based on multiple-trajectory data, rather than on single-trajectory data.
Multiple-trajectory data arises in two broad situations: (1) episodic tasks
where a system is reset to an initial state after a finite run time, as
encountered in iterative learning control and reinforcement learning (Matni et
al.,, 2019); and (2) data collected from multiple identical systems in
parallel, for example, robotic-grasping dataset collected by Google running
several robot arms concurrently (Gu et al.,, 2017; Levine et al.,, 2018). For
multiple-trajectory data, the length of each trajectory may be small, but the
number of trajectories can be large. However, the classic literature of system
identification mainly focuses on studying online estimation over a single
trajectory, so there is a need to study how to identify systems based on
multiple-trajectory data. In addition, system identification based on multiple
trajectories can be a pre-step of conducting other tasks such as control
design of LQR (Dean et al.,, 2019). Thus, studying the performance of
identification algorithms based on multiple trajectories is necessary for
obtaining performance guarantees of later tasks.

###  1.1 Related Work

For identification of a nominal linear system, recursive algorithms, such as
the recursive least-squares algorithm, have been developed in the control
literature (Lai and Wei,, 1982; Ljung,, 1986; Chen and Guo,, 2012). These
algorithms can be applied to identification of linear systems with
multiplicative noise, provided that certain conditions of system stability and
noise hold. Non-asymptotic performance analysis of identification methods can
be found in Weyer and Campi, (2002); Campi and Weyer, (2002, 2005). It has
once again attracted attention from different domains and been investigated
more extensively, because of recent development of random matrix theories,
self-normalized martingales, and so on (see Dean et al., (2019); Matni and Tu,
(2019); Zheng and Li, (2020) and references therein).

For estimation of noise covariance, both recursive and batch methods have been
proposed over the last few decades (Duník et al.,, 2017), but most of these
methods focus on the additive-noise case. In order to estimate multiplicative
noise covariance, Schön et al., (2011) introduces a maximum-likelihood
approach, and Kitagawa, (1998); Kantas et al., (2015) utilize Bayesian
frameworks. These methods, however, require prior assumptions on the noise
distributions, whose incorrectness may worsen algorithm performance. Coppens
and Patrinos, (2020); Coppens et al., (2020) study stochastic LQR design for a
special case of linear systems with multiplicative noise. It is assumed that
the multiplicative noise is observed directly so that a concentration
inequality can be obtained for estimates of the noise covariance. The most
relevant work to our paper is Di and Lamperski, (2021), which studies
simultaneously estimating the nominal system parameters and noise covariance
matrix based on single-trajectory data. In that paper, a self-normalizing
(ellipsoidal) bound and a Euclidean (box) bound are provided for least-squares
estimates, but it is not clear whether the bounds converge to zero under the
setting of linear systems with multiplicative noise.

There is a growing interest in system identification based on multiple-
trajectory data, along with their applications in data-driven control (Dean et
al.,, 2019; Matni and Tu,, 2019), due to the powerful and convenient estimator
schemes facilitated by resetting the system. This framework can be applied to
both stable and unstable systems, because of the finite duration of each
trajectory. The authors in Tu and Recht, (2018); Sun et al., (2020) introduce
the procedure of collecting multiple trajectories, to identify finite impulse
response systems. In Dean et al., (2019), the authors develop a framework
called coarse-ID control to solve the problem of LQR with unknown linear
dynamics. The first step of this framework is to learn a coarse model of the
unknown linear system, by observing multiple independent trajectories with
finite length. However, only the last input-state pairs of the trajectories
are used in the theoretical analysis of the learning algorithm. The
performance of a least-squares algorithm, using all samples of every
trajectory, is studied in Zheng and Li, (2020), for partially observed,
possibly open-loop unstable, linear systems.

###  1.2 Contributions

This paper considers identification of linear systems with multiplicative
noise from multiple-trajectory data. The contributions are three-fold:

  * 1.

An algorithm (Algorithm 1) based on the least-squares method and multiple-
trajectory data is proposed for joint identification of the nominal system
matrices and the multiplicative noise covariance from multiple-trajectory
data. The algorithm does not need prior knowledge of the noise or stability of
the system, but requires only independent inputs with pre-designed first and
second moments, relatively small length for each trajectory, and the
assumption of independent and identically distributed (i.i.d.) noise with
finite first and second moments. It is theoretically shown that, under the
preceding conditions, the algorithm solves the identification problem.

  * 2.

Identifiability of the noise covariance matrix is investigated (Propositions 1
and 2). It is shown that there exists an equivalent class of covariance
matrices that generate the same second-moment dynamic of system states. In
addition, it is studied when such equivalent class has a unique element,
meaning that the covariance matrix can be uniquely determined. An explicit
expression of the equivalent class is provided for the recovery of the noise
covariance based on estimates given by the proposed algorithm.

  * 3.

Asymptotic consistency of the proposed algorithm is verified (Theorem 3.10),
under sufficiently exciting inputs and system controllability conditions. Non-
asymptotic estimation performance is also analyzed under the assumption that
the system is bounded. This analysis provides high-probability error bounds,
which vanish as the number of trajectories grows to infinity (Theorems 3.16
and 3.17).

Compared with Di and Lamperski, (2021), the current paper provides high-
probability error bounds, for the proposed algorithm, that converge to zero as
the number of trajectories increases. In addition, identifiability of the
noise covariance matrix is thoroughly studied, and conditions, under which the
covariance matrix is uniquely determined, are provided. In our problem,
because of the complicated structure of the second-moment dynamic of system
states, both analysis of the error bounds and study of the identifiability
require more elaborate use of tools from linear algebra and high-dimensional
probability theory. The differences between this paper and its conference
version (Xing et al.,, 2020) are as follows. This paper studies
identifiability of the noise covariance matrix in detail, demonstrating a
framework to recover the equivalent class of covariance matrices. Moreover,
sharper bounds for the required length of each trajectory are obtained.
Finally, finite sample analysis of the proposed algorithm is provided.

###  1.3 Outline

The remainder of the paper is organized as follows. The problem is formulated
in Section 2. In Section 3 the algorithm is introduced and theoretical results
are given. Numerical simulation results are presented in Section 4. Section 5
concludes the paper. Some proofs are postponed to Appendix.

Notation.
Denote the $n$-dimensional Euclidean space by $\mathds{R}^{n}$, and the set of
$n\times m$ real matrices by $\mathds{R}^{n\times m}$. Let $\mathds{N}$ stand
for the set of nonnegative integers, and
$\mathds{N}^{+}:=\mathds{N}\setminus\\{0\\}$. Let $[k]:=\\{1,2,\dots,k\\}$,
$k\in\mathds{N}^{+}$. We use $\|\cdot\|$ to denote the Euclidean norm for
vectors, and use $\|\cdot\|_{F}$ and $\|\cdot\|_{2}$ to denote the Frobenius
and spectral norm for matrices. The probability of an event $E$ is denoted by
$\mathds{P}\\{E\\}$, and the expectation of a random vector $x$ is represented
by $\mathds{E}\\{x\\}$. An event happening almost surely (a.s.) means that it
happens with probability one. Let $A\times B$ be the Cartesian product of sets
$A$ and $B$, namely, $A\times B=\\{(a,b):a\in A,b\in B\\}$. For two sequences
of real numbers $a_{k}$ and $b_{k}\not=0$, $k\in\mathds{N}^{+}$, denote
$a_{k}=\mathcal{O}(b_{k})$, if there exists a positive constant $C$ such that
$|a_{k}/b_{k}|\leq C$ for all $k\in\mathds{N}^{+}$.

Let $a_{ij}$ or $[A]_{ij}$ represent the $(i,j)$-th entry of
$A\in\mathds{R}^{n\times m}$. Denote the $n$-dimensional all-one vector and
all-zero vector by $\mathbf{1}_{n}$ and $\mathbf{0}_{n}$, respectively. The
$n$-dimensional unit vector with $i$-th component being one is represented by
$\mathbf{e}_{i}^{n}$. $I_{n}$ is the $n$-dimensional identity matrix. For two
symmetric matrices $A,B\in\mathds{R}^{n\times n}$, $A\succeq 0$ ($A\succ 0$)
means that $A$ is positive semidefinite (positive definite), and $A\succeq B$
($A\succ B$) means that $A-B\succeq 0$ ($A-B\succ 0$). For a matrix
$A\in\mathds{R}^{n\times n}$, $\rho(A)$ represents the spectral radius of $A$.
For a symmetric matrix $A\in\mathds{R}^{n}$, denote its smallest and largest
eigenvalue by $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ respectively. A
block diagonal matrix $A$ with $A_{1}$, $\dots$, $A_{k}$ on its diagonal is
denoted by $\text{blockdiag}(A_{1},\dots,A_{k})$.

The Kronecker product of two matrices $A\in\mathds{R}^{m\times n}$ and
$B\in\mathds{R}^{p\times q}$ is represented by $A\otimes B$. The full
vectorization of $A=[a_{ij}]\in\mathds{R}^{m\times n}$ is found by stacking
the columns of $A$ (i.e.,
${\operatorname{vec}(A)=[a_{11}~{}a_{21}~{}\cdots~{}a_{m1}~{}a_{12}~{}a_{22}~{}\cdots~{}a_{mn}]^{\intercal}}$).
The symmetric vectorization (also called half-vectorization) of a symmetric
matrix $A\in\mathds{R}^{n\times n}$ is found by stacking the upper triangular
part of the columns of $A$ (i.e.,
${\operatorname{svec}(A)=[a_{11}~{}a_{12}~{}a_{22}~{}\cdots~{}a_{1n}~{}a_{2n}~{}\cdots~{}a_{nn}]^{\intercal}}$).
The inverse operations of $\operatorname{vec}(\cdot)$ and
$\operatorname{svec}(\cdot)$, given $p,q\in\mathds{N}$, are the full
matricization $\operatorname{mat}_{p\times
q}(x):=(\operatorname{vec}(I_{q})^{\intercal}\otimes I_{p})(I_{q}\otimes x)$
for a vector $x\in\mathds{R}^{pq}$ and symmetric matricization
$\operatorname{smat}_{p}(y)$ for a vector $y\in\mathds{R}^{p(p+1)/2}$,
respectively. To generalize the vectorization and matricization operations to
a block matrix

$\displaystyle B=\begin{bmatrix}B_{11}&B_{12}&\cdots&B_{1n}\\\
\vdots&\vdots&&\vdots\\\
B_{m1}&B_{m2}&\cdots&B_{mn}\end{bmatrix}\in\mathds{R}^{mp\times nq},$

where $B_{ij}\in\mathds{R}^{p\times q}$, define the following matrix reshaping
operator $F:\mathds{R}^{mp\times nq}\rightarrow\mathds{R}^{mn\times pq}$,

$\displaystyle
F(B,m,n,p,q):=[\operatorname{vec}(B_{11})~{}\operatorname{vec}(B_{21})~{}\cdots~{}\operatorname{vec}(B_{m1})\cdots~{}\operatorname{vec}(B_{12})~{}\operatorname{vec}(B_{22})~{}\cdots~{}\operatorname{vec}(B_{mn})]^{\intercal}.$

Then it holds that $F(A\otimes
A,m,n,m,n)=\operatorname{vec}(A)\operatorname{vec}(A)^{\intercal}$ for
$A\in\mathds{R}^{m\times n}$, which demonstrates the correspondence between
the entries of $A\otimes A$ and those of
$\operatorname{vec}(A)\operatorname{vec}(A)^{\intercal}$. Note when $p=q=1$,
$F(\cdot)$ degenerates to $\operatorname{vec}(\cdot)$. Define the inverse
reshaping operator $G:\mathds{R}^{mn\times pq}\to\mathds{R}^{mp\times nq}$ as

$\displaystyle G(B,m,n,p,q):=\begin{bmatrix}\operatorname{mat}_{p\times
q}(B_{1})&\cdots&\operatorname{mat}_{p\times q}(B_{(n-1)m+1})\\\
\operatorname{mat}_{p\times q}(B_{2})&\cdots&\operatorname{mat}_{p\times
q}(B_{(n-1)m+2})\\\ \vdots&&\vdots\\\ \operatorname{mat}_{p\times
q}(B_{m})&\cdots&\operatorname{mat}_{p\times q}(B_{mn})\end{bmatrix},$

where $B\in\mathds{R}^{mn\times pq}$, $B_{i}^{\intercal}$ is the $i$-th row of
$B$. Thus $F$ and $G$ are inverses of each other in the sense that

$\displaystyle F(G(A,m,n,p,q),m,n,p,q)$ $\displaystyle=A,$ $\displaystyle
G(F(B,m,n,p,q),m,n,p,q)$ $\displaystyle=B,$

for any $A\in\mathds{R}^{mn\times pq}$ and $B\in\mathds{R}^{mp\times nq}$. In
this way,
$G(\operatorname{vec}(A)\operatorname{vec}(A)^{\intercal},m,n,m,n)=A\otimes A$
for $A\in\mathds{R}^{m\times n}$. Note that both $F$ and $G$ are linear:
$F(A+B,m,n,p,q)=F(A,m,n,p,q)+F(B,m,n,p,q)$ for $A,B\in\mathds{R}^{mp\times
nq}$, and $G(A+B,m,n,p,q)=G(A,m,n,p,q)+G(B,m,n,p,q)$ for
$A,B\in\mathds{R}^{mn\times pq}$.

##  2 Problem Formulation

Consider the linear system with multiplicative noise

$\displaystyle x_{t+1}$
$\displaystyle=(A+\bar{A}_{t})x_{t}+(B+\bar{B}_{t})u_{t},~{}t\in\mathds{N},$
(1)

where $x_{t}\in\mathds{R}^{n}$ is the system state, and
$u_{t}\in\mathds{R}^{m}$ is the control input, $m\leq n$. The system is
described by the nominal dynamic matrix $A\in\mathds{R}^{n\times n}$ and the
nominal input matrix $B\in\mathds{R}^{n\times m}$, and incorporates
multiplicative noise terms modeled by i.i.d. and mutually independent random
matrices $\bar{A}_{t}$ and $\bar{B}_{t}$, which have zero mean and covariance
matrices
$\Sigma_{A}:=\mathds{E}\\{\operatorname{vec}(\bar{A}_{t})\operatorname{vec}(\bar{A}_{t})^{T}\\}\in\mathds{R}^{n^{2}\times
n^{2}}$ and
$\Sigma_{B}:=\mathds{E}\\{\operatorname{vec}(\bar{B}_{t})\operatorname{vec}(\bar{B}_{t})^{T}\\}\in\mathds{R}^{nm\times
nm}$, respectively. The multiplicative noise is assumed to be independent of
the inputs. Note that if $\bar{A}_{t}$ and $\bar{B}_{t}$ have non-zero means
$\bar{A}$ and $\bar{B}$, respectively, then we can consider a system with
nominal matrix $[A+\bar{A}~{}B+\bar{B}]$, as well as noise terms
$\bar{A}_{t}-\bar{A}$ and $\bar{B}_{t}-\bar{B}$, which satisfies the preceding
zero-mean assumption. The term multiplicative noise refers to that noise,
$\bar{A}_{t}$ and $\bar{B}_{t}$, enters the system as multipliers of $x_{t}$
and $u_{t}$, rather than as additions. The independence of $\bar{A}_{t}$ and
$\bar{B}_{t}$ is assumed for simplicity, and under this assumption the
covariance matrix of the entire multiplicative noise is a block diagonal
matrix
$\mathds{E}\\{\operatorname{vec}([\bar{A}_{t}~{}\bar{B}_{t}])\operatorname{vec}([\bar{A}_{t}~{}\bar{B}_{t}])^{\intercal}\\}=\text{blockdiag}(\Sigma_{A},\Sigma_{B})$.
Throughout the paper, we use
$(\Sigma_{A},\Sigma_{B})\in\mathds{R}^{n^{2}\times
n^{2}}\times\mathds{R}^{nm\times nm}$ to represent this matrix. If
$\bar{A}_{t}$ and $\bar{B}_{t}$ are dependent, there is an extra but amenable
term on their correlations,
$\mathds{E}\\{\operatorname{vec}(\bar{A}_{t})\operatorname{vec}(\bar{B}_{t})^{\intercal}\\}$.

An example of System (1) is the following system studied in the optimal
control literature (Boyd et al.,, 1994; Gravell et al.,, 2021),

$\displaystyle x_{t+1}$
$\displaystyle=\Big{(}A+\sum_{i=1}^{r}A_{i}p_{i,t}\Big{)}x_{t}+\Big{(}B+\sum_{j=1}^{s}B_{j}q_{j,t}\Big{)}u_{t},$
(2)

where $\\{p_{i,t}\\}$ and $\\{q_{i,t}\\}$ are mutually independent scalar
random variables, with $\mathds{E}\\{p_{i,t}\\}=\mathds{E}\\{q_{j,t}\\}=0$,
$\mathds{E}\\{p_{i,t}^{2}\\}=\sigma^{2}_{i}$, and
$\mathds{E}\\{q_{j,t}^{2}\\}=\delta^{2}_{j}$, $\forall
i\in[r],j\in[s],t\in\mathds{N}$. It can be seen that
$\bar{A}_{t}=\sum_{i=1}^{r}A_{i}p_{i,t}$ and
$\bar{B}_{t}=\sum_{j=1}^{s}B_{j}q_{j,t}$, where $\sigma_{i}$ and $\delta_{j}$
are the eigenvalues of $\Sigma_{A}$ and $\Sigma_{B}$, and $A_{i}$ and $B_{j}$
are the reshaped eigenvectors of $\Sigma_{A}$ and $\Sigma_{B}$. These
parameters are necessary for optimal controller design (Gravell et al.,,
2021). It is also possible to use System (2) to model cyber-physical systems
in which fault signals appear as multiplicative noise (Wang et al.,, 2020).
For new systems with unknown parameters, the key problem is to identify the
parameters in the first place. Another example of System (1) is interconnected
systems, where the nominal part captures relationships between different
subsystems, and multiplicative noise characterizes randomly varying topologies
(Haber and Verhaegen,, 2014).

In the rest of the paper, a trajectory sample is referred to as a _rollout_.
Suppose that multiple rollouts consisting of system states and inputs (i.e.,
$\\{[x_{0}^{(k)},u_{0}^{(k)},\dots,x_{\ell-1}^{(k)},u_{\ell-1}^{(k)},x_{\ell}^{(k)}],k\in[n_{r}]\\}$)
are available, where
$[x_{0}^{(k)},u_{0}^{(k)},\dots,x_{\ell-1}^{(k)},u_{\ell-1}^{(k)},x_{\ell}^{(k)}]$
is the $k$-th trajectory, $\ell$ is the length (index of the final time-step)
of every rollout, and $n_{r}$ is the number of rollouts. The problem
considered in this paper is as follows.

Problem. Given multiple-trajectory data $\\{[x_{0}^{(k)},u_{0}^{(k)},$
$\dots,$ $x_{\ell-1}^{(k)},u_{\ell-1}^{(k)},x_{\ell}^{(k)}],k\in[n_{r}]\\}$,
estimate the nominal system matrix $[A~{}B]$ and the noise covariance matrix
$(\Sigma_{A},\Sigma_{B})$.

##  3 Identification Algorithm Based on Least-Squares and Multiple-Trajectory
Data

In this section, we propose and study an identification algorithm solving the
considered problem. Section 3.1 studies identifiability of the noise
covariance matrix, paving the way to algorithm design. Consistency of the
algorithm is given by Theorem 3.10 in Section 3.2. Finally, sample complexity
of the algorithm is studied in Section 3.3, and the results are provided in
Theorems 3.16 and 3.17.

###  3.1 Moment Dynamics and Algorithm Design

In this subsection, we propose an algorithm based on multiple trajectories
collected independently to estimate system parameters. Before algorithm
design, the effect of multiplicative noise on moment dynamics is studied, and
identifiability of the noise covariance matrix is clarified.

Taking the expectation of both sides of System (1) and denoting
$\mu_{t}:=\mathds{E}\\{x_{t}\\}$ and $\nu_{t}:=\mathds{E}\\{u_{t}\\}$ yield
the first-moment dynamic of system states (i.e., the dynamic of
$\mathds{E}\\{x_{t}\\}$) as follows,

$\displaystyle\mu_{t+1}=A\mu_{t}+B\nu_{t},~{}t\in\mathds{N}.$ (3)

Denote the vectorization of the second-moment matrices of state, state-input,
and input at time $t$ by
$X_{t}:=\operatorname{vec}(\mathds{E}\\{x_{t}x_{t}^{\intercal}\\})$,
$W_{t}:=\operatorname{vec}(\mathds{E}\\{x_{t}u_{t}^{\intercal}\\})$,
$W_{t}^{\prime}:=\operatorname{vec}(\mathds{E}\\{u_{t}x_{t}^{\intercal}\\})$,
and $U_{t}:=\operatorname{vec}(\mathds{E}\\{u_{t}u_{t}^{\intercal}\\})$. From
the independence of $\bar{A}_{t}$ and $\bar{B}_{t}$, as well as vectorization,
the second-moment dynamic of system states is

$\displaystyle X_{t+1}$ $\displaystyle=(A\otimes A)X_{t}+(B\otimes
A)W_{t}+(A\otimes B)W_{t}^{\prime}+(B\otimes
B)U_{t}+\mathds{E}\\{(\bar{A}_{t}\otimes\bar{A}_{t})\operatorname{vec}(x_{t}x_{t}^{\intercal})\\}+\mathds{E}\\{(\bar{B}_{t}\otimes\bar{B}_{t})\operatorname{vec}(u_{t}u_{t}^{\intercal})\\}$
$\displaystyle=(A\otimes A+\Sigma_{A}^{\prime})X_{t}+(B\otimes
B+\Sigma_{B}^{\prime})U_{t}+(B\otimes A)W_{t}+(A\otimes
B)W_{t}^{\prime},~{}t\in\mathds{N},$ (4)

where
$\Sigma_{A}^{\prime}=\mathds{E}\\{\bar{A}_{t}\otimes\bar{A}_{t}\\}\in\mathds{R}^{n^{2}\times
n^{2}}$ and
$\Sigma_{B}^{\prime}=\mathds{E}\\{\bar{B}_{t}\otimes\bar{B}_{t}\\}\in\mathds{R}^{n^{2}\times
m^{2}}$. The relation between $(\Sigma_{A},\Sigma_{B})$ and
$(\Sigma_{A}^{\prime},\Sigma_{B}^{\prime})$ can be illustrated by
$F(\Sigma_{A}^{\prime},n,n,n,n)=\Sigma_{A}$ and
$F(\Sigma_{B}^{\prime},n,m,n,m)=\Sigma_{B}$, where the reshaping operator
$F(\cdot)$ is defined in the notation section.

An intrinsic identifiability issue arises in the second-moment dynamic (4).
Since $\mathds{E}\\{x_{t}x_{t}^{\intercal}\\}$ is symmetric, $X_{t}$ has
$n(n-1)/2$ pairs of identical entries corresponding to the off-diagonal
entries of $\mathds{E}\\{x_{t}x_{t}^{\intercal}\\}$ (i.e.,
$\mathds{E}\\{x_{t,i}x_{t,j}\\}=\mathds{E}\\{x_{t,j}x_{t,i}\\}$ for all
$i,j\in[n]$). To remove the redundant terms, introduce binary row- and column-
selection matrices, which are also called elimination and duplication matrices
(Magnus and Neudecker,, 1980).

To begin, notice that the redundant entries of $X_{t}$ are associated with the
index set $\\{(j-1)n+i:i,j\in[n],i<j\\}$. Define matrix
$T_{1}\in\mathds{R}^{n^{2}\times n^{2}}$ by replacing the $[(j-1)n+i]$-th row
of $I_{n^{2}}$ by $(\mathbf{e}_{(i-1)n+j}^{n^{2}})^{\intercal}$ for all
$i,j\in[n]$ with $i<j$. Note that $\mathds{E}\\{x_{t,i}x_{t,j}\\}$ is the
$[(j-1)n+i]$-th entry of $X_{t}$, so $X_{t}$ is invariant under $T_{1}$ (i.e.,
$X_{t}=T_{1}X_{t}$). Furthermore, define a binary elimination matrix $P_{1}$
that picks out only the unique entries of $X_{t}$, and a complementary binary
duplication matrix $Q_{1}$ which in turn reconstructs $X_{t}$ from the unique
representation, by repeating the redundant entries in the proper order. These
matrices are defined explicitly as $P_{1}\in\mathds{R}^{[n(n+1)/2]\times
n^{2}}$ by removing the $[(j-1)n+i]$-th row of $I_{n^{2}}$, $i,j\in[n]$ with
$i<j$, and $Q_{1}\in\mathds{R}^{n^{2}\times[n(n+1)/2]}$ by removing the
$[(j-1)n+i]$-th column of $T_{1}$, $i,j\in[n]$ with $i<j$. Then one is able to
freely convert between the full vectorization (with redundant entries) $X_{t}$
and the symmetric vectorization (without redundant entries)
$\tilde{X}_{t}:=\operatorname{svec}(X_{t})$, by employing the linear
transformations defined by the matrices $P_{1}$ and $Q_{1}$:

$\displaystyle\tilde{X}_{t}=P_{1}X_{t},\qquad X_{t}=Q_{1}\tilde{X}_{t}.$

Now apply the same arguments to the second moment of input $U_{t}$: $U_{t}$
has $m(m-1)/2$ pairs of identical entries corresponding to the off-diagonal
entries of $\mathds{E}\\{u_{t}u_{t}^{\intercal}\\}$, so define
$T_{2}\in\mathds{R}^{m^{2}\times m^{2}}$,
$P_{2}\in\mathds{R}^{[m(m+1)/2]\times m^{2}}$, and
$Q_{2}\in\mathds{R}^{m^{2}\times[m(m+1)/2]}$ by replacing $n$ by $m$ in the
definitions of $T_{1},P_{1}$, and $Q_{1}$, respectively.

Applying the symmetric vectorization transformations
$\tilde{X}_{t}=P_{1}X_{t}$ and $\tilde{U}_{t}=P_{2}U_{t}$ yields the second-
moment dynamic with unique entries,

$\displaystyle\tilde{X}_{t+1}$ $\displaystyle=P_{1}X_{t+1}$
$\displaystyle=P_{1}(A\otimes A+\Sigma_{A}^{\prime})X_{t}+P_{1}(B\otimes
B+\Sigma_{B}^{\prime})U_{t}+P_{1}(B\otimes A)W_{t}+P_{1}(A\otimes
B)W_{t}^{\prime}$ $\displaystyle=P_{1}(A\otimes
A+\Sigma_{A}^{\prime})Q_{1}P_{1}X_{t}+P_{1}(B\otimes
B+\Sigma_{B}^{\prime})Q_{2}P_{2}U_{t}+P_{1}(B\otimes A)W_{t}+P_{1}(A\otimes
B)W_{t}^{\prime}$
$\displaystyle=(\tilde{A}+\tilde{\Sigma}_{A}^{\prime})\tilde{X}_{t}+(\tilde{B}+\tilde{\Sigma}_{B}^{\prime})\tilde{U}_{t}+K_{BA}W_{t}+K_{AB}W_{t}^{\prime},$
(5)

where the penultimate equation follows from $T_{1}=Q_{1}P_{1}$ and
$T_{2}=Q_{2}P_{2}$. In the last equation the following notations are
introduced:

$\displaystyle\tilde{A}$ $\displaystyle:=P_{1}(A\otimes
A)Q_{1}\in\mathds{R}^{[n(n+1)/2]\times[n(n+1)/2]},$
$\displaystyle\tilde{\Sigma}_{A}^{\prime}$
$\displaystyle:=P_{1}\Sigma_{A}^{\prime}Q_{1}\in\mathds{R}^{[n(n+1)/2]\times[n(n+1)/2]},$
$\displaystyle\tilde{B}$ $\displaystyle:=P_{1}(B\otimes
B)Q_{2}\in\mathds{R}^{[n(n+1)/2]\times[m(m+1)/2]},$
$\displaystyle\tilde{\Sigma}_{B}^{\prime}$
$\displaystyle:=P_{1}\Sigma_{B}^{\prime}Q_{2}\in\mathds{R}^{[n(n+1)/2]\times[m(m+1)/2]},$
$\displaystyle K_{BA}$ $\displaystyle:=P_{1}(B\otimes
A),~{}K_{AB}:=P_{1}(A\otimes B).$

Note that $\tilde{X}_{t}$ and $\tilde{U}_{t}$ have no redundant entries but
are able to capture the second-moment dynamic of system states. By the
definition of Kronecker product, $\Sigma_{A}^{\prime}$ and
$\Sigma_{B}^{\prime}$ have the following structures.

$\displaystyle\begin{array}[]{lc}\mbox{}&\begin{array}[]{ccccc}&(k-1)n+l&&(l-1)n+k\end{array}\\\
\begin{array}[]{c}{}\hfil\\\ (i-1)n+j\\\ {}\hfil\\\ (j-1)n+i\\\
{}\hfil\end{array}&\left[\begin{array}[]{ccccc}&\vdots&&\vdots&\\\
\cdots&\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jl}\\}&\cdots&\mathds{E}\\{[\bar{A}_{t}]_{il}[\bar{A}_{t}]_{jk}\\}&\cdots\\\
&\vdots&&\vdots&\\\
\cdots&\mathds{E}\\{[\bar{A}_{t}]_{jk}[\bar{A}_{t}]_{il}\\}&\cdots&\mathds{E}\\{[\bar{A}_{t}]_{jl}[\bar{A}_{t}]_{ik}\\}&\cdots\\\
&\vdots&&\vdots&\end{array}\right]\end{array}\begin{array}[]{lc}\mbox{}&\begin{array}[]{ccccc}&(p-1)m+q&&(q-1)m+p\end{array}\\\
\begin{array}[]{c}{}\hfil\\\ {}\hfil\\\ {}\hfil\\\ {}\hfil\\\
{}\hfil\end{array}&\left[\begin{array}[]{ccccc}&\vdots&&\vdots&\\\
\cdots&\mathds{E}\\{[\bar{B}_{t}]_{ip}[\bar{B}_{t}]_{jq}\\}&\cdots&\mathds{E}\\{[\bar{B}_{t}]_{iq}[\bar{B}_{t}]_{jp}\\}&\cdots\\\
&\vdots&&\vdots&\\\
\cdots&\mathds{E}\\{[\bar{B}_{t}]_{jp}[\bar{B}_{t}]_{iq}\\}&\cdots&\mathds{E}\\{[\bar{B}_{t}]_{jq}[\bar{B}_{t}]_{ip}\\}&\cdots\\\
&\vdots&&\vdots&\end{array}\right]\end{array},$ (32)

where $i,j,k,l\in[n]$, $p,q\in[m]$, and $[\bar{A}_{t}]_{ij}$
($[\bar{B}_{t}]_{ip}$) is the $(i,j)$-th entry of $\bar{A}_{t}$ ($(i,p)$-th
entry of $\bar{B}_{t}$). If $i=j$ ($k=l$), the corresponding two rows (two
columns) coincide. The following proposition demonstrates the correspondences
between the entries of $\tilde{\Sigma}_{A}^{\prime}$ and
$\tilde{\Sigma}_{B}^{\prime}$ and those of $\Sigma_{A}^{\prime}$ and
$\Sigma_{B}^{\prime}$, respectively.

######  Proposition 1.

Denote the $(i,j)$-th entry of $\tilde{\Sigma}_{A}^{\prime}$ by
$[\tilde{\Sigma}_{A}^{\prime}]_{ij}$. It holds for $i,j,k,l\in[n]$ with $i<j$
and $k<l$ that

$\displaystyle[\tilde{\Sigma}_{A}^{\prime}]_{(i-1)(n-i/2)+i,(k-1)(n-k/2)+k}=\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{ik}\\},$
$\displaystyle[\tilde{\Sigma}_{A}^{\prime}]_{(i-1)(n-i/2)+i,(k-1)(n-k/2)+l}=2\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{il}\\},$
$\displaystyle[\tilde{\Sigma}_{A}^{\prime}]_{(i-1)(n-i/2)+j,(k-1)(n-k/2)+k}=\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jk}\\},$
$\displaystyle[\tilde{\Sigma}_{A}^{\prime}]_{(i-1)(n-i/2)+j,(k-1)(n-k/2)+l}=\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jl}\\}+\mathds{E}\\{[\bar{A}_{t}]_{il}[\bar{A}_{t}]_{jk}\\}.$

Denote the $(i,j)$-th entry of $\tilde{\Sigma}_{B}^{\prime}$ by
$[\tilde{\Sigma}_{B}^{\prime}]_{ij}$. It holds for $i,j\in[n]$ with $i<j$ and
$p,q\in[m]$ with $p<q$ that

$\displaystyle[\tilde{\Sigma}_{B}^{\prime}]_{(i-1)(n-i/2)+i,(p-1)(m-p/2)+p}=\mathds{E}\\{[\bar{B}_{t}]_{ip}[\bar{B}_{t}]_{ip}\\},$
$\displaystyle[\tilde{\Sigma}_{B}^{\prime}]_{(i-1)(n-i/2)+i,(p-1)(m-p/2)+q}=2\mathds{E}\\{[\bar{B}_{t}]_{ip}[\bar{B}_{t}]_{iq}\\},$
$\displaystyle[\tilde{\Sigma}_{B}^{\prime}]_{(i-1)(n-i/2)+j,(p-1)(m-p/2)+p}=\mathds{E}\\{[\bar{B}_{t}]_{ip}[\bar{B}_{t}]_{jp}\\},$
$\displaystyle[\tilde{\Sigma}_{B}^{\prime}]_{(i-1)(n-i/2)+j,(p-1)(m-p/2)+q}=\mathds{E}\\{[\bar{B}_{t}]_{ip}[\bar{B}_{t}]_{jq}\\}+\mathds{E}\\{[\bar{B}_{t}]_{iq}[\bar{B}_{t}]_{jp}\\}.$

By observing the definitions of $P_{i}$ and $Q_{i}$, $i=1,2$, and the
structures of $\Sigma_{A}^{\prime}$ and $\Sigma_{B}^{\prime}$ shown in (32),
we can get the expressions of the entries of $\tilde{\Sigma}_{A}^{\prime}$ and
$\tilde{\Sigma}_{B}^{\prime}$ as in the proposition. To determine their
positions, note from the definition of $P_{1}$ that all of the $[(j-1)n+i]$-th
rows of $\Sigma_{A}^{\prime}$ are removed during the transformation
$P_{1}\Sigma_{A}^{\prime}$, where $j>i$, $i,j\in[n]$. This means that the
following rows above the $[(i-1)n+j]$-th row of $\Sigma_{A}^{\prime}$, $i\leq
j$, $i,j\in[n]$, are removed: $(i-1)n+1$, $\dots$, $(i-1)n+i-1$, $(i-2)n+1$,
$\dots$, $(i-2)n+i-2$, $\dots$, $n+1$, whose total number is $i(i-1)/2$. Thus,
the $[(i-1)n+j]$-th rows of $\Sigma_{A}^{\prime}$ becomes the
$[(i-1)n+j-i(i-1)/2]$-th row of $\tilde{\Sigma}_{A}^{\prime}$, i.e., the
$[(i-1)(n-i/2)+j]$-th row, where $i\leq j$, $i,j\in[n]$. Applying the same
argument to the columns of $\Sigma_{A}^{\prime}$ and to $\Sigma_{B}^{\prime}$,
we obtain the correspondence given in the proposition.

######  Remark 3.1.

The preceding discussion indicates that $X_{t}$ is determined by $[A~{}B]$ and
$[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]$, and the
proposition shows that there exists a set of equivalent covariance matrices in
the sense that they generate the same second-moment dynamic of system states,
given the nominal matrix $[A~{}B]$. This fact results from that the dynamic of
$X_{t}=Q_{1}\tilde{X}_{t}$ only depends on $[A~{}B]$ and
$[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]$, and is the same
under all $(\Sigma_{1}^{\prime},\Sigma_{2}^{\prime})$ satisfying
$P_{1}\Sigma_{1}^{\prime}Q_{1}=\tilde{\Sigma}_{A}^{\prime}$ and
$P_{2}\Sigma_{2}^{\prime}Q_{2}=\tilde{\Sigma}_{B}^{\prime}$.

From an entry-wise point of view,
$\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jl}\\}$ and
$\mathds{E}\\{[\bar{A}_{t}]_{il}[\bar{A}_{t}]_{jk}\\}$, $i\not=j$ and
$k\not=l$, have a coupled effect on the second-moment dynamic of system
states. We may only estimate the sum of these two entries out of $X_{t}$,
rather than their exact values, since realizations of $\bar{A}_{t}$ and
$\bar{B}_{t}$ are not observed directly but indirectly through their effect on
system states. Fortunately, some entries of $\Sigma_{A}^{\prime}$ and
$\Sigma_{B}^{\prime}$ are identifiable, such as
$\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{ik}\\}$, the variance of
$[\bar{A}_{t}]_{ik}$, and
$\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jk}\\}$, the covariance between
entries in the same column. Similar issues also appear, when estimating
covariance matrices, in topics such as Kalman filtering (Mehra,, 1970; Moghe
et al.,, 2019). Critically, since these identifiable quantities uniquely
generate the second-moment dynamic of system states, it suffices to estimate
$\tilde{\Sigma}_{A}^{\prime}$ and $\tilde{\Sigma}_{B}^{\prime}$ for LQR
design. This fact can be verified by expanding the Bellman equation; we omit
the details to keep the paper concise.

Given $(\Sigma_{A},\Sigma_{B})$ with $\Sigma_{A}\succeq 0$ and
$\Sigma_{B}\succeq 0$ (then
$\tilde{\Sigma}_{A}^{\prime}=P_{1}\Sigma_{A}^{\prime}Q_{1}$ and
$\tilde{\Sigma}_{B}^{\prime}=P_{2}\Sigma_{B}^{\prime}Q_{2}$), the set of
equivalent matrices discussed in Remark 3.1 can be written explicitly as
follows, where positive semidefinite conditions are imposed because
$\Sigma_{A}$ and $\Sigma_{B}$ are covariance matrices,

$\displaystyle S^{*}(\tilde{\Sigma}_{A}^{\prime})$
$\displaystyle:=\Big{\\{}\Sigma_{A}(\alpha)\in\mathds{R}^{n^{2}\times
n^{2}}:\Sigma_{A}(\alpha)\succeq
0,\alpha\in\mathds{R}^{n^{2}(n-1)^{2}/4}\Big{\\}},$ $\displaystyle
S^{*}(\tilde{\Sigma}_{B}^{\prime})$
$\displaystyle:=\Big{\\{}\Sigma_{B}(\beta)\in\mathds{R}^{nm\times
nm}:\Sigma_{B}(\beta)\succeq 0,\beta\in\mathds{R}^{nm(n-1)(m-1)/4}\Big{\\}},$
$\displaystyle S^{*}_{\Sigma}$
$\displaystyle:=S^{*}(\tilde{\Sigma}_{A}^{\prime})\times
S^{*}(\tilde{\Sigma}_{B}^{\prime}),$ (33)

with
$\Sigma_{A}(\alpha):=F(Q_{1}\tilde{\Sigma}_{A}^{\prime}Q_{1}^{\intercal}D_{n}+E_{\alpha},n,n,n,n)$
and
$\Sigma_{B}(\beta):=F(Q_{1}\tilde{\Sigma}_{B}^{\prime}Q_{2}^{\intercal}D_{m}+E_{\beta},n,m,n,m)$.
Here

$\displaystyle E_{\alpha}$
$\displaystyle=\underset{i<j,k<l}{\underset{i,j,k,l\in[n]}{\sum}}\bigg{[}\alpha_{ij,kl}\Big{(}\mathbf{e}_{(i-1)n+j}^{n^{2}}-\mathbf{e}_{(j-1)n+i}^{n^{2}}\Big{)}\Big{(}\mathbf{e}_{(k-1)n+l}^{n^{2}}-\mathbf{e}_{(l-1)n+k}^{n^{2}}\Big{)}^{\intercal}\bigg{]},$
$\displaystyle E_{\beta}$
$\displaystyle=\underset{p,q\in[m],p<q}{\underset{i,j\in[n],i<j}{\sum}}\bigg{[}\beta_{ij,pq}\Big{(}\mathbf{e}_{(i-1)n+j}^{n^{2}}-\mathbf{e}_{(j-1)n+i}^{n^{2}}\Big{)}\Big{(}\mathbf{e}_{(p-1)m+q}^{m^{2}}-\mathbf{e}_{(q-1)m+p}^{m^{2}}\Big{)}^{\intercal}\bigg{]},$

where $\alpha=[\alpha_{ij,kl}]\in\mathds{R}^{n^{2}(n-1)^{2}/4}$,
$\beta=[\beta_{ij,pq}]\in\mathds{R}^{nm(n-1)(m-1)/4}$, $i,j,k,l\in[n]$,
$p,q\in[m]$, $i<j$, $k<l$, $p<q$, $Q_{1}$ and $Q_{2}$ are given before (5),
$D_{n}$ is an $n^{2}$-dimensional diagonal matrix with $[(i-1)n+i]$-th
diagonal entry being $1$ and the rest being $1/2$, $i\in[n]$, and $D_{m}$ is
an $m^{2}$-dimensional diagonal matrix with $[(p-1)m+p]$-th diagonal entry
being $1$ and the rest being $1/2$, $p\in[m]$. Note that $S^{*}_{\Sigma}$ is
given by two inequalities which respectively depend on $\alpha$ and $\beta$.
These two inequalities are linear matrix inequalities (Boyd et al.,, 1994),
since the reshaping operator $F$ is linear. Obviously $S^{*}_{\Sigma}$ is not
empty, because $(\Sigma_{A},\Sigma_{B})$ is one of its elements. The following
example provides an intuitive idea of previous discussions.

######  Example 3.2.

Consider System (1) with $n=2$ and $m=1$, where
$X_{t}=[\mathds{E}\\{x_{t,1}x_{t,1}\\}~{}\mathds{E}\\{x_{t,2}x_{t,1}\\}~{}\mathds{E}\\{x_{t,1}x_{t,2}\\}$
$\mathds{E}\\{x_{t,2}x_{t,2}\\}]^{T}$. So $\mathds{E}\\{X_{t,2}X_{t,1}\\}$ and
$\mathds{E}\\{X_{t,1}X_{t,2}\\}$ are identical and have the same dynamic from
(4). Thus,

$\displaystyle\tilde{X}_{t}=\left[\mathds{E}\\{x_{t,1}x_{t,1}\\}~{}~{}\mathds{E}\\{x_{t,2}x_{t,1}\\}~{}~{}\mathds{E}\\{x_{t,2}x_{t,2}\\}\right]^{T},$
$\displaystyle P_{1}=\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&0&1\end{bmatrix},~{}Q_{1}=\begin{bmatrix}1&0&0\\\ 0&1&0\\\ 0&1&0\\\
0&0&1\end{bmatrix},~{}T_{1}=\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&1&0&0\\\
0&0&0&1\end{bmatrix},~{}P_{2}=Q_{2}=T_{2}=1.$

According to the previously discussed simplification, from

$\displaystyle\Sigma_{A}^{\prime}=\begin{bmatrix}\sigma_{a,11,11}&\sigma_{a,11,12}&\sigma_{a,12,11}&\sigma_{a,12,12}\\\
\sigma_{a,11,21}&\sigma_{a,11,22}&\sigma_{a,12,21}&\sigma_{a,12,22}\\\
\sigma_{a,21,11}&\sigma_{a,21,12}&\sigma_{a,22,11}&\sigma_{a,22,12}\\\
\sigma_{a,21,21}&\sigma_{a,21,22}&\sigma_{a,22,21}&\sigma_{a,22,22}\end{bmatrix},~{}\Sigma_{B}^{\prime}=\begin{bmatrix}\sigma_{b,11}&\sigma_{b,12}&\sigma_{b,21}&\sigma_{b,22}\end{bmatrix}^{T},$

we have that

$\displaystyle\tilde{\Sigma}_{A}^{\prime}=\begin{bmatrix}\sigma_{a,11,11}&2\sigma_{a,11,12}&\sigma_{a,12,12}\\\
\sigma_{a,11,21}&\sigma_{a,11,22}+\sigma_{a,12,21}&\sigma_{a,12,22}\\\
\sigma_{a,21,21}&2\sigma_{a,21,22}&\sigma_{a,22,22}\end{bmatrix},~{}\tilde{\Sigma}_{B}^{\prime}=\begin{bmatrix}\sigma_{b,11}&\sigma_{b,12}&\sigma_{b,22}\end{bmatrix}^{T},$

where $\sigma_{a,ij,kl}=\mathds{E}\\{[\bar{A}_{t}]_{ij}[\bar{A}_{t}]_{kl}\\}$
and $\sigma_{b,ij}=\mathds{E}\\{[\bar{B}_{t}]_{i},[\bar{B}_{t}]_{j}\\}$, and
that

$\displaystyle\tilde{A}=\begin{bmatrix}a_{11}a_{11}&a_{11}a_{12}+a_{12}a_{11}&a_{12}a_{12}\\\
a_{11}a_{21}&a_{11}a_{22}+a_{12}a_{21}&a_{12}a_{22}\\\
a_{21}a_{21}&a_{21}a_{22}+a_{22}a_{21}&a_{22}a_{22}\end{bmatrix},~{}\tilde{B}=\begin{bmatrix}b_{1}b_{1}&b_{1}b_{2}&b_{2}b_{2}\end{bmatrix}^{T},$
$\displaystyle K_{AB}=\begin{bmatrix}a_{11}b_{1}&a_{12}b_{1}\\\
a_{21}b_{1}&a_{22}b_{2}\\\
a_{21}b_{2}&a_{22}b_{2}\end{bmatrix},~{}K_{BA}=\begin{bmatrix}a_{11}b_{1}&a_{12}b_{1}\\\
a_{11}b_{2}&a_{12}b_{2}\\\ a_{21}b_{2}&a_{22}b_{2}\end{bmatrix}.$

In this example, $\Sigma_{B}$ is unique, but based on (33) the covariance
matrix $\Sigma_{A}(\alpha)$, equivalent to $\Sigma_{A}$, is given by

$\displaystyle\Sigma_{A}(\alpha)=\begin{bmatrix}\sigma_{a,11,11}&\sigma_{a,11,21}&\sigma_{a,11,12}&\sigma_{a,11,22}+\alpha\\\
\sigma_{a,21,11}&\sigma_{a,21,21}&\sigma_{a,21,12}-\alpha&\sigma_{a,21,22}\\\
\sigma_{a,12,11}&\sigma_{a,12,21}-\alpha&\sigma_{a,12,12}&\sigma_{a,12,22}\\\
\sigma_{a,22,11}+\alpha&\sigma_{a,22,21}&\sigma_{a,22,12}&\sigma_{a,22,22}\end{bmatrix},$

where $\alpha\in\mathds{R}$ is such that $\Sigma_{A}(\alpha)\succeq 0$.

######  Example 3.3.

Consider System (2) with $\bar{A}_{t}=\sum_{i=1}^{r}A_{i}p_{i,t}$,
$\bar{B}_{t}=\sum_{j=1}^{s}B_{j}q_{j,t}$. Hence,

$\displaystyle\Sigma_{A}=\sum_{i=1}^{r}\mathds{E}\\{p_{i,t}^{2}\\}\operatorname{vec}(A_{i})\operatorname{vec}(A_{i})^{\intercal},~{}\Sigma_{A}^{\prime}=\sum_{i=1}^{r}\mathds{E}\\{p_{i,t}^{2}\\}A_{i}\otimes
A_{i},~{}\tilde{\Sigma}_{A}^{\prime}=\sum_{i=1}^{r}\mathds{E}\\{p_{i,t}^{2}\\}P_{1}(A_{i}\otimes
A_{i})Q_{1},$
$\displaystyle\Sigma_{B}=\sum_{j=1}^{s}\mathds{E}\\{q_{j,t}^{2}\\}\operatorname{vec}(B_{j})\operatorname{vec}(B_{j})^{\intercal},~{}\Sigma_{B}^{\prime}=\sum_{j=1}^{s}\mathds{E}\\{q_{j,t}^{2}\\}B_{j}\otimes
B_{j},~{}\tilde{\Sigma}_{B}^{\prime}=\sum_{j=1}^{s}\mathds{E}\\{q_{j,t}^{2}\\}P_{1}(B_{j}\otimes
B_{j})Q_{2}.$

Suppose that for $A_{i}$, $i\in[r]$, there exist $k_{i}$, $l_{i}\in[n]$ such
that $[A_{i}]_{k_{i},l_{i}}\not=0$ and $[A_{j}]_{k_{i},l_{i}}$ for all
$j\in[r]\setminus\\{i\\}$. That is, the $(k_{i},l_{i})$-th entry of $A_{i}$ is
nonzero but the $(k_{i},l_{i})$-th entry of $A_{j}$ is zero for all $j\not=i$.
Then
$\mathds{E}\\{[\bar{A}_{t}]_{k_{i},l_{i}}[\bar{A}_{t}]_{k_{i},l_{i}}\\}=[A_{i}]_{k_{i},l_{i}}^{2}\mathds{E}\\{p_{i,t}^{2}\\}=[A_{i}]_{k_{i},l_{i}}^{2}\sigma_{i}^{2}$.
From Proposition 1 we know that $\sigma_{i}^{2}=\mathds{E}\\{p_{i,t}^{2}\\}$
can be uniquely determined if second-moment dynamic (4), or
$\tilde{\Sigma}_{A}^{\prime}$, is given. A similar conclusion holds for
$\\{q_{j,t}\\}$. However there are also situations where $\sigma_{i}^{2}$
cannot be uniquely determined. For instance, assume that $r\geq 2$ and for all
$i\in[r]$, $[A_{i}]_{11}\not=0$ but all other entries of $A_{i}$ are zero.
Then we only have a single equation
$\sum_{i=1}^{r}[A_{i}]_{11}^{2}\sigma_{i}^{2}=[\tilde{\Sigma}_{A}^{\prime}]_{11}$
for $\\{\sigma_{i}^{2}\\}$.

As shown in Example 3.2, given $(\Sigma_{A},\Sigma_{B})$ with
$\Sigma_{A}\succeq 0$ and $\Sigma_{B}\succeq 0$, the set $S^{*}_{\Sigma}$ is
not empty but may have infinitely many elements, resulting in unidentifiable
entries $\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jl}\\}$ and
$\mathds{E}\\{[\bar{B}_{t}]_{ip}[\bar{B}_{t}]_{jq}\\}$, $i\not=j$, $k\not=l$,
$p\not=q$, $i,j,k,l\in[n]$, $p,q\in[m]$. The following proposition gives
several conditions under which the covariance matrix of the multiplicative
noise can or cannot be uniquely determined from
$[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]$.

######  Proposition 2.

Given $(\Sigma_{A},\Sigma_{B})$ with $\Sigma_{A}\succeq 0$ and
$\Sigma_{B}\succeq 0$,
$\tilde{\Sigma}_{A}^{\prime}=P_{1}\Sigma_{A}^{\prime}Q_{1}$ and
$\tilde{\Sigma}_{B}^{\prime}=P_{2}\Sigma_{B}^{\prime}Q_{2}$, the following
results hold.
(i) If $n=m=1$, then $S^{*}_{\Sigma}$ has a unique element. If $m=1$, then
$S^{*}(\tilde{\Sigma}_{B}^{\prime})$ has a unique element. If $n\geq 2$ and
$\Sigma_{A}\succ 0$ (resp. $m\geq 2$ and $\Sigma_{B}\succ 0$), then
$S^{*}(\tilde{\Sigma}_{A}^{\prime})$ (resp.
$S^{*}(\tilde{\Sigma}_{B}^{\prime})$) has infinitely many elements. As a
result, under either condition, $S^{*}_{\Sigma}$ has infinitely many elements.
(ii) If $S^{*}(\tilde{\Sigma}_{A}^{\prime})$ has infinitely many elements,
then $S^{*}(\tilde{\Sigma}_{A}^{\prime})\cap T_{A}$ has a unique element,
where

$\displaystyle T_{A}$
$\displaystyle:=\Big{\\{}\Sigma\in\mathds{R}^{n^{2}\times
n^{2}}:\gamma_{ij,kl}\Sigma_{(k-1)n+i,(l-1)n+j}+\delta_{ij,kl}\Sigma_{(l-1)n+i,(k-1)n+j}=\tau_{ij,kl},i<j,k<l,i,j,k,l\in[n]\Big{\\}},$

with constants $\gamma_{ij,kl},\delta_{ij,kl},\tau_{ij,kl}\in\mathds{R}$ and
$\gamma_{ij,kl}\not=\delta_{ij,kl}$ for all $i<j$, $k<l$, $i,j,k,l\in[n]$. The
same result holds for $S^{*}(\tilde{\Sigma}_{B}^{\prime})$ by modifying the
definition of $T_{A}$ according to the dimension of $\Sigma_{B}$.

The first two conclusions of (i) are trivial. If $n\geq 2$, then
$\tilde{\Sigma}_{A}^{\prime}$ has entries of the form
$\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jl}\\}+\mathds{E}\\{[\bar{A}_{t}]_{il}[\bar{A}_{t}]_{jk}\\}$.
Since $\Sigma_{A}\succ 0$, its minimum eigenvalue is larger than zero. Note
that from the definition of $S^{*}(\tilde{\Sigma}_{A}^{\prime})$ there exists
$\alpha^{*}$ such that $\Sigma_{A}=\Sigma_{A}(\alpha^{*})$. Because the
eigenvalues of a matrix depend continuously on its entries (Theorem 2.4.9.2 of
Horn and Johnson, (2012)), $\Sigma_{A}(\alpha^{*}+\varepsilon)$ is still a
positive definite matrix for small enough $\varepsilon>0$. This proves the
last result in (i). From (i), we know that if
$S^{*}(\tilde{\Sigma}_{A}^{\prime})$ has infinitely many elements, then $n\geq
2$. To show (ii), just note that if
$\Sigma=\operatorname{vec}(\bar{A})\operatorname{vec}(\bar{A})^{\intercal}$
for some $A\in\mathds{R}^{n\times n}$, then the $[(k-1)n+i,(l-1)n+j]$-th entry
of $\Sigma$ is $[\bar{A}]_{ik}[\bar{A}]_{jl}$ and the $[(l-1)n+i,(k-1)n+j]$-th
entry is $[\bar{A}]_{il}[\bar{A}]_{jk}$, $i\not=j$, $k\not=l$. Hence if
$\gamma_{ij,kl}\not=\delta_{ij,kl}$ then we have two linearly independent
equations for $[\bar{A}]_{ik}[\bar{A}]_{jl}$ and
$[\bar{A}]_{il}[\bar{A}]_{jk}$ (the other one from Proposition 1 is
$[\bar{A}]_{ik}[\bar{A}]_{jl}+[\bar{A}]_{il}[\bar{A}]_{jk}=[\tilde{\Sigma}_{A}^{\prime}]_{(i-1)(n-i/2)+j,(k-1)(n-k/2)+l}$).
So these entries can be uniquely determined, and the conclusion follows.

######  Remark 3.4.

The first part of the proposition shows that if $\Sigma_{A}\succ 0$ or
$\Sigma_{B}\succ 0$ and $n\geq m\geq 2$, then it is impossible to uniquely
determine $(\Sigma_{A},\Sigma_{B})$ only based on second-moment dynamic (5).
However the second part indicates that more conditions imposed on the
covariance matrix can make all entries of $\Sigma_{A}$ and $\Sigma_{B}$
identifiable. The set $T_{A}$ introduces additional constraints for
$\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jl}\\}$, $i\not=j$, $k\not=l$.
For example, if entries in $\bar{A}_{t}$ are mutually independent, then
$\Sigma_{A}$ is diagonal. In this case, it holds that
$\mathds{E}\\{[\bar{A}_{t}]_{ik}[\bar{A}_{t}]_{jl}\\}-\mathds{E}\\{[\bar{A}_{t}]_{il}[\bar{A}_{t}]_{jk}\\}=0$,
$i\not=j$, $k\not=l$, and hence the covariance matrix of $\bar{A}_{t}$ is
uniquely determined.

Now we are ready to propose our estimation algorithm. Following the previous
discussion, we introduce an algorithm based on the first- and second-moment
dynamics (3) and (5). Since the exact moment dynamics are unavailable, we
average over multiple independent rollouts to obtain their estimates. To get
persistently exciting inputs, it is necessary to design their first and second
moments in advance, in either a deterministic or a stochastic way. For
example, generate the two moments from standard Gaussian and Wishart
distributions (Gupta and Nagar,, 2018), respectively, or set them
periodically. The initial states of different rollouts are assumed to be
i.i.d. subject to a same distribution $\mathcal{X}_{0}$ with finite second
moment (see Section 3.2.2). The overall algorithm is shown in Algorithm 1,
where the superscript $(k)$ represents the $k$-th rollout. Note that Algorithm
1 is different from classic recursive identification algorithms. The recursive
least-squares algorithm (Lai and Wei,, 1982; Chen and Guo,, 2012), for
example, uses only one trajectory of a system. In contrast, Algorithm 1 is
based on multiple trajectories with finite length.

Based on the estimates $\hat{\tilde{\Sigma}}_{A}^{\prime}$ and
$\hat{\tilde{\Sigma}}_{B}^{\prime}$, it is able to obtain an estimate
$\hat{S}^{*}_{\Sigma}$ of the equivalent class (33), via replacing
$\tilde{\Sigma}_{A}^{\prime}$ and $\tilde{\Sigma}_{B}^{\prime}$ in the
definition (33) by their estimates. If the linear matrix inequalities are
infeasible (i.e., $\hat{S}^{*}_{\Sigma}=\emptyset$), then project the
estimates onto the positive semidefinite cone. However this situation is
unlikely to happen when $n_{r}$ is large, because of the consistency of
Algorithm 1 given in the next section.

Algorithm 1
Multiple-trajectory averaging least-squares (MALS)

1:Input: Rollout length $\ell$ and the number of rollouts $n_{r}$.

2:Output: $[\hat{A}~{}\hat{B}]$,
$[\hat{\tilde{\Sigma}}_{A}^{\prime}~{}\hat{\tilde{\Sigma}}_{B}^{\prime}]$.

3:// Control-input design

4:for $t$ from $0$ to $\ell-1$ do

5: Generate $\nu_{t}\in\mathds{R}^{m}$ and $\bar{U}_{t}\in\mathds{R}^{m\times
m}$ with $\bar{U}_{t}\succeq 0$.

6:end for

7:// Multiple-trajectory collection

8:for $k$ from $1$ to $n_{r}$ do

9: Generate $x_{0}^{(k)}$ independently from the initial multivariate
distribution $\mathcal{X}_{0}$.

10: for $t$ from $0$ to $\ell-1$ do

11: Generate $u_{t}^{(k)}$ independently from a multivariate distribution with
first moment $\nu_{t}$ and second central

12: moment $\bar{U}_{t}$,

13:
$x_{t+1}^{(k)}=(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}$.

14: end for

15:end for

16:// Least-squares estimation

17:for $t$ from $0$ to $\ell$ do

18: Compute  $\displaystyle\hat{\mu}_{t}$
$\displaystyle:=\frac{1}{n_{r}}\sum_{k=1}^{n_{r}}x_{t}^{(k)},$
$\displaystyle\hat{\tilde{X}}_{t}$
$\displaystyle:=\frac{1}{n_{r}}P_{1}\operatorname{vec}\left(\sum_{k=1}^{n_{r}}x_{t}^{(k)}(x_{t}^{(k)})^{\intercal}\right),$
$\displaystyle\hat{W}_{t}$
$\displaystyle:=\frac{1}{n_{r}}\operatorname{vec}\left(\sum_{k=1}^{n_{r}}x_{t}^{(k)}\nu_{t}^{\intercal}\right)=\operatorname{vec}(\hat{\mu}_{t}\nu_{t}^{\intercal}),$
$\displaystyle\hat{W}_{t}^{\prime}$
$\displaystyle:=\frac{1}{n_{r}}\operatorname{vec}\left(\sum_{k=1}^{n_{r}}\nu_{t}{x_{t}^{(k)}}^{\intercal}\right)=\operatorname{vec}(\nu_{t}\hat{\mu}_{t}^{\intercal}),$
$\displaystyle\tilde{U}_{t}$
$\displaystyle:=P_{2}\operatorname{vec}(\bar{U}_{t}+\nu_{t}\nu_{t}^{\intercal}).$

19:end for

20:$[\hat{A}~{}\hat{B}]=\underset{{[A~{}B]}}{\text{argmin}}\\{\sum_{t=0}^{\ell-1}\|\hat{\mu}_{t+1}-(A\hat{\mu}_{t}+B\nu_{t})\|_{2}^{2}\\}$,

21:Compute $\hat{\tilde{A}}=P_{1}(\hat{A}\otimes\hat{A})Q_{1}$,
$\hat{\tilde{B}}=P_{1}(\hat{B}\otimes\hat{B})Q_{2}$,
$\hat{K}_{BA}=P_{1}(\hat{B}\otimes\hat{A})$, and
$\hat{K}_{AB}=P_{1}(\hat{A}\otimes\hat{B})$, where $P_{1},P_{2},Q_{1}$, and
$Q_{2}$ are given before (5),

22:$[\hat{\tilde{\Sigma}}_{A}^{\prime}~{}\hat{\tilde{\Sigma}}_{B}^{\prime}]=\underset{[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]}{\text{argmin}}\\{\sum_{t=0}^{\ell-1}\|\hat{\tilde{X}}_{t+1}$
$-[\tilde{A}\hat{\tilde{X}}_{t}+K_{BA}\hat{W}_{t}$
$+K_{AB}\hat{W}_{t}^{\prime}+\tilde{B}\tilde{U}_{t}+\tilde{\Sigma}_{A}^{\prime}\hat{\tilde{X}}_{t}+\tilde{\Sigma}_{B}^{\prime}\tilde{U}_{t}]\|_{2}^{2}\\}$.

###  3.2 Performance of Algorithm 1

This section analyzes performance of Algorithm 1 by investigating the moment
dynamics (3) and (5).

####  3.2.1 Moment Dynamics and Input Design

Provided that $\mu_{t}$ and $\tilde{X}_{t}$ are known, it is possible to
recover the parameters via least-squares as in lines $14$-$16$ in Algorithm 1.
Denote

$\displaystyle\mathbf{Y}:=[\mu_{\ell}~{}\cdots~{}\mu_{1}],~{}\mathbf{Z}:=\begin{bmatrix}\mu_{\ell-1}&\cdots&\mu_{0}\\\
\nu_{\ell-1}&\cdots&\nu_{0}\end{bmatrix},~{}\mathbf{C}:=[C_{\ell}~{}\cdots~{}C_{1}],~{}\mathbf{D}:=\begin{bmatrix}\tilde{X}_{\ell-1}&\cdots&\tilde{X}_{0}\\\
\tilde{U}_{\ell-1}&\cdots&\tilde{U}_{0}\end{bmatrix},$ (34)

where
$C_{t}=\tilde{X}_{t}-\big{(}\tilde{A}\tilde{X}_{t-1}+K_{BA}W_{t-1}+K_{AB}W_{t-1}^{\prime}+\tilde{B}\tilde{U}_{t-1}\big{)}$,
$1\leq t\leq\ell$. Then closed-form solutions of the least-squares problems
are

$\displaystyle\big{[}\hat{A}~{}\hat{B}\big{]}=\mathbf{Y}\mathbf{Z}^{\intercal}(\mathbf{Z}\mathbf{Z}^{\intercal})^{\dagger},~{}\Big{[}\hat{\tilde{\Sigma}}_{A}^{\prime}~{}\hat{\tilde{\Sigma}}_{B}^{\prime}\Big{]}=\mathbf{C}\mathbf{D}^{\intercal}(\mathbf{D}\mathbf{D}^{\intercal})^{\dagger},$

where $\dagger$ represents the pseudoinverse. When the inverse matrices exist,
the solutions are identical to true values; that is,
$[\hat{A}~{}\hat{B}]=[A~{}B]$ and
$[\hat{\tilde{\Sigma}}_{A}^{\prime}~{}\hat{\tilde{\Sigma}}_{B}^{\prime}]=[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]$.
Hence, the first question towards the consistency of Algorithm 1 is whether
the matrices $\mathbf{Z}\mathbf{Z}^{\intercal}$ and
$\mathbf{D}\mathbf{D}^{\intercal}$ are invertible. As to be shown, designing a
proper input sequence ensures this invertibility, if systems $(A,B)$ and
$(\tilde{A}+\tilde{\Sigma}_{A}^{\prime},\tilde{B}+\tilde{\Sigma}_{B}^{\prime})$
are controllable, and the rollout length $\ell$ is large enough.

######  Proposition 3.

Suppose that $\ell\geq n+m$ and $(A,B)$ is controllable. For fixed
$\mu_{0}\in\mathds{R}^{n}$, the matrix $\mathbf{Z}$ has full row rank, and
consequently $\mathbf{Z}\mathbf{Z}^{\intercal}$ is invertible, for almost all
$[\nu_{0}^{\intercal}~{}\cdots~{}\nu_{\ell-1}^{\intercal}]^{\intercal}\in\mathds{R}^{m\ell}$.

See Appendix A.

######  Remark 3.5.

The proposition shows that for large enough rollout length, the full row
rankness of $\mathbf{Z}$ can be guaranteed for almost all
$[\nu_{0}^{\intercal}~{}\cdots~{}\nu_{\ell-1}^{\intercal}]^{\intercal}\in\mathds{R}^{m\ell}$.
The controllability of $(A,B)$ plays a key role in the proof, similar to
classic results on identification of linear systems (Chen and Guo,, 2012). The
condition $\ell\geq n+m$ is necessary for the invertibility of
$\mathbf{Z}\mathbf{Z}^{\intercal}$. This lower bound is much smaller than that
given in Xing et al., (2020). According to the proposition,
$\mathbf{Z}\mathbf{Z}^{\intercal}$ is invertible with probability one if the
first moments of inputs are generated i.i.d. from a distribution absolutely
continuous with respect to Lebesgue measure (e.g., Gaussian distribution or
uniform distribution). This proposition can be seen as a generalization of the
single-input case studied in Schmidt et al., (2005).

######  Proposition 4.

Suppose that $\ell\geq[n(n+1)+m(m+1)]/2$ and
$(\tilde{A}+\tilde{\Sigma}_{A}^{\prime},\tilde{B}+\tilde{\Sigma}_{B}^{\prime})$
is controllable. For fixed $\mu_{0}\in\mathds{R}^{n}$ and
$\tilde{X}_{0}\in\mathds{R}^{n(n+1)/2}$, the matrix $\mathbf{D}$ has full row
rank, and consequently $\mathbf{D}\mathbf{D}^{\intercal}$ is invertible, for
almost all
$[\nu_{0}^{\intercal}~{}\cdots~{}\nu_{\ell-1}^{\intercal}~{}\operatorname{svec}(\bar{U}_{0})^{\intercal}$
$\cdots$
$\operatorname{svec}(\bar{U}_{\ell-1}))^{\intercal}]^{\intercal}\in\mathds{R}^{\ell
m(m+3)/2}$, where $\bar{U}_{t}$ is defined in line 2 of Algorithm 1.

See Appendix B.

######  Remark 3.6.

The controllability condition in Proposition 4 reflects the nature of the
multiplicative noise (i.e., coupling between $\bar{A}_{t}$ and $x_{t}$, and
that between $\bar{B}_{t}$ and $u_{t}$). The result indicates that a
controllability condition on (5) may be necessary to ensure successful
identification. The lower bound for $\ell$ is necessary for the invertibility
of $\mathbf{D}\mathbf{D}^{\intercal}$, and is much smaller than that given in
Xing et al., (2020). As in Algorithm 1,
$\tilde{U}_{t}=\operatorname{svec}(\bar{U}_{t}+\nu_{t}\nu_{t}^{\intercal})$,
so random generation of $\nu_{t}$ and $\bar{U}_{t}$ ensures
$\mathbf{D}\mathbf{D}^{\intercal}$ is invertible with probability one.

We summarize the preceding two results in the following corollary.

######  Corollary 3.7.

Suppose that $\ell\geq[n(n+1)+m(m+1)]/2$, and both $(A,B)$ and
$(\tilde{A}+\tilde{\Sigma}_{A}^{\prime},\tilde{B}+\tilde{\Sigma}_{B}^{\prime})$
are controllable. For fixed $\mu_{0}\in\mathds{R}^{n}$ and
$\tilde{X}_{0}\in\mathds{R}^{n(n+1)/2}$, the matrices
$\mathbf{Z}\mathbf{Z}^{\intercal}$ and $\mathbf{D}\mathbf{D}^{\intercal}$ are
invertible, for almost all
$[\nu_{0}^{\intercal}~{}\cdots~{}\nu_{\ell-1}^{\intercal}~{}\operatorname{svec}(\bar{U}_{0})^{\intercal}$
$\cdots$
$\operatorname{svec}(\bar{U}_{\ell-1}))^{\intercal}]^{\intercal}\in\mathds{R}^{\ell
m(m+3)/2}$, where $\bar{U}_{t}$ is defined in line 2 of Algorithm 1.

######  Remark 3.8.

The corollary implies that the existence of
$(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}$ and
$(\mathbf{D}\mathbf{D}^{\intercal})^{-1}$ can be guaranteed with probability
one, as long as both $\nu_{t}$ and $\bar{U}_{t}$ are independently generated
from distributions that is absolutely continuous with respect to Lebesgue
measure. For example, the entries of $\nu_{t}$ are generated i.i.d. from a
non-degenerate Gaussian distribution and then $\bar{U}_{t}$ is generated
i.i.d. from a non-degenerate Wishart distribution, $0\leq t\leq\ell-1$.

####  3.2.2 Asymptotic Consistency

In this subsection, we assume that the expectations and covariance matrices of
inputs have been generated as discussed in the previous section, and that both
$\mathbf{Z}\mathbf{Z}^{\intercal}$ and $\mathbf{D}\mathbf{D}^{\intercal}$ have
been designed to be invertible. The closed-form estimates generated by
Algorithm 1 are

$\displaystyle\big{[}\hat{A}~{}\hat{B}\big{]}$
$\displaystyle=\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger},$
(35)
$\displaystyle\Big{[}\hat{\tilde{\Sigma}}_{A}^{\prime}~{}\hat{\tilde{\Sigma}}_{B}^{\prime}\Big{]}$
$\displaystyle=\hat{\mathbf{C}}\hat{\mathbf{D}}^{\intercal}(\hat{\mathbf{D}}\hat{\mathbf{D}}^{\intercal})^{\dagger},$
(36)

where

$\displaystyle\hat{\mathbf{Y}}:=\big{[}\hat{\mu}_{\ell}~{}\cdots~{}\hat{\mu}_{1}\big{]},\quad\hat{\mathbf{Z}}:=\begin{bmatrix}\hat{\mu}_{\ell-1}&\cdots&\hat{\mu}_{0}\\\
\nu_{\ell-1}&\cdots&\nu_{0}\end{bmatrix},$ (37)
$\displaystyle\hat{\mathbf{C}}:=\big{[}\hat{C}_{\ell}~{}\cdots~{}\hat{C}_{1}\big{]},\quad\hat{\mathbf{D}}:=\begin{bmatrix}\hat{\tilde{X}}_{\ell-1}&\cdots&\hat{\tilde{X}}_{0}\\\
\tilde{U}_{\ell-1}&\cdots&\tilde{U}_{0}\end{bmatrix},$ (38)

and
$\hat{C}_{t}=\hat{\tilde{X}}_{t}-\big{(}\hat{\tilde{A}}\hat{\tilde{X}}_{t-1}+\hat{K}_{BA}\hat{W}_{t-1}+\hat{K}_{AB}\hat{W}_{t-1}^{\prime}+\hat{\tilde{B}}\tilde{U}_{t-1}\big{)}$,
$1\leq t\leq\ell$. Here $\hat{\tilde{A}}$, $\hat{\tilde{B}}$, $\hat{K}_{AB}$,
and $\hat{K}_{BA}$ are estimates of $\tilde{A}$, $\tilde{B}$, $K_{AB}$, and
$K_{BA}$, obtained from $\hat{A}$ and $\hat{B}$ given by Algorithm 1. The
estimates depend on the number of rollouts $n_{r}$, which is omitted for
convenience. For the convergence result, we present the following assumptions.

######  Assumption 1.

For all rollouts indexed by $k\in[n_{r}]$, the below conditions hold.
(i) The rollout length is $\ell\geq[n(n+1)+m(m+1)]/2$.
(ii) The initial states $x_{0}^{(k)},k\in[n_{r}]$, are i.i.d. subject to the
same distribution $\mathcal{X}_{0}$ with finite second moment, and are
independent of the multiplicative noise and inputs.
(iii) $\\{\bar{A}_{t}^{(k)},0\leq t\leq\ell,k\in[n_{r}]\\}$ and
$\\{\bar{B}_{t}^{(k)},0\leq t\leq\ell,k\in[n_{r}]\\}$, are i.i.d. sequences
respectively and are mutually independent, both with zero mean and finite
second moments (i.e., $\mathds{E}\\{\bar{A}_{t}^{(k)}\\}$ and
$\mathds{E}\\{\bar{B}_{t}^{(k)}\\}$ are zero matrices, and
$\|\Sigma_{A}\|_{2},\|\Sigma_{B}\|_{2}<\infty$).
(iv) The parameters of inputs are given by lines $1$-$3$ of Algorithm 1, and
the inputs are generated, according to line $7$ of Algorithm 1. The inputs and
noise are independent.
(v) Both $\mathbf{Z}\mathbf{Z}^{\intercal}$ and
$\mathbf{D}\mathbf{D}^{\intercal}$ are invertible.

######  Remark 3.9.

From Corollary 3.7, the lower bound of the rollout length in Assumption 1 (i)
is necessary for estimating the noise covariance matrix, whereas, from
Proposition 3, trajectories with length $\ell\geq n+m$ may be enough for
estimating the nominal system matrix. The initial states of different
trajectories need not start with the same value, but it is required that they
have the same first and second moments (Assumption 1 (ii)). The mutual
independence of noise at different time steps in one trajectory is a standard
assumption (Assumption 1 (iii)), but the results in this paper still hold, if
the noise sequence in the same trajectory is dependent, but the noise
sequences in different trajectories are mutually independent and the noise has
zero mean and the same second moment. The physical meaning of the independence
between the noise and the inputs in Assumption 1 (iv) is that the former is an
intrinsic part of the system and cannot be influenced by inputs. To keep the
analysis concise, we separately discuss the input design (Section 3.2.1) and
the performance of Algorithm 1. Assumption 1 (v) indicates that the input
design yields invertible $\mathbf{Z}\mathbf{Z}^{\intercal}$ and
$\mathbf{D}\mathbf{D}^{\intercal}$, but note that it implicitly assumes the
controllability of the first- and second-moment dynamics of system states.

Under Assumption 1 the rollouts $[x_{0}^{(k)},\dots,x_{l}^{(k)}]$,
$k\in[n_{r}]$, are i.i.d., so the following consistency result can be obtained
from strong law of large numbers.

######  Theorem 3.10.

(Consistency) Suppose that Assumption 1 holds, then the estimators (35)-(36)
are asymptotically consistent, namely,

$\displaystyle\big{[}\hat{A}~{}\hat{B}\big{]}\to[A~{}B]\ \ \text{ and }\ \
\Big{[}\hat{\tilde{\Sigma}}_{A}^{\prime}~{}\hat{\tilde{\Sigma}}_{B}^{\prime}\Big{]}\to\big{[}\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}\big{]},$

with probability one as the number of rollouts $n_{r}\to\infty$.

See Appendix C.

######  Remark 3.11.

This theorem indicates that consistency of Algorithm 1 may hold even when the
rollout length is relatively small. In Di and Lamperski, (2021), the
estimation of the first and second moments of multiplicative noise is
decoupled, whereas here the estimate of
$[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]$ relies on
$[\hat{A}~{}\hat{B}]$. The coupling exists because here the noise covariance
matrix, which from definition depends on the mean of the noise, is estimated.
Note that $\ell$ is assumed to be fixed and we do not consider the case where
$\ell\to\infty$, since an averaging step is used in Algorithm 1. Study of the
case with increasing rollout length is left to future work.

###  3.3 Finite-Sample Analysis

This subsection studies finite-sample performance of Algorithm 1,
demonstrating its non-asymptotic behavior. The existence of multiplicative
noise complicates the analysis, so the following assumptions, ensuring that
the system is bounded a.s., are introduced.

######  Assumption 2.

For all rollouts indexed by $k\in[n_{r}]$, the following conditions hold.
(i) The initial state is bounded a.s. for all $k\in[n_{r}]$ as

$\displaystyle\|x_{0}^{(k)}\|\leq c_{X}<\infty.$

(ii) The inputs are bounded a.s. for all $0\leq t\leq\ell-1$ and $k\in[n_{r}]$
as

$\displaystyle\|u_{t}^{(k)}\|\leq c_{U}<\infty.$

(iii) The multiplicative noise, $\bar{A}_{t}^{(k)}$ and $\bar{B}_{t}^{(k)}$,
is bounded a.s. for all $0\leq t\leq\ell-1$ and $k\in[n_{r}]$ as

$\displaystyle\|\bar{A}_{t}^{(k)}\|_{2}\leq
c_{\bar{A}}<\infty,~{}\|\bar{B}_{t}^{(k)}\|_{2}\leq c_{\bar{B}}<\infty.$

######  Remark 3.12.

The assumption of bounded multiplicative noise is reasonable for physical
systems, which cannot have infinite variations. For example, in interconnected
systems, the noise represents randomly varying topologies of subsystems, and
is naturally bounded.

Introduce the state- and input-deviation quantities

$\displaystyle e_{t}^{(k)}$
$\displaystyle:=x_{t}^{(k)}-\mathds{E}[x_{t}^{(k)}]=x_{t}^{(k)}-\mu_{t},$
$\displaystyle d_{t}^{(k)}$
$\displaystyle:=u_{t}^{(k)}-\mathds{E}[u_{t}^{(k)}]=u_{t}^{(k)}-\nu_{t}.$

The next proposition is a natural consequence of Assumption 2.

######  Proposition 5.

Under Assumption 2, the following results hold.
(i) The initial state-deviation is bounded a.s. for all rollouts $k\in[n_{r}]$
as

$\displaystyle\|e_{0}^{(k)}\|\leq c_{\mu}<\infty.$

(ii) The outer product initial state deviation is bounded a.s. for all
rollouts $k\in[n_{r}]$ as

$\displaystyle\left\|\operatorname{vec}\left(x_{0}^{(k)}(x_{0}^{(k)})^{\intercal}-\mathds{E}\left[x_{0}^{(k)}(x_{0}^{(k)})^{\intercal}\right]\right)\right\|\leq
c_{\Delta X}.$

(iii) The input-deviations are bounded a.s. for all $0\leq t\leq\ell-1$ and
$k\in[n_{r}]$ as

$\displaystyle\|d_{t}^{(k)}\|\leq c_{\nu}<\infty.$

(iv) The Kronecker products of $\bar{A}_{t}^{(k)}$ and $\bar{B}_{t}^{(k)}$ are
bounded a.s. for all $0\leq t\leq\ell-1$ and $k\in[n_{r}]$ as

$\displaystyle\|\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)}-\Sigma_{A}^{\prime}\|_{2}\leq
c_{\Sigma_{A}^{\prime}},$
$\displaystyle\|\bar{B}_{t}^{(k)}\otimes\bar{B}_{t}^{(k)}-\Sigma_{B}^{\prime}\|_{2}\leq
c_{\Sigma_{B}^{\prime}}.$

######  Remark 3.13.

This proposition captures the deviations of random components of System (1)
from their expectations. Using the bounds in Assumption 2 one could upper-
bound these deviations, for instance,

$\displaystyle\left\|\operatorname{vec}\left(x_{0}^{(k)}(x_{0}^{(k)})^{\intercal}-\mathds{E}\left[x_{0}^{(k)}(x_{0}^{(k)})^{\intercal}\right]\right)\right\|\leq\left\|x_{0}^{(k)}({x_{0}^{(k)}})^{\intercal}\right\|_{F}+\left\|\mathds{E}\big{\\{}x_{0}^{(k)}({x_{0}^{(k)}})^{\intercal}\big{\\}}\right\|_{F}\leq
2c_{X}^{2},$

and

$\displaystyle\|\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)}-\Sigma_{A}^{\prime}\|_{2}$
$\displaystyle\leq\|\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)}\|_{2}+\|\Sigma_{A}^{\prime}\|_{2}=\|\bar{A}_{t}^{(k)}\|_{2}^{2}+\|\Sigma_{A}^{\prime}\|_{2}\leq
c_{\bar{A}}^{2}+\|\Sigma_{A}^{\prime}\|_{2},$
$\displaystyle\|\bar{B}_{t}^{(k)}\otimes\bar{B}_{t}^{(k)}-\Sigma_{B}^{\prime}\|_{2}$
$\displaystyle\leq\|\bar{B}_{t}^{(k)}\otimes\bar{B}_{t}^{(k)}\|_{2}+\|\Sigma_{B}^{\prime}\|_{2}=\|\bar{B}_{t}^{(k)}\|_{2}^{2}+\|\Sigma_{B}^{\prime}\|_{2}\leq
c_{\bar{B}}^{2}+\|\Sigma_{B}^{\prime}\|_{2}.$

However, these bounds may not depend on those in Assumption 2. For example,
when $x_{0}^{(k)}$ is a nonzero constant, $c_{\mu}=0$ but $c_{X}$ is positive.

The boundedness of the states and state-deviations follows from Assumptions 1
and 2 according to the following statement.

######  Lemma 3.14.

Suppose that Assumptions 1 and 2 hold, then for all $k\in[n_{r}]$ and $0\leq
t\leq\ell$ we have that

$\displaystyle\big{\|}x_{t}^{(k)}\big{\|}\leq
c_{M},\quad\big{\|}e_{t}^{(k)}\big{\|}\leq
c_{N},\quad\big{\|}x_{t}^{(k)}(x_{t}^{(k)})^{\intercal}\big{\|}_{2}\leq
c_{M}^{2},\quad\big{\|}e_{t}^{(k)}(e_{t}^{(k)})^{\intercal}\big{\|}_{2}\leq
c_{N}^{2},$

and

$\displaystyle\left\|\operatorname{vec}\left(x_{t}^{(k)}(x_{t}^{(k)})^{\intercal}-\mathds{E}\left\\{x_{t}^{(k)}(x_{t}^{(k)})^{\intercal}\right\\}\right)\right\|\leq
c_{F},\quad\left\|\operatorname{vec}\left(x_{t}^{(k)}(u_{t}^{(k)})^{\intercal}-\mathds{E}\left\\{x_{t}^{(k)}(u_{t}^{(k)})^{\intercal}\right\\}\right)\right\|\leq
c_{W},$

where

$\displaystyle c_{M}$ $\displaystyle:=\max_{0\leq
t\leq\ell}\left\\{c_{A}^{t}c_{X}+\sum_{i=0}^{t-1}c_{A}^{i}c_{B}c_{U}\right\\},$
$\displaystyle c_{N}$ $\displaystyle:=\max_{0\leq
t\leq\ell}\left\\{\|A\|_{2}^{t}c_{\mu}+\sum_{i=0}^{t-1}\|A\|_{2}^{i}(\|B\|_{2}c_{\nu}+c_{\bar{A}}c_{M}+c_{\bar{B}}c_{U})\right\\},$
$\displaystyle c_{F}$ $\displaystyle:=\max_{0\leq
t\leq\ell}\left\\{(\|A\|_{2}^{2}+\|\Sigma_{A}^{\prime}\|_{2})^{t}c_{\Delta
X}+\sum_{i=0}^{t-1}(\|A\|_{2}^{2}+\|\Sigma_{A}^{\prime}\|_{2})^{i}(c_{FX}+c_{FU}+c_{FXU})\right\\},$
$\displaystyle c_{W}$ $\displaystyle:=c_{N}c_{U}+c_{M}c_{\nu},$

and

$\displaystyle c_{A}$ $\displaystyle:=\|A\|_{2}+c_{\bar{A}},$ $\displaystyle
c_{B}$ $\displaystyle:=\|B\|_{2}+c_{\bar{B}},$ $\displaystyle c_{FX}$
$\displaystyle:=(2\|A\|_{2}c_{\bar{A}}+c_{\Sigma_{A}^{\prime}})c_{M}^{2},$
$\displaystyle c_{FU}$
$\displaystyle:=3(\|B\|_{2}^{2}+\|\Sigma_{B}^{\prime}\|_{2})c_{U}c_{\nu}+(2\|B\|_{2}c_{\bar{B}}+c_{\Sigma_{B}^{\prime}})c_{U}^{2},$
$\displaystyle c_{FXU}$
$\displaystyle:=2\|A\|_{2}\|B\|_{2}c_{W}+(2\|A\|_{2}c_{\bar{B}}+2\|B\|_{2}c_{\bar{A}}+2c_{\bar{A}}c_{\bar{B}})c_{M}c_{U}.$

See Appendix D.

######  Remark 3.15.

The quantity $c_{M}$ can be interpreted as a bound on the radius from the
origin to the outer boundary of the set of reachable states from any valid
$x_{0}$ over $\ell$ time steps. If the system is not robustly stable in the
sense that $c_{A}>1$, then the limit as $\ell\to\infty$ of $c_{M}$ could be
infinite. However, since we consider only finite-length rollouts, $c_{M}$ is
finite regardless of the stability properties of the system.

Analogous interpretations follow for the quantity $c_{N}$ and the reachable
state-deviations. Notice that the constants $c_{N}$ grows with increasing
maximum initial state and input deviations $c_{\mu}$ and $c_{\nu}$, and
maximum noise magnitudes $c_{\bar{A}}$ and $c_{\bar{B}}$. Conversely, $c_{N}$
vanishes as those quantities become smaller, i.e. in the case that the initial
state $x_{0}$ is a fixed deterministic value, the inputs $u_{t}$ follow a
deterministic sequence, and there is no multiplicative noise. Likewise,
$c_{F}$ vanishes in such a scenario, so that $c_{\Delta
X}=c_{FX}=c_{FU}=c_{FXU}=0$.

The following theorems state finite-sample results for the estimates of
$[A~{}B]$ and $[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]$,
whose proofs are given in Appendices F and G, respectively.

######  Theorem 3.16.

Suppose that Assumptions 1 and 2 hold. Fix a failure probability
$\delta\in(0,1)$. It holds with probability at least $1-\delta$ that

$\displaystyle\big{\|}\big{[}\hat{A}~{}\hat{B}\big{]}-[A~{}B]\big{\|}_{2}\leq\mathcal{O}\left(\sqrt{\frac{\ell\log(\ell/\delta)}{n_{r}}}\right).$

######  Theorem 3.17.

Under the same condition of Theorem 3.16, with probability at least
$1-\delta$, it holds that

$\displaystyle\left\|\Big{[}\hat{\tilde{\Sigma}}^{\prime}_{A}~{}\hat{\tilde{\Sigma}}^{\prime}_{B}\Big{]}-\big{[}\tilde{\Sigma}^{\prime}_{A}~{}\tilde{\Sigma}^{\prime}_{B}\big{]}\right\|_{2}\leq\mathcal{O}\left(\sqrt{\frac{\ell\log(\ell/\delta)}{n_{r}}}\right).$

######  Remark 3.18.

In Theorems 3.16 and 3.17, high-probability upper bounds are given for the
estimates of $[A~{}B]$ and
$[\tilde{\Sigma}^{\prime}_{A}~{}\tilde{\Sigma}^{\prime}_{B}]$. It can be
observed that these bounds shrink as $\mathcal{O}\left(1/\sqrt{n_{r}}\right)$
with the number of rollouts, and converge to zero as the number of rollouts
grows to infinity, indicating the consistency of the estimators. Note that the
bounds are deterministic, although they depend on the failure probability
$\delta$. The theorems also indicate that the probability of the estimation
error exceeding an arbitrary positive constant decays exponentially fast with
the number of rollouts, which is illustrated in Section 4.1.

The $\mathcal{O}(\cdot)$ notation hides the coefficients of the error bounds,
and the polynomial and exponential factors of $n$ and $m$ in the logarithm
term. Their explicit forms are given in Appendices F and G, respectively. The
coefficient of the estimation error of $[A~{}B]$ increases with
$\|\mathbf{Y}\|_{2}$, $\|\mathbf{Z}\|_{2}$, and the bound of the system, but
decreases with the minimum eigenvalue of $\mathbf{Z}\mathbf{Z}^{\intercal}$.
Similarly, the coefficient of the estimation error of
$[\tilde{\Sigma}^{\prime}_{A}~{}\tilde{\Sigma}^{\prime}_{B}]$ decreases with
the minimum eigenvalue of $\mathbf{D}\mathbf{D}^{\intercal}$, but increases
with $\|\mathbf{C}\|_{2}$, $\|\mathbf{D}\|_{2}$, and the bound of the system.
It also increases with $\|A\|_{2}$, $\|B\|_{2}$, and quantities related to the
second-moment dynamic of system states, because of the dependence of
$[\hat{\tilde{\Sigma}}^{\prime}_{A}~{}\hat{\tilde{\Sigma}}^{\prime}_{B}]$ on
$[\hat{A}~{}\hat{B}]$. From definition, $\mathbf{Y}$, $\mathbf{Z}$,
$\mathbf{C}$, and $\mathbf{D}$ depend on system parameters and inputs, so
proper input design could reduce the estimation error. It remains for future
study how to design the moments of inputs so that the coefficients of the
bounds can achieve their smallest values, and how to obtain data-dependent
bounds, because the nominal system matrix is unknown.

In Di and Lamperski, (2021), the authors study identification of System (1)
from single-trajectory data, by developing error bounds for a least-squares
algorithm, but it is unclear under what conditions of System (1) these error
bounds converge to zero. In contrast, our analysis provides sufficient
conditions under which the error bounds for estimates given by Algorithm 1
vanish. The results show that a relatively small rollout length is enough to
guarantee consistency, but the current bounds imply that longer rollout length
$\ell$ may lead to worse performance, which seems to be contrary to the
intuition that longer trajectory provides more information. This could result
from the averaging step which eliminates some excitation. Future work will
consider how to use the data more efficiently.

##  4 Numerical Simulations

In this section we empirically validate the theoretical results for Algorithm
1, and compare its performance with the recursive least-squares algorithm
based on single-trajectory data (Chen and Guo,, 2012; Lai and Wei,, 1982).

###  4.1 Consistency and Finite-Sample Result

This subsection considers identification of the $2$-dimensional system
discussed in Example 3.2 with parameters

$\displaystyle A=\begin{bmatrix}1&0.2\\\
0&1\end{bmatrix},~{}B=\begin{bmatrix}0.8\\\
1\end{bmatrix},~{}\Sigma_{A}=\frac{1}{40}\begin{bmatrix}8&-2&0&0\\\
-2&16&2&0\\\ 0&2&2&0\\\
0&0&0&8\end{bmatrix},~{}\Sigma_{B}=\frac{1}{40}\begin{bmatrix}5&-2\\\
-2&20\end{bmatrix}.$

According to the reshaping operator $G$ defined in the notation section and
the discussion in Example 3.2, it holds that

$\displaystyle\tilde{\Sigma}_{A}^{\prime}=\frac{1}{40}\begin{bmatrix}8&0&2\\\
-2&2&0\\\
16&0&8\end{bmatrix},~{}\tilde{\Sigma}_{B}^{\prime}=\frac{1}{40}\begin{bmatrix}5&-2&20\end{bmatrix}^{\intercal}.$
(39)

A simulated experiment is conducted with rollout data of length $\ell=4$. For
$0\leq t\leq 3$, $\nu_{t}$ is generated independently from uniform
distribution $\mathcal{U}([0,1])$ and then fixed. Three types of inputs are
considered: Gaussian, uniform, and deterministic inputs. An identical sequence
of input covariances, independently generated from $1$-dimensional Wishart
distribution $W_{p}(0.1,1)$ and then fixed, is used in the former two cases.
For the case of deterministic inputs, the covariances are set to be zero
(i.e., $\bar{U}_{t}=0$). In this setting $\mathbf{D}\mathbf{D}^{\intercal}$
can be invertible because the second moment of the input at time $t$ satisfies
that $U_{t}=\bar{U}_{t}+\nu_{t}\nu_{t}^{\intercal}$, and the generation of
$\nu_{t}$ provides randomness. For each case, Algorithm 1 is run for $50$
times. The mean of estimation error in each case is shown in Fig. 1. It can be
seen that Algorithm 1 converges with convergence rate
$\mathcal{O}(1/\sqrt{n_{r}})$, and performs similarly under all three types of
inputs. The algorithm fluctuates when the number of rollouts is small, which
may result from the averaging step.

Fig. 2 provides the relative frequency of the normalized estimation errors,
$\|[\hat{A}~{}\hat{B}]-[A~{}B]\|_{2}/\|[A~{}B]\|_{2}$ and
$\|[\hat{\tilde{\Sigma}}_{A}^{\prime}~{}\hat{\tilde{\Sigma}}_{B}^{\prime}]$
$-[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]\|_{2}/$
$\|[\tilde{\Sigma}_{A}^{\prime}~{}\tilde{\Sigma}_{B}^{\prime}]\|_{2}$,
exceeding a given constant, under the uniform-input case. This result shows an
exponential decay of the frequency and validates the finite-sample results.
The relative frequency of $n_{r}$ rollouts is denoted by $p_{n_{r}}$.

From Remark 3.1 and (33), it follows that (39) defines an equivalent class of
covariance matrices that generates the same second-moment dynamic of system
states. In the current example, $\Sigma_{B}$ is unique, but the following
covariance matrix is equivalent to $\Sigma_{A}$,

$\displaystyle\Sigma_{A}(\alpha)=\frac{1}{40}\begin{bmatrix}8&-2&0&1+\alpha\\\
-2&16&1-\alpha&0\\\ 0&1-\alpha&2&0\\\ 1+\alpha&0&0&8\end{bmatrix},$

with $\alpha\in\mathds{R}$ such that $\Sigma_{A}(\alpha)\succeq 0$. Fig. 3
illustrates the dynamic (4), starting with the same initial condition
$\mu_{0}=\mathbf{0}_{2}$ and $X_{0}=\mathbf{0}_{4}$, and with the noise
covariance matrix given by $(\Sigma_{A},\Sigma_{B})$,
$(\Sigma_{A}(1),\Sigma_{B})$, and estimates from Algorithm 1, respectively.
The parameters of inputs ($\nu_{t}$ and $\bar{U}_{t}$) are the same as the
uniform-input case. Note that $\Sigma_{A}(-1)=\Sigma_{A}$, and
$\Sigma_{A}(1)\succ 0$. It can be observed that the dynamics defined by
$(\Sigma_{A},\Sigma_{B})$ and $(\Sigma_{A}(1),\Sigma_{B})$ are identical, and
the dynamic defined by the estimates from Algorithm 1 is close to the former.

Figure 1: Consistency of Algorithm 1. Figure 2: Finite-sample result of
Algorithm 1. Figure 3:  The second-moment dynamic of system states defined by
several noise covariance matrices.

It is assumed that there is no additive noise in System (1), but Algorithm 1
can also be applied to identifying linear systems with both multiplicative and
additive noise. If additive noise $w_{t}$, independent of the inputs and the
multiplicative noise, exists, then write the system as

$\displaystyle x_{t+1}$
$\displaystyle=(A+\bar{A}_{t})x_{t}+(B+\bar{B}_{t})u_{t}+w_{t}$
$\displaystyle=(A+\bar{A}_{t})x_{t}+\big{[}B+\bar{B}_{t}~{}w_{t}\big{]}\begin{bmatrix}u_{t}\\\
1\end{bmatrix}.$ (40)

In other words, $w_{t}$ can be considered as a part of multiplicative noise
corresponding to a constant input equal to one. Consider the above
$2$-dimensional system with Gaussian noise
$w_{t}\sim\mathcal{N}(\mathbf{0}_{2},\sigma^{2}I_{2})$ and previously designed
Gaussian inputs. Note that in this case $\ell=6$ is needed because the
dimension of inputs increases by one in (40), compared with the original
system. Fig. 4 shows the consistency of Algorithm 1 under the presence of
additive noise.

Figure 4: Consistency of Algorithm 1 under both multiplicative and additive
noise.

###  4.2 Performance Comparison

The recursive form of the ordinary least-squares (OLS), namely, the recursive
least-squares (RLS), is widely used in identification of dynamic systems (Chen
and Guo,, 2012; Lai and Wei,, 1982). It is possible to apply RLS to identify
System (1) if certain conditions hold. Note that from System (1), we have that

$\displaystyle x_{t+1}$
$\displaystyle=Ax_{t}+Bu_{t}+(\bar{A}_{t}x_{t}+\bar{B}_{t}u_{t})$
$\displaystyle=Ax_{t}+Bu_{t}+w^{(1)}_{t},$

where $w^{(1)}_{t}:=\bar{A}_{t}x_{t}+\bar{B}_{t}u_{t}$ is considered to be
noise. Under Assumption 1, $\\{w^{(1)}_{t},\mathcal{F}_{t}\\}$ is a martingale
difference sequence, i.e., $\mathds{E}\\{w^{(1)}_{t}|\mathcal{F}_{t-1}\\}=0$,
where $\mathcal{F}_{t}:=\sigma(\bar{A}_{k},\bar{B}_{k},u_{k},0\leq k\leq t)$.
A mild condition for $w^{(1)}_{t}$ to ensure convergence of RLS in literature
(Chen and Guo,, 2012; Lai and Wei,, 1982) is that
$\sup_{t}\mathds{E}\\{\|w^{(1)}_{t}\|^{\beta}|\mathcal{F}_{t-1}\\}<\infty$
holds a.s. for some $\beta>2$. However in our case $w^{(1)}_{t}$ is state-
dependent, so certain stability assumption is needed to ensure this
boundedness condition. This fact means that RLS could fail if the nominal part
of System (1) is marginally stable ($\rho(A)=1$) or unstable ($\rho(A)>1$). In
contrast, Algorithm 1 can handle this situation with the help of multiple-
trajectory data. Similarly, the noise covariance matrix of System (1) may be
estimated using the following dynamic

$\displaystyle P_{1}\operatorname{vec}(x_{t+1}x_{t+1}^{\intercal})$
$\displaystyle=P_{1}\big{(}(A+\bar{A}_{t})\otimes(A+\bar{A}_{t})\big{)}Q_{1}P_{1}\operatorname{vec}(x_{t}x_{t}^{\intercal})+P_{1}\big{(}(B+\bar{B}_{t})\otimes(B+\bar{B}_{t})\big{)}Q_{2}P_{2}\operatorname{vec}(u_{t}u_{t}^{\intercal})$
$\displaystyle\quad+P_{1}\big{(}(B+\bar{B}_{t})\otimes(A+\bar{A}_{t})\big{)}\operatorname{vec}(x_{t}u_{t}^{\intercal})+P_{1}\big{(}(A+\bar{A}_{t})\otimes(B+\bar{B}_{t})\big{)}\operatorname{vec}(u_{t}x_{t}^{\intercal})$
$\displaystyle=P_{1}\mathds{E}\big{\\{}(A+\bar{A}_{t})\otimes(A+\bar{A}_{t})\big{\\}}Q_{1}P_{1}\operatorname{vec}(x_{t}x_{t}^{\intercal})+P_{1}\mathds{E}\big{\\{}(B+\bar{B}_{t})\otimes(B+\bar{B}_{t})\big{\\}}Q_{2}P_{2}\operatorname{vec}(u_{t}u_{t}^{\intercal})$
$\displaystyle\quad+P_{1}(B\otimes
A)\operatorname{vec}(x_{t}u_{t}^{\intercal})+P_{1}(A\otimes
B)\operatorname{vec}(u_{t}x_{t}^{\intercal})+w^{(2)}_{t},$

where

$\displaystyle w^{(2)}_{t}$
$\displaystyle:=P_{1}\big{(}(A+\bar{A}_{t})\otimes(A+\bar{A}_{t})-\mathds{E}\\{(A+\bar{A}_{t})\otimes(A+\bar{A}_{t})\\}\big{)}Q_{1}P_{1}\operatorname{vec}(x_{t}x_{t}^{\intercal})$
$\displaystyle\quad+P_{1}\big{(}(B+\bar{B}_{t})\otimes(B+\bar{B}_{t})-\mathds{E}\\{(B+\bar{B}_{t})\otimes(B+\bar{B}_{t})\\}\big{)}Q_{2}P_{2}\operatorname{vec}(u_{t}u_{t}^{\intercal})$
$\displaystyle\quad+P_{1}(B\otimes\bar{A}_{t}+\bar{B}_{t}\otimes
A+\bar{B}_{t}\otimes\bar{A}_{t})\operatorname{vec}(x_{t}u_{t}^{\intercal})$
$\displaystyle\quad+P_{1}(A\otimes\bar{B}_{t}+\bar{A}_{t}\otimes
B+\bar{A}_{t}\otimes\bar{B}_{t})\operatorname{vec}(u_{t}x_{t}^{\intercal}).$

It can be verified that, under Assumption 1, despite state-dependent,
$\\{w^{(2)}_{t},\mathcal{F}_{t}\\}$ is also a martingale difference sequence.
To estimate the covariance matrix of the multiplicative noise, Di and
Lamperski, (2021) apply OLS, which is equivalent to RLS. Note that, when using
OLS or RLS, one estimates the second moments of $A+\bar{A}_{t}$ and
$B+\bar{B}_{t}$, rather than their covariance matrices, which are $\Sigma_{A}$
and $\Sigma_{B}$ in our context. The estimation of noise covariance is still
coupled with the estimation of the nominal system, since
$\Sigma_{A}^{\prime}=\mathds{E}\\{(A+\bar{A}_{t})\otimes(A+\bar{A}_{t})\\}-A\otimes
A$ and
$\Sigma_{B}^{\prime}=\mathds{E}\\{(B+\bar{B}_{t})\otimes(B+\bar{B}_{t})\\}-B\otimes
B$.

Figure 5: Performance comparison of RLS, RLSp, and Algorithm 1.

To compare the performance of RLS and Algorithm 1, we consider four systems.
In the first case, the nominal system matrices are

$\displaystyle A=\begin{bmatrix}0.6&0.2\\\
0&0.6\end{bmatrix},~{}B=\begin{bmatrix}0.8\\\ 1\end{bmatrix},$

and both $\Sigma_{A}$ and $\Sigma_{B}$ are zero matrices. That is, a linear
system without noise and $\rho(A)=0.6$, where $\rho(A)$ is the spectral radius
of $A$. We use this case to show the consistency of RLS. In the other three
cases, the matrix $A$ is set to be

$\begin{bmatrix}0.6&0.2\\\ 0&0.6\end{bmatrix},~{}\begin{bmatrix}0.8&0.2\\\
0&0.8\end{bmatrix},\text{ and }\begin{bmatrix}1&0.2\\\ 0&1\end{bmatrix},$

respectively. $B$ is the same as the first case, while $\Sigma_{A}$ and
$\Sigma_{B}$ in Section 4.1 are adopted to be the covariance matrices. The
implementation of Algorithm 1 is the same as in Section 4.1. That is,
$\nu_{t}$ and $\bar{U}_{t}$ are randomly generated, and then fixed in all runs
of the entire numerical experiment. The input $u_{t}$ at time $0\leq
t\leq\ell-1$ in each rollout is generated from Gaussian distribution
$\mathcal{N}(\nu_{t},\bar{U}_{t})$, and $\ell=4$. Since RLS is based on
single-trajectory data, the length of the trajectory is set to be $\ell
n_{r}$, so that the number of samples that RLS uses is the same as that of
Algorithm 1. RLS with independent standard Gaussian inputs is considered as a
baseline. In order to rule out the effect of different input design, we also
run RLS with periodic inputs (RLSp) satisfying that, in each period, the
inputs are generated in the same way as those in a rollout of Algorithm 1.

For each system, the three algorithms, RLS, RLSp, and Algorithm 1 are run for
$50$ times, respectively. The mean of estimation error in each case is
presented in Fig. 5. It can be observed that RLS and RLSp perform similarly in
all cases. When multiplicative noise is absent, they converge slightly faster
than Algorithm 1. They are also a little better than Algorithm 1, in the case
$\rho(A)=0.6$ with noise, for the estimation of $[A~{}B]$, indicating OLS
could be applied to Algorithm 1 as a way to estimate $[A~{}B]$. However,
Algorithm 1 surpasses RLS and RLSp when identifying the noise covariance
matrix. Moreover, the performance of RLS gets worse as $\rho(A)$ grows.
Interestingly, in the case of $\rho(A)=0.8$, although the nominal system is
stable, the second-moment dynamic of system states is not. This instability
leads to degraded performance of RLS estimating $[A~{}B]$ and divergence of
RLS estimating the covariance matrix. In the marginally stable case, namely
$\rho(A)=1$, RLS and RLSp explode in finite time. In contrast, Algorithm 1
behaves almost identically for all cases (the consistency of Algorithm 1 in
the marginally stable case is shown in Fig. 1). To sum up, Algorithm 1 can
deal with the estimation of noise covariance matrix better and relies less on
the stability of both the nominal system and the second-moment dynamic of
system states.

##  5 Conclusion and Future Work

In this paper an identification algorithm based on multiple-trajectory data
was proposed for linear systems with multiplicative noise. With appropriately
designed exciting inputs, the proposed algorithm is able to jointly estimate
the nominal system and the multiplicative noise covariance. The asymptotic and
non-asymptotic performance of the algorithm was analyzed theoretically, and
illustrated by numerical experiments. Future work include studying more
efficient algorithms that can be used in online settings, optimal and adaptive
input design, sparsity-promoting regularization for identification of
networked systems, end-to-end finite-sample performance guarantees for
identification-based optimal control, and applications to identification of
cyber-physical systems with coupling between noise and inputs.

Appendix

##  Appendix A Proof of Proposition 3

We begin with a standard result from real analysis (Caron and Traynor,, 2005;
Federer,, 2014) regarding the zero set of a polynomial.

######  Lemma A.19.

A polynomial function $\mathds{R}^{n}$ to $\mathds{R}$ is either identically
$0$ or non-zero almost everywhere.

It suffices to consider the case with $\ell=n+m$, implying $\mathbf{Z}$ is a
square matrix. Note that when $(A,B)$ and $\mu_{0}$ are fixed, $|\mathbf{Z}|$
is a polynomial of
$\nu:=[\nu_{0}^{\intercal}~{}\cdots~{}\nu_{\ell-1}^{\intercal}]^{\intercal}\in\mathds{R}^{m\ell}$.
Hence the existence of a vector in $\mathds{R}^{m\ell}$ such that
$|\mathbf{Z}|\not=0$ implies that $|\mathbf{Z}|\not=0$ almost everywhere from
Lemma A.19.

First we verify the conclusion for $\mu_{0}=\mathbf{0}_{n}$. It follows from
the definition of controllability and the assumption that
$[A^{n-1}B~{}A^{n-2}B~{}\cdots~{}B]$ has full row rank. Without loss of
generality, let $B_{1}$, $AB_{1}$, $\dots$, $A^{r_{1}}B_{1}$, $B_{2}$,
$\dots$, $A^{r_{2}}B_{2}$, $\dots$, $B_{p}$, $\dots$, $A^{r_{p}}B_{p}$ be a
basis of $\mathds{R}^{n}$, where $B_{i}$ is the $i$-th column of $B$, $1\leq
p\leq m$, $0\leq r_{i}\leq n-1$, $1\leq i\leq p$, and
$p+\sum_{i=1}^{p}r_{i}=n$. Moreover, $A^{k}B_{i}$, $k>r_{i}$ and $1\leq i\leq
p$, can be written as a linear combination of $B_{1}$, $\dots$,
$A^{r_{1}}B_{1}$, $\cdots$, $A^{r_{i}}B_{i}$.

Let
$\tilde{\nu}:=[\tilde{\nu}_{0}^{\intercal}~{}\cdots~{}\tilde{\nu}_{n-1}^{\intercal}]^{\intercal}$
be such that $\tilde{\nu}_{0}=\mathbf{e}_{1}^{m}$,
$\tilde{\nu}_{q_{1}}=\mathbf{e}_{2}^{m}$, $\dots$,
$\tilde{\nu}_{q_{k}}=\mathbf{e}_{k+1}^{m}$, $\dots$,
$\tilde{\nu}_{q_{p-1}}=\mathbf{e}_{p}^{m}$ (or any nonzero multiplier of
respective unit vectors), where $q_{k}=k+\sum_{i=1}^{k}r_{i}$, $1\leq k\leq
p$, and $\tilde{\nu}_{i}=\mathbf{0}_{m}$ for other $1\leq i\leq n-1$. Then for
$\tilde{\mu}_{t+1}=A\tilde{\mu}_{t}+B\tilde{\nu}_{t}$, $0\leq t\leq n-1$, and
$\tilde{\mu}_{0}=\mathbf{0}_{n}$, it holds that

$\displaystyle\tilde{\mu}_{t}$ $\displaystyle=A^{t-1}B_{1},$
$\displaystyle~{}1\leq t\leq q_{1},$ $\displaystyle\tilde{\mu}_{t}$
$\displaystyle=A^{t-q_{1}-1}B_{2}+A^{t-1}B_{1},$ $\displaystyle~{}q_{1}+1\leq
t\leq q_{1}+q_{2},$ $\displaystyle\vdots$ $\displaystyle\tilde{\mu}_{t}$
$\displaystyle=A^{t-q_{p-1}-1}B_{p}+\cdots+A^{t-q_{1}-1}B_{2}+A^{t-1}B_{1},$
$\displaystyle~{}q_{p-1}+1\leq t\leq n.$

The definitions of $B_{1}$, $\dots$, $A^{r_{p}}B_{p}$, and $q_{k}$ imply that
$\tilde{\mu}_{1}$, $\dots$, $\tilde{\mu}_{n}$ are linearly independent. Hence
it is shown that there exists a vector
$[\nu_{0}^{\intercal}~{}\cdots~{}\nu_{n-1}^{\intercal}]^{\intercal}$ such that
$[\mu_{1}~{}\cdots~{}\mu_{n}]$ has full rank. If $m=1$, let
$\nu=[\tilde{\nu}^{\intercal}~{}0]^{\intercal}$ and

$\displaystyle\mathbf{Z}=\begin{bmatrix}\tilde{\mu}_{n}&\cdots&\tilde{\mu}_{1}&\mathbf{0}_{n}\\\
0&\cdots&0&1\end{bmatrix}$

has full rank. In the case of $m\geq 2$, fix $\nu_{n}$ to be zero. Set
$\nu_{n+1}=c_{1}\mathbf{e}_{2}^{m}$, $c_{1}\in\mathds{R}$, and

$\displaystyle\left|\begin{array}[]{cccccccc}\mu_{n+1}&\tilde{\mu}_{n}&\cdots&\cdots&\tilde{\mu}_{q_{1}}&\cdots&\cdots&\mathbf{0}_{n}\\\
0&0&\cdots&\cdots&0&\cdots&0&1\\\
c_{1}&0&\cdots&0&d_{1}&0&\cdots&0\end{array}\right|=c_{1}\left|\begin{array}[]{cccc}\tilde{\mu}_{n}&\cdots&\tilde{\mu}_{1}&\mathbf{0}_{n}\\\
0&\cdots&0&1\end{array}\right|+(-1)^{2n+4-q_{1}}d_{1}\left|\begin{array}[]{cccccccc}\mu_{n+1}&\tilde{\mu}_{n}&\cdots&\tilde{\mu}_{q_{1}+1}&\tilde{\mu}_{q_{1}-1}&\cdots&\mathbf{0}_{n}\\\
0&0&\cdots&\cdots&\cdots&0&1\end{array}\right|,$

where $d_{1}=1$ if $p>1$ and $d_{1}=0$ if $p=1$. Hence there must exist
$\tilde{c}_{1}\in\mathds{R}$ such that the above determinant is nonzero.
Inductively, set $\nu_{n+k}=c_{k}\mathbf{e}_{k+1}^{m}$, $2\leq k\leq m-1$, and
it is able to find $\tilde{c}_{k}$, $2\leq k\leq m-1$, (consequently
$\tilde{\nu}_{0}$, $\dots$, $\tilde{\nu}_{n+m-1}$) such that
$|\mathbf{Z}|\not=0$. Now suppose that the system starts with
$\mu_{0}\not=\mathbf{0}_{n}$. Since $[A^{n-1}B~{}A^{n-2}B~{}\cdots~{}B]$ has
full row rank, so does the matrix $[A^{n-1}B$ $A^{n-2}B$
$~{}\cdots~{}B~{}A^{n-1}\mu_{0}~{}\cdots~{}\mu_{0}]$. Thus, without loss of
generality, assume that $\mu_{0}$, $A\mu_{0}$, $\dots$, $A^{r_{0}}\mu_{0}$,
$B_{1}$, $\dots$, $A^{r_{1}}B_{1}$, $\dots$, $B_{p}$, $\dots$,
$A^{r_{p}}B_{p}$ is a basis of $\mathds{R}^{n}$, where $0\leq p\leq m$ ($p=0$
means that there is no $B_{i}$), $0\leq r_{i}\leq n-1$, $0\leq i\leq p$, and
$p+1+\sum_{i=0}^{p}r_{i}=n$. Moreover, $r_{i}\leq n-2$ for all $1\leq i\leq
p$. It can be verified that there exists a vector
$[\nu_{0}^{\intercal}~{}\cdots~{}\nu_{n-1}^{\intercal}]^{\intercal}$ such that
$[\mu_{0}~{}\cdots~{}\mu_{n-1}]$ has full rank. Therefore in a similar way the
conclusion can be obtained.

######  Remark A.20.

From the proof we know that even if $(A,B)$ is not controllable, $\mathbf{Z}$
can still have full row rank as long as $[A^{n-1}B$ $A^{n-2}B$
$~{}\cdots~{}B~{}A^{n-1}\mu_{0}~{}\cdots~{}\mu_{0}]$ has full row rank for
some $\mu_{0}\in\mathds{R}^{n}$.

##  Appendix B Proof of Proposition 4

Write (5) as

$\displaystyle\tilde{X}_{t+1}$
$\displaystyle=(\tilde{A}+\tilde{\Sigma}_{A}^{\prime})\tilde{X}_{t}+(\tilde{B}+\tilde{\Sigma}_{B}^{\prime})\tilde{U}_{t}+(K_{BA}W_{t}+K_{AB}W_{t}^{\prime})$
$\displaystyle:=\breve{A}\tilde{X}_{t}+\breve{B}\tilde{U}_{t}+\eta_{t}.$

Note that after setting $\nu_{t}=0$ for all $0\leq t\leq\ell-1$,
$\tilde{U}_{t}=\operatorname{svec}(\bar{U}_{t})$ and the above dynamic becomes

$\displaystyle\tilde{X}_{t+1}=\breve{A}\tilde{X}_{t}+\breve{B}\tilde{U}_{t}.$

Similar to the proof of Proposition 3, first examine if $\tilde{X}_{0}$,
$\breve{A}\tilde{X}_{0}$, $\dots$, $\breve{A}^{n(n+1)/2-1}\tilde{X}_{0}$ are
linearly independent. If not, select some columns of
$[\breve{A}^{n(n+1)/2-1}\breve{B},\dots,\breve{B}]$ to together form a basis
of $\mathds{R}^{n(n+1)/2}$. Considering $\tilde{U}_{t}$ as an input makes the
rest of the proof essentially the same as the proof of Proposition 3.

##  Appendix C Proof of Theorem 3.10

Consider each rollout
$[(x_{0}^{(k)})^{\intercal},\dots,(x_{\ell}^{(k)})^{\intercal}]^{\intercal}$
as an independent sample of the random vector
$\mathbf{x}_{\ell}:=[x_{0}^{\intercal},\dots,x_{\ell}^{\intercal}]^{\intercal}$,
and from Assumption 1 (ii) and (iii) we know that the random vector
$\mathbf{x}_{\ell}$ has finite first and second moments. So it follows from
the Kolmogorov’s strong law of large numbers that
$\hat{\mathbf{Y}}\to\mathbf{Y}$ a.s., and similarly
$\hat{\mathbf{Z}}\to\mathbf{Z}$ a.s., as $n_{r}\to\infty$. Hence
$\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}\to\mathbf{Y}\mathbf{Z}^{\intercal}$
and
$\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}\to\mathbf{Z}\mathbf{Z}^{\intercal}$
a.s. From the assumption that $\mathbf{Z}\mathbf{Z}^{\intercal}$ is invertible
and the continuous mapping theorem (Theorem 2.3 of (Van der Vaart,, 2000)), it
can be obtained that as $n_{r}\to\infty$

$\displaystyle(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}\to(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1},\text{
a.s.}$

When $(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}$ does not exist, in
Algorithm 1 we replace it by
$(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger}$. Thus
$(\hat{A},\hat{B})\to(A,B)$. Combining the above convergence with the
Kolmogorov’s strong law of large numbers, the convergence of
$\hat{\mathbf{C}}$ and $\hat{\mathbf{D}}$ follows. Therefore, applying the
continuous mapping theorem again, we obtain the consistency of the estimator
$(\hat{\tilde{\Sigma}}_{A}^{\prime},\hat{\tilde{\Sigma}}_{B}^{\prime})$.

##  Appendix D Proof of Lemma 3.14

For the first claim, regarding states, taking the norm of both sides of System
(1), at time step $t$ we have

$\displaystyle\big{\|}x_{t+1}^{(k)}\big{\|}$
$\displaystyle=\big{\|}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}\big{\|}.$
(41)

Using the triangle inequality and the fact that the spectral norm is
compatible with the Euclidean norm, we have

$\displaystyle\big{\|}x_{t+1}^{(k)}\big{\|}$
$\displaystyle\leq\big{\|}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}\big{\|}+\big{\|}(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}\big{\|}$
$\displaystyle\leq\big{\|}A+\bar{A}_{t}^{(k)}\big{\|}_{2}\big{\|}x_{t}^{(k)}\big{\|}+\big{\|}B+\bar{B}_{t}^{(k)}\big{\|}_{2}\big{\|}u_{t}^{(k)}\big{\|}.$

Using the triangle inequality and Assumption 2 (iii) we have

$\displaystyle\big{\|}A+\bar{A}_{t}^{(k)}\big{\|}_{2}\leq\max_{0\leq
t\leq\ell}\big{\|}A+\bar{A}_{t}^{(k)}\big{\|}_{2}\leq\|A\|_{2}+\max_{0\leq
t\leq\ell}\big{\|}\bar{A}_{t}^{(k)}\big{\|}_{2}$ $\displaystyle=:c_{A},$
$\displaystyle\big{\|}B+\bar{B}_{t}^{(k)}\big{\|}_{2}\leq\max_{0\leq
t\leq\ell}\big{\|}B+\bar{B}_{t}^{(k)}\big{\|}_{2}\leq\|B\|_{2}+\max_{0\leq
t\leq\ell}\big{\|}\bar{B}_{t}^{(k)}\big{\|}_{2}$ $\displaystyle=:c_{B},$

so from Assumption 2 (ii)

$\displaystyle\big{\|}x_{t+1}^{(k)}\big{\|}$ $\displaystyle\leq
c_{A}\big{\|}x_{t}^{(k)}\big{\|}+c_{B}c_{U}.$ (42)

For the base case when $t=0$, we have by Assumption 2 (i) that
$\big{\|}x_{0}^{(k)}\big{\|}\leq c_{X}$. Applying (42) inductively with the
base case proves the first claim.

For the second claim, regarding the state-deviations, we have

$\displaystyle e_{t+1}^{(k)}$
$\displaystyle=x_{t+1}^{(k)}-\mathds{E}\big{\\{}x_{t+1}^{(k)}\big{\\}}$
$\displaystyle=\big{(}A+\bar{A}_{t}^{(k)}\big{)}x_{t}^{(k)}+\big{(}B+\bar{B}_{t}^{(k)}\big{)}u_{t}^{(k)}-\mathds{E}\big{\\{}\big{(}A+\bar{A}_{t}^{(k)}\big{)}x_{t}^{(k)}+\big{(}B+\bar{B}_{t}^{(k)}\big{)}u_{t}^{(k)}\big{\\}}$
$\displaystyle=\big{(}A+\bar{A}_{t}^{(k)}\big{)}x_{t}^{(k)}+\big{(}B+\bar{B}_{t}^{(k)}\big{)}u_{t}^{(k)}-A\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}-B\mathds{E}\big{\\{}u_{t}^{(k)}\big{\\}}$
$\displaystyle=A\big{(}x_{t}^{(k)}-\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\big{)}+B\big{(}u_{t}^{(k)}-\mathds{E}\big{\\{}u_{t}^{(k)}\big{\\}}\big{)}+\bar{A}_{t}^{(k)}x_{t}^{(k)}+\bar{B}_{t}^{(k)}u_{t}^{(k)}$
$\displaystyle=Ae_{t}^{(k)}+Bd_{t}^{(k)}+\bar{A}_{t}^{(k)}x_{t}^{(k)}+\bar{B}_{t}^{(k)}u_{t}^{(k)}.$

Taking the norm of both sides, using submultiplicativity and triangle
inequality, we have

$\displaystyle\|e_{t+1}^{(k)}\|$
$\displaystyle=\|Ae_{t}^{(k)}+Bd_{t}^{(k)}+\bar{A}_{t}^{(k)}x_{t}^{(k)}+\bar{B}_{t}^{(k)}u_{t}^{(k)}\|$
$\displaystyle\leq\|A\|_{2}\|e_{t}^{(k)}\|+\|B\|_{2}\|d_{t}^{(k)}\|+\|\bar{A}_{t}^{(k)}\|_{2}\|x_{t}^{(k)}\|+\|\bar{B}_{t}^{(k)}\|_{2}\|u_{t}^{(k)}\|$
$\displaystyle\leq\|A\|_{2}\|e_{t}^{(k)}\|+\|B\|_{2}c_{\nu}+c_{\bar{A}}c_{M}+c_{\bar{B}}c_{U},$
(43)

where the final inequality follows from the first part of Lemma 3.14 and
Assumptions 2 (ii) and (iii). For the base case when $t=0$, we have by
Assumption 2 (i) that $\big{\|}e_{0}^{(k)}\big{\|}\leq c_{\mu}$. Applying (43)
inductively with the base case proves the second claim.

By the definition of the spectral norm and the first and second claims we have

$\displaystyle\big{\|}x_{t}^{(k)}{x_{t}^{(k)}}^{\intercal}\big{\|}_{2}$
$\displaystyle=\big{\|}x_{t}^{(k)}\big{\|}_{2}^{2}\leq c_{M}^{2},$
$\displaystyle\big{\|}e_{t}^{(k)}{e_{t}^{(k)}}^{\intercal}\big{\|}_{2}$
$\displaystyle=\big{\|}e_{t}^{(k)}\big{\|}_{2}^{2}\leq c_{N}^{2},$

proving the third and fourth claims.

For the fifth and sixth claims, define the quantities

$\displaystyle\Delta
X_{t}:=\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\big{\\}}\right)=\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\right)-X_{t},$
$\displaystyle\Delta
U_{t}:=\operatorname{vec}\left(u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\big{\\}}\right)=\operatorname{vec}\left(u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right)-U_{t},$
$\displaystyle\Delta
W_{t}:=\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\big{\\}}\right)=\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right)-W_{t}.$

We can bound $\|\Delta U_{t}\|$ as

$\displaystyle\|\Delta U_{t}\|$
$\displaystyle=\left\|\operatorname{vec}\left(u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\big{\\}}\right)\right\|$
$\displaystyle=\left\|\operatorname{vec}\left(u_{t}^{(k)}\big{(}{u_{t}^{(k)}}-\mathds{E}\\{u_{t}^{(k)}\\}\big{)}^{\intercal}+\big{(}u_{t}^{(k)}-\mathds{E}\\{u_{t}^{(k)}\\}\big{)}\mathds{E}\\{u_{t}^{(k)}\\}^{\intercal}-\mathds{E}\big{\\{}u_{t}^{(k)}\big{(}{u_{t}^{(k)}}-\mathds{E}{u_{t}^{(k)}}\big{)}^{\intercal}\big{\\}}\right)\right\|$
$\displaystyle=\left\|u_{t}^{(k)}\big{(}{u_{t}^{(k)}}-\mathds{E}\\{u_{t}^{(k)}\\}\big{)}^{\intercal}+\big{(}u_{t}^{(k)}-\mathds{E}\\{u_{t}^{(k)}\\}\big{)}\mathds{E}\\{u_{t}^{(k)}\\}^{\intercal}-\mathds{E}\big{\\{}u_{t}^{(k)}\big{(}{u_{t}^{(k)}}-\mathds{E}{u_{t}^{(k)}}\big{)}^{\intercal}\big{\\}}\right\|_{F}$
$\displaystyle\leq\left\|u_{t}^{(k)}\big{(}{u_{t}^{(k)}}-\mathds{E}\\{u_{t}^{(k)}\\}\big{)}^{\intercal}\right\|_{F}+\left\|\big{(}u_{t}^{(k)}-\mathds{E}\\{u_{t}^{(k)}\\}\big{)}\mathds{E}\\{u_{t}^{(k)}\\}^{\intercal}\right\|_{F}+\mathds{E}\left\\{\left\|u_{t}^{(k)}\big{(}{u_{t}^{(k)}}-\mathds{E}{u_{t}^{(k)}}\big{)}^{\intercal}\right\|_{F}\right\\}$
$\displaystyle\leq 3c_{U}c_{\nu}.$

For the sixth claim, we can bound $\|\Delta W_{t}\|$ as

$\displaystyle\|\Delta W_{t}\|$
$\displaystyle=\left\|\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\big{\\}}\right)\right\|$
$\displaystyle=\left\|\operatorname{vec}\left(\big{(}x_{t}^{(k)}-\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\big{)}({u_{t}^{(k)}})^{\intercal}+\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\big{(}u_{t}^{(k)}-\mathds{E}\big{\\{}u_{t}^{(k)}\big{\\}}\big{)}^{\intercal}\right)\right\|$
$\displaystyle=\Big{\|}\big{(}x_{t}^{(k)}-\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\big{)}({u_{t}^{(k)}})^{\intercal}+\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\big{(}u_{t}^{(k)}-\mathds{E}\big{\\{}u_{t}^{(k)}\big{\\}}\big{)}^{\intercal}\Big{\|}_{F}$
$\displaystyle\leq\big{\|}x_{t}^{(k)}-\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\big{\|}\big{\|}{u_{t}^{(k)}}\big{\|}+\big{\|}\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\big{\|}\big{\|}u_{t}^{(k)}-\mathds{E}\big{\\{}u_{t}^{(k)}\big{\\}}\big{\|}$
$\displaystyle\leq c_{N}c_{U}+c_{M}c_{\nu}=:c_{W}.$

For the fifth claim, substituting the dynamics and expanding the products we
have

$\displaystyle\Delta X_{t+1}$
$\displaystyle=\operatorname{vec}\left(x_{t+1}^{(k)}({x_{t+1}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t+1}^{(k)}({x_{t+1}^{(k)}})^{\intercal}\big{\\}}\right)$
$\displaystyle=\operatorname{vec}\left(\left[(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}\right]\left[(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}\right]^{\intercal}\right.$
$\displaystyle\quad\left.-\mathds{E}\left\\{\left[(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}\right]\left[(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}\right]^{\intercal}\right\\}\right)$
$\displaystyle=\operatorname{vec}\bigg{(}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}+(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}$
$\displaystyle\quad+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}$
$\displaystyle-\mathds{E}\big{\\{}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}+(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}$
$\displaystyle\quad+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}+(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}\big{\\}}\bigg{)}$
$\displaystyle=\operatorname{vec}\bigg{(}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}-\mathds{E}\big{\\{}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}\big{\\}}$
$\displaystyle\quad(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}-\mathds{E}\big{\\{}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}\big{\\}}$
$\displaystyle\quad(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}-\mathds{E}\big{\\{}(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}\big{\\}}$
$\displaystyle\quad(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}-\mathds{E}\big{\\{}(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}\big{\\}}\bigg{)}.$
(44)

Considering the first of the four terms of (44), we have

$\displaystyle\boxed{1}$
$\displaystyle:=\operatorname{vec}\left((A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}-\mathds{E}\big{\\{}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(A+\bar{A}_{t}^{(k)})^{\intercal}\big{\\}}\right)$
$\displaystyle=\operatorname{vec}\bigg{(}Ax_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}A^{\intercal}+Ax_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(\bar{A}_{t}^{(k)})^{\intercal}+\bar{A}_{t}^{(k)}x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}A^{\intercal}+\bar{A}_{t}^{(k)}x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(\bar{A}_{t}^{(k)})^{\intercal}$
$\displaystyle-A\mathds{E}\big{\\{}x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\big{\\}}A^{\intercal}-0-0-\mathds{E}\big{\\{}\bar{A}_{t}^{(k)}x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}(\bar{A}_{t}^{(k)})^{\intercal}\big{\\}}\bigg{)}$
$\displaystyle=(A\otimes
A)\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\big{\\}}\right)$
$\displaystyle\quad+(A\otimes\bar{A}_{t}^{(k)})\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\right)+(\bar{A}_{t}^{(k)}\otimes
A)\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\right)$
$\displaystyle\quad+(\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\right)-\mathds{E}\big{\\{}(\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})\big{\\}}\mathds{E}\big{\\{}\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\right)\big{\\}}$
$\displaystyle=(A\otimes
A)\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\big{\\}}\right)$
$\displaystyle\quad+(A\otimes\bar{A}_{t}^{(k)})\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\right)+(\bar{A}_{t}^{(k)}\otimes
A)\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\right)$
$\displaystyle\quad+\left[(\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})-\mathds{E}\big{\\{}(\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})\big{\\}}\right]\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\right)$
$\displaystyle\quad+\mathds{E}\big{\\{}(\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})\big{\\}}\operatorname{vec}\left(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal}\big{\\}}\right).$

Taking norms, and substituting notated quantities, we have

$\displaystyle\left\|\boxed{1}\right\|$ $\displaystyle\leq\|A\otimes
A\|_{2}\|\Delta
X_{t}\|+\|A\otimes\bar{A}_{t}^{(k)}\|_{2}\|\operatorname{vec}(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal})\|+\|\bar{A}_{t}^{(k)}\otimes
A\|_{2}\|\operatorname{vec}(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal})\|$
$\displaystyle\quad+\|(\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})-\Sigma_{A}^{\prime}\|_{2}\|\operatorname{vec}(x_{t}^{(k)}({x_{t}^{(k)}})^{\intercal})\|+\|\Sigma_{A}^{\prime}\|_{2}\|\Delta
X_{t}\|$ $\displaystyle\leq(\|A\|_{2}^{2}+\|\Sigma_{A}^{\prime}\|_{2})\|\Delta
X_{t}\|+\left(2\|A\|_{2}\|\bar{A}_{t}^{(k)}\|_{2}+\|(\bar{A}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})-\Sigma_{A}^{\prime}\|_{2}\right)\|x_{t}^{(k)}\|^{2}$
$\displaystyle\leq(\|A\|_{2}^{2}+\|\Sigma_{A}^{\prime}\|_{2})\|\Delta
X_{t}\|+(2\|A\|_{2}c_{\bar{A}}+c_{\Sigma_{A}^{\prime}})c_{M}^{2}.$

Applying identical arguments to the fourth term of (44)

$\displaystyle\boxed{4}$
$\displaystyle:=\operatorname{vec}\left((B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}-\mathds{E}\big{\\{}(B+\bar{B}_{t}^{(k)})u_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}\big{\\}}\right),$

we obtain the norm bound

$\displaystyle\left\|\boxed{4}\right\|$
$\displaystyle\leq(\|B\|_{2}^{2}+\|\Sigma_{B}^{\prime}\|_{2})\|\Delta
U_{t}\|+(2\|B\|_{2}c_{\bar{B}}+c_{\Sigma_{B}^{\prime}})c_{U}^{2}$
$\displaystyle\leq
3(\|B\|_{2}^{2}+\|\Sigma_{B}^{\prime}\|_{2})c_{U}c_{\nu}+(2\|B\|_{2}c_{\bar{B}}+c_{\Sigma_{B}^{\prime}})c_{U}^{2}.$

Likewise, for the second term of (44), we have

$\displaystyle\boxed{2}$
$\displaystyle:=\operatorname{vec}\left((A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}-\mathds{E}\big{\\{}(A+\bar{A}_{t}^{(k)})x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(B+\bar{B}_{t}^{(k)})^{\intercal}\big{\\}}\right)$
$\displaystyle=\operatorname{vec}\bigg{(}Ax_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}B^{\intercal}+Ax_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(\bar{B}_{t}^{(k)})^{\intercal}+\bar{A}_{t}^{(k)}x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}B^{\intercal}+\bar{A}_{t}^{(k)}x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(\bar{B}_{t}^{(k)})^{\intercal}$
$\displaystyle-A\mathds{E}\big{\\{}x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\big{\\}}B^{\intercal}-0-0-\mathds{E}\big{\\{}\bar{A}_{t}^{(k)}x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}(\bar{B}_{t}^{(k)})^{\intercal}\big{\\}}\bigg{)}$
$\displaystyle=(B\otimes
A)\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\big{\\}}\right)$
$\displaystyle\quad+(B\otimes\bar{A}_{t}^{(k)})\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right)+(\bar{B}_{t}^{(k)}\otimes
A)\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right)$
$\displaystyle\quad+(\bar{B}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right)-\mathds{E}\big{\\{}(\bar{B}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})\big{\\}}\mathds{E}\big{\\{}\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right)\big{\\}}$
$\displaystyle=(B\otimes
A)\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}-\mathds{E}\big{\\{}x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\big{\\}}\right)$
$\displaystyle\quad+(B\otimes\bar{A}_{t}^{(k)})\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right)+(\bar{B}_{t}^{(k)}\otimes
A)\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right)$
$\displaystyle\quad+(\bar{B}_{t}^{(k)}\otimes\bar{A}_{t}^{(k)})\operatorname{vec}\left(x_{t}^{(k)}({u_{t}^{(k)}})^{\intercal}\right).$

Taking norms, and substituting notated quantities, we have

$\displaystyle\left\|\boxed{2}\right\|$
$\displaystyle\leq\|A\|_{2}\|B\|_{2}\|\Delta
W_{t}\|+(\|A\|_{2}c_{\bar{B}}+\|B\|_{2}c_{\bar{A}}+c_{\bar{A}}c_{\bar{B}})c_{M}c_{U}.$

The third term of (44) is simply the transpose of the second term, so an
identical norm bound holds.

Putting together the four terms of (44), we obtain

$\displaystyle\|\Delta X_{t+1}\|$
$\displaystyle\leq(\|A\|_{2}^{2}+\|\Sigma_{A}^{\prime}\|_{2})\|\Delta
X_{t}\|+(2\|A\|_{2}c_{\bar{A}}+c_{\Sigma_{A}^{\prime}})c_{M}^{2}$
$\displaystyle\
+3(\|B\|_{2}^{2}+\|\Sigma_{B}^{\prime}\|_{2})c_{U}c_{\nu}+(2\|B\|_{2}c_{\bar{B}}+c_{\Sigma_{B}^{\prime}})c_{U}^{2}$
$\displaystyle\
+2\|A\|_{2}\|B\|_{2}c_{W}+(2\|A\|_{2}c_{\bar{B}}+2\|B\|_{2}c_{\bar{A}}+2c_{\bar{A}}c_{\bar{B}})c_{M}c_{U}.$
(45)

For the base case when $t=0$, we have by Assumption 2 (v) that $\|\Delta
X_{0}\|\leq c_{\Delta X}$. Applying (45) inductively with the base case proves
the fifth claim.

##  Appendix E Basic identities and inequalities

In the proofs of Theorems 3.16 and 3.17, the following facts will be used.

Submultiplicativity of spectral norm

For $A\in\mathds{R}^{m\times n}$ and $B\in\mathds{R}^{n\times p}$,

$\displaystyle\|AB\|_{2}\leq\|A\|_{2}\|B\|_{2}.$

Norm of Kronecker product

For $A\in\mathds{R}^{m\times n}$ and $B\in\mathds{R}^{p\times q}$,

$\displaystyle\|A\otimes B\|_{2}=\|A\|_{2}\|B\|_{2}.$

Inverse of spectral norm

For any invertible matrix $A\in\mathds{R}^{n\times n}$ we have

$\displaystyle\|A^{-1}\|_{2}=\frac{1}{\sqrt{\lambda_{\min}(AA^{\intercal})}}.$
(46)

Difference of matrix inverses

Suppose $A,E\in\mathds{R}^{n\times n}$ are invertible square matrices. Then

$\displaystyle
A^{-1}-E^{-1}=E^{-1}(((E-A)+A)A^{-1})-E^{-1}=E^{-1}(E-A)A^{-1}.$ (47)

Matrix inverse perturbation bound

(Equation (5.8.1) of (Horn and Johnson,, 2012))
Suppose $A,A+F\in\mathds{R}^{n\times n}$ are invertible square matrices. Then

$\displaystyle\|A^{-1}-(A+F)^{-1}\|\leq\|A^{-1}\|\cdot\|F\|\cdot\|(A+F)^{-1}\|.$
(48)

This follows from taking $E=A+F$ in (47) and using submultiplicativity of
spectral norm.

Probability bound on the sum of random variables

Consider $k$ random variables, $X_{1}$, $\dots$, $X_{k}$, and a positive
number $\varepsilon$. Note that $X_{i}<\varepsilon/k$ for all $i\in[k]$,
implies $\sum_{i=1}^{k}X_{i}<\varepsilon$, so it follows from the union bound
that

$\displaystyle\mathbb{P}\left\\{\sum_{i=1}^{k}X_{i}\geq\varepsilon\right\\}\leq\sum_{i=1}^{k}\mathbb{P}\\{X_{i}\geq\varepsilon/k\\}.$
(49)

Probability bound on the product of nonnegative random variables

Consider two nonnegative random variables, $X_{1}$ and $X_{2}$, and a positive
number $\varepsilon$. Since $X_{1}<\sqrt{\varepsilon}$ and
$X_{2}<\sqrt{\varepsilon}$ implies $X_{1}X_{2}<\varepsilon$ we have

$\displaystyle\mathbb{P}\\{X_{1}X_{2}\geq\varepsilon\\}\leq\mathbb{P}\\{X_{1}\geq\sqrt{\varepsilon}\\}+\mathbb{P}\\{X_{2}\geq\sqrt{\varepsilon}\\}.$
(50)

We will need the following geometrical result later in the use of covering
arguments.

######  Lemma E.21 (Covering numbers of the Euclidean Sphere).

Consider a minimal $\gamma$-net $\\{w_{k},k\in[M_{\gamma}]\\}$ of the
$n$-dimensional sphere surface
$\mathcal{S}_{n-1}:=\\{w\in\mathds{R}^{n}:\|w\|=1\\}$. That is, for all
$w\in\mathcal{S}_{n-1}$ there exists $w_{i}\in\\{w_{k},k\in[M_{\gamma}]\\}$
such that $\|w-w_{i}\|\leq\gamma$, and $M_{\gamma}$ is the smallest number
satisfies this condition. Then for any $\gamma>0$, the covering number
$M_{\gamma}$, i.e., the cardinality of the $\gamma$-net satisfies

$\displaystyle\left(\frac{1}{\gamma}\right)^{n}\leq
M_{\gamma}\leq\left(\frac{2}{\gamma}+1\right)^{n}.$ (51)

A standard volume comparison involving Euclidean balls, e.g. Corollary 4.2.13
of (Vershynin,, 2018), yields the result.

We also need the following matrix concentration inequality.

######  Lemma E.22 (Matrix Bernstein inequality (Tropp,, 2015)).

Consider a finite sequence of independent random matrices $\\{X_{k},$
$k\in[N]\\}$ with common dimension $m\times n$. Assume that
$\mathds{E}\\{X_{k}\\}=0$ and $\|X_{k}\|_{2}\leq L$, $k\in[N]$. Introduce
$S:=\sum_{k=1}^{N}X_{k}$ and let
$v\geq\max\\{\|\mathds{E}\\{SS^{\intercal}\\}\|_{2},\|\mathds{E}\\{S^{\intercal}S\\}\|_{2}\\}$.
Then, for all $\varepsilon\geq 0$,

$\displaystyle\mathds{P}\\{\|S\|_{2}\geq\varepsilon\\}\leq(n+m)\exp\left\\{-\frac{3}{2}\cdot\frac{\varepsilon^{2}}{3v+L\varepsilon}\right\\}.$
(52)

We obtain the following corollary from Lemma E.22.

######  Corollary E.23.

Consider a finite sequence of independent random matrices
$\\{Y_{k},k\in[N]\\}$ with common dimension $m\times n$. Assume that
$\mathds{E}\\{Y_{k}\\}=0$ and $\|Y_{k}\|_{2}\leq M$, $k\in[N]$. Then, for all
$\varepsilon\geq 0$,

$\displaystyle\mathds{P}\left\\{\left\|\frac{1}{N}\sum_{k=1}^{N}Y_{k}\right\|_{2}\geq\varepsilon\right\\}\leq\delta(\varepsilon),$

where

$\displaystyle\delta(\varepsilon):=(n+m)\exp\left\\{-\frac{3}{2}\cdot\frac{N\varepsilon^{2}}{3M^{2}+M\varepsilon}\right\\}.$

Towards application of Lemma E.22, assign $X_{k}=Y_{k}/N$ and thus $L=M/N$.
Now we get a crude bound on $v$ as

$\displaystyle\|\mathds{E}\\{SS^{\intercal}\\}\|_{2}$
$\displaystyle=\left\|\mathds{E}\left\\{\left(\sum_{k=1}^{N}X_{k}\right)\left(\sum_{j=1}^{N}X_{j}^{\intercal}\right)\right\\}\right\|_{2}$
$\displaystyle=\left\|\sum_{k=1}^{N}\sum_{j=1}^{N}\mathds{E}\left\\{X_{k}X_{j}^{\intercal}\right\\}\right\|_{2}$
(linearity of $\mathds{E}\\{\cdot\\}$)
$\displaystyle=\left\|\sum_{k=1}^{N}\mathds{E}\left\\{X_{k}X_{k}^{\intercal}\right\\}\right\|_{2}$
(mutual independence of $X_{k}$)
$\displaystyle\leq\sum_{k=1}^{N}\left\|\mathds{E}\left\\{X_{k}X_{k}^{\intercal}\right\\}\right\|_{2}$
(triangle inequality)
$\displaystyle\leq\sum_{k=1}^{N}\mathds{E}\left\\{\|X_{k}X_{k}^{\intercal}\|_{2}\right\\}$
(Jensen’s inequality)
$\displaystyle=\sum_{k=1}^{N}\mathds{E}\left\\{\|X_{k}\|_{2}^{2}\right\\}$
(property of the spectral norm)
$\displaystyle\leq\sum_{k=1}^{N}(M/N)^{2}=M^{2}/N.$

An identical argument shows $\|\mathds{E}\\{SS^{\intercal}\\}\|_{2}\leq
M^{2}/N$ so we can take $v=M^{2}/N$. Applying Lemma E.22 with $X_{k}=Y_{k}/N$,
$L=M/N$, and $v=M^{2}/N$ yields the claim.

##  Appendix F Proof of Theorem 3.16

In this section we obtain bounds for system parameter error matrix
$[\hat{A}~{}\hat{B}]-[A~{}B]$ by decomposing the difference using their
representation in the least-squares estimators as

$\displaystyle\big{[}\hat{A}~{}\hat{B}\big{]}-[A~{}B]$
$\displaystyle=\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger}-\mathbf{Y}\mathbf{Z}^{\intercal}(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}$
$\displaystyle=\big{[}\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Y}\mathbf{Z}^{\intercal}\big{]}(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}+\mathbf{Y}\mathbf{Z}^{\intercal}\big{[}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger}-(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}\big{]}+\big{[}\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Y}\mathbf{Z}^{\intercal}\big{]}\big{[}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger}-(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}\big{]}.$

In this form it is obvious that there are four unique terms, which fall into
two groups. The first group is
$\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Y}\mathbf{Z}^{\intercal}$
and
$(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger}-(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}$,
which represent error terms amenable to analysis. The second group is
$(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}$ and
$\mathbf{Y}\mathbf{Z}^{\intercal}$, which are inherent to the system and do
not depend on the estimator quality. The terms
$\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Y}\mathbf{Z}^{\intercal}$
and
$(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger}-(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}$
are treated first, then the bound on $[\hat{A}~{}\hat{B}]-[A~{}B]$ is
obtained.

Throughout this section, small probability bounds are denoted by
$\delta_{[\cdot]}$, where $[\cdot]$ are various subscripts, and each of these
bounds decreases monotonically towards $0$ with increasing number of rollouts
$n_{r}$.

The following bounds for $\hat{\mathbf{Y}}-\mathbf{Y}$ and
$\hat{\mathbf{Z}}-\mathbf{Z}$ follow from Corollary E.23.

######  Lemma F.24.

Suppose that Assumptions 1 and 2 hold. Then for all $\varepsilon>0$,

$\displaystyle\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}\geq\varepsilon\big{\\}}=\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}\geq\varepsilon\big{\\}}=\mathds{P}\Big{\\{}\big{[}\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}\geq\varepsilon\big{]}\bigcup\big{[}\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}\geq\varepsilon\big{]}\Big{\\}}$
$\displaystyle\leq\delta_{Y}(\varepsilon),$

where

$\displaystyle\delta_{Y}(\varepsilon)$
$\displaystyle:=(n+\ell)\exp\left\\{-\frac{3}{2}\cdot\frac{n_{r}\varepsilon^{2}}{3\ell
c_{N}^{2}+\varepsilon\sqrt{\ell c_{N}^{2}}}\right\\},$

Using the bound
$\|e_{t}^{(k)}\|=\|x_{t}^{(k)}-\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\|\leq
c_{N}$ from Lemma 3.14, and denoting
$\hat{\mathbf{Y}}_{k}:=\big{[}x_{\ell}^{(k)}~{}\cdots~{}x_{1}^{(k)}\big{]}$ so
$\hat{\mathbf{Y}}=(\sum\nolimits_{k=1}^{n_{r}}\hat{\mathbf{Y}}_{k})/n_{r}$ and
$\mathds{E}\big{\\{}\hat{\mathbf{Y}}_{k}\big{\\}}=\mathbf{Y}$, we obtain

$\displaystyle\big{\|}\hat{\mathbf{Y}}_{k}-\mathbf{Y}\big{\|}_{2}\leq\big{\|}\hat{\mathbf{Y}}_{k}-\mathbf{Y}\big{\|}_{F}=\sqrt{\sum_{t=1}^{\ell}\big{\|}x_{t}^{(k)}-\mathds{E}\big{\\{}x_{t}^{(k)}\big{\\}}\big{\|}^{2}}\leq\sqrt{\ell
c_{N}^{2}}.$

Applying Corollary E.23 with $Y_{k}=\hat{\mathbf{Y}}_{k}-\mathbf{Y}$,
$N=n_{r}$, and $M=\sqrt{\ell c_{N}^{2}}$, we conclude

$\displaystyle\mathds{P}\\{\|\hat{\mathbf{Y}}-\mathbf{Y}\|_{2}\geq\varepsilon\\}$
$\displaystyle\leq\delta_{Y}(\varepsilon).$

Denote
$\hat{\mathbf{Z}}_{k}=\begin{bmatrix}x_{\ell-1}^{(k)}&\cdots&x_{0}^{(k)}\\\
\nu_{\ell-1}&\cdots&\nu_{0}\end{bmatrix}$ so
$\hat{\mathbf{Z}}=(\sum\nolimits_{k=1}^{n_{r}}\hat{\mathbf{Z}}_{k})/n_{r}$ and
$\mathds{E}\big{\\{}\hat{\mathbf{Z}}_{k}\big{\\}}=\mathbf{Z}$. Noticing that
the last $m$ rows of $\hat{\mathbf{Z}}_{k}-\mathbf{Z}$ all have zero entries,
we have that

$\displaystyle\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}=\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}.$

Hence the events
$\big{\\{}\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}\geq\varepsilon\big{\\}}$
and
$\big{\\{}\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}\geq\varepsilon\big{\\}}$
are precisely the same, concluding the proof.

######  Remark F.25.

The reason that the last $m$ rows of $\hat{\mathbf{Z}}_{k}-\mathbf{Z}$ all
have zero entries is that the first moments of the inputs are known in
Algorithm 1, and appear identically in both $\hat{\mathbf{Z}}_{k}$ and
$\mathbf{Z}$. Hence the probability bounds are independent of the input
dimension $m$.

######  Lemma F.26.

Suppose that Assumptions 1 and 2 hold. Then for all $\varepsilon>0$,

$\displaystyle\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Y}\mathbf{Z}^{\intercal}\big{\|}_{2}\geq\varepsilon\big{\\}}\leq\delta_{YZ}(\varepsilon)$

where

$\displaystyle\delta_{YZ}(\varepsilon)$
$\displaystyle:=\delta_{Y}\left(\sqrt{\varepsilon+\left(\frac{\|\mathbf{Y}\|_{2}+\|\mathbf{Z}\|_{2}}{2}\right)^{2}}-\frac{\|\mathbf{Y}\|_{2}+\|\mathbf{Z}\|_{2}}{2}\right).$

We begin with the decomposition

$\displaystyle\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Y}\mathbf{Z}^{\intercal}=(\hat{\mathbf{Y}}-\mathbf{Y})(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}+(\hat{\mathbf{Y}}-\mathbf{Y})\mathbf{Z}^{\intercal}+\mathbf{Y}(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}.$

By the triangle inequality and submultiplicativity we have

$\displaystyle\big{\|}\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Y}\mathbf{Z}^{\intercal}\big{\|}_{2}$
$\displaystyle=\big{\|}(\hat{\mathbf{Y}}-\mathbf{Y})(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}+(\hat{\mathbf{Y}}-\mathbf{Y})\mathbf{Z}^{\intercal}+\mathbf{Y}(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}\big{\|}_{2}$
$\displaystyle\leq\big{\|}(\hat{\mathbf{Y}}-\mathbf{Y})(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}\big{\|}_{2}+\big{\|}(\hat{\mathbf{Y}}-\mathbf{Y})\mathbf{Z}^{\intercal}\big{\|}_{2}+\big{\|}\mathbf{Y}(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}\big{\|}_{2}$
$\displaystyle\leq\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}+\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}\|\mathbf{Z}\|_{2}+\|\mathbf{Y}\|_{2}\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}$
$\displaystyle=\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}^{2}+\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}\left(\left\|\mathbf{Y}\right\|_{2}+\left\|\mathbf{Z}\right\|_{2}\right).$

Considering a probability bound, solving the quadratic inequality in
$\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}$, and applying Lemma F.24 we
have

$\displaystyle\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Y}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Y}\mathbf{Z}^{\intercal}\big{\|}_{2}\geq\varepsilon\big{\\}}$
$\displaystyle\leq\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}^{2}+\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}(\|\mathbf{Y}\|_{2}+\|\mathbf{Z}\|_{2})\geq\varepsilon\big{\\}}$
$\displaystyle=\mathds{P}\left\\{\big{\|}\hat{\mathbf{Y}}-\mathbf{Y}\big{\|}_{2}\geq\frac{1}{2}(\|\mathbf{Y}\|_{2}+\|\mathbf{Z}\|_{2})\left(\sqrt{1+\frac{4\varepsilon}{(\|\mathbf{Y}\|_{2}+\|\mathbf{Z}\|_{2})^{2}}}-1\right)\right\\}$
$\displaystyle\leq\delta_{Y}\left(\frac{\|\mathbf{Y}\|_{2}+\|\mathbf{Z}\|_{2}}{2}\left(\sqrt{1+\frac{4\varepsilon}{(\|\mathbf{Y}\|_{2}+\|\mathbf{Z}\|_{2})^{2}}}-1\right)\right),$

which was the claimed inequality.

######  Lemma F.27.

Suppose Assumptions 1 and 2 hold. Given a positive constant
$\varepsilon_{\max}$, then for all $0<\varepsilon<\varepsilon_{\max}$,

$\displaystyle\mathds{P}\big{\\{}\big{\|}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger}-(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}\big{\|}_{2}\geq\varepsilon\big{\\}}\leq\delta_{ZZ}(\varepsilon,\varepsilon_{\max}),$

where

$\displaystyle\delta_{ZZ}(\varepsilon,\varepsilon_{\max})$
$\displaystyle:=\delta_{0}\left(\frac{1}{2}\lambda_{\min}^{2}(\mathbf{Z}\mathbf{Z}^{\intercal})\left(1-\frac{\varepsilon}{\varepsilon_{\max}}\right)\varepsilon\right)+\delta_{m}\left(\frac{\varepsilon\lambda_{\min}(\mathbf{Z}\mathbf{Z}^{\intercal})}{\varepsilon_{\max}(2+\lambda_{\min}(\mathbf{Z}\mathbf{Z}^{\intercal})/\lambda_{\max}(\mathbf{Z}\mathbf{Z}^{\intercal}))}\right),$
$\displaystyle\delta_{0}(\varepsilon)$
$\displaystyle:=\delta_{Y}\left(\sqrt{\lambda_{\max}(\mathbf{Z}\mathbf{Z}^{\intercal})+\varepsilon}-\sqrt{\lambda_{\max}(\mathbf{Z}\mathbf{Z}^{\intercal})}\right),$
$\displaystyle\delta_{m}(\varepsilon)$
$\displaystyle:=\left(9^{n+m}+\left(\frac{16\lambda_{\max}(\mathbf{Z}\mathbf{Z}^{\intercal})}{\lambda_{\min}(\mathbf{Z}\mathbf{Z}^{\intercal})}+1\right)^{n+m}\right)\delta_{0}(\varepsilon).$

######  Remark F.28.

The additional parameter $\varepsilon_{\max}$ arises when bounding
$\lambda_{\min}^{2}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})$, and is in
fact independent of the estimation bound.

Later we will show that $\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}$ is
invertible with high probability, so we now assume the existence of
$(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}$, i.e.,
$(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{\dagger}=(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}$.
Hence we may apply (48) to obtain the decomposition

$\displaystyle\big{\|}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}-(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}\big{\|}_{2}$
$\displaystyle=\big{\|}[(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}\big{(}\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}\big{)}\big{\|}_{2}$
$\displaystyle\leq\big{\|}(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}\big{\|}_{2}\big{\|}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}\big{\|}_{2}\big{\|}\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}\big{\|}_{2}.$

Considering a probability bound, rearranging, and using (46) we obtain

$\displaystyle\mathds{P}\big{\\{}\big{\|}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}-(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}\big{\|}_{2}\geq\varepsilon\big{\\}}$
$\displaystyle\leq\mathds{P}\big{\\{}\big{\|}(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}\big{\|}_{2}\big{\|}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}\big{\|}_{2}\big{\|}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})-(\mathbf{Z}\mathbf{Z}^{\intercal})\big{\|}_{2}\geq\varepsilon\big{\\}}$
$\displaystyle=\mathds{P}\left\\{\big{\|}\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}\big{\|}_{2}\geq\frac{\varepsilon}{\big{\|}(\mathbf{Z}\mathbf{Z}^{\intercal})^{-1}\big{\|}_{2}\big{\|}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})^{-1}\big{\|}_{2}}\right\\}$
$\displaystyle=\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}\big{\|}_{2}\geq\varepsilon\lambda_{\min}(\mathbf{Z}\mathbf{Z}^{\intercal})\lambda_{\min}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})\big{\\}}.$
(53)

We are now faced with providing bounds on both
$\big{\|}\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}\big{\|}_{2}$
on the left side and
$\lambda_{\min}(\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal})$ on the right
side of the inequality inside the probability.

First, we consider the bound of
$\big{\|}\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}\big{\|}_{2}$
by beginning with the decomposition

$\displaystyle\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}=(\hat{\mathbf{Z}}-\mathbf{Z})(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}+(\hat{\mathbf{Z}}-\mathbf{Z})\mathbf{Z}^{\intercal}+\mathbf{Z}(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}.$

Using the triangle inequality, submultiplicativity, and solving the quadratic
inequality in $\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}$, we obtain

$\displaystyle\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}\big{\|}_{2}\geq\varepsilon\big{\\}}$
$\displaystyle=\mathds{P}\big{\\{}\big{\|}(\hat{\mathbf{Z}}-\mathbf{Z})(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}+(\hat{\mathbf{Z}}-\mathbf{Z})\mathbf{Z}^{\intercal}+\mathbf{Z}(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}\big{\|}_{2}\geq\varepsilon\big{\\}}$
$\displaystyle\leq\mathds{P}\big{\\{}\big{\|}(\hat{\mathbf{Z}}-\mathbf{Z})(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}\big{\|}_{2}+\big{\|}(\hat{\mathbf{Z}}-\mathbf{Z})\mathbf{Z}^{\intercal}\big{\|}_{2}+\big{\|}\mathbf{Z}(\hat{\mathbf{Z}}-\mathbf{Z})^{\intercal}\big{\|}_{2}\geq\varepsilon\big{\\}}$
$\displaystyle\leq\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}^{2}+2\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}\big{\|}\mathbf{Z}\big{\|}_{2}\geq\varepsilon\big{\\}}$
$\displaystyle=\mathds{P}\left\\{\big{\|}\hat{\mathbf{Z}}-\mathbf{Z}\big{\|}_{2}\geq\|\mathbf{Z}\|_{2}\left(\sqrt{1+\frac{\varepsilon}{\|\mathbf{Z}\|_{2}^{2}}}-1\right)\right\\}.$

Applying Lemma F.24 with the appropriate settings of $\varepsilon$ yields

$\displaystyle\mathds{P}\big{\\{}\big{\|}\hat{\mathbf{Z}}\hat{\mathbf{Z}}^{\intercal}-\mathbf{Z}\mathbf{Z}^{\intercal}\big{\|}_{2}\geq\varepsilon\big{\\}}$
$\displaystyle\leq\delta_{Y}\left(\|\mathbf{Z}\|_{2}\left(\sqrt{1+\frac{\varepsilon}{\|\mathbf{Z}\|_{2}^{2}}}-1\right)\right)=:\delta_{0}(\varepsilon).$
(54)

Now we seek a lower bound of