# Finite Sample Frequency Domain Identification

Anastasios Tsiamis  ETH Zürich, Switzerland  Mohamed Abdalmoaty  ETH Zürich,
Switzerland  Roy S. Smith  ETH Zürich, Switzerland  John Lygeros  ETH Zürich,
Switzerland

###### Abstract

We study non-parametric frequency-domain system identification from a finite-
sample perspective. We assume an open loop scenario where the excitation input
is periodic and consider the Empirical Transfer Function Estimate (ETFE),
where the goal is to estimate the frequency response at certain desired
(evenly-spaced) frequencies, given input-output samples. We show that under
sub-Gaussian colored noise (in time-domain) and stability assumptions, the
ETFE estimates are concentrated around the true values. The error rate is of
the order of
$\mathcal{O}((d_{{u}}+\sqrt{d_{{u}}d_{{y}}})\sqrt{M/N_{\mathrm{tot}}})$, where
$N_{\mathrm{tot}}$ is the total number of samples, $M$ is the number of
desired frequencies, and $d_{{u}},\,d_{{y}}$ are the dimensions of the input
and output signals respectively. This rate remains valid for general
irrational transfer functions and does not require a finite order state-space
representation. By tuning $M$, we obtain a $N_{\mathrm{tot}}^{-1/3}$ finite-
sample rate for learning the frequency response over all frequencies in the
$\mathcal{H}_{\infty}$ norm. Our result draws upon an extension of the Hanson-
Wright inequality to semi-infinite matrices. We study the finite-sample
behavior of ETFE in simulations.

##  1 Introduction

We consider the identification of _unknown_ linear, discrete-time, time-
invariant systems of the form

$\displaystyle y_{t}$ $\displaystyle=\bar{y}_{t}+v_{t},$ (1)
$\displaystyle\bar{y}_{t}$ $\displaystyle=G(\mathrm{\mathstrut
q})u_{t},\,G(\mathrm{\mathstrut
q})\triangleq\sum_{s=0}^{\infty}g_{s}\mathrm{\mathstrut q}^{-s},$

where $t\in\mathbb{Z}$ is the time, $u_{t}\in\mathbb{R}^{d_{{u}}}$ is the
input, $y_{t}\in\mathbb{R}^{d_{{y}}}$ is the output, $\mathrm{\mathstrut
q}^{-s}u_{t}=u_{t-s}$ is the backward shift operator, and
$g_{t}\in\mathbb{R}^{d_{{y}}\times d_{{u}}}$ is the impulse response. The
noiseless output $\bar{y}_{t}$ is perturbed by some random noise process
$v_{t}\in\mathbb{R}^{d_{{y}}}$. We are interested in estimating the frequency
response $G(e^{j\omega})$ from finite input-output data.

Frequency domain identification has been extensively studied [1, 2, 3]. The
estimation error guarantees (on its distribution) are typically asymptotic,
e.g. see Central Limit Theorem in [3, Ch. 16], and, thus, are valid when the
number of samples grows to infinity. Here, we adopt a finite-sample point of
view, motivated by advances in modern statistics [4, 5] and statistical
learning theory. Asymptotic methods are sharp asymptotically but are often
heuristically applied for finite samples. Finite-sample bounds, on the other
hand, are valid for any number of samples, but suffer from looser bounding
constants. Nonetheless, they can provide a more detailed qualitative
characterization of the statistical difficulty of learning [6].

While finite-sample system identification has been studied before, most
results are focused on time domain identification [7, 8, 9, 10, 11, 12, 13,
14, 15, 16]. Detailed related work and a tutorial on the subject can be found
in [6, 17]. Frequency domain and time domain identification have many
similarities–ignoring initial conditions, transients, or leakage effects, the
two domains are equivalent from a prediction error framework perspective [2].
Still, working in one domain may offer some advantages over the other [2]. For
example, the frequency domain approach allows a unified treatment of discrete
and continuous time systems, simplifies the analysis of systems with delays,
and offers a more explicit way of designing the input excitation.

Our contributions are the following:

Finite-sample guarantees for the ETFE. We provide finite sample guarantees for
the well-established Empirical Transfer Function Estimate (ETFE) [1], a non-
parametric method for frequency domain identification, under open-loop
periodic excitation. While the mean and variance of the ETFE have been
characterized before, we provide guarantees on the distribution of the
estimation error, the tail probabilities in particular. Under certain
stability conditions, we prove that the estimation error decays with a rate of
$\sqrt{M/N_{\mathrm{tot}}}$, where $N_{\mathrm{tot}}$ is the total number of
samples. The parameter $M$ is the number of selected frequencies at which we
estimate the frequency response; it controls the frequency resolution. The
rate holds for general irrational transfer functions and does not require a
finite order state-space representation, unlike prior non-asymptotic bounds
[11].

Guarantees in the $\mathcal{H}_{\infty}$ norm. Based on our finite-sample
bound, we tune the number of frequencies $M$ to provide guarantees for
learning the frequency response across all frequencies. We provide a non-
asymptotic rate of $N_{\mathrm{tot}}^{-1/3}$ in the $\mathcal{H}_{\infty}$
norm of the estimation error, which reflects optimal rates for non-parametric
learning of Lipschitz functions [18].

Extension of the Hanson-Wright inequality. To prove our main result we have to
deal with quadratic forms of a (countably) infinite number of sub-Gaussian
variables. To achieve this, we extend the celebrated Hanson-Wright inequality
[4, 19] to semi-infinite matrices; that is, bounded operators mapping
sequences to finite vector spaces.

Our paper is related to non-parametric system identification, which includes
works on both time [20] and frequency domain [21]. Using Gaussian Processes as
in [21], where the unknown frequency response follows a Gaussian prior, we can
also obtain finite sample guarantees. Here, we follow a different approach and
we do not consider Gaussian priors. Note that in this work we focus on
qualitative data-independent bounds linking sample requirements to system
theoretic properties. Application-oriented, data-dependent bounds have also
been studied before [22].

Notation. Let $(H,\mathbb{F},\langle\cdot,\cdot\rangle_{H})$ be a Hilbert
space with field $\mathbb{F}=$ or $\mathbb{C}$ and inner product
$\langle\cdot,\cdot\rangle_{H}$. For any vector $x\in H$, let
$\|x\|\triangleq\sqrt{\langle x,x\rangle_{H}}$ denote the inner product norm.
Let $H,V$ be Hilbert spaces with $\mathbb{F}=$ or $\mathbb{C}$ and let
$\mathcal{A}:H\rightarrow V$ be any linear map. Let
$\|\mathcal{A}\|_{\mathrm{op}}\triangleq\sup_{\|x\|=1}\|\mathcal{A}(x)\|$
denote the operator norm and $\mathcal{A}^{*}$ denote the adjoint operator. If
$\left\\{b_{i}\right\\}_{i\in\mathcal{I}}$ is an orthonormal basis, the
Hilbert-Schmidt or Frobenius norm is defined as
$\|\mathcal{A}\|_{\mathrm{F}}^{2}\triangleq\sum_{i\in\mathcal{I}}\|\mathcal{A}(b_{i})\|^{2}$.
Let
$\mathcal{L}_{2}(p)\triangleq\\{\left\\{x_{k}\right\\}^{\infty}_{k=1}:x_{k}\in^{p},\,\sum_{k=1}^{\infty}\|x_{k}\|_{2}^{2}<\infty\\}$
be the Hilbert space of $p$-dimensional square summable sequences. A universal
constant is a constant that is independent of the problem at hand, e.g., the
system or the algorithm. For any integer $M$, let
$[M]\triangleq{0,\dots,M-1}$. The $\mathcal{H}_{\infty}$ norm of
$G(e^{j\omega}):[0,2\pi)\rightarrow\mathbb{C}^{d_{1}\times d_{2}}$ is given by
$\sup_{\omega\in[0,2\pi)}\|G(\omega)\|_{\mathrm{op}}$; it is denoted by
$\|G\|_{\mathcal{H}_{\infty}}$.

##  2 Problem formulation

Consider the input-output system (1). We make the following assumption about
the noise process $v_{t}$.

######  Assumption 1 (Noise).

The noise process $v_{t}$ is filtered sub-Gaussian white noise, that is,

$v_{t}=H(\mathrm{\mathstrut q})e_{t},\quad H(\mathrm{\mathstrut
q})\triangleq\sum_{s=0}^{\infty}h_{s}\mathrm{\mathstrut q}^{-s},$ (2)

