# A Tutorial on the Non-Asymptotic Theory of System Identification

Ingvar Ziemann  University of Pennsylvania  Anastasios Tsiamis  ETH Zürich
Bruce Lee  University of Pennsylvania  Yassir Jedra  KTH Royal Institute of
Technology  Nikolai Matni  University of Pennsylvania  George J. Pappas
University of Pennsylvania

###### Abstract

This tutorial serves as an introduction to recently developed non-asymptotic
methods in the theory of—mainly linear—system identification. We emphasize
tools we deem particularly useful for a range of problems in this domain, such
as the covering technique, the Hanson-Wright Inequality and the method of
self-normalized martingales. We then employ these tools to give streamlined
proofs of the performance of various least-squares based estimators for
identifying the parameters in autoregressive models. We conclude by sketching
out how the ideas presented herein can be extended to certain nonlinear
identification problems.

###### Contents

  1. 1 Introduction
    1. 1.1 Problem Formulation
    2. 1.2 Least Squares Regression and the Path Ahead
    3. 1.3 Overview
  2. 2 Preliminaries: Concentration Inequalities, Packing and Covering
    1. 2.1 Sub-Gaussian Concentration and the Hanson-Wright Inequality
    2. 2.2 Covering and Discretization Arguments
    3. 2.3 Concentration of the Covariance Matrix of Linear Systems
    4. 2.4 Notes
  3. 3 The Lower Spectrum of the Empirical Covariance
    1. 3.1 A Decoupling Inequality for sub-Gaussian Quadratic Forms
    2. 3.2 The Lower Tail of the Empirical Covariance of Causal sub-Gaussian Processes
    3. 3.3 Notes
  4. 4 Self-Normalized Martingale Bounds
    1. 4.1 Exponential Inequalities via Pseudo-maximization
    2. 4.2 Self-Normalized Martingales Satisfy the Canonical Assumption
    3. 4.3 Notes
  5. 5 System Identification
    1. 5.1 ARX Systems
      1. 5.1.1 Persistency of Excitation in ARX Models
      2. 5.1.2 Dealing with the Noise Term
    2. 5.2 State-Space Systems
      1. 5.2.1 Reduction to ARX learning with Bias
      2. 5.2.2 Non-Asymptotic Guarantees
    3. 5.3 Notes
  6. 6 An Alternative Viewpoint: the Basic Inequality
    1. 6.1 Sparse Autoregressions
    2. 6.2 Notes
  7. 7 Beyond Linear Models
    1. 7.1 Many Trajectories and Finite Hypothesis Classes
    2. 7.2 Notes
  8. A Proof of The Hanson-Wright Inequality
    1. A.1 Gaussian Comparison Inequalities for sub-Gaussian Quadratic Forms
    2. A.2 Finishing the proof of Theorem 2.1
  9. B Proof of Theorem 2.2
    1. B.1 A Variant of the Hanson-Wright Inequality
    2. B.2 The Main Concentration Inequality
    3. B.3 Approximate Isometries
    4. B.4 Finishing the proof of Theorem 2.2
  10. C Proofs Relating to the Lower Tail of the Empirical Covariance
  11. D Proof of the Self-Normalized Martingale Theorem
    1. D.1 Proof of Theorem 4.1
  12. E Proofs for System Identification
    1. E.1 Proof of Theorem 5.1
    2. E.2 Proof of Theorem 5.2
    3. E.3 Proof of Theorem 5.4
    4. E.4 Proof of Theorem 5.3
  13. F Proofs for Section 6 and Section 7
    1. F.1 Proofs for Sparse Identification
    2. F.2 Proofs for Nonlinear Identification

## Notation

Maxima (resp. minima) of two numbers $a,b\in\mathbb{R}$ are denoted by $a\vee
b=\max(a,b)$ ($a\wedge b=\min(a,b)$). For two sequences
$\\{a_{t}\\}_{t\in\mathbb{Z}}$ and $\\{b_{t}\\}_{t\in\mathbb{Z}}$ we introduce
the shorthand $a_{t}\lesssim b_{t}$ if there exists a universal constant $C>0$
and an integer $t_{0}$ such that $a_{t}\leq Cb_{t}$ for every $t\geq t_{0}$.
If $a_{t}\lesssim b_{t}$ and $b_{t}\lesssim a_{t}$ we write $a_{t}\asymp
b_{t}$. Let $\mathsf{X}\subset\mathbb{R}^{d}$ and let $f,g\in\mathsf{X}\to R$.
We write $f=O(g)$ if $\limsup_{x\to x_{0}}|f(x)/g(x)|<\infty$, where the limit
point $x_{0}$ is typically understood from the context. We use $\tilde{O}$ to
hide logarithmic factors and write $f=o(g)$ if $\limsup_{x\to
x_{0}}|f(x)/g(x)|=0$. We write $f=\Omega(g)$ if $\limsup_{x\to
x_{0}}|f(x)/g(x)|>0$. For an integer $N$, we also define the shorthand
$[N]\triangleq\\{1,\dots,N\\}$.

Expectation (resp. probability) with respect to all the randomness of the
underlying probability space is denoted by $\mathbf{E}$ (resp. $\mathbf{P}$).

The Euclidean norm on $\mathbb{R}^{d}$ is denoted $\|\cdot\|_{2}$, and the
unit sphere in $\mathbb{R}^{d}$ is denoted $\mathbb{S}^{d-1}$. The standard
inner product on $\mathbb{R}^{d}$ is denoted $\langle\cdot,\cdot\rangle$. We
embed matrices $M\in\mathbb{R}^{d_{1}\times d_{2}}$ in Euclidean space by
vectorization: $\operatorname{\mathsf{vec}}M\in\mathbb{R}^{d_{1}d_{2}}$, where
$\operatorname{\mathsf{vec}}$ is the operator that vertically stacks the
columns of $M$ (from left to right and from top to bottom). For a matrix $M$
the Euclidean norm is the Frobenius norm, i.e.,
$\|M\|_{F}\triangleq\|\operatorname{\mathsf{vec}}M\|_{2}$. We similarly define
the inner product of two matrices $M,N$ by $\langle
M,N\rangle\triangleq\langle\operatorname{\mathsf{vec}}M,\operatorname{\mathsf{vec}}N\rangle$.
The transpose of a matrix $M$ is denoted by $M^{\mathsf{T}}$ and
$\operatorname{tr}M$ denotes its trace. For a matrix
$M\in\mathbb{R}^{d_{1}\times d_{2}}$, we order its singular values
$\sigma_{1}(M),\dots,\sigma_{d_{1}\wedge d_{2}}(M)$ in descending order by
magnitude. We also write $\|M\|_{\mathsf{op}}$ for its largest singular value:
$\|M\|_{\mathsf{op}}\triangleq\sigma_{1}(M)$. To not carry dimensional
notation, we will also use $\sigma_{\min}(M)$ for the smallest nonzero
singular value. For square matrices $M\in\mathbb{R}^{d\times d}$ with real
eigenvalues, we similarly order the eigenvalues of $M$ in descending order as
$\lambda_{1}(M),\dots,\lambda_{d}(M)$. In this case, $\lambda_{\min}(M)$ will
also be used to denote the minimum (possibly zero) eigenvalue of $M$. For two
symmetric matrices $M,N$, we write $M\succ N$ ($M\succeq N)$ if $M-N$ is
positive (semi-)definite.

##  1 Introduction

Machine learning methods are at an ever increasing pace being integrated into
domains that have classically been within the purview of controls. There is a
wide range of examples, including perception-based control, agile robotics,
and autonomous driving and racing. As exciting as these developments may be,
they have been most pronounced on the experimental and empirical sides. To
deploy these systems safely, stably, and robustly into the real world, we
argue that a principled and integrated theoretical understanding of a)
fundamental limitations and b) statistical optimality is needed. Under the
past few years, a host of new techniques have been introduced to our field.
Unfortunately, existing results in this area are relatively inaccessible to a
typical first or second year graduate student in control theory, as they
require both sophisticated mathematical tools not typically included in a
control theorist’s training (e.g., high-dimensional statistics and learning
theory).

This tutorial seeks to provide a streamlined exposition of some of these
recent advances that are most relevant to the non-asymptotic theory of linear
system identification. Our aim is not to be encyclopedic but rather to give
simple proofs of the main developments and to highlight and collect the key
technical tools to arrive at these results. For a broader—and less
technical—overview of the literature we point the reader to our recent survey
(Tsiamis et al., 2023). It is also worth to point out that the classical
literature on system identification has done a formidable job at—often very
accurately—characterizing the asymptotic performance of identification
algorithms (Ljung, 1999). Our aim is not to supplant this literature but
rather to complement the asymptotic picture with finite sample guarantees by
relaying recently developed technical tools drawn from high-dimensional
probability, statistics and learning theory (Vershynin, 2018; Wainwright,
2019).

###  1.1 Problem Formulation

Let us now fix ideas. We are concerned with linear time-series models of the
form:

$\displaystyle Y_{t}$ $\displaystyle=\theta^{\star}X_{t}+V_{t}$ $\displaystyle
t=1,2,\dots,T$ (1.1)

where $Y_{1:T}$ is a sequence of outputs (or targets) assuming values in
$\mathbb{R}^{d_{\mathsf{Y}}}$ and $X_{1:T}$ is a sequence of inputs (or
covariates) assuming values in $\mathbb{R}^{d_{\mathsf{X}}}$. The goal of the
user (or learner) is to recover the a priori unknown linear map
$\theta^{\star}\in\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$ using only
the observations $X_{1:T}$ and $Y_{1:T}$. The linear relationship in the
regression model (1.1) is perturbed by a stochastic noise sequence $V_{1:T}$
assuming values in $\mathbb{R}^{d_{\mathsf{Y}}}$. We refer to the regression
model (1.1) as a time-series to emphasize the fact that the observations
$X_{1:T}$ and $Y_{1:T}$ may arrive sequentially and in particular that past
$X_{t}$ and $Y_{t}$ may influence future $X_{t^{\prime}}$ and $Y_{t^{\prime}}$
(i.e. with $t^{\prime}>t)$.

##### Example: Autoregressive Models.

For instance, a model class of particular interest to us which is subsumed by
(1.1) are the (vector) autoregressive exogenous models of order $p$ and $q$
(briefly ARX$(p,q)$):

$\displaystyle Y_{t}$
$\displaystyle=\sum_{i=1}^{p}A^{\star}_{i}Y_{t-i}+\sum_{j=1}^{q}B^{\star}_{i}U_{t-j}+W_{t}$
(1.2)

where typically $U_{1:T-1}$ is a sequence of user specified inputs taking
values in $\mathbb{R}^{d_{\mathsf{U}}}$ and $W_{1:T}$ is an iid sequence of
noise variables taking values in $\mathbb{R}^{d_{\mathsf{W}}}$. If we are only
interested in the parameters
$\begin{bmatrix}A_{1:p}^{\star}&B_{1:q}^{\star}\end{bmatrix}$, we obtain the
model (1.2) by setting

$\displaystyle X_{t}$
$\displaystyle=\begin{bmatrix}Y_{t-1:t-p}^{\mathsf{T}}&U_{t-1:t-q}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}};$
$\displaystyle\theta^{\star}=\begin{bmatrix}A_{1:p}^{\star}&B_{1:q}^{\star}\end{bmatrix};$
$\displaystyle V_{t}=W_{t}.$ (1.3)

We point out that that the above discussion presupposes that the order of the
model, $(p,q)$, is known (there are ways around this).

In this tutorial we will provide the necessary tools to tackle the following
problem.

######  Problem 1.1.

Fix $\varepsilon>0$, $\delta\in(0,1)$, and a norm $\|\cdot\|$. Fix also a
‘reasonable’ estimator $\widehat{\theta}$ of $\theta_{\star}$ using a sample
$(X,Y)_{1:T}$ from (1.1). We seek to establish finite sample guarantees of the
form

$\|\widehat{\theta}-\theta^{\star}\|\leq\varepsilon\qquad\textnormal{ with
probability at least }1-\delta$ (1.4)

where $\varepsilon$ controls the accuracy (or rate) and the failure parameter
$\delta$ controls the confidence.

In the sequel, ‘reasonable’ estimator will typically mean some form of least
squares estimator (1.7). These are introduced in Section 1.2 below. A bound of
the form (1.4) is typically thought of as follows. We fix a priori the failure
parameter $\delta$ and then provide guarantees of the form
$\|\widehat{\theta}-\theta^{\star}\|\leq\varepsilon(T,\delta,\mathsf{P}_{XY})$
where $\mathsf{P}_{XY}$ is the joint distribution of $(X,Y)_{1:T}$. Hence, the
sample size $T$, the failure probability $\delta$ and the distribution of the
samples all impact the performance guarantee $\varepsilon$ we are able to
establish. To be more specific, $\varepsilon$ will typically be of the form

$\varepsilon\propto(\textnormal{Noise
Scale})\times\sqrt{\frac{\textnormal{problem
dimension}+\log(1/\delta)}{\textnormal{sample size}}}.$ (1.5)

Thus in principle, the best possible choice of $\varepsilon^{2}$ can be
thought of as a high probability version of the (inverse) signal-to-noise
ratio of the problem at hand. The fact that the confidence parameter $\delta$
typically affects (1.5) additively in $\log(1/\delta)$ is consistent with
classical asymptotic normality theory of estimators. One often expects the
normalized difference $T^{-1/2}(\widehat{\theta}-\theta^{\star})$ to converge
in law to a normal distribution (van der Vaart, 2000). In this tutorial we
will provide tools that allow us to match such classical asymptotics but with
a finite sample twist. Let us also remark that there often is a minimal
requirement on the sample size necessary for a bound of the form (1.4)-(1.5)
to hold. Such requirements are typically of the form

$\textnormal{sample size}\>\gtrsim\>\textnormal{problem
dimension}+\log(1/\delta).$ (1.6)

Requirements such as (1.6) are called burn-in times and are related to the
notion of persistence of excitation. They correspond to the rather minimal
requirement that the parameter identification problem is feasible in the
complete absence of observation noise.

###  1.2 Least Squares Regression and the Path Ahead

Let us now return to the general setting of (1.1). Fix a subset $\mathsf{M}$
of $\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$, called the model class.
The estimator

$\widehat{\theta}\in\operatorname*{argmin}_{\theta\in\mathsf{M}}\frac{1}{T}\sum_{t=1}^{T}\|Y_{t}-\theta
X_{t}\|_{2}^{2}$ (1.7)

is the least squares estimator (LSE) of $\theta^{\star}$ (with respect to
$\mathsf{M}$). Often we simply set
$\mathsf{M}=\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$. In this case,
equivalently:

$\widehat{\theta}=\left(\sum_{t=1}^{T}Y_{t}X_{t}^{\mathsf{T}}\right)\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{\dagger}$
(1.8)

and the LSE reduces to the (minimum norm) ordinary least squares (OLS)
estimator (1.8).

For simplicity, let us further assume that the (normalized) empirical
covariance matrix:

$\widehat{\Sigma}\triangleq\frac{1}{T}\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}};$
(1.9)

is full rank almost surely.

##### The Path Ahead.

Let us now briefly sketch the path ahead to solve 1.1. If (1.9) is full
rank—as required above—the estimator (1.8) admits the convenient error
representation:

$\widehat{\theta}-\theta^{\star}=\left[\left(\sum_{t=1}^{T}V_{t}X_{t}^{\mathsf{T}}\right)\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{-1/2}\right]\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{-1/2}.$
(1.10)

The leftmost term of (1.10) (in square brackets) can be shown to be (almost)
time-scale invariant in many situations. For instance, if the noise $V_{1:T}$
is a sub-Gaussian martingale difference sequence with respect to the
filtration generated by the covariates $X_{1:T}$, one can invoke methods from
the theory of self-normalized processes to show this (Peña et al., 2009;
Abbasi-Yadkori, 2013). These methods are the topic of Section 4.

Whenever this is the case, the dominant term in the rate of convergence of the
least squares estimator is
$\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{-1/2}$. In other words,
providing control of the smallest eigenvalue of (1.9) effectively yields
control of the rate of convergence of the least squares estimator in many
situations. Thus, to analyze the rate of convergence of (1.7) when
$\mathsf{M}=\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$ it suffices to:

  * •

Analyze the smallest eigenvalue (or lower tail) of (1.9). We provide such
analyses in Section 3

  * •

Analyze the scale invariant term (in square brackets) of (1.10). This can in
many situations be handled for instance by the self-normalized martingale
method described in Section 4.

###  1.3 Overview

Before covering these more technical topics in Section 3 and Section 4, we
also briefly review some preliminaries from probability theory in Section 2.
We then demonstrate how to apply these ideas in the setting of identifying the
parameters of an ARX$(p,q)$ model of the form (1.2) in Section 5. An
alternative perspective not based on the decomposition (1.10) for more general
least squares algorithms is given in Section 6. We conclude with a brief
discussion on how the tools in Section 6 can be extended to study more general
nonlinear phenomena in Section 7.

##  2 Preliminaries: Concentration Inequalities, Packing and Covering

Before we proceed to tackle the more advanced question of analyzing the LSE
(1.7), let us discuss a few preliminary inequalities that control the tail of
a random variable. Our first inequality is Markov’s.

######  Lemma 2.1.

Let $X$ be a nonnegative random variable. For every $s>0$ we have that

$\mathbf{P}(X\geq s)\leq s^{-1}\mathbf{E}[X].$ (2.1)

###### Proof.

We have that $\mathbf{E}[X]\geq\mathbf{E}[\mathbf{1}_{X\geq s}X]\geq
s\mathbf{E}[\mathbf{1}_{X\geq s}]$. Since $\mathbf{E}[\mathbf{1}_{X\geq
s}]=\mathbf{P}(X\geq s)$ the result follows by rearranging. ∎

Typically, Markov’s inequality itself is insufficient for our goals: we seek
deviation inequalities that taper of exponentially fast in $s$ and not as
$s^{-1}$. Such scaling is for instance predicted asymptotically by the central
limit theorem by the asymptotic normality of renormalized sums of square
integrable iid random variables; that is, sums of the form
$S_{n}/\sqrt{n}=(X_{1}+X_{2}+\dots+X_{n})/\sqrt{n}$ where the $X_{i},i\in[n]$
are independent and square integrable. For random variables possessing a
moment generating function, Markov’s inequality can be "boosted" by the so-
called "Chernoff trick". Namely, we apply Markovs inequality to the moment
generating function of the random variable instead of applying it directly to
the random variable itself.

######  Corollary 2.1 (Chernoff).

Fix $s>0$ and suppose that $\mathbf{E}\exp\left(\lambda X\right)$ exists. Then

$\mathbf{P}\left(X\geq s\right)\leq\min_{\lambda\geq 0}e^{-\lambda
s}\mathbf{E}\exp\left(\lambda X\right).$ (2.2)

###### Proof.

Fix $\lambda\geq 0$. We have:

