# Efficient Learning of a Linear Dynamical System
with Stability Guarantees

Wouter Jongeneel  Tobias Sutter  Daniel Kuhn

( Risk Analytics and Optimization Chair, École Polytechnique Fédérale de
Lausanne,
{wouter.jongeneel, tobias.sutter<EMAIL_ADDRESS>)

###### Abstract

We propose a principled method for projecting an arbitrary square matrix to
the non-convex set of asymptotically stable matrices. Leveraging ideas from
large deviations theory, we show that this projection is optimal in an
information-theoretic sense and that it simply amounts to shifting the initial
matrix by an optimal linear quadratic feedback gain, which can be computed
exactly and highly efficiently by solving a standard linear quadratic
regulator problem. The proposed approach allows us to learn the system matrix
of a stable linear dynamical system from a single trajectory of correlated
state observations. The resulting estimator is guaranteed to be stable and
offers explicit statistical bounds on the estimation error.

##  1 Introduction

We study the problem of learning a stable linear dynamical system from a
single trajectory of correlated state observations. This problem is of
fundamental importance in various disciplines such as adaptive control [2],
system identification [34, 66], reinforcement learning [57, 10, 43, 51] and
approximate dynamic programming [7, 50]. Specifically, we consider a discrete-
time linear time-invariant system of the form

$x_{t+1}=\theta x_{t}+w_{t},\quad x_{0}\sim\nu,$ (1.1)

where $x_{t}\in\mathbb{R}^{n}$ and $w_{t}\in\mathbb{R}^{n}$ denote the state
and the exogenous noise at time $t\in\mathbb{N}$, respectively, while $\theta$
represents the system matrix, and $\nu$ stands for the marginal distribution
of the initial state $x_{0}$. We assume that $\theta$ is asymptotically
stable, that is, it belongs to $\Theta=\\{{\theta}\in\mathbb{R}^{n\times
n}:\rho({\theta})<1\\}$, where $\rho({\theta})$ denotes the largest absolute
eigenvalue of ${\theta}$. For ease of terminology, we will usually refer to
$\Theta$ as the set of stable matrices and to its complement in
$\mathbb{R}^{n\times n}$ as the set of unstable matrices. We assume that
nothing is known about $\theta$ except for its membership in $\Theta$, and we
aim to learn $\theta$ from a single-trajectory of data
$\\{\widehat{x}_{t}\\}_{t=0}^{T}$ generated by (1.1). To this end, one can use
the least squares estimator

$\widehat{\theta}_{T}=\left(\textstyle\sum_{t=1}^{T}\widehat{x}_{t}\widehat{x}_{t-1}^{\mathsf{T}}\right)\left(\textstyle\sum_{t=1}^{T}\widehat{x}_{t-1}\widehat{x}_{t-1}^{\mathsf{T}}\right)^{-1},$
(1.2)

which may take any value in $\Theta^{\prime}=\mathbb{R}^{n\times n}$ under
standard assumptions on the noise distribution. It is therefore possible that
$\widehat{\theta}_{T}\notin\Theta$ even though $\theta\in\Theta$. This is
troubling because stability is important in many applications, for example,
when the estimated model is used for prediction, filtering or control.
Estimating a stable model is also crucial for assessing the performance of a
stable system or, in a simulation context, for generating useful exploration
data.

Given the prior structural information that $\theta$ is stable, we thus seek
an estimator that is guaranteed to preserve stability. A natural approach to
achieve this goal would be to project the least squares estimator
$\widehat{\theta}_{T}$ to the nearest stable matrix with respect to some
discrepancy function on $\mathbb{R}^{n\times n}$. This seems challenging,
however, because $\Theta$ is open, unbounded and non-convex; see Figure
1.1($a$). To circumvent this difficulty, we introduce a new discrepancy
function that adapts to the geometry of $\Theta$ and is thus ideally suited
for projecting unstable matrices onto $\Theta$. We will characterize the
statistical properties of this projection when applied to the least squares
estimator, and we will show that it can be computed efficiently even for
systems with $O(10^{3})$ states.

The following example shows that naïve heuristics to project $\theta^{\prime}$
into the interior of $\Theta$ could spectacularly fail.

######  Example 1.1 (Projection by eigenvalue scaling).

A naïve method to stabilize a matrix $\theta^{\prime}\notin\Theta$ would be to
scale its unstable eigenvalues into the complex unit circle. To see that the
output of this transformation may not retain much similarity with the input
$\theta^{\prime}$, consider the matrices

${\theta}^{\prime}\\!=\\!\begin{bmatrix}1.01&10\\\
.01&1\end{bmatrix},\;{\theta}^{\prime}_{a}\\!=\\!\begin{bmatrix}.84&4.77\\\
.005&.84\end{bmatrix},\;{\theta}^{\prime}_{b}\\!=\\!\begin{bmatrix}.99&10\\\
0&.99\end{bmatrix}.$

Clipping off the unstable eigenvalues of $\theta^{\prime}$ at $|\lambda|=.99$
yields ${\theta}^{\prime}_{a}$ with $\rho(\theta^{\prime}_{a})=.99$ and
$\|{\theta}^{\prime}-{\theta}^{\prime}_{a}\|_{2}\gtrsim 5$. However, the
matrix $\theta^{\prime}_{b}$ also has spectral radius
$\rho(\theta^{\prime}_{b})=.99$ but is much closer to $\theta^{\prime}$.
Indeed, we have $\|{\theta}^{\prime}-{\theta}^{\prime}_{b}\|_{2}\approx 0.02$.

Figure 1.1: $(a)$ Visualization of a two-dimensional slice of the set $\Theta$
of stable $2\times 2$-matrices with diagonal elements
$\theta_{11}=\theta_{22}=\frac{1}{2}$. The set of all corresponding off-
diagonal elements is open, non-convex and unbounded along the coordinate axes.
$(b)$ Some eigenvalues and the trace of the covariance matrix $S_{\theta}$ of
the invariant state distribution explode as $\theta$ approaches the boundary
of $\Theta$.

###  1.1 Related work

The problem of learning a stable dynamical system is widely studied in
dynamical systems theory, while the problem of projecting an unstable matrix
onto $\Theta$ with respect to some norm has attracted considerable interest in
matrix analysis.

In dynamical systems theory, [40] proposed one of the first methods to project
a possibly unstable estimator onto $\Theta$ by using subspace methods. This
pioneering approach has significant practical merits [64] but may also
significantly distort the original estimator. To overcome this deficiency,
[35] approximate $\Theta$ by the set of contractive matrices whose operator
norm is at most $1$. While this set is convex, it offers but a conservative
approximation of $\Theta$. Several related methods have since been proposed to
enforce stability [36, 13, 59], which are all either conservative or
computationally expensive. Moreover, these methods do not provide any
statistical guarantees. [62, 63] regularize the least squares objective and
show that the spectral radius of the resulting estimator is bounded by a
function of the regularization weights. As $\Theta$ is an open set, however,
the tuning of these weights remains a matter of taste. More recently, [60]
propose a maximum likelihood approach that is attractive from a statistical
point of view but can be computationally challenging in certain applications.
On the other hand, several authors use Lyapunov theory to provide stability
guarantees for deterministic vector fields; see, e.g., [45, 5, 32, 61]. A more
recent approach by [12] learns stability certificates from i.i.d.
trajectories. There is also a substantial body of literature on (sub-)optimal
finite-sample concentration bounds for linear systems identified via least
squares estimation [54, 27, 52, 28, 53]. These approaches offer fast learning
rates but cannot guarantee stability of the identified systems for finite
sample sizes.

Much like in dynamical systems theory, in matrix analysis one seeks algorithms
for projecting an unstable deterministic matrix $\theta^{\prime}$ onto
$\Theta$, which is equivalent to finding the smallest additive perturbation
that stabilizes $\theta^{\prime}$. More specifically, matrix analysis studies
the nearest stable matrix problem

$\Pi_{\Theta}(\theta^{\prime})\in\arg\min_{\theta\in\mathop{\mathsf{cl}}\Theta}\|\theta^{\prime}-\theta\|^{2},$
(1.3)

where $\|\cdot\|$ represents a prescribed norm on $\mathbb{R}^{n\times n}$.
Note that optimizing over the closure of $\Theta$ is necessary for (1.3) to be
well-defined because any minimizer lies on the boundary of the open set
$\Theta$. Solving (1.3) is challenging because $\Theta$ is non-convex.
Existing numerical solution procedures rely on successive convex
approximations [47], on local optimization schemes based on the solution of
low-rank matrix differential equations [23] or on an elegant reparametrization
of the set of stable matrices, which simplifies the numerics of the projection
operation [21, 15]. The latter approach was recently used for learning stable
systems [42, 41]. [46] solve (1.3) for certain polyhedral norms and non-
negative matrices $\theta^{\prime}$, which allows them to find exact
solutions. See [25] for a general discussion on matrix nearness problems.

Optimal control offers a promising alternative perspective on problem (1.3),
which is closely related to the approach advocated in this paper: one could
try to design a linear quadratic regulator (LQR) problem whose optimal
feedback gain $K^{\star}\in\mathbb{R}^{n\times n}$ renders
$\theta^{\prime}+K^{\star}$ stable. By proposing an LQR objective that is
inversely proportional to the sample covariance matrix of the measurement
noise, [58] show that this idea is indeed valid, but they provide no error
analysis or statistical guarantees. Using optimal control techniques, [29]
prove that one can find matrices $K$ that not only render $\theta^{\prime}+K$
stable but are also structurally equivalent to $\theta^{\prime}$ (e.g.,
$\theta^{\prime}+K$ preserves the null space of $\theta^{\prime}$). Such a
structure-preserving approach seems preferable over the plain nearest stable
matrix problem (1.3), which merely seeks stability at minimal ‘cost.’
Appealing to the theory of large deviations, we will give such approaches a
statistical underpinning.

Notation. For a matrix $A\in\mathbb{C}^{n\times n}$, we denote by $\rho(A)$
the largest absolute eigenvalue and by $\kappa(A)$ the condition number of
$A$. For a set $\mathcal{D}\subset\mathbb{R}^{n}$, we denote by
$\mathcal{D}^{\mathsf{c}}$ the complement, by $\mathsf{cl}\,\mathcal{D}$ the
closure and by $\mathsf{int}\,\mathcal{D}$ the interior of $\mathcal{D}$. For
a real sequence $\\{a_{T}\\}_{T\in\mathbb{N}}$ we use $1\ll a_{T}\ll T$ to
express that $a_{T}/T\to 0$ and $a_{T}\to\infty$ as $T\to\infty$.

###  1.2 Contributions

Even though the least squares estimator $\widehat{\theta}_{T}$ is strongly
consistent and thus converges $\mathbb{P}_{\theta}$-almost surely to $\theta$
[14], it differs $\mathbb{P}_{\theta}$-almost surely from $\theta$ for any
finite $T$. In order to quantify estimation errors, we introduce a discrepancy
function $I:\Theta^{\prime}\times\Theta\rightarrow[0,\infty]$ defined through

$I(\theta^{\prime},\theta)=\frac{1}{2}\mathrm{tr}\left(S_{w}^{-1}(\theta^{\prime}-\theta)S_{\theta}(\theta^{\prime}-\theta)^{\mathsf{T}}\right).$
(1.4)

Here, $S_{\theta}$ denotes the covariance matrix of $x_{t}$ under the
stationary state distribution, which exists for $\theta\in\Theta$ but diverges
as $\theta$ approaches the boundary of $\Theta$; see Figure 1.1($b$). Note
that since $S_{w}\succ 0$ and hence $S_{\theta}\succ 0$,
$I(\theta^{\prime},\theta)$ vanishes if and only if $\theta^{\prime}=\theta$.
In this sense $I$ behaves like a distance. Note, however, that
$I(\theta^{\prime},\theta)$ is not symmetric in $\theta$ and
$\theta^{\prime}$. In this paper we propose to use the discrepancy function
(1.4) for projecting an unstable matrix $\theta^{\prime}$ onto $\Theta$.
Specifically, we define the reverse $I$-projection of any
$\theta^{\prime}\in\mathbb{R}^{n\times n}$ as

$\mathcal{P}(\theta^{\prime})\in\arg\min_{\theta\in\Theta}I(\theta^{\prime},{\theta}).$
(1.5)

We emphasize that the minimum in (1.5) is always attained even though $\Theta$
is open. The reason for this is that $S_{\theta}$, and thus also
$I(\theta^{\prime},\theta)$, diverges as $\theta$ approaches the boundary of
$\Theta$; see Proposition 3.6 below. Thus, the minimum must be attained inside
$\Theta$. In fact, as $I(\theta^{\prime},{\theta})$ trades off distance
against stability, $\mathcal{P}(\theta^{\prime})$ may not even be close to the
boundary of $\Theta$; see Figure 1.2. Moreover, we will see that the
discrepancy function (1.4) has a natural statistical interpretation, which
enables us to derive strong statistical guarantees for the reverse
$I$-projection of the least squares estimator.

We will actually show that the discrepancy function (1.4) determines the speed
at which the probabilities of untypical realizations of the least squares
estimator $\widehat{\theta}_{T}$ decay with $T$. More precisely, if
$\\{a_{T}\\}_{T\in\mathbb{N}}$ is a real sequence with $1\ll a_{T}\ll T$ and
$I(\widehat{\theta}_{T},\theta)\leq r(a_{T}/T)$ for some $r>0$, then

$\mathbb{P}_{\theta}(\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)\in U)\geq
e^{-ra_{T}+o(a_{T})}$ (1.6)