where $h_{t}\in\mathbb{R}^{d_{{y}}\times d_{{e}}}$ are the _unknown_ filter
coefficients. Let $e_{t}\in\mathbb{R}^{d_{{e}}}$ be i.i.d. zero mean, with
covariance $\mathbb{E}{e_{t}e_{t}^{\top}}=\sigma_{e}^{2}I_{d_{{e}}}$, and
$K^{2}$-sub-Gaussian, i.e., for any $\xi\in^{d_{{e}}}$

$\mathbb{E}\exp\left(\xi^{\top}e_{t}\right)\leq\exp\left(\frac{K^{2}\|\xi\|^{2}}{2}\right),$
(3)

for some $K>0$.

The noise process $v_{t}$ is colored. It is used to model measurement noise as
well as any stochastic disturbances acting on the dynamical system.

We assume throughout that the input is bounded. This guarantees that any
transient phenomena have a limited effect on the estimation problem.

######  Assumption 2 (Input Bound).

All inputs are bounded

$\left\|u_{t}\right\|_{2}\leq D_{{u}},\,\text{ for all }t\in\mathbb{Z}$

for some $D_{u}>0$ independent of $t$.

We start all identification experiments at time $t=0$. Hence, the initial
conditions are determined by all past signals $u_{-1},u_{-2},\dots$ and
$e_{-1},e_{-2},\dots,$, which are nonzero in general, and unknown. Note that
our formulation allows general irrational transfer functions and does not
assume a state-space representation of finite dimension.

###  2.1 Empirical Transfer Function Estimate

The goal of non-parametric frequency domain identification is to estimate the
frequency response $G(e^{j\omega})$, given input-output data. We assume access
to $d_{{u}}$ experiments of length $N$, that is, data
$(u^{(i)}_{0},y^{(i)}_{0},\dots,u^{(i)}_{N-1},y^{(i)}_{N-1})$, for
$i=1,\dots,d_{{u}}$. This brings the total number of samples to
$N_{\mathrm{tot}}\triangleq d_{{u}}N$. For simplicity, we assume that the
trajectories are statistically independent. Our analysis can be extended to
the single trajectory case, where we gather all samples sequentially.

We are interested in the performance of the ETFE, which we review here. Given
any signal $z=\\{z_{t}\\}_{t\in[N]}$, let

$Z_{k}\triangleq\mathcal{F}^{N}_{k}(z)\triangleq\frac{1}{\sqrt{N}}\sum_{t=0}^{N-1}z_{t}e^{-j\frac{2\pi
k}{N}t},\,k\in[N]$

denote its $N$-point Discrete Fourier Transform (DFT), evaluated at $k$. Let
$Y^{(i)}_{k},U^{(i)}_{k}$ be the $N$ point DFTs of $y^{(i)}_{t}$ and
$u^{(i)}_{t}$ respectively for the $i-$th experiment, $i=1,\dots,d_{{u}}$. Let
$Y_{k}\in\mathbb{C}^{d_{{y}}\times d_{{u}}},U_{k}\in\mathbb{C}^{d_{{u}}\times
d_{{u}}}$ denote the stacked DFTs for all experiments

$\\!\\!\\!Y_{k}\triangleq\begin{bmatrix}Y^{(1)}_{k}&\cdots&Y^{(d_{{u}})}_{k}\end{bmatrix},\,\\!U_{k}\triangleq\begin{bmatrix}U^{(1)}_{k}&\cdots&U^{(d_{{u}})}_{k}\end{bmatrix}\\!.$
(4)

Then, an estimate of $G(e^{j\omega})$ at frequency $\omega_{k}=2\pi k/N$, for
$k=0,\dots,N-1$, can be obtained using the _ETFE_

$\hat{G}_{k}\triangleq Y_{k}U^{-1}_{k},$ (5)

provided that $U_{k}$ is invertible; the estimate is undefined if not. Since
the number of frequencies $N$ scales with the number of data, it is generally
impossible to estimate the responses at all frequencies consistently (without
assuming structure) [1]. Instead, we can learn the responses at a smaller
frequency set. Given a frequency-resolution parameter $M<N$, we focus on
estimating $G(e^{j\omega})$ at $\omega=2\pi\ell/M$, for $\ell\in[M]$. Based on
the ETFE and under some additional Lipschitz assumptions on the frequency
response, we can extend the estimation over all frequencies
$\omega\in[0,2\pi)$.

###  2.2 Excitation Method

The estimation performance also depends on the excitation method. Since we
only need to estimate the frequency responses at $2\pi\ell/M$, $\ell\in[M]$,
it is sufficient to excite the system at only these frequencies [1]. Assuming
that $M$ divides $N$, the DFT $U_{k}$ of the input can be non-zero at only
$2\pi k/N=2\pi\ell/M$ or $k=\ell N/M$. The latter condition is satisfied if
and only if the excitation input is periodic with a period equal to $M$. Note
that we also need invertibility of $U_{k}$ at $k=\ell N/M$. To achieve this,
we assume the following.

######  Assumption 3 (Excitation).

Let the input signals be periodic with period $M$ such that
$u^{(i)}_{t+M}=u^{(i)}_{t}$, for $t\geq 0$ and every experiment
$i=1,\dots,d_{{u}}$. Assume that $M$ divides $N$ with $N_{p}\triangleq N/M$.
Consider _one period_ of the input signals and let the respective $M$-point
DFTs be

$\tilde{U}^{(i)}_{\ell}=\mathcal{F}^{M}_{\ell}(u^{(i)})=\frac{1}{\sqrt{M}}\sum_{t=0}^{M-1}u^{(i)}_{t}e^{-\frac{2\pi\ell}{M}t},i=1,\dots,d_{{u}}$

for $\ell\in[M]$, with respective stacked DFTs

$\tilde{U}_{\ell}\triangleq\begin{bmatrix}\tilde{U}^{(1)}_{\ell}&\cdots&\tilde{U}^{(d_{{u}})}_{\ell}\end{bmatrix}.$

Assume that for all $\ell\in[M]$ the stacked DFTs satisfy

$\sigma_{u,\ell}^{2}I_{d_{{u}}}\preceq\tilde{U}_{\ell}\tilde{U}^{*}_{\ell},$
(6)

for some $0<\sigma^{2}_{u,\ell}$ such that
$\sum_{\ell=0}^{M-1}\sigma^{2}_{u,\ell}\leq MD_{{u}}^{2}$, where $D_{{u}}$ is
the input upper bound of Assumption 2.

Such assumptions are standard when dealing with experiment design in frequency
domain. For example, Assumption 3 is satisfied by design (with uniform
$\sigma^{2}_{u,\ell}$ across $\ell\in[M]-\\{0\\}$) when pseudorandom binary
sequence (PRBS) signals are used and we excite one input at a time [1, Ch.
13]. Another choice could be multisine signals [23], where the user simply
designs the input to have sinusoids with non-zero amplitudes at the required
frequencies. Another option is to design the input spectrum and generate the
input by passing a white noise realization through the spectral factor [24].

###  2.3 Objective

We can now state our objective, which is providing finite-sample guarantees
for estimating the frequency responses. We focus on $\epsilon-\delta$
probabilistic guarantees, where $\epsilon$ controls the estimation accuracy
and $\delta$ controls the confidence.

######  Problem 1 (Finite-Sample ETFE).

Fix a frequency resolution $M<N$ such that $M$ divides $N$ and denote their
ratio by $N_{p}=N/M$. Consider $d_{{u}}$ independent input-ouput trajectories
of length $N$ $\\{u^{(i)}_{t},y^{(i)}_{t}\\}_{t\in[N]}$, for
$i=1,\dots,d_{{u}}$, generated by system (1) with excitation inputs as in
Assumption 3. Determine $\epsilon_{\ell}>0$, $\ell\in[M]$ and $\delta>0$ such
that $\mathbb{P}\left(\|G(e^{j\frac{2\pi\ell}{M}})-\hat{G}_{\ell
N_{p}}\|_{\mathrm{op}}\geq\epsilon_{\ell},\,\forall\ell\in[M]\right)\leq\delta,$
where the ETFE $\hat{G}_{\ell N_{p}}$ is defined in (5). Problem 1 only
focuses on the desired discretized frequency grid
$\\{2\pi\ell/M\\}_{\ell\in[M]}$. In Section 4, we also study uniform
guarantees over all frequencies in the $\mathcal{H}_{\infty}$ norm.

To guarantee a well-defined estimation problem, we consider the following
stability conditions.

######  Assumption 4 (Strict Stability).

The input-output impulse response is strictly stable [1], that is,

$\|G\|_{\star}\triangleq\sum_{t=0}^{\infty}t\|g_{t}\|_{\mathrm{op}}<\infty.$
(7)

The auto-correlation function of the noise
$R_{t}\triangleq\mathbb{E}v_{s}v_{s-t}^{\top}$ is also strictly stable