$\displaystyle\mathbf{P}\left(X\geq s\right)$
$\displaystyle=\mathbf{P}\left(\exp\left(\lambda X\right)\geq\exp\left(\lambda
s\right)\right)$ $\displaystyle(\textnormal{monotonicity of }x\mapsto
e^{\lambda x})$ $\displaystyle\leq e^{-\lambda s}\mathbf{E}\exp\left(\lambda
X\right)$ $\displaystyle(\textnormal{Markov's inequality}).$

The result follows by optimizing. ∎

Recall that the function
$\psi_{X}(\lambda)\triangleq\mathbf{E}\exp\left(\lambda X\right)$ is the
moment generating function of $X$. For instance, if $X$ has univariate
Gaussian distribution with mean zero and variance $\sigma^{2}$, the moment
generating function appearing in (2.2) is just $\mathbf{E}\exp\left(\lambda
X\right)=\exp\left(\lambda^{2}\sigma^{2}/2\right)$. Hence the probability that
said Gaussian exceeds $s$ is upper-bounded:

$\mathbf{P}\left(X>s\right)\leq\min_{\lambda\geq 0}e^{-\lambda
s}\exp\left(\lambda^{2}\sigma^{2}/2\right)=\exp\left(\frac{-s^{2}}{2\sigma^{2}}\right)$
(2.3)

which (almost) exhibits the correct Gaussian tails as compared to (2.1).111We
write almost because $\exp(-s^{2}/2\sigma^{2})\approx\mathbf{P}(V>s)$ where
$V\sim N(0,\sigma^{2})$ but the expression is not exact. It should be pointed
out that assumptions stronger than those of the Central Limit Theorem (finite
variance) are indeed needed for a non-asymptotic theory with sub-Gaussian
tails as in (2.3). An assumption of this kind which is relatively standard in
the literature is introduced next.

###  2.1 Sub-Gaussian Concentration and the Hanson-Wright Inequality

In the sequel, we will not want to impose the Gaussian assumption. Instead, we
define a class of random variables that admit reasoning analogous to (2.3).

######  Definition 2.1.

We say that a random vector $W$ taking values in $\mathbb{R}^{d}$ is
$\sigma^{2}$-sub-Gaussian ($\sigma^{2}$-subG) if for every
$v\in\mathbb{R}^{d}$ we have that:

$\mathbf{E}\exp\left(\langle
v,W\rangle\right)\leq\exp\left(\frac{\sigma^{2}\|v\|^{2}}{2}\right).$ (2.4)

Similarly, we say that $W$ is $\sigma^{2}$-conditionally sub-Gaussian with
respect to a $\sigma$-field $\mathcal{F}$ if (2.4) holds with
$\mathbf{E}[\cdot]$ replaced by $\mathbf{E}[\cdot|\mathcal{F}]$.

The term $\sigma^{2}$ appearing in (2.4) is called the variance proxy of a
sub-Gaussian random variable. The significance of this definition is that the
one-dimensional projections $X=\langle v,W\rangle$ (with $\|v\|=1$) satisfy
the tail inequality (2.3). While obviously Gaussian random variables are sub-
Gaussian with their variance as variance-proxy, there are many examples beyond
Gaussians that fit into this framework. It is for instance straightforward to
show that bounded random variables have variance proxy proportional to the
square of their width (see eg. Wainwright, 2019, Examples 2.3 and 2.4).
Moreover, it is readily verified that the normalized sum mentioned
above—$S_{n}/\sqrt{n}=(X_{1}+\dots+X_{n})/\sqrt{n}$—satisfies the same bound
(2.3) provided that the entries of $X_{1:n}$ are independent, mean zero and
$\sigma^{2}$-sub-Gaussian. To see this, notice that the moment generating
function "tensorizes" across products. Namely, for every
$\lambda\in\mathbb{R}$:

$\mathbf{E}\exp\left(\frac{\lambda}{\sqrt{n}}\sum_{i=1}^{n}X_{i}\right)=\prod_{i=1}^{n}\mathbf{E}\exp\left(\frac{\lambda}{\sqrt{n}}X_{i}\right)\leq\prod_{i=1}^{n}\exp\left(\frac{\lambda^{2}\sigma^{2}}{2n}\right)=\exp\left(\frac{\lambda^{2}\sigma^{2}}{2}\right).$
(2.5)

Hence, by the exact same reasoning leading up to (2.3) such normalized sub-
Gaussian sums satisfy the same tail bound (2.3).

When analyzing linear regression models, most variables of interest are
typically either linear or quadratic in the variables of interest (cf.
(1.10)). Hence, we also need to understand how squares of sub-Gaussian random
variables behave. The next result shows that sub-Gaussian quadratic forms
exhibit similar tail behavior to the Chi-squared distribution (often in the
literature referred to as sub-exponential tails). It is known as the Hanson-
Wright Inequality.

######  Theorem 2.1 (Hanson and Wright (1971); Rudelson and Vershynin (2013)).

Let $M\in\mathbb{R}^{d\times d}$. Fix a random variable $W=W_{1:d}$ where each
$W_{i},i\in[d]$ is a scalar, mean zero and independent $\sigma^{2}$-sub-
Gaussian random variable. Then for every $s\in[0,\infty)$:

$\mathbf{P}\left(|W^{\mathsf{T}}MW-\mathbf{E}W^{\mathsf{T}}MW|>s\right)\leq
2\exp\left(-\min\left(\frac{s^{2}}{144\sigma^{4}\|M\|_{F}^{2}},\frac{s}{16\sqrt{2}\sigma^{2}\|M\|_{\mathsf{op}}}\right)\right).$
(2.6)

The proof of Theorem 2.1 is rather long and technical and thus relegated to
Appendix A. There, the reader may also find further useful concentration
inequalities for quadratic forms in sub-Gaussian variables. In fact, there are
plethora of useful concentration inequalities not covered here and the
interested reader is urged to consult the first few chapters of Vershynin
(2018).

###  2.2 Covering and Discretization Arguments

We will often find ourselves in a situation where it is possible to obtain a
scalar concentration bound but need this to hold uniformly for many random
variables at once. The $\varepsilon$-net argument, which proceeds via the
notion of covering numbers, is a relatively straightforward way of converting
concentration inequalities for scalars into their counterparts for vectors,
matrices and functions more generally.

The reader will for instance notice that the quantity being controlled by
Theorem 2.1 is a scalar quadratic form in sub-Gaussian random variables. By
contrast, the empirical covariance matrix (1.9) is a matrix and so a
conversion step is needed. This idea will be used frequently and in various
forms throughout the manuscript, so we review it briefly here for the
particular case of controlling the operator norm of a random matrix. To this
end, we notice that for any matrix $M\in\mathbb{R}^{m\times d}$:

$\|M\|_{\mathsf{op}}^{2}=\max_{v\in\mathbb{S}^{d-1}}\langle Mv,Mv\rangle.$
(2.7)

Hence, the operator norm of a random matrix is a maximum of scalar random
variables indexed by the unit sphere $\mathbb{S}^{d-1}$.

Recall now that the union bound states that the probability that the maximum
of a _finite collection_ ($|S|<\infty$) $\\{X_{i}\\}_{i\in S}$ of random
variables exceeds a certain threshhold can be bounded by the sum of their
probabilities:

$\mathbf{P}\left(\max_{i\in S}X_{i}>t\right)\leq\sum_{i\in
S}\mathbf{P}\left(X_{i}>t\right).$ (2.8)

Unfortunately, the unit sphere appearing (2.7) is not a finite set and so the
union bound (2.8) cannot be directly applied. However, when the domain of
optimization has geometric structure, one can often exploit this to leverage
the union bound not directly but rather in combination with a discretization
argument. Returning to our example of the operator norm of a matrix, the set
$S$ appearing in (2.8) will be a discretized version of the unit sphere
$\mathbb{S}^{d-1}$.

The following notion is key.

######  Definition 2.2.

Let $(\mathsf{X},d)$ be a compact metric space and fix $\varepsilon>0$. A
subset $\mathcal{N}$ of $\mathsf{X}$ is called an $\varepsilon$-net of
$\mathsf{X}$ if every point of $\mathsf{X}$ is within radius $\varepsilon$ of
a point of $\mathcal{N}$:

$\sup_{x\in\mathsf{X}}\inf_{x^{\prime}\in\mathcal{N}}d(x,x^{\prime})\leq\varepsilon.$
(2.9)

Moreover, the minimal cardinality of $\mathcal{N}$ necessary such that (2.9)
holds is called the covering number at resolution $\varepsilon$ of
$(\mathsf{X},d)$ and is denoted $\mathcal{N}(\varepsilon,\mathsf{X},d)$.

We will not explore this notion in full, but simply content ourselves to note
that it plays very nicely with the notion of operator norm.

######  Lemma 2.2 (Lemma 4.4.1 in Vershynin (2018)).

Let $M\in\mathbb{R}^{m\times d}$ and let $\varepsilon\in(0,1)$. Then for any
$\varepsilon$-net $\mathcal{N}$ of $(\mathbb{S}^{d-1},\|\cdot\|_{2})$ we have
that:

$\|M\|_{\mathsf{op}}\leq\frac{1}{1-\varepsilon}\sup_{v\in\mathcal{N}}\|Mv\|_{2}.$
(2.10)

Hence at a small multiplicative cost, the computation of the operator norm can
be restricted to the discretized sphere $\mathcal{N}$. Our intention is now to
apply the union bound (2.8) to the right hand side of (2.10). To do so, we
also need control of the size (cardinality) of the $\varepsilon$-net.

######  Lemma 2.3 (Corollary 4.2.13 in Vershynin (2018)).

For any $\varepsilon>0$ the covering numbers of $\mathbb{S}^{d-1}$ satisfy

$\mathcal{N}(\varepsilon,\mathbb{S}^{d-1},\|\cdot\|)\leq\left(1+\frac{1}{2\varepsilon}\right)^{d}.$
(2.11)

We now provide two instances of this covering argument combined with the union
bound. The second of these uses an alternative variational characterization of
the operator norm but otherwise similar ideas.

######  Lemma 2.4.

Let $M$ be an $m\times d$ random matrix, and $\epsilon\in(0,1)$. Furthermore,
let $\mathcal{N}$ be an $\epsilon$-net of $\mathbb{S}^{d-1}$ of minimal
cardinality. Then for all $\rho>0$, we have

$\displaystyle\mathbf{P}\left(\|M\|_{\mathsf{op}}>\rho\right)$
$\displaystyle\leq\left(\frac{2}{\epsilon}+1\right)^{d}\max_{v\in\mathcal{N}}\mathbf{P}\left(\|Mv\|_{2}>(1-\epsilon)\rho\right).$

######  Lemma 2.5.

Let $M$ be an $d\times d$ symmetric random matrix, and let
$\epsilon\in(0,1/2)$. Furthermore, let $\mathcal{N}$ be an $\epsilon$-net of
$\mathbb{S}^{d-1}$ with minimal cardinality. Then for all $\rho>0$, we have

$\displaystyle\mathbf{P}\left(\|M\|_{\mathsf{op}}>\rho\right)$
$\displaystyle\leq\left(\frac{2}{\epsilon}+1\right)^{d}\max_{v\in\mathcal{N}}\mathbf{P}\left(|v^{\top}Mv|>(1-2\epsilon)\rho\right).$

Lemma 2.4 and Lemma 2.5 exploit two different variational forms of the
operator norm. Namely for any $M$ we have that
$\|M\|_{\mathsf{op}}^{2}=\sup_{v\in\mathbf{S}^{d-1}}\|Mv\|^{2}$ and in
addition, when $M$ is symmetric we also have, $\|M\|_{\mathsf{op}}=\sup_{v\in
S^{d-1}}|v^{\top}Mv|$. The proof of these last two lemmas are standard and can
be found for example in (Vershynin, 2018, Chapter 4).

###  2.3 Concentration of the Covariance Matrix of Linear Systems

To not get lost in the weeds, let us provide an example showcasing the use of
Theorem 2.1 due to Jedra and Proutiere (2022). Recall that the matrix
$\widehat{\Sigma}$ appearing in (1.9) is crucial to the performance of the
least squares estimator. We will now see that this matrix is well-conditioned
when we consider stable first order auto-regressions of the form:

$X_{t+1}=A^{\star}X_{t}+W_{t}\qquad t=1,\dots,T\qquad W_{1:T}\textnormal{ iid\
isotropic and $K^{2}$-subG}$ (2.12)

taking values in $\mathbb{R}^{d_{\mathsf{X}}}$. By stable we mean that the
largest eigenvalue of $A^{\star}$ has module strictly smaller than 1.

The following result is a consequence of the Hanson-Wright inequality together
with the discretization strategy outlined in Section 2.2. The full proof is
given in Appendix B.

######  Theorem 2.2.

Let $\varepsilon>0$ and set
$M\triangleq\left(\sum_{t=1}^{T}\sum_{k=0}^{t-1}(A^{\star})^{k}(A^{\star,\mathsf{T}})^{k}\right)^{-\frac{1}{2}}$.
Let also $\mathbf{L}$ be the linear operator such that
$X_{1:T}=\mathbf{L}W_{1:T}$. Then simultaneously for every
$i\in[d_{\mathsf{X}}]$:

$(1-\varepsilon)^{2}\lambda_{\min}\left(\sum_{t=1}^{T}\sum_{k=0}^{t-1}(A^{\star})^{k}(A^{\star,\mathsf{T}})^{k}\right)\leq\lambda_{i}\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)\leq(1+\varepsilon)^{2}\lambda_{\max}\left(\sum_{t=1}^{T}\sum_{k=0}^{t-1}(A^{\star})^{k}(A^{\star,\mathsf{T}})^{k}\right)$

holds with probability at least

$1-\exp\left(-\frac{\varepsilon^{2}}{576\,K^{2}\|M\|_{\mathsf{op}}^{2}\|\mathbf{L}\|_{\mathsf{op}}^{2}}+d_{\mathsf{X}}\log(18)\right).$
(2.13)

Put differently, on the same event as in Theorem 2.2, the spectrum of

$\widehat{\Sigma}=\frac{1}{T}\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}$ (2.14)

is sandwiched by that of its population counterpart
($\mathbf{E}\widehat{\Sigma}$) within a small multiplicative factor. The
result holds with high probability for strictly stable systems.

The quantity $\|L\|_{\mathsf{op}}$ in (2.13) grows very quickly as the
spectral radius of $A^{\star}$ tends to $1$; Theorem 2.2 becomes vacuous in
the marginally stable regime. It turns out that requirement of two-sided
concentration—the sandwiching of the entire spectrum—is too stringent a
requirement to obtain bounds that degrade gracefully with the stability of the
system. Fortunately, we only need sharp control of the lower half of the
spectrum to control the error (1.10). This motivates Section 3 below, in which
we will see how to relax the stability assumption and analyze more general
linear systems.

###  2.4 Notes

The basic program carried out in Section 2.3 can be summarized as follows: (1)
introduce a discretization of the problem considered—for matrices this is
typically a discretization of the unit sphere; (2) prove an exponential
inequality for a family of scalar random variables corresponding to one-
dimensional projection of the discretization—in our case: prove bounds on the
moment generating function of quadratic forms in random matrices; and (3)
conclude to obtain a uniform bound by using the union bound across the
discretization. This roughly summarizes the proof of Theorem 2.2. These tools
are thematic throughout this manuscript.

##  3 The Lower Spectrum of the Empirical Covariance

Recall that our outline of the analysis of the least squares estimator in
Section 1.2 consists of two main components, one of which being the lower tail
of the empirical covariance matrix (1.9). In this section we provide a self-
contained analysis of this random matrix for a class of "causal" systems.
Moreover, we will emphasize only the lower tail of this random matrix as to
sidestep issues with bounds degrading with the stability of the system
considered. This allows us to quantitatively separate the notions of
persistence of excitation and stability.

Let us now carry out this program. Fix two integers $T$ and $k$ such that
$T/k\in\mathbb{N}$. We consider causal processes of the form
$X_{1:T}=(X_{1}^{\mathsf{T}},\dots,X_{T}^{\mathsf{T}})^{\mathsf{T}}$ evolving
on $\mathbb{R}^{d}$. More precisely, we assume the existence of an isotropic
sub-Gaussian process evolving on $\mathbb{R}^{p}$, $W_{1:T}$ with
$\mathbf{E}W_{1:T}W_{1:T}^{\mathsf{T}}=I_{pT}$ and a (block-) lower-triangular
matrix $\mathbf{L}\in\mathbb{R}^{dT\times pT}$ such that

$X_{1:T}=\mathbf{L}W_{1:T}.$ (3.1)

We will assume that all the $pT$-many entries of $W_{1:T}$ are independent
$K^{2}$-sub-Gaussian for some positive $K\in\mathbb{R}$.

We say that $X_{1:T}$ is $k$-causal if the matrix $\mathbf{L}$ has the block
lower-triangular form:

$\displaystyle\mathbf{L}=\begin{bmatrix}\mathbf{L}_{1,1}&0&0&0&0\\\
\mathbf{L}_{2,1}&\mathbf{L}_{2,2}&0&0&0\\\
\mathbf{L}_{3,1}&\mathbf{L}_{3,2}&\mathbf{L}_{3,3}&0&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
\mathbf{L}_{T/k,1}&\dots&\dots&\dots&\dots\mathbf{L}_{T/k,T/k}\end{bmatrix}=\begin{bmatrix}\mathbf{L}_{1}\\\
\mathbf{L}_{2}\\\ \mathbf{L}_{3}\\\ \vdots\\\ \mathbf{L}_{T/k}\end{bmatrix}$
(3.2)

where each $\mathbf{L}_{ij}\in\mathbb{R}^{dk\times
pk},i,j\in[T/k]\triangleq\\{1,2,\dots,T/k\\}$. In brief, we say that
$X_{0:T-1}$ satisfying the above construction is $k$-causal with independent
$K^{2}$-sub-Gaussian increments.

Obviously, every $1$-causal process is $k$-causal for every $k\in\mathbb{N}$
as long as the divisibility condition holds. To analyze the lower tail of the
empirical covariance of $X_{0:T-1}$ we will also associate a decoupled random
process

$\tilde{X}_{1:T}=\mathrm{blkdiag}(\mathbf{L}_{11},\dots,\mathbf{L}_{T/k,T/k})W_{1:T}.$

Hence, the process $\tilde{X}_{1:T}$ is generated in much the same way as
$X_{1:T}$ but by removing the sub-diagonal entries of $\mathbf{L}$:

$\tilde{\mathbf{L}}\triangleq\begin{bmatrix}\mathbf{L}_{1,1}&0&0&0\\\
0&\mathbf{L}_{2,2}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&0\\\
0&\dots&0&\mathbf{L}_{T/k,T/k}\end{bmatrix}\quad\Longrightarrow\quad\tilde{X}_{1:T}=\tilde{\mathbf{L}}W_{1:T}.$

We emphasize that by our assumptions on $W_{1:T}$ and the block-diagonal
structure of $\tilde{\mathbf{L}}$ the variables
$\tilde{X}_{1:k},\tilde{X}_{k+1:2k},\dots,\tilde{X}_{T-k+1:T}$ are all
independent of each other; they have been decoupled. This decoupled process
will effectively dictate our lower bound, and we will show under relatively
mild assumptions that

$\displaystyle\lambda_{\min}\left(\frac{1}{T}\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)\gtrsim\lambda_{\min}\left(\frac{1}{T}\sum_{t=1}^{T}\mathbf{E}\tilde{X}_{t}\tilde{X}_{t}^{\mathsf{T}}\right)$
(3.3)

with probability that approaches $1$ at an exponential rate in the sample size
$T$. More precisely, the following statement is the main result of this
section.