for any open neighbourhood $U\subset\mathbb{R}^{n\times n}$ of $0$. Hence, the
smaller $r$, the more unlikely it is that $\widehat{\theta}_{T}$ will rapidly
change, and $r$ characterizes the rate of change. Adopting an optimistic
perspective, assume for the moment that the underlying model $\theta\in\Theta$
maximizes the lower-bound in (1.6), that is, $r\geq 0$ is made as small as
possible for the given value of $\widehat{\theta}_{T}$. This is indeed
accomplished by setting $\theta=\mathcal{P}(\widehat{\theta}_{T})$. Note that
$\mathcal{P}(\widehat{\theta}_{T})=\widehat{\theta}_{T}$ whenever
$\widehat{\theta}_{T}\in\Theta$ by the definition of the reverse
$I$-projection. If $\widehat{\theta}_{T}\notin\Theta$, on the other hand, then
selecting $\theta$ optimistically as described corresponds to minimizing the
rate of change of $\widehat{\theta}_{T}$. Hence, the reverse $I$-projection of
$\widehat{\theta}_{T}$ can be viewed as being optimistic and has the flavor of
a maximum likelihood estimator.

Figure 1.2: Schematic visualization of $\Pi_{\Theta}(\theta^{\prime})$ and
$\mathcal{P}(\theta^{\prime})$ for different estimator realizations
$\theta^{\prime}$ inside and outside of $\Theta$.

Our main contributions can be summarized as follows.

  1. (i)

We prove that the discrepancy function (1.4) has a natural statistical
interpretation as the rate function of a moderate deviation principle for the
transformed least squares estimators
$\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)+\theta$, $T\in\mathbb{N}$.

  2. (ii)

We derive finite-sample and asymptotic statistical error bounds on the
operator norm distance between the reverse $I$-projection
$\mathcal{P}(\widehat{\theta}_{T})$ of the least squares estimator
$\widehat{\theta}_{T}$ and the unknown true system matrix $\theta$.

  3. (iii)

We show that the reverse $I$-projection $\mathcal{P}(\theta^{\prime})$ can be
computed highly efficiently to within any desired accuracy by solving a
standard LQR problem, e.g., via numerical routines that are readily available
in MATLAB or Julia. This method finds the ‘cheapest’ feedback gain matrix
$K^{\star}$ that renders $\theta^{\prime}+K^{\star}$ stable, and it can
evaluate $\mathcal{P}(\theta^{\prime})$ in under a minute even if $n\approx
10^{3}$.

In addition, numerical experiments corroborate our theoretical results and
showcase the statistical and computational merits of using the reverse
$I$-projection of $\widehat{\theta}_{T}$ to estimate $\theta$. To our best
knowledge, we present the first method for the identification of a linear
dynamical system with stability guarantees that is both computationally
efficient and offers asymptotic consistency and tight statistical error
bounds.

##  2 Main results

Throughout the paper we assume that all random objects are defined on a
measurable space $(\Omega,\mathcal{F})$ equipped with a probability measure
$\mathbb{P}_{\theta}$ that depends parametrically on $\theta$. We denote the
expectation operator with respect to $\mathbb{P}_{\theta}$ by
$\mathbb{E}_{\theta}[\cdot]$. From now on we impose the following assumption.

######  Assumption 2.1 (Linear system).

The following hold.

  1. (i)

The linear system (1.1) is stable, i.e., $\theta\in\Theta$.

  2. (ii)

For each $\theta\in\Theta$ the disturbances $\\{w_{t}\\}_{t\in\mathbb{N}}$ are
independent and identically distributed (i.i.d.) and independent of $x_{0}$
under $\mathbb{P}_{\theta}$. The marginal noise distributions are unbiased
($\mathbb{E}_{\theta}[w_{t}]=0$), non-degenerate
($S_{w}=\mathbb{E}_{\theta}[w_{t}w_{t}^{\mathsf{T}}]\succ 0$ is finite) and
have an everywhere positive probability density function.

Assumption 2.1 ensures that the linear system (1.1) admits an invariant
distribution $\nu_{\theta}$ [44, § 10.5.4]. This means that
$x_{t}\sim\nu_{\theta}$ implies $x_{t+1}\sim\nu_{\theta}$ for any
$t\in\mathbb{N}$. Moreover, as the probability density function of $w_{t}$ is
everywhere positive, $\\{x_{t}\\}_{t\in\mathbb{N}}$ represents a uniformly
ergodic Markov process, which implies that the marginal distribution of
$x_{t}$ under $\mathbb{P}_{\theta}$ converges weakly to $\nu_{\theta}$ as $t$
tends to infinity [44, Theorems 16.5.1 and 16.2.1]. Assumption 2.1 then
implies that the mean vector of $\nu_{\theta}$ vanishes and that the
covariance matrix $S_{\theta}$ of $\nu_{\theta}$ coincides with the unique
solution of the discrete Lyapunov equation

$S_{\theta}=\theta S_{\theta}\theta^{\mathsf{T}}+S_{w},$ (2.1)

which provides for a convenient way to compute $S_{\theta}$; see, e.g., [1, §
6.10 E]. Recall that $S_{\theta}$ critically enters the discrepancy function
$I(\theta^{\prime},\theta)$ defined in (1.4) and thus also the reverse
$I$-projection defined in (1.5). The following main theorem summarizes the key
statistical and computational properties of the reverse $I$-projection that
will be proved in the remainder of the paper. This theorem involves the
function $\mathsf{dlqr}(A,B,Q,R)$, which outputs the optimal feedback gain
matrix of an infinite-horizon deterministic LQR problem in discrete time. Such
problems are described by two system matrices $A$ and $B$, a state cost matrix
$Q\succeq 0$ and an input cost matrix $R\succ 0$ of compatible dimensions that
satisfy standard stabilizability and detectability conditions [8, § 4].

######  Theorem 2.2 (Efficient guaranteed stable identification).

Suppose that Assumption 2.1 holds and that $\widehat{\theta}_{T}$ is the least
squares estimator (1.2). Then, for any fixed $\theta\in\Theta$ the reverse
$I$-projection defined in (1.5) displays the following properties.

  1. (i)

Asymptotic consistency.

$\lim_{T\to\infty}\mathcal{P}(\widehat{\theta}_{T})=\theta\quad\mathbb{P}_{\theta}\text{-a.s.}$

  2. (ii)

Finite sample guarantees. For any $\beta\in(0,1)$ and sample size
$T\geq\kappa(S_{w}){\log(1/\beta)}/({\varepsilon^{2}C(n)})$ we have

$\mathbb{P}_{\theta}\left(\|\theta-\mathcal{P}(\widehat{\theta}_{T})\|_{2}\leq
2\varepsilon\right)\geq 1-\beta,$

where $\varepsilon=[2\kappa(S_{w})\cdot
I(\widehat{\theta}_{T},\mathcal{P}(\widehat{\theta}_{T}))]^{1/2}$, and $C(n)$
is a constant that depends only on the dimension $n$.

  3. (iii)

Efficient computation. For any $\theta^{\prime}\notin\Theta$ and $S_{w},Q\succ
0$ there is a $p\geq 1$, such that for all $\delta>0$ we have that

$\theta^{\star}_{\delta}=\theta^{\prime}+\mathsf{dlqr}(\theta^{\prime},I_{n},Q,(2\delta
S_{w})^{-1})$

is stable and satisfies
$\|\mathcal{P}(\theta^{\prime})-\theta^{\star}_{\delta}\|_{2}\leq
O(\delta^{p})$.

The asymptotic consistency (i) formalizes the intuitive requirement that more
data is preferable to less data. We emphasize that the reverse $I$-projection
does not introduce unnecessary bias because
$\mathcal{P}(\theta^{\prime})=\theta^{\prime}$ if $\theta^{\prime}$ is already
stable. The finite sample guarantees (ii) stipulate that the projected least
squares estimator $\mathcal{P}(\widehat{\theta}_{T})$ is guaranteed to be
close to the (unknown) true stable matrix $\theta$ with high probability
$1-\beta$. Note that if the observed state trajectory
$\\{\widehat{x}_{t}\\}_{t=0}^{T}$ is generated under $\mathbb{P}_{\theta}$,
then the inverse matrix appearing in (1.2) exists $\mathbb{P}_{\theta}$-almost
surely for any sample size $T\geq n$ thanks to Assumption 2.1 (ii). The
efficient computability property (iii), finally, shows that computing the
reverse $I$-projection to within high accuracy is no harder than solving a
standard LQR problem. The function $\mathsf{dlqr}$ is readily available as a
standard routine in software packages such as MATLAB or Julia. We also
emphasize that setting $Q=I_{n}$ works well in practice, that is, no tuning is
required to compute $\mathcal{P}(\theta^{\prime})$. Recall from Section 1.2
that the reverse $I$-projection exhibits optimism in the face of uncertainty,
a decision-making paradigm that is used with great success in various
reinforcement learning applications [38]. In general, however, optimism in the
face of uncertainty leads to computational intractability [14]. Thus, the
tractability result of Theorem 2.2 (iii) is a perhaps unexpected exception to
this rule; see Proposition 3.11 below for further details.

In the remainder of the paper we will prove Theorem 2.2. The proofs of
auxiliary results are relegated to the appendix (the supplementary material).

##  3 Reverse $I$-projection

We now demonstrate that the discrepancy function (1.4) underlying the reverse
$I$-projection has a natural statistical interpretation, which is crucial for
the proof of Theorem 2.2.

###  3.1 Moderate Deviations Theory

We leverage recent results from moderate deviations theory to show that the
discrepancy function (1.4) is intimately related to the least squares
estimator (1.2). To this end, we first introduce the basic notions of a rate
function and a moderate deviation principle. For a comprehensive introduction
to moderate deviations theory we refer to [26, 19].

######  Definition 3.1 (Rate function).

An extended real-valued function
$I:\Theta^{\prime}\times\Theta\rightarrow[0,\infty]$ is called a rate function
if it is lower semi-continuous in its first argument.

######  Definition 3.2 (Moderate deviation principle).

A sequence of estimators $\\{\widehat{\theta}_{T}\\}_{T\in\mathbb{N}}$ is said
to satisfy a moderate deviation principle with rate function $I$ if for every
sequence $\\{a_{T}\\}_{T\in\mathbb{N}}$ of real numbers with $1\ll a_{T}\ll
T$, for every Borel set $\mathcal{D}\subset\Theta^{\prime}$ and for every
$\theta\in\Theta$ we have

$\displaystyle-\inf_{\theta^{\prime}\in\mathsf{int}{\mathcal{D}}}\,I(\theta^{\prime},\theta)\leq$
$\displaystyle\liminf_{T\to\infty}\leavevmode\nobreak\
\frac{1}{a_{T}}\log\mathbb{P}_{\theta}\left(\widehat{\theta}_{T}\in\mathcal{D}\right)$
(3.1a) $\displaystyle\leq$
$\displaystyle\limsup_{T\to\infty}\leavevmode\nobreak\
\frac{1}{a_{T}}\log\mathbb{P}_{\theta}\left(\widehat{\theta}_{T}\in\mathcal{D}\right)$
(3.1b) $\displaystyle\leq$
$\displaystyle-\inf_{\theta^{\prime}\in\mathsf{cl}{\mathcal{D}}}\,I(\theta^{\prime},\theta).$
(3.1c)

If the rate function $I(\theta^{\prime},\theta)$ is continuous in
$\theta^{\prime}$ and the interior of $\mathcal{D}$ is dense in $\mathcal{D}$,
then the infima in (3.1a) and (3.1c) coincide, which implies that all
inequalities in (3.1) collapse to equalities. In this case, (3.1) can be
paraphrased as
$\mathbb{P}_{\theta}(\widehat{\theta}_{T}\in\mathcal{D})=e^{-ra_{T}+o(a_{T})}$,
where $r=\inf_{\theta^{\prime}\in\mathcal{D}}I(\theta^{\prime},\theta)$
represents the $I$-distance between the system matrix $\theta$ and the set
$\mathcal{D}$ of estimator realizations. Thus, $r$ represents the decay rate
of the probability $\mathbb{P}_{\theta}(\widehat{\theta}_{T}\in\mathcal{D})$,
while $\\{a_{T}\\}_{T\in\mathbb{N}}$ can be viewed as the speed of
convergence. The condition $1\ll a_{T}\ll T$ is satisfied, for example, if
$a_{T}=\sqrt{T}$, $T\in\mathbb{N}$. However, many other choices are possible.
It is perhaps surprising that if a sequence of estimators satisfies a moderate
deviations principle, then the choice of the speed
$\\{a_{T}\\}_{T\in\mathbb{N}}$ has no impact on the decay rate $r$ but may
only influence the coefficients of the higher-order terms hidden in
$o(a_{T})$. We also remark that if the inequalities in (3.1) hold for
$a_{T}=T$, $T\in\mathbb{N}$ (in which case the speed of convergence violates
the condition $1\ll a_{T}\ll T$), then
$\\{\widehat{\theta}_{T}\\}_{T\in\mathbb{N}}$ is said to satisfy a large
deviation principle [19]. It is also customary to talk about a moderate
deviation principle as being a large deviation principle with reduced speed
$\\{a_{T}\\}_{T\in\mathbb{N}}$ such that $1\ll a_{T}\ll T$.