$\|R\|_{\star}\triangleq\sum_{t=0}^{\infty}t\|R_{t}\|_{\mathrm{op}}<\infty.$
(8)

Strict stability guarantees that the derivative of the frequency response
$\partial G(e^{j\omega})/\partial\omega$ is uniformly bounded over all
frequencies. This, in turn, implies that the response $G(e^{j\omega})$ is
Lipschitz. Strict stability also guarantees that the transient phenomena have
a limited effect on the estimation procedure.

##  3 Finite-sample guarantees for the ETFE

In this section, we focus on estimating the frequency response at the selected
frequencies $\\{2\pi\ell/M\\}_{\ell\in[M]}$. Following the convention of (4),
we define the stacked DFTs of the noises and the noiseless outputs as

$V_{k}\triangleq\begin{bmatrix}V^{(1)}_{k}&\cdots&V^{(d_{{u}})}_{k}\end{bmatrix},\,\bar{Y}_{k}\triangleq\begin{bmatrix}\bar{Y}^{(1)}_{k}&\cdots&\bar{Y}^{(d_{{u}})}_{k}\end{bmatrix}$

Then, for every frequency $\omega_{k}=2\pi k/N$ we have

$Y_{k}=\bar{Y}_{k}+V_{k}=G(e^{j\omega_{k}})U_{k}+V_{k}+T_{k,N},$ (9)

where $T_{k,N}=\bar{Y}_{k}-G(e^{j\omega_{k}})U_{k}$ accounts for transient and
time-aliasing phenomena since the DFT of $\\{\bar{y}^{(i)}_{t}\\}_{t=0}^{N-1}$
is different from $\\{G(e^{j\omega_{k}})U^{(i)}_{k}\\}_{k=0}^{N-1}$ for finite
$N$. This term vanishes as $N$ grows to infinity.

######  Remark 1.

The above relation fits the framework of non-parametric function estimation.
However, there are some notable differences with standard formulations [18,
5]. First, we have the presence of the input $U_{k}$, which affects the
signal-to-noise ratio (SNR) and is an additional degree of freedom. For
example, if the input matrix is not invertible at some $k$, we do not get a
well-defined sample of $G(e^{j\omega_{k}})$. Second, the noise $V_{k}$ is
heteroscedastic since its variance depends on the frequency $k$. Moreover, the
sequence $V_{k},\,k\in[N]$ is non-Gaussian and non-independent across
frequencies for finite samples $N$ (only asymptotically as $N$ goes to
infinity). Hence, the non-asymptotic techniques of [5, Ch. 13] do not apply
directly.

The estimation error is equal to

$\displaystyle\hat{G}_{k}-G(e^{j\omega_{k}})$
$\displaystyle=V_{k}U^{-1}_{k}+T_{k,N}U^{-1}_{k},$ (10)

where the input matrix $U_{k}$ is invertible, and we only look at the
frequencies $k=\ell N_{p}$, $\ell\in[M]$. Let
$\Phi_{v,N}(k)\triangleq\mathbb{E}V^{(i)}_{k}(V^{(i)}_{k})^{*}$ be the aliased
power spectrum of the process $v_{t}$ at frequency $k$, where due to
independence, the experiment index $i$ does not affect the definition. Define
the signal-to-noise ratio (SNR) at frequency $k=\ell N_{p},$ $\ell\in[M]$ as

$\mathrm{SNR}_{k,N}\triangleq\frac{\sigma_{u,\ell}}{\sqrt{\|\Phi_{v,N}(k)\|_{\mathrm{op}}}},$
(11)

where $\|\Phi_{v,N}(k)\|_{\mathrm{op}}$ is interpreted as the matrix norm for
fixed $k$. We obtain the following finite-sample guarantees.

######  Theorem 1 (ETFE Finite-Sample).

Consider Problem 1 and fix a failure probability $\delta>0$. Under Assumptions
1-4, with probability at least $1-\delta$ for all $k=\ell N_{p}$, $\ell\in[M]$

$\displaystyle\|G(e^{j\omega_{k}})-\hat{G}_{k}\|_{\mathrm{op}}\leq\frac{2\|G\|_{\star}D_{u}\sqrt{M}}{\sigma_{u,\ell}N}$
(12)
$\displaystyle+\frac{\sqrt{M}}{\sqrt{N}}\mathrm{SNR}^{-1}_{k,N}\left(\sqrt{d_{{y}}}+c\frac{K^{2}}{\sigma_{{e}}^{2}}\sqrt{d_{{u}}+\log
M/\delta}\right)$

where $c$ is a universal constant, $\|G\|_{\star}$ is defined in (7), and
$D_{{u}}$ is the maximum input norm.

The first term of the right-hand side captures the transient error
$T_{k,N}U^{-1}_{k}$, while the second one captures the error $V_{k}U^{-1}_{k}$
due to stochastic noise. Recall that the total number of samples is equal to
$N_{\mathrm{tot}}=d_{{u}}N$. As we increase the number of samples
$N_{\mathrm{tot}}$ while keeping $M$ constant, the former term decays at a
faster rate of $1/N_{\mathrm{tot}}$ compared to the latter’s
$1/\sqrt{N_{\mathrm{tot}}}$. Hence, the non-asymptotic rate is

$\mathcal{O}\left(\frac{\sqrt{M}}{\sqrt{N_{\mathrm{tot}}}}\sqrt{d_{{u}}}\big{(}\sqrt{d_{{u}}+\log
M/\delta}+\sqrt{d_{{y}}}\big{)}\right).$

The rate is similar to the ones for non-asymptotic parametric identification
in time-domain [17]; the optimal rate in that line of work is typically of the
order $\sqrt{d}$ for some $d$ scaling with the number of unknown parameters.
Here, we have a similar scaling of $\sqrt{M}(\sqrt{d_{{u}}}+\sqrt{d_{{y}}})$
times an additional $\sqrt{d_{{u}}}$ dimensional dependence. This is an
artifact of imposing a strict input norm bound in Assumption 2. If we allow
$\sigma_{{u}},D_{{u}}$ to scale with $\sqrt{d_{{u}}}$ (as is the case for
white-noise inputs in the time-domain [17]), we can remove this extra term.

A benefit of frequency-domain identification is that it provides specialized
guarantees for every frequency of interest by breaking down the SNR into SNRs
for every frequency. This offers direct insights on which frequencies to focus
on and how to design the excitation inputs. Note that the inverse
$\mathrm{SNR}^{-1}_{k,N}$ is upper bounded and converges to a limit as $N$
grows to infinity. This is a consequence of the following result which
exploits the strict stability condition (8).

######  Lemma 1 (Stochastic Transient [1]).

Denote the power spectrum of the noise at frequency $\omega_{k}=2\pi k/N$, for
some $k=0,\dots,N-1$, by
$\Phi_{v}(k)\triangleq\sum_{t=-\infty}^{\infty}R_{t}e^{-j\omega_{k}t}$, with
$R_{-t}=R_{t}^{\top}$. We have

$\|\Phi_{v}(k)-\Phi_{v,N}(k)\|_{\mathrm{op}}\leq 2\|R\|_{\star}/N.$

In the remainder of the section, we provide a sketch of the proof of Theorem
1. We simply bound every term that appears in (10) separately.

###  3.1 Deterministic transient

Even in the absence of any stochastic noise, the ETFE suffers from estimation
errors due to transient phenomena (e.g. aliasing, leakage) [2]. Fortunately,
the transient error $T_{k,N}=\bar{Y}_{k}-G(e^{j\omega_{k}})U_{k}$ decays
uniformly to zero as the DFT horizon $N$ goes to infinity.

######  Lemma 2 (Deterministic Transient [1]).

We have

$\sup_{k=0,\dots,N-1}\|T_{k,N}\|_{\mathrm{op}}\leq
2\|G\|_{\star}D_{{u}}/\sqrt{N}.$

The above result is a consequence of bounded inputs and strict stability
(Assumptions 2, 4). Together they guarantee that the deterministic transient
phenomena have a vanishing effect on the estimation procedure.

###  3.2 Input energy

Next, we review a standard result for periodic inputs. The input $U_{k}$ at
frequencies $k=\ell N_{p}$ is equal to $\sqrt{N_{p}}\tilde{U}_{\ell}$, where
$\tilde{U}_{\ell}$ is the $M-$point DFT based on one period of the input
signals (see Assumption 3). As a result, $U^{-1}_{k}$ is well-defined and
vanishes to zero with $N^{-1/2}_{p}$. This phenomenon is a direct consequence
of periodicity and the properties of DFT.

######  Lemma 3 (Input energy).

Let Assumption 3 be in effect. Then, for all $k=N_{p}\ell$, $\ell\in[M]$