######  Theorem 3.1.

Fix an integer $k\in\mathbb{N}$, let $T\in N$ be divisible by $k$ and suppose
$X_{1:T}$ is a $k$-causal process taking values in $\mathbb{R}^{d}$ with
$K^{2}$-sub-Gaussian increments. Suppose further that the diagonal blocks are
all equal: $\mathbf{L}_{j,j}=\mathbf{L}_{1,1}$ for all $j\in[T/k]$. Suppose
$\lambda_{\min}\left(\sum_{t=1}^{T}\mathbf{E}\tilde{X}_{t}\tilde{X}_{t}^{\mathsf{T}}\right)>0$.
We have that:

$\mathbf{P}\left(\frac{1}{T}\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\nsucceq\frac{1}{8T}\sum_{t=1}^{T}\mathbf{E}\tilde{X}_{t}\tilde{X}_{t}^{\mathsf{T}}\right)\leq\left(C_{\mathsf{sys}}\right)^{d}\exp\left(-\frac{T}{576K^{2}k}\right)$
(3.4)

where

$C_{\mathsf{sys}}\triangleq
1+2\sqrt{2}\frac{\left(\frac{T\|\mathbf{L}\mathbf{L}^{\mathsf{T}}\|_{\mathsf{op}}}{18k\lambda_{\min}\left(\sum_{t=1}^{T}\mathbf{E}X_{t}X_{t}^{\mathsf{T}}\right)}+9\right)\lambda_{\max}\left(\sum_{t=1}^{T}\mathbf{E}X_{t}X_{t}^{\mathsf{T}}\right)}{\lambda_{\min}\left(\sum_{t=1}^{T}\mathbf{E}\tilde{X}_{t}\tilde{X}_{t}^{\mathsf{T}}\right)}.$
(3.5)

To parse Theorem 3.1, note that it simply informs us that there exist a a
system-dependent constant $C_{\mathsf{sys}}$—which itself has no more than
polynomial dependence on relevant quantities—such that if

$T/k\geq 576K^{2}(d\log C_{\mathsf{sys}}+\log(1/\delta))$ (3.6)

then on an event with probability mass at least $1-\delta$:

$\frac{1}{T}\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\succeq\frac{1}{8T}\sum_{t=1}^{T}\mathbf{E}\tilde{X}_{t}\tilde{X}_{t}^{\mathsf{T}}.$

######  Remark 3.1.

Since the blocks of $\mathbf{L}$ can be regarded to specify the noise-to-
output map, the assumption that the diagonal blocks are constant is for
instance satisfied by linear time-invariant (LTI) systems. The assumption can
be removed at the cost of a more complicated expression.

The next example serves as the archetype for the reduction from $\mathbf{L}$
to $\tilde{\mathbf{L}}$.

######  Example 3.1.

Suppose that (3.1) is specified via

$X_{t}=A_{\star}X_{t-1}+B_{\star}W_{t}$ (3.7)

for $t\in[T]$ and where
$(A_{\star},B_{\star})\in\mathbb{R}^{d_{\mathsf{X}}\times
d_{\mathsf{X}}+d_{\mathsf{X}}\times d_{\mathsf{W}}}$. We set
$d=d_{\mathsf{X}}$ and $p=d_{\mathsf{W}}$ in the theorem above. The reduction
from $X_{1:T}=\mathbf{L}W_{1:T}$ to
$\tilde{X}_{1:T}=\mathrm{blkdiag}(\mathbf{L}_{11},\dots,\mathbf{L}_{T/k,T/k})W_{1:T}$
corresponds to replacing a single trajectory from the linear system (3.7) of
length $T$ by $T/k$ trajectories of length $k$ each and sampled independently
of each other. The price we pay for decoupling these systems is that our lower
bound is dictated by the gramians up to range $k$:

$\frac{1}{T}\sum_{t=1}^{T}\mathbf{E}\tilde{X}_{t}\tilde{X}_{t}^{\mathsf{T}}=\frac{1}{k}\sum_{t=1}^{k}\mathbf{E}\tilde{X}_{t}\tilde{X}_{t}^{\mathsf{T}}=\frac{1}{k}\sum_{t=1}^{k}\sum_{j=0}^{t-1}(A^{\star})^{j}B^{\star}B^{\star,\mathsf{T}}(A^{\star,\mathsf{T}})^{j}$
(3.8)

instead of the gramians up to range $T$:

$\frac{1}{T}\sum_{t=1}^{T}\mathbf{E}X_{t}X_{t}^{\mathsf{T}}=\frac{1}{T}\sum_{t=1}^{T}\sum_{j=0}^{t-1}(A^{\star})^{j}B^{\star}B^{\star,\mathsf{T}}(A^{\star,\mathsf{T}})^{j}.$
(3.9)

Put differently, the reduction from $\mathbf{L}$ to $\tilde{\mathbf{L}}$ can
be thought of as restarting the system every $k$ steps.

Comparing with Theorem 2.2, the advantage of Theorem 3.1 is that it allows us
to provide persistence-of-excitation type guarantees that do not rely strongly
on the stability of the underlying system. While Theorem 2.2 gives in
principle stronger two-sided concentration results, it comes at the cost of
the guarantees becoming vacuous as the spectral radius of $A^{\star}$ in
Example 3.1 tends to marginal stability (tends to $1$). By contrast, Theorem
3.1 does not exhibit such a blow-up since the dependence on $C_{\mathsf{sys}}$
in (3.6) is logarithmic (instead of polynomial). The distinction might seem
small, but it is qualitatively important as it (almost) decouples the
phenomena of stability and persistence of excitation.

###  3.1 A Decoupling Inequality for sub-Gaussian Quadratic Forms

Our proof of Theorem 3.1 will make heavy use of Proposition 3.1 below. This is
the crucial probalistic inequality that allows us to decouple—or restart as
discussed in Example 3.1.

######  Proposition 3.1.

Fix $K\geq 1$, $x\in\mathbb{R}^{n}$ and a symmetric positive semidefinite
$Q\in\mathbb{R}^{(n+m)\times(n+m)}$ of the form $\displaystyle
Q=\begin{bmatrix}Q_{11}&Q_{12}\\\ Q_{21}&Q_{22}\end{bmatrix}$ with
$Q_{22}\succ 0$. Let $W$ be an $m$-dimensional mean zero, isotropic and
$K^{2}$-sub-Gaussian random vector with independent entries. Then for every
$\lambda\in\left[0,\frac{1}{8\sqrt{2}K^{2}\|Q_{22}\|_{\mathsf{op}}}\right]$ it
holds true that:

$\displaystyle\mathbf{E}\exp\left(-\lambda\begin{bmatrix}x\\\
W\end{bmatrix}^{\mathsf{T}}\begin{bmatrix}Q_{11}&Q_{12}\\\
Q_{21}&Q_{22}\end{bmatrix}\begin{bmatrix}x\\\
W\end{bmatrix}\right)\leq\exp\left(-\lambda\operatorname{tr}Q_{22}+36K^{4}\lambda^{2}\operatorname{tr}Q_{22}^{2}\right).$
(3.10)

By combining Lemma 3.1 below with the exponential form of Hanson-Wright
(Proposition A.1) we obtain the exponential inequality (3.10), which in the
sequel will allow us to control the lower tail of the conditionally random
quadratic form

$\begin{bmatrix}x\\\
W\end{bmatrix}^{\mathsf{T}}\begin{bmatrix}Q_{11}&Q_{12}\\\
Q_{21}&Q_{22}\end{bmatrix}\begin{bmatrix}x\\\ W\end{bmatrix}.$

We point out that (3.10) is not the best possible if the entries of are $W$
independent and Gaussian as opposed to just isotropic and sub-Gaussian. In
this case, the factor $36K^{4}\lambda^{2}(\operatorname{tr}Q_{22})^{2}$ in
(3.10) can be improved to $\frac{\lambda^{2}}{2}\operatorname{tr}Q_{22}^{2}$
and the inequality can be shown to hold for the entire range of non-negative
$\lambda$ (Ziemann, 2022, Lemma 2.1). Irrespectively, we will see in the
sequel that it captures the correct qualitative behavior.

######  Lemma 3.1 (sub-Gaussian Decoupling).

Fix $K\geq 1$, $x\in\mathbb{R}^{n}$ and a symmetric positive semidefinite
$Q\in\mathbb{R}^{(n+m)\times(n+m)}$ of the form $\displaystyle
Q=\begin{bmatrix}Q_{11}&Q_{12}\\\ Q_{21}&Q_{22}\end{bmatrix}$. Let $W$ be an
$m$-dimensional mean zero and $K^{2}$-sub-Gaussian random vector. Then for
every $\lambda\in\left[0,\frac{1}{4K^{2}\|Q_{22}\|_{\mathsf{op}}}\right]$ it
holds true that:

$\displaystyle\mathbf{E}\exp\left(-\lambda\begin{bmatrix}x\\\
W\end{bmatrix}^{\mathsf{T}}\begin{bmatrix}Q_{11}&Q_{12}\\\
Q_{21}&Q_{22}\end{bmatrix}\begin{bmatrix}x\\\
W\end{bmatrix}\right)\leq\sqrt{\mathbf{E}\exp\left(-2\lambda
W^{\mathsf{T}}Q_{22}W\right)}.$ (3.11)

Once equipped with (3.11), Proposition 3.1 follows immediately by Proposition
A.1. The proof of Lemma 3.1 is given in Appendix C.

###  3.2 The Lower Tail of the Empirical Covariance of Causal sub-Gaussian
Processes

Repeated application of Proposition 3.1 to the process
$X_{1:T}=\mathbf{L}W_{1:T}$ in combination with the tower property of
conditional expectation yields the following exponential inequality that
controls the lower tail of (1.9) in any fixed direction.

######  Theorem 3.2.

Fix an integer $k\in\mathbb{N}$, let $T\in N$ be divisible by $k$ and suppose
$X_{1:T}$ is a $k$-causal process driven by independent $K^{2}$-sub-Gaussian
increments as described in Section 3. Fix also a matrix
$\Delta\in\mathbb{R}^{d^{\prime}\times d}$. Let
$Q_{\max}\triangleq\max_{j\in[T/k]}\|\mathbf{L}_{j,j}^{\mathsf{T}}\mathrm{blkdiag}(\Delta^{\mathsf{T}}\Delta)\mathbf{L}_{j,j}\|_{\mathsf{op}}$
Then for every $\lambda\in\left[0,\frac{1}{8\sqrt{2}K^{2}Q_{\max}}\right]$:

$\mathbf{E}\exp\left(-\lambda\sum_{t=1}^{T}\|\Delta X_{t}\|_{2}^{2}\right)\\\
\leq\exp\Bigg{(}-\lambda\sum_{j=1}^{T/k}\operatorname{tr}\left(\mathbf{L}_{j,j}^{\mathsf{T}}\mathrm{blkdiag}(\Delta^{\mathsf{T}}\Delta)\mathbf{L}_{j,j}\right)+36K^{4}\lambda^{2}\sum_{j=1}^{T/k}\operatorname{tr}\left(\mathbf{L}_{j,j}^{\mathsf{T}}\mathrm{blkdiag}(\Delta^{\mathsf{T}}\Delta)\mathbf{L}_{j,j}\right)^{2}\Bigg{)}.$

To appreciate the terms appearing in Theorem 3.2, it is worth to point out
that

$\sum_{j=1}^{T/k}\operatorname{tr}\left(\mathbf{L}_{j,j}^{\mathsf{T}}\mathrm{blkdiag}(\Delta^{\mathsf{T}}\Delta)\mathbf{L}_{j,j}\right)=\sum_{t=1}^{T}\mathbf{E}\|\Delta\tilde{X}_{t}\|_{2}^{2}.$

Hence Theorem 3.2 effectively passes the expectation inside the exponential at
the cost of working with the possibly less excited process $\tilde{X}_{1:T}$
and a quadratic correction term. Note also that the assumption that $T$ is
divisible by $k$ is not particularly important. If not, let $T^{\prime}$ be
the largest integer such that $T^{\prime}/k\in\mathbb{N}$ and $T^{\prime}\leq
T$ and apply the result with $T^{\prime}$ in place of $T$.

The significance of Theorem 3.2 is demonstrated by the following simple
observation, which is just the Chernoff approach applied to the exponential
inequality in Theorem 3.2.

######  Lemma 3.2.

Fix an integer $k\in\mathbb{N}$, let $T\in N$ be divisible by $k$ and suppose
$X_{1:T}$ is a $k$-causal process with independent $K^{2}$-sub-Gaussian
increments. Suppose further that the diagonal blocks are all equal:
$\mathbf{L}_{j,j}=\mathbf{L}_{1,1}$ for all $j\in[T/k]$. For every size-
conforming matrix $\Delta$ we have that:

$\mathbf{P}\left(\sum_{t=1}^{T}\|\Delta
X_{t}\|^{2}_{2}\leq\frac{1}{2}\sum_{t=1}^{T}\mathbf{E}\|\Delta\tilde{X}_{t}\|^{2}_{2}\right)\leq\exp\left(-\frac{T}{576K^{2}k}\right).$
(3.12)

Note that Lemma 3.2 only yields _pointwise_ control of the empirical
covariance—i.e. pointwise on the sphere $\mathbb{S}^{d-1}$ . By setting
$\Delta=v\in\mathbb{S}^{d-1}$, the result holds for a fixed vector on the
sphere, but not uniformly for all such vectors at once. Thus, returning to our
over-arching goal of providing control of the smallest eigenvalue of the
empirical covariance matrix (1.9), we now combine (3.12) (using
$d^{\prime}=1$) with a union bound. This approach yields Theorem 3.1, of which
the proof—along with that of its supporting lemmas—is given in full in
Appendix C.222Similar results can also be obtained for restricted eigenvalues.

###  3.3 Notes

In this manuscript we have chosen a perhaps less well-known but conceptually
simpler approach to establishing lower bounds on the empirical covariance
matrix Equation 3.3. The first proof of a statement similar to Theorem 3.1 is
due to Simchowitz et al. (2018) which in turn relies on a more advanced notion
from probability theory known as the small-ball method, due to Mendelson
(2014). The emphasis therein is on anti-concentration—which can hold under
milder moment assumptions—rather than concentration. However, the introduction
of this tool is not necessary for Gaussian (or sub-Gaussian) system
identification. For instance, Sarkar and Rakhlin (2019) leverage the method of
self-normalized martingales introduced in Section 4 below.

Our motivation for providing a different proof is to streamline the exposition
as to fit control of the lower tail into the "standard machinery", which
roughly consists of: (1) prove a family of scalar exponential inequalities,
(2) invoke the Chernoff method, and (3) conclude by a discretization argument
and a union bound to port the result from scalars to matrices. Our proof here
follows this outline and emphasizes the exponential inequality in Theorem 3.2.
We finally remark that the proof presented here is new to the literature and
extends a result in Ziemann (2022) from the Gaussian setting to the sub-
Gaussian setting.

##  4 Self-Normalized Martingale Bounds

The objective in this section is to bound the operator and Frobenius norms of
the self-normalized term of (1.10):

$\displaystyle\left(\sum_{t=1}^{T}V_{t}X_{t}^{\mathsf{T}}\right)\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{-1/2}.$
(4.1)

This object has special structure. Firstly, in many cases of interest, e.g.
the autoregressive model in (1.2), the noise $V_{t}$ is independent of $X_{k}$
for all $k\leq t$. This is what provides martingale structure, as will be made
precise shortly. Secondly, it is self-normalized: if the covariates $X_{t}$
are large for some $t$, then any increase in the left sum will be compensated
by an increase in the sum in the term on the right. Together, these properties
make the object above a self-normalized martingale term.

To express results generally and compactly, several definitions are in order.

######  Definition 4.1.

(Filtration and Adapted Process) A sequence of sub-$\sigma$-algebras
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=1}^{T}$ is said to
be a _filtration_ if $\mathcal{F}_{t}\subseteq\mathcal{F}_{k}$ for $t\leq k$.
A stochastic process $\mathopen{}\left\\{W_{t}\right\\}\mathclose{}_{t=1}^{T}$
is said to be _adapted_ to the filtration
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=1}^{T}$ if for all
$t\geq 1$, $W_{t}$ is $\mathcal{F}_{t}$-measurable.

Conditioning on a sub-$\sigma$-algebra provides partial information about the
total randomness. Therefore, the requirement that a filtration is non-
decreasing captures the fact that information is not forgotten. An adapted
process is one in which all the randomness at a particular time is explained
by the information in the filtration up to that time.

######  Definition 4.2.

(Martingale) Consider a stochastic process
$\mathopen{}\left\\{W_{t}\right\\}\mathclose{}_{t=1}^{T}$ which is adapted to
a filtration
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=1}^{T}$. This
process is called a _martingale_ if for all $1\leq t\leq T$, $W_{t}$ is
integrable and for all $1\leq t<T$,
$\mathbf{E}[W_{t+1}|\mathcal{F}_{t}]=W_{t}$.

Martingales model causal or non-anticipative processes. To better appreciate
this, note that the increments $W_{t+1}-W_{t}$ are mean zero and conditionally
orthogonal to the past; they can be thought of as the "next step" in a random
walk whose path is traced out by $W_{t}$.

In the context of the linear time-series model in (1.1), we may define the
sub-$\sigma$-algebras $\mathcal{F}_{t}$ as those induced by the randomness up
to time $t$: $\mathcal{F}_{t}=\sigma(X_{1},\dots,X_{t+1},V_{1},\dots,V_{t})$.
In this case, the process
$\mathopen{}\left\\{X_{t}\right\\}\mathclose{}_{t=1}^{T}$ is adapted to the
filtration
$\mathopen{}\left\\{\mathcal{F}_{t-1}\right\\}\mathclose{}_{t=1}^{T}$ and the
process $\mathopen{}\left\\{V_{t}\right\\}\mathclose{}_{t=1}^{T}$ is adapted
to the filtration
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=1}^{T}$.

Recall now again that the "numerator" in the ordinary least squares error is
(4.1). We see that if we define the sum,
$S_{t}\triangleq\sum_{s=1}^{t}V_{s}X_{s}^{\top}$, then the process
$\mathopen{}\left\\{S_{t}\right\\}\mathclose{}_{t=1}^{T}$ is adapted to
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=1}^{T}$.
Furthermore,
$\mathbf{E}\left(S_{t+1}|\mathcal{F}_{t}\right)=S_{t}+\mathbf{E}\left(V_{t+1}|\mathcal{F}_{t}\right)X_{t+1}^{\top}$.
In particular, as long as the noise has conditional mean zero
($\mathbf{E}\left(V_{t+1}|\mathcal{F}_{t}\right)=0$), the process
$\mathopen{}\left\\{S_{t}\right\\}\mathclose{}_{t=1}^{T}$ is a
martingale.333Indeed, $S_{t}$ is a so-called martingale transform of
$X_{1:t}$. Normalizing the sum $S_{t}$ by the covariates as
$S_{t}\left(\sum_{s=1}^{t}X_{s}X_{s}^{\top}\right)^{-1/2}$ _almost_ preserves
the martingale structure. Expressions of this type are called self-normalized
martingales—although we stress that they are not strictly speaking martingales
but only constructed from them.