We now show that the transformed least squares estimators

$\widehat{\vartheta}_{T}=\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)+\theta$
(3.2)

satisfy a moderate deviation principle, where the discrepancy function (1.4)
plays the role of the rate function. This result relies on another standard
regularity condition.

######  Assumption 3.3 (Light-tailed noise and stationarity).

The following hold for every $\theta\in\Theta$.

  1. (i)

The disturbances $\\{w_{t}\\}_{t\in\mathbb{N}}$ are light-tailed, i.e., there
exists $\alpha>0$ with $\mathbb{E}_{\theta}[e^{\alpha\|w_{t}\|^{2}}]<\infty$
for all $t\in\mathbb{N}$.

  2. (ii)

The initial distribution $\nu$ coincides with the invariant distribution
$\nu_{\theta}$ of the linear system (1.1).

Assumption 3.3 (i) essentially requires the tails of the noise to have no
heavier tails than a normal distribution, while Assumption 3.3 (ii) stipulates
that the linear system is in the stationary regime already at time $t=0$.

######  Proposition 3.4 (Moderate deviation principle).

If Assumptions 2.1 and 3.3 hold, $\\{\widehat{\theta}_{T}\\}_{T\in\mathbb{N}}$
denote the least squares estimators defined in (1.2) and
$\\{a_{T}\\}_{T\in\mathbb{N}}$ is a real sequence with $1\ll a_{T}\ll T$, then
the transformed least squares estimators
$\\{\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)+\theta\\}_{T\in\mathbb{N}}$
satisfy a moderate deviation principle with rate function (1.4).

Unlike the standard least squares estimators (1.2), the transformed estimators
of Proposition 3.4 depend on the unknown parameter $\theta$. However, as we
will explain below, they are useful for theoretical considerations.
Proposition 3.4 can be viewed as a corollary of [68], which uses ideas from
[67] to show that the transformed least squares estimators satisfy a moderate
deviation principle with a rate function that is defined implicitly in
variational form. Proposition 3.4 shows that this rate function admits the
explicit representation (1.4) and thus confirms (1.6). In addition, it also
relaxes the restrictive condition $\|\theta\|_{2}<1$ from [68, Proposition
2.2] to $\rho(\theta)<1$.

By identifying the discrepancy function (1.4) with the rate function of a
moderate deviation principle, Proposition 3.4 justifies our terminology,
whereby $\mathcal{P}(\theta^{\prime})$ is called the reverse $I$-projection of
$\theta^{\prime}$. Indeed, [16] use this term to denote any projection with
respect to an information divergence $I(\theta^{\prime},\theta)$. Note that
swapping the arguments $\theta^{\prime}$ and $\theta$ of the (asymmetric)
function $I(\theta^{\prime},\theta)$ would give rise to an ordinary
$I$-projection [18]. Proposition 3.4 also suggests that the reverse
$I$-projection is intimately related to maximum likelihood estimation, as
already alluded to in the introduction. Indeed, for i.i.d. training data it is
well-known that every maximum likelihood estimator can be regarded as a
reverse $I$-projection with respect to the rate function of some large
deviation principle [17, Lemma 3.1].

The power of Proposition 3.4 lies in its generality. Indeed, a moderate
deviation principle provides tight bounds on the probability of any Borel set
of estimator realizations. A simple direct application of the moderate
deviation principle established in Proposition 3.4 is described below.

######  Example 3.5 (System identification).

Consider a scalar system with $S_{w}=1$ that satisfies Assumptions 2.1 and
3.3. In this case $\Theta=(-1,1)$ with the rate function (1.4) reducing to
$I(\theta^{\prime},\theta)=\frac{1}{2}(\theta^{\prime}-\theta)^{2}/(1-\theta^{2})$.
Using the least squares estimators (1.2) to identify $\theta$, Proposition 3.4
reveals that

$\displaystyle\mathbb{P}_{\theta}(|\widehat{\theta}_{T}-\theta|\\!>\varepsilon\sqrt{a_{T}/T})$
$\displaystyle=\mathbb{P}_{\theta}(\theta+\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)\in\mathcal{D})$
$\displaystyle=\textstyle\exp\left(-\inf_{\theta^{\prime}\in\mathcal{D}}I(\theta^{\prime},\theta)\cdot
a_{T}+o(a_{T})\right)$
$\displaystyle=\textstyle\exp\left(-\frac{1}{2}\varepsilon^{2}\,a_{T}/(1-\theta^{2})+o(a_{T})\right)$

for any $\varepsilon>0$ and $T\in\mathbb{N}$, where
$\mathcal{D}=\\{\theta^{\prime}\in\mathbb{R}:|\theta^{\prime}-\theta|>\varepsilon\\}$.
This result confirms the counterintuitive insight of [54] whereby stable
systems with $|\theta|\approx 1$ are easier to identify than systems with
$|\theta|\approx 0$.

The moderate deviation principle established in Proposition 3.4 also enables
us to find statistically optimal data-driven decisions for stochastic
optimization problems, where the underlying probability measure is only
indirectly observable through finitely many training samples. Indeed, [65] and
[56] show that such optimal decisions can be found by solving data-driven
distributionally robust optimization problems.

Next, we establish several structural properties of the rate function (1.4).

######  Proposition 3.6 (Properties of $I(\theta^{\prime},\theta)$).

The rate function $I(\theta^{\prime},\theta)$ defined in (1.4) has the
following properties.

  1. (i)

$I(\theta^{\prime},\theta)$ is analytic in
$(\theta^{\prime},\theta)\in\Theta^{\prime}\times\Theta$.

  2. (ii)

If $\theta^{\prime}\in\Theta$, then the sublevel set
$\\{\theta\in\Theta:I(\theta^{\prime},\theta)\leq r\\}$ is compact for every
$r\geq 0$.

  3. (iii)

If $\theta^{\prime}\in\Theta$, then $I(\theta^{\prime},\theta)$ tends to
infinity as $\theta$ approaches the boundary of $\Theta$.

Proposition 3.6 (ii) guarantees that the minimum in (1.5) is indeed attained
and that the reverse $I$-projection is well-defined. To close this section, we
present a useful relation between the rate function $I$ and the operator norm.

######  Lemma 3.7 (Pinsker-type inequality).

For any $\theta^{\prime}\in\Theta^{\prime}$ and $\theta\in\Theta$ we have
$\|\theta^{\prime}-\theta\|_{2}^{2}\leq{2\kappa(S_{w})\cdot
I(\theta^{\prime},\theta)}$.

Lemma 3.7 provides a direct link between the nearest stable matrix problem
(1.3) and the reverse $I$-projection (1.5) because

$\inf_{{\theta}\in\Theta}\|\theta^{\prime}-{\theta}\|_{2}^{2}\leq
2\kappa(S_{w})\cdot I(\theta^{\prime},\mathcal{P}(\theta^{\prime})).$

###  3.2 Statistics of the reverse $I$-projection

In the following we apply the reverse $I$-projection to the least squares
estimator $\widehat{\theta}_{T}$ and not to the transformed least squares
estimator $\widehat{\vartheta}_{T}$ defined in (3.2) (which is anyway
unaccessible because $\theta$ is unknown), even though Proposition 3.4 relates
$I$ to $\widehat{\vartheta}_{T}$. However, an elementary calculation shows
that
$I(\widehat{\theta}_{T},\theta)=(a_{T}/T)I(\widehat{\vartheta}_{T},\theta)$,
and thus $I(\widehat{\theta}_{T},\theta)$ inherits any statistical
interpretations from $I(\widehat{\vartheta}_{T},\theta)$.

Using the results of Section 3.1, we are now ready to establish probabilistic
bounds on the operator norm distance between the projected least squares
estimator $\mathcal{P}(\widehat{\theta}_{T})$ and and the unknown true system
matrix $\theta$ with respect to the data-generating probability measure
$\mathbb{P}_{\theta}$. These bounds are random because they depend on the
$I$-distance between $\widehat{\theta}_{T}$ and
$\mathcal{P}(\widehat{\theta}_{T})$, and thus they have an a posteriori
flavour. Specifically, the following proposition provides both asymptotic as
well as finite-sample bounds, which display different strengths and weaknesses
and are both structurally identical to existing (and in some cases optimal)
finite-sample bounds for $\widehat{\theta}_{T}$; see, e.g., [52, § 6].

######  Proposition 3.8 (Statistical guarantees).

Suppose that Assumption 2.1 holds and that $\widehat{\theta}_{T}$ and
$\mathcal{P}(\widehat{\theta}_{T})$ represent the least squares estimator and
its reverse $I$-projection, respectively. Set $\varepsilon=[\kappa(S_{w})\cdot
I(\widehat{\theta}_{T},\mathcal{P}(\widehat{\theta}_{T}))]^{1/2}$, and fix any
$\theta\in\Theta$ and $\beta\in(0,1)$ . Then, we have
$\|\widehat{\theta}_{T}-\mathcal{P}(\widehat{\theta}_{T})\|_{2}\leq\varepsilon$
$\mathbb{P}_{\theta}$-almost surely, and the following statistical guarantees
hold.

  1. (i)

There is a constant $C(n)$ that depends only on $n$ such that for all
$T\geq\kappa(S_{w}){\log(1/\beta)}/({\varepsilon^{2}C(n)})$ we have

$\mathbb{P}_{\theta}\left(\|\theta-\mathcal{P}(\widehat{\theta}_{T})\|_{2}\leq
2\varepsilon\right)\geq 1-\beta.$ (3.3a)

  2. (ii)

If Assumption 3.3 holds and $\\{a_{T}\\}_{T\in\mathbb{N}}$ is a real sequence
satisfying $1\ll a_{T}\ll T$, then for all $T$ with $a_{T}\geq
2\kappa(S_{w})\big{(}\log(1/\beta)+o(a_{T})\big{)}/\varepsilon^{2}$ we have

$\mathbb{P}_{\theta}\left(\|\theta-\mathcal{P}(\widehat{\theta}_{T})\|_{2}\leq\varepsilon\left(\sqrt{{a_{T}}/{T}}+1\right)\right)\geq
1-\beta.$ (3.3b)

As $\widehat{\theta}_{T}$ and $\mathcal{P}(\widehat{\theta}_{T})$ both
converge $\mathbb{P}_{\theta}$-almost surely to $\theta\in\Theta$, and as $I$
is continuous in both of its arguments, it is easy to show that the random
bound $\varepsilon$ as defined in Proposition 3.8 converges
$\mathbb{P}_{\theta}$-almost surely to $0$ as $T$ grows. Therefore, the bounds
(3.3a) and (3.3b) improve with $T$. We further emphasize that the finite-
sample bound (3.3a) and the asymptotic bound (3.3b) are remarkably similar.
Yet the bound (3.3b), which follows from the moderate deviation principle
established in Proposition 3.4, involves no unknown constants. As the
inequalities in (3.1) are asymptotically tight, we conjecture that (3.3b) is
the best possible asymptotic bound.

###  3.3 Computation of the reverse $I$-projection

We now address the numerical computation of $\mathcal{P}(\theta^{\prime})$ as
defined in (1.5) for any given estimator realization
$\theta^{\prime}\in\Theta^{\prime}$. To this end, we fix $Q\succ 0$ and show
that solving (1.5) is equivalent to finding a minimizer of the optimization
problem

$\min_{\theta\in\Theta}\\{\mathsf{tr}(QS_{\theta}):I(\theta^{\prime},\theta)\leq
r\\}$ (3.4)

for the smallest radius $r=\underline{r}$ that renders (3.4) feasible. Note
that $\underline{r}$ exists because the optimal value of (3.4) is lower semi-
continuous in $r$. In addition, problem (3.4) admits a minimizer for any
$r\geq\underline{r}$ due to Proposition 3.6 (ii). The proposed procedure works
because if $\theta^{\prime}$ is unstable, then any $\theta$ feasible in (3.4)
is stable, and its $I$-distance to $\theta^{\prime}$ is at most $r$. Setting
$r=\underline{r}$ thus ensures that any minimizer of (3.4) is a reverse
$I$-projection of $\theta^{\prime}$ and that
$I(\theta^{\prime},\mathcal{P}(\theta^{\prime}))=\underline{r}$. Moreover, the
proposed procedure is computationally attractive because we will prove below
that (3.4) is equivalent to a standard LQR problem. We emphasize that the
exact choice of $Q$ has no effect on the validity and hardly any effect on the
numerical performance of this procedure.

######  Proposition 3.9 (Reformulation of (1.5)).

If $\theta^{\prime}\in\Theta^{\prime}$, $Q\succ 0$ and $\underline{r}$ is the
smallest $r\geq 0$ for which (3.4) is feasible, then any minimizer of (3.4) at
$r=\underline{r}$ is a reverse $I$-projection.