$\sigma_{u,\ell}^{2}N_{p}I_{d_{{u}}}\preceq U_{k}U^{*}_{k},\text{ for }k=\ell
N_{p}.$ (13)

Lemma 3 is key to achieving consistency. While the known input is periodic,
the noise is not. Hence, at the selected frequencies, the noise is averaged
out.

###  3.3 Noise concentration

Finally, we bound the noise term $V_{k}$ by showing that its norm concentrates
around $\sqrt{\mathrm{tr}(\Phi_{v,N}(k))}$.

######  Lemma 4 (Concentration of noise).

Under Assumption 1, there exists a universal constant $c$ such that for any
$s>0$

$\mathbb{P}\left\\{\left|\|V_{k}\|_{\mathrm{op}}-\sqrt{\mathrm{tr}(\Phi_{v,N}(k))}\right|>2s\sqrt{\mathrm{tr}(\Phi_{v,N}(k))}\right\\}\leq
9^{2d_{{u}}}2\exp\left(-cs^{2}\frac{\sigma_{{e}}^{4}}{K^{4}}\frac{\mathrm{tr}(\Phi_{v,N}(k))}{\|\Phi_{v,N}(k)\|_{\mathrm{op}}}\right).$
(14)

Prior results show that $V^{(i)}_{k}$ is asymptotically circular Gaussian as
$N$ goes to infinity [3, Ch. 16]. We show that even under finite $N$, the
tails of the distribution decay exponentially reflecting the properties of the
Gaussian distribution.

The proof of Lemma 4 is based on a novel extension of the Hanson-Wright
inequality [19], a standard tool for proving concentration of quadratic forms
involving random variables. Note that $V_{k}$ is a linear function of an
infinite number of random noises $e^{(i)}_{t}$, $i=1,\dots,d_{{u}}$,
$-\infty<t\leq N-1$. Hence, we need to extend the Hanson-Wright inequality
[19, 25] to semi-infinite matrices, that is, bounded operators from square
summable (real-valued) sequences $\mathcal{L}_{2}$ to finite vector spaces.

######  Theorem 2 (Semi-infinite Hanson-Wright).

Consider a sequence $z\triangleq\\{z_{t}\\}_{t=1}^{\infty},\,z_{t}\in^{p}$ of
independent, zero-mean, $K^{2}$-sub-Gaussian random variables with covariance
$\mathbb{E}z_{t}z^{\top}_{t}=\sigma_{{e}}^{2}I_{p}$. Let
$\mathcal{A}:\mathcal{L}_{2}(p)\rightarrow^{d}$ be a linear map with bounded
Frobenius norm $\|\mathcal{A}\|_{\mathrm{F}}<\infty$. For any $\alpha>0$

$\mathbb{P}\left\\{\left|\|\mathcal{A}(z)\|^{2}-\mathbb{E}\|\mathcal{A}(z)\|^{2}\right|>\alpha\|\mathcal{A}\|_{\mathrm{F}}^{2}\right\\}\leq\\\
2\exp\left(-\min\Bigg{\\{}\frac{\alpha^{2}\|\mathcal{A}\|_{\mathrm{F}}^{2}}{144K^{4}\|\mathcal{A}\|_{\mathrm{op}}^{2}},\frac{\alpha\|\mathcal{A}\|_{\mathrm{F}}^{2}}{16\sqrt{2}K^{2}\|\mathcal{A}\|_{\mathrm{op}}^{2}}\Bigg{\\}}\right),$
(15)

where the extension $\mathcal{A}(z)$ to $z$ is defined almost surely and
$\|\cdot\|$ denotes the inner product norm.

The constants that appear in the statement are similar to [17]; we extend
their proof to semi-infinite matrices. The Hanson-Wright inequality can be
used to prove concentration of quadratic forms of the form
$\xi^{*}V^{*}_{k}V_{k}\xi$. This, in turn, implies concentration of the norm
$\|V_{k}\|_{\mathrm{op}}$ via the variational representation of the norm
$\|V_{k}\|_{\mathrm{op}}^{2}=\sup_{\|\xi\|=1}\xi^{*}V^{*}_{k}V_{k}\xi$.

##  4 Estimation over all frequencies

In the previous section, we derived finite-sample guarantees for estimating
the frequency responses at fixed selected frequencies $2\pi\ell/M$,
$\ell\in[M]$. Here, we derive guarantees for estimating the function
$G(e^{j\omega})$ uniformly over all $\omega\in[0,2\pi)$ in the
$\mathcal{H}_{\infty}$ norm. We consider a naive estimator where to compute
$\hat{G}(e^{j\omega})$ we use the closest frequency
$\hat{G}(e^{j2\pi\ell/M})$, for some $\ell\in[M]$.

Let $\chi(t)$ be the indicator function of the half-open interval
$[-1/2,1/2)$. Define the naive estimator

$\hat{G}^{N}(e^{j\omega})\triangleq\sum_{\ell=0}^{M-1}\hat{G}_{\ell
N_{p}}\chi\left(\frac{\omega-\omega_{\ell N_{p}}}{2\pi/M}\right).$ (16)

Due to strict stability, it follows that the frequency response is smooth with
Lipschitz constant upper bounded by $\|G\|_{\star}$. This follows from the
fact that

$\left\|\frac{\partial
G(e^{j\omega})}{\partial\omega}\right\|_{\mathrm{op}}=\left\|\sum_{t=0}^{\infty}jtg_{t}e^{-j\omega
t}\right\|_{\mathrm{op}}\leq\|G\|_{\star}.$

Hence, for any $\omega\in[0,2\pi)$, the error is bounded by

$\|G(e^{j\omega})-\hat{G}^{N}(e^{j\omega})\|_{\mathrm{op}}\leq\frac{\pi}{M}\|G\|_{\star}+\|G(e^{j\omega_{\ell
N_{p}}})-\hat{G}_{\ell N_{p}}\|_{\mathrm{op}},$

for some $\ell\in[M]$. The first term scales with $1/M$, while the latter
scales with $\sqrt{M/N}$, excluding logarithmic terms. Balancing the two
terms, we obtain the following guarantees.

######  Theorem 3 (Guarantees in $\mathcal{H}_{\infty}$ norm).

Let $M=c_{1}N^{1/3}$, for some $c_{1}>0$ such that $M$, $N/M$ are integers.
Consider the naive estimator (16) and fix a failure probability $\delta$.
Under Assumptions 1, 2, 3, 4, with probability at least $1-\delta$:

$\displaystyle\|G-\hat{G}^{N}\|_{\mathcal{H}_{\infty}}\leq\pi
c^{-1}_{1}N^{-1/3}\|G\|_{\star}+c^{1/2}_{1}N^{-5/6}\frac{2\|G\|_{\star}D_{u}}{\underline{\sigma}_{u}}$
$\displaystyle+c^{1/2}_{1}cN^{-1/3}\underline{\mathrm{SNR}}^{-1}_{N}\left(\sqrt{d_{{y}}}+\frac{K^{2}}{\sigma_{{e}}^{2}}\sqrt{d_{{u}}+\log
N/\delta}\right)$ (17)

where $c$ is a universal constant,

$\underline{\mathrm{SNR}}_{N}=\min_{\ell\in[M]}\mathrm{SNR}_{\ell N_{p},N}$

is the worst case SNR, and
$\underline{\sigma}_{u}=\min_{\ell\in[M]}\sigma_{u,\ell}$ is the worst case
excitation among the frequencies of interest.

We can tune the constant $c_{1}$ to guarantee $N_{p}$ is an integer and trade
between the Lipschitz constant and the SNR. Excluding logarithmic factors, we
obtain a rate of $N_{\mathrm{tot}}^{-1/3}$ which is the optimal one for non-
parametric estimation of Lipschitz functions [18, 5]. The rate is suboptimal
when higher order derivatives exist, which would imply stricter stability
conditions than Assumption 4. For example, for twice differentiable functions,
the rate can be improved to $N_{\mathrm{tot}}^{-2/5}$ [1]. In the setting of
rational functions, strict stability is equivalent to exponential stability,
which, in turn, implies the existence of all high-order derivatives of $G$. In
this case, it would be suboptimal to employ the naive estimator (16) without
additional smoothing. We leave this for future work.

The $\mathcal{H}_{\infty}$ norm bound picks up the frequency which is the
hardest to learn as it depends on the worst-case SNR. Assume that the input
excites uniformly all frequencies, that is, $\sigma_{u,\ell}=\sigma_{u}$, for
all $\ell\in[M]$. Then, the worst case SNR scales inversely with
$\sup_{\ell}\sqrt{\|\Phi_{v,N}(\ell N_{p})\|_{\mathrm{op}}}$; as we increase
the frequency resolution $M$, this quantity scales, in turn, with the
$\mathcal{H}_{\infty}$ norm of the noise filter $H(\mathrm{\mathstrut q})$.