We now state bounds on the operator and Frobenius norms of the self-normalized
martingale. The main idea behind the result—the technique of pseudo-
maximization—is due to Robbins and Siegmund (1970). The formulations presented
here are a consequence of Theorem 3.4 in Abbasi-Yadkori (2013).

######  Theorem 4.1.

(Special cases of Theorem 3.4 in Abbasi-Yadkori (2013)) Let
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=0}^{T}$ be a
filtration such that $\mathopen{}\left\\{X_{t}\right\\}\mathclose{}_{t=1}^{T}$
is adapted to
$\mathopen{}\left\\{\mathcal{F}_{t-1}\right\\}\mathclose{}_{t=1}^{T}$ and
$\mathopen{}\left\\{V_{t}\right\\}\mathclose{}_{t=1}^{T}$ is adapted to
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=1}^{T}$.
Additionally, suppose that for all $1\leq t\leq T$, $V_{t}$ is
$\sigma^{2}$-conditionally sub-Gaussian with respect to $\mathcal{F}_{t}$. Let
$\Sigma$ be a positive definite matrix in $\mathbb{R}^{d_{\mathsf{X}}\times
d_{\mathsf{X}}}$. For a fixed $T\in\mathbb{N}_{+}$ and $\delta\in(0,1)$, with
probability at least $1-\delta$,

$\displaystyle\left\lVert\sum_{t=1}^{T}V_{t}X_{t}^{\top}\left(\Sigma+\sum_{t=1}^{T}X_{t}X_{t}^{\top}\right)^{-1/2}\right\rVert_{F}^{2}\leq
d_{\mathsf{Y}}\sigma^{2}\log\left(\frac{\det\left(\Sigma+\sum_{t=1}^{T}X_{t}X_{t}^{\top}\right)}{\det(\Sigma)}\right)+2\sigma^{2}\log\frac{1}{\delta}.$

Additionally, for a fixed $T\in\mathbb{N}_{+}$ and $\delta\in(0,1)$, with
probability at least $1-\delta$,

$\displaystyle\left\|\sum_{t=1}^{T}V_{t}X_{t}^{\top}\left(\Sigma+\sum_{t=1}^{T}X_{t}X_{t}^{\top}\right)^{-1/2}\right\|_{\mathsf{op}}^{2}$
$\displaystyle\qquad\leq
4\sigma^{2}\log\left(\frac{\det\left(\Sigma+\sum_{t=1}^{T}X_{t}X_{t}^{\top}\right)}{\det(\Sigma)}\right)+8d_{\mathsf{Y}}\sigma^{2}\log
5+8\sigma^{2}\log\frac{1}{\delta}.$

Note that the quantities bounded above have a positive definite matrix
$\Sigma$ added to the normalization quantities that was not present in the
original term of interest, (4.1). Furthermore, the covariates
$\sum_{t=1}^{T}X_{t}X_{t}$ appear in the upper bound. Hence, one typically
combines the self-normalized martingale bound with some weak form of
concentration.444Alternatively, in the analysis of ridge regression, $\Sigma$
takes the role of the penalizing matrix which can be tuned. This is done in
Section 5.

In the sequel, we prove the above bounds. To obtain the bound on the Frobenius
norm, we directly consider the object

$\displaystyle\left\lVert\left(\sum_{t=1}^{T}V_{t}X_{t}^{\mathsf{T}}\right)\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{-1/2}\right\rVert_{F},$
(4.2)

while to obtain the bound on the operator norm we consider the following
vector norm for an arbitrary unit vector $w$ as an intermediate step:

$\displaystyle\left\lVert\left(w^{\top}\sum_{t=1}^{T}V_{t}X_{t}^{\mathsf{T}}\right)\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{-1/2}\right\rVert_{2}$
(4.3)

and combine with a covering argument (recall Section 2.2).

Let us also briefly consider the dimensional dependencies of the Frobenius and
operator norm bounds in Theorem 4.1. The leading term in the Frobenius norm
bound is $d_{\mathsf{Y}}$ multiplied by the $\log\det$ term, which scales with
$d_{\mathsf{X}}\log T$ when the empirical covariance is well-conditioned. In
particular, the leading term scales with $d_{\mathsf{X}}d_{\mathsf{Y}}\log T$.
The factor of $d_{\mathsf{Y}}$ is no longer present on the $\log\det$ term for
the operator norm. The term therefore scales as $d_{\mathsf{X}}\log T$ when
the the empirical covariance matrix is well-conditioned. There is, however, a
term $8d_{\mathsf{Y}}\sigma^{2}\log 5$ which results from the covering
argument. The operator norm bound therefore scales as
$\max\mathopen{}\left\\{d_{\mathsf{X}}\log
T,d_{\mathsf{Y}}\right\\}\mathclose{}$.

###  4.1 Exponential Inequalities via Pseudo-maximization

We begin by neglecting the details of the process that generated the data in
(4.1). In particular, consider a random matrix $P$ assuming values in
$\mathbb{R}^{d_{\mathsf{\eta}}\times d_{\mathsf{X}}}$ for
$d_{\mathsf{\eta}}\in\mathbb{N}_{+}$ and a random matrix $Q$ assuming values
in $\mathbb{R}^{d_{\mathsf{X}}\times d_{\mathsf{X}}}$ with $Q$ almost surely
nonsingular. Bounding the quantities in (4.2) and (4.3) are special cases of
bounding $\lVert PQ^{-\frac{1}{2}}\rVert_{F}$. A naive first approach is to
apply a Chernoff bound (2.2). Doing so results in the inequality

$\displaystyle\mathbf{P}\left(\lVert PQ^{-1/2}\rVert_{F}\geq
x\right)\leq\min_{\lambda\geq
0}\exp\left(-\frac{\lambda}{2}x^{2}\right)\mathbf{E}\exp\left(\frac{\lambda}{2}\lVert
PQ^{-1/2}\rVert_{F}^{2}\right).$

If it is possible to bound the moment generating function
$\mathbf{E}\exp\left(\frac{\lambda}{2}\lVert PQ^{-1/2}\rVert_{F}^{2}\right)$
by one for some $\lambda>0$, then the above bound provides an exponential
inequality. Obtaining a bound of the form
$\mathbf{E}\exp\left(\frac{\lambda}{2}\lVert
PQ^{-1/2}\rVert_{F}^{2}\right)\leq 1$ requires very strong assumptions on $P$
and $Q$ which would not be suitable for our purposes. However, we may observe
that $\frac{1}{2}\lVert
PQ^{-1/2}\rVert_{F}^{2}=\max_{\Lambda}\operatorname{tr}(P\Lambda-\frac{1}{2}\Lambda^{\top}Q\Lambda)$.
This motivates the following canonical assumption in self-normalized process
theory:

$\displaystyle\max_{\Lambda\in\mathbb{R}^{d_{\mathsf{X}}\times
d_{\mathsf{\eta}}}}\mathbf{E}\exp\operatorname{tr}\left(P\Lambda-\frac{1}{2}\Lambda^{\top}Q\Lambda\right)\leq
1.$ (4.4)

This inequality is called the canonical assumption because a wide variety of
self-normalized processes satisfy it. We will demonstrate that it is satisfied
for (4.2) and (4.3) in Section 4.2. If we could exchange the order of the
maximization with the expectation in (4.4), then the bound
$\mathbf{E}\exp\left(\frac{1}{2}\lVert PQ^{-1/2}\rVert_{F}^{2}\right)\leq 1$
would be satisfied, and the Chernoff bound above would provide a valuable
exponential inequality. As this exchange is not possible, we instead lower
bound the maximization over $\Lambda$ by assigning a probability distribution
to a random variable $\Psi$ which takes values in
$\mathbb{R}^{d_{\mathsf{X}}\times d_{\mathsf{\eta}}}$, and taking the
expectation over this distribution. Doing so preserves the inequality (4.4):

$\displaystyle\mathbf{E}\mathbf{E}\left[\exp\operatorname{tr}\left(P\Psi-\frac{1}{2}\Psi^{\top}Q\Psi\right)\Big{|}\Psi\right]\leq
1.$

The order of expectation over $\Psi$ and over the random variables $P$ and $Q$
may then be exchanged by an appeal to Fubini’s theorem:

$\displaystyle
1\geq\mathbf{E}\mathbf{E}\left[\exp\operatorname{tr}\left(P\Psi-\frac{1}{2}\Psi^{\top}Q\Psi\right)\Big{|}\Psi\right]=\mathbf{E}\mathbf{E}\left[\exp\operatorname{tr}\left(P\Psi-\frac{1}{2}\Psi^{\top}Q\Psi\right)\Big{|}P,Q\right].$
(4.5)

By selecting the distribution over $\Psi$ appropriately, the result is a so-
called pseudo-maximization. In particular, by completing the square, the inner
conditional expectation on the right may be expressed as

$\displaystyle\mathbf{E}\left[\exp\operatorname{tr}\left(P\Psi-\frac{1}{2}\Psi^{\top}Q\Psi\right)\Big{|}P,Q\right]$
(4.6)
$\displaystyle=\exp\operatorname{tr}\left(PQ^{-1}P^{\top}/2\right)\mathbf{E}\left[\exp\operatorname{tr}\left(-\frac{1}{2}(\Psi-Q^{-1}P^{\top})^{\top}Q(\Psi-Q^{-1}P^{\top})\right)\Big{|}P,Q\right].$

For particular choices of the distribution of $\Psi$, the right side of the
above expression approximates the maximimum value,
$\exp\operatorname{tr}\left(PQ^{-1}P^{\top}/2\right)$, of
$\exp\operatorname{tr}(P\Lambda-\frac{1}{2}\Lambda^{\top}Q\Lambda)$. This
allows us to apply a Chernoff argument similar to the one sketched above to
obtain an exponential bound on a quantity related to $\lVert
PQ^{-1/2}\rVert_{F}$. The following lemma demonstrates one such bound that
results by selecting the distribution of $\Psi$ as a matrix normal
distribution.

######  Lemma 4.1 (Extension of Theorem 14.7 in Peña et al. (2009)).

Suppose that (4.4) is satisfied. Let $\Sigma$ be a positive definite matrix in
$\mathbb{R}^{d_{\mathsf{X}}\times d_{\mathsf{X}}}$. Then, for $\delta>0$, with
probability at least $1-\delta$,

$\displaystyle\lVert P(Q+\Sigma)^{-1/2}\rVert_{F}^{2}\leq
2\log\left(\frac{\det(Q+\Sigma)^{d_{\mathsf{\eta}}/2}\det(\Sigma)^{-d_{\mathsf{\eta}}/2}}{\delta}\right).$

###  4.2 Self-Normalized Martingales Satisfy the Canonical Assumption

In order to make use of Lemma 4.1 to bound (4.2) or (4.3), we must ensure that
the condition (4.4) holds for

$\displaystyle P=\sum_{t=1}^{T}\frac{\eta_{t}X_{t}^{\top}}{\sigma},\quad
Q=\sum_{t=1}^{T}X_{t}X_{t}^{\top},$

where $\eta_{t}\in\mathbb{R}^{d_{\mathsf{\eta}}}$ is either the noise process
$V_{t}$ or the scalar process $w^{\top}V_{t}$ for some fixed unit vector $w$.
The following lemma shows that it is sufficient for $\eta_{t}$ to be
$\sigma^{2}$-conditionally sub-Gaussian.

######  Lemma 4.2.

Fix $T\in\mathbb{N}_{+}$. Let
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=0}^{T}$ be a
filtration such that $\mathopen{}\left\\{X_{t}\right\\}\mathclose{}_{t=1}^{T}$
is adapted to
$\mathopen{}\left\\{\mathcal{F}_{t-1}\right\\}\mathclose{}_{t=1}^{T}$ and
$\mathopen{}\left\\{\eta_{t}\right\\}\mathclose{}_{t=1}^{T}$ is adapted to
$\mathopen{}\left\\{\mathcal{F}_{t}\right\\}\mathclose{}_{t=1}^{T}$.
Additionally, suppose that for all $t\geq 1$, $\eta_{t}$ is
$\sigma^{2}$-conditionally sub-Gaussian with respect to $\mathcal{F}_{t}$. Let
$\Lambda\in\mathbb{R}^{d_{\mathsf{X}}\times d_{\mathsf{\eta}}}$ be arbitrary
and consider for $t\in\mathopen{}\left\\{1,\dots,T\right\\}\mathclose{}$

$\displaystyle
M_{t}(\Lambda)\triangleq\exp\operatorname{tr}\left(\sum_{s=1}^{t}\mathopen{}\left[\frac{\eta_{s}X_{s}^{\top}\Lambda}{\sigma}-\frac{1}{2}\Lambda^{\top}X_{s}X_{s}^{\top}\Lambda\right]\mathclose{}\right).$

Then $\mathbf{E}M_{T}(\Lambda)\leq 1$.

Combining the above exponential inequality with the ideas outlined in Section
4.1 along with a covering argument yields Theorem 4.1.

###  4.3 Notes

Theorem 4.1 holds for a fixed $T\in\mathbb{N}_{+}$, which is sufficient for
analyzing the system identification error. In contrast, the self-normalized
margtingale bound in Abbasi-Yadkori (2013) holds for an arbitrary stopping
time and thus uniformly for all $T\in\mathbb{N}_{+}$ by a stopping time
construction. This uniform bound may be required in some settings, e.g. in
error bounds for adapative control.

##  5 System Identification

In this section, we analyze well-known linear system identification algorithms
that rely on the least squares algorithm. Note that the problem formulation
changes with the system parameterization (e.g., state space, ARMAX, etc.).
However, a nice property of linear systems is that under certain conditions,
we can obtain a linear non-parametric ARX model by regressing the system
output to past outputs and inputs. Then, depending on the parameterization, we
can recover a particular realization. In the following, we first review ARX
identification, which can be seen as a fundamental building block for many
linear system identification algorithms. Then, we analyze identification of
Markov parameters in the case of state-space systems. We focus exclusively on
the case of single trajectory data.

###  5.1 ARX Systems

Consider an unknown vector autoregressive system with exogenous inputs (ARX)

$\displaystyle Y_{t}$
$\displaystyle=\sum_{i=1}^{p}A^{\star}_{i}Y_{t-i}+\sum_{j=1}^{q}B^{\star}_{j}U_{t-j}+\Sigma_{W}^{1/2}W_{t},$
(5.1)

where $Y_{t}\in\mathbb{R}^{d_{\mathsf{Y}}}$ are the system outputs,
$U_{t}\in\mathbb{R}^{d_{\mathsf{U}}}$ are the control (exogenous) inputs, and
$W_{t}\in\mathbb{R}^{d_{\mathsf{Y}}}$ is the normalized process noise with
$\Sigma_{W}\in\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{Y}}}$ capturing the
(non-normalized) noise covariance. Matrices $A^{\star}_{i},$ $i\leq p$ and
$B^{\star}_{j}$, $j\leq q$ contain the unknown ARX coefficients. For the
initial conditions, we assume $Y_{-1}=\cdots=Y_{-p}=0$,
$U_{-1}=\cdots=U_{-q}=0$.

######  Assumption 5.1 (System and Noise model).

Let the noise covariance $\Sigma_{W}\succ 0$ be full rank. Let the normalized
process noise $W_{t},t\geq 0$ be independent and identically distributed,
$K^{2}$-sub-Gaussian (see Definition 2.1), with zero mean and unit covariance
$\mathbf{E}W_{t}W^{\top}_{t}=I_{d_{\mathsf{Y}}}$. The orders $p,q$ are known.
System (5.1) is non-explosive, that is, the eigenvalues of matrix

$\mathcal{A}_{11}\triangleq\begin{bmatrix}A^{\star}_{1}&A^{\star}_{2}&\cdots&A^{\star}_{p-1}&A^{\star}_{p}\\\
I&0&\cdots&0&0\\\ 0&I&\cdots&0&0\\\ \vdots&&\ddots&&\vdots\\\
0&0&\cdots&I&0\end{bmatrix},\,$ (5.2)

lie strictly on or inside the unit circle $\rho(\mathcal{A}_{11})\leq 1$.

The techniques below can provide meaningful finite-sample bounds only when the
system is non-explosive. Deriving finite sample guarantees for identifying
explosively unstable partially-observed systems from single trajectory data is
open to the best of our knowledge (Tsiamis et al., 2023).

In this tutorial, we focus solely on white-noise excitation inputs. Excitation
strategies—experiment designs—beyond this simple structure form a vast
literature in system identification and statistics, see for instance Ljung
(1999); Gevers (2005); Pukelsheim (2006) and Bombois et al. (2011) and the
references therein. For a recent non-asymptotic treatment see also Wagenmaker
and Jamieson (2020).

######  Assumption 5.2 (White-noise excitation policy).

We assume that the control input is generated by a random i.i.d. Gaussian
process, that is, $U_{t}\sim\mathcal{N}(0,\sigma^{2}_{u}I)$.

Grouping all covariates into one vector and defining

$X_{t}=\begin{bmatrix}Y^{\top}_{t-1:t-p}&U^{\top}_{t-1:t-q}\end{bmatrix}^{\top},\,\theta^{\star}=\begin{bmatrix}A^{\star}_{1:p}&B^{\star}_{1:q}\end{bmatrix}$
(5.3)

we can re-write (5.1) in terms of (1.1)

$Y_{t}=\theta^{\star}X_{t}+\Sigma_{W}^{1/2}W_{t},$

where $W_{t}$ is independent from $X_{t}$, but $X_{t}$ has the special time-
dependent structure induced by (5.1). Given samples $(Y_{1:T},U_{0:T-1})$, the
least-squares estimate is given by

$\widehat{\theta}_{T}\triangleq\sum_{t=1}^{T}Y_{t}X^{\top}_{t}\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{\dagger},$
(5.4)

where we purposely highlight the dependence of the estimate on the number of
samples with the subscript $T$. Before we present the main result, let us
define some quantities which are related to the quality of system estimates.
The covariance at time $t\geq 0$ is defined as

$\Sigma_{t}\triangleq\mathbf{E}X_{t}X^{\top}_{t}.$ (5.5)

It captures the expected richness of the data, i.e., how excited the modes of
the system are on average. In particular, the relative excitation of the data
compared to the noise magnitude significantly affects the quality of system
identification. This motivates the definition of signal-to-noise (SNR) as the
ratio of the (directionally) worst-case excitation over the worst-case noise
magnitude

$\mathsf{SNR}_{t}\triangleq\frac{\lambda_{\min}(\Sigma_{t})}{\|\Sigma_{W}\|_{\mathsf{op}}K^{2}}.$
(5.6)

The following theorem provides a finite-sample upper bound on the performance
of the least-square estimator.

######  Theorem 5.1 (ARX Finite-Sample Bound).

Let $(Y_{1:T},U_{0:T-1})$ be single trajectory input-output samples generated
by system (5.1) under Assumptions 5.1, 5.2 for some horizon $T$. Fix a failure
probability $0<\delta<1$ and a time index $\tau\geq\max\\{p,q\\}$. Let
$T_{\mathsf{pe}}(\delta,\tau)\triangleq\min\\{t:t\geq
T_{0}(t,\delta/3,\tau)\\}$, where $T_{0}$ is defined in (5.9). If $T\geq
T_{\mathsf{pe}}(\delta,\tau)$, then with probability at least $1-\delta$