Note that if $r\geq\overline{r}=I(\theta^{\prime},0)$, then problem (3.4) has
the trivial solution $\theta=0$, and its optimal value reduces to
$\mathsf{tr}(QS_{w})$.111One readily verifies that
$\overline{r}=\frac{1}{2}\|S_{w}^{-1/2}\theta^{\prime}S_{w}^{1/2}\|_{F}^{2}$,
where $\|\cdot\|_{F}$ stands for the Frobenius norm. In this case, the rate
constraint is not binding at optimality. If $r<\overline{r}$, on the other
hand, then problem (3.4) is infeasible for $r<\underline{r}$ and admits a
quasi-closed form solution for $r>\underline{r}$ as explained in the following
proposition.

######  Proposition 3.10 (Optimal solution of (3.4)).

Suppose that Assumption 2.1 holds. Then, for every
$\theta^{\prime}\notin\Theta$ there exists an analytic function
$\varphi:(\underline{r},\overline{r})\rightarrow(0,\infty)$ that is increasing
and bijective such that the following hold for all
$r\in(\underline{r},\overline{r})$.

  1. (i)

For any $\delta\in(0,\infty)$ the matrix $P_{\delta}\in\mathcal{S}^{n}$ is the
unique positive definite solution of the Riccati equation

$\displaystyle P=$
$\displaystyle\;Q+{\theta^{\prime}}^{\mathsf{T}}P\left(I_{n}+2\delta
S_{w}P\right)^{-1}\theta^{\prime}.$ (3.5)

  2. (ii)

The matrix $\theta_{\delta}^{\star}=(I_{n}+2\delta
S_{w}P_{\delta})^{-1}\theta^{\prime}$ is the unique solution of problem (3.4)
at $r=\varphi^{-1}(\delta)$, and the rate constraint is binding at optimality,
i.e., $I(\theta^{\prime},\theta^{\star}_{\delta})=r$.

We have seen that evaluating $\mathcal{P}(\theta^{\prime})$ is equivalent to
solving (3.4) at $r=\underline{r}$. Unfortunately, $\underline{r}$ is unknown,
and Proposition 3.10 only characterizes solutions of (3.4) for
$r>\underline{r}$. However, by the properties of $\varphi$ established in
Proposition 3.10, we also have $\lim_{r\downarrow\underline{r}}\varphi(r)=0$,
which is equivalent to $\lim_{\delta\downarrow
0}\varphi^{-1}(\delta)=\underline{r}$. A standard continuity argument
therefore implies that $\lim_{\delta\downarrow 0}\theta^{\star}_{\delta}$
solves (3.4) at $r=\underline{r}$. In practice, we may simply set $\delta$ to
a small positive number and compute $\theta^{\star}_{\delta}$ by solving (3.5)
to find a high-accuracy approximation for the reverse $I$-projection
$\mathcal{P}(\theta^{\prime})$.

######  Proposition 3.11 (Computing the reverse $I$-projection).

If Assumption 2.1 holds, $\theta^{\prime}\notin\Theta$ and $Q\succ 0$, then
there exists a $p\geq 1$ such that for all $\delta>0$ the matrix
$\theta^{\star}_{\delta}$ from Proposition 3.10 (ii) is stable and satisfies

$\|\mathcal{P}(\theta^{\prime})-\theta^{\star}_{\delta}\|_{2}=O(\delta^{p}).$
(3.6)

In addition, $\theta^{\star}_{\delta}$ can be computed as

$\theta^{\star}_{\delta}=\theta^{\prime}+\mathsf{dlqr}(\theta^{\prime},I_{n},Q,(2\delta
S_{w})^{-1}),$

where the standard LQR routine222For example, see
https://juliacontrol.github.io/ControlSystems.jl/latest/examples/example/#LQR-
design $\mathsf{dlqr}(\cdot)$ has time and memory complexity of the order
$O(n^{3})$ and $O(n^{2})$, respectively.

######  Corollary 3.12 (The (numerical) reverse $I$-projection preserves the
structure of $\theta^{\prime}$).

For any $\theta^{\prime}\in\Theta^{\prime}$ there exist invertible matrices
$\Lambda,\Lambda_{\delta}\in\mathbb{R}^{n\times n}$ such that
$\mathcal{P}(\theta^{\prime})=\Lambda^{-1}\theta^{\prime}$ and
$\theta^{\star}_{\delta}=\Lambda^{-1}_{\delta}\theta^{\prime}$.

Corollary 3.12 implies, among other things, that the reverse $I$-projection
preserves the kernel of $\theta^{\prime}$, see [29] for more information. [24]
have shown that if $\theta^{\prime}=\alpha\mathbbold{1}_{n\times n}$ with
$\alpha\in[\frac{1}{n},\frac{2}{n}]$, then the solution of the closest stable
matrix problem (1.3) with respect to the Frobenius norm is
$\Pi_{\Theta}(\theta^{\prime})=\frac{1}{n}\mathbbold{1}_{n\times n}$, which
lies on the boundary of $\Theta$. A simple calculation further shows that
$\mathcal{P}(\theta^{\prime})=\frac{1}{2n}\mathbbold{1}_{n\times n}$, which
lies in the interior of $\Theta$ and has the same structure as
$\theta^{\prime}$ and $\Pi_{\Theta}(\theta^{\prime})$, thus exemplifying
Corollary 3.12. For $\alpha>\frac{2}{n}$ problem (1.3) appears to have many
local minima, while the reverse $I$-projection remains unique as well as
structurally similar to $\theta^{\prime}$.

(a) Approximation in operator norm for $m=1$.

(b) Eigenvalue spectra for $m=1$.

(c) Historgrams of the spectral radii for $m=1$.

(d) Historgrams of the spectral radii for $m=9$.

Figure 3.1: Comparison of the least squares estimator $\widehat{\theta}_{T}$
and its reverse $I$-projection $\mathcal{P}(\widehat{\theta}_{T})$ against the
CG and FG estimators $\widehat{\theta}_{\rm CG}$ and $\widehat{\theta}_{\rm
FG}$, respectively, based on 250 independent simulation runs with
$\widehat{\theta}_{T}\notin\Theta$. To aid visibility, Figure 1(a) only shows
the first 50 experiments.

(a) Convergence in operator norm for $n=1$.

(b) Convergence in operator norm for $n=10$.

(c) Convergence of spectral radii for $n=1$.

(d) Convergence of spectral radii for $n=10$.

Figure 3.2: Convergence behavior of $\widehat{\theta}_{T}$ and
$\mathcal{P}(\widehat{\theta}_{T})$. Solid lines represent averages and shaded
areas represent ranges across $10^{4}$ simulations (corresponding to 100
randomly generated system matrices $\theta$ and 100 randomly generated state
trajectories per system matrix).

Now we have all the tools in place to prove Theorem 2.2.

###### Proof of Theorem 2.2.

As for Assertion (i), recall from [14] that
$\lim_{T\to\infty}\widehat{\theta}_{T}=\theta$ $\mathbb{P}_{\theta}$-almost
surely. Therefore, we have $\mathbb{P}_{\theta}$-almost surely that

$\displaystyle\textstyle\lim_{T\to\infty}\mathcal{P}(\widehat{\theta}_{T})$
$\displaystyle=\textstyle\lim_{T\to\infty}\arg\min_{\bar{\theta}\in\Theta}I(\widehat{\theta}_{T},\bar{\theta})$
$\displaystyle=\arg\min_{\bar{\theta}\in\Theta}\textstyle\lim_{T\to\infty}I(\widehat{\theta}_{T},\bar{\theta})$
$\displaystyle=\arg\min_{\bar{\theta}\in\Theta}I(\textstyle\lim_{T\to\infty}\widehat{\theta}_{T},\bar{\theta})$
$\displaystyle=\arg\min_{\bar{\theta}\in\Theta}I(\theta,\bar{\theta})$
$\displaystyle=\theta,$

where the first equality exploits the definition of
$\mathcal{P}(\widehat{\theta}_{T})$ in (1.5). The second equality follows from
the strict convexity of the rate function in its first argument and [55,
Theorem 9.17], which imply that the reverse $I$-projection is continuous. The
third equality follows from the continuity of the rate function established in
Proposition 3.6 (i), and the last equality holds because the rate function
vanishes if and only if its arguments coincide. This proves Assertion (i).
Assertion (ii) follows from Proposition 3.8, and Assertion (iii) from
Proposition 3.11. ∎

##  4 Numerical simulations

The two subsequent examples showcase the statistical and computational
properties of the reverse $I$-projection.

######  Example 4.1 (Spectral approximation quality).

Select $m\in\\{1,9,64\\}$, and set $Q=I_{3m}$,
$w_{t}\sim\mathcal{N}(0,I_{3m})$ and

$\theta=(Y\otimes I_{m}),\quad\text{where}\quad Y=\begin{bmatrix}0.95&0.1&1\\\
-0.1&0.95&0\\\ 0&0&0.9\end{bmatrix}.$

The eigenvalues of $Y$ are $0.9$ and $0.95\pm i\cdot 0.1$, where $i$ is the
imaginary unit, and thus $\theta$ is almost unstable. Set now $T=25\sqrt{m}$,
and generate 250 independent state trajectories for which
$\widehat{\theta}_{T}\notin\Theta$. Next, compute
$\mathcal{P}(\widehat{\theta}_{T})$ approximately as described in Proposition
3.11 for $\delta=10^{-9}$. In addition, compute
$\Pi_{\Theta}(\widehat{\theta}_{T})$ with respect to the Frobenius norm by
using the approximate constraint generation (CG) method of [13] and the exact
fast gradient (FG) method of [21].333sites.google.com/site/nicolasgillis/code
Figure 1(a) shows that for $m=1$ all methods succeed in approximating $\theta$
reasonably closely, with the FG method having a slight edge. However, from
Figure 1(b) it becomes apparent that $\mathcal{P}(\widehat{\theta}_{T})$
approximates the eigenvalue spectrum of $\theta$ best. All of its eigenvalues
reside within the complex unit circle and concentrate near the true spectrum
of $\theta$, which might be explained by the structure-preserving property of
the reverse $I$-projection established in Corollary 3.12. In contrast, the CG
method often produces unstable estimators, and the FG method generates
estimators that reside on the boundary of $\Theta$. These observations are
consistent with Figure 1(c), which displays the empirical distribution of the
spectral radii corresponding to the different estimators. Indeed, the
histogram corresponding to the reverse $I$-projection is confined to $[0,1]$
and centered around $\rho(\theta)$. The FG method, on the other hand, is
designed to generate estimators with unit spectral radius, which could be
undesirable in applications. In higher dimensions it is less likely that
$\widehat{\theta}_{T}\notin\Theta$, but if this happens, then
$\widehat{\theta}_{T}$ is likely to be close to the boundary of $\Theta$, see
Figure 1(d). Thus, $\mathcal{P}(\widehat{\theta}_{T})$ is more likely to have
a spectral radius close to $1$. This phenomenon is further accentuated for
$m=64$, see Figure 1(d) in the appendix.

As both the reverse $I$-projection and the FG method have complexity
$O(n^{3})$, we compare their runtimes for $\theta^{\prime}=(Y\otimes
2I_{m})\notin\Theta$ as a function of $n=3m$, see Figure 1(c) in the appendix.
We observe that the reverse $I$-projection is faster for $n\lesssim 500$,
while the FG method dominates for higher dimensions. We remark that the
reverse $I$-projection is computed using off-the-shelf software but could be
sped up by using dedicated large-scale algorithms [20, 4]. All simulations
were implemented in Julia [11] and run on a 4GHz CPU with 16Gb RAM.

######  Example 4.2 (Statistical guarantees).

The second experiment is designed to validate the statistical guarantees of
Proposition 3.8. To this end, choose $n\in\\{1,10,100\\}$, and sample 100
stable matrices from a standard normal distribution on $\mathbb{R}^{n\times
n}$ restricted to $\Theta$. For each such matrix $\theta$, generate 100 state
trajectories of length $\overline{T}=10^{2}(n+1)$, and compute
$\mathcal{P}(\widehat{\theta}_{T})$ for every $T=1,\ldots,\overline{T}$
approximately as described in Proposition 3.11 for $\delta=10^{-9}$. Figures
2(a) and 2(b) visualize the convergence of the estimators
$\widehat{\theta}_{T}$ and $\mathcal{P}(\widehat{\theta}_{T})$ to $\theta$
with respect to the operator norm for $n=1$ and $n=10$, respectively. Both
figures are consistent with the $1/\sqrt{T}$ scaling law anticipated by
Proposition 3.8 (ii). Although Example 4.1 revealed that
$\rho(\mathcal{P}(\widehat{\theta}_{T}))$ can concentrate away from
$\rho(\theta)$ in high dimensions, Figures 2(c) and 2(d) show that
$\rho(\mathcal{P}(\widehat{\theta}_{T}))$ converges to $\rho(\theta)$ on
average. Figures 1(a) and 1(b) in the appendix extend these results to
$n=100$. Note that the finite sample bound (3.3a) is more difficult to
validate numerically due to the unknown constant $C(n)$. Although asymptotic,
the bound (3.3b) derived from the moderate deviation principle of Proposition
3.4 is fully transparent.

##  5 Appendix

In this appendix we collect all proofs not contained in the main body of the
paper plus we provide some auxiliary results.