Figure 1: The (normalized) empirical maximum error of the ETFE over the fixed
frequency grid $2\pi\ell/M$, $\ell\in[M]$, for fixed $M$. The shaded areas
show one (empirical) standard deviation. It decays with a rate of $N^{-1/2}$.
Moreover, the error increases as we require more resolution, i.e., larger $M$.
The error for $M=2047$ is roughly $\sqrt{2}$ times larger than for $M=1023$,
verifying the result of Theorem 1. Figure 2: The (normalized) empirical
$\mathcal{H}_{\infty}$ norm of the ETFE, based on the naive estimator (16).
The shaded areas show one (empirical) standard deviation. We optimize the
value of $M$ for every choice of $N$. Unlike the error at fixed frequencies,
it decays slower, with a rate of $N^{-1/3}$.

##  5 Simulations

We study the performance of the ETFE and the naive estimator via a numerical
example. Consider the system

$G(\mathrm{\mathstrut q})=\frac{0.12\mathrm{\mathstrut
q}^{-1}+0.18\mathrm{\mathstrut q}^{-2}}{1-1.4\mathrm{\mathstrut
q}^{-1}+1.443\mathrm{\mathstrut q}^{-2}-1.123\mathrm{\mathstrut
q}^{-3}+0.7729q^{-4}}$

with noise filter $H(\mathrm{\mathstrut q})=(1-0.2\mathrm{\mathstrut
q}^{-1})^{-1}$. Let all past inputs be zero $u_{t}=0,t<0$. We generate the
excitation signal based on PRBS [1] with an additional offset to excite the
zero frequency. Note that PRBS maximal length signals require $M=2^{d}-1$, for
some $d\in\mathbb{N}$.

In the first simulation, we study the maximum error of the ETFE over the fixed
grid $\\{2\pi\ell/M\\}_{\ell\in[M]}$

$\|G-\hat{G}^{N}\|_{\widehat{\mathcal{H}_{\infty}}}\triangleq\sup_{\ell\in[M]}\|G(e^{j2\pi\frac{\ell}{M}})-\hat{G}_{\ell
N_{p}}\|_{\mathrm{op}},$

which can be thought as the “discretized” $\mathcal{H}_{\infty}$ norm of the
error. We keep $M$ fixed, $\sigma_{{e}}^{2}=0.1$, and we vary the total number
of samples $N_{\mathrm{tot}}=N$. To visualize the results, we perform $100$
Monte Carlo iterations for every number of samples $N$ and we present the
empirical mean along with one empirical standard deviation. As shown in Fig.
1, the error decays with a rate of $N^{-1/2}$ validating Theorem 1.

In the second simulation, we study the maximum error of the naive estimator
(16) over all frequencies, that is, the “true” $\mathcal{H}_{\infty}$ norm of
the error. For every number of samples $N$, we tune $d$ in $M=2^{d}-1$ to be
the optimal one (empirically), and we perform $100$ Monte Carlo simulations.
As shown in Fig. 2, the error decays with a rate of $N^{-1/3}$ reflecting the
result of Theorem 3.

##  6 Conclusion and Future Work

We provide finite-sample guarantees for the ETFE over a selected frequency
grid, in the case of open-loop periodic excitation and under strict stability
assumptions. By tuning the frequency resolution and exploiting Lipschitz
continuity, we also obtain estimation guarantees in the $\mathcal{H}_{\infty}$
norm. An interesting direction for future work is studying finite-sample non-
parametric least squares [5, 16] in the frequency domain. This approach could
lead to interesting connections between function class complexity and
experiment design. Moreover, adding more structure, beyond Lipschitz
continuity, will lead to faster rates. Other topics that are left for future
work include the estimation of the noise statistics, the prediction problem,
and obtaining finite-sample guarantees when we use random inputs instead of
deterministic ones. Finally, extending minimax lower bounds to this setting is
also open [26].

## Acknowledgment

This work has been supported by the Swiss National Science Foundation under
NCCR Automation (grant agreement 51NF40 180545), and by the European Research
Council under the ERC Advanced grant agreement 787845 (OCAL).

## References

  * [1] L. Ljung, _System Identification: Theory for the User_. Prentice Hall, 1999.
  * [2] J. Schoukens, R. Pintelon, and Y. Rolain, “Time domain identification, frequency domain identification. Equivalencies! Differences?” in _Proceedings of the 2004 American Control Conference_ , vol. 1. IEEE, 2004, pp. 661–666.
  * [3] R. Pintelon and J. Schoukens, _System identification: a frequency domain approach_. John Wiley & Sons, 2012.
  * [4] R. Vershynin, _High-dimensional probability: An introduction with applications in data science_. Cambridge university press, 2018, vol. 47.
  * [5] M. J. Wainwright, _High-dimensional statistics: A non-asymptotic viewpoint_. Cambridge university press, 2019, vol. 48.
  * [6] A. Tsiamis, I. Ziemann, N. Matni, and G. J. Pappas, “Statistical learning theory for control: A finite-sample perspective,” _IEEE Control Systems Magazine_ , vol. 43, no. 6, pp. 67–97, 2023.
  * [7] A. Goldenshluger, “Nonparametric estimation of transfer functions: rates of convergence and adaptation,” _IEEE Transactions on Information Theory_ , vol. 44, no. 2, pp. 644–658, 1998.
  * [8] M. K. S. Faradonbeh, A. Tewari, and G. Michailidis, “Finite Time Identification in Unstable Linear Systems,” _Automatica_ , vol. 96, pp. 342–353, 2018.
  * [9] M. Simchowitz, H. Mania, S. Tu, M. I. Jordan, and B. Recht, “Learning Without Mixing: Towards A Sharp Analysis of Linear System Identification,” _arXiv preprint arXiv:1802.08334_ , 2018.
  * [10] S. Oymak and N. Ozay, “Revisiting Ho-Kalman based system identification: robustness and finite-sample analysis,” _IEEE Transactions on Automatic Control_ , 2021.
  * [11] T. Sarkar, A. Rakhlin, and M. A. Dahleh, “Finite time LTI system identification,” _Journal of Machine Learning Research_ , vol. 22, no. 26, pp. 1–61, 2021.
  * [12] A. Tsiamis and G. J. Pappas, “Finite Sample Analysis of Stochastic System Identification,” in _IEEE 58th Conference on Decision and Control (CDC)_ , 2019.
  * [13] A. Wagenmaker and K. Jamieson, “Active learning for identification of linear dynamical systems,” in _Conference on Learning Theory_. PMLR, 2020, pp. 3487–3582.
  * [14] S. Tu, R. Frostig, and M. Soltanolkotabi, “Learning from many trajectories,” _arXiv preprint arXiv:2203.17193_ , 2022.
  * [15] I. Ziemann and S. Tu, “Learning with little mixing,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 4626–4637, 2022.
  * [16] I. M. Ziemann, H. Sandberg, and N. Matni, “Single trajectory nonparametric learning of nonlinear dynamics,” in _conference on Learning Theory_. PMLR, 2022, pp. 3333–3364.
  * [17] I. Ziemann, A. Tsiamis, B. Lee, Y. Jedra, N. Matni, and G. J. Pappas, “A Tutorial on the Non-Asymptotic Theory of System Identification,” in _2023 62nd IEEE Conference on Decision and Control (CDC)_. IEEE, 2023, pp. 8921–8939.
  * [18] A. Tsybakov, _Introduction to Nonparametric Estimation_ , ser. Springer Series in Statistics. Springer New York, 2008.
  * [19] D. L. Hanson and F. T. Wright, “A bound on tail probabilities for quadratic forms in independent random variables,” _The Annals of Mathematical Statistics_ , vol. 42, no. 3, pp. 1079–1083, 1971.
  * [20] A. Carè, R. Carli, A. Dalla Libera, D. Romeres, and G. Pillonetto, “Kernel methods and gaussian processes for system identification and control: A road map on regularized kernel-based learning for control,” _IEEE Control Systems Magazine_ , vol. 43, no. 5, pp. 69–110, 2023.
  * [21] A. Devonport, P. Seiler, and M. Arcak, “Frequency domain gaussian process models for $\mathcal{H}_{\infty}$ uncertainties,” in _Learning for Dynamics and Control Conference_. PMLR, 2023, pp. 1046–1057.
  * [22] G. Baggio, A. Carè, A. Scampicchio, and G. Pillonetto, “Bayesian frequentist bounds for machine learning and system identification,” _Automatica_ , vol. 146, p. 110599, 2022.
  * [23] T. P. Dobrowiecki, J. Schoukens, and P. Guillaume, “Optimized excitation signals for MIMO frequency response function measurements,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 55, no. 6, pp. 2072–2079, 2006.
  * [24] M. Gevers, X. Bombois, R. Hildebrand, and G. Solari, “Optimal experiment design for open and closed-loop system identification,” _Communications in Information and Systems_ , vol. 11, no. 3, pp. 197–224, 2011.
  * [25] M. Rudelson and R. Vershynin, “Hanson-Wright inequality and sub-gaussian concentration,” _Electronic Communications in Probability_ , vol. 18, 2013.
  * [26] S. Tu, R. Boczar, and B. Recht, “Minimax lower bounds for $\mathcal{H}_{\infty}$-norm estimation,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 3538–3543.
  * [27] R. Durrett, _Probability: theory and examples_. Cambridge University Press, 2010.