$\|\widehat{\theta}_{T}-\theta^{\star}\|_{\mathsf{op}}^{2}\leq\frac{C}{\mathsf{SNR}_{\tau}T}\left((pd_{\mathsf{Y}}+qd_{\mathsf{U}})\log\frac{pd_{\mathsf{Y}}+qd_{\mathsf{U}}}{\delta}+\log\det\left(\Sigma_{T}\Sigma^{-1}_{\tau}\right)\right),$
(5.7)

where $C$ is a universal constant, i.e., it is independent of system,
confidence $\delta$ and index $\tau$.

For non-explosive systems, matrix $\Sigma_{T}\Sigma^{-1}_{\tau}$ increases at
most polynomially with $T$ in norm (in view of Lemma 5.1). Hence, the
identification error decays with a rate of $\tilde{O}(1/\sqrt{T})$.

Dimensional dependence. Ignoring logarithmic terms, the bound implies that the
number of samples $T$ should scale linearly with the dimension
$pd_{\mathsf{Y}}+qd_{\mathsf{U}}$ of the covariates $X_{t}$. Since every
sample $Y_{t}$ contains at least $d_{\mathsf{Y}}$ measurements, this implies
that the total number of measurements should be linear with
$d_{\mathsf{Y}}\times(pd_{\mathsf{Y}}+qd_{\mathsf{U}})$. This scaling is
qualitatively correct since $\theta^{\star}$ has
$d_{\mathsf{Y}}\times(pd_{\mathsf{Y}}+qd_{\mathsf{U}})$ unknowns, requiring at
least as many independent equations.

Logarithmic dependence on confidence. The error norm scales linearly with
$\sqrt{\log 1/\delta}$. In the asymptotic regime we also recover the same
order of $\sqrt{\log 1/\delta}$ by applying the Central Limit Theorem (CLT).
However, in the regime of finite samples, obtaining the rate is non-trivial,
see Tsiamis et al. (2023), and requires the analysis presented in this
tutorial.

System theoretic constants. The identification error is directly affected by
the SNR of the system. The more the system is excited and the smaller the
noise, the better the SNR becomes. However, excitability varies heavily
depending on the system and the choice of excitation policy. In particular,
the system’s controllability structure can affect the degree of excitation
dramatically. Systems with poor controllability structure can exhibit SNR
which suffers from curse of dimensionality, i.e., the smallest eigenvalue of
$\Sigma_{\tau}$ degrades exponentially with the system dimension Tsiamis and
Pappas (2021).

The upper bound also increases with the logarithm of the “condition number"
$\det(\Sigma_{T}\Sigma^{-1}_{\tau})$. For stable systems, this condition
number is bounded since $\Sigma_{T}$ converges to a steady-state covariance as
$T$ increases; we can neglect it in this case. On the other hand, the term
might be significant in the case of general non-explosive systems. Let
$\kappa$ be the size of the largest Jordan block of $\mathcal{A}_{11}$ with
eigenvalues on the unit circle. Then, this term can be as large as $\kappa\log
T$.

Burn-in time. The upper bound holds as soon as the number of samples exceeds a
“burn-in" time $T_{\mathsf{pe}}(\delta,\tau)$. If the system is non-explosive,
$T_{\mathsf{pe}}(\delta,\tau)$ is always finite for fixed $\tau$. Exceeding
the burn-in time guarantees that we have persistency of excitation, that is,
all modes of the system are excited. The burn-in time increases as we require
more confidence $\delta$ and we chose larger time indices $\tau$. On the other
had, larger $\tau$ leads to larger $\Sigma_{\tau}$, which improves the
$\mathsf{SNR}_{\tau}$. In other words, there is a tradeoff between improving
the SNR and deteriorating the burn-in time. We analyze persistency of
excitation in more detail, in the next subsection.

Proof outline. We outline the proof below, the full proofs can be found in
Appendix E. To analyze the least squares error, observe that

$\widehat{\theta}_{T}-\theta^{\star}=\underbrace{\sum_{t=1}^{T}\Sigma_{W}^{1/2}W_{t}X^{\top}_{t}\left(\sum_{t=1}^{T}X_{t}X_{t}^{\top}\right)^{-1/2}}_{\text{noise}}\times\underbrace{\left(\sum_{t=1}^{T}X_{t}X_{t}^{\top}\right)^{-1/2}}_{\text{excitation}}$
(5.8)

where we assumed that the inverse exists. To deal with the second term, we
will prove persistency of excitation in finite time leveraging the techniques
of Section 3, which requires most of the work. To deal with the noise part we
will apply the self-normalized martingale methods, which are reviewed in
Section 4. We study both terms in the following subsections.

####  5.1.1 Persistency of Excitation in ARX Models

In this subsection, we leverage the result of Theorem 3.1 to prove persistency
of excitation. By persistency of excitation, we refer to the case when we have
rich input-output data, that is, data which characterize all possible
behaviors of the system. Recall the definition (1.9) of the empirical
covariance matrix

$\widehat{\Sigma}_{T}\triangleq\frac{1}{T}\sum_{t=1}^{T}X_{t}X_{t}^{\top}.$

Using this definition, the excitation term in the least squares error can be
re-written as $(T\widehat{\Sigma}_{T})^{-1/2}$. We say that persistency of
excitation holds if and only if the empirical covariance matrix is strictly
positive definite (full rank). In the following, we show that the full rank
condition holds with high probability, provided that the number of samples
exceeds a certain threshold, i.e., the burn-in time.

######  Theorem 5.2 (ARX PE).

Let $(Y_{1:T},U_{0:T-1})$ be input-output samples generated by system (5.1)
under Assumptions 5.1, 5.2 for some fixed horizon $T$. Fix a failure
probability $0<\delta<1$ and a time index $\tau\geq\max\\{p,q\\}$. Then,
$\lambda_{\min}(\Sigma_{\tau})>0$. Moreover, if $T$ is large enough

$T\geq T_{0}(T,\delta,\tau)\triangleq
1152\tau\max\\{K^{2},1\\}\Big{(}(pd_{\mathsf{Y}}+qd_{\mathsf{U}})\log
C_{\mathsf{sys}}(T,\tau)+\log(1/\delta)\Big{)}$ (5.9)

where

$C_{\mathsf{sys}}(T,\tau)\triangleq\frac{T}{3\tau}\frac{\|\Sigma_{T}\|_{\mathsf{op}}^{2}}{\lambda^{2}_{\min}(\Sigma_{\tau})},$

then,

$\mathbf{P}\left(\widehat{\Sigma}_{T}\succeq\frac{1}{16}\Sigma_{\tau}\right)\geq
1-\delta.$

The detailed proof can be found in Appendix E, we only sketch the main ideas
here. The first step is to express the covariates $X_{1:T}$ as a causal linear
combination of the noises and the inputs, mimicking (3.1). The evolution of
the covariates follows a state-space recursion

$X_{t+1}=\mathcal{A}X_{t}+\mathcal{B}V_{t},$ (5.10)

where we concatenate noises and inputs
$V_{t}\triangleq\begin{bmatrix}W^{\top}_{t}&U^{\top}_{t}\end{bmatrix}^{\top}$.
Matrices $\mathcal{A},\,\mathcal{B}$ are given by

$\displaystyle\mathcal{A}\triangleq\begin{bmatrix}\mathcal{A}_{11}&\mathcal{A}_{12}\\\
0&\mathcal{A}_{22}\end{bmatrix},\quad\mathcal{B}\triangleq\begin{bmatrix}\mathcal{B}_{1}&\mathcal{B}_{2}\end{bmatrix},\text{
where }$ $\displaystyle\mathcal{A}_{11}\text{ is defined
in~{}\eqref{eq:arx_A_matrix}},\,$
$\displaystyle\mathcal{A}_{12}\triangleq\begin{bmatrix}1\\\ 0\\\ \vdots\\\
0\end{bmatrix}\otimes\begin{bmatrix}B^{\star}_{1}&\cdots&B^{\star}_{q}\end{bmatrix}\in\mathbb{R}^{pd_{\mathsf{Y}}\times
qd_{\mathsf{U}}},\,$
$\displaystyle\mathcal{A}_{22}\triangleq\begin{bmatrix}0&\cdots&0&0&0\\\
I_{d_{\mathsf{U}}}&\cdots&0&0&0\\\ \vdots&\ddots&&\cdots&\\\
0&\cdots&I_{d_{\mathsf{U}}}&0&0\\\
0&\cdots&0&I_{d_{\mathsf{U}}}&0\end{bmatrix}\in\mathbb{R}^{qd_{\mathsf{U}}\times
qd_{\mathsf{U}}}\qquad\qquad\qquad$
$\displaystyle\mathcal{B}_{1}\triangleq\begin{bmatrix}\Sigma_{W}^{1/2}&0_{((p-1)d_{\mathsf{Y}}+qd_{\mathsf{U}})\times
d_{\mathsf{Y}}}\end{bmatrix}^{\top},\,$
$\displaystyle\mathcal{B}_{2}\triangleq\begin{bmatrix}0_{pd_{\mathsf{Y}}\times
d_{\mathsf{U}}}&I_{d_{\mathsf{U}}}&0_{(q-1)d_{\mathsf{U}}\times
d_{\mathsf{U}}}\end{bmatrix}^{\top}.$

The vector $X_{1:T}$ of all covariates satisfies the following causal linear
relation

$X_{1:T}=\underbrace{\begin{bmatrix}\mathcal{B}&0&\cdots&0\\\
\mathcal{A}\mathcal{B}&\mathcal{B}&\cdots&0\\\ \vdots&&\ddots&\\\
\mathcal{A}^{T-1}\mathcal{B}&\mathcal{A}^{T-2}\mathcal{B}&\cdots&\mathcal{B}\end{bmatrix}}_{\mathbf{L}}V_{0:T-1}.$
(5.11)

where the lower-block triangular matrix is the Toeplitz matrix generated by
the Markov parameters matrices $\mathcal{A},\,\mathcal{B}$. The second step is
to apply Theorem 3.1. The details can be found in the Appendix.

######  Remark 5.1 (Existence of burn-in time.).

For the above result to be meaningful, we need inequality (5.9) to be
feasible. For non-explosive systems, the system theoretic term $\log
C_{\mathsf{sys}}(T,\tau)$ increases at most logarithmically with $T$, since
$\Sigma_{t}$ increases polynomially with $t$ in view of Lemma 5.1. Hence, for
any fixed $\tau$, or, in general, any $\tau$ that increases mildly
(sublinearly) with $T$, e.g. $O(\sqrt{T})$, it is possible to satisfy (5.9).
Note that $\rho(\mathcal{A})=\rho(\mathcal{A}_{11})$ due to the triangular
structure of $\mathcal{A}$. Hence, by Assumption 5.1, system $\mathcal{A}$ is
also non-explosive.

######  Remark 5.2 (Unknown system orders $p,q$).

The result of Theorem 5.2 still holds if the orders $p,q$ are unknown and we
use the wrong orders $\hat{p},\hat{q}$ in the covariates $X_{t}$. We just need
to replace $p,q$ with $\hat{p},\hat{q}$ with $\hat{p},\hat{q}$ and revise the
size of $\Sigma$ accordingly in (5.9). The finite-sample bounds of Theorem 5.1
also hold (by revising accordingly), but only if we overestimate $p,q$, that
is $\hat{p}\geq p$, $\hat{q}\geq q$. This also generalizes the single
trajectory result of Du et al. (2022) to non-explosive systems.

The following supporting lemma proves that the $k$-th powers of non-explosive
matrices increase at most polynomially with $k$.

######  Lemma 5.1 (Lemma 1 in Tsiamis and Pappas (2021)).

Let $\mathcal{A}\in\mathbb{R}^{d\times d}$ have all eigenvalues inside or on
the unit circle, with $\|\mathcal{A}\|_{\mathsf{op}}\leq M$, for some $M>0$.
Then,

$\|\mathcal{A}^{k}\|_{\mathsf{op}}\leq(ek)^{d-1}\max\left\\{M^{d},1\right\\}.$
(5.12)

As a corollary, the covariance matrices $\Sigma_{t}$ also grow at most
polynomially with the time $t$.

####  5.1.2 Dealing with the Noise Term

In this subsection, we modify the noise term so that we can leverage Theorem
4.1, which cannot be applied directly. We first manipulate the inverse of
$T\widehat{\Sigma}_{T}$ to relate it to the inverse of
$\Sigma+T\widehat{\Sigma}_{T}$, for some carefully selected $\Sigma$. Inspired
by Sarkar and Rakhlin (2019), we leverage the result of Theorem 5.2. Under the
event that persistency of excitation holds we have
$\widehat{\Sigma}_{T}\succeq T\Sigma_{\tau}/16$. Thus, selecting
$\Sigma=T\Sigma_{\tau}/16$ guarantees that

$\left(T\widehat{\Sigma}_{T}\right)^{-1}\preceq
2\left(\Sigma+T\widehat{\Sigma}_{T}\right)^{-1}.$

We can now apply Theorem 4.1. To finish the proof we need to upper-bound the
determinant of $\log\det(\Sigma+T\widehat{\Sigma}_{T})$. It is sufficient to
establish a crude upper-bound on the empirical covariance
$T\widehat{\Sigma}_{T}$ as in the following lemma.

######  Lemma 5.2 (Matrix Markov’s inequality).

Fix a failure probability $\delta>0$. With probability at least $1-\delta$

$\widehat{\Sigma}_{T}\preceq\frac{pd_{\mathsf{Y}}+qd_{\mathsf{U}}}{\delta}\Sigma_{T}.$
(5.13)

A more refined upper bound can also be applied (see e.g. the proof of
Proposition 6.1 below or the results in Jedra and Proutiere (2022)).

###  5.2 State-Space Systems

In this subsection, we derive finite-sample guarantees for learning Markov
parameters of linear systems in state-space form. Consider the following
state-space system in the so-called innovation form:

$\displaystyle X_{t+1}$
$\displaystyle=A^{\star}X_{t}+B^{\star}U_{t}+F^{\star}\Sigma_{E}^{1/2}E_{t}$
(5.14) $\displaystyle Y_{t}$
$\displaystyle=C^{\star}X_{t}+\Sigma_{E}^{1/2}E_{t},$

where $A^{\star}\in\mathbb{R}^{d_{\mathsf{X}}\times d_{\mathsf{X}}}$,
$B^{\star}\in\mathbb{R}^{d_{\mathsf{X}}\times d_{\mathsf{U}}}$,
$F^{\star}\in\mathbb{R}^{d_{\mathsf{X}}\times d_{\mathsf{Y}}}$, and
$C^{\star}\in\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$ are _unknown_
state-space parameters. For the initial condition, we assume $X_{0}=0$. We
call the normalized noise process $E_{t}$ the innovation error process.
Similar to the ARX case, we focus on white-noise excitation inputs, namely
Assumption 5.2 also holds here. Moreover, we assume the following.

######  Assumption 5.3 (System and Noise model).

Let the noise covariance $\Sigma_{E}\succ 0$ be full rank. Let the normalized
innovation process $E_{t}$ be independent, identically distributed,
$K^{2}$-sub-Gaussian (see Definition 2.1), with zero mean and unit covariance
$\mathbf{E}E_{t}E^{\top}_{t}=I_{d_{\mathsf{Y}}}$. The order $d_{\mathsf{X}}$
is unknown. System (5.14) is non-explosive, that is, the eigenvalues of matrix
$A^{\star}$ lie strictly on or inside the unit circle $\rho(A)\leq 1$. The
system is also minimum-phase, i.e., the closed loop matrix

$A^{\star}_{\mathsf{cl}}\triangleq A^{\star}-F^{\star}C^{\star}$ (5.15)

has all eigenvalues inside the unit circle $\rho(A^{\star}_{\mathsf{cl}})<1$.

The innovation form (5.14) might seem puzzling at first. In particular, the
correlation between process and measurement noise via $F^{\star}$, and the
requirement $\rho(A^{\star}_{\mathsf{cl}})<1$ seem restrictive. However, the
representation (5.14) is standard in the system identification literature
Verhaegen and Verdult (2007). Moreover, as we review below, standard state-
space models have input-output second-order statistics, which are equivalent
to the ones generated by system (5.14) (for appropriate $F^{\star}$,
$\Sigma_{E}$).

######  Remark 5.3 (Generality of model).

System class (5.14) captures general state-space systems driven by Gaussian
noise. Consider the following state-space model

$\displaystyle S_{t+1}$ $\displaystyle=A^{\star}S_{t}+B^{\star}U_{t}+W_{t}$
(5.16) $\displaystyle Y_{t}$ $\displaystyle=C^{\star}S_{t}+V_{t},$

where $W_{t},V_{t}$ are i.i.d., independent of each other, mean-zero Gaussian,
with covariances $\Sigma_{W}$ and $\Sigma_{V}$ respectively. Assume that
$\Sigma_{V}\succ 0$ is full rank, the pair $(C^{\star},A^{\star})$ is
detectable, and the pair $(A^{\star},\Sigma_{W})$ is stabilizable. These three
assumptions imply that the Kalman filter of system (5.16) is well-defined
(Anderson and Moore, 2012). In particular, define the Riccati operator as

$\mathsf{RIC}(P)\triangleq
A^{\star}P(A^{\star})^{\top}+\Sigma_{W}-A^{\star}P(C^{\star})^{\top}(C^{\star}P(C^{\star})^{\top}+\Sigma_{V})^{-1}C^{\star}P(A^{\star})^{\top}$
(5.17)

and let $P^{\star}$ be the unique positive semidefinite solution of
$P^{\star}=\mathsf{RIC}(P^{\star})$. Then the Kalman filter gain is equal to

$F^{\star}=-A^{\star}P(C^{\star})^{\top}(C^{\star}P(C^{\star})^{\top}+\Sigma_{V})^{-1}.$
(5.18)

Assume that the initial state is also mean-zero Gaussian with covariance
$P^{\star}$ and independent of the noises. Finally set

$\Sigma_{E}=C^{\star}P^{\star}(C^{\star})^{\top}+\Sigma_{V}.$ (5.19)

Under the above assumptions and selection of $F^{\star}$, $\Sigma_{E}$ systems
(5.14) and (5.16) are statistically equivalent from an input-output
perspective, see Qin (2006). Both system descriptions lead to input-output
trajectories with identical statistics. Moreover, due to the properties of
Kalman filter, stability of $A^{\star}_{\mathsf{cl}}$ (minimum phase property)
and independence of $E_{t}$ are satisfied automatically (Anderson and Moore,
2012).

In this tutorial we will only focus on recovering the first few (logarithmic
in $T$-many) Markov parameters
$C^{\star}(A^{\star}_{\mathsf{cl}})^{i}B^{\star}$, $i\geq 0$ and
$C^{\star}(A^{\star}_{\mathsf{cl}})^{j}F^{\star}$, $j\geq 0$ of system (5.14).
From a learning theory point of view, this is also known as improper learning,
since the search space (finitely many Markov parameters) does not exactly, but
only approximately, coincide with the hypothesis class (state space models).
In principle, this forms the backbone of the SSARX method introduced by
Jansson (2003). One can then proceed to recover the original state-space
parameters (up to similarity transformation) from the Markov parameters by
employing some realization method. We refer to Oymak and Ozay (2021); Tsiamis
et al. (2023) for a discussion on this approach from a finite sample
perspective.