###  5.1 Proofs of Section 3.1

###### Proof of Proposition 3.4.

Fix any $\theta\in\mathbb{R}^{n\times n}$ and assume that $\|\theta\|_{2}<1$.
This condition is stronger than Assumption 2.1 (i) because the spectral radius
$\rho(\theta)$ is bounded above by the spectral norm $\|\theta\|_{2}$.
Together with Assumptions 2.1 (ii) and 3.3, this condition implies via [68,
Proposition 2.2], that the transformed least squares estimators
$\\{\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)+\theta\\}_{T\in\mathbb{N}}$
satisfy a moderate deviation principle with rate function

$\displaystyle I(\theta^{\prime},\theta)=\sup_{L\in\mathbb{R}^{n\times
n}}\left\\{\langle
L,\theta^{\prime}-\theta\rangle-\frac{1}{2}\mathbb{E}_{\theta}\left[{\langle
L,w_{1}x_{0}^{\mathsf{T}}S_{\theta}^{-1}}\rangle^{2}\right]\right\\},$ (5.1)

where the inner product of two matrices $A,B\in\mathbb{R}^{n\times n}$ is
defined as $\langle A,B\rangle=\mathsf{tr}(A^{\mathsf{T}}B)$. This rate
function, which is defined implicitly as the optimal value of an optimization
problem, captures the speed at which the transformed least squares estimators
(and indirectly also the standard least squares estimators) converge to
$\theta$.

Next, we will demonstrate that (5.1) is equivalent to (1.4). As a preparation,
we derive the analytical solution of the following unconstrained convex
quadratic maximization problem over the matrix space $\mathbb{R}^{n\times n}$,

$\max_{X\in\mathbb{R}^{n\times n}}\left\\{\langle
C,X\rangle-\frac{1}{2}\mathsf{tr}\left(XB_{1}X^{\mathsf{T}}B_{2}\right)\right\\},$
(5.2)

which is parameterized by $B_{1},B_{2}\in\mathcal{S}^{n}_{\succ 0}$ and
$C\in\mathbb{R}^{n\times n}$. As the trace term
$\mathsf{tr}(XB_{1}X^{\mathsf{T}}B_{2})$ is convex in $X$ by virtue of [39,
Corollary 1.1], we can solve (5.2) by setting the gradient of the objective
function to zero. Specifically, using [6, Propositions 10.7.2 & 10.7.4], we
find

$\displaystyle{\nabla_{X}}\left(\langle{C},{X}\rangle-\frac{1}{2}\mathsf{tr}(XB_{1}X^{\mathsf{T}}B_{2})\right)=C^{\mathsf{T}}-B_{1}X^{\mathsf{T}}B_{2}.$

As $B_{1},B_{2}\succ 0$, it is easy to verify that this gradient vanishes at
$X^{\star}=B_{2}^{-1}CB_{1}^{-1}$, which implies that the optimal value of
problem (5.2) amounts to
$\frac{1}{2}\mathsf{tr}(B_{2}^{-1}CB_{1}^{-1}C^{\mathsf{T}})$.

We next rewrite the expectation term in (5.1) as

$\displaystyle\mathbb{E}_{\theta}\left[{\langle{L},{w_{1}x_{0}^{\mathsf{T}}S_{\theta}^{-1}}}\rangle^{2}\right]$
$\displaystyle=\mathbb{E}_{\theta}\left[{(w_{1}^{\mathsf{T}}LS_{\theta}^{-1}x_{0})^{2}}\right]$
$\displaystyle=\mathbb{E}_{\theta}\left[{w_{1}^{\mathsf{T}}LS_{\theta}^{-1}S_{\theta}S_{\theta}^{-1}L^{\mathsf{T}}w_{1}}\right]$
$\displaystyle=\mathbb{E}_{\theta}{\mathsf{tr}({LS_{\theta}^{-1}L^{\mathsf{T}}w_{1}w_{1}^{\mathsf{T}}}})$
$\displaystyle=\mathsf{tr}({LS_{\theta}^{-1}L^{\mathsf{T}}S_{w}}),$

where the second equality follows from the independence of $x_{0}$ and $w_{1}$
and from Assumption 3.3 (ii), which implies that the initial state $x_{0}$ is
governed by the invariant state distribution $\nu_{\theta}$ and thus has zero
mean and covariance matrix $S_{\theta}$. Substituting the resulting trace term
into (5.1) yields

$\displaystyle I(\theta^{\prime},\theta)$
$\displaystyle=\max_{L\in\mathbb{R}^{d\times
d}}\left\\{\langle{\theta^{\prime}-\theta},{L}\rangle-\frac{1}{2}\mathsf{tr}\left({LS_{\theta}^{-1}L^{\mathsf{T}}S_{w}}\right)\right\\}$
$\displaystyle=\frac{1}{2}\mathsf{tr}\left({S_{w}^{-1}(\theta^{\prime}-\theta)S_{\theta}(\theta^{\prime}-\theta)^{\mathsf{T}}}\right),$

where the second equality follows from our analytical solution of problem
(5.2) in the special case where $B_{1}=S_{\theta}^{-1}$, $B_{2}=S_{2}$ and
$C=\theta^{\prime}-\theta$. Thus, the rate function coincides indeed with the
discrepancy function (1.4).

At last, we show that the moderate deviations principle established for
$\|\theta\|_{2}<1$ remains valid for all asymptotically stable system
matrices. To this end, fix any $\theta$ with $\rho(\theta)<1$. By standard
Lyapunov stability theory, e.g., Theorem 5.3.5 from [37], there exists a
$P\in\mathcal{S}^{n}_{\succ 0}$ with $P-\theta^{\mathsf{T}}P\theta\succ 0$.
Using the matrix $P$, we can apply the change of variables
$\bar{x}_{t}=P^{\frac{1}{2}}x_{t}$ and $\bar{w}_{t}=P^{\frac{1}{2}}w_{t}$ to
obtain the auxiliary linear dynamical system

$\bar{x}_{t+1}=\bar{\theta}\,\bar{x}_{t}+\bar{w}_{t},\quad\bar{x}_{0}\sim\bar{\nu},$

with $\bar{\theta}=P^{\frac{1}{2}}\theta P^{-\frac{1}{2}}$, where the noise
$\bar{w}_{t}$ has mean zero and covariance matrix
$S_{\bar{w}}=P^{\frac{1}{2}}S_{w}P^{\frac{1}{2}}$ for all $t\in\mathbb{N}$,
and $\bar{\nu}=\nu\circ P^{-\frac{1}{2}}$ is the pushforward distribution of
$\nu$ under the coordinate transformation $P^{\frac{1}{2}}$. Note also that
the invariante state covariance matrix is given by
$S_{\bar{\theta}}=P^{\frac{1}{2}}S_{\theta}P^{\frac{1}{2}}$. By construction,
the auxiliary linear system is equivalent to (1.1) and satisfies Assumptions
2.1 (ii) and 3.3. Moreover, multiplying $P-\theta^{\mathsf{T}}P\theta\succ 0$
from both sides with $P^{-\frac{1}{2}}$ yields
$I_{n}-\bar{\theta}^{\mathsf{T}}\bar{\theta}\succ 0$, which means that the
largest eigenvalue of $\bar{\theta}^{\mathsf{T}}\bar{\theta}$ is strictly
smaller than 1 or, equivalently, that $\|\bar{\theta}\|_{2}<1$. If we denote
by $\widehat{\bar{\theta}}_{T}$ the least squares estimator for $\bar{\theta}$
based on $T$ state observations of the auxiliary linear system, we may then
conclude from the first part of the proof that the estimators
$\\{\sqrt{T/a_{T}}(\widehat{\bar{\theta}}_{T}-\theta)+\theta\\}_{T\in\mathbb{N}}$
satisfy a moderate deviations principle with rate function

$\bar{I}(\bar{\theta}^{\prime},\bar{\theta})=\mathrm{tr}\left(S_{\bar{w}}^{-1}(\bar{\theta}^{\prime}-\bar{\theta})S_{\bar{\theta}}(\bar{\theta}^{\prime}-\bar{\theta})^{\mathsf{T}}\right).$

One also readily verifies from (1.2) that the least squares estimators
pertaining to the original and the auxiliary linear systems are related
through the continuous transformation
$\widehat{\bar{\theta}}_{T}=P^{\frac{1}{2}}\widehat{\theta}_{T}P^{-\frac{1}{2}}$.
The corresponding transformed estimators evidently obey the same relation. By
the contraction principle [19, Theorem 4.2.1], the estimators
$\\{\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)+\theta\\}_{T\in\mathbb{N}}$
thus satisfy a moderate deviations principle with rate function

$\displaystyle\bar{I}(P^{-\frac{1}{2}}\theta^{\prime}P^{\frac{1}{2}},P^{-\frac{1}{2}}\theta
P^{\frac{1}{2}})=I(\theta^{\prime},\theta).$

This observation completes the proof. ∎

To prove that $I(\theta^{\prime},\theta)$ is analytic, we recall that the
stationary state covariance matrix $S_{\theta}$ is analytic in
$\theta\in\Theta$.

######  Lemma 5.1 (Analyticity of $S_{\theta}$ [49], Lemma 3.2).

If $S_{w}\succ 0$, then the solution $S_{\theta}$ to the Lyapunov equation
(2.1) is analytic in $\theta\in\Theta$.

###### Proof of Proposition 3.6.

As for Assertion (i), note that the quadratic function
$(\theta^{\prime}-\theta)^{\mathsf{T}}(\theta^{\prime}-\theta)S_{w}^{-1}$ is
manifestly analytic in $(\theta^{\prime},\theta)$. Moreover, the stationary
state covariance matrix $S_{\theta}$ is analytic in $\theta$ (and thus also in
$(\theta^{\prime},\theta)$) by virtue of Lemma 5.1. The rate function
$I(\theta^{\prime},\theta)$ defined in (1.4) can therefore be viewed as an
inner product of two matrix-valued analytic functions and is thus analytic
thanks to [33, Proposition 2.2.2].

The proof of Assertion (ii) consists of two steps. We first prove that if a
sequence $\\{\theta_{k}\\}_{k\in\mathbb{N}}$ in $\Theta$ has an unstable limit
$\theta$ (i.e., $\rho(\theta)=1$), then there exists a subsequence
$\\{\theta_{k_{l}}\\}_{l\in\mathbb{N}}$ with
$\lim_{l\to\infty}I(\theta^{\prime},\theta_{k_{l}})=\infty$ for all
$\theta^{\prime}\in\Theta$ (Step 1). We then use this result to show that the
set $\\{\theta\in\Theta:I(\theta^{\prime},\theta)\leq r\\}$ is compact for all
$r\geq 0$ (Step 2).

Step 1: We first derive an easily computable lower bound on the rate function
$I(\theta^{\prime},\theta)$ for any asymptotically stable matrices
$\theta^{\prime},\theta\in\Theta$. To this end, we denote by
$\lambda\in\mathbb{C}$ an eigenvalue of $\theta$ whose modulus $|\lambda|$
matches the spectral radius $\rho(\theta)<1$. We further denote by
$v\in\mathbb{C}^{n}$ a normalized eigenvector corresponding to the eigenvalue
$\lambda$, that is, $\|v\|=1$ and $\theta v=\lambda v$. We also use
$\beta=\lambda_{\min}(S_{w})/\lambda_{\max}(S_{w})>0$ as a shorthand for the
inverse condition number of the noise covariance matrix $S_{w}\succ 0$.
Recalling that for any $A,B,C\in\mathcal{S}^{n}_{\succeq 0}$ the semidefinite
inequality $A\succeq B$ implies $\mathsf{tr}({AC})\geq\mathsf{tr}({BC})$, we
find the following estimate.

$\displaystyle 2I(\theta^{\prime},\theta)$
$\displaystyle=\mathsf{tr}\left({S_{w}^{-1}(\theta^{\prime}-\theta)S_{\theta}(\theta^{\prime}-\theta)^{\mathsf{T}}}\right)$
$\displaystyle\geq\lambda^{-1}_{\max}(S_{w})\,\mathsf{tr}\left({(\theta^{\prime}-\theta)S_{\theta}(\theta^{\prime}-\theta)^{\mathsf{T}}}\right)$
$\displaystyle\geq\beta\sum_{k=0}^{\infty}\mathsf{tr}\left({(\theta^{\prime}-\theta)\theta^{k}(\theta^{k})^{\mathsf{T}}(\theta^{\prime}-\theta)^{\mathsf{T}}}\right)$
$\displaystyle=\beta\sum_{k=0}^{\infty}\mathsf{tr}\left({(\theta^{k})^{\mathsf{T}}(\theta^{\prime}-\theta)^{\mathsf{T}}(\theta^{\prime}-\theta)\theta^{k}}\right)$
$\displaystyle\geq\beta\sum_{k=0}^{\infty}v^{\mathsf{H}}(\theta^{k})^{\mathsf{T}}(\theta^{\prime}-\theta)^{\mathsf{T}}(\theta^{\prime}-\theta)\theta^{k}v$
$\displaystyle=\beta\sum_{k=0}^{\infty}|\lambda|^{2k}v^{\mathsf{H}}(\theta^{\prime}-\theta)^{\mathsf{T}}(\theta^{\prime}-\theta)v$
$\displaystyle=\beta\|(\theta^{\prime}-\theta)v\|^{2}\frac{1}{1-|\lambda|^{2}}$