##  Appendix A Proof of Theorem 2

Recall that $\mathcal{L}_{2}(p)$ represents the Hilbert space of
$p$-dimensional real-valued square summable sequences. For any sequences
$\alpha,\beta\in\mathcal{L}_{2}(p)$, their inner product is defined as

$\langle\alpha,\beta\rangle\triangleq\sum_{k=1}^{\infty}\alpha^{\top}_{k}\beta_{k},$

with the respective inner product norm

$\|\alpha\|=\sqrt{\langle\alpha,\alpha\rangle}.$

Let $\mathcal{A}:\mathcal{L}_{2}(p)\rightarrow^{d}$ be a linear map from
sequences to finite-dimensional vectors. Consider the standard orthonormal
basis for $\mathcal{L}_{2}(p)$, that is, the set $\\{b_{i}\\}_{i=1}^{\infty}$,
where $b_{i}=(0,0,\dots,1,0,\dots)$ is the unit norm sequence with all
elements $0$ except for the $i$-th one. Then, we can represent $\mathcal{A}$
with a semi-infinite block matrix

$\begin{bmatrix}\mathcal{A}_{11}&\mathcal{A}_{12}&\cdots\\\
\mathcal{A}_{21}&\mathcal{A}_{22}&\cdots\\\ \vdots\\\
\mathcal{A}_{d1}&\mathcal{A}_{d2}&\cdots\end{bmatrix},$

where $\mathcal{A}_{ij}\in^{1\times p}$, $i=1,\dots,d$, $j\geq 1$. Then, the
Frobenius norm is given by

$\|\mathcal{A}\|_{\mathrm{F}}^{2}\triangleq\sum_{i=1}^{d}\sum_{k=0}^{\infty}\left\|\mathcal{A}_{ij}\right\|^{2}.$

In the following, we will study operators with bounded Frobenius norm
$\|\mathcal{A}\|_{\mathrm{F}}<\infty$; the operator norm is upper bounded by
the Frobenius norm, and, thus, also bounded.

Consider now a sequence $z\triangleq\\{z_{t}\\}_{t=1}^{\infty}$ of independent
zero mean $K-$sub-Gaussian random variables, with $z_{t}\in^{p}$. We will
study quadratic forms $\|\mathcal{A}(z)\|^{2}$. Since $z$ is not square
summable, $\mathcal{A}(z)$ is interpreted in terms of the extension

$[\mathcal{A}(z)]_{i}=\lim_{k\rightarrow\infty}\sum_{j=1}^{k}{\mathcal{A}_{ij}z_{j}},\,i=1,\dots,d$
(18)

when the above limit exists. We note that in our setting $\mathcal{A}(z)$ is
well-defined.

######  Lemma 5 (Well-posedness).

Let $\mathcal{A}:\mathcal{L}_{2}(p)\rightarrow^{d}$ be a linear map with
bounded Frobenius norm. Let $z=\\{z_{t}\\}_{t=1}^{\infty}$ be zero-mean,
independent, and $K-$sub-Gaussian. Then, the extension $\mathcal{A}(z)\in^{d}$
is well-defined almost surely.

###### Proof.

Since $z_{t}$, $t\geq 1$ are $K^{2}-$sub-Gaussian, we have

$\sup_{t\geq 1}\mathbb{E}z^{\top}_{t}z_{t}=cpK^{2}<\infty,$

for some universal constant $c$. This follows from the fact that every
coordinate $[z_{t}]_{s}$, $s=1,\dots,p$ is also $K^{2}-$sub-gaussian and the
bound $\mathbb{E}([z_{t}]_{s})^{2}\leq cK^{2}$ (see [4, Prop. 2.5.2]). Fix an
$i$, $i=1,\dots,d$. We have

$\sum_{j=1}^{k}\mathrm{var}(\mathcal{A}_{ij}z_{j})=\sum_{j=1}^{k}\mathbb{E}z_{j}^{\top}\mathcal{A}^{\top}_{ij}\mathcal{A}_{ij}z_{j}=\sum_{j=1}^{k}\|\mathcal{A}_{ij}\|^{2}\mathbb{E}\|z_{j}\|^{2}\leq\sum_{j=1}^{k}\|\mathcal{A}_{ij}\|^{2}cpK^{2}$

Hence, we have

$\sum_{j=1}^{\infty}\mathrm{var}(\mathcal{A}_{ij}z_{j})\leq\|\mathcal{A}\|_{\mathrm{F}}^{2}cpK^{2}<\infty.$

Hence, by the Kolmogorov maximal inequality [27, Th. 2.5.3], the series
$\lim_{k\rightarrow\infty}\sum_{j=1}^{k}{\mathcal{A}_{ij}z_{j}}$ converges
almost surely for any $i=1,\dots,d$ and, thus, $\mathcal{A}(z)$ is almost
surely well-defined. ∎

The main idea behind proving Theorem 2 is using the Chernoff bound method. Let
$X=\|\mathcal{A}(z)\|^{2}-\mathbb{E}\|\mathcal{A}(z)\|^{2}$. Then, for any
$\lambda\geq 0$

$\mathbb{P}(X\geq t)\leq\min_{\lambda\geq 0}e^{-\lambda t}\mathbb{E}e^{\lambda
X}.$

To control the right-hand side, we have to upper bound the moment-generating
function $\mathbb{E}e^{\lambda X}$, which requires most of the work.

######  Lemma 6 (MGF bound).

Let $z=\\{z_{t}\\}_{t=1}^{\infty}$ be a sequence of independent, zero-mean,
$K^{2}$-sub-Gaussian random variables. Let
$\mathcal{A}:\mathcal{L}_{2}(p)\rightarrow^{d}$ be a linear map with bounded
Frobenius norm $\|\mathcal{A}\|_{\mathrm{F}}<\infty$. For every $\lambda\in$
such that
$\left|\lambda\right|\leq(8\sqrt{2}K^{2}\|\mathcal{A}\|_{\mathrm{op}}^{2})^{-1}$,
we have

$\mathbb{E}\exp(\lambda\|\mathcal{A}(z)\|^{2}-\lambda\mathbb{E}\|\mathcal{A}(z)\|^{2})\leq
36\lambda^{2}K^{4}\|\mathcal{A}\|_{\mathrm{op}}^{2}\|\mathcal{A}\|_{\mathrm{F}}^{2},$

where the extension $\mathcal{A}(z)$ is interpreted as in (18).

###### Proof.

We will consider finite-dimensional truncations of the quadratic form
$\|\mathcal{A}(z)\|^{2}_{2}$, apply the result of [17], and take limits. Let
$\hat{\mathcal{A}}^{(k)}:\mathcal{L}_{2}(p)\rightarrow^{d}$ be the truncated
map with semi-infinite matrix representation

$\hat{\mathcal{A}}^{(k)}_{ij}=\left.\begin{aligned}
&\mathcal{A}_{ij},\,\text{for }i=1,\dots,d,\,j\leq k\\\ &0\text{
elsewhere}\end{aligned}\right\\}.$ (19)

Define the random variables

$\displaystyle X_{k}$
$\displaystyle\triangleq\exp(\lambda\|\hat{\mathcal{A}}^{(k)}(z)\|^{2}-\lambda\mathbb{E}\|\hat{\mathcal{A}}^{(k)}(z)\|^{2}),\,k\geq
1$ $\displaystyle X$
$\displaystyle\triangleq\exp(\lambda\|\mathcal{A}(z)\|^{2}-\lambda\mathbb{E}\|\mathcal{A}(z)\|^{2}).$

Step a): truncation. First, we will establish the desired inequality for all
$X_{k}$. Notice that $\|\hat{\mathcal{A}}^{(k)}(z)\|^{2}_{2}$ is a finite
dimensional quadratic form. Hence, Proposition A.1 of [17] applies and
provides the bound