####  5.2.1 Reduction to ARX learning with Bias

Let $p>0$ be a past horizon. Denote the Markov parameters up to time $p$ by

$\theta^{\star}_{p}\triangleq[\begin{array}[]{cccccc}C^{\star}B^{\star}&\cdots&C^{\star}(A^{\star}_{\mathsf{cl}})^{p-1}B^{\star}&C^{\star}F^{\star}&\cdots&C^{\star}(A^{\star}_{\mathsf{cl}})^{p-1}F^{\star}\end{array}].$
(5.20)

Note that the innovation errors are equal to
$\Sigma_{E}^{1/2}E_{t}=Y_{t}-C^{\star}X_{t}$. Replacing this expression into
the state equation (5.14), we obtain

$X_{t}=A^{\star}_{\mathsf{cl}}X_{t-1}+B^{\star}U_{t-1}+F^{\star}Y_{t-1}.$

Unrolling the state equation $p$ times, we get

$Y_{t}=\underbrace{\theta^{\star}_{p}Z_{t}+\Sigma_{E}^{1/2}E_{t}}_{\text{ARX}}+\underbrace{C^{\star}(A^{\star}_{\mathsf{cl}})^{p}X_{t-p}}_{\text{bias}},$
(5.21)

where $Z_{t}$ includes the past $p$ covariates

$Z_{t}=\begin{bmatrix}Y^{\top}_{t-1:t-p}&U^{\top}_{t-1:t-p}\end{bmatrix}^{\top}.$
(5.22)

The above recursion is an approximate ARX equation. There is an additive bias
error term on top of the statistical noise. The least-squares solution is
given by

$\widehat{\theta}_{p,T}\triangleq\sum_{t=1}^{T}Y_{t}Z^{\top}_{t}\left(\sum_{t=1}^{T}Z_{t}Z_{t}^{\mathsf{T}}\right)^{\dagger},$
(5.23)

where we also highlight the dependence on the past $p$. By the minimum phase
assumption, the bias term decays exponentially with the past horizon $p$. This
follows from the fact that $A^{\star}_{\mathsf{cl}}$ is asymptotically stable,
while $X_{t}$ scales at most polynomially with $t$ (in view of Lemma 5.1). By
selecting $p=\Omega(\log T)$, we can make the bias term decay very fast,
making its contribution to the error
$\theta^{\star}_{p}-\widehat{\theta}_{p,T}$ negligible. On the other hand,
increasing the past horizon $p$ increases the statistical error since the
search space is larger.

####  5.2.2 Non-Asymptotic Guarantees

To derive finite-time guarantees for state space systems of the form (5.14),
we follow the same steps as in the case of ARX systems. However, we have to
account for the bias term and the fact that $p$ grows with $\log T$. Let us
define again the covariance at time $t\geq 0$

$\Sigma_{p,t}\triangleq\mathbf{E}Z_{t}Z^{\top}_{t},$ (5.24)

where we highlight the dependence on both the past horizon $p$ and the time
$t$. The covariance of the state is defined similarly

$\Sigma_{X,t}\triangleq\mathbf{E}X_{t}X^{\top}_{t}.$ (5.25)

Define the SNR as

$\mathsf{SNR}_{p,t}\triangleq\frac{\lambda_{\min}(\Sigma_{p,t})}{\|\Sigma_{E}\|_{\mathsf{op}}K^{2}}.$
(5.26)

Unlike the ARX case, here the SNR might degrade since we allow $p$ to grow
with $\log T$. For this reason, we require the following additional
assumption.

######  Assumption 5.4 (Non-degenerate SNR).

We assume that the SNR is uniformly lower bounded for all possible past
horizons

$\liminf_{t\geq 0}\mathsf{SNR}_{t,t}>0.$

Later on, in Theorem 5.4, we show that the above condition is non-vacuous and
is satisfied for quite general systems.

######  Theorem 5.3 (State Space Finite-Sample Bound).

Let $(Y_{1:T},U_{0:T-1})$ be single trajectory input-output samples generated
by system (5.14) under Assumptions 5.2, 5.3, 5.4, for some horizon $T$. Fix a
failure probability $0<\delta<1$ and select $p=\beta\log T$, for $\beta$ large
enough such that

$\|C^{\star}(A^{\star}_{\mathsf{cl}})^{p}\|_{\mathsf{op}}\|\Sigma_{X,T}\|_{\mathsf{op}}\leq
T^{-3}.$ (5.27)

Let $T^{\mathsf{ss}}_{\mathsf{pe}}(\delta,\beta)\triangleq\min\\{t:t\geq
T_{0}(t,\delta,\beta\log t)\\}$, where $T_{0}$ is defined in (5.9). If $T\geq
T^{\mathsf{ss}}_{\mathsf{pe}}(\delta,\beta)$, then with probability at least
$1-2\delta$

$\|\widehat{\theta}_{p,T}-\theta^{\star}_{p}\|_{\mathsf{op}}^{2}\leq\frac{C_{1}}{\mathsf{SNR}_{p,p}T}\left(p(d_{\mathsf{Y}}+d_{\mathsf{U}})\log\frac{p(d_{\mathsf{Y}}+d_{\mathsf{U}})}{\delta}+\log\det\left(\Sigma_{p,T}\Sigma^{-1}_{p,p}\right)\right),$
(5.28)

where $C_{1}$ is a universal constant, i.e., it is independent of system,
confidence $\delta$ and past horizon $p$.

For non-explosive systems, matrix $\Sigma_{p,T}\Sigma^{-1}_{p,p}$ increases at
most polynomially with $T$ in norm. Since the SNR is uniformly lower bounded,
the identification error decays with a rate of $\tilde{O}(1/\sqrt{T})$. The
bound seems similar to the one for ARX systems for $\tau=p$. However, since
$p\asymp\log T$, we have an extra logarithmic term.

Role of $\beta$. Recall the approximate ARX relation (5.21). For the bias term
to be small, the exponentially decaying $(A^{\star}_{\mathsf{cl}})^{p}$ should
counteract the magnitude of the state $\|X_{t-p}\|_{\mathsf{op}}$.
Intuitively, the state grows as fast as
$\|\Sigma_{X,t}\|_{\mathsf{op}}^{1/2}$, where
$\Sigma_{X,t}=\mathbf{E}X_{t}X^{\top}_{t}$. Hence the state norm grows at most
polynomially with $T$. Meanwhile,
$\|(A^{\star}_{\mathsf{cl}})^{p}\|_{\mathsf{op}}=O(\rho^{p})$ for some
$\rho>\rho(A^{\star}_{\mathsf{cl}})$. With the choice $p=\beta\log T$, we get
$\|(A^{\star}_{\mathsf{cl}})^{p}\|_{\mathsf{op}}=O(T^{-\beta/\log(1/\rho)})$.
Hence, if we select large enough $\beta$, we can make the bias term very
small, even smaller than the dominant $\tilde{O}(1/\sqrt{T})$ term.

Burn-in time. Since the system is non-explosive,
$T^{\mathsf{ss}}_{\mathsf{pe}}(\delta,\beta)$ is always finite under
Assumption 5.4, for any $\beta$. As before, exceeding the burn-in time
guarantees that we have persistency of excitation. Naturally, larger $\beta$
lead to larger past horizons $p$, which, in turn, increase the burn-in time.

Finally, we prove that Assumption 5.4 is non-vacuous. It is sufficient for
$F^{\star}$ and $\Sigma_{W}$ to be generated by a Kalman filter as in (5.18),
(5.19).

######  Theorem 5.4.

Consider system (5.14) and the definition of $\mathsf{SNR}_{p,t}$ in (5.26).
If matrices $F^{\star}$, $\Sigma_{E}$ are generated as in (5.18), (5.19) with
$(A^{\star},\Sigma_{W}^{1/2})$ stabilizable, $(C^{\star},A^{\star})$
detectable and $\Sigma_{V}\succ 0$, then the SNR is uniformly lower bounded
$\liminf_{t\geq 0}\mathsf{SNR}_{t,t}>0$.

Both conditions are sufficient. It is subject of future work to extend the
result to more general non-exposive systems.

###  5.3 Notes

The exposition above is inspired by prior work on identifying fully-observed
systems (Faradonbeh et al., 2018; Simchowitz et al., 2018; Sarkar and Rakhlin,
2019) and partially-observed systems (Oymak and Ozay, 2019; Simchowitz et al.,
2019; Sarkar et al., 2021; Tsiamis and Pappas, 2019; Lee and Lamperski, 2020;
Lale et al., 2021; Lee, 2022). For a wider overview of the literature, we
refer the reader to Tsiamis et al. (2023).

Let us further remark that the guarantee for the ARX model in Theorem 5.1 is
almost optimal. The use of Matrix Markov’s inequality yields extraneous
dependency on the problem dimension multiplying the deviation term
$\log(1/\delta)$. This can in principle be removed by a more refined analysis
(see e.g. the proof of Proposition 6.1 below or the results in Jedra and
Proutiere (2022)). Indeed, the signal-to-noise term (5.6) is closely related
to the Fisher Information Matrix appearing in the classical asymptotic
optimality theory. Let us also point out that the question of optimality in
identifying partially observed state-space systems is more subtle, and while
consistent, the bounds presented here are not (asymptotically) optimal.

##  6 An Alternative Viewpoint: the Basic Inequality

In many situations, the choice of the model class
$\mathsf{M}=\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$ leading to (1.8)
is not appropriate. For instance physical or other modelling considerations
might have already informed us that the true $\theta^{\star}$ belongs to some
smaller model class such as the family of low rank or sparse matrices which
are strict subsets of $\mathsf{M}$. Other properties one might wish to enforce
include, stable, low norm, or even passivity-type properties. In either of the
above examples no error expression of the form (1.10) is directly available.
Instead, we observe by optimality of $\widehat{M}$ to the optimization program
(1.7) that

$\frac{1}{T}\sum_{t=1}^{T}\|Y_{t}-\widehat{\theta}X_{t}\|_{2}^{2}\leq\frac{1}{T}\sum_{t=1}^{T}\|Y_{t}-\theta^{\star}X_{t}\|_{2}^{2}.$
(6.1)

Expanding the squares and re-arranging terms we arrive at the so-called basic
inequality of least squares:

$\frac{1}{T}\sum_{t=1}^{T}\|(\widehat{\theta}-\theta^{\star})X_{t}\|_{2}^{2}\leq\frac{2}{T}\sum_{t=1}^{T}\langle
V_{t},(\widehat{\theta}-\theta^{\star})X_{t}\rangle.$ (6.2)

The inequality (6.2) serves as an alternative to the explicit error equation
(1.10). To drive home this point, let us first re-arrange (6.2) slightly:

$\frac{1}{T}\sum_{t=1}^{T}\|(\widehat{\theta}-\theta^{\star})X_{t}\|_{2}^{2}\leq\frac{4}{T}\sum_{t=1}^{T}\langle
V_{t},(\widehat{\theta}-\theta^{\star})X_{t}\rangle-\frac{1}{T}\sum_{t=1}^{T}\|(\widehat{\theta}-\theta^{\star})X_{t}\|_{2}^{2}.$
(6.3)

Note now that $\widehat{\theta}-\theta^{\star}$ are elements of
$\mathsf{M}_{\star}\triangleq\mathsf{M}-\theta^{\star}$. Hence—by considering
the worst-case (supremum) right hand side of (6.3)—we obtain:

$\frac{1}{T}\sum_{t=1}^{T}\|(\widehat{\theta}-\theta^{\star})X_{t}\|_{2}^{2}\leq\sup_{\theta\in\mathsf{M}_{\star}}\left\\{\frac{4}{T}\sum_{t=1}^{T}\langle
V_{t},\theta X_{t}\rangle-\frac{1}{T}\sum_{t=1}^{T}\|\theta
X_{t}\|_{2}^{2}\right\\}.$ (6.4)

In fact, if $\mathsf{M}=\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$, the
optimization on the right hand side of (6.4) has an explicit solution. This
implies that we always have the following upper-bound on the event that the
design is nondegenerate:

$\displaystyle\sup_{\theta\in\mathsf{M}_{\star}}\left\\{\frac{4}{T}\sum_{t=1}^{T}\langle
V_{t},\theta X_{t}\rangle-\frac{1}{T}\sum_{t=1}^{T}\|\theta
X_{t}\|_{2}^{2}\right\\}$ (6.5)
$\displaystyle\leq\sup_{\theta\in\mathbb{R}^{d_{\mathsf{Y}}\times
d_{\mathsf{X}}}}\left\\{\frac{4}{T}\sum_{t=1}^{T}\langle V_{t},\theta
X_{t}\rangle-\frac{1}{T}\sum_{t=1}^{T}\|\theta X_{t}\|_{2}^{2}\right\\}$
$\displaystyle(\mathsf{M}_{\star}\subset\mathbb{R}^{d_{\mathsf{Y}}\times
d_{\mathsf{X}}})$
$\displaystyle=\frac{4}{T}\left\|\left(\sum_{t=1}^{T}V_{t}X_{t}^{\mathsf{T}}\right)\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{-1/2}\right\|^{2}_{F}.$
$\displaystyle(\textnormal{direct calculation})$

Hence, we have in principle recovered an in-norm version of (1.10) with
slightly worse constants. Put differently, we may regard (6.4) as a
variational (or dual) form of the explicit error (1.10). Now, the advantage of
(6.4) is twofold:

  1. 1.

(6.4) and (6.5) hold for any
$\mathsf{M}_{\star}\subset\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$
and hence allows us to analyze the LSE (1.7) beyond OLS
($\mathsf{M}_{\star}=\mathbb{R}^{d_{\mathsf{Y}}\times d_{\mathsf{X}}}$). This
is important in identification problems where the parameter space is
restricted.

  2. 2.

We do not have to rely on (6.5) to control (6.4). In fact, for many reasonable
classes of $\mathsf{M}_{\star}\subset\mathbb{R}^{d_{\mathsf{Y}}\times
d_{\mathsf{X}}}$ we are able to give alternative arguments that are much
sharper (in terms of e.g. dimensional scaling) than the naive bound (6.5). See
Section 6.1 below.

A third advantage of the variational form (6.4) is that it generalizes
straightforwardly beyond linear least squares. In fact, none of the steps
(6.1),(6.2), (6.3) and (6.4) relied on the linearity of
$x\mapsto\widehat{\theta}x$ or that of $x\mapsto\theta^{\star}x$
($x\in\mathbb{R}^{d_{\mathsf{X}}}$). We will explore this theme further in
Section 6.1 and Section 7.

###  6.1 Sparse Autoregressions

Before we proceed to sketch out how the basic inequality above extends to
nonlinear problems in Section 7, let us use it to analyze a simple variation
of the autoregression already encountered in Section 5. Namely, the
autoregressive model (5.1) which—for simplicity—is further assumed one-
dimensional:

$Y_{t}=\sum_{i=1}^{p}A^{\star}_{i}Y_{t-i}+W_{t}\qquad(W_{1:T}\textnormal{
iid\, mean zero and $\sigma^{2}$-subG})$ (6.6)

and assume in addition that it is known that only $s\in\mathbb{N}$ of the $p$
entries of $\theta^{\star}=[A^{\star}_{1},\dots,A^{\star}_{p}]$ are nonzero.
Put differently, the vector $\theta^{\star}$ is known to be $s$-sparse and we
write $\theta^{\star}\in\\{\theta\in\mathbb{R}^{p}:\|\theta\|_{0}\leq
s\\}\triangleq\mathsf{M}$. Hence, in this case the model class $\mathsf{M}$ is
the union of ${p}\choose s$ subspaces. Clearly, we could use OLS (1.8) but
this estimator does not take advantage of the additional information that
$A^{\star}=\theta^{\star}$ lies in the $s$-dimensional submanifold
$\mathsf{M}$. Intuitively, if $s\ll p$ this set should be much smaller than
$\mathbb{R}^{p}$ and so one expects that identification occurs at a faster
rate.

In this section we demonstrate that the least squares estimator (1.7) in which
the search is restricted to the low-dimensional manifold $\mathsf{M}$
outperforms the OLS. We stress that this is _not_ a computationally efficient
estimator and the results in this section should be thought of as little more
than an illustrative example.

Returning to the problem of controlling the error of this estimator, we note
that in this case there is no closed form for the LSE and we do not have
direct access to the error equation (1.10).555Although, in this particular
case an alternative analysis based on this equation is possible. Hence, we
instead use the offset basic inequality approach from Section 6. As before, it
is convenient to set $X_{t}=[Y_{t-1},\dots,Y_{t-p}]^{\mathsf{T}}$. With this
additional bit of notation in place, we recall from (6.4) that:

$\frac{1}{T}\sum_{t=1}^{T}\|(\widehat{\theta}-\theta^{\star})X_{t}\|_{2}^{2}\leq\max_{\theta\in\mathsf{M}_{\star}}\left\\{\frac{4}{T}\sum_{t=1}^{T}W_{t}\theta
X_{t}-\frac{1}{T}\sum_{t=1}^{T}|\theta X_{t}|_{2}^{2}\right\\}$ (6.7)

where $\mathsf{M}_{\star}$ is the translation $\mathsf{M}-\theta^{\star}$.
Since $\mathsf{M}$ is the union of ${p}\choose s$-many linear $s$-dimensional
subspaces $S\subset\mathbb{R}^{d_{\mathsf{X}}\times d_{\mathsf{X}}}$,
$\mathsf{M}_{\star}$ is the union of ${p}\choose s$ affine subspaces
$s$-dimensional affine subspaces of the form $S-\theta^{\star}$. Let us also
note that
$\mathsf{M}_{\star}\subset\mathsf{M}-\mathsf{M}=\\{\theta\in\mathbb{R}^{p}:\|\theta\|_{0}\leq
2s\\}$. Consequently:

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}\|(\widehat{\theta}-\theta^{\star})X_{t}\|_{2}^{2}$
$\displaystyle\leq\max_{\theta\in\mathsf{M}_{\star}}\left\\{\frac{4}{T}\sum_{t=1}^{T}W_{t}\theta
X_{t}-\frac{1}{T}\sum_{t=1}^{T}|\theta X_{t}|_{2}^{2}\right\\}$ (6.8)
$\displaystyle\leq\max_{S}\max_{\theta\in
S}\left\\{\frac{4}{T}\sum_{t=1}^{T}W_{t}\theta
X_{t}-\frac{1}{T}\sum_{t=1}^{T}|\theta X_{t}|_{2}^{2}\right\\}.$

where maximization over $S$ occurs over the $p\choose 2s$-many sparse
subspaces. Notice now that since $\theta$ in $\eqref{eq:sparseoffset}$ is
$s$-sparse, the products $\theta X_{t}$ are just $\theta X_{t}=\sum_{i\in
S}\theta_{i}(X_{t})_{i}$ where we have abused notation and identified $S$ with
its support set. Hence, by the same direct calculation as in (6.5), if we
denote $(X_{t})_{S}$ the $s$-dimensional vector obtained by coordinate
projection onto part of $S$ not constrained to be identically zero (i.e. the
image of the projection onto $S$ represented as the $s$-dimensional Euclidean
space) we find that:

$\frac{1}{T}\sum_{t=1}^{T}\|(\widehat{\theta}-\theta^{\star})X_{t}\|_{2}^{2}\leq\frac{4}{T}\max_{S}\left\|\left(\sum_{t=1}^{T}W_{t}(X_{t})_{S}\right)\left(\sum_{t=1}^{T}(X_{t})_{S}(X_{t})_{S}^{\mathsf{T}}\right)^{-1/2}\right\|^{2}_{2}.$
(6.9)

The right hand side of (6.9) can be controlled by the self-normalized
inequality in Theorem 4.1 for each fixed $S$. Moreover, there are only
$p\choose 2s$ such subspaces, so we can apply a union bound to control the
maximum over these subspaces. Note also that the left hand side of (6.9) can
be controlled by the tools developed in Section 3. Carrying out these steps
leads to the following guarantee.

######  Proposition 6.1.

Fix $T,k\in\mathbb{N}$ with $T$ divisible by $k$ and let $\mathbf{L}$ be the
linear operator defined in (5.11). Let $\widehat{\theta}$ be the LSE (1.7)
over the set $\mathsf{M}=\\{\theta\in\mathbb{R}^{p}:\|\theta\|_{0}\leq s\\}$
for the system (6.6). Define
$\Sigma_{j}\triangleq\frac{1}{j}\sum_{t=1}^{j}\mathbf{E}X_{t}X_{t}^{\mathsf{T}}$
for $j\in[T]$ and

$\mathrm{cond}_{\mathsf{sys}}(T,k)\triangleq\left(1+\frac{\|\mathbf{L}\mathbf{L}^{\mathsf{T}}\|_{\mathsf{op}}}{k\lambda_{\min}\left(\Sigma_{T}\right)}\right)\frac{\lambda_{\max}\left(\Sigma_{T}\right)}{\lambda_{\min}\left(\Sigma_{k}\right)}.$

There exist univeral positive constants $c,c^{\prime}$ such that for any
$\delta\in(0,1)$ it holds with probability at least $1-\delta$ that:

$\|(\widehat{\theta}-\theta^{\star})\sqrt{\Sigma_{k}}\|_{2}^{2}\leq
c\sigma^{2}\times\frac{s\log\left(\frac{p\times\mathrm{cond}_{\mathsf{sys}}(T,k)}{s}\right)+\log(1/\delta)}{T}$
(6.10)

as long as

$T/k\geq
c^{\prime}\sigma^{2}\left(s\left[\log\left(\mathrm{cond}_{\mathsf{sys}}(T,k)\right)+\log(p/s)\right]+\log(1/\delta)\right).$
(6.11)

A few remarks are in order. The guarantee (6.10) depends on the dimension $s$
of $\mathsf{M}$, and not the total parameter dimension $p$. Similarly, the
burn-in (6.11) exhibits a similar win, by depending linearly on $s$ and only
logarithmically on $p$. There is also the difference that the left hand side
of (6.10) is given in the problem-dependent Mahalanobis norm induced by
$\Sigma_{k}$ and opposed to just the standard Euclidean $2$-norm. This implies
that if we actually want parameter identification in the sense of the previous
section, a restricted eigenvalue condition on $\Sigma_{k}$ is needed.666Note
that $\widehat{\theta}-\theta^{\star}$ is $2s$-sparse. Indeed, for some
positive number $\lambda_{\mathrm{restricted}}$, one requires that
$v^{\mathsf{T}}\Sigma_{k}v\geq\lambda_{\mathrm{restricted}}$ for all
$2s$-sparse vectors $v$ on the unit sphere: $v\in\mathbb{S}^{p-1}$ and
$\|v\|_{0}\leq 2s$. Obviously the requirements on
$\lambda_{\mathrm{restricted}}$ are much milder than the corresponding ones on
$\lambda_{\min}(\Sigma_{k})$ and we always have
$\lambda_{\mathrm{restricted}}\geq\lambda_{\min}(\Sigma_{k})$.

The following lemma is central. Namely, we begin the proof of Proposition 6.1
by restricting to an event in which the designs
$\sum_{t=1}^{T}(X_{t})_{S}(X_{t})_{S}^{\mathsf{T}}$ are sufficiently well-
conditioned for all the subspaces $S$ at once. The requirements on this event
are relatively milder than the corresponding one over $\mathbb{R}^{p}$ and
explains the "dimensional win" (when $s\ll p$) of the sparse estimator over
OLS.

######  Lemma 6.1.

Let $\mathbf{L}$ be the linear operator defined in (5.11). Fix
$\delta\in(0,1)$ and let $T$ be divisible by $k\in\mathbb{N}$. There exist
universal positive constants $c_{1},c_{2},c_{3}\in\mathbb{R}$ such that the
following two-sided control holds uniformly in $S$ with probability
$1-\delta$:

$\frac{c_{1}}{k}\sum_{t=1}^{k}\mathbf{E}\left[(X_{t})_{S}(X_{t})_{S}^{\mathsf{T}}\right]\preceq\frac{1}{T}\sum_{t=1}^{T}(X_{t})_{S}(X_{t})_{S}^{\mathsf{T}}\\\
\preceq
c_{2}\left(1+\frac{T\|\mathbf{L}\mathbf{L}^{\mathsf{T}}\|_{\mathsf{op}}}{k\lambda_{\min}\left(\sum_{t=0}^{T-1}\mathbf{E}X_{t}X_{t}^{\mathsf{T}}\right)}\right)\left(\frac{1}{T}\sum_{t=1}^{T}\mathbf{E}\left[(X_{t})_{S}(X_{t})_{S}^{\mathsf{T}}\right]\right)$
(6.12)

as long as

$T\geq c_{3}\sigma^{2}\left(s\left[\log
C_{\mathsf{sys}}+\log(p/s)\right]+\log(1/\delta)\right).$ (6.13)

Equation 6.13 is revealing about the advantage of using the sparse estimator
searching over $\mathsf{M}=\\{\theta\in\mathbb{R}^{p}:\|\theta\|_{0}\leq
s\\}$. The burn-in period in (6.13) is proportional to the dimension of the
low-dimensional parameter manifold $\mathsf{M}$ instead of that of the latent
space $\mathbb{R}^{p}$. Finally, as usual we have relegated the full proof of
Proposition 6.1 to the appendix, see Section F.1.

###  6.2 Notes

The variational formulation of the least squares error—the basic inequality
(6.2)—is standard in the nonparametric statistics literature (see e.g.
Wainwright, 2019, Chapters 13 and 14). The idea to rewrite the basic
inequality (6.2) as (6.3) was introduced to the statistical literature by
Liang et al. (2015), but has its roots in online learning (Rakhlin and
Sridharan, 2014).

##  7 Beyond Linear Models

Let us now make another gradual shift of perspective. Instead of considering
the linear model (1.1) introduced in Section 1.1 we consider the following
_nonlinear_ regression model:

$Y_{t}=f^{\star}(X_{t})+V_{t},\qquad t\in[T].$ (7.1)

As before, $Y_{1:T}$,$X_{1:T}$ and $V_{1:T}$ are stochastic processes taking
values in $\mathbb{R}^{d_{\mathsf{Y}}}$ and $\mathbb{R}^{d_{\mathsf{X}}}$
respectively. However, this time $f^{\star}$ is no longer constrained to be a
linear map of the form $x\mapsto Ax$ for matrix $A$. Rather, we suppose that
$f^{\star}$ in (7.1) belongs to some (square integrable) space of functions
$\mathscr{F}$ such that $\mathscr{F}\ni f:x\mapsto f(x)$. It is perhaps now
that the motivation behind the change of perspective from Section 6 becomes
most apparent: the basic inequality (6.3) remains valid. To be precise, let us
define the _nonlinear_ least squares estimator

$\widehat{f}\in\operatorname*{argmin}_{f\in\mathscr{F}}\left\\{\frac{1}{T}\sum_{t=1}^{T}\|Y_{t}-f(X_{t})\|_{2}^{2}\right\\}.$
(7.2)

Let $\mathscr{F}_{\star}\triangleq\mathscr{F}-f^{\star}$. By the exact same
optimality argument as in Section 6, the reader can now readily verify that:

$\frac{1}{T}\sum_{t=1}^{T}\|\widehat{f}(X_{t})-f^{\star}(X_{t})\|^{2}_{2}\leq\sup_{f\in\mathscr{F}_{\star}}\frac{1}{T}\left(\sum_{t=1}^{T}4\langle
V_{t},f(X_{t})\rangle-\sum_{t=1}^{T}\|f(X_{t})\|^{2}_{2}\right).$ (7.3)

What does (7.3) entail in terms of estimating the unknown function
$f^{\star}$? To answer this, we first need to define a performance criterion.
The simplest one is small average $L^{2}$-norm-error, where

$f\in\mathscr{F}:\quad\|f\|_{L^{2}}^{2}\triangleq\frac{1}{T}\sum_{t=1}^{T}\mathbf{E}\|f(X_{t})\|^{2}_{2}.$
(7.4)

The program we have carried out in the previous sections now generalizes as
follows:

  * •

First, prove a so-called lower uniform law. That is to say, we wish to show
that with overwhelming probability

$\|f-f_{\star}\|_{L^{2}}^{2}\leq\frac{C}{T}\sum_{t=1}^{T}\|f(X_{t})-f_{\star}(X_{t})\|^{2}\quad(\textnormal{simultaneously
}\forall f\in\mathscr{F}).$ (7.5)

for some universal positive constant $C$.

  * •

Second, control the supremum of the _empirical process_ :

$f\mapsto\left(\sum_{t=1}^{T}4\langle
V_{t},f(X_{t})\rangle-\sum_{t=1}^{T}\|f(X_{t})\|^{2}_{2}\right)$ (7.6)

in terms of the noise level $\sigma$ and the complexity of the class
$\mathscr{F}$.

By combining (7.5) and (7.6) we arrive at a high probability bound of the
form:

$\|\widehat{f}-f^{\star}\|^{2}_{L^{2}}\leq\frac{C}{T}\sum_{t=1}^{T}\|f(X_{t})-f_{\star}(X_{t})\|^{2}\leq\frac{C\times\mathrm{comp}(\mathscr{F},\sigma^{2})+\textnormal{deviation
term}}{T}.$ (7.7)

A statement of this form is given as Theorem 7.1 below.

######  Remark 7.1.

It is worth to take pause and appreciate the analogy to the analysis of linear
regression models. The first step (7.5) exactly corresponds to controlling the
lower spectrum of the empirical covariance matrix. Suppose for simplicity that
$d_{\mathsf{Y}}=1$. Then for a linear map $\mathbb{S}^{d_{\mathsf{X}}-1}\ni
f\mapsto\langle f,x\rangle$ we have:

$\frac{1}{T}\sum_{t=1}^{T}\|f(X_{t})\|^{2}_{2}=\frac{1}{T}\sum_{t=1}^{T}\langle
f,(X_{t}X_{t}^{\mathsf{T}})f\rangle=\left\langle
f,\left[\frac{1}{T}\sum_{t=1}^{T}(X_{t}X_{t}^{\mathsf{T}})\right]f\right\rangle$
(7.8)

which are just the one-dimensional projections of the empirical covariance
matrix (1.9). In the context of linear models, establishing (7.5) was the
topic of Section 3. Analogously, for a linear predictor, the $L^{2}$-norm
(7.4) becomes a Mahalanobis norm:
$f\in\mathbb{R}^{d_{\mathsf{X}}}\Rightarrow\|f\|_{L^{2}}^{2}=\langle
f,\Sigma_{X}f\rangle$ for
$\Sigma_{X}=\frac{1}{T}\sum_{t=1}^{T}\mathbf{E}X_{t}X_{t}^{\mathsf{T}}$.

Moreover, For linear models, we had:

$\sup_{f\in\mathbb{R}^{d_{\mathsf{X}}}}\left(\sum_{t=1}^{T}4\langle
V_{t},f(X_{t})\rangle-\sum_{t=1}^{T}\|f(X_{t})\|^{2}_{2}\right)=4\left\|\left(\sum_{t=1}^{T}V_{t}X_{t}^{\mathsf{T}}\right)\left(\sum_{t=1}^{T}X_{t}X_{t}^{\mathsf{T}}\right)^{-1/2}\right\|^{2}_{F}.$
(7.9)

Analyzing terms of this form was the topic of Section 4.

In other words, the approach outlined above is very much in the same spirit as
that in the rest of the manuscript. There are a few changes that need to be
made since we less access to linearity in our argument, but in principle the
key difference is that we will have to replace the indexing set
$\mathbb{S}^{d-1}$ with a more general function class $\mathscr{F}_{\star}$.

###  7.1 Many Trajectories and Finite Hypothesis Classes

In order to make the exposition self-contained, we will now make two
simplifying assumptions relating to the finiteness of the hypothesis class and
the dependence structure of the covariate process $X_{1:T}$. A more general
treatment without these can be found in Ziemann and Tu (2022). Here, we impose
the following:

  1. A1.

The hypothesis class $\mathscr{F}$ is finite; $|\mathscr{F}|<\infty$.

  2. A2.

We have access to $T/k$-many independent trajectories from the same process:
there exists an integer $k\in\mathbb{N}$ dividing $T$ such that
$X_{1:k},X_{k+1:2k},\dots$ are drawn iid.

We will also impose the following rather minimal integrability condition:

  1. A3.

All functions $f\in\mathscr{F}$ are such that
$\mathbf{E}\|f(X_{t})\|_{2}^{4}<\infty$ for all $t\in[T]$.

Moreover, as in Section 5, we require the noise to be a sub-Gaussian
martingale difference sequence:

  1. A4.

For each $t\in[T]$, $V_{t}|X_{1:t}$ is $\sigma^{2}$ conditionally-sub-Gaussian
and mean zero.

######  Remark 7.2.

Note that A4. above entails that $f_{\star}(x)=\mathbf{E}[Y_{t}|X_{t}=x]$ for
every time instance $t$ and so the setup is akin to the study of "predictor
models" from system identification (see e.g. Davis and Vinter, 1985, Chapter
2.6).

Under these assumptions, the main result of Ziemann and Tu (2022) essentially
simplifies to the following theorem.

######  Theorem 7.1.

Impose A1-A4, fix $\delta\in(0,1)$ and define

$\mathrm{cond}_{\mathscr{F}}\triangleq\max_{f\in\mathscr{F}_{\star}}\max_{t\in
T}\frac{\sqrt{\mathbf{E}\|f(X_{t})\|^{4}_{2}}}{\mathbf{E}\|f(X_{t})\|^{2}_{2}}.$
(7.10)

Suppose further that

$T/k\geq
4\mathrm{cond}_{\mathscr{F}}^{2}\left(\log|\mathscr{F}|+\log(2/\delta)\right)$

then we have that:

$\|\widehat{f}-f^{\star}\|_{L^{2}}^{2}\leq
16\sigma^{2}\left(\frac{\log(|\mathscr{F}|)+\log(2/\delta)}{T}\right).$ (7.11)

A few remarks are in order. The structure of Theorem 7.1 is by now familiar
and it is very much of the same structure as our previous results, cf. (1.5).
The key differences are that: (1) we now control the $L^{2}$ norm of our
estimator instead of the Euclidean or spectral norm; and (2) that the
dimensional dependency has been replaced by the complexity term
$\log|\mathscr{F}_{\star}|$. The proof is also structurally similar, as noted
in Remark 7.1. We also caution the reader that (7.11) is strictly a
statistical guarntee; we have said nothing—and will say nothing more—about the
computational feasibility of the estimator (7.2).

Let us now discuss A1-A4. Assumption A1 informs us that the search space for
the LSE (7.2) is finite. This is mainly imposed to avoid the introduction of
the chaining technique which is the standard alternative to the bounds from
Section 4. Using this technique, similar statements can for instance be
derived for compact subsets of bounded function classes (Ziemann and Tu,
2022). Assumption A2 controls the dependence structure of the process. Here,
we assume that we are able to restart the process every $k$ time steps. Again,
a more general statement relying on stochastic stability can be found in
Ziemann and Tu (2022). Assumption A3 is relatively standard. Arguably the
strongest assumption is A4, which in principle yields that the conditional
expectation (given all past data) is a function in the search space
$\mathscr{F}$. It is a so-called realizability assumption—the model (7.1) is
well-specified—and it is not currently known how to remove it and still obtain
sharp bounds beyond linear classes (for an analysis of linear misspecified
models, see Ziemann et al., 2023).

###  7.2 Notes

As noted in the previous section, the idea of using the “offset” basic
inequality relied on here is due to Rakhlin and Sridharan (2014); Liang et al.
(2015). The “many trajectores”-style of analysis used here is due to Tu et al.
(2022) who introduced it in the linear setting. Here, we have extended their
style of analysis to simplify the exposition of Ziemann and Tu (2022) who
consider the single trajectory setting, but rely on a rather more advanced
exponential inequality due to Samson (2000). Alternatively, one can also
easily extend the lower uniform law in Proposition F.1 to certain classes of
mixing processes by invoking the blocking technique of Yu (1994) combined with
a truncation style of argument such as that used in the proof of Theorem 14.12
of Wainwright (2019), see also Ziemann et al. (2023, the proof of Theorem
4.3). Note however, that all the analyses above and in this section
necessitate some degree of stability (mixing). This should be contrasted with
the system identification bounds of Section 5, which work even in the
marginally regime. In principle, the consequence of this is that while the
convergence rates for bounds such as Theorem 7.1 are correct, the burn-ins are
deflated by various dependency measures.

There have also been other, more algorithmically focused, approaches to
nonlinear identification problems in the recent literature. Noteably, gradient
based methods in generalized linear models of the form
$X_{t+1}=\phi(A^{\star}X_{t})+V_{t}$ (with $\phi$ a known nonlinearity) have
been the topic of a number of recent papers (see e.g. Foster et al., 2020;
Sattar and Oymak, 2022). The sharpest bounds for parameter recovery in this
setting are due to Kowshik et al. (2021).

##### Acknowledgements.

Ingvar Ziemann is supported by a Swedish Research Council international
postdoc grant. Nikolai Matni is funded by NSF awards CPS-2038873, CAREER award
ECCS-2045834, and ECCS-2231349.