Here, the first equality follows from the definition of the rate function in
(1.4), and the first inequality exploits the bound
$\lambda_{\max}(S_{w})I_{n}\succeq S_{w}$. The second inequality holds due to
the series representation
$S_{\theta}=\sum_{t=0}^{\infty}\theta^{t}S_{w}(\theta^{t})^{\mathsf{T}}$, the
bound $S_{w}\succeq\lambda_{\min}(S_{w})I_{n}$ and the definition of $\beta$.
The second equality exploits the cyclicity property of the trace, and the
third inequality holds because any (real) matrix $C\in\mathcal{S}^{n}_{\succeq
0}$ satisfies

$\mathsf{tr}({C})\geq w^{\mathsf{H}}Cw\quad\forall
w\in\mathbb{C}^{n}:\;\|w\|=1.$

The third equality then uses the eigenvalue equation $\theta v=\lambda v$, and
the last equality holds because $|\lambda|=\rho(\theta)<1$. We thus conclude
that the rate function admits the lower bound

$I(\theta^{\prime},\theta)\geq\frac{\beta}{2}\|(\theta^{\prime}-\theta)v\|^{2}\frac{1}{1-|\lambda|^{2}}.$
(5.3)

Consider now a converging sequence $\\{\theta_{k}\\}_{k\in\mathbb{N}}$ in
$\Theta$ whose limit $\theta$ satisfies $\rho(\theta)=1$. Define
$\lambda_{k}\in\mathbb{C}$ as an eigenvalue of $\theta_{k}$ with
$|\lambda_{k}|=\rho(\theta_{k})<1$ and let $v_{k}\in\mathbb{C}^{n}$ be a
normalized eigenvector corresponding to $\lambda_{k}$, that is, $\|v_{k}\|=1$
and $\theta_{k}v_{k}=\lambda_{k}v_{k}$. As the spectral radius is a continuous
function, we then have

$\lim_{k\to\infty}|\lambda_{k}|=\lim_{k\to\infty}\rho(\theta_{k})=\rho(\lim_{k\to\infty}\theta_{k})=\rho(\theta)=1.$

In addition, as the unit spheres in $\mathbb{C}$ and in $\mathbb{C}^{n}$ are
both compact, there exists a subsequence
$\\{(\lambda_{k_{l}},v_{k_{l}})\\}_{l\in\mathbb{N}}$ converging to a point
$(\lambda,v)\in\mathbb{C}\times\mathbb{C}^{n}$ with $|\lambda|=1$ and
$\|v\|=1$. This limit satisfies the eigenvalue equation

$\theta
v=\lim_{l\to\infty}\theta_{k_{l}}v_{k_{l}}=\lim_{l\to\infty}\lambda_{k_{l}}v_{k_{l}}=\lambda
v,$ (5.4)

which implies that $v$ is an eigenvector of $\theta$ corresponding to the
eigenvalue $\lambda$ with $|\lambda|=1=\rho(\theta)$.

The above reasoning allows us to conclude that

$\displaystyle\lim_{l\to\infty}I(\theta^{\prime},\theta_{k_{l}})$
$\displaystyle\geq\lim_{l\to\infty}\frac{\beta}{2}\|(\theta_{k_{l}}-\theta^{\prime})v_{k_{l}}\|^{2}\frac{1}{1-|\lambda_{k_{l}}|^{2}}$
$\displaystyle=\lim_{l\to\infty}\frac{\beta}{2}\|\lambda_{k_{l}}v_{k_{l}}-\theta^{\prime}v_{k_{l}}\|^{2}\frac{1}{1-|\lambda_{k_{l}}|^{2}}$
$\displaystyle=\frac{\beta}{2}\|\lambda
v-\theta^{\prime}v\|^{2}\lim_{l\to\infty}\frac{1}{1-|\lambda_{k_{l}}|^{2}}=\infty,$

where the inequality follows from (5.3), the first equality holds because
$\theta_{k}v_{k}=\lambda_{k}v_{k}$, and the second equality exploits (5.4).
Finally, the last equality holds because $\lim_{k\to\infty}|\lambda_{k}|=1$
and because the term $\frac{\beta}{2}\|\lambda v-\theta^{\prime}v\|^{2}$ is
strictly positive. Indeed, this non-negative term can only vanish if
$\theta^{\prime}v=\lambda v$, which would imply that $\theta^{\prime}$ is
unstable (as $|\lambda|=1)$ and thus contradict the assumption that
$\theta^{\prime}\in\Theta$. This observation completes Step 1.

Step 2: Select now any $\theta^{\prime}\in\Theta$ and $r\geq 0$, and define
$\mathcal{A}=\\{\theta\in\Theta:I(\theta^{\prime},\theta)\leq r\\}$. In order
to prove that $\mathcal{A}$ is compact, we need to show that it is bounded and
closed. This is potentially difficult because $\Theta$ itself is unbounded and
open. In order to prove boundedness of $\mathcal{A}$, note that every
$\theta\in\mathcal{A}$ satisfies

$\displaystyle r\geq I(\theta^{\prime},\theta)$
$\displaystyle=\frac{1}{2}\mathsf{tr}\left({S_{w}^{-1}(\theta^{\prime}-\theta)S_{\theta}(\theta^{\prime}-\theta)^{\mathsf{T}}}\right)$
$\displaystyle\geq\frac{1}{2}\mathsf{tr}\left({S_{w}^{-1}(\theta^{\prime}-\theta)S_{w}(\theta^{\prime}-\theta)^{\mathsf{T}}}\right),$

where the second inequality follows from the trivial bound $S_{\theta}\succeq
S_{w}$, which is implied by the Lyapunov equation (2.1). Thus, the sublevel
set $\mathcal{A}$ is contained in a bounded ellipsoid,

$\mathcal{A}\subset\left\\{\theta\in\mathbb{R}^{n\times
n}:\frac{1}{2}\mathsf{tr}\left({S_{w}^{-1}(\theta^{\prime}-\theta)S_{w}(\theta^{\prime}-\theta)^{\mathsf{T}}}\right)\leq
r\right\\},$

and thus $\mathcal{A}$ is bounded. To show that $\mathcal{A}$ is closed,
consider a converging sequence $\\{\theta_{k}\\}_{k\in\mathbb{N}}$ in
$\mathcal{A}$ with limit $\theta$. We first prove that $\theta\in\Theta$.
Suppose for the sake of argument that $\theta\notin\Theta$. As $\theta$ is the
limit of a sequence in $\mathcal{A}\subset\Theta$, this implies that $\theta$
must reside on the boundary of $\Theta$ ( i.e., $\rho(\theta)=1$). By the
results of Step 1, we may thus conclude that there exists a subsequence
$\\{\theta_{k_{l}}\\}_{l\in\mathbb{N}}$ with
$\lim_{l\to\infty}I(\theta^{\prime},\theta_{k_{l}})=\infty$. Clearly, we then
have $I(\theta^{\prime},\theta_{k_{l}})>r$ for all sufficiently large $l$,
which contradicts the assumption that $\theta_{k_{l}}\in\mathcal{A}$ for all
$l\in\mathbb{N}$. Thus, our initial hypothesis was wrong, and we may conclude
that $\theta\in\Theta$. In addition, we have

$\displaystyle
r\geq\lim_{k\to\infty}I(\theta^{\prime},\theta_{k})=I(\theta^{\prime},\lim_{k\to\infty}\theta_{k})=I(\theta^{\prime},\theta),$

where the inequality holds because $\theta_{k}\in\mathcal{A}$ for all
$k\in\mathbb{N}$. Here, the first equality follows from Assertion (i), which
ensures that the rate function is analytic and thus continuous. Hence, we find
that $\theta\in\mathcal{A}$. As the sequence
$\\{\theta_{k}\\}_{k\in\mathbb{N}}$ was chosen arbitrarily, we conclude that
$\mathcal{A}$ is closed. In summary, we have shown that $\mathcal{A}$ is
bounded and closed and thus compact. This observation completes Step 2. Hence,
Assertion (ii) follows.

As for Assertion (iii), fix $\theta^{\prime}\in\Theta$ and consider a sequence
$\\{\theta_{k}\\}_{k\in\mathbb{N}}$ in $\Theta$ whose limit $\theta$ resides
on the boundary of the open set $\Theta$. This implies that
$\theta\notin\Theta$. Next, choose any $r\geq 0$. We know from Assertion (ii)
that $\mathcal{A}=\\{\theta\in\Theta:I(\theta^{\prime},\theta)\leq r\\}$ is a
compact subset of $\Theta$, and thus $\theta\notin\mathcal{A}$. Hence, the
complement of $\mathcal{A}$ represents an open neighborhood of $\theta$, and
thus there exists $k(r)\in\mathbb{N}$ such that $\theta_{k}\notin\mathcal{A}$
and $I(\theta^{\prime},\theta_{k})\geq r$ for all $k\geq k(r)$. As $r$ was
chosen freely, this means that
$\lim_{k\rightarrow\infty}I(\theta^{\prime},\theta_{k})=\infty$. ∎

###### Proof of Lemma 3.7.

By the definition of the rate function we have

$\displaystyle I(\theta^{\prime},\theta)$
$\displaystyle=\frac{1}{2}\mathrm{tr}\left(S_{w}^{-1}(\theta^{\prime}-\theta)S_{\theta}(\theta^{\prime}-\theta)^{\mathsf{T}}\right)$
$\displaystyle\geq\frac{1}{2}\sigma_{\mathrm{min}}(S_{w}^{-1})\sigma_{\mathrm{min}}(S_{\theta})\|\theta^{\prime}-\theta\|_{F}^{2}$
$\displaystyle\geq\frac{1}{2\kappa(S_{w})}\|\theta^{\prime}-\theta\|_{2}^{2},$

where the third inequality holds because $S_{\theta}\succeq S_{w}$ and
$\sigma_{\mathrm{min}}(S_{w}^{-1})=1/\sigma_{\mathrm{max}}(S_{w})$. The claim
then follows by multiplying the above inequality with $2\kappa(S_{w})$ and
taking square roots on both sides. ∎

###  5.2 Proofs of Section 3.2

###### Proof of Proposition 3.8.

First, given some $\varepsilon_{I}>0$, it follows from Lemma 3.7 that
$\|\theta^{\prime}-\theta\|^{2}_{2}>2\varepsilon_{I}\kappa(S_{w})$, implies
$I(\theta^{\prime},\theta)>\varepsilon_{I}$. So, if we define the set
$\mathcal{D}=\\{\theta^{\prime}\in\mathbb{R}^{n\times
n}:\|\theta^{\prime}-\theta\|_{2}>\sqrt{2\varepsilon_{I}\kappa(S_{w})}\\}$,
then,
$\inf_{\theta^{\prime}\in\mathcal{D}}I(\theta^{\prime},\theta)\geq\varepsilon_{I}$.

Now recall that the transformed estimator from Proposition 3.4 equals
$\widehat{\vartheta}_{T}(\widehat{\theta}_{T},\theta)=\sqrt{T/a_{T}}(\widehat{\theta}_{T}-\theta)+\theta$.
Indeed,
$\|\widehat{\theta}_{T}-\theta\|_{2}\leq(T/a_{T})^{-1/2}\sqrt{2\varepsilon_{I}\kappa(S_{w})}$
is equivalent to
$\widehat{\vartheta}_{T}(\widehat{\theta}_{T},\theta)\notin\mathcal{D}$. We
can plug these ideas into our MDP framework

$\limsup\limits_{T\to\infty}\frac{1}{a_{T}}\log\mathbb{P}_{\theta}(\widehat{\vartheta}_{T}(\widehat{\theta}_{T},\theta)\in\mathcal{D})\leq-\inf_{\theta^{\prime}\in\mathsf{cl}\,\mathcal{D}}I(\theta^{\prime},\theta)$

and obtain

$\mathbb{P}_{\theta}(\widehat{\vartheta}_{T}(\widehat{\theta}_{T},\theta)\notin\mathcal{D})\geq
1-e^{-\varepsilon_{I}a_{T}+o(a_{T})}.$

Therefore, if the product of $\varepsilon_{I}$ and $a_{T}$ is sufficiently
large such that $o(a_{T})$ becomes negligible, that is, if
$a_{T}\geq(1/{\varepsilon_{I}})\big{(}\log(1/\beta)+o(a_{T})\big{)}$ we have
with probability at least $1-\beta$ that
$\|\widehat{\theta}_{T}-\theta\|_{2}\leq(T/a_{T})^{-1/2}\sqrt{2\varepsilon_{I}\kappa(S_{w})}$.

Then, the last part of the proof follows from the triangle inequality
$\|\theta-\mathcal{P}(\widehat{\theta}_{T})\|_{2}\leq\|\theta-\widehat{\theta}_{T}\|_{2}+\|\widehat{\theta}_{T}-\mathcal{P}(\widehat{\theta}_{T})\|_{2}$.
Specifically, set
$\varepsilon_{I}=I(\widehat{\theta}_{T},\mathcal{P}(\widehat{\theta}_{T}))$,
then, the first term can be bounded using the result from the prior paragraph
whereas the second term can be bounded using Lemma 3.7. As $T$, and thereby
$\varepsilon$, are fixed, this yields both
$\|\widehat{\theta}_{T}-\mathcal{P}(\widehat{\theta}_{T})\|_{2}\leq\varepsilon$
$\mathbb{P}_{\theta}$-almost surely and (3.3b).