$\mathbb{E}X_{k}\leq
36\lambda^{2}K^{4}\|(\hat{\mathcal{A}}^{(k)})^{*}\hat{\mathcal{A}}^{(k)}\|_{\mathrm{F}}^{2},$

for
$\left|\lambda\right|\leq(8\sqrt{2}K^{2}\|(\hat{\mathcal{A}}^{(k)})^{*}\hat{\mathcal{A}}^{(k)}\|_{\mathrm{op}})$
where $(\hat{\mathcal{A}}^{(k)})^{*}$ denotes the adjoint. Truncation does not
increase the norm, i.e.,
$\|\mathcal{A}\|_{\mathrm{F}}\geq\|\hat{\mathcal{A}}^{(k)}\|_{\mathrm{F}}$ and
$\|\mathcal{A}\|_{\mathrm{op}}\geq\|\hat{\mathcal{A}}^{(k)}\|_{\mathrm{op}}$.
Moreover

$\|(\hat{\mathcal{A}}^{(k)})^{*}\hat{\mathcal{A}}^{(k)}\|_{\mathrm{F}}^{2}\leq\|\hat{\mathcal{A}}^{(k)}\|_{\mathrm{op}}^{2}\|\hat{\mathcal{A}}^{(k)}\|_{\mathrm{F}}^{2}.$

Hence, we can write

$\mathbb{E}X_{k}\leq
36\lambda^{2}K^{4}\|\mathcal{A}\|_{\mathrm{op}}^{2}\|\mathcal{A}\|_{\mathrm{F}}^{2},$
(20)

for
$\left|\lambda\right|\leq(8\sqrt{2}K^{2}\|\mathcal{A}\|_{\mathrm{op}}^{2})$.

Step b): convergence. We prove that $X_{k}$ converges almost surely to $X$.
Convergence of $\hat{\mathcal{A}}^{(k)}(z)$ to $\mathcal{A}(z)$ follows from
Lemma 5. To show convergence of $\mathbb{E}\|\hat{\mathcal{A}}^{(k)}(z)\|^{2}$
observe that

$|\|\mathcal{A}(z)\|^{2}-\|\hat{\mathcal{A}}^{(k)}(z)\|_{2}^{2}|\leq
2\left|\left(\hat{\mathcal{A}}^{(k)}(z)\right)^{\top}(\mathcal{A}-\hat{\mathcal{A}}^{k})(z)\right|+\|(\mathcal{A}-\hat{\mathcal{A}}^{k})(z)\|^{2}.$

Let $C_{z}=\sup_{t\geq 1}\mathbb{E}z^{\top}_{t}z_{t}$; due to sub-Gaussianity
$C_{z}=cpK^{2}$ for some universal constant $c$–see proof of Lemma 5. Invoking
the Cauchy-Schwartz inequality

$\mathbb{E}|\|\mathcal{A}(z)\|^{2}-\|\hat{\mathcal{A}}^{(k)}(z)\|_{2}^{2}|\leq
2C_{z}\|\mathcal{A}\|_{\mathrm{F}}\|\mathcal{A}-\hat{A}^{(k)}\|_{\mathrm{F}}+C_{z}\|\mathcal{A}-\hat{A}^{(k)}\|_{\mathrm{F}}^{2}\rightarrow
0,$

as $k\rightarrow\infty$. Convergence of $X_{k}$ to $X$ follows by continuity.

Step c): limit. Since $X_{k}>0$, by Fatou’s lemma

$\mathbb{E}X\leq\liminf\mathbb{E}X_{k}.$ (21)

The result follows from (20), (21). ∎

The proof of Theorem 2 now follows from Lemma 6 and the Chernoff bound method.
The remaining details are omitted since they are identical to [17, proof of
Th. 2.1 ]. ∎

##  Appendix B Proof of Lemma 4

###### Proof.

Due to the independence between the experiments

$\mathbb{E}V^{*}_{k}V_{k}=\mathrm{diag}(\mathbb{E}(V^{(i)})^{*}_{k}V^{(i)}_{k})=\mathrm{tr}(\Phi_{v,N}(k))I_{d_{{u}}}.$

By the variational representation of the operator norm for symmetric matrices

$\displaystyle\|V^{*}_{k}V_{k}-\mathrm{tr}(\Phi_{v,N}(k))I_{d_{{u}}}\|_{\mathrm{op}}$
$\displaystyle=\sup_{\xi\in\mathbb{C}^{d_{{u}}},\|\xi\|=1}\left|\xi^{*}(V^{*}_{k}V_{k}-\mathrm{tr}(\Phi_{v,N}(k))I_{d_{{u}}})\xi\right|$
(22)
$\displaystyle=\sup_{\xi\in\mathbb{C}^{d_{{u}}},\|\xi\|=1}\left|\|V_{k}\xi\|^{2}-\mathbb{E}\|V_{k}\xi\|^{2}\right|.$

We need to control the quantity inside the supremum over the whole unit
sphere. First, we will bound it for a fixed $\xi$ using the Hanson-Wright
inequality. Then, we will discretize the unit sphere and use a covering
argument.

Step a: fixed $\xi$. Observe that $V^{(i)}_{k}=\mathcal{F}_{k}((h*e^{(i)}))$
is the output of a linear map applied to the infinite sequence
$e^{(i)}\triangleq\\{e^{(i)}_{N-1},e^{(i)}_{N-2},\dots\\}$, where $*$ denotes
convolution. We can rewrite the above linear relation using the notation
$V^{(i)}_{k}=H_{k}(e^{(i)})$, where $H_{k}$ represents the map
$\mathcal{F}_{k}\circ((h*\cdot))$. For convenience, let us group all noise
sequences $e^{(i)}$ into one single sequence $e$. We use the enumeration

$\displaystyle e$
$\displaystyle\triangleq\left\\{e_{t}\right\\}_{t=1}^{\infty}\triangleq\left\\{e^{(1)}_{N-1},\dots,e^{(d_{{u}})}_{N-1},e^{(1)}_{N-2},\dots,e^{(d_{{u}})}_{N-2},\dots\right\\}.$

Let $\xi\in\mathbb{C}^{d_{{u}}}$ be any vector with unit norm $\|\xi\|=1$.
Then,

$V_{k}\xi=H_{k}(\sum_{i=1}^{d_{{u}}}e^{(i)}\xi_{i})=(H_{k}\circ\Xi)(e),$

where $\Xi$ is the linear map from $e$ to
$\sum_{i=1}^{d_{{u}}}e^{(i)}\xi_{i}$. We can represent $\Xi$ with the infinite
matrix

$\Xi=\begin{bmatrix}\xi^{*}&0&0&\cdots\\\ 0&\xi^{*}&0&\cdots\\\ 0&0&\xi^{*}\\\
\vdots&&&\ddots\end{bmatrix}.$

Since $\xi^{*}\xi=1$, $\Xi\Xi^{*}$ is the identity operator and has unit norm.
Since $V_{k}\xi\in\mathbb{C}^{d_{{y}}}$ is complex-valued, we need to embed
the vector to a real vector space to apply the Hanson-Wright inequality. Let
$\Re(H_{k}\circ\Xi)$ denote the real part of the operator $H_{k}\circ\Xi$ with
$\Im(H_{k}\circ\Xi)$ the respective imaginary part. Then, lift the space to
${}^{2d_{{y}}}$

$\hat{V}_{k}\xi=\mathcal{A}(e)\triangleq\begin{bmatrix}\Re(H_{k}\circ\Xi)(e)\\\
\Im(H_{k}\circ\Xi)(e)\end{bmatrix}\in^{2d_{{y}}}.$

Since the sequence $e$ is real-valued

$\xi^{*}\hat{V}^{*}_{k}\hat{V}_{k}\xi=\xi^{*}V^{*}_{k}V_{k}\xi.$

Hence, we have

$\sigma_{{e}}^{2}\|\mathcal{A}\|_{\mathrm{F}}^{2}=\mathbb{E}\|\mathcal{A}(e)\|_{2}^{2}=\mathbb{E}\xi^{\top}V^{*}_{k}V_{k}\xi=\mathrm{tr}(\Phi_{v,N}(k)).$

For the operator norm, we use the bound
$\|\mathcal{A}\|_{\mathrm{op}}\leq\|H_{k}\|_{\mathrm{op}}=\sqrt{\|\Phi_{v,N}(k)\|_{\mathrm{op}}}/\sigma_{{e}}$,
where the first inequality follows from the fact that $\Xi$ has unit operator
norm. The equality follows from the fact that
$\Phi_{v,N}(k)=\mathbb{E}(V^{(i)}_{k}(V^{(i)}_{k})^{*})=\sigma_{{e}}^{2}H_{k}H^{*}_{k}$,
where $H^{*}_{k}$ denotes the adjoint of $H_{k}$. Applying the Hanson-Wright
inequality to $\mathcal{A}(e)$, we obtain