## References

  * Abbasi-Yadkori [2013] Yasin Abbasi-Yadkori.  Online learning for linearly parametrized control problems.  2013\.
  * Anderson and Moore [2012] Brian DO Anderson and John B Moore.  _Optimal Filtering_.  Courier Corporation, 2012.
  * Bombois et al. [2011] Xavier Bombois, Michel Gevers, Roland Hildebrand, and Gabriel Solari.  Optimal experiment design for open and closed-loop system identification.  _Communications in Information and Systems_ , 11(3):197–224, 2011.
  * Davis and Vinter [1985] M.H.A Davis and R.B. Vinter.  _Stochastic modelling and control_.  Chapman and Hall, 1985.
  * Du et al. [2022] Zhe Du, Zexiang Liu, Jack Weitze, and Necmiye Ozay.  Sample complexity analysis and self-regularization in identification of over-parameterized ARX models.  In _2022 IEEE 61st Conference on Decision and Control (CDC)_ , pages 6026–6033. IEEE, 2022.
  * Faradonbeh et al. [2018] Mohamad Kazem Shirani Faradonbeh, Ambuj Tewari, and George Michailidis.  Finite time identification in unstable linear systems.  _Automatica_ , 96:342–353, 2018.
  * Foster et al. [2020] Dylan Foster, Tuhin Sarkar, and Alexander Rakhlin.  Learning nonlinear dynamical systems from a single trajectory.  In _Learning for Dynamics and Control_ , pages 851–861. PMLR, 2020\.
  * Gevers [2005] Michel Gevers.  Identification for Control: From the Early Achievements to the Revival of Experiment Design.  _European journal of control_ , 11(4-5):335–352, 2005.
  * Hanson and Wright [1971] David Lee Hanson and Farroll Tim Wright.  A bound on tail probabilities for quadratic forms in independent random variables.  _The Annals of Mathematical Statistics_ , 42(3):1079–1083, 1971.
  * Jansson [2003] Magnus Jansson.  Subspace identification and ARX modeling.  _IFAC Proceedings Volumes_ , 36(16):1585–1590, 2003.
  * Jedra and Proutiere [2022] Yassir Jedra and Alexandre Proutiere.  Finite-time identification of linear systems: Fundamental limits and optimal algorithms.  _IEEE Transactions on Automatic Control_ , 2022.
  * Kowshik et al. [2021] Suhas Kowshik, Dheeraj Nagaraj, Prateek Jain, and Praneeth Netrapalli.  Near-optimal offline and streaming algorithms for learning non-linear dynamical systems.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Lale et al. [2021] Sahin Lale, Kamyar Azizzadenesheli, Babak Hassibi, and Anima Anandkumar.  Finite-time system identification and adaptive control in autoregressive exogenous systems.  In _Learning for Dynamics and Control_ , pages 967–979. PMLR, 2021\.
  * Lee and Lamperski [2020] Bruce Lee and Andrew Lamperski.  Non-asymptotic Closed-Loop System Identification using Autoregressive Processes and Hankel Model Reduction.  In _IEEE 59th Conference on Decision and Control (CDC)_ , 2020.
  * Lee [2022] Holden Lee.  Improved rates for prediction and identification of partially observed linear dynamical systems.  In _International Conference on Algorithmic Learning Theory_ , pages 668–698. PMLR, 2022.
  * Liang et al. [2015] Tengyuan Liang, Alexander Rakhlin, and Karthik Sridharan.  Learning with square loss: Localization through offset rademacher complexity.  In _Conference on Learning Theory_ , pages 1260–1285. PMLR, 2015\.
  * Ljung [1999] Lennart Ljung.  System identification: theory for the user.  _PTR Prentice Hall, Upper Saddle River, NJ_ , 28, 1999.
  * Mendelson [2014] Shahar Mendelson.  Learning without concentration.  In _Conference on Learning Theory_ , pages 25–39. PMLR, 2014.
  * Oymak and Ozay [2019] Samet Oymak and Necmiye Ozay.  Non-asymptotic identification of LTI systems from a single trajectory.  In _2019 American control conference (ACC)_ , pages 5655–5661. IEEE, 2019.
  * Oymak and Ozay [2021] Samet Oymak and Necmiye Ozay.  Revisiting Ho–Kalman-based system identification: Robustness and finite-sample analysis.  _IEEE Transactions on Automatic Control_ , 67(4):1914–1928, 2021.
  * Peña et al. [2009] Victor H Peña, Tze Leung Lai, and Qi-Man Shao.  _Self-normalized processes: Limit theory and Statistical Applications_.  Springer, 2009.
  * Pukelsheim [2006] Friedrich Pukelsheim.  _Optimal Design of Experiments_.  SIAM, 2006.
  * Qin [2006] S Joe Qin.  An overview of subspace identification.  _Computers & chemical engineering_, 30(10-12):1502–1513, 2006.
  * Rakhlin and Sridharan [2014] Alexander Rakhlin and Karthik Sridharan.  Online non-parametric regression.  In _Conference on Learning Theory_ , pages 1232–1264. PMLR, 2014\.
  * Robbins and Siegmund [1970] Herbert Robbins and David Siegmund.  Boundary crossing probabilities for the Wiener process and sample sums.  _The Annals of Mathematical Statistics_ , pages 1410–1429, 1970.
  * Rudelson and Vershynin [2013] Mark Rudelson and Roman Vershynin.  Hanson-Wright inequality and sub-gaussian concentration.  _Electronic Communications in Probability_ , 18, 2013.
  * Samson [2000] Paul-Marie Samson.  Concentration of measure inequalities for Markov chains and $\phi$-mixing processes.  _The Annals of Probability_ , 28(1):416–461, 2000\.
  * Sarkar and Rakhlin [2019] Tuhin Sarkar and Alexander Rakhlin.  Near optimal finite time identification of arbitrary linear dynamical systems.  In _International Conference on Machine Learning_ , pages 5610–5618, 2019.
  * Sarkar et al. [2021] Tuhin Sarkar, Alexander Rakhlin, and Munther A Dahleh.  Finite time LTI system identification.  _Journal of Machine Learning Research_ , 22(26):1–61, 2021.
  * Sattar and Oymak [2022] Yahya Sattar and Samet Oymak.  Non-asymptotic and accurate learning of nonlinear dynamical systems.  _The Journal of Machine Learning Research_ , 23(1):6248–6296, 2022.
  * Simchowitz et al. [2018] Max Simchowitz, Horia Mania, Stephen Tu, Michael I. Jordan, and Benjamin Recht.  Learning without mixing: Towards a sharp analysis of linear system identification.  In _Conference On Learning Theory_ , pages 439–473. PMLR, 2018.
  * Simchowitz et al. [2019] Max Simchowitz, Ross Boczar, and Benjamin Recht.  Learning linear dynamical systems with semi-parametric least squares.  In _Conference on Learning Theory_ , pages 2714–2802. PMLR, 2019\.
  * Tsiamis and Pappas [2019] Anastasios Tsiamis and George J. Pappas.  Finite sample analysis of stochastic system identification.  In _2019 IEEE 58th Conference on Decision and Control (CDC)_ , pages 3648–3654. IEEE, 2019.
  * Tsiamis and Pappas [2021] Anastasios Tsiamis and George J Pappas.  Linear systems can be hard to learn.  _arXiv preprint arXiv:2104.01120_ , 2021.
  * Tsiamis et al. [2023] Anastasios Tsiamis, Ingvar Ziemann, Nikolai Matni, and George J Pappas.  Statistical learning theory for control: A finite sample perspective.  _to appear: IEEE Control Systems Magazine_ , 2023.
  * Tu et al. [2022] Stephen Tu, Roy Frostig, and Mahdi Soltanolkotabi.  Learning from many trajectories.  _arXiv preprint arXiv:2203.17193_ , 2022.
  * van der Vaart [2000] Aad W van der Vaart.  _Asymptotic Statistics_.  Cambridge university press, 2000.
  * Verhaegen and Verdult [2007] Michel Verhaegen and Vincent Verdult.  _Filtering and system identification: a least squares approach_.  Cambridge university press, 2007.
  * Vershynin [2018] Roman Vershynin.  _High-Dimensional Probability: An Introduction with Applications in Data Science_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2018.  doi: 10.1017/9781108231596.
  * Wagenmaker and Jamieson [2020] Andrew Wagenmaker and Kevin Jamieson.  Active learning for identification of linear dynamical systems.  In _Conference on Learning Theory_ , pages 3487–3582. PMLR, 2020\.
  * Wainwright [2019] Martin J. Wainwright.  _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48.  Cambridge University Press, 2019.
  * Yu [1994] Bin Yu.  Rates of convergence for empirical processes of stationary mixing sequences.  _The Annals of Probability_ , 22(1):94–116, 1994\.
  * Ziemann [2022] Ingvar Ziemann.  A note on the smallest eigenvalue of the empirical covariance of causal gaussian processes.  _arXiv preprint arXiv:2212.09508_ , 2022.
  * Ziemann and Tu [2022] Ingvar Ziemann and Stephen Tu.  Learning with little mixing.  _arXiv preprint arXiv:2206.08269. NeurIPS’22_ , 2022.
  * Ziemann et al. [2023] Ingvar Ziemann, Stephen Tu, George J Pappas, and Nikolai Matni.  The noise level in linear regression with dependent data.  _arXiv preprint arXiv:2305.11165_ , 2023.

##  Appendix A Proof of The Hanson-Wright Inequality

Let us now proceed with the proof of the Hanson-Wright Inequality (Theorem
2.1). We will need a few intermediate exponential inequalities relating sub-
Gaussian variables to Gaussian counterparts.

###  A.1 Gaussian Comparison Inequalities for sub-Gaussian Quadratic Forms

######  Lemma A.1.

Fix symmetric positive semidefinite $M\in\mathbb{R}^{d\times d}$ and let $W$
be a $\sigma^{2}$-sub-Gaussian assuming values in $\mathbb{R}^{d}$. Let
further $G$ be an independent $N(0,I_{d})$-distributed random variable. For
every $\lambda\in\mathbb{R}$ with $|\lambda|\leq
1/(\sqrt{2}\sigma^{2}\|M\|_{\mathsf{op}})$:

$\mathbf{E}\exp\left(\lambda
W^{\mathsf{T}}MW\right)\leq\mathbf{E}\exp\left(\lambda\sigma^{2}G^{\mathsf{T}}MG\right)\leq\exp\left(\lambda^{2}\sigma^{4}\|M\|_{F}^{2}\right)$
(A.1)

###### Proof.

For the first inequality note that:

$\displaystyle\mathbf{E}\exp\left(\lambda W^{\mathsf{T}}MW\right)$
$\displaystyle=\mathbf{E}\exp\left(\sqrt{2\lambda}W^{\mathsf{T}}\sqrt{M}G\right)$
$\displaystyle(\textnormal{Gaussian MGF})$ (A.2)
$\displaystyle\leq\mathbf{E}\exp\left(\lambda\sigma^{2}G^{\mathsf{T}}MG\right).$
$\displaystyle(\textnormal{sub-Gaussian MGF})$

To obtain the second inequality, let $VSV^{\mathsf{T}}$ be the
eigendecomposition of $M$ with orthonormal $V$. Let $v_{i},i\in[d]$ be the
columns of $V$. Then

$G^{\mathsf{T}}MG=\sum_{i=1}^{d}s_{i}^{2}(v_{i}^{\mathsf{T}}G)^{2}$ (A.3)

Note by Gaussian rotational invariance that the vector $G^{\prime}$ with
entries $v_{i}^{\mathsf{T}}G$ is equal in distribution to $G$. Hence:

$\displaystyle\mathbf{E}\exp\left(\lambda\sigma^{2}G^{\mathsf{T}}MG\right)$
$\displaystyle=\prod_{i=1}^{d}\mathbf{E}\exp\left(\lambda\sigma^{2}s_{i}^{2}(v_{i}^{\mathsf{T}}G)^{2}\right)$
(A.4)
$\displaystyle=\prod_{i=1}^{d}\mathbf{E}\exp\left(\lambda\sigma^{2}s_{i}^{2}G_{i}^{2}\right)$
$\displaystyle(\textnormal{Gaussian Rotation Invariance})$
$\displaystyle=\prod_{i=1}^{d}\left(1-\lambda^{2}\sigma^{4}s_{i}^{2}\right)^{-1/2}$
$\displaystyle(\textnormal{Gaussian-Squared MGF; $\lambda$ in our range})$
$\displaystyle=\prod_{i=1}^{d}\exp\left(\frac{-1}{2}\log(1-\lambda^{2}\sigma^{4}s_{i}^{2})\right)$
$\displaystyle\leq\prod_{i=1}^{d}\exp\left(\lambda^{2}\sigma^{4}s_{i}^{2}\right)$
$\displaystyle\left(-\log(1-x)\leq 2x\textnormal{ if }x\in[0,1/2]\right)$
$\displaystyle=\exp\left(\lambda^{2}\sigma^{4}\|M\|_{F}^{2}\right)$
$\displaystyle\left(\sum_{i=1}^{d}s_{i}^{2}=\|M\|_{F}^{2}\right)$

as per requirement. ∎

######  Lemma A.2 (Comparison to Gaussians).

Fix $M\in\mathbb{R}^{d\times d}$ and let $W,W^{\prime}$ be independent
$\sigma^{2}$-sub-Gaussian assuming values in $\mathbb{R}^{d}$. Let further
$G,G^{\prime}$ be two independent $N(0,I_{d})$-distributed random variables.
For every $\lambda\in\mathbb{R}$:

$\mathbf{E}\exp\left(\lambda
W^{\mathsf{T}}MW^{\prime}\right)\leq\mathbf{E}\exp\left(\lambda{\sigma^{2}}G^{\mathsf{T}}MG^{\prime}\right).$

###### Proof.

The result follows by straightforward computation, alternatingly using the
sub-Gaussian property and the closed form of the Gaussian MGF. Namely:

$\displaystyle\mathbf{E}\exp\left(\lambda W^{\mathsf{T}}MW^{\prime}\right)$
$\displaystyle\leq\mathbf{E}\exp\left(\frac{\sigma^{2}\lambda^{2}\|MW^{\prime}\|_{2}^{2}}{2}\right)$
$\displaystyle(W\textnormal{ is }\sigma^{2}-\textnormal{subG})$ (A.5)
$\displaystyle=\mathbf{E}\exp\left(\lambda\sigma
G^{\mathsf{T}}MW^{\prime}\right)$ $\displaystyle(\textnormal{MGF of }G)$
$\displaystyle\leq\mathbf{E}\exp\left(\frac{\lambda^{2}\sigma^{4}\|MG\|_{2}^{2}}{2}\right)$
$\displaystyle(W^{\prime}\textnormal{ is }\sigma^{2}-\textnormal{subG})$
$\displaystyle=\mathbf{E}\exp\left(\lambda{\sigma^{2}}G^{\mathsf{T}}MG^{\prime}\right).$
$\displaystyle(\textnormal{MGF of }G^{\prime})$

The desired result has been established. ∎

######  Lemma A.3 (MGF of Gaussian Chaos).

Let $W,W^{\prime}\sim N(0,I_{d})$ be independent and let
$M\in\mathbb{R}^{d\times d}$. Then for every $\lambda\in\mathbb{R}$ with
$|\lambda|\leq 1/(\sqrt{2}\|M\|_{\mathsf{op}})$:

$\mathbf{E}\exp\left(\lambda
W^{\mathsf{T}}MW^{\prime}\right)\leq\exp\left(\lambda^{2}\|M\|_{F}^{2}\right)$
(A.6)

###### Proof.

By the singular value decomposition we may write $M=USV^{\mathsf{T}}$ where
the columns $u_{1:d}$ ($v_{1:d}$) of $U$ (of $V$) form orthonormal bases of
$\mathbb{R}^{d}$. Hence

$W^{\mathsf{T}}MW=\sum_{i=1}^{d}s_{i}(W^{\mathsf{T}}u_{i})((W^{\prime})^{\mathsf{T}}v_{i})$
(A.7)

where $s_{1:d}$ are the singular values of $M$. Note next that
$G_{i}\triangleq(W^{\mathsf{T}}u_{i})$ and
$G_{i}^{\prime}\triangleq((W^{\prime})^{\mathsf{T}}v_{i})$ are again
independent and standard normal by Gaussian rotational invariance. Hence:

$\displaystyle\mathbf{E}\exp\left(\lambda W^{\mathsf{T}}MW\right)$
$\displaystyle=\mathbf{E}\exp\left(\lambda\sum_{i=1}^{d}s_{i}G_{i}G_{i}^{\prime}\right)$
(A.8) $\displaystyle=\prod_{i=1}^{d}\mathbf{E}\exp\left(\lambda
s_{i}G_{i}G_{i}^{\prime}\right)$ $\displaystyle(\textnormal{independence})$
$\displaystyle=\prod_{i=1}^{d}\mathbf{E}\exp\left(\lambda^{2}s_{i}^{2}G_{i}^{2}/2\right)$
$\displaystyle(\textnormal{Gaussian MGF})$
$\displaystyle=\prod_{i=1}^{d}\left(1-\lambda^{2}s_{i}^{2}\right)^{-1/2}$
$\displaystyle(\textnormal{Gaussian-Squared MGF; $\lambda$ in our range})$
$\displaystyle=\prod_{i=1}^{d}\exp\left(\frac{-1}{2}\log(1-\lambda^{2}s_{i}^{2})\right)$
$\displaystyle(\exp\circ\log\textnormal{is the identity function})$
$\displaystyle\leq\prod_{i=1}^{d}\exp\left(\lambda^{2}s_{i}^{2}\right)$
$\displaystyle\left(-\log(1-x)\leq 2x\textnormal{ if }x\in[0,1/2]\right)$
$\displaystyle=\exp\left(\lambda^{2}\|M\|_{F}^{2}\right)$
$\displaystyle\left(\sum_{i=1}^{d}s_{i}^{2}=\|M\|_{F}^{2}\right)$

as per requirement. ∎

###  A.2 Finishing the proof of Theorem 2.1

We need one more preliminary tool before we arrive at the proof of Theorem
2.1. The next Theorem constitutes a useful decoupling inequality which allows
us to treat the mixed terms in the quadratic form $W^{\mathsf{T}}MW$ as
independent.

######  Theorem A.1 (Theorem 6.1.1 in of Vershynin [2018]).

Let $W$ be a $d$-dimensional random vector with mean zero and independent
entries. For every convex function $f$ and every
$M=(m_{ij})_{i,j=1}^{d}\in\mathbb{R}^{d\times d}$ it holds true that:

$\mathbf{E}f\left(\sum_{i,j=1,i\neq
j}^{d}m_{ij}W_{i}W_{j}\right)\leq\mathbf{E}f\left(4\sum_{i,j=1}^{d}m_{ij}W_{i}W_{j}^{\prime}\right)$
(A.9)

where $W^{\prime}$ is an independent copy of $W$ (i.e., equal to $W$ in
distribution but independent of $W$).

######  Proposition A.1 (Hanson-Wright Exponential Inequality Form).

Let $M\in\mathbb{R}^{d\times d}$. Fix a random variable $W=W_{1:d}$ where each
$W_{i},i\in[d]$ is an independent scalar $\sigma^{2}$-sub-Gaussian random
variable. For every $\lambda\in\mathbb{R}$ with
$|\lambda|\leq\frac{1}{8\sqrt{2}\sigma^{2}\|M\|_{\mathsf{op}}}$ we have that:

$\max\left\\{\mathbf{E}\exp\left(\lambda
W^{\mathsf{T}}MW-\lambda\mathbf{E}W^{\mathsf{T}}MW\right),\mathbf{E}\exp\left(\lambda
W^{\mathsf{T}}MW\right)\right\\}\leq\exp\left(36\lambda^{2}\sigma^{4}\|M\|_{F}^{2}\right)$
(A.10)

###### Proof.

By rescaling $\lambda$ we may assume without loss of generality that
$\sigma^{2}=1$. Let now $m_{ij}$ with $i,j\in[d]$ denote the entries of $M$.
We begin by writing

$W^{\mathsf{T}}MW-\mathbf{E}W^{\mathsf{T}}MW=\sum_{i=1}^{d}m_{ii}(W_{i}^{2}-\mathbf{E}W_{i}^{2})+\sum_{i\neq
j}m_{ij}W_{i}W_{j}.$ (A.11)

Hence by the Cauchy-Schwarz inequality:

$\mathbf{E}\exp\left(\lambda