Regarding the finite-sample bound, we appeal to [52, Section 6]. First, to
deal with the non-isotropic case, recall that if (1.1), for some $\theta$,
would be driven by noise with an identity covariance matrix, then a change of
coordinates $z=S_{w}^{\frac{1}{2}}x$ would bring us back to the general case
with the new system matrix being $S_{w}^{\frac{1}{2}}\theta
S_{w}^{-\frac{1}{2}}$. Now observe that

$\begin{array}[]{lll}\|\theta-\widehat{\theta}_{T}\|_{2}&\leq\frac{\varepsilon}{\kappa(S_{w})^{1/2}}&\iff\\\
\|\theta-\widehat{\theta}_{T}\|_{2}\kappa(S_{w})^{\frac{1}{2}}&\leq\varepsilon&\implies\\\
\|S_{w}^{\frac{1}{2}}(\theta-\widehat{\theta}_{T})S_{w}^{-\frac{1}{2}}\|_{2}&\leq\varepsilon&\end{array}$

since
$\|S_{w}^{\frac{1}{2}}\|_{2}\|S_{w}^{-\frac{1}{2}}\|_{2}=\kappa(S_{w}^{\frac{1}{2}})=\kappa(S_{w})^{\frac{1}{2}}$.
Using this observation we recover bound (3.3a) by setting $\varepsilon$ in
[52, § 6] to $\varepsilon$ defined in Proposition 3.8. ∎

###  5.3 Proofs of Section 3.3

###### Proof of Proposition 3.9.

The follows immediately from the fact that for any $Q\in\mathcal{S}^{n}_{\succ
0}$, (3.4) is only finite when $\bar{\theta}\in\Theta$ indeed. Let
$\theta_{r}^{\star}$ be a minimizing argument of (3.4), then, if
$\theta^{\star}_{\underline{r}}$ was not the reverse $I$-projection it would
contradict the definition of $\underline{r}$. ∎

The numerical reverse $I$-projection relies on a result from infinite-horizon
dynamic programming, e.g., see [3, Chapter 3] or [9], and the following exact
constraint relaxation result borrowed from [30, Lemma A-0.1]; see also [31].
We repeat this to keep the paper self-contained.

######  Lemma 5.2 (Exact constraint relaxation).

Let $f$ and $g$ be arbitrary functions from an arbitrary set $\mathcal{X}$ to
$(-\infty,\infty]$, and consider the two closely related maximization problems

$\displaystyle\mathcal{P}_{1}(r)\leavevmode\nobreak\ :\leavevmode\nobreak\
\sup_{x\in\mathcal{X}}\left\\{f(x):g(x)\leq
r\right\\},\quad\mathcal{P}_{2}(\delta)\leavevmode\nobreak\
:\leavevmode\nobreak\ \sup_{x\in\mathcal{X}}f(x)-\frac{g(x)}{\delta}$

parametrized by $r\in\mathbb{R}$ and $\delta\in\mathbb{R}_{>0}$, respectively.
If the penalty-based maximization problem $\mathcal{P}_{2}(\delta)$ admits an
optimal solution $x^{\star}_{2}(\delta)$ for the parameter values $\delta$
within some set $\Delta\subset\mathbb{R}_{>0}$, then the following hold.

  1. (i)

The function $h(\delta)=g(x^{\star}_{2}(\delta))$ is non-decreasing in the
parameter $\delta\in\Delta$.

  2. (ii)

If there exits $\delta\in\Delta$ with $h(\delta)=r$, then the constrained
maximization problem $\mathcal{P}_{1}(r)$ is solved by
$x^{\star}_{2}(\delta)$.

###### Proof.

In order to prove Assertion (i), choose any parameters
$\delta_{1},\delta_{2}\in\Delta$ with $\delta_{1}>\delta_{2}$. As
$x^{\star}_{2}(\delta_{1})$ is optimal in $\mathcal{P}_{2}(\delta_{1})$ and
$x^{\star}_{2}(\delta_{2})$ is optimal in $\mathcal{P}_{2}(\delta_{2})$, one
can readily verify that

$\displaystyle
f\big{(}x^{\star}_{2}(\delta_{1})\big{)}-\frac{g(x^{\star}_{2}(\delta_{1}))}{\delta_{1}}$
$\displaystyle\geq
f(x^{\star}_{2}(\delta_{2}))-\frac{g(x^{\star}_{2}(\delta_{2}))}{\delta_{1}}$

and

$\displaystyle
f(x^{\star}_{2}(\delta_{2}))-\frac{g(x^{\star}_{2}(\delta_{2}))}{\delta_{2}}$
$\displaystyle\geq
f\big{(}x^{\star}_{2}(\delta_{1})\big{)}-\frac{g(x^{\star}_{2}(\delta_{1}))}{\delta_{2}}.$

Summing up these two inequalities yields

$\displaystyle\left(\frac{1}{\delta_{2}}-\frac{1}{\delta_{1}}\right)g\big{(}x^{\star}_{2}(\delta_{2})\big{)}\leq\left(\frac{1}{\delta_{2}}-\frac{1}{\delta_{1}}\right)g\big{(}x^{\star}_{2}(\delta_{1})\big{)}\iff
h(\delta_{2})=g\big{(}x^{\star}_{2}(\delta_{2})\big{)}\leq
g\big{(}x^{\star}_{2}(\delta_{1})\big{)}=h(\delta_{1}),$

where the equivalence holds because $\delta_{1}>\delta_{2}$. This completes
the proof of Assertion (i).

As for Assertion (ii), fix any $r\in\mathbb{R}$ and assume that there exists
$\delta\in\Delta$ with $r=h(\delta)$. We need to show that the optimizer
$x^{\star}_{2}(\delta)$ of $\mathcal{P}_{2}(\delta)$ is also optimal in
$\mathcal{P}_{1}(r)$. To this end, observe that $x^{\star}_{2}(\delta)$ is
feasible in $\mathcal{P}_{1}(\gamma)$ because $r=h(\delta)$. It then suffices
to prove optimality. Assume for the sake of contradiction that there exists
$x^{\prime}_{1}\in\mathcal{X}$ with
$f(x^{\prime}_{1})>f(x^{\star}_{2}(\delta))$ and $g(x^{\prime}_{1})\leq
g(x^{\star}_{2}(\delta))=r$. In this case, we have

$\displaystyle
f(x^{\prime}_{1})-\frac{g(x^{\prime}_{1})}{\delta}>f\big{(}x^{\star}_{2}(\delta)\big{)}-\frac{g(x^{\star}_{2}(\delta))}{\delta}\,,$

which contradicts the optimality of $x^{\star}_{2}(\delta)$ in
$\mathcal{P}_{2}(\delta)$. We thus conclude that $x^{\star}_{2}(\delta)$ must
indeed solve $\mathcal{P}_{1}(r)$. ∎

Another useful result is the following.

######  Lemma 5.3 (Matrix Inversion [66, page 19]).

Let $A$ and $C$ be invertible matrices, then,

$\displaystyle(A$ $\displaystyle+BCD)^{-1}=A^{-1}$
$\displaystyle-A^{-1}B(C^{-1}+DA^{-1}B)^{-1}DA^{-1}.$ (5.5)

###### Proof of Proposition 3.10.

Let $\theta^{\prime}\notin\Theta$. The proof relies on transforming (3.4) into
a unconstrained problem. To that end, we appeal to Lemma 5.2, yet, we change
the programs to

$\displaystyle\mathcal{P}_{1}(r)=\leavevmode\nobreak\
\min_{x\in\mathcal{X}}\left\\{f(x):g(x)\leq
r\right\\},\quad\mathcal{P}_{2}(\delta)=\leavevmode\nobreak\
\min_{x\in\mathcal{X}}f(x)+\frac{g(x)}{\delta}.$

and under the assumption that $\mathcal{P}_{2}(\delta)$ can be solved, we
obtain that a solution to $\mathcal{P}_{1}(r)$ is given by
$x_{2}^{\star}(\delta)$ if $g\big{(}x_{2}^{\star}(\delta)\big{)}=r$. Hence,
focusing on $\mathcal{P}_{2}(\delta)$, we appeal to a classical dynamic
programming formulation and recover the classical average cost LQR problem
since if $\mathcal{P}_{1}(r)$ is of the form (3.4), then, problem
$\mathcal{P}_{2}(\delta)$ becomes

$\begin{array}[]{cl}\min\limits_{L\in\mathbb{R}^{n\times
n}}&\displaystyle\lim_{T\to\infty}\frac{1}{T}\mathbb{E}_{\theta^{\prime}+L}\left[\textstyle\sum^{T-1}_{k=0}c(x_{k},Lx_{k})\right]\\\
\mathrm{s.t.}&x_{k+1}=(\theta^{\prime}+L)x_{k}+w_{k},\quad
x_{0}\sim\nu\end{array}$ (5.6)

for $c(x,Lx)=x^{\mathsf{T}}\left(Q+L^{\mathsf{T}}(2\delta
S_{w})^{-1}L\right)x$. Since we trivially satisfy standard stabilizability and
detectability assumptions [8, Chapter 4], the solution to (5.6) is given by
${L}_{\delta}=-({P}_{\delta}+(2\delta
S_{w})^{-1})^{-1}{P}_{\delta}\theta^{\prime}$ for ${P}_{\delta}$ the unique
symmetric positive definite solution to

$P=Q+\theta^{\prime\mathsf{T}}P\theta^{\prime}-\theta^{\prime\mathsf{T}}P(P+(2\delta
S_{w})^{-1})^{-1}P\theta^{\prime}.$ (5.7)

Moreover, using [49, Lemma 3.2] one can show that $P_{\delta}$ is real-
analytic in $\delta$ over the open interval $(0,\infty)$. Combining this with
the monotonicty from Lemma 5.2 yields that the real-analytic map
$I(\theta^{\prime},\theta^{\prime}+L_{\delta})\equiv\varphi^{-1}(\delta)$ is
strictly monotone in $\delta$ and hence bijective. Then, regarding the domain
of $\varphi$, the key observation here is the following. Consider e.g., (5.6),
then clearly
$\lim_{\delta\uparrow\infty}\varphi^{-1}(\delta)=\overline{r}=I(\theta^{\prime},0)$
and in a similar vain $\lim_{\delta\downarrow
0}\varphi^{-1}(\delta)=\underline{r}$, a different limit would contradict
optimality of $L_{\delta}$ in (5.6).

Hence, for each $r\in(\underline{r},\overline{r})$ there is some
$\delta=\varphi(r)$, such that
$I(\theta^{\prime},\theta^{\prime}+L_{\delta})=r$. Given such a radius $r$,
the main claim now follows simply by setting
$\theta_{\delta}^{\star}=\theta^{\prime}+L_{\varphi(r)}$, parametrized by
$P_{\delta}=P_{\varphi(r)}$.

Note that the matrix inversion formula (5.5) then enables us to reformulate
the above Riccati equation as (3.5). ∎

###### Proof of Proposition 3.11.

Given some $\theta^{\prime}\notin\Theta$, it follows from Proposition 3.10
that $\lim_{\delta\downarrow
0}\theta^{\star}_{\delta}=\mathcal{P}(\theta^{\prime})$. However, one cannot
appeal to $\theta^{\star}_{\delta}$ for $\delta\equiv 0$ since (3.5) would
have no solution in $\mathcal{S}^{n}_{\succ 0}$. Nevertheless, combining Lemma
3.7 with the fact that $I(\theta^{\prime},\theta^{\star}_{\delta})$ is real-
analytic in $\delta$ over $(0,\infty)$ implies the error complexity bound
(3.6).

Then, the fact that $\theta^{\star}_{\delta}$ can be computed using a standard
(discrete-time) LQR routine follows directly the proof of Proposition 3.10 as
well, e.g., see (5.6).

At last, the computational complexity mainly follows from the fact that the
numerically optimal way to solve (3.5) utilizes a QZ algorithm and requires
$O(n^{3})$ flops in general, see [48] and [22, Algorithm 7.7.3]. Although we
generically need ${O}(n^{3})$ flops and ${O}(n^{2})$ memory (see the
references below), for large scale problems one should resort to different
schemes like proposed in [20, 4]. ∎

###### Proof of Corollary 3.12.

Let $\theta^{\prime}\in\Theta^{\prime}$, then, given Proposition 3.10 it
follows directly from [29] that for any $\delta\in(0,\infty)$ there is a
$\Lambda_{\delta}$ with $\mathrm{det}(\Lambda_{\delta})>0$ such that
$\theta^{\star}_{\delta}=\Lambda_{\delta}^{-1}\theta^{\prime}$. Now, as
$\Lambda_{\delta}=(I_{n}+2\delta S_{w}P_{\delta})$, $S_{w}\succ 0$,
$P_{\delta}\succeq Q\succ 0$, $P_{\delta}$ is real-analytic in $\delta$ over
$(0,\infty)$ and $\lambda_{i}(S_{w}P_{\delta})>0$ for $i\in[n]$ we have that
$\lim_{\delta\downarrow 0}\mathrm{det}(\Lambda_{\delta})>0$, which concludes
the proof. ∎