$\mathbb{P}\left\\{\big{|}\|V_{k}\xi\|^{2}-\mathrm{tr}(\Phi_{v,N}(k))\big{|}>\beta\>\mathrm{tr}(\Phi_{v,N}(k))\right\\}\leq
2\exp\left(-c\min\\{\beta^{2},\beta\\}\frac{\sigma_{{e}}^{4}}{K^{4}}\frac{\mathrm{tr}(\Phi_{v,N}(k))}{\|\Phi_{v,N}(k)\|_{\mathrm{op}}}\right),$
(23)

where we select $\alpha=\beta\mathrm{tr}(\Phi_{v,N}(k))$. We simplify the
expression by choosing $c=1/144$ and exploiting the fact that the sub-Gaussian
parameter upper-bounds the variance that $K\geq\sigma_{{e}}$.

Step b: covering argument. To control the supremum in (22), we need to
discretize the unit sphere. In particular, we consider an $\epsilon-$net and
apply the following result.

######  Lemma 7 (Operator Norm on a net [4, 17]).

Let $A\in\mathbb{C}^{d\times d}$ be a hermitian random matrix and let
$\epsilon\in(0,1/2)$. Let $\mathcal{N}$ be an $\epsilon-$net of the complex
unit sphere with minimal cardinality. Then, for any $t>0$

$\mathbb{P}(\|A\|_{\mathrm{op}}\geq
t)\leq\left|\mathcal{N}\right|\max_{z\in\mathcal{N}}\mathbb{P}(\left|z^{*}Az\right|\geq(1-2\epsilon)t),$

where the cardinality $\left|\mathcal{N}\right|$ is upper bounded

$\left|\mathcal{N}\right|\leq\left(1+\frac{2}{\epsilon}\right)^{2d}.$

Lemma 7 follows directly by Lemma 2.5 in [17]. The exponent scales with $2d$
instead of $d$ since the complex unit sphere in $\mathbb{C}^{d}$ is isomorphic
to the unit sphere $S^{2d-1}$ in 2d.

Applying Lemma 7 with $\epsilon=1/4$, $t$ replaced with $2\beta$, and $A$
replaced with $V^{*}_{k}V_{k}-\mathrm{tr}(\Phi_{v,N}(k))I_{d_{{u}}}$, we
obtain

$\mathbb{P}\left\\{\|V^{*}_{k}V_{k}-\mathrm{tr}(\Phi_{v,N}(k))I_{d_{{u}}}\|_{\mathrm{op}}>2\beta\>\mathrm{tr}(\Phi_{v,N}(k))\right\\}\leq\\\
9^{2d_{{u}}}2\exp\left(-c\min\left\\{\beta^{2},\beta\right\\}\frac{\sigma_{{e}}^{4}}{K^{4}}\frac{\mathrm{tr}(\Phi_{v,N}(k))}{\|\Phi_{v,N}(k)\|_{\mathrm{op}}}\right)$
(24)

We only need to simplify the final expression. We invoke the following lemma.

######  Lemma 8 ((3.2) in [4]).

Let $z,s\geq 0$ be any two positive real numbers. If $\left|z-1\right|\geq s$
then $\left|z^{2}-1\right|\geq\max\\{s,s^{2}\\}$.

Let $s=\min\\{\beta,\beta^{2}\\}$, which implies that
$\beta=\max\\{s,s^{2}\\}$. Then, the event

$\mathcal{E}_{1}\triangleq\left\\{\Big{|}\|V_{k}\|_{\mathrm{op}}/\sqrt{\mathrm{tr}(\Phi_{v,N}(k))}-1\Big{|}\geq
2s\right\\}$

implies that there exists a unit norm $\xi\in^{d_{{u}}}$ such that

$\Big{|}\|V_{k}\xi\|/\sqrt{\mathrm{tr}(\Phi_{v,N}(k))}-1\Big{|}\geq 2s.$

By Lemma 8, this implies that

$\Big{|}\|V_{k}\xi\|^{2}/\mathrm{tr}(\Phi_{v,N}(k))-1\Big{|}\\!\geq\\!\max\\{2s,4s^{2}\\}\geq
2\max\\{s,s^{2}\\}.$

This together with (22) imply the following event holds

$\mathcal{E}_{2}\triangleq\left\\{\|V^{*}_{k}V_{k}-\mathrm{tr}(\Phi_{v,N}(k))I_{d_{{u}}}\|_{\mathrm{op}}>2\beta\>\mathrm{tr}(\Phi_{v,N}(k))\right\\}.$

Since $\mathcal{E}_{1}\subset\mathcal{E}_{2}$ the result follows from (24). ∎

##  Appendix C Proof of Lemma 3

Fix an experiment index $i$. Since $u^{(i)}_{t}$ is periodic with period $M$
and $N_{p}$ is an integer, we have

$\displaystyle\sum_{t=0}^{N-1}u^{(i)}_{t}e^{-j\omega_{k}t}$
$\displaystyle=\sum_{\tau=0}^{N_{p}-1}\sum_{t=0}^{M-1}u^{(i)}_{t}e^{-j\omega_{k}(M\tau+t)}$
$\displaystyle=\left(\sum_{\tau=0}^{N_{p}-1}e^{-j\frac{2\pi
k}{N_{p}}\tau}\right)\left(\sum_{t=0}^{M-1}u^{(i)}_{t}e^{-j\omega_{k}t}\right)$

Therefore, we obtain

$\displaystyle U^{(i)}_{N_{p}\ell}$
$\displaystyle=\sqrt{N_{p}}\tilde{U}^{(i)}_{l},\,\text{for
}N_{p}\ell,\ell\in[M]$ $\displaystyle U^{(i)}_{k}$
$\displaystyle=0,\,\text{elsewhere. }$

For $k=N_{p}\ell$ we have

$U_{k}U^{*}_{k}=N_{p}\tilde{U}_{\ell}\tilde{U}^{*}_{\ell}\succeq
N_{p}\sigma^{2}_{u,\ell}I_{d_{{u}}}.$

##  Appendix D Proof of Theorem 1

By (10), the estimation error is upper bounded by

$\|\hat{G}_{k}-G(e^{j\omega_{k}})\|_{\mathrm{op}}\leq\|V_{k}\|_{\mathrm{op}}\|U_{k}^{-1}\|_{\mathrm{op}}+\|T_{k,N}\|_{\mathrm{op}}\|U^{-1}_{k}\|_{\mathrm{op}}.$

Invoking Lemma 4, we have with probability at least $1-\tilde{\delta}$

$\displaystyle\|V_{k}\|_{\mathrm{op}}$
$\displaystyle\leq\sqrt{\mathrm{tr}(\Phi_{v,N}(k))}+2\tilde{s}\sqrt{\mathrm{tr}(\Phi_{v,N}(k))}$
$\displaystyle\tilde{s}^{2}\mathrm{tr}(\Phi_{v,N}(k))$
$\displaystyle=c^{\prime}\frac{K^{4}}{\sigma_{{e}}^{4}}\|\Phi_{v,N}(k))\|_{\mathrm{op}}(\log(1/\tilde{\delta})+d_{{u}}),$

for some universal constant $c^{\prime}$. Set $\delta/M=\tilde{\delta}$ and
consider a union bound over all $\ell\in[M]$. Then, with probability at least
$1-\delta$, for all $k=N_{p}\ell$, $\ell\in[M]$

$\|V_{k}\|_{\mathrm{op}}\leq\sqrt{d_{{y}}}\sqrt{\|\Phi_{v,N}(k)\|_{\mathrm{op}}}+\sqrt{c^{\prime}}\frac{K^{2}}{\sigma_{{e}}^{2}}\sqrt{\|\Phi_{v,N}(k))\|_{\mathrm{op}}}\sqrt{\log(M/\delta)+d_{{u}}},$
(25)

where we used $\mathrm{tr}(\Phi_{v,N}(k))\leq
d_{{y}}\|\Phi_{v,N}(k)\|_{\mathrm{op}}$. The proof now follows from Lemma 2
and Lemma 3. ∎

##  Appendix E Proof of Theorem 3

It follows directly from Theorem 1 and

$\|G(e^{j\omega})-\hat{G}^{N}(e^{j\omega})\|_{\mathrm{op}}\leq\frac{\pi}{M}\|G\|_{\star}+\|G(e^{j\omega_{\ell
N_{p}}})-\hat{G}_{\ell N_{p}}\|_{\mathrm{op}}.$