###  5.4 Additional numerical results to Section 4

Figure 5.1 provides additional numerical results for the examples in Section
4.

(a) Convergence in operator norm for $n=100$.

(b) Convergence of spectral radii for $n=100$.

(c) Runtime comparison on $\theta^{\prime}=(Y\otimes 2I_{3})$.

(d) Historgrams of the spectral radii for $m=9$.

Figure 5.1: Additional material related to Section 4.

#### Acknowledgements

This research was supported by the Swiss National Science Foundation as part
of the National Centre of Competence in Research on Dependable, Ubiquitous
Automation (NCCR Automation).

### Bibliography

## References

  * [1] Panos J. Antsaklis and Anthony N. Michel  “Linear Systems”  Birkhäuser, 2006
  * [2] Karl Johan Åström and Björn Wittenmark  “On self tuning regulators”  In _Automatica_ 9.2 Elsevier, 1973, pp. 185–199
  * [3] Tamer Başar and Pierre Bernhard  “${H_{\infty}}$-Optimal Control and Related Minimax Design Problems A Dynamic Game Approach”  Birkhäuser, 1995
  * [4] Peter Benner and Heike Fassbender  “On the numerical solution of large-scale sparse discrete-time Riccati equations”  In _Advances in Computational Mathematics_ 35, 2011, pp. 119–147
  * [5] Felix Berkenkamp, Matteo Turchetta, Angela P. Schoellig and Andreas Krause  “Safe Model-based Reinforcement Learning with Stability Guarantees”  In _Neural Information Processing Systems_ , 2017
  * [6] D. S. Bernstein  “Matrix Mathematics”  Princeton University Press, 2009
  * [7] D. P. Bertsekas and J. Tsitsiklis  “Neuro-Dynamic Programming”  Athena Scientific, 1996
  * [8] Dimitri P. Bertsekas  “Dynamic Programming and Optimal Control”, Volume 1  Athena Scientific, 2005
  * [9] Dimitri P. Bertsekas  “Dynamic Programming and Optimal Control”, Volume 2  Athena Scientific, 2007
  * [10] Dimitri P Bertsekas  “Reinforcement Learning and Optimal Control”  Athena Scientific, 2019
  * [11] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B Shah  “Julia: A fresh approach to numerical computing”  In _SIAM Review_ 59.1 SIAM, 2017, pp. 65–98
  * [12] Nicholas M. Boffi et al.  “Learning Stability Certificates from Data”  In _arXiv preprint arXiv:2008.05952_ , 2020
  * [13] Byron Boots, Geoffrey J Gordon and Sajid M. Siddiqi  “A Constraint Generation Approach to Learning Stable Linear Dynamical Systems”  In _Advances in Neural Information Processing Systems 20_ , 2008, pp. 1329–1336
  * [14] Marco C Campi and P. R. Kumar  “Adaptive linear quadratic Gaussian control: The cost-biased approach revisited”  In _SIAM Journal on Control and Optimization_ 36.6 SIAM, 1998, pp. 1890–1907
  * [15] Neelam Choudhary, Nicolas Gillis and Punit Sharma  “On approximating the nearest ${\Omega}$-stable matrix”  In _Numerical Linear Algebra with Applications_ 27.3, 2020
  * [16] I. Csiszar and F. Matus  “Information projections revisited”  In _IEEE Transactions on Information Theory_ 49.6, 2003, pp. 1474–1490
  * [17] I. Csiszar and P.C. Shields  “Information Theory and Statistics: A Tutorial”  In _Foundations and Trends in Communications and Information Theory_ 1.4, 2004, pp. 417–528
  * [18] Imre Csiszar  “Sanov Property, Generalized $I$-Projection and a Conditional Limit Theorem”  In _Annals of Probability_ 12.3 The Institute of Mathematical Statistics, 1984, pp. 768–793
  * [19] A. Dembo and O. Zeitouni  “Large Deviations Techniques and Applications”  Springer, 2009
  * [20] Judith D. Gardiner and Alan J. Laub  “Parallel algorithms for algebraic Riccati equations”  In _International Journal of Control_ 54.6, 1991, pp. 1317–1333
  * [21] Nicolas Gillis, Michael Karow and Punit Sharma  “Approximating the nearest stable discrete-time system”  In _Linear Algebra and its Applications_ 573, 2019, pp. 37–53
  * [22] Gene H. Golub and Charles F. Loan  “Matrix Computations”  John Hopkins University Press, 2013
  * [23] Nicola Guglielmi and Christian Lubich  “Matrix stabilization using differential equations”  In _SIAM Journal on Numerical Analysis_ 55.6, 2017, pp. 3097–3119
  * [24] Nicola Guglielmi and Vladimir Yu Protasov  “On the closest stable/unstable nonnegative matrix and related stability radii”  In _SIAM Journal on Matrix Analysis and Applications_ 39.4 SIAM, 2018, pp. 1642–1669
  * [25] Nicholas J Higham  “Matrix nearness problems and applications”  In _Applications of Matrix Theory_ Oxford University Press, 1989, pp. 1–27
  * [26] F. Hollander  “Large Deviations”  American Mathematical Society, 2008
  * [27] Yassir Jedra and Alexandre Proutiere  “Sample Complexity Lower Bounds for Linear System Identification”  In _arXiv preprint arXiv:1903.10343_ , 2019
  * [28] Yassir Jedra and Alexandre Proutiere  “Finite-time Identification of Stable Linear Systems: Optimality of the Least-Squares Estimator”  In _arXiv preprint arXiv:2003.07937_ , 2020
  * [29] W. Jongeneel and D. Kuhn  “On Topological Equivalence in Linear Quadratic Optimal Control”  In _submitted_ , 2020
  * [30] Wouter Jongeneel  “Controlling the Unknown: A Game Theoretic Perspective”, 2019
  * [31] Wouter Jongeneel, Tyler Summers and Peyman Mohajerin Esfahani  “Robust Linear Quadratic Regulator: Exact Tractable Reformulation”  In _IEEE Conference on Decision and Control_ , 2019, pp. 6742–6747
  * [32] J. Zico Kolter and Gaurav Manek  “Learning Stable Deep Dynamics Models”  In _Advances in Neural Information Processing Systems 32_ , 2019, pp. 11128–11136
  * [33] Steven G. Krantz and Harold R. Parks  “A Primer of Real Analytic Functions”  Birkhäuser, 2002
  * [34] P. R. Kumar and Pravin Varaiya  “Stochastic Systems: Estimation, Identification and Adaptive Control”  USA: Prentice-Hall, Inc., 1986
  * [35] S. L. Lacy and D. S. Bernstein  “Subspace identification with guaranteed stability using constrained optimization”  In _Proceedings of the 2002 American Control Conference_ 4, 2002, pp. 3307–3312 vol.4
  * [36] S. L. Lacy and D. S. Bernstein  “Subspace identification with guaranteed stability using constrained optimization”  In _IEEE Transactions on Automatic Control_ 48.7, 2003, pp. 1259–1263
  * [37] Peter Lancaster and Leiba Rodman  “Algebraic Riccati Equations”  Oxford University Press, 1995
  * [38] Tor Lattimore and Csaba Szepesvári  “Bandit Algorithms”  Cambridge University Press, 2020
  * [39] Elliott Lieb  “Convex trace functions and the Wigner-Yanase-Dyson conjecture”  In _Advances in Mathematics_ 11.3, 1973, pp. 267–288
  * [40] J.M. Maciejowski  “Guaranteed stability with subspace methods”  In _Systems & Control Letters_ 26.2, 1995, pp. 153 –156
  * [41] Giorgos Mamakoukas, Ian Abraham and Todd D Murphey  “Learning Data-Driven Stable Koopman Operators”  In _arXiv preprint arXiv:2005.04291_ , 2020
  * [42] Giorgos Mamakoukas, Orest Xherija and T. D. Murphey  “Memory-Efficient Learning of Stable Linear Dynamical Systems for Prediction and Control”  In _arXiv e-prints arXiv:2006.03937_ , 2020
  * [43] Nikolai Matni, Alexandre Proutiere, Anders Rantzer and Stephen Tu  “From self-tuning regulators to reinforcement learning and back again”  In _arXiv preprint arXiv:1906.11392_ , 2019
  * [44] Sean Meyn and Richard L Tweedie  “Markov Chains and Stochastic Stability”  Cambridge University Press, 2009
  * [45] S. Mohammad Khansari-Zadeh and Aude Billard  “Learning control Lyapunov function to ensure stability of dynamical system-based robot reaching motions”  In _Robotics and Autonomous Systems_ 62.6, 2014, pp. 752 –765
  * [46] Yu. Nesterov and V. Yu. Protasov  “Computing Closest Stable Nonnegative Matrix”  In _SIAM Journal on Matrix Analysis and Applications_ 41.1, 2020, pp. 1–28
  * [47] Francois-Xavier Orbandexivry, Yurii Nesterov and Paul Van Dooren  “Nearest stable system using successive convex approximations”  In _Automatica_ 49.5, 2013, pp. 1195 –1203
  * [48] T. Pappas, A. Laub and N. Sandell  “On the numerical solution of the discrete-time algebraic Riccati equation”  In _IEEE Transactions on Automatic Control_ 25.4, 1980, pp. 631–641
  * [49] Jan Willem Polderman  “A note on the structure of two subsets of the parameter space in adaptive control problems”  In _Systems & Control Letters_ 7.1, 1986, pp. 25–34
  * [50] Warren B Powell  “Approximate Dynamic Programming: Solving the Curses of Dimensionality”  John Wiley & Sons, 2007
  * [51] Benjamin Recht  “A tour of reinforcement learning: The view from continuous control”  In _Annual Review of Control, Robotics, and Autonomous Systems_ 2 Annual Reviews, 2019, pp. 253–279
  * [52] Tuhin Sarkar and Alexander Rakhlin  “Near optimal finite time identification of arbitrary linear dynamical systems”  In _International Conference on Machine Learning_ , 2019, pp. 5610–5618
  * [53] Tuhin Sarkar, Alexander Rakhlin and Munther A. Dahleh  “Nonparametric Finite Time LTI System Identification”  In _arXiv e-prints arXiv:1902.01848_ , 2020
  * [54] Max Simchowitz et al.  “Learning Without Mixing: Towards A Sharp Analysis of Linear System Identification”  In _Conference On Learning Theory_ , 2018, pp. 439–473
  * [55] Rangarajan K. Sundaram  “A First Course in Optimization Theory”  Cambridge University Press, 1996
  * [56] Tobias Sutter, Bart P. G. Van Parys and Daniel Kuhn  “A General Framework for Optimal Data-Driven Optimization”  In _arXiv preprint, arXiv:2010.06606_ , 2020
  * [57] Richard S. Sutton and Andrew G. Barto  “Reinforcement Learning: An Introduction”  MIT Press, 2018
  * [58] Hideyuki Tanaka and Tohru Katayama  “Stochastic subspace identification guaranteeing stability and minimum phase”  In _IFAC Proceedings Volumes_ 38.1, 2005, pp. 910 –915
  * [59] K. Turksoy et al.  “Guaranteed stability of recursive multi-input-single-output time series models”  In _2013 American Control Conference_ , 2013, pp. 77–82
  * [60] Jack Umenberger, Johan Wågberg, Ian R. Manchester and Thomas B. Schön  “Maximum likelihood identification of stable linear dynamical systems”  In _Automatica_ 96, 2018, pp. 280 –292
  * [61] Jonas Umlauft and Sandra Hirche  “Learning stochastically stable Gaussian process state–space models”  In _IFAC Journal of Systems and Control_ 12, 2020
  * [62] T. Van Gestel, J. A. K. Suykens, P. Van Dooren and B. De Moor  “Imposing stability in subspace identification by regularization”  In _Proceedings of the 39th IEEE Conference on Decision and Control_ , 2000, pp. 1555–1560 vol.2
  * [63] T. Van Gestel, J. A. K. Suykens, P. Van Dooren and B. De Moor  “Identification of stable models in subspace identification by using regularization”  In _IEEE Transactions on Automatic Control_ 46.9, 2001, pp. 1416–1420
  * [64] Peter Van Overschee and Bart De Moor  “Subspace identification for linear systems: Theory-Implementation-Applications”  Kluwer Academic Publishers, 1996
  * [65] Bart Van Parys, Peyman Mohajerin Esfahani and Daniel Kuhn  “From Data to Decisions: Distributionally Robust Optimization Is Optimal” Articles in Advance  In _Management Science_ , 2021
  * [66] Michel Verhaegen and Vincent Verdult  “Filtering and System Identification”  Cambridge University Press, 2007
  * [67] Julien Worms  “Moderate deviations for stable Markov chains and regression models”  In _Electronic Journal of Probability_ 4 The Institute of Mathematical Statisticsthe Bernoulli Society, 1999, pp. 28 pp.
  * [68] Miao Yu and Shen Si  “Moderate deviation principle for autoregressive processes”  In _Journal of Multivariate Analysis_ 100.9 Elsevier, 2009, pp. 1952–1961