<figure><img src="image_1.png" /></figure>

HAL
open science

# Minimax Estimation of Partially-Observed Vector AutoRegressions

Guillaume Dalle, Yohann de Castro

► To cite this version:

Guillaume Dalle, Yohann de Castro. Minimax Estimation of Partially-Observed Vector AutoRegressions. 2022. hal-03263275v3

HAL Id: hal-03263275

https://hal.science/hal-03263275v3

Preprint submitted on 5 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Minimax estimation of partially-observed vector autoregressions

Guillaume Dalle¹ and Yohann De Castro²

¹LVMT, ENPC, Institut Polytechnique de Paris, Univ Gustave Eiffel, Marne-la-Vallée, France,
e-mail: guillaume.dalle@enpc.fr

²Univ. Lyon, École Centrale de Lyon, CNRS UMR 5208, Institut Camille Jordan, France,
e-mail: yohann.de-castro@ec-lyon.fr

**Abstract:** Multi-dimensional time series are a core ingredient of the statistical modeling toolkit, for which numerous estimation methods are known. But when observations are scarce or corrupted, the learning task becomes much harder. The question is: how much harder?

In this paper, we study the properties of a partially-observed Vector AutoRegressive process, which is a state-space model endowed with a stochastic observation mechanism. Our goal is to estimate its transition matrix, but we only have access to a small and noisy subsample of the state components.

We start by describing an estimator based on the Yule-Walker equation and the Dantzig selector, and we give an upper bound on its non-asymptotic error. Then, we provide a matching minimax lower bound, thus proving near-optimality of our estimator. The convergence rate we obtain sheds light on the role of several key parameters such as the sampling ratio, the amount of noise and the number of non-zero coefficients in the transition matrix, which can be dense in low dimension but is usually sparse in high dimension. These theoretical findings are commented and illustrated by numerical experiments on simulated data.

**MSC2020 subject classifications:** Primary 62M10, 62F12; secondary 62D10.

**Keywords and phrases:** Vector autoregression, partial observation, sparsity, minimax lower bound, railway modeling.

Received December 2022.

## 1. Introduction

Time series provide a natural representation for periodic measurements of a stochastic process. In particular, those defined by linear Gaussian recursions may be the most widely used and the easiest to study. Well-known examples include the AutoRegressive (AR) process and its multivariate counterpart, the Vector AutoRegressive (VAR) process, two staples of mathematical statistics.

Industrial applications of these models encounter two main challenges. First, they often involve signals in high dimension, which means sparsity assumptions play an important role. Second, the variables of interest are rarely measured exactly or entirely. Indeed, physical constraints such as the cost of sensors can make it impossible to capture every component of the system's state at all times. It is therefore natural to ask: *how much harder does high-dimensional learning become when one only observes a fraction of the relevant values?*
# MINIMAX ESTIMATION OF PARTIALLY-OBSERVED VECTOR AUTOREGRESSIONS

Guillaume Dalle

CERMICS, École des Ponts
Marne-la-Vallée, France
guillaume.dalle@enpc.fr

Yohann De Castro

Institut Camille Jordan, École Centrale Lyon
Écully, France
yohann.de-castro@ec-lyon.fr

April 2022

## Abstract

High-dimensional time series are a core ingredient of the statistical modeling toolkit, for which numerous estimation methods are known. But when observations are scarce or corrupted, the learning task becomes much harder. The question is: how much harder?

In this paper, we study the properties of a partially-observed Vector AutoRegressive process, which is a state-space model endowed with a stochastic observation mechanism. Our goal is to estimate its sparse transition matrix, but we only have access to a small and noisy subsample of the state components. Interestingly, the sampling process itself is random and can exhibit temporal correlations, a feature shared by many realistic data acquisition scenarios.

We start by describing an estimator based on the Yule-Walker equation and the Dantzig selector, and we give an upper bound on its non-asymptotic error. Then, we provide a matching minimax lower bound, thus proving near-optimality of our estimator. The convergence rate we obtain sheds light on the role of several key parameters such as the sampling ratio, the amount of noise and the number of non-zero coefficients in the transition matrix. These theoretical findings are commented and illustrated by numerical experiments on simulated data.

# 1 INTRODUCTION

Time series provide a natural representation for periodic measurements of a stochastic process. In particular, those defined by linear Gaussian recursions may be the most widely used and the easiest to study. Well-known examples include the AutoRegressive (AR) process and its multivariate counterpart, the Vector AutoRegressive (VAR) process.

Industrial applications of these models encounter two main challenges. First, they often involve signals in high dimension, which means sparsity assumptions play an important role. Second, the variables of interest are rarely measured exactly or entirely. Indeed, physical constraints such as the cost of sensors can make it impossible to capture every component of the system's state at all times. It is therefore natural to ask: *how much harder does high-dimensional learning become when one only observes a fraction of the relevant values?*

## 1.1 CONTEXT OF THE STUDY

To answer this question, we study a state-space model where the state $X_t \in \mathbb{R}^D$ follows a VAR process of order 1 over a period of length $T$. Since the dimension $D$ of $X_t$ is high, we assume that its transition matrix $\theta \in \mathbb{R}^{D \times D}$ is $s$-sparse (there are no more than $s$ non-zero coefficients in each row). However, we do not observe the state itself: our observations $Y_t$ only involve the subset of components $X_{t,d}$ for which $\pi_{t,d} = 1$, where $\pi_t$ is a vector of Bernoulli variables. To make matters worse, this subset is corrupted with noise, which leads to the following generative procedure:

$$X_t = \theta X_{t-1} + \mathcal{N}(0, \sigma^2 I) \quad \pi_{t,d} \sim \mathcal{B}(p) \quad Y_t = \text{diag}(\pi_t) X_t + \mathcal{N}(0, \omega^2 I). \quad (1)$$
Among numerous practical applications of VAR models with high societal impact, we can cite econometrics and causality analysis [28]. In such scenarios, sparsity is a key asset for interpreting mutually correlated time series with many components, especially when the data are noisy. The lens of partial observations also offers a direct connection with state-space models [11], heavily used in control of dynamical systems and industrial reliability applications.

## 1.1. Context of the study

We study a state-space model indexed by time $t$, where the state $X_t \in \mathbb{R}^D$ follows a VAR process of order 1 over a period of length $T$. We account for the case where the dimension $D$ of $X_t$ is high, by assuming that its transition matrix $\theta \in \mathbb{R}^{D \times D}$ is $s$-sparse: there are no more than $s$ non-zero coefficients in each row. Note that we recover the dense case by taking $s = D$. The distribution of innovations is an isotropic multivariate Gaussian $N(0, \sigma^2 I)$ with variance $\sigma^2$. However, we do not observe the state itself: our observations $Y_t \in \mathbb{R}^D$ only involve the subset of components $X_{t,d}$ for which $\pi_{t,d} = 1$, where $\pi_t$ is a vector of independent Bernoulli variables $\mathcal{B}(p)$ with the same parameter $p$. To make matters worse, this subset is corrupted by an isotropic Gaussian noise $N(0, \omega^2 I)$ with variance $\omega^2$. This leads to the following generative procedure:

$$X_t = \theta X_{t-1} + N(0, \sigma^2 I) \quad \pi_{t,d} \sim \mathcal{B}(p) \quad Y_t = \text{diag}(\pi_t) X_t + N(0, \omega^2 I). \quad (1)$$

To shed light on the properties of our model, we start by constructing a sparsity-aware estimator for $\theta$, whose non-asymptotic error we upper-bound. We complement this finding with a lower bound on the minimax error that does not depend on the choice of estimator. Upper and lower bound match in most regards, which proves their optimality. A rough summary of our analysis is that the best possible estimator $\hat{\theta}$ satisfies

$$\|\hat{\theta} - \theta\|_{\infty} \lesssim \left(1 + \frac{\omega^2}{\sigma^2}\right) \frac{s}{p\sqrt{T}} \quad (2)$$

with high probability, where $\|\cdot\|_{\infty}$ denotes the operator $\ell_{\infty}$ norm. We observe that the error does not crucially depend on the state dimension $D$ (except for logarithmic terms, see below), but rather on the sparsity $s$ of the transition matrix. As expected, it decreases linearly as $p$ grows, since more information becomes available. Lastly, it is a function of $\omega^2/\sigma^2$, which means that precise recovery of $\theta$ is only possible when the noise $\omega$ is not too much larger than the signal amplitude $\sigma$.

Novel features of our work include the first proof of a minimax lower bound in this setting, the combination of discrete and continuous concentration inequalities to obtain error estimates, as well as detailed numerical experiments on simulated data. We highlight that the analysis covers both sparse and dense transition matrices, bringing new results to both cases.

## 1.2. Example of application

Our study was inspired by concrete questions related to delay propagation on railway networks, which came up during a collaboration with a leading railway company. When external factors (weather, passenger behavior, mechanical failures) trigger a primary delay, resource
When we write $\pi_{t,d} \sim B(p)$, we mean that the marginals of the sampling variables are identical, which requires that every state component be sampled with equal probability $p$. However, we reject the standard independence assumption in favor of temporal dependencies between the Bernoulli variables $\pi_{t,d}$ (see Section 1.2 for a practical justification).

To shed light on the properties of our model, we start by constructing a sparse estimator for $\theta$, whose non-asymptotic error we upper bound. We complement this finding with a lower bound on the minimax error that does not depend on the choice of estimator. Upper and lower bound match in most regards, which proves their optimality. A rough summary of our analysis is that the best possible estimator $\hat{\theta}$ satisfies

$$ \|\hat{\theta} - \theta\|_{\infty} \lesssim \left(1 + \frac{\omega^2}{\sigma^2}\right) \frac{s}{p\sqrt{T}} \qquad (2) $$

with high probability. We observe that the error does not depend on the state dimension $D$, but only on the sparsity $s$ of the transition matrix. As expected, it decreases linearly as $p$ grows, since more information becomes available. Lastly, it is a function of $\omega^2/\sigma^2$, which means that precise recovery of $\theta$ is only possible when the noise is not too much larger than the signal.

Novel features of our work include the first proof of a minimax lower bound in this setting (to the best of our knowledge), the investigation of temporal correlations within the sampling process, the combination of discrete and continuous concentration inequalities to obtain error estimates, as well as detailed numerical experiments on simulated data.

## 1.2 EXAMPLE OF APPLICATION

Our study was inspired by concrete questions related to delay propagation on railway networks, which came up during a collaboration with a leading railway company. When external factors (weather, passenger behavior, mechanical failures) trigger a primary delay, resource conflicts between trains can amplify the initial incident and send ripple effects through the whole network. Understanding and predicting this propagation phenomenon is a crucial task for traffic management and robust scheduling.

To model it, we construct a network graph $G = (V, E)$ linking the railway stations, and we assume the existence of a hidden congestion variable $X_{t,d}$ that lives on the edges $d \in E$. This congestion evolves according to a VAR process, whose transition matrix $\theta$ represents pairwise interactions between edges. The sparsity structure of $\theta$ expresses the local nature of delay propagation, which is why it is closely related to the adjacency structure of $G$. Indeed, between times $t$ and $t+1$, edges are expected to transmit congestion to their close neighbors, and not to regions of the network that are very far away.

Unfortunately for us, $X_t$ is never observed directly. The only information we have is collected by the trains whenever they cross an edge of the network. The crossing time of a train is influenced by the congestion, but also by other individual factors: in this sense, our observations $Y_t$ are a noisy version of the underlying process $X_t$. Furthermore, the observations are limited in size: the dimension of $X_t$ is the number of edges $D = |E|$, while the dimension of $Y_t$ is linked to the number of trains on the timetable and the length of their respective journeys. We can thus define a random variable $\pi_{t,d}$ equal to 1 if a train crosses edge $d$ between $t$ and $t+1$, and 0 otherwise. A more realistic model would account for the possibility of multiple trains crossing an edge in the same time step, especially if the discretization interval is large. However, our binary assumption greatly simplifies exposition without betraying the qualitative behavior of the system. Crucially, this sampling mechanism exhibits temporal correlations: periods of dense traffic are likely to be followed by dense traffic, which means that the sequence of sampling variables $\pi_{t,d}$ is not independently distributed.

We recognize the framework of Equation (1), and can therefore apply the theoretical result of Equation (2). This error quantification provides useful insight on the estimation of $\theta$, which is essential to help railway operatives dimension their data sets or evaluate prediction uncertainty.

## 1.3 RELATED WORKS

The theory of VAR processes has been known for a long time: the book of Lütkepohl [2005] provides a detailed account. If we have full and noiseless observations of the process $X_t$, we can use conditional Least Squares to
conflicts between trains can amplify the initial incident and send ripple effects through the whole network. Understanding and predicting this propagation phenomenon is a crucial task for traffic management and robust scheduling.

To model it, we construct a network graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ linking the railway stations, and we assume the existence of a hidden congestion variable $X_{t,d}$ that lives on the edges $d \in \mathcal{E}$. This congestion evolves according to a VAR process, whose transition matrix $\theta$ represents pairwise interactions between edges. If $\theta_{d_1,d_2} > 0$, it means that congestion on edge $d_1$ at time $t$ can lead to congestion on edge $d_2$ at time $t+1$. A typical scenario is when edge $d_2$ lies upstream of $d_1$: trains arriving towards a congested part of the network are often forced to slow down and eventually stop, which is why the traffic jam moves in the opposite direction of traffic.

The sparsity structure of $\theta$ expresses the local nature of delay propagation, thus it is closely related to the adjacency structure of $\mathcal{G}$. Indeed, between times $t$ and $t+1$, edges are expected to transmit congestion to their close neighbors, and not to regions of the network that are very far away. Such a locality prior massively reduces the search space for $\theta$, yet it still hinges on the speed at which congestion propagates. This propagation speed is not a known quantity, and may vary geographically for reasons that are not captured by our simple model (like track configuration). In other words, the sparsity mask of $\theta$ is not available a priori.

To make matters worse, $X_t$ is never observed directly. The only information we have is collected by the trains whenever they cross an edge of the network. The crossing time of a train is influenced by the congestion, but also by other individual factors: in this sense, our observations $Y_t$ are a noisy version of the underlying process $X_t$. Furthermore, the observations are limited in size: the dimension of $X_t$ is the number of edges $D = |\mathcal{E}|$, while the dimension of $Y_t$ is linked to the number of trains on the timetable and the length of their respective journeys.

Railway timetables result from a lengthy scheduling process, and they are constantly modified in real time following stochastic incidents. To abstract away this complexity, we make a hypothesis of random sampling: a random mask $\pi_{t,d}$ equal to 1 if a train crosses edge $d$ between $t$ and $t+1$, and 0 otherwise. A more realistic model would account for the possibility of multiple trains crossing an edge in the same time step, especially if the discretization interval is large. However, our binary assumption greatly simplifies exposition without betraying the qualitative behavior of the system. Similarly, some edges are busier than others, which is why the probability of sampling $P(\pi_{t,d} = 1)$ should actually be heterogeneous within the network.

Allowing some model simplifications, we recognize the framework of Equation (1), and can therefore apply the theoretical result of Equation (2). This error quantification provides useful insight on the estimation of $\theta$, which is essential to help railway operatives dimension their data sets or evaluate prediction uncertainty. Details on the concrete model, as well as an application of partially-observed VAR estimation to a dataset of tramway arrival times, can be found in Dalle [8, Chapter 10]. Since this paper is concerned with the underlying theory, we only present experiments on synthetic data.

## 1.3. Related works

The theory of VAR processes has been known for a long time: the book of Lütkepohl [27] provides a detailed account. If we have full and noiseless observations of the process $X_t$, we
estimate $\theta$ by minimizing the quadratic error $\sum_t ||X_t - \theta X_{t-1}||_2^2$. This is equivalent to solving the Yule-Walker equation $\Gamma_h = \theta\Gamma_{h-1}$, where we replace the autocovariance matrix $\Gamma_h = \text{Cov}[X_{t+h}, X_t]$ with its empirical counterpart $\hat{\Gamma}_h$. In the case of Gaussian innovations, both approaches coincide with the Maximum Likelihood Estimator (MLE).

Neither of these methods was initially designed for missing or noisy data. Luckily, statistical estimation with imprecise measurements has been thoroughly studied [Buonaccorsi, 2010]. The same goes for incomplete data sets; an extensive survey was recently published by Little and Rubin [2019]. According to their terminology, our work deals with data that is missing completely at random (MCAR), which means that the projection $\pi_t$ is independent from the underlying process $X_t$. We also assume to know the distribution of the missingness indicators $1 - \pi_{t,d}$, which is not necessarily true for other applications (e.g. clinical trials).

A principled approach to deal with missing data would require extending the MLE to partially-observed time series, also known as state-space models [Cappé et al., 2006]. Most of the time, exact or approximate inference is achievable using some version of the Kalman filter [Kalman, 1960] or particle methods [Doucet et al., 2000], whereas parameter estimation typically involves the Expectation-Maximization (EM) algorithm [Shumway and Stoffer, 1982]. Unfortunately, the EM algorithm is hard to analyze explicitly in terms of statistical error, which is why other methods are sometimes preferred in theoretical studies. In particular, plug-in methods that use covariance estimates within the Yule-Walker equation have been quite popular in the machine learning community.

In this line of work, the core challenge is the high dimension $D$ of the VAR process $X_t$. To address it, many authors use sparsity-inducing penalties as a way to reduce data requirements and computational workload. In the last ten years, the LASSO [Tibshirani, 1996] has been increasingly applied to random designs exhibiting correlations or missing data. This trend started with the seminal work of Loh and Wainwright [2012], and numerous other papers followed [see for example Basu and Michailidis, 2015, Kock and Callot, 2015, Melnyk and Banerjee, 2016, Jalali and Willett, 2018].

As an alternative to the LASSO, the Dantzig selector [Candes and Tao, 2007] enforces sparsity in the objective and data fidelity in the constraints. While the LASSO requires solving a Quadratic Program (QP), for instance with proximal methods, the Dantzig selector gives rise to a Linear Program (LP) which can be parallelized across dimensions. Han et al. [2015] studied its application to VAR estimation, obtaining finite-sample error bounds with very natural hypotheses. A little later, Rao et al. [2017a] extended these results to the more general scenario in which a hidden VAR process is randomly sampled or projected, and then corrupted with noise. This last work is quite similar to ours, but we think that the proof they present to control the non-asymptotic error is incomplete at best<fnref n="1" />.

Another salient feature of our paper is the search for a minimax lower bound, which allows us to prove the optimality of our convergence rates. To the best of our knowledge, this was only attempted once for partially-observed VAR processes. Rao et al. [2017b] presented a lower bound on the minimax error in a setting very similar to ours, but their result is less generic in several regards. Indeed, we account for the possibility of temporal correlations within sampling, as well as observation noise. Moreover, unlike the one proposed by Rao et al. [2017b], our proof focuses on geometric properties and doesn't make use of the admissible set of transition matrices until the very end. This makes it easy to handle many different types of structured transitions without additional work: sparse, Toeplitz, banded, etc.

Finally, the error bounds we obtain are backed up by detailed numerical experiments on simulated data, which allow us to visualize the influence of every parameter of interest.

## 1.4 OUTLINE OF THE PAPER

In Section 2, we define the generative procedure behind the partially-observed VAR process, and we present a sparse estimator of the transition matrix. We then state both of our theoretical results in Section 3: an upper-bound on the error of our specific estimator, complemented by a minimax lower bound on the error of any estimation algorithm. Section 4 contains numerical experiments demonstrating the impact of various parameters, which lead to the conclusion in Section 5.

<footnote n="1">¹Indeed, the combination of discrete and Gaussian concentration inequalities as performed on page 2 (middle of right column) of the supplementary material for Rao et al. [2017a] glosses over the fact that $L_F$ is itself a random variable. As we will discover during our own proof, this introduces an additional difficulty and forces us to use a more complex Gaussian concentration result (Lemma 37). See https://web.stanford.edu/~milind/papers/system_id_icassp_proof.pdf for the supplementary material in question.</footnote>
can use conditional Least Squares to estimate $\theta$ by minimizing the quadratic error $\sum_t \|X_t - \theta X_{t-1}\|_2^2$. This is equivalent to solving the Yule-Walker equation $\Gamma_h = \theta \Gamma_{h-1}$, where we replace the autocovariance matrix $\Gamma_h = \text{Cov}[X_{t+h}, X_t]$ with its empirical counterpart $\hat{\Gamma}_h$. In the case of Gaussian innovations, both approaches coincide with the Maximum Likelihood Estimator (MLE).

Neither of these methods was initially designed for missing or noisy data. Luckily, statistical estimation with imprecise measurements has been thoroughly studied [4]. The same goes for incomplete data sets; an extensive survey was recently published by Little and Rubin [25]. According to their terminology, our work deals with data that is missing completely at random (MCAR), which means that the projection $\pi_t$ is independent of the underlying process $X_t$. We also assume to know the distribution of the missingness indicators $1 - \pi_{t,d}$, which is not necessarily true for other applications (e.g. clinical trials).

A principled approach to deal with missing data would require extending the MLE to partially-observed time series, also known as state-space models [6]. Most of the time, exact or approximate inference is achievable using some version of the Kalman filter [22] or particle methods [12], whereas parameter estimation typically involves the Expectation-Maximization (EM) algorithm [34]. Unfortunately, the EM algorithm is hard to analyze explicitly in terms of statistical error, which is why other methods are sometimes preferred in theoretical studies. In particular, plug-in methods that use covariance estimates within the Yule-Walker equation have been quite popular in the machine learning community.

In this line of work, the core challenge is the high dimension $D$ of the VAR process $X_t$. To address it, many authors use sparsity-inducing penalties as a way to reduce data requirements and computational workload. In the last ten years, the LASSO [36] has been increasingly applied to random designs exhibiting correlations or missing data. This trend started with the seminal work of Loh and Wainwright [26], and numerous other papers followed [see for example 1, 23, 30, 21].

As an alternative to the LASSO, the Dantzig selector [5] enforces sparsity in the objective and data fidelity in the constraints. While the LASSO requires solving a Quadratic Program (QP), for instance with proximal methods, the Dantzig selector gives rise to a Linear Program (LP) which can be parallelized across dimensions. Han, Lu and Liu [17] studied its application to VAR estimation, obtaining finite-sample error bounds with very natural hypotheses. Note that unlike standard proofs for the LASSO in high-dimension, their argument does not require any hypotheses on the design matrix (for instance restricted eigenvalues) and it is taylored for dealing with correlated design matrices from a VAR process. In that latter case, the design matrix is squared and given by the lag-$h_0$ covariance matrix. A little later, Rao et al. [31] extended the results of Han, Lu and Liu [17] to the more general scenario in which a hidden VAR process is randomly sampled or projected, and then corrupted with noise. This last work is quite similar to ours, but we think that the proof they present to control the non-asymptotic error is incomplete (see below).

A related model is the so-called Factor-Augmented VAR (FAVAR), introduced by Bernanke, Boivin and Eliasz [2]. In their framework, a set of latent factors $X_t$ and observed covariates $C_t$ jointly follow a VAR process, and jointly influence the observed sequence $Y_t$. It can be seen as a variant of our partially-observed VAR, where the observation mask $\pi_t$ would be deterministic and always hide the first block of the state vector $(X_t, C_t)$. Our stochastic
Appendix A is dedicated to proving the convergence rate of the sparse estimator, while Appendix B contains the derivation of the minimax lower bound. A number of useful results from linear algebra and probability are presented in Appendix C to make the paper as self-contained as possible. Most of them are well-known, some were obtained or adapted specifically for our proof. Appendix D contains a summary of the main notations and symbols.

# 2 THE PARTIALLY-OBSERVED VAR PROCESS AND ITS SPARSE ESTIMATOR

Before stating our theoretical results, we introduce our statistical model and the estimator we use.

## 2.1 MODEL DEFINITION

The model we study was described approximately in the introduction. We now fill the gaps of the generative procedure it relies on.

THE UNDERLYING STATE $X = (X_t)_{t \in [T]}$ follows a stationary VAR process of order 1. This process has dimension $D$ and the following recursive definition:

$$X_t = \theta X_{t-1} + \varepsilon_t \quad \text{with} \quad \varepsilon_t \sim \mathcal{N}(0, \Sigma). \quad (3)$$

Here $\theta \in \mathbb{R}^{D \times D}$ is the transition matrix and $\Sigma \in \mathbb{R}^{D \times D}$ is the covariance matrix of the innovations (in the introduction, we assumed $\Sigma = \sigma^2 I$).

To ensure stationarity of the VAR process, we must constrain the spectral radius of $\theta$ to satisfy $\rho(\theta) < 1$. Throughout the paper, we actually make the following (slightly stronger) assumption on the spectral norm of $\theta$: there exists $\vartheta \in (0, 1)$ such that for all the values of $\theta$ we consider, $\|\theta\|_2 \le \vartheta < 1$. Furthermore, we only study row-sparse transition matrices, having at most $s$ nonzero coefficients in each row. In other words, we restrict our choice of parameters to

$$\theta \in \Theta_s \quad \text{where} \quad \Theta_s = \{\theta \in \mathbb{R}^{D \times D} : \|\theta\|_2 \le \vartheta < 1 \quad \text{and} \quad \forall i, \|\theta_{i, \cdot}\|_0 \le s\}. \quad (4)$$

We denote by $\sigma_{\min}^2 = \lambda_{\min}(\Sigma)$ and $\sigma_{\max}^2 = \lambda_{\max}(\Sigma)$ the minimum and maximum eigenvalues of the covariance matrix $\Sigma$.

THE OBSERVATION MECHANISM we chose implies that we do not have direct access to the latent process $X_t$. To construct the observations $Y_t$, we sample a subset of state components according to the binary vectors $\pi_t$. Then, independent Gaussian noise with variance $\omega^2$ is added to these selected components, and we observe the result. If we denote by $\Pi_t = \text{diag}(\pi_t)$ the diagonal projection matrix, we have

$$Y_t = \Pi_t X_t + \eta_t \quad \text{with} \quad \eta_t \sim \mathcal{N}(0, \omega^2 I). \quad (5)$$

An essential hypothesis we make is the mutual independence between our three sources of randomness: the innovations $\varepsilon_t$, the projections $\pi_t$ and the observation noise $\eta_t$.

A major feature of the present work is the non-deterministic selection of observed state components, that is, the fact that $\pi_t$ is a random sequence of Bernoulli vectors following a known distribution. In order to sum up the amount of information available using one parameter $p \in (0, 1)$, we want this distribution to satisfy the following condition: each component $X_{t,d}$ of the latent state must be sampled with the same marginal probability $p = \mathbb{P}(\pi_{t,d} = 1)$.

On the other hand, we also want to introduce temporal dependencies between the projections. The simplest way to achieve that is through a Markovian hypothesis: independently along each dimension $d$, time indices $t$ are selected for observation according to a binary-valued Markov chain with transition matrix $T = \begin{pmatrix} 1-a & a \\ b & 1-b \end{pmatrix}$. Its coefficients are chosen to make the chain stationary with invariant measure $(\frac{b}{a+b}, \frac{a}{a+b}) = (1-p, p)$. Note that when $a = 1-b = p$, this reduces to independent sampling of each component with probability $p$. We also assume there exists a universal constant $\chi$ such that $0 < \chi \le a, b \le 1-\chi < 1$: this means that the chain does not transition too fast nor too slowly.

Our data set is built from $N$ independent realizations of this process. For the sake of simplicity however, we will prove all convergence theorems in the case $N = 1$: extending those results to the general case simply amounts to replacing $T$ with $NT$ in the resulting error bounds.
analysis does not generalize though, because we assume the mask marginals are all identical
Bernoullis $\pi_{t,d} \sim \mathcal{B}(p)$. A further divergence lies in the dimensionality assumptions: even in
recent works like Lin and Michailidis [24], the size of the latent vector is assumed to remain
small, whereas we allow for a high-dimensional latent space thanks to transition sparsity.

Another salient feature of our paper is the search for a minimax lower bound, which allows us to prove the optimality of our convergence rates. To the best of our knowledge, this was only attempted once for partially-observed VAR processes. Rao et al. [32] presented a lower bound on the minimax error in a setting very similar to ours, but their result is less generic in several regards. Indeed, we account for the possibility of temporal correlations within sampling, as well as observation noise. Moreover, unlike the one proposed by Rao et al. [32], our proof focuses on geometric properties and doesn't make use of the admissible set of transition matrices until the very end. This makes it easy to handle many types of structured transitions without much additional work: sparse, Toeplitz, banded, etc.

Finally, the error bounds we obtain are backed up by detailed numerical experiments on
simulated data, which allow us to visualize the influence of every parameter of interest.

**Remark** (Correcting incomplete previous proofs). As mentioned above, Rao et al. [31, 32] previously claimed a similar result, but we think their upper bound proof is incomplete. To obtain a concentration result, we work with a combination of continuous randomness (from the innovations and noise) and discrete randomness (from the sampling mask). In our proof, we first control the Bernoulli variables, and then control the Gaussian variables conditionally. For this, we leverage a conditional version of the Hanson-Wright inequality (Lemma C.12) which is new (to the best of our knowledge) and of independent interest. Meanwhile, their proof amounts to using the unconditional Hanson-Wright inequality instead (Lemma C.11), which is missing a term compared to its conditional counterpart. We give more details in Appendix A.6.

**1.4. Outline of the paper**

In Section 2, we define the generative procedure behind the partially-observed VAR process, and we present a sparse estimator of the transition matrix. We then state both of our theoretical results in Section 3: an upper-bound on the error of our specific estimator, complemented by a minimax lower bound on the error of any estimation algorithm. Section 4 contains numerical experiments demonstrating the impact of various parameters, which lead to the conclusion in Section 5.

Appendix A is dedicated to proving the convergence rate of the sparse estimator, while Appendix B contains the derivation of the minimax lower bound. A number of useful results from linear algebra and probability are presented in Appendix C to make the paper as self-contained as possible. Most of them are well-known, some were obtained or adapted specifically for our proof.

**2. The partially-observed VAR process and its sparse estimator**

Before stating our theoretical results, we introduce our statistical model and the estimator we
use.
## 2.2 SPARSE ESTIMATOR FOR THE TRANSITION MATRIX

We now introduce the estimation method chosen for this problem.

THE TRANSITION ESTIMATOR presented here is a straightforward generalization of the one used by Rao et al. [2017a]. The lag-h covariance matrix of the VAR process $X_t$ is given by the Yule-Walker recursion (see Lemma 1):

$$ \Gamma_h(\theta) = \mathrm{Cov}_{\theta}[X_{t+h}, X_t] = \theta \Gamma_{h-1}(\theta) = \theta^h \Gamma_0(\theta) \quad (6) $$

We can use it to define a simple two-step procedure:

1. For a given $h_0$, build estimators $\hat{\Gamma}_{h_0}$ and $\hat{\Gamma}_{h_0+1}$ of the covariances $\Gamma_{h_0}$ and $\Gamma_{h_0+1}$.

2. Use them to approximate the transition matrix by inverting Equation (6).

A simple inversion technique uses the Moore-Penrose pseudoinverse (just in case $\hat{\Gamma}_{h_0}$ is singular):

$$ \hat{\theta}^{\mathrm{dense}} = \hat{\Gamma}_{h_0+1} \hat{\Gamma}_{h_0}^{\dagger}. \quad (7) $$

The problem with this procedure is that it does not guarantee sparsity of $\hat{\theta}$. To obtain a sparse result, we follow Han et al. [2015] and cast Equation (6) as a soft constraint enforcing proximity between $\hat{\Gamma}_{h_0+1}$ and $\hat{\theta}\hat{\Gamma}_{h_0}$. This amounts to solving the following constrained optimization problem:

$$ \hat{\theta} \in \underset{M \in \mathbb{R}^{D \times D}}{\operatorname{argmin}} \|\mathrm{vec}(M)\|_1 \quad \text{subject to} \quad \|M\hat{\Gamma}_{h_0} - \hat{\Gamma}_{h_0+1}\|_{\max} \le \lambda. \quad (8) $$

Here $\|\mathrm{vec}(\cdot)\|_1$ denotes the sum of the absolute values of all the coefficients of a matrix, while $\|\cdot\|_{\max}$ is the maximum of these absolute values. Given that both of these norms are piecewise linear, the problem of Equation (8) can be reformulated as an LP. It can even be decomposed along each dimension, which allows for an efficient and parallel solution procedure. The only thing left to do is decide how to estimate the covariance matrices $\Gamma_h$.

THE COVARIANCE ESTIMATOR we use is a variant of the empirical covariance. Since $Y_t = \Pi_t X_t + \eta_t$ where $\eta_t$ is zero-mean, a natural proxy for $X_t$ is obtained by inverting the sampling operator: $\hat{X}_t = \Pi_t^\dagger Y_t$. It would therefore seem logical to build an estimator of $\Gamma_h$ by plugging this proxy into the empirical covariance between $X_{t+h}$ and $X_t$. However, in order for this idea to work, we must make two small adjustments.

To account for the random sampling, the plug-in empirical covariance must be scaled elementwise by a matrix $S(h) = \mathbb{E}[\pi_{t+h} \pi_t'$. Intuitively, since $\hat{X}_{t+h} \hat{X}_t'$ has a fraction $p^2$ of nonzero coefficients, we need to divide it by something close to $p^2$ to get an unbiased covariance estimator. Furthermore, to account for the observation noise, we must incorporate an additive correction $-\omega^2 I$. This correction becomes unnecessary for $h \ge 1$ since the observation noise $\eta_t$ is independent across time.

In conclusion, we obtain the following covariance estimator:

$$ \hat{\Gamma}_h = \frac{1}{S(h)} \odot \frac{1}{T-h} \sum_{t=1}^{T-h} \left( \Pi_{t+h}^{\dagger} Y_{t+h} \right) \left( \Pi_t^{\dagger} Y_t \right)' - \mathbf{1}_{\{h=0\}} \omega^2 I. \quad (9) $$

The coefficients of the scaling matrix $S(h)$ are computed in Lemma 4.

## 3 LOWER AND UPPER BOUND ON THE ESTIMATION ERROR

We now have the necessary background to formulate our theoretical results. In all the following statements (and their proofs), the letter $c$ denotes a universal positive constant, which may change from one line to the next but never depends on any varying problem parameters. More specifically, statements involving it should always be understood as “there exists $c > 0$ such that”...
## 2.1. Model definition

The model we study was described approximately in the introduction. We now fill the gaps of the generative procedure it relies on.

**The underlying state** $X = (X_t)_{t \in [T]}$ follows a stationary VAR process of order 1. This process has dimension $D$ and the following recursive definition:

$$X_t = \theta X_{t-1} + \varepsilon_t \quad \text{with} \quad \varepsilon_t \sim N(0, \sigma^2 I). \quad (3)$$

Here $\theta \in \mathbb{R}^{D \times D}$ is the transition matrix and $\sigma^2$ is the variance of the innovations.

To ensure stationarity of the VAR process, we must constrain the spectral radius of $\theta$ to satisfy $\rho(\theta) < 1$. Throughout the paper, we actually make the following (slightly stronger) assumption on the spectral norm (or operator $\ell_2$ norm) of $\theta$, denoted by $\|\theta\|_2$: there exists $\vartheta \in (0, 1)$ such that for all the values of $\theta$ we consider, $\|\theta\|_2 \le \vartheta < 1$. Furthermore, we only study row-sparse transition matrices, having at most $s$ nonzero coefficients in each row. This is equivalent to constraining the $\ell_0$ norm of each row $i$, denoted by $\|\theta_{i, \cdot}\|_0$, to be smaller than $s$. In other words, we restrict our choice of parameters to

$$\theta \in \Theta_s \quad \text{where} \quad \Theta_s = \{\theta \in \mathbb{R}^{D \times D} : \|\theta\|_2 \le \vartheta < 1 \quad \text{and} \quad \forall i, \|\theta_{i, \cdot}\|_0 \le s\}. \quad (4)$$

**The observation mechanism** we chose implies that we do not have direct access to the latent process $X_t$. To construct the observations $Y_t$, we sample a subset of state components according to the binary vectors $\pi_t$. Then, independent Gaussian noise with variance $\omega^2$ is added to these selected components, and we observe the result. If we denote by $\Pi_t = \text{diag}(\pi_t)$ the diagonal projection matrix, we have

$$Y_t = \Pi_t X_t + \eta_t \quad \text{with} \quad \eta_t \sim N(0, \omega^2 I). \quad (5)$$

An essential hypothesis we make is the mutual independence between our three sources of randomness: the innovations $\varepsilon_t$, the projections $\pi_t$ and the observation noise $\eta_t$.

A major feature of the present work is the non-deterministic selection of observed state components, that is, the fact that $\pi_t$ is a random sequence of Bernoulli vectors following a known distribution. In order to sum up the amount of information available using one parameter $p \in (0, 1)$, we want this distribution to satisfy the following condition: each component $X_{t,d}$ of the latent state must be sampled independently with the same marginal probability $p = \mathbb{P}(\pi_{t,d} = 1)$.

Our data set is built from $N$ independent realizations of this process. For the sake of simplicity however, we will prove all convergence theorems in the case $N = 1$: extending those results to the general case simply amounts to replacing $T$ with $NT$ in the resulting error bounds.

## 2.2. Sparse estimator for the transition matrix

We now introduce the estimation method chosen for this problem.
## 3.1 MAIN THEOREMS

We start by bounding the non-asymptotic error of the estimator we just introduced.

**Theorem 1** (Error upper bound). Consider the partially-observed VAR model defined in Section 2.1. We use the estimator $\hat{\theta}$ of Section 2.2 with $h_0 = 0$, and we suppose that $T$ is “large enough”, as specified by Equations (20) and (23). Let us define

$$ \gamma_u(\theta) = \frac{\|\theta\|_{\infty} + 1}{(1 - \|\theta\|_2)^2} \frac{\sigma_{\max}^2 + \omega^2}{\|\Gamma_0(\theta)^{-1}\|_1^{-1}} \quad \text{and} \quad q_u = \min\{p, 1-b\} \le p. \quad (10) $$

Then there is a value of $\lambda$ such that the following upper bound holds with probability at least $1 - \delta$:

$$ \|\hat{\theta} - \theta\|_{\infty} \le c \frac{\gamma_u(\theta)s}{\sqrt{Tpq_u}} \sqrt{\log(D/\delta)}. \quad (11) $$

*Proof.* The argument combines discrete and continuous concentration inequalities, to account for both the Bernoulli sampling and the Gaussian noise. More precisely, we exploit a recent Chernoff bound that applies to non-reversible Markov chains, and we plug it into a conditional version of the Hanson-Wright inequality that we derived specifically for our purposes. See Appendix A for more details. □

We now move on to a minimax lower bound which is estimator-independent, and quantifies the intrinsic difficulty of our statistical problem. The term minimax means that we study the probability of making an error of magnitude $\zeta$, when we pick the best possible estimator $\hat{\theta}$ and nature replies by choosing the worst possible parameter $\theta$:

$$ \mathfrak{P}(\zeta) = \inf_{\hat{\theta}} \sup_{\theta \in \Theta_s} \mathbb{P}_{\theta} \left[ \|\hat{\theta} - \theta\|_{\infty} \ge \zeta \right]. \quad (12) $$

More precisely, we want to find a threshold $\zeta$ such that the probability of exceeding it is non-negligible, for instance $\mathfrak{P}(\zeta) \ge \frac{1}{2}$. The evolution of this threshold will tell us how the error behaves with respect to the various problem parameters.

**Theorem 2** (Error lower bound). Consider the partially-observed VAR model defined in Section 2.1. We suppose that $T$ is “large enough”, as specified by Equations (25) and (27). Let us define

$$ \gamma_{\ell} = (1 - \vartheta)^{3/2} \frac{\sigma_{\min}^2 + \omega^2}{\sigma_{\max}^2} \quad \text{and} \quad q_{\ell} = \max\{1 - b, 2p - (1 - b)\} \ge p. \quad (13) $$

Then the following minimax lower bound holds:

$$ \inf_{\hat{\theta}} \sup_{\theta \in \Theta_s} \mathbb{P}_{\theta} \left[ \|\hat{\theta} - \theta\|_{\infty} \ge c \frac{\gamma_{\ell} s}{\sqrt{Tpq_{\ell}}} \right] \ge \frac{1}{2}. \quad (14) $$

*Proof.* The argument is based on an information-theoretical result known as Fano's inequality. To apply it, we need to upper bound the Kullback-Leibler (KL) divergence between the distributions $\mathbb{P}_{\theta_0}(\Pi, Y)$ and $\mathbb{P}_{\theta_1}(\Pi, Y)$, where $\theta_0$ and $\theta_1$ are sufficiently far apart. See Appendix B for more details. □

## 3.2 INFLUENCE OF THE PROBLEM PARAMETERS

Let us now compare the error bounds of Theorems 1 and 2. Our first remark is that $s$ and $T$ play exactly the same roles in both bounds (up to a logarithmic factor), which shows that the dependency of the error in $s/\sqrt{T}$ is optimal.

THE SAMPLING PARAMETERS appear as $1/\sqrt{pq_u}$ in the upper bound, whereas the lower bound scales as $1/\sqrt{pq_{\ell}}$ instead. This means that we have not proven the optimality of either bound with respect to $p$ or $b$. However, it is reassuring to note that there is no conflict between them since $q_{\ell} \ge p \ge q_u$. Furthermore, when $a = 1 - b = p$ (that is, when Markov sampling boils down to independent sampling), both bounds simplify into the $1/p$ dependency we would expect (since $q_u = q_{\ell} = p$). So in the case of independent sampling, $1/p$ is indeed the optimal rate.
**The transition estimator** presented here is a straightforward generalization of the one used by Rao et al. [31]. We recall that the covariance matrix between two random vectors is defined as $\text{Cov}[U, V] = \mathbb{E}[(U - \mathbb{E}[U])(V - \mathbb{E}[V])']$, where we use an apostrophe to denote matrix transposition. The lag-$h$ autocovariance matrix of the VAR process $X_t$ is given by the Yule-Walker recursion (see Lemma A.1):

$$ \Gamma_h(\theta) = \text{Cov}_{\theta}[X_{t+h}, X_t] = \theta\Gamma_{h-1}(\theta) = \theta^h\Gamma_0(\theta) \quad (6) $$

We can use it to define a simple two-step procedure:

1. For a given $h_0$, build estimators $\hat{\Gamma}_{h_0}$ and $\hat{\Gamma}_{h_0+1}$ of the covariances $\Gamma_{h_0}$ and $\Gamma_{h_0+1}$.

2. Use them to approximate the transition matrix by inverting Equation (6).

A simple inversion technique uses the Moore-Penrose pseudoinverse (just in case $\hat{\Gamma}_{h_0}$ is singular):

$$ \hat{\theta}^{\text{dense}} = \hat{\Gamma}_{h_0+1} \hat{\Gamma}_{h_0}^{\dagger}. \quad (7) $$

This is the method of choice when the transition matrix is small and dense. In high dimension however, computing the pseudoinverse requires a matrix factorization, which can be costly. Moreover, Equation (7) does not guarantee sparsity of $\hat{\theta}$. To obtain a sparse result, we follow Han, Lu and Liu [17] and cast Equation (6) as a soft constraint enforcing proximity between $\hat{\Gamma}_{h_0+1}$ and $\hat{\theta}\hat{\Gamma}_{h_0}$. This amounts to solving the following constrained optimization problem, known as the Dantzig selector:

$$ \hat{\theta} \in \underset{M \in \mathbb{R}^{D \times D}}{\operatorname{argmin}} \|\operatorname{vec}(M)\|_1 \quad \text{subject to} \quad \|M\hat{\Gamma}_{h_0} - \hat{\Gamma}_{h_0+1}\|_{\max} \le \lambda. \quad (8) $$

Here $\|\operatorname{vec}(\cdot)\|_1$ denotes the sum of the absolute values of all the coefficients of a matrix, while $\|\cdot\|_{\max}$ is the maximum of these absolute values. Given that both of these norms are piecewise linear, the problem of Equation (8) can be reformulated as an LP. It can even be decomposed along each dimension, which allows for an efficient and parallel solution procedure. The only thing left to do is decide how to estimate the covariance matrices $\Gamma_h$. Note that for $\lambda = 0$, we recover $\hat{\theta} = \hat{\theta}_{\text{dense}}$ (at least when the covariance estimate $\hat{\Gamma}_0$ is invertible).

**The covariance estimator** we use is a variant of the empirical covariance. Since $Y_t = \Pi_t X_t + \eta_t$ where $\eta_t$ is zero-mean, a natural proxy for $X_t$ is obtained by inverting the sampling operator: $\hat{X}_t = \Pi_t^\dagger Y_t$. It would therefore seem logical to build an estimator of $\Gamma_h$ by plugging this proxy into the empirical covariance between $X_{t+h}$ and $X_t$. However, in order for this idea to work, we must make two small adjustments.

To account for the random sampling, the plug-in empirical covariance must be scaled elementwise by a matrix $S(h) = \mathbb{E}[\pi_{t+h}\pi_t']$. Intuitively, since $\hat{X}_{t+h}\hat{X}_t'$ has a fraction $p^2$ of nonzero coefficients, we need to divide it by something close to $p^2$ to get an unbiased covariance estimator. Furthermore, to account for the observation noise, we must incorporate an additive correction $-\omega^2 I$. This correction becomes unnecessary for $h \ge 1$ since the observation noise $\eta_t$ is independent across time. We denote this special case with an indicator function $1_{\{h=0\}}$.
THE $\ell_2$ NORM OF THE TRANSITION MATRIX plays opposite roles on each side. In the lower bound, $1 - \vartheta = 1 - \max_{\theta \in \Theta_s} \|\theta\|_2$ appears in the numerator, whereas in the upper bound, $1 - \|\theta\|_2$ appears in the denominator. It is likely that these dependencies are suboptimal, but at least they are compatible with one another: as $\|\theta\|_2 \to 1$, that is, as the VAR process becomes unstable, the lower bound tends to 0 and the upper bound to $+\infty$. This is a reflection of the fact that our proofs make heavy use of the distance between $\theta$ and the unit sphere, which means they become meaningless when $\theta$ gets too large.

THE VARIANCES $\Sigma$ AND $\omega^2$ are involved in $\gamma_\ell$ for the lower bound, and in $\gamma_u(\theta)$ for the upper bound. In both cases, the ratio $\gamma$ tells us whether the underlying process is large enough to be detected among the noise. Roughly speaking, the magnitude of $X_t$ is related to the spectrum of $\Sigma$, while the magnitude of $Y_t$ is related to the spectrum of $\Sigma + \omega^2 I$. If the latter is significantly larger than the former, recovering $X_t$ (and thus $\theta$) is a hopeless endeavor.

To simplify the comparison, let us assume in this discussion that $\Sigma = \sigma^2 I$, and that $\theta$ commutes with its transpose. Then we have $\|\Gamma_0^{-1}(\theta)\|_1^{-1} = \|(\sigma^2(I - \theta\theta')^{-1})^{-1}\|_1^{-1} = \sigma^2\|I - \theta\theta'\|_1^{-1}$, and we can give a simpler expression of $\gamma_\ell$ and $\gamma_u(\theta)$:

$$ \gamma_u(\theta) = \frac{(\|\theta'\|_1 + 1)\|I - \theta\theta'\|_1 \sigma^2 + \omega^2}{(1 - \vartheta)^2 \sigma^2} \qquad \gamma_\ell = (1 - \vartheta)^{3/2} \frac{\sigma^2 + \omega^2}{\sigma^2}. $$

We recognize the same dependency in both bounds, namely $\gamma \propto 1 + \frac{\sigma^2}{\omega^2}$. Lemma 38 gives a heuristic argument linking this functional form to the asymptotic behavior of the MLE.

### 3.3 EXTENSION TO VAR PROCESSES OF HIGHER ORDER

Although our results only apply to state-space models based on an underlying VAR process of order 1, we could try to extend them to the more general case of VAR($K$) processes. Just for this Section, suppose $X_t$ is no longer given by Equation (3), but instead satisfies:

$$ X_t = \theta_1 X_{t-1} + \theta_2 X_{t-2} + \dots + \theta_K X_{t-K} + \varepsilon_t. $$

Then we can represent this as a VAR(1) process using augmented variables [Lütkepohl, 2005]. Indeed, observe that defining $\tilde{X}_t = (X_t \ X_{t-1} \ \dots \ X_{t-K+1})'$ and $\tilde{\varepsilon}_t = (\varepsilon_t \ 0 \ \dots \ 0)'$ yields

$$ \tilde{X}_t = \tilde{\theta} \tilde{X}_{t-1} + \tilde{\varepsilon}_t \quad \text{with} \quad \tilde{\theta} = \begin{bmatrix} \theta_1 & \theta_2 & \cdots & \theta_{K-1} & \theta_K \\ I_D & 0 & \cdots & 0 & 0 \\ 0 & I_D & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & I_D & 0 \end{bmatrix}. $$

Unfortunately, by this reasoning, the Markov sampling mechanism that generates $\tilde{\Pi}_t$ gives rise to a new distribution for $\tilde{\Pi}_t$ which is no longer part of the same family. Indeed, the augmented sampling process $\tilde{\Pi}_t$ is still Markovian but with a memory of size $K$ instead of 1. Therefore, the adaptation would require new arguments and we leave it for future work.

## 4 NUMERICAL ILLUSTRATIONS

We now illustrate our results on simulated data. All experiments were performed on a Dell Precision 5530 mobile workstation with Intel Core i7-8850H CPU (2.60GHz × 12) and 31 GiB of RAM, running under Ubuntu 20.04. Our code was written in Julia [Bezanson et al., 2017], linear optimization problems were modeled using JuMP [Dunning et al., 2017] and solved with the COIN-OR Clp solver [Forrest et al., 2022]. The reproducible Pluto notebook used to generate all the plots will be made available on GitHub as soon as the review procedure is complete and anonymity is no longer required.
In conclusion, we obtain the following covariance estimator, where $\odot$ is the elementwise multiplication and $\frac{1}{S(h)}$ is the elementwise inverse:

$$ \hat{\Gamma}_h = \frac{1}{S(h)} \odot \frac{1}{T-h} \sum_{t=1}^{T-h} (\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' - \mathbf{1}_{\{h=0\}} \omega^2 I. \quad (9) $$

The coefficients of the scaling matrix $S(h)$ are computed in Lemma A.4.

**Remark (Practicality of the estimator).** The estimator we study here has some drawbacks which could hinder its use in practice. Not only does it require a linear programming solver, it also demands the knowledge of the noise level $\omega$ to estimate the lag-0 covariance $\hat{\Gamma}_0$. While this latter limitation could be lifted by picking a higher lag $h_0$, we do not address it here because our main focus is the estimation of $\theta$. In this setting, the Dantzig selector allows for a very straightforward convergence analysis [5, 17], which justifies our choice.

### 3. Lower and upper bound on the estimation error

We now have the necessary background to formulate our theoretical results. In all the following statements (and their proofs), the letter $c$ denotes a universal positive constant, which may change from one line to the next but never depends on any varying problem parameters. More specifically, statements involving it should always be understood as “there exists $c > 0$ such that”...

#### 3.1. Main theorems

We start by bounding the non-asymptotic error of the estimator we just introduced.

**Theorem 3.1 (Error upper bound).** Consider the partially-observed VAR model defined in Section 2.1. We use the estimator $\hat{\theta}$ of Section 2.2 with $h_0 = 0$, and we suppose that $T$ is “large enough” (see Remark 3.1). Let us assume that $\theta$ commutes with its transpose, and define

$$ \gamma_u(\theta) = \frac{(\|\theta\|_\infty + 1)\|I - \theta\theta'\|_1}{(1 - \|\theta\|_2)^2} \left(1 + \frac{\omega^2}{\sigma^2}\right). \quad (10) $$

where $\|\cdot\|_1$ is the operator $\ell_1$ norm. Then there is a value of $\lambda$ such that the following upper bound holds with probability at least $1 - \delta$:

$$ \|\hat{\theta} - \theta\|_{\infty} \le c \frac{\gamma_u(\theta)s}{\sqrt{T}p} \sqrt{\log(D/\delta)} \quad (11) $$

*Proof.* The argument combines discrete and continuous concentration inequalities, to account for both the Bernoulli sampling and the Gaussian noise. More precisely, we exploit a Chernoff bound for discrete variables, and we plug it into a conditional version of the Hanson-Wright inequality that we derived specifically for our purposes. The very end of the proof follows the exact same steps as Han, Lu and Liu [17], but with our own concentration bounds for the covariance estimates. It is in these bounds that part of the novelty and difficulty lies, because partial observation requires combining two sources of randomness: Gaussian variables from the VAR process, and Bernoulli variables from the sampling mechanism. See Appendix A for more details. □
## 4.1 DATA GENERATION

Simulating a partially-observed VAR process with known transition matrix $\theta$ allows us to compute the estimation error $\|\hat{\theta} - \theta\|_{\infty}$ and study the influence of parameters such as $T$, $D$, $s$, $p$, $\omega$, etc. Real values for $\theta$ were drawn using independent standard Gaussian distributions for each coefficient, and then normalized to satisfy $\|\theta\|_{2} = \vartheta = \frac{1}{2}$. To simplify comparison with the theoretical bounds, we used a diagonal innovation covariance $\Sigma = \sigma^2 I$ and set the sampling parameters to $a = 1-b = p$, which amounts to independent sampling (except for the experiment that focuses specifically on the influence of $b$). When not mentioned explicitly, all other parameters are equal to their default values given below (we assume $\omega$ is known):

$$T = 10000 \quad D = 5 \quad \sigma = 1.0 \quad \omega = 0.1 \quad p = 1.0.$$

Most of the simulations are run in a dense estimation scenario. For those that require the sparse procedure, selecting a good regularization parameter $\lambda$ is paramount: indeed, Theorem 1 is only valid for a specific value of $\lambda$ (which is not known in practice, but we can hope to approximate this near-optimal choice).

A standard way to tune $\lambda$ would be cross-validation. However, evaluating a choice of $\lambda$ (and the resulting estimate $\hat{\theta}$) requires inferring the hidden state sequence $X_t$ from the observations $Y_t$. If the projection matrices $\Pi_t$ were deterministic, the inference could be performed with Kalman filtering [Kalman, 1960], but since they are stochastic, the distribution of $(X, Y)$ is no longer jointly Gaussian and the justification behind the Kalman filter breaks down. Finding an appropriate inference method in our setting will be the topic of future studies.

In the meantime, to tune $\lambda$, we suppose that the sparsity level of the real transition matrix $\theta$ is known. We then use this target sparsity $\hat{s}$ to guide a dichotomy search on $\lambda$, until we find a transition matrix estimate $\hat{\theta}$ whose row sparsity level is sufficiently close to $\hat{s}$.

## 4.2 RESULTS

The main results are presented on Figure 1. With the exception of 1f, all plots have the estimation error $\|\hat{\theta} - \theta\|_{\infty}$ on their y-axis, and some parameter of interest on their x-axis. The axes are displayed with logarithmic scaling, in order to highlight the exponent of the dependencies. Each point corresponds to one run of the algorithm, aimed at estimating a single random value of $\theta$. When a straight line is added to a scatter plot, it is the result of a Theil-Sen regression [Sen, 1968] applied to the points of the same color: its slope is denoted by $\alpha$ in the legend.

Figure 1a confirms that the error decreases as $1/\sqrt{T}$. This is only true because the sampling probability $p$ remains constant. If instead we had a limited observation budget but an increasing temporal precision, we would have $p \propto 1/T$, in which case the error would increase as $\sqrt{T}$ instead of decreasing.

Figure 1b exhibits three clearly identifiable regimes with respect to the noise variance. In the first one, corresponding to $\omega/\sigma \ll 1$, the error remains small and constant. Then, the error increases when $\omega/\sigma \simeq 1$. In the third phase, corresponding to $\omega/\sigma \gg 1$, the error remains high and volatile. This is consistent with the theoretical dependency in $1 + \omega^2/\sigma^2$.

Figure 1c compares the respective benefits of sparse and dense estimation by increasing the ambient dimension $D$ while keeping the true sparsity level $s$ constant. The error for $\hat{\theta}^{\text{dense}}$ scales linearly with $D$, while its sparse counterpart $\hat{\theta}$ achieves a much slower error growth. As a side note, the fact that the error grows with $D$ is not surprising. Indeed, we measure it with the $\ell_{\infty}$ operator norm, which scales with the dimension of the matrix.

Figure 1d takes the opposite perspective by increasing the number of nonzero coefficients in a space of fixed dimension. In this case, the theory predicts that the error should scale linearly with $s$, but the slope we observe is below 1. Our interpretation is that the function $\gamma_u(\theta)$ also depends on the sparsity level in complicated ways through $\theta$, especially since the real values are renormalized to satisfy $\|\theta\|_{2} = \frac{1}{2}$.

Figure 1e shows that the error evolves as $1/p$, which is consistent with our upper bound. It is also informative w.r.t. the choice of $h_0$. Choosing $h_0 = 0$ means we need to know $\omega$ to perform estimation. If this parameter is unknown, we can choose $h_0 \ge 1$, which leads to a much higher variance of the estimator (this is not visible in our results since we wrote the proof in the case where $h_0 = 0$). An alternate solution would be to keep $h_0 = 0$ and plug in a guess such as $\omega = 0$, effectively trading lower variance for a higher bias.
**Remark (Parameter regime).** The required size of $T$ for Theorem 3.1 is specified by Equations (19) and (22) in the Appendix. A sufficient condition for both to hold is

$$T = \Omega\left(\frac{\log(1/\delta) \max\{1, (\sigma^2 + \omega^2)^{-2}\}}{p^2(1-\vartheta)^4}\right)$$

We now move on to a minimax lower bound which is estimator-independent, and quantifies the intrinsic difficulty of our statistical problem. The term minimax means that we study the probability of making an error of magnitude $\zeta$, when we pick the best possible estimator $\hat{\theta}$ and nature replies by choosing the worst possible parameter $\theta$:

$$\mathfrak{P}(\zeta) = \inf_{\hat{\theta}} \sup_{\theta \in \Theta_s} \mathbb{P}_{\theta} \left[ \|\hat{\theta} - \theta\|_{\infty} \ge \zeta \right]. \quad (12)$$

More precisely, we want to find a threshold $\zeta$ such that the probability of exceeding it is non-negligible, for instance $\mathfrak{P}(\zeta) \ge \frac{1}{2}$. The evolution of this threshold will tell us how the error behaves with respect to the various problem parameters.

**Theorem 3.2 (Error lower bound).** Consider the partially-observed VAR model defined in Section 2.1. We suppose that $T$ is “large enough” (see Remark 3.1). Let us define

$$\gamma_{\ell} = (1 - \vartheta)^{3/2} \left(1 + \frac{\omega^2}{\sigma^2}\right) \quad (13)$$

Then the following minimax lower bound holds:

$$\inf_{\hat{\theta}} \sup_{\theta \in \Theta_s} \mathbb{P}_{\theta} \left[ \|\hat{\theta} - \theta\|_{\infty} \ge c \frac{\gamma_{\ell} s}{\sqrt{T} p} \right] \ge \frac{1}{2}. \quad (14)$$

*Proof.* The argument is based on an information-theoretical result known as Fano's inequality. To apply it, we need to upper-bound the Kullback-Leibler (KL) divergence between the distributions $\mathbb{P}_{\theta_0}(\Pi, Y)$ and $\mathbb{P}_{\theta_1}(\Pi, Y)$, where $\theta_0$ and $\theta_1$ are sufficiently far apart. See Appendix B for more details. □

**Remark (Parameter regime).** The required size of $T$ for Theorem 3.2 is specified by Equations (24) and (26) in the Appendix. A sufficient condition for both to hold is

$$T = \Omega\left(\frac{\gamma_{\ell}^2 s^2}{p^3 \min\{\vartheta^2, \gamma_{\ell}^2(1-\vartheta)\}}\right)$$

## 3.2. Influence of the problem parameters

Let us now compare the error bounds of Theorems 3.1 and 3.2.

**The sparsity $s$, duration $T$ and sampling probability $p$** play exactly the same roles in both bounds up to a logarithmic factor in $D$, which shows that the dependency of the error in $s/p\sqrt{T}$ is essentially optimal. Note that we only obtain this dependency for sufficiently large values of $T$, with thresholds that could probably be improved.
<figure><img src="image_2.png" /><figcaption>(a) Influence of $T$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(c) Influence of $D$ with fixed $s$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(e) Influence of $p$ and $h_0$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Influence of $\omega$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) Influence of $s$ with fixed $D$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(f) Influence of $(p, b)$</figcaption></figure>

<caption>Figure 1: Impact of model parameters on the estimation error</caption>
**The $\ell_2$ norm of the transition matrix** plays opposite roles on each side. In the lower bound, $1 - \theta = 1 - \max_{\theta \in \Theta_s} \|\theta\|_2$ appears in the numerator, whereas in the upper bound, $1 - \|\theta\|_2$ appears in the denominator. It is likely that these dependencies are suboptimal, but at least they are compatible with one another: as $\|\theta\|_2 \to 1$, that is, as the VAR process becomes unstable, the lower bound tends to 0 and the upper bound to $+\infty$. This is a reflection of the fact that our proofs make heavy use of the distance between $\theta$ and the unit sphere, which means they become meaningless when $\theta$ gets too large. Note that many previous works on VAR estimation make use of the spectral radius or spectral norm of $\theta$. While Simchowitz et al. [35] is among the few exceptions, their assumptions are much stronger than ours (no transition sparsity, full observations). They also use a different estimator, and their bound contains some terms that ours does not need.

**The variances $\sigma^2$ and $\omega^2$** are involved in $\gamma_\ell$ for the lower bound, and in $\gamma_u(\theta)$ for the upper bound. In both cases, the quantity $1 + \omega^2/\sigma^2$ tells us whether the underlying process is large enough to be detected among the noise. Roughly speaking, the magnitude of $X_t$ is related to $\sigma^2$, while the magnitude of $Y_t$ is related to $\sigma^2 + \omega^2$. If the latter is significantly larger than the former, recovering $X_t$ (and thus $\theta$) is a hopeless endeavor. Lemma C.13 gives a heuristic argument linking this dependency to the asymptotic behavior of the MLE.

**Remark (Choice of error norm).** Our theorems both provide results for the error in $\ell_\infty$ operator norm. While the $\ell_2$ operator norm would have been more usual, our choice is dictated by the proof technique for the Dantzig selector, which we adapted from [17]. Note that their paper gives error bounds in $\ell_1$ operator norm, but that is because their VAR equation reads $X_t = \theta' X_{t-1}$, and we have $\|\theta'\|_1 = \|\theta\|_\infty$.

### 3.3. Possible extensions

**Heterogeneous sampling** A key assumption in the previous discussion is the uniformity of sampling marginals: every state component has the same probability $p$ of being observed. This is a useful simplification, which condenses the influence of sampling into a single parameter $p$. To handle the case of heterogeneous sampling, where each state component gets a different probability $p_d$, the proofs can be adapted by replacing $p$ with $p_{\min} = \min_d p_d$ or $p_{\max} = \max p_d$, and performing other similar changes. We thus end up with a lower bound scaling with $1/p_{\max}$ and an upper bound scaling with $1/p_{\min}$. It is likely that the gap between both can be reduced, but we leave better approaches for future research.

**Non-isotropic variances** An initial version of this manuscript handled the more general case of a non-isotropic innovation variance $\Sigma \neq \sigma^2 I$. This was done by leveraging bounds on its extremal eigenvalues $\sigma_{\min}$ and $\sigma_{\max}$, so that the lower bound for example scaled as $\frac{\sigma_{\min}^2 + \omega^2}{\sigma_{\max}^2}$. Again, more refined spectral considerations could lead to smarter bounds, but we sought to simplify exposition instead of pursuing them.

**Weak sparsity** While sparsity is a useful hypothesis for theoretical analysis, many realistic models are only weakly sparse. Adapting our results to this more generic context would only require two minor changes, which we do not include for the sake of conciseness. For the upper bound, the proof technique of Han, Lu and Liu [17] generalizes directly to the weakly
Figure 1f takes a closer look at the role of the Markov sampling parameter *b*. The white region corresponds to values of *b* for which there is no *a* such that $p = a/(a+b)$. On this logarithmic heatmap, we see regularly-spaced and nearly vertical contour lines, which is consistent with a convergence rate of 1/*p* that does not depend on *b*. We conjecture that 1/*p* is the true order of magnitude for the optimal error, and that the dependencies 1/√*pq*_u and 1/√*pq*_l from our Theorems could be refined and brought together with a more careful theoretical analysis.

# 5 CONCLUSION AND PERSPECTIVES

In this paper, we studied a partially-observed VAR process, whose latent state components are randomly projected and corrupted with noise before being observed. The temporal correlations within the sampling process are a novel feature, and combining both sources of randomness (discrete and continuous) required the use of tailored probabilistic methods. We provided upper and lower bounds for the optimal estimation error on the transition matrix, and found that these bounds roughly match. Our analysis, supported by empirical results, sheds light on the intrinsic difficulty of such statistical problems, which arise naturally when analyzing several types of network processes.

However, our study leaves many questions open for future work. On the theoretical side, bridging the gap between our bounds will probably require more sophisticated tools to capture the precise behavior of Markov sampling. Going from the uniform case, where the sampling probability equals *p* everywhere, to more realistic heterogeneous settings, is also a worthy avenue to explore. On the practical side, this linear Gaussian model may not perform well when applied to real prediction problems. Finding ways to enhance it will be necessary if we want to gather insights on complex high-dimensional dynamics, especially for graph-structured data.

# ACKNOWLEDGEMENTS

The authors would like to thank their colleague Axel Parmentier for his collaboration and careful proofreading. Clément Mantoux, Éloïse Berthier, Maxime Godin and Pierre Marion (by alphabetical order of first names) provided more support and advice than can be described in such a constrained space. We also thank Emeline Luirard for her help with a critical Lemma, and Mathieu Besançon for his last-minute look at the draft.

We are very grateful to the SNCF, especially its departments DGEX Solutions (SNCF Réseau) and Transilien (SNCF Voyageurs) for providing us with the inspiration behind this work.

# REFERENCES

S. Basu and G. Michailidis. Regularized estimation in sparse high-dimensional time series models. *Annals of Statistics*, 43(4):1535–1567, Aug. 2015. ISSN 0090-5364, 2168-8966. doi: 10.1214/15-AOS1315. URL https://projecteuclid.org/euclid.aos/1434546214.

J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. Julia: A Fresh Approach to Numerical Computing. *SIAM Review*, 59(1):65–98, Jan. 2017. ISSN 0036-1445. doi: 10.1137/141000671. URL https://epubs.siam.org/doi/10.1137/141000671.

J. P. Buonaccorsi. *Measurement Error: Models, Methods, and Applications*. CRC Press, Mar. 2010. ISBN 978-1-4200-6658-6.

E. Candes and T. Tao. The Dantzig selector: Statistical estimation when *p* is much larger than *n*. *The Annals of Statistics*, 35(6):2313–2351, Dec. 2007. ISSN 0090-5364, 2168-8966. doi: 10/b4rfq6. URL https://projecteuclid.org/journals/annals-of-statistics/volume-35/issue-6/The-Dantzig-selector--Statistical-estimation-when-p-is-much/10.1214/009053606000001523.full.

O. Cappé, É. Moulines, and T. Ryden. *Inference in Hidden Markov Models*. Springer Science & Business Media, Apr. 2006. ISBN 978-0-387-28982-3.
sparse case, as they only assume a bound on the sum of the q-th powers of all coefficients in a single column. As for the lower bound, one would need to adjust the packing in Fano's method (Section B.9) depending on the exponent q.

**Support recovery** The Dantzig selector is an implicit method for choosing a desirable subset of nonzero entries. There is a huge literature on support recovery using $l_1$-minimization programs, and many procedures motivated by a wide array of criteria have been proposed over the years. By and large, the most commonly discussed approach is the LASSO [36]. The LASSO is known to have a support recovery property under quite restrictive conditions on the design or the noise level, such as the Irrepresentability condition and the $\beta$-min condition. The Dantzig selector, on the other hand, is a generalization of the LASSO that can handle a wider range of design matrices and still guarantees interesting error bounds on the coefficients, as done in this paper. It is particularly well-suited to our problem because our constraints does not involve a design matrix *per se* but estimates of the Yule-Walker equations.

**VAR processes of higher order** Although our results only apply to state-space models based on an underlying VAR process of order 1, we could try to extend them to the more general case of VAR(K) processes. However, as we show here, this does not work. Just for this Section, suppose $X_t$ is no longer given by Equation (3), but instead satisfies:

$$X_t = \theta_1 X_{t-1} + \theta_2 X_{t-2} + \dots + \theta_K X_{t-K} + \varepsilon_t.$$

Then we can represent this as a VAR(1) process using augmented variables [27]. Observe that if we define $\tilde{X}_t = (X_t \ X_{t-1} \ \dots \ X_{t-K+1})'$ and also $\tilde{\varepsilon}_t = (\varepsilon_t \ 0 \ \dots \ 0)'$, we get

$$\tilde{X}_t = \tilde{\theta} \tilde{X}_{t-1} + \tilde{\varepsilon}_t \quad \text{with} \quad \tilde{\theta} = \begin{bmatrix} \theta_1 & \theta_2 & \cdots & \theta_{K-1} & \theta_K \\ I_D & 0 & \cdots & 0 & 0 \\ 0 & I_D & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & I_D & 0 \end{bmatrix}.$$

There are two major issues preventing generalization of our theorems. The first is the assumption bounding the spectral radius of $\theta$, which would no longer hold for $\tilde{\theta}$. The second is the change of sampling mechanism. Indeed, the process that generates $\Pi_t$ gives rise to a new distribution for $\tilde{\Pi}_t$ which is no longer part of the same family as it exhibits more internal correlations. In fact, the augmented sampling process $\tilde{\Pi}_t$ is still Markovian but with a memory of size K instead of 1. Therefore, the adaptation would require new arguments, and we leave it for future work.

# 4. Numerical illustrations

We now illustrate our results on simulated data. All experiments were performed on a Mac-Book Pro (2023) with an Apple Silicon M3 processor and 36 GiB of RAM, running under Sequoia 15.3. Our code was written in Julia [3], linear optimization problems were modeled using JuMP.jl [16] and solved with the HiGHS LP solver [20]. The code used to generate all the plots is available as a supplementary material [10] and on Zenodo [9].
K.-M. Chung, H. Lam, Z. Liu, and M. Mitzenmacher. Chernoff-Hoeffding Bounds for Markov Chains: Generalized and Simplified. In C. Dürr and T. Wilke, editors, *29th International Symposium on Theoretical Aspects of Computer Science (STACS 2012)*, volume 14 of *Leibniz International Proceedings in Informatics (LIPIcs)*, pages 124–135, Dagstuhl, Germany, 2012. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik. ISBN 978-3-939897-35-4. doi: 10.4230/LIPIcs.STACS.2012.124. URL http://drops.dagstuhl.de/opus/volltexte/2012/3437.

T. M. Cover and J. A. Thomas. *Elements of Information Theory*. John Wiley & Sons, Nov. 2012. ISBN 978-1-118-58577-1.

R. Douc, É. Moulines, and D. Stoffer. *Nonlinear Time Series: Theory, Methods and Applications with R Examples*. CRC Press, Jan. 2014. ISBN 978-1-4665-0225-3.

A. Doucet, S. Godsill, and C. Andrieu. On sequential Monte Carlo sampling methods for Bayesian filtering. *Statistics and Computing*, 10(3):197–208, July 2000. ISSN 1573-1375. doi: 10.1023/A:1008935410038. URL https://doi.org/10.1023/A:1008935410038.

D. P. Dubhashi and A. Panconesi. *Concentration of Measure for the Analysis of Randomized Algorithms*. Cambridge University Press, June 2009. ISBN 978-0-521-88427-3.

J. Duchi. Derivations for linear algebra and optimization. 2007. URL https://web.stanford.edu/~jduchi/projects/general_notes.pdf.

J. Duchi. Information Theory and Statistics. 2019. URL https://stanford.edu/class/stats311/lecture-notes.pdf.

I. Dunning, J. Huchette, and M. Lubin. JuMP: A Modeling Language for Mathematical Optimization. *SIAM Review*, 59(2):295–320, Jan. 2017. ISSN 0036-1445. doi: 10/gftshn. URL https://epubs.siam.org/doi/abs/10.1137/15M1020575.

J. Forrest, S. Vigerske, T. Ralphs, L. Hafer, J. Forrest, jpfasano, H. G. Santos, M. Saltzman, Jan-Willem, B. Kristjanson, h-i-gassmann, A. King, pobonomo, S. Brito, and to-st. Coin-or/Clp: Release releases/1.17.7. Zenodo, Jan. 2022. URL https://zenodo.org/record/5839302.

F. Han, H. Lu, and H. Liu. A Direct Estimation of High Dimensional Stationary Vector Autoregressions. *Journal of Machine Learning Research*, 16(97):3115–3150, 2015. ISSN 1533-7928. URL http://jmlr.org/papers/v16/han15a.html.

R. A. Horn and C. R. Johnson. *Topics in Matrix Analysis*. Cambridge University Press, June 1994. ISBN 978-0-521-46713-1.

R. A. Horn and C. R. Johnson. *Matrix Analysis*. Cambridge University Press, Oct. 2012. ISBN 978-1-139-78888-5.

A. Jalali and R. Willett. Missing Data in Sparse Transition Matrix Estimation for Sub-Gaussian Vector Autoregressive Processes. *arXiv:1802.09511 [cs, stat]*, Feb. 2018. URL http://arxiv.org/abs/1802.09511.

R. E. Kalman. A New Approach to Linear Filtering and Prediction Problems. *Journal of Basic Engineering*, 82(1):35–45, Mar. 1960. ISSN 0021-9223. doi: 10.1115/1.3662552. URL https://asmedigitalcollection.asme.org/fluidsengineering/article/82/1/35/397706/A-New-Approach-to-Linear-Filtering-and-Prediction.

A. B. Kock and L. Callot. Oracle inequalities for high dimensional vector autoregressions. *Journal of Econometrics*, 186(2):325–344, June 2015. ISSN 0304-4076. doi: 10.1016/j.jeconom.2015.02.013. URL http://www.sciencedirect.com/science/article/pii/S0304407615000378.

D. A. Levin and Y. Peres. *Markov Chains and Mixing Times*. American Mathematical Soc., Oct. 2017. ISBN 978-1-4704-2962-1.
## 4.1. Setup

Simulating a partially-observed VAR process with known transition matrix $\theta$ allows us to compute the estimation error $||\hat{\theta} - \theta||_{\infty}$ and study the influence of parameters such as $T$, $D$, $s$, $p$, $\sigma$ and $\omega$. The non-varying parameters are specified in the title of each plot. Each error value (marker) corresponds to the median of 10 samples, with 80% confidence intervals represented as shaded bands. For each sample, we draw a new $\theta$ with every coefficient following an independent standard Gaussian distribution, and we normalize it to satisfy $||\theta||_2 = \theta = \frac{1}{2}$. We then sample a training dataset $(X, \pi, Y)$ based on the true $\theta$ and use it to compute $\hat{\theta}$.

Some of the simulations are run in a dense estimation scenario. For those that require the sparse procedure, selecting a good regularization parameter $\lambda$ is paramount: indeed, Theorem 3.1 is only valid for a specific value of $\lambda$ (which is not known in practice, but we can hope to approximate this near-optimal choice). We tune $\lambda$ using cross-validation on a held-out validation set, for which we compare the observations $Y_t$ to the predictions $\hat{Y}_t = \hat{\theta}X_{t-1}$ (assuming the previous state had been observed). The best value for $\lambda$ is selected through a grid search to minimize the mean squared prediction error.

Axes are displayed with logarithmic scaling, in order to highlight the exponent of the dependencies. When a straight line is added to a scatter plot, it is the result of a Theil-Sen regression [33] applied to the points of the same color: its slope is given in the legend.

## 4.2. Results

Figure 1 validates the dependency of the error in $1/p$, at least for large enough values of $p$. Judging by Remarks 3.1 and 3.1, the smaller $p$ becomes, the larger $T$ must be for the asymptotic scaling behavior to kick in.

Figure 2 shows that the error also decreases with $1/\sqrt{T}$, for various values of $p$. Note the regular intervals between the colored curves, which mirror the logarithmic intervals between

<figure><img src="image_1.png" /><figcaption>Fig 1. Influence of $\omega$ on the estimation error</figcaption></figure>
R. J. A. Little and D. B. Rubin. *Statistical Analysis with Missing Data*. John Wiley & Sons, Apr. 2019. ISBN 978-0-470-52679-8.

P.-L. Loh and M. J. Wainwright. High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity. *The Annals of Statistics*, 40(3):1637–1664, June 2012. ISSN 0090-5364, 2168-8966. doi: 10/ggx5bj. URL https://projecteuclid.org/journals/annals-of-statistics/volume-40/issue-3/High-dimensional-regression-with-noisy-and-missing-data--Provable/10.1214/12-AOS1018.full.

H. Lütkepohl. *New Introduction to Multiple Time Series Analysis*. Springer Science & Business Media, Dec. 2005. ISBN 978-3-540-27752-1.

L. Malagò and G. Pistone. Information Geometry of the Gaussian Distribution in View of Stochastic Optimization. In *Proceedings of the 2015 ACM Conference on Foundations of Genetic Algorithms XIII*, FOGA '15, pages 150–162, Aberystwyth, United Kingdom, Jan. 2015. Association for Computing Machinery. ISBN 978-1-4503-3434-1. doi: 10.1145/2725494.2725510. URL https://doi.org/10.1145/2725494.2725510.

I. Melnyk and A. Banerjee. Estimating Structured Vector Autoregressive Models. In *International Conference on Machine Learning*, pages 830–839, June 2016. URL http://proceedings.mlr.press/v48/melnyk16.html.

K. B. Petersen and M. S. Pedersen. The matrix cookbook. Technical report, 2012. URL https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf.

M. Rao, T. Javidi, Y. C. Eldar, and A. Goldsmith. Estimation in autoregressive processes with partial observations. In *2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, pages 4212–4216, Mar. 2017a. doi: 10.1109/ICASSP.2017.7952950.

M. Rao, T. Javidi, Y. C. Eldar, and A. Goldsmith. Fundamental estimation limits in autoregressive processes with compressive measurements. In *2017 IEEE International Symposium on Information Theory (ISIT)*, pages 2895–2899, June 2017b. doi: 10.1109/ISIT.2017.8007059.

P. K. Sen. Estimates of the Regression Coefficient Based on Kendall’s Tau. *Journal of the American Statistical Association*, 63(324):1379–1389, Dec. 1968. ISSN 0162-1459. doi: 10/gfxz87. URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1968.10480934.

R. H. Shumway and D. S. Stoffer. An Approach to Time Series Smoothing and Forecasting Using the EM Algorithm. *Journal of Time Series Analysis*, 3(4):253–264, 1982. ISSN 1467-9892. doi: 10.1111/j.1467-9892.1982.tb00349.x. URL http://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9892.1982.tb00349.x.

R. Tibshirani. Regression Shrinkage and Selection via the Lasso. *Journal of the Royal Statistical Society. Series B (Methodological)*, 58(1):267–288, 1996. ISSN 0035-9246. doi: 10.1111/j.2517-6161.1996.tb02080.x. URL https://www.jstor.org/stable/2346178.

A. B. Tsybakov. *Introduction to Nonparametric Estimation*. Springer Science & Business Media, Oct. 2008. ISBN 978-0-387-79052-7.

R. Vershynin. *High-Dimensional Probability: An Introduction with Applications in Data Science*. Cambridge University Press, Sept. 2018. ISBN 978-1-108-24454-1.

M. J. Wainwright. *High-Dimensional Statistics: A Non-Asymptotic Viewpoint*. Cambridge University Press, Feb. 2019. ISBN 978-1-108-49802-9.
<figure><img src="image_2.png" /><figcaption>Fig 2. Influence of T on the estimation error</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig 3. Influence of $\omega$ on the estimation error</figcaption></figure>

the corresponding values of $p$. This dependence on $T$ only holds because the sampling probability $p$ remains constant. If instead we had a limited observation budget but an increasing temporal precision, we would have $p \propto 1/T$, which would render our analysis invalid.

Figure 3 exhibits three clearly identifiable regimes with respect to the noise variance. In the first one, corresponding to $\omega/\sigma \ll 1$, the error remains small and constant. Then, the error increases when $\omega/\sigma \approx 1$. In the third phase, corresponding to $\omega/\sigma \gg 1$, the error remains high and volatile. This is consistent with the theoretical dependency in $1 + \omega^2/\sigma^2$.

Figure 4 explores the role of $D$ in the dense case, or $s$ in the sparse case. Unlike the other figures, it shows a divergence from the theoretical bound. While the dense error does scale with $D$, the sparse error does not have the same dependency on $s$ (the slope is < 1). Furthermore, if we plot its dependency as a function of both $s$ and $D$ in a heatmap, we find that both parameters have a non-negligible effect. This is in stark contrast to our theoretical
# A PROOF OF THE ESTIMATOR'S CONVERGENCE RATE

Here we present the detailed proof of Theorem 1.

## A.1 OVERVIEW

The main steps of the argument are the following:

1. Prove the Yule-Walker Equation (6) and deduce an expression for the covariance matrix of X (Lemmas 1 and 2).

2. Justify the formula of Equation (9) for $\hat{\Gamma}_h$ by showing that it defines an unbiased estimator of $\Gamma_h$ (Lemmas 3 and 4).

3. Fixing two indices $d_1$ and $d_2$, rewrite $(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}$ using quadratic forms $g'_a \Psi'_a L \Psi_b g_b$ of standard Gaussian vectors (Lemma 5).

4. Control the deviation of the matrix L using discrete concentration inequalities (Lemmas 6, 7, 8, 9 and 10).

5. Apply a conditional version of the Hanson-Wright inequality (Lemma 37) to the quadratic forms $g'_a \Psi'_a L \Psi_b g_b$ (Lemma 11).

6. Obtain a high-probability control on $\|\hat{\Gamma}_h - \Gamma_h\|_{\max}$ with a union bound (Lemma 13).

7. Deduce the error of $\hat{\theta}$ from the error of $\hat{\Gamma}_{h_0}$ and $\hat{\Gamma}_{h_0+1}$ by drawing inspiration from Han et al. [2015] (Lemmas 14 and 15).

## A.2 COVARIANCE MATRICES

The Yule-Walker equation is a direct consequence of the VAR recursion, as can be seen from this Lemma.

**Lemma 1 (VAR covariance matrices).** The autocovariance matrices of the stationary VAR process defined by Equation (3) have the following expressions:

$$ \Gamma_0(\theta) = \mathrm{Cov}_\theta[X_t] = \sum_{k=0}^{\infty} \theta^k \Sigma \theta'^k $$

$$ \Gamma_h(\theta) = \mathrm{Cov}_\theta[X_{t+h}, X_t] = \theta^h \Gamma_0(\theta). $$

*Proof.* We start by noting that according to Equation (3), the stacked vector $X = (X_t)_{t \in [T]}$ follows a TD-dimensional centered multivariate Gaussian distribution. The covariance matrix of $X_t$ can be deduced from the recursion:

$$ \Gamma_0(\theta) = \mathrm{Cov}_\theta[X_t] = \mathrm{Cov}_\theta[\theta X_{t-1} + \varepsilon_t] = \theta \mathrm{Cov}_\theta[X_{t-1}]\theta' + \Sigma = \theta \Gamma_0(\theta)\theta' + \Sigma. $$

There is a unique stationary solution:

$$ \Gamma_0(\theta) = \sum_{k=0}^{\infty} \theta^k \Sigma \theta'^k. $$

The covariance matrix between $X_{t+h}$ and $X_t$ is obtained similarly:

$$ \Gamma_h(\theta) = \mathrm{Cov}_\theta[X_{t+h}, X_t] = \mathbb{E}[X_{t+h}X_t'] = \mathbb{E}[(\theta X_{t+h-1} + \varepsilon_{t+h})X_t'] \\ = \theta \mathrm{Cov}_\theta[X_{t+h-1}, X_t] = \theta^h \mathrm{Cov}_\theta[X_t, X_t] = \theta^h \Gamma_0(\theta). $$

And $\mathrm{Cov}_\theta[X_t, X_{t+h}] = \mathrm{Cov}_\theta[X_{t+h}, X_t]'$. In other words, we just proved that

$$ \mathrm{Cov}_\theta[X] = \begin{bmatrix} \Gamma_0(\theta) & \Gamma_0(\theta)\theta'^1 & \Gamma_0(\theta)\theta'^2 & \cdots & \Gamma_0(\theta)\theta'^{T-1} \\ \theta^1\Gamma_0(\theta) & \Gamma_0(\theta) & \Gamma_0(\theta)\theta'^1 & & \Gamma_0(\theta) \\ \theta^2\Gamma_0(\theta) & \theta^1\Gamma_0(\theta) & \Gamma_0(\theta) & & \Gamma_0(\theta) \\ \vdots & & & & \vdots \\ \theta^{T-1}\Gamma_0(\theta) & & & & \Gamma_0(\theta) \end{bmatrix} $$
<figure><img src="image_4.png" /><figcaption>Fig 4. Influence of *D* and *s* on the estimation error</figcaption></figure>

expectations, according to which *D* only appears in a logarithmic term, whereas *s* is the leading influence on the error.

We do not yet fully make sense of this behavior. One possible interpretation is that the function $\gamma_u(\theta)$ also depends on the sparsity level in complicated ways through $\theta$, especially since the real values are renormalized to satisfy $\|\theta\|_2 = \frac{1}{2}$. Thus, there may be an effect of *s* and *D* hidden inside $\gamma_u(\theta)$, which is hard to disentangle from the predicted scaling when $\|\theta\|_2$ is close to 1.

# 5. Conclusion and perspectives

In this paper, we studied a partially-observed VAR process, whose latent state components are randomly projected and corrupted with noise before being observed. The partial observations
□

The following result will come in handy later.

**Lemma 2 (Norm of $\Gamma_0(\theta)$).** *The covariance matrix $\Gamma_0(\theta)$ satisfies*

$$ \|\Gamma_0(\theta)\|_2 \le \frac{\sigma_{\max}^2}{1 - \vartheta^2} $$

*Proof.* By Lemma 1,

$$ \|\Gamma_0(\theta)\|_2 \le \sum_{k=0}^{\infty} \|\theta^k \Sigma \theta'^k\|_2 \le \sum_{k=0}^{\infty} \|\theta\|_2^k \|\Sigma\|_2 \|\theta\|_2^k = \frac{\|\Sigma\|_2}{1 - \|\theta\|_2^2} \le \frac{\sigma_{\max}^2}{1 - \vartheta^2}. $$

□

## A.3 CONSTRUCTION OF THE COVARIANCE ESTIMATOR

Now we justify the construction of our covariance estimator. Let $h_0 = 0$: for most of the proof, we fix a lag value $h \in \{h_0, h_0 + 1\} = \{0, 1\}$.

**Lemma 3 (Bias of the covariance estimator).** *The estimator $\hat{\Gamma}_h$ given by Equation (9) for the covariance matrix $\Gamma_h$ is unbiased.*

*Proof.* First, let us remember that since $\Pi_t = \text{diag}(\pi_t)$ is diagonal and binary, we also have $\Pi_t^\dagger = \Pi_t' = \Pi_t$. By Equation (5),

$$ \begin{aligned} (\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_t^{\dagger} Y_t)' &= \Pi_{t+h}^{\dagger} (\Pi_{t+h} X_{t+h} + \eta_{t+h}) (X_t' \Pi_t' + \eta_t') \Pi_t^{\dagger'} \\ &= \text{diag}(\pi_{t+h}) (X_{t+h} X_t' + X_{t+h} \eta_t' + \eta_{t+h} \eta_t' + \eta_{t+h} \eta_t') \text{diag}(\pi_t). \end{aligned} \quad (15) $$

Taking the conditional expectation and removing the cross-product terms (by independence of $X$ and $\Pi$), we get:

$$ \mathbb{E}[(\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_t^{\dagger} Y_t)' | \Pi] = \text{diag}(\pi_{t+h}) (\mathbb{E}[X_{t+h} X_t'] + \mathbb{E}[\eta_{t+h} \eta_t']) \text{diag}(\pi_t). $$

Since $\mathbb{E}[X_{t+h} X_t'] = \Gamma_h$ and $\mathbb{E}[\eta_{t+h} \eta_t] = \mathbf{1}_{\{h=0\}} \omega^2 I$, we are left with:

$$ \mathbb{E}[(\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_t^{\dagger} Y_t)' | \Pi] = (\pi_{t+h} \pi_t') \odot \Gamma_h + \mathbf{1}_{\{h=0\}} \omega^2 \text{diag}(\pi_t). $$

where $\odot$ is the elementwise Hadamard product. We now take the expectation w.r.t. $\Pi$:

$$ \mathbb{E}[(\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_t^{\dagger} Y_t)'] = \mathbb{E}[\pi_{t+h} \pi_t'] \odot \Gamma_h + \mathbf{1}_{\{h=0\}} \omega^2 \mathbb{E}[\text{diag}(\pi_t)]. $$

Dividing elementwise by the scaling matrix $S(h) = \mathbb{E}[\pi_{t+h} \pi_t']$, we get

$$ \begin{aligned} \mathbb{E} \left[ \frac{1}{\mathbb{E}[\pi_{t+h} \pi_t'] \odot (\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_t^{\dagger} Y_t)'} \right] &= \Gamma_h + \mathbf{1}_{\{h=0\}} \omega^2 \mathbb{E}[\text{diag}(\pi_t)] \odot \frac{1}{\mathbb{E}[\pi_t \pi_t']'} \\ &= \Gamma_h + \mathbf{1}_{\{h=0\}} \omega^2 \text{diag} \left( \frac{\mathbb{E}[\pi_t]}{\mathbb{E}[\pi_t^2]} \right) \\ &= \Gamma_h + \mathbf{1}_{\{h=0\}} \omega^2 I \end{aligned} $$

which shows that our estimator

$$ \hat{\Gamma}_h = \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{1}{\mathbb{E}[\pi_{t+h} \pi_t'] \odot (\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_t^{\dagger} Y_t)'} - \mathbf{1}_{\{h=0\}} \omega^2 I $$

is unbiased.

□
are a novel feature, and combining both sources of randomness (discrete and continuous) required the use of tailored probabilistic methods. We provided upper and lower bounds for the optimal estimation error on the transition matrix, and found that these bounds roughly match. Our analysis, supported by empirical results, sheds light on the intrinsic difficulty of such statistical problems, which arise naturally when analyzing several types of network processes.

However, our study leaves many questions open for future work. On the theoretical side, richer behaviors may be uncovered by allowing for non-isotropic variances, as well as non-uniform or correlated sampling. The discrepancy between numerical experiments and theory with respect to the role of the sparsity level also deserves further investigation. On the practical side, this linear Gaussian model may not perform well when applied to real prediction problems. Finding ways to enhance it will be necessary if we want to gather insights on complex high-dimensional dynamics, especially for graph-structured data.

# Appendix A: Proof of the estimator's convergence rate

Here we present the detailed proof of Theorem 3.1.

## A.1. Overview

The main steps of the argument are the following:

1. Prove the Yule-Walker Equation (6) and deduce an expression for the covariance matrix of $X$ (Lemmas A.1 and A.2).

2. Justify the formula of Equation (9) for $\hat{\Gamma}_h$ by showing that it defines an unbiased estimator of $\Gamma_h$ (Lemmas A.3 and A.4).

3. Fixing two indices $d_1$ and $d_2$, rewrite $(\hat{\Gamma}_h - \Gamma_h)_{d_1,d_2}$ thanks to quadratic forms $g'_a \Psi'_a L \Psi_b g_b$ of standard Gaussian vectors (Lemma A.5).

4. Control the deviation of the matrix $L$ using discrete concentration inequalities (Lemmas A.6, A.7, A.8, A.9 and A.10).

5. Apply a conditional version of the Hanson-Wright inequality (Lemma C.12) to the quadratic forms $g'_a \Psi'_a L \Psi_b g_b$ (Lemma A.11).

6. Obtain a high-probability control on $\|\hat{\Gamma}_h - \Gamma_h\|_{\max}$ (where $\|\cdot\|_{\max}$ is the maximum absolute value of the coefficients) with a union bound (Lemma A.13).

7. Deduce the error of $\hat{\theta}$ from the error of $\hat{\Gamma}_{h_0}$ and $\hat{\Gamma}_{h_0+1}$ by drawing inspiration from Han, Lu and Liu [17] (Lemmas A.14 and A.15).

## A.2. Covariance matrices

The Yule-Walker equation is a direct consequence of the VAR recursion, as can be seen from this Lemma.

**Lemma A.1 (VAR covariance matrices).** The autocovariance matrices of the stationary VAR process defined by Equation (3) have the following expressions:

$$ \Gamma_0(\theta) = \mathrm{Cov}_{\theta}[X_t] = \sigma^2 \sum_{k=0}^{\infty} \theta^k \theta'^k \quad \text{and} \quad \Gamma_h(\theta) = \mathrm{Cov}_{\theta}[X_{t+h}, X_t] = \theta^h \Gamma_0(\theta). $$
Note that since the process $(\Pi_t)$ is stationary, the coefficients of $S(h)$ do not depend on $t$. They are computed in the next Lemma.

**Lemma 4.** The second-order moments of $\pi$ are given by

$$S(h)_{d_1, d_2} = \mathbb{E}[\pi_{t+h, d_1} \pi_{t, d_2}] = \begin{cases} p^2 & \text{if } d_1 \neq d_2 \\ p & \text{if } d_1 = d_2 \text{ and } h = 0 \\ p^2 + p(1-p)(1-a-b)^h & \text{if } d_1 = d_2 \text{ and } h \ge 1 \end{cases}$$

In particular, every coefficient of the scaling matrix $S(h)$ is lower-bounded by

$$\min_{d_1, d_2, h} S(h)_{d_1, d_2} = \min\{p^2, p(1-b)\} = pq_u \quad \text{where} \quad q_u = \min\{p, 1-b\}.$$

*Proof.* Let $i = (t+h, d_1)$ and $j = (t, d_2)$ be two indices in $[T] \times [D]$. We have $\mathbb{E}[\pi_i] = \mathbb{E}[\pi_i^2] = p$. If $d_1 \neq d_2$, then the variables $\pi_i$ and $\pi_j$ belong to independent Markov chains, and thus $\mathbb{E}[\pi_i \pi_j] = p^2$. Otherwise, we have $i = (t+h, d)$ and $j = (t, d)$, which means these two variables are part of the same Markov chain. Stationarity yields

$$\mathbb{E}[\pi_i \pi_j] = \mathbb{P}(\pi_{t,d} = 1) \times \mathbb{P}(\pi_{t+h,d} = 1 | \pi_{t,d} = 1) = p(\mathcal{T}^h)_{11}.$$

When diagonalizing the transition matrix $\mathcal{T}$, we see that the bottom-right coefficient of $\mathcal{T}^h$ is

$$(\mathcal{T}^h)_{11} = \frac{a + b(1 - a - b)^h}{a + b} = p + (1 - p)(1 - a - b)^h.$$

Plugging this in, we get

$$\mathbb{E}[\pi_i \pi_j] = p^2 + p(1-p)(1-a-b)^h.$$

Among all the possible values of $S(h)_{d_1, d_2}$, the smallest one is $p^2$ if $1 - a - b \ge 0$, and $p^2 + p(1 - p)(1 - a - b)$ otherwise. But since

$$\begin{aligned} p + (1 - p)(1 - a - b) &= \frac{a}{a + b} + \frac{b}{a + b}(1 - a - b) \\ &= \frac{a + b - ab - b^2}{a + b} = \frac{a(1 - b) + b(1 - b)}{a + b} \\ &= 1 - b, \end{aligned}$$

we conclude

$$\min_{d_1, d_2, h} S(h)_{d_1, d_2} = \min\{p^2, p^2 + p(1-p)(1-a-b)\} = \min\{p^2, p(1-b)\}.$$

□

## A.4 GAUSSIAN CONCENTRATION, EPISODE 1

From now on, we will study the concentration of $\hat{\Gamma}_h$, coefficient by coefficient. Let us fix two indices $d_1$ and $d_2$: our goal is to control the deviation of $(\hat{\Gamma}_h)_{d_1, d_2}$ around its mean.

**Lemma 5 (Deviation of $(\hat{\Gamma}_h)_{d_1, d_2}$).** The deviation probability for $(\hat{\Gamma}_h)_{d_1, d_2}$ can be decomposed as follows:

$$\begin{aligned} \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \ge u) \le & \mathbb{P}(|g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon - \mathbb{E}[g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon]| \ge u/4) \\ & + \mathbb{P}(|g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon - \mathbb{E}[g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon]| \ge u/4) \\ & + \mathbb{P}(|g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta - \mathbb{E}[g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta]| \ge u/4) \\ & + \mathbb{P}(|g'_\eta \Psi'_\eta L \Psi_\eta g_\eta - \mathbb{E}[g'_\eta \Psi'_\eta L \Psi_\eta g_\eta]| \ge u/4) \end{aligned}$$

where the random matrix $L$ is defined in Equation (16), $\Psi_\epsilon$ and $\Psi_\eta$ are defined in Equation (18), and $g_\epsilon$ and $g_\eta$ are standard Gaussian vectors.
**Proof.** We start by noting that according to Equation (3), the stacked vector $X = (X_t)_{t \in [T]}$ follows a TD-dimensional centered multivariate Gaussian distribution. The covariance matrix of $X_t$ can be deduced from the recursion:

$$ \Gamma_0(\theta) = \mathrm{Cov}_{\theta}[X_t] = \mathrm{Cov}_{\theta}[\theta X_{t-1} + \varepsilon_t] = \theta \mathrm{Cov}_{\theta}[X_{t-1}]\theta' + \sigma^2 I = \theta \Gamma_0(\theta)\theta' + \sigma^2 I. $$

There is a unique stationary solution $\Gamma_0(\theta) = \sigma^2 \sum_{k=0}^{\infty} \theta^k \theta'^k$. The covariance matrix between $X_{t+h}$ and $X_t$ is obtained similarly:

$$ \begin{aligned} \Gamma_h(\theta) &= \mathrm{Cov}_{\theta}[X_{t+h}, X_t] = \mathbb{E}[X_{t+h} X_t'] = \mathbb{E}[(\theta X_{t+h-1} + \varepsilon_{t+h}) X_t'] \\ &= \theta \mathrm{Cov}_{\theta}[X_{t+h-1}, X_t] = \theta^h \mathrm{Cov}_{\theta}[X_t, X_t] = \theta^h \Gamma_0(\theta). \end{aligned} $$

And $\mathrm{Cov}_{\theta}[X_t, X_{t+h}] = \mathrm{Cov}_{\theta}[X_{t+h}, X_t]'$. In other words, we just proved that

$$ \mathrm{Cov}_{\theta}[X] = \left[ \begin{array}{ccccc} \Gamma_0(\theta) & \Gamma_0(\theta)\theta'^1 & \Gamma_0(\theta)\theta'^2 & \cdots & \Gamma_0(\theta)\theta'^{T-1} \\ \theta^1\Gamma_0(\theta) & \Gamma_0(\theta) & \Gamma_0(\theta)\theta'^1 & & \\ \theta^2\Gamma_0(\theta) & \theta^1\Gamma_0(\theta) & \Gamma_0(\theta) & & \\ \vdots & & \ddots & & \\ \theta^{T-1}\Gamma_0(\theta) & & & & \Gamma_0(\theta) \end{array} \right] \quad \square $$

The following result will come in handy later.

**Lemma A.2 (Norm of $\Gamma_0(\theta)$).** The covariance matrix $\Gamma_0(\theta)$ satisfies $\|\Gamma_0(\theta)\|_2 \le \frac{\sigma^2}{1-\vartheta^2}$.

*Proof.* By Lemma A.1,

$$ \|\Gamma_0(\theta)\|_2 \le \sigma^2 \sum_{k=0}^{\infty} \|\theta^k \theta'^k\|_2 \le \sigma^2 \sum_{k=0}^{\infty} \|\theta\|_2^k \|\theta\|_2^k = \frac{\sigma^2}{1 - \|\theta\|_2^2} \le \frac{\sigma^2}{1 - \vartheta^2}. \quad \square $$

### A.3. Construction of the covariance estimator

Now we justify the construction of our covariance estimator. Let $h_0 = 0$: for most of the proof, we fix a lag value $h \in \{h_0, h_0+1\} = \{0, 1\}$.

**Lemma A.3 (Bias of the covariance estimator).** The estimator $\hat{\Gamma}_h$ given by Equation (9) for the covariance matrix $\Gamma_h$ is unbiased.

*Proof.* First, let us remember that since $\Pi_t = \mathrm{diag}(\pi_t)$ is diagonal and binary, we also have $\Pi_t^\dagger = \Pi_t' = \Pi_t$. By Equation (5),

$$ \begin{aligned} (\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' &= \Pi_{t+h}^{\dagger} (\Pi_{t+h} X_{t+h} + \eta_{t+h}) (X_t' \Pi_t' + \eta_t') \Pi_t^{\dagger} \\ &= \mathrm{diag}(\pi_{t+h}) (X_{t+h} X_t' + X_{t+h} \eta_t' + \eta_{t+h} X_t' + \eta_{t+h} \eta_t') \mathrm{diag}(\pi_t). \end{aligned} \quad (15) $$

Taking the conditional expectation and removing the cross-product terms (by independence of $X$ and $\Pi$), we get:

$$ \mathbb{E}[(\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' | \Pi] = \mathrm{diag}(\pi_{t+h}) (\mathbb{E}[X_{t+h} X_t'] + \mathbb{E}[\eta_{t+h} \eta_t']) \mathrm{diag}(\pi_t). $$

Since $\mathbb{E}[X_{t+h} X_t'] = \Gamma_h$ and $\mathbb{E}[\eta_{t+h} \eta_t] = \mathbf{1}_{\{h=0\}} \omega^2 I$, we are left with:

$$ \mathbb{E}[(\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' | \Pi] = (\pi_{t+h} \pi_t') \odot \Gamma_h + \mathbf{1}_{\{h=0\}} \omega^2 \mathrm{diag}(\pi_t) $$
Proof. We denote by $e_d$ the basis vector filled with zeros except for a 1 in position $d$. By Equation (9),

$$ (\hat{\Gamma}_h + \mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} = \frac{1}{T-h} \sum_{t=1}^{T-h} \left( \frac{1}{S(h)} \odot (\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' \right)_{d_1, d_2} \\ = \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{1}{S(h)_{d_1, d_2}} e'_{d_1} (\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' e_{d_2} \\ = \frac{1}{T-h} \sum_{t=1}^{T-h} \mathrm{Tr} \left[ \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} (\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' \right] $$

Equation (15) allows us to rewrite $(\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_t^{\dagger} Y_t)'$:

$$ (\hat{\Gamma}_h + \mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} = \frac{1}{T-h} \sum_{t=1}^{T-h} X'_t \mathrm{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \mathrm{diag}(\pi_{t+h}) X_{t+h} \\ + \frac{1}{T-h} \sum_{t=1}^{T-h} \eta'_t \mathrm{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \mathrm{diag}(\pi_{t+h}) X_{t+h} \\ + \frac{1}{T-h} \sum_{t=1}^{T-h} X'_t \mathrm{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \mathrm{diag}(\pi_{t+h}) \eta_{t+h} \\ + \frac{1}{T-h} \sum_{t=1}^{T-h} \eta'_t \mathrm{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \mathrm{diag}(\pi_{t+h}) \eta_{t+h} $$

Let us denote by $P_t$ the projection of $\mathbb{R}^{TD}$ keeping only the components associated with time $t$, i.e. such that $X_t = P_t X$ and $\eta_t = P_t \eta$. We recognize the following matrix $L$ in all four lines of the expression above:

$$ L = \frac{1}{T-h} \sum_{t=1}^{T-h} P'_t \mathrm{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \mathrm{diag}(\pi_{t+h}) P_{t+h} \quad (16) \\ = \frac{1}{T-h} \sum_{t=1}^{T-h} P'_t \frac{\pi_{t+h, d_1} \pi_{t, d_2} e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} P_{t+h} $$

This leads to:

$$ (\hat{\Gamma}_h + \mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} = X' LX + \eta' LX + X' L \eta + \eta' L \eta $$

Since $X$ and $\eta$ both follow centered multivariate Gaussian distributions, we can express them as linear combinations of standard Gaussian vectors $g_\epsilon$ and $g_\eta$ of dimension $TD$ (indexed by the source of randomness):

$$ X = \Psi_{\epsilon} g_{\epsilon} \quad \text{and} \quad \eta = \Psi_{\eta} g_{\eta} \quad (17) $$

where $\Psi_\epsilon$ and $\Psi_\eta$ are the square roots of the respective covariance matrices

$$ \Psi_{\epsilon} = \mathrm{Cov}[X]^{1/2} \quad \text{and} \quad \Psi_{\eta} = \mathrm{Cov}[\eta]^{1/2} = \omega I. \quad (18) $$

We substitute $X$ and $\eta$ to get:

$$ (\hat{\Gamma}_h + \mathbf{1}_{\{h=0\}} \omega^2 I)_{d_2, d_1} = g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon + g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon + g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta + g'_\eta \Psi'_\eta L \Psi_\eta g_\eta, $$

which implies

$$ (\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2} = g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon - \mathrm{E}[g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon] \\ + g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon - \mathrm{E}[g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon] \\ + g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta - \mathrm{E}[g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta] \\ + g'_\eta \Psi'_\eta L \Psi_\eta g_\eta - \mathrm{E}[g'_\eta \Psi'_\eta L \Psi_\eta g_\eta]. $$

The union bound gives us the expected result. $\square$
where $\odot$ is the elementwise Hadamard product. We now take the expectation w.r.t. $\Pi$:

$$ \mathbb{E}[(\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_{t}^{\dagger} Y_{t})'] = \mathbb{E}[\pi_{t+h} \pi'_{t}] \odot \Gamma_{h} + \mathbf{1}_{\{h=0\}} \omega^{2} \mathbb{E}[\mathrm{diag}(\pi_{t})]. $$

Dividing elementwise by the scaling matrix $S(h) = \mathbb{E}[\pi_{t+h} \pi'_{t}]$, we get

$$ \begin{aligned} \mathbb{E} \left[ \frac{1}{\mathbb{E}[\pi_{t+h} \pi'_{t}]} \odot (\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_{t}^{\dagger} Y_{t})' \right] &= \Gamma_{h} + \mathbf{1}_{\{h=0\}} \omega^{2} \mathbb{E}[\mathrm{diag}(\pi_{t})] \odot \frac{1}{\mathbb{E}[\pi_{t} \pi'_{t}]} \\ &= \Gamma_{h} + \mathbf{1}_{\{h=0\}} \omega^{2} \mathrm{diag} \left( \frac{\mathbb{E}[\pi_{t}]}{\mathbb{E}[\pi_{t}^{2}]} \right) \\ &= \Gamma_{h} + \mathbf{1}_{\{h=0\}} \omega^{2} I \end{aligned} $$

which shows that our estimator

$$ \hat{\Gamma}_h = \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{1}{\mathbb{E}[\pi_{t+h} \pi'_{t}]} \odot (\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_{t}^{\dagger} Y_{t})' - \mathbf{1}_{\{h=0\}} \omega^2 I $$

is unbiased.

<figure><img src="image_5.png" /></figure>

Note that since the process $(\Pi_t)$ is stationary, the coefficients of $S(h)$ do not depend on $t$. They are computed in the next Lemma.

**Lemma A.4.** The second-order moments $S(h)_{d_1, d_2} = \mathbb{E}[\pi_{t+h, d_1} \pi_{t, d_2}]$ of $\pi$ are given by $S(h)_{d_1, d_2} = p$ if $d_1 = d_2$ and $h = 0$ and $S(h)_{d_1, d_2} = p^2$ otherwise. In particular, every coefficient of the scaling matrix $S(h)$ is lower-bounded by $\min_{d_1, d_2, h} S(h)_{d_1, d_2} = p^2$.

*Proof.* Let $i = (t+h, d_1)$ and $j = (t, d_2)$ be two indices in $[T] \times [D]$. We have $\mathbb{E}[\pi_i] = \mathbb{E}[\pi_t^2] = p$. If $d_1 \neq d_2$ or $t+h \neq t$, then the variables $\pi_i$ and $\pi_j$ are independent, and thus $\mathbb{E}[\pi_t \pi_j] = p^2$. $\square$

### A.4. Gaussian concentration, episode 1

From now on, we will study the concentration of $\hat{\Gamma}_h$, coefficient by coefficient. Let us fix two indices $d_1$ and $d_2$: our goal is to control the deviation of $(\hat{\Gamma}_h)_{d_1, d_2}$ around its mean.

**Lemma A.5 (Deviation of $(\hat{\Gamma}_h)_{d_1, d_2}$).** The deviation probability for $(\hat{\Gamma}_h)_{d_1, d_2}$ can be decomposed as follows:

$$ \begin{aligned} \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \ge u) \le & \mathbb{P}(|g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon - \mathbb{E}[g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon]| \ge u/4) \\ & + \mathbb{P}(|g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon - \mathbb{E}[g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon]| \ge u/4) \\ & + \mathbb{P}(|g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta - \mathbb{E}[g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta]| \ge u/4) \\ & + \mathbb{P}(|g'_\eta \Psi'_\eta L \Psi_\eta g_\eta - \mathbb{E}[g'_\eta \Psi'_\eta L \Psi_\eta g_\eta]| \ge u/4) \end{aligned} $$

where the random matrix $L$ is defined in Equation (16), $\Psi_\epsilon$ and $\Psi_\eta$ are defined in Equation (17), and $g_\epsilon$ and $g_\eta$ are standard Gaussian vectors.

*Proof.* We denote by $e_d$ the basis vector filled with zeros except for a 1 in position $d$. By Equation (9),

$$ (\hat{\Gamma}_h + \mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} = \frac{1}{T-h} \sum_{t=1}^{T-h} \left( \frac{1}{S(h)} \odot (\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_{t}^{\dagger} Y_{t})' \right)_{d_1, d_2} $$
Now, our goal is to apply a Gaussian concentration inequality to these deviation probabilities. However, since $L$ is generated by the discrete sampling process $\pi$, it is random, and so are the products $\Psi'_a L \Psi_b$ (where $a, b \in \{\varepsilon, \eta\}$). We thus need a conditional version of the Hanson-Wright inequality (Lemma 37), in which the following random variables will come into play:

* The spectral norm $\|\Psi'_a L \Psi_b\|_2$

* The Frobenius norm $\|\Psi'_a L \Psi_b\|_F^2$

* The shifted trace $\mathrm{Tr}(\Psi'_a L \Psi_b - \mathbb{E}[\Psi'_a L \Psi_b])$

## A.5 INTERLUDE: DISCRETE CONCENTRATION

We exploit discrete concentration results to bound the deviations of the three quantities we just mentioned, starting with the norms.

**Lemma 6 (Norm reformulation for $L$).** The spectral and Frobenius norms of $L$ are given by

$$ \|L\|_2 = \frac{\max_{t \in [T-h]} \pi_{t+h, d_1} \pi_{t, d_2}}{(T-h)S(h)_{d_1, d_2}} \quad \text{and} \quad \|L\|_F^2 = \frac{1}{(T-h)^2 S(h)_{d_1, d_2}} \sum_{t=1}^{T-h} \pi_{t+h, d_1} \pi_{t, d_2}. $$

*Proof.* We first notice that $L$ has a block-superdiagonal structure of rank $h$:

$$ L = \frac{1}{T-h} \sum_{t=1}^{T-h} P'_t L_{[t,t+h]} P_{t+h} \quad \text{with} \quad L_{[t,t+h]} = \frac{\pi_{t+h,d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} e_{d_2} e'_{d_1}. \qquad (19) $$

The spectral and Frobenius norms of such a matrix can easily be deduced from those of its blocks. Since $\|e_{d_2} e'_{d_1}\|_2 = \|e_{d_2} e'_{d_1}\|_F = 1$ and the $\pi_t$ are binary-valued, this leads to the following formulas:

$$ \|L\|_2 = \frac{1}{T-h} \max_{t \in [T-h]} \|L_{[t,t+h]}\|_2 = \frac{1}{(T-h)S(h)_{d_1,d_2}} \max_{t \in [T-h]} \pi_{t+h,d_1} \pi_{t, d_2} $$

$$ \|L\|_F^2 = \frac{1}{(T-h)^2} \sum_{t=1}^{T-h} \|L_{[t,t+h]}\|_F^2 = \frac{1}{(T-h)^2 S(h)_{d_1,d_2}^2} \sum_{t=1}^{T-h} \pi_{t+h,d_1} \pi_{t, d_2}. $$

□

We can bound the spectral norm for free.

**Lemma 7 (Spectral norm bound for $L$).** With probability 1, the spectral norm $\|L\|_2$ satisfies

$$ \|L\|_2 \le \frac{c}{Tpq_u} $$

*Proof.* Note that $S(h)_{d_1, d_2} \ge pq_u$, and since $h \in \{0, 1\}$, we can state that $T-h \ge cT$. By Lemma 6, we deduce

$$ \|L\|_2 = \frac{\max_{t \in [T-h]} \pi_{t+h, d_1} \pi_{t, d_2}}{(T-h)S(h)_{d_1, d_2}} \le \frac{1}{(T-h)S(h)_{d_1, d_2}} \le \frac{1}{(T-h)pq_u} \le \frac{1}{cTpq_u}. $$

□

The Frobenius norm requires a little more work because of the sum it contains.

**Lemma 8 (Concentration of the sampling Bernoulli).** For all $u \in [0, 1]$,

$$ \mathbb{P}\left(\left|\frac{1}{T-h} \sum_{t=1}^{T-h} \pi_{t+h, d_1} \pi_{t, d_2} - S(h)_{d_1, d_2}\right| \ge uS(h)_{d_1, d_2}\right) \le c_1 \exp(-c_2 u^2 TS(h)_{d_1, d_2}). $$
$$ 
\begin{aligned}
&= \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{1}{S(h)_{d_1, d_2}} e'_{d_1} (\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' e_{d_2} \\
&= \frac{1}{T-h} \sum_{t=1}^{T-h} \text{Tr} \left[ \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} (\Pi_{t+h}^{\dagger} Y_{t+h}) (\Pi_t^{\dagger} Y_t)' \right]
\end{aligned}
$$

Equation (15) allows us to rewrite $(\Pi_{t+h}^{\dagger} Y_{t+h})(\Pi_t^{\dagger} Y_t)'$:

$$ 
\begin{aligned}
(\hat{\Gamma}_h + \mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} = & \frac{1}{T-h} \sum_{t=1}^{T-h} X_t' \text{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \text{diag}(\pi_{t+h}) X_{t+h} \\
& + \frac{1}{T-h} \sum_{t=1}^{T-h} \eta_t' \text{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \text{diag}(\pi_{t+h}) X_{t+h} \\
& + \frac{1}{T-h} \sum_{t=1}^{T-h} X_t' \text{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \text{diag}(\pi_{t+h}) \eta_{t+h} \\
& + \frac{1}{T-h} \sum_{t=1}^{T-h} \eta_t' \text{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \text{diag}(\pi_{t+h}) \eta_{t+h}
\end{aligned}
$$

Let us denote by $P_t$ the projection of $\mathbb{R}^{TD}$ keeping only the components associated with time $t$, i.e. such that $X_t = P_t X$ and $\eta_t = P_t \eta$. We recognize the following matrix $L$ in all four lines of the expression above:

$$ 
\begin{aligned}
L &= \frac{1}{T-h} \sum_{t=1}^{T-h} P_t' \text{diag}(\pi_t) \frac{e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} \text{diag}(\pi_{t+h}) P_{t+h} \\
&= \frac{1}{T-h} \sum_{t=1}^{T-h} P_t' \frac{\pi_{t+h, d_1} \pi_{t, d_2} e_{d_2} e'_{d_1}}{S(h)_{d_1, d_2}} P_{t+h}
\end{aligned}
\quad (16) $$

This leads to:

$$ (\hat{\Gamma}_h + \mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} = X' LX + \eta' LX + X' L\eta + \eta' L\eta $$

Since $X$ and $\eta$ both follow centered multivariate Gaussian distributions, we can express them as linear combinations of standard Gaussian vectors $g_\epsilon$ and $g_\eta$ of dimension $TD$ (indexed by the source of randomness): $X = \Psi_\epsilon g_\epsilon$ and $\eta = \Psi_\eta g_\eta$. Here $\Psi_\epsilon$ and $\Psi_\eta$ are the square roots of the respective covariance matrices

$$ \Psi_{\epsilon} = \text{Cov}[X]^{1/2} \quad \text{and} \quad \Psi_{\eta} = \text{Cov}[\eta]^{1/2} = \omega I. \quad (17) $$

We substitute $X$ and $\eta$ to get:

$$ (\hat{\Gamma}_h + \mathbf{1}_{\{h=0\}} \omega^2 I)_{d_2, d_1} = g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon + g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon + g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta + g'_\eta \Psi'_\eta L \Psi_\eta g_\eta, $$

which implies

$$ 
\begin{aligned}
(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2} = & g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon - \mathbb{E}[g'_\epsilon \Psi'_\epsilon L \Psi_\epsilon g_\epsilon] \\
& + g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon - \mathbb{E}[g'_\eta \Psi'_\eta L \Psi_\epsilon g_\epsilon] \\
& + g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta - \mathbb{E}[g'_\epsilon \Psi'_\epsilon L \Psi_\eta g_\eta]
\end{aligned}
$$
Proof. We distinguish three cases:

• When $d_1 = d_2$ and $h = 0$, we have $\pi_{t+h,d_1} = \pi_{t,d_2}$, which is a 2-state Markov chain with transition matrix $\mathcal{T} \otimes \mathcal{I}$, depicted on Figure 2a.

• When $d_1 \neq d_2$, the couple $(\pi_{t,d_2}, \pi_{t+h,d_1})$ is a 4-state Markov chain with transition matrix $\mathcal{T} \otimes \mathcal{T}$ since the chains $\pi_{t+h,d_1}$ and $\pi_{t,d_2}$ evolve along independent dimensions. It is shown on Figure 2b.

• When $d_1 = d_2$ and $h \ge 1$, we must study the $(h+1)$-tuple $(\pi_{t,d_1}, \pi_{t+1,d_1}, \dots, \pi_{t+h,d_1})$. It is a $2^{h+1}$-state Markov chain with transition matrix $S(h)$, whose non-reversible transition diagram can be seen on Figure 2c.

<figure><img src="image_8.png" /><figcaption>Figure 2: State space and transitions for the Markov chains used in the discrete concentration result</figcaption></figure>

In all of these cases, our variable of interest $\pi_{t+h,d_1} \pi_{t,d_2}$ is a function of the underlying Markov chain. The relevant functions are:

$$f_1 : x \mapsto x \qquad f_2 : (x, y) \mapsto yx \qquad f_3 : (x_0, \dots, x_h) \mapsto x_h x_0.$$

We note that since $\chi \le a, b \le 1 - \chi$, all the coefficients of $\mathcal{T}$ are greater than $\chi$. Furthermore, all the coefficients of $\mathcal{T} \otimes \mathcal{T}$ are greater than $\chi^2$. Finally, all the coefficients of $S(h)^{h+1}$ are greater than $\chi^{h+1}$, because all pairs of states are connected after $h+1$ steps. Let us illustrate this with $h=1$:

$$S(1) = \begin{pmatrix} 1-a & a & 0 & 0 \\ 0 & 0 & b & 1-b \\ 1-a & a & 0 & 0 \\ 0 & 0 & b & 1-b \end{pmatrix} \qquad S(1)^2 = \begin{pmatrix} (1-a)^2 & a(1-a) & ab & a(1-b) \\ (1-a)b & ab & (1-b)b & (1-b)^2 \\ (1-a)^2 & a(1-a) & ab & a(1-b) \\ (1-a)b & ab & (1-b)b & (1-b)^2 \end{pmatrix}.$$

Subsequently, all the transition matrices $\mathcal{R}$ we are interested in, namely $\mathcal{R} \in \{\mathcal{T}, \mathcal{T} \otimes \mathcal{T}, S(h)^{h+1}\}$, satisfy the Doeblin condition with $r = h+1$ and $\delta = \chi^{h+1}$:

$$\mathcal{R}^{h+1} \ge \chi^{h+1} \begin{pmatrix} 1 & \cdots & 1 \\ \vdots & \ddots & \vdots \\ 1 & \cdots & 1 \end{pmatrix}.$$

Since we only consider $h \in \{0, 1\}$ and since $\chi$ is fixed for our purposes, these chains fulfill the assumptions of Lemma 34. We thus conclude:

$$\mathbb{P}\left(\left| \frac{1}{T-h} \sum_{t=1}^{T-h} \pi_{t+h, d_1} \pi_{t, d_2} - S(h)_{d_1, d_2} \right| \ge uS(h)_{d_1, d_2} \right) \le c_1 \exp(-c_2 u^2 (T-h) S(h)_{d_1, d_2}).$$

We finally replace $T-h$ with $cT$ in the exponential, leading to the result we announced. $\square$
$$+ g'_{\eta} \Psi'_{\eta} L \Psi_{\eta} g_{\eta} - \mathbb{E}[g'_{\eta} \Psi'_{\eta} L \Psi_{\eta} g_{\eta}].$$

The union bound gives us the expected result. □

Now, our goal is to apply a Gaussian concentration inequality to these deviation probabilities. However, since $L$ is generated by the discrete sampling process $\pi$, it is random, and so are the products $\Psi'_a L \Psi_b$ (where $a, b \in \{\varepsilon, \eta\}$). We thus need a conditional version of the Hanson-Wright inequality (Lemma C.12), in which the following random variables will come into play:

* The spectral norm $\|\Psi'_a L \Psi_b\|_2$

* The Frobenius norm $\|\Psi'_a L \Psi_b\|_F^2$

* The shifted trace $\mathrm{Tr}(\Psi'_a L \Psi_b - \mathbb{E}[\Psi'_a L \Psi_b])$

### A.5. Interlude: Discrete concentration

We exploit discrete concentration results to bound the deviations of the three quantities we just mentioned, starting with the norms.

**Lemma A.6 (Norm reformulation for $L$).** The spectral and Frobenius norms of $L$ are given by

$$\|L\|_2 = \frac{\max_{t \in [T-h]} \pi_{t+h, d_1} \pi_{t, d_2}}{(T-h)S(h)_{d_1, d_2}}$$

$$\|L\|_F^2 = \frac{1}{(T-h)^2 S(h)_{d_1, d_2}} \sum_{t=1}^{T-h} \pi_{t+h, d_1} \pi_{t, d_2}.$$

*Proof.* We first notice that $L$ has a block-superdiagonal structure of rank $h$:

$$L = \frac{1}{T-h} \sum_{t=1}^{T-h} P'_{t} L_{[t,t+h]} P_{t+h} \quad \text{with blocks} \quad L_{[t,t+h]} = \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} e_{d_2} e'_{d_1}. \quad (18)$$

The spectral and Frobenius norms of such a matrix can easily be deduced from those of its blocks. Since $\|e_{d_2} e'_{d_1}\|_2 = \|e_{d_2} e'_{d_1}\|_F = 1$ and the $\pi_t$ are binary-valued, this leads to the following formulas:

$$\|L\|_2 = \frac{1}{T-h} \max_{t \in [T-h]} \|L_{[t,t+h]}\|_2 = \frac{1}{(T-h)S(h)_{d_1,d_2}} \max_{t \in [T-h]} \pi_{t+h, d_1} \pi_{t, d_2}$$

$$\|L\|_F^2 = \frac{1}{(T-h)^2} \sum_{t=1}^{T-h} \|L_{[t,t+h]}\|_F^2 = \frac{1}{(T-h)^2 S(h)_{d_1,d_2}^2} \sum_{t=1}^{T-h} \pi_{t+h, d_1} \pi_{t, d_2}. \quad \square$$

We can bound the spectral norm for free.

**Lemma A.7 (Spectral norm bound for $L$).** With probability 1, the spectral norm $\|L\|_2$ satisfies $\|L\|_2 \le \frac{c}{Tp^2}$.

*Proof.* Note that $S(h)_{d_1, d_2} \ge p^2$, and since $h \in \{0, 1\}$, we can state that $T-h \ge cT$. By Lemma A.6, we deduce

$$\|L\|_2 = \frac{\max_{t \in [T-h]} \pi_{t+h, d_1} \pi_{t, d_2}}{(T-h)S(h)_{d_1, d_2}} \le \frac{1}{(T-h)S(h)_{d_1, d_2}} \le \frac{1}{cTp^2}. \quad \square$$
Based on this concentration property, we can now bound the norms of the random matrix L with high probability.

**Lemma 9** (Frobenious norm bound for L). For any $\delta$ such that Equation (20) holds, with probability at least $1 - \delta$, the Frobenius norm $\|L\|_F^2$ satisfies

$$ \|L\|_F^2 \le \frac{c}{Tpq_u}. $$

*Proof.* By Lemma 8: for all $u \in [0, 1]$,

$$ \mathbb{P} \left( \frac{1}{T-h} \sum_{t=1}^{T-h} \pi_{t+h, d_1} \pi_{t, d_2} \ge (1+u)S(h)_{d_1, d_2} \right) \le c_1 \exp(-c_2 u^2 T S(h)_{d_1, d_2}). $$

We remember the expression of Lemma 6 for $\|L\|_F^2$ and notice that:

$$ \begin{aligned} & \mathbb{P} \left( \|L\|_F^2 \ge \frac{1+u}{(T-h)S(h)_{d_1,d_2}} \right) \\ &= \mathbb{P} \left( \frac{1}{(T-h)S(h)_{d_1,d_2}} \left( \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h,d_1} \pi_{t,d_2}}{S(h)_{d_1,d_2}} \right) \ge \frac{1}{(T-h)S(h)_{d_1,d_2}} (1+u) \right) \\ &\le c_1 \exp(-c_2 u^2 T S(h)_{d_1,d_2}) \end{aligned} $$

We finally recall that $S(h)_{d_1,d_2} \ge pq_u$ and $T-h \ge cT$, so that

$$ \begin{aligned} \mathbb{P} \left( \|L\|_F^2 \ge \frac{1+u}{cTpq_u} \right) &\le \mathbb{P} \left( \|L\|_F^2 \ge \frac{1+u}{(T-h)pq_u} \right) \\ &\le \mathbb{P} \left( \|L\|_F^2 \ge \frac{1+u}{(T-h)S(h)_{d_1,d_2}} \right) \\ &\le c_1 \exp(-c_2 u^2 T S(h)_{d_1,d_2}) \\ &\le c_1 \exp(-c_2 u^2 T pq_u). \end{aligned} $$

All we need to make sure that $\mathbb{P} \left( \|L\|_F^2 \ge \frac{1+u}{cTpq_u} \right) \le \delta$ is to choose $u$ such that

$$ c_1 \exp(-c_2 u^2 T pq_u) \le \delta \iff u \ge \sqrt{\frac{\log(c_1/\delta)}{c_2 T pq_u}} $$

Note that we can replace $\log(c_1/\delta)$ by a constant times $\log(1/\delta)$ to simplify expressions: this is possible as long as $\delta$ is chosen “small enough” (i.e. smaller than some universal constant). We will assume this fairly often in the rest of the proof.

For Lemma 8 to apply, we must ensure that our choice of $u$ is smaller than 1. With the previous discussion in mind, $u \le 1$ is implied by

$$ \sqrt{\frac{\log(1/\delta)}{Tpq_u}} \le c. \qquad (20) $$

If this holds, then we have

$$ \mathbb{P} \left( \|L\|_F^2 \ge \frac{2}{cTpq_u} \right) \le \mathbb{P} \left( \|L\|_F^2 \ge \frac{1+u}{cTpq_u} \right) \le \delta. $$

This yields the result we wanted. □

We now move on to studying the shifted trace of $\Psi'_a L \Psi_b$, which is the last ingredient we need for our application of Lemma 37.
The Frobenius norm requires a little more work because of the sum it contains.

**Lemma A.8 (Concentration of the sampling Bernoulli).** For all $u \in [0, 1]$,

$$ \mathbb{P}\left(\left|\frac{1}{T-h}\sum_{t=1}^{T-h} \pi_{t+h,d_1}\pi_{t,d_2} - S(h)_{d_1,d_2}\right| \ge uS(h)_{d_1,d_2}\right) \le c_1 \exp(-c_2 u^2 T p^2). $$

*Proof.* Recall that $S(h)_{d_1,d_2} = \mathbb{E}[\pi_{t+h,d_1}\pi_{t,d_2}]$. We distinguish three cases:

* When $d_1 = d_2$ and $h = 0$, we have $\pi_{t+h,d_1}$ equal to $\pi_{t,d_2}$, so their product is a Bernoulli $\mathcal{B}(p)$, and the sum is made of independent terms. The Chernoff bound of Lemma C.9 directly applies.

* When $d_1 \neq d_2 = d$ and $h = 1$, we have $\pi_{t+h,d_1}$ independent from $\pi_{t,d_2}$, so their product is a Bernoulli $\mathcal{B}(pp)$, and the sum is made of independent terms. Lemma C.9 directly applies here as well.

* When $d_1 = d_2$ and $h = 1$, we have $\pi_{t+h,d_1}$ independent from $\pi_{t,d_2}$, so their product is a Bernoulli $\mathcal{B}(p^2)$, but the sum is no longer made of independent terms. However, the terms become independent if we split the sum between even and odd indices. Notice that the event

$$ \left| \frac{1}{T-h} \sum_{t=1}^{T-h} \pi_{t+h,d_1} \pi_{t,d_2} - S(h)_{d_1,d_2} \right| \ge uS(h)_{d_1,d_2} $$

can only occur if at least one of the two following events occurs:

$$ \left| \frac{1}{T-h} \sum_{\substack{t=1 \\ \text{even}}}^{T-h} \pi_{t+h,d_1} \pi_{t,d_2} - \frac{S(h)_{d_1,d_2}}{2} \right| \ge u \frac{S(h)_{d_1,d_2}}{2} $$

$$ \left| \frac{1}{T-h} \sum_{\substack{t=1 \\ \text{odd}}}^{T-h} \pi_{t+h,d_1} \pi_{t,d_2} - \frac{S(h)_{d_1,d_2}}{2} \right| \ge u \frac{S(h)_{d_1,d_2}}{2} $$

By the union bound,

$$ \mathbb{P}\left(\left|\frac{1}{T-h}\sum_{t=1}^{T-h} \pi_{t+h,d_1}\pi_{t,d_2} - S(h)_{d_1,d_2}\right| \ge uS(h)_{d_1,d_2}\right) $$

$$ \le \mathbb{P}\left(\left|\frac{1}{T-h}\sum_{\substack{t=1 \\ \text{even}}}^{T-h} \pi_{t+h,d_1}\pi_{t,d_2} - \frac{S(h)_{d_1,d_2}}{2}\right| \ge u\frac{S(h)_{d_1,d_2}}{2}\right) $$

$$ + \mathbb{P}\left(\left|\frac{1}{T-h}\sum_{\substack{t=1 \\ \text{odd}}}^{T-h} \pi_{t+h,d_1}\pi_{t,d_2} - \frac{S(h)_{d_1,d_2}}{2}\right| \ge u\frac{S(h)_{d_1,d_2}}{2}\right) $$

And we can apply Lemma C.9 to each of the two terms separately.

In all three cases, we thus get

$$ \mathbb{P}\left(\left|\frac{1}{T-h}\sum_{t=1}^{T-h} \pi_{t+h,d_1}\pi_{t,d_2} - S(h)_{d_1,d_2}\right| \ge uS(h)_{d_1,d_2}\right) \le c_1 \exp(-c_2 u^2 T S(h)_{d_1,d_2}). $$
**Lemma 10** (Trace bound for the L matrices). For all $u \in [0, 1]$,

$$ \mathbb{P}(|\operatorname{Tr}(\Psi'_{\varepsilon} L \Psi_{\varepsilon} - \mathbb{E}[\Psi'_{\varepsilon} L \Psi_{\varepsilon}])| \ge u) \le c_1 \exp \left(-\frac{c_2 u^2 T p q_u}{\|\Gamma_h\|_2^2}\right) $$

$$ \mathbb{P}(|\operatorname{Tr}(\Psi'_{\eta} L \Psi_{\eta} - \mathbb{E}[\Psi'_{\eta} L \Psi_{\eta}])| \ge u) \le c_1 \exp \left(-\frac{c_2 u^2 T p q_u}{\omega^4}\right). $$

*Proof.* We can compute an explicit formula thanks to Equation (19): if $a \in \{\varepsilon, \eta\}$ then

$$ \begin{aligned} \operatorname{Tr}(\Psi'_{a} L \Psi_{a}) &= \operatorname{Tr} \left( \frac{1}{T-h} \sum_{t=1}^{T-h} \Psi'_{a} P'_{t} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} e_{d_2} e'_{d_1} P_{t+h} \Psi_{a} \right) \\ &= \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} \operatorname{Tr} (\Psi'_{a} P'_{t} e_{d_2} e'_{d_1} P_{t+h} \Psi_{a}) \\ &= \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} (e'_{d_1} P_{t+h} \Psi_{a} \Psi'_{a} P'_{t} e_{d_2}) \\ &= \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} ((\Psi_{a} \Psi'_{a})_{[t+h, t]})_{d_1, d_2} \end{aligned} $$

where $((\Psi_a \Psi'_a)_{[t+h, t]})_{d_1, d_2}$ denotes the $(d_1, d_2)$ coefficient of the $(t+h, t)$ block of $\Psi_a \Psi'_a$. Now is the time to look back on Equation (18), which tells us that both $\Psi_\varepsilon \Psi'_\varepsilon$ and $\Psi_\eta \Psi'_\eta$ are constant along their superdiagonal of rank $h$. We thus find that

$$ \operatorname{Tr}(\Psi'_{\varepsilon} L \Psi_{\varepsilon} - \mathbb{E}[\Psi'_{\varepsilon} L \Psi_{\varepsilon}]) = (\Gamma_h)_{d_1, d_2} \left( \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} - 1 \right) $$

$$ \operatorname{Tr}(\Psi'_{\eta} L \Psi_{\eta} - \mathbb{E}[\Psi'_{\eta} L \Psi_{\eta}]) = (\mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} \left( \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} - 1 \right) $$

Like before, we can apply Lemma 8: for all $u \in [0, 1]$,

$$ \mathbb{P}\left(\left|\frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} - 1\right| \ge u\right) \le c_1 \exp(-c_2 u^2 T S(h)_{d_1, d_2}) \le c_1 \exp(-c_2 u^2 T p q_u). $$

Since $|(\Gamma_h)_{d_1, d_2}| \le \|\Gamma_h\|_2$ and $(\mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} \le \omega^2$, we can deduce

$$ \mathbb{P}(|\operatorname{Tr}(\Psi'_{\varepsilon} L \Psi_{\varepsilon} - \mathbb{E}[\Psi'_{\varepsilon} L \Psi_{\varepsilon}])| \ge u \|\Gamma_h\|_2) \le c_1 \exp(-c_2 u^2 T p q_u) $$

$$ \mathbb{P}(|\operatorname{Tr}(\Psi'_{\eta} L \Psi_{\eta} - \mathbb{E}[\Psi'_{\eta} L \Psi_{\eta}])| \ge u \omega^2) \le c_1 \exp(-c_2 u^2 T p q_u) $$

which, after rescaling, yields the result we announced. □

## A.6 GAUSSIAN CONCENTRATION, EPISODE 2

We are now ready to apply our conditional concentration result.

**Lemma 11** (Applying Hanson-Wright). Let $\delta > 0$ and $u \in [0, 1]$. Assume that Equations (20) and (21) hold. Then the deviation probability for $(\hat{\Gamma}_h)_{d_1, d_2}$ satisfies

$$ \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \ge u) \le 4\delta + c_1 \exp \left( - \frac{c_2 u^2 T p q_u}{\max \{ (\|\Psi_\varepsilon\|_2^2 + \omega^2)^2, \|\Gamma_h\|_2^2, \omega^4 \}} \right). $$
We finally replace $T-h$ with $cT$ and $S(h)_{d_1,d_2}$ with $p^2$ in the exponential, leading to the result we announced. □

Based on this concentration property, we can now bound the norms of the random matrix $L$ with high probability.

**Lemma A.9** (Frobenious norm bound for $L$). For any $\delta$ such that Equation (19) holds, with probability at least $1 - \delta$, the Frobenius norm $\|L\|_F^2$ satisfies $\|L\|_F^2 \le \frac{c}{Tp^2}$.

*Proof.* By Lemma A.8: for all $u \in [0, 1]$,

$$ \mathbb{P}\left(\frac{1}{T-h} \sum_{t=1}^{T-h} \pi_{t+h, d_1} \pi_{t, d_2} \ge (1+u)S(h)_{d_1, d_2}\right) \le c_1 \exp(-c_2 u^2 T p^2). $$

We remember the expression of Lemma A.6 for $\|L\|_F^2$ and notice that:

$$ \begin{aligned} & \mathbb{P}\left(\|L\|_F^2 \ge \frac{1+u}{(T-h)S(h)_{d_1,d_2}}\right) \\ &= \mathbb{P}\left(\frac{1}{(T-h)S(h)_{d_1,d_2}} \left(\frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h,d_1}\pi_{t, d_2}}{S(h)_{d_1,d_2}}\right) \ge \frac{1}{(T-h)S(h)_{d_1,d_2}}(1+u)\right) \\ &\le c_1 \exp(-c_2 u^2 T q^2 p^2) \end{aligned} $$

We finally recall that $S(h)_{d_1,d_2} \ge p^2$ and $T-h \ge cT$, so that

$$ \begin{aligned} \mathbb{P}\left(\|L\|_F^2 \ge \frac{1+u}{cTp^2}\right) &\le \mathbb{P}\left(\|L\|_F^2 \ge \frac{1+u}{(T-h)S(h)_{d_1,d_2}}\right) \\ &\le c_1 \exp(-c_2 u^2 T p^2) \end{aligned} $$

All we need to make sure that $\mathbb{P}\left(\|L\|_F^2 \ge \frac{1+u}{cTp^2}\right) \le \delta$ is to choose $u$ such that

$$ c_1 \exp(-c_2 u^2 T p^2) \le \delta \iff u \ge \sqrt{\frac{\log(c_1/\delta)}{c_2 T p^2}} $$

Note that we can replace $\log(c_1/\delta)$ by a constant times $\log(1/\delta)$ to simplify expressions: this is possible as long as $\delta$ is chosen “small enough” (i.e. smaller than some universal constant). We will assume this fairly often in the rest of the proof.

For Lemma A.8 to apply, we must ensure that our choice of $u$ is smaller than 1. With the previous discussion in mind, $u \le 1$ is implied by

$$ \sqrt{\frac{\log(1/\delta)}{Tp^2}} \le c. \qquad (19) $$

If this holds, then we have

$$ \mathbb{P}\left(\|L\|_F^2 \ge \frac{2}{cTp^2}\right) \le \mathbb{P}\left(\|L\|_F^2 \ge \frac{1+u}{cTp^2}\right) \le \delta. $$

This yields the result we wanted. □
Proof. The conclusion we had reached before our discrete interlude is given by Lemma 5, and we can rewrite it as

$$ \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \ge u) \le p_{\varepsilon\varepsilon} + p_{\eta\varepsilon} + p_{\varepsilon\eta} + p_{\eta\eta}, $$

where each $p_{ab}$ represents a deviation probability for a specific quadratic form $g'_a \Psi'_a L \Psi_b g_b$. Let us choose $\delta$ such that Equation (20) holds. By Lemmas 7 and 9, with probability at least $1 - \delta$, the following eight inequalities occur at the same time (we use Lemma 25 to split the products):

$$ \begin{aligned} \|\Psi'_\epsilon L \Psi_\epsilon\|_F^2 &\le \frac{c \|\Psi_\epsilon\|_2^4}{Tpq_u} & \|\Psi'_\epsilon L \Psi_\epsilon\|_2^2 &\le \frac{c \|\Psi_\epsilon\|_2^2}{Tpq_u} \\ \|\Psi'_\eta L \Psi_\epsilon\|_F^2 &\le \frac{c \|\Psi_\eta\|_2^2 \|\Psi_\epsilon\|_2^2}{Tpq_u} & \|\Psi'_\eta L \Psi_\epsilon\|_2^2 &\le \frac{c \|\Psi_\eta\|_2 \|\Psi_\epsilon\|_2}{Tpq_u} \\ \|\Psi'_\epsilon L \Psi_\eta\|_F^2 &\le \frac{c \|\Psi_\epsilon\|_2^2 \|\Psi_\eta\|_2^2}{Tpq_u} & \|\Psi'_\epsilon L \Psi_\eta\|_2^2 &\le \frac{c \|\Psi_\epsilon\|_2 \|\Psi_\eta\|_2}{Tpq_u} \\ \|\Psi'_\eta L \Psi_\eta\|_F^2 &\le \frac{c \|\Psi_\eta\|_2^4}{Tpq_u} & \|\Psi'_\eta L \Psi_\eta\|_2^2 &\le \frac{c \|\Psi_\eta\|_2^2}{Tpq_u}. \end{aligned} $$

The spectral norm of $\Psi_\eta$ is easily seen to equal $\|\Psi_\eta\|_2 = \|\omega^2 I\|_2^{1/2} = \omega$, which allows us to lighten these expressions. From there, Lemma 37 (applied with $X = g_a$, $Y = g_b$ and $A = \Psi'_a L \Psi_b$) provides the concentration bounds we need<fnref n="2" />:

$$ \begin{aligned} p_{\varepsilon\varepsilon} &\le \delta + 2 \exp \left( -cTpq_u \min \left\{ \frac{(u/4)^2}{\|\Psi_\epsilon\|_2^4}, \frac{(u/4)}{\|\Psi_\epsilon\|_2^2} \right\} \right) + \mathbb{P} \left( |\mathrm{Tr}(\Psi'_\eta L \Psi_\epsilon) - \mathbb{E}[\Psi'_\epsilon L \Psi_\epsilon]| \ge u/8 \right) \\ p_{\eta\varepsilon} &\le \delta + 2 \exp \left( -cTpq_u \min \left\{ \frac{(u/4)^2}{\omega^2 \|\Psi_\epsilon\|_2^2}, \frac{(u/4)}{\omega \|\Psi_\epsilon\|_2} \right\} \right) \\ p_{\varepsilon\eta} &\le \delta + 2 \exp \left( -cTpq_u \min \left\{ \frac{(u/4)^2}{\|\Psi_\epsilon\|_2^2 \omega^2}, \frac{(u/4)}{\|\Psi_\epsilon\|_2 \omega} \right\} \right) \\ p_{\eta\eta} &\le \delta + 2 \exp \left( -cTpq_u \min \left\{ \frac{(u/4)^2}{\omega^4}, \frac{(u/4)}{\omega^2} \right\} \right) + \mathbb{P} \left( |\mathrm{Tr}(\Psi'_\eta L \Psi_\epsilon) - \mathbb{E}[\Psi'_\eta L \Psi_\epsilon]| \ge u/8 \right). \end{aligned} $$

The denominators inside the minima can be unified: for the left column,

$$ \max \{ \|\Psi_\epsilon\|_2^4, \|\Psi_\epsilon\|_2^2 \omega^2, \omega^4 \} \le (\|\Psi_\epsilon\|_2^2 + \omega^2)^2, $$

and for the right column,

$$ \max \{ \|\Psi_\epsilon\|_2^2, \|\Psi_\epsilon\|_2 \omega, \omega^2 \} \le (\|\Psi_\epsilon\|_2 + \omega)^2 \le 2 (\|\Psi_\epsilon\|_2^2 + \omega^2). $$

This means we can upper bound each of the four minima by

$$ \min \left\{ \left( \frac{u/4}{\|\Psi_\epsilon\|_2^2 + \omega^2} \right)^2, \frac{u/8}{\|\Psi_\epsilon\|_2^2 + \omega^2} \right\}. $$

From now on, we additionally suppose that

$$ \frac{u/4}{\|\Psi_\epsilon\|_2^2 + \omega^2} \le \frac{1}{2} \qquad (21) $$

<footnote n="2">²The additional trace terms that appear when applying Lemma 37 (as opposed to the non-conditional version of Lemma 36) are absent from the papers by Rao et al. [2017a,b], which is why we think their upper bound proofs are incomplete.</footnote>
We now move on to studying the shifted trace of $\Psi'_a L \Psi_b$, which is the last ingredient we need for our application of Lemma C.12.

**Lemma A.10** (Trace bound for the $L$ matrices). For all $u \in [0, 1]$,

$$ \mathbb{P}(|\mathrm{Tr}(\Psi'_\epsilon L \Psi_\epsilon - \mathbb{E}[\Psi'_\epsilon L \Psi_\epsilon])| \ge u) \le c_1 \exp\left(-\frac{c_2 u^2 T p^2}{\|\Gamma_h\|_2^2}\right) $$

$$ \mathbb{P}(|\mathrm{Tr}(\Psi'_\eta L \Psi_\eta - \mathbb{E}[\Psi'_\eta L \Psi_\eta])| \ge u) \le c_1 \exp\left(-\frac{c_2 u^2 T p^2}{\omega^4}\right). $$

*Proof.* We can compute an explicit formula thanks to Equation (18): if $a \in \{\epsilon, \eta\}$ then

$$ \begin{aligned} \mathrm{Tr}(\Psi'_a L \Psi_a) &= \mathrm{Tr} \left( \frac{1}{T-h} \sum_{t=1}^{T-h} \Psi'_a P'_t \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} e_{d_2} e'_{d_1} P_{t+h} \Psi_a \right) \\ &= \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} \mathrm{Tr} \left( \Psi'_a P'_t e_{d_2} e'_{d_1} P_{t+h} \Psi_a \right) \\ &= \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} \left( e'_{d_1} P_{t+h} \Psi_a \Psi'_a P'_t e_{d_2} \right) \\ &= \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} \left( (\Psi_a \Psi'_a)_{[t+h, t]} \right)_{d_1, d_2} \end{aligned} $$

where $((\Psi_a \Psi'_a)_{[t+h,t]})_{d_1,d_2}$ denotes the $(d_1, d_2)$ coefficient of the $(t+h, t)$ block of $\Psi_a \Psi'_a$. Now is the time to look back on Equation (17), which tells us that both $\Psi_\epsilon \Psi'_\epsilon$ and $\Psi_\eta \Psi'_\eta$ are constant along their superdiagonal of rank $h$. We thus find that

$$ \mathrm{Tr}(\Psi'_\epsilon L \Psi_\epsilon - \mathbb{E}[\Psi'_\epsilon L \Psi_\epsilon]) = (\Gamma_h)_{d_1, d_2} \left( \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} - 1 \right) $$

$$ \mathrm{Tr}(\Psi'_\eta L \Psi_\eta - \mathbb{E}[\Psi'_\eta L \Psi_\eta]) = (\mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1, d_2} \left( \frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} - 1 \right) $$

Like before, we can apply Lemma A.8: for all $u \in [0, 1]$,

$$ \mathbb{P}\left(\left|\frac{1}{T-h} \sum_{t=1}^{T-h} \frac{\pi_{t+h, d_1} \pi_{t, d_2}}{S(h)_{d_1, d_2}} - 1\right| \ge u\right) \le c_1 \exp(-c_2 u^2 T p^2) $$

Since $|(\Gamma_h)_{d_1,d_2}| \le \|\Gamma_h\|_2$ and $(\mathbf{1}_{\{h=0\}} \omega^2 I)_{d_1,d_2} \le \omega^2$, we can deduce

$$ \mathbb{P}(|\mathrm{Tr}(\Psi'_\epsilon L \Psi_\epsilon - \mathbb{E}[\Psi'_\epsilon L \Psi_\epsilon])| \ge u \|\Gamma_h\|_2) \le c_1 \exp(-c_2 u^2 T p^2) $$

$$ \mathbb{P}(|\mathrm{Tr}(\Psi'_\eta L \Psi_\eta - \mathbb{E}[\Psi'_\eta L \Psi_\eta])| \ge u \omega^2) \le c_1 \exp(-c_2 u^2 T p^2) $$

which, after rescaling, yields the result we announced. □

## A.6. Gaussian concentration, episode 2

We are now ready to apply our conditional concentration result.
This enables us to get rid of these minima by reducing them to the (smaller) quadratic term on the left. We end up with

$$
\begin{align*}
p_{\varepsilon\varepsilon} &\le \delta + 2 \exp\left(-\frac{cu^2Tpq_u}{(\|\Psi_\varepsilon\|_2^2 + \omega^2)^2}\right) + \mathbb{P}\left(|\operatorname{Tr}(\Psi'_\varepsilon L\Psi_\varepsilon) - \mathbb{E}[\Psi'_\varepsilon L\Psi_\varepsilon]| \ge u/8\right) \\
p_{\eta\varepsilon} &\le \delta + 2 \exp\left(-\frac{cu^2Tpq_u}{(\|\Psi_\varepsilon\|_2^2 + \omega^2)^2}\right) \\
p_{\varepsilon\eta} &\le \delta + 2 \exp\left(-\frac{cu^2Tpq_u}{(\|\Psi_\varepsilon\|_2^2 + \omega^2)^2}\right) \\
p_{\eta\eta} &\le \delta + 2 \exp\left(-\frac{cu^2Tpq_u}{(\|\Psi_\eta\|_2^2 + \omega^2)^2}\right) + \mathbb{P}\left(|\operatorname{Tr}(\Psi'_\eta L\Psi_\eta) - \mathbb{E}[\Psi'_\eta L\Psi_\eta]| \ge u/8\right).
\end{align*}
$$

As for the trace terms, they are taken care of by Lemma 10:

$$
\begin{align*}
\mathbb{P} \left( |\operatorname{Tr}(\Psi'_{\eta} L \Psi_{\varepsilon}) - \mathbb{E}[\Psi'_{\varepsilon} L \Psi_{\varepsilon}]| \ge u/8 \right) &\le c_3 \exp \left( -c_4 \frac{(u/8)^2 T pq_u}{\|\Gamma_h\|_2^2} \right) \\
\mathbb{P} \left( |\operatorname{Tr}(\Psi'_{\eta} L \Psi_{\eta}) - \mathbb{E}[\Psi'_{\eta} L \Psi_{\eta}]| \ge u/8 \right) &\le c_3 \exp \left( -c_4 \frac{(u/8)^2 T pq_u}{\omega^4} \right)
\end{align*}
$$

We plug this in and rearrange to get:

$$
p_{\varepsilon\varepsilon} + p_{\eta\varepsilon} + p_{\varepsilon\eta} + p_{\eta\eta} \le 4\delta + c_1 \exp \left( - \frac{c_2 u^2 T pq_u}{\max \{ (\|\Psi_\varepsilon\|_2^2 + \omega^2)^2, \|\Gamma_h\|_2^2, \omega^4 \}} \right).
$$

<figure><img src="image_9.png" /></figure>

The following result will simplify the denominator inside the exponential.

**Lemma 12 (Spectral norms of Ψ<sub>ε</sub> and Γ<sub>h</sub>).* The matrices Ψ<sub>ε</sub> and Γ<sub>h</sub> satisfy:

$$
\|\Psi_{\varepsilon}\|_{2}^{2} \le \frac{\sigma_{\max}^{2}}{(1 - \vartheta)^{2}} \quad \text{and} \quad \|\Gamma_{h}\|_{2} \le \frac{\vartheta^{h} \sigma_{\max}^{2}}{1 - \vartheta}.
$$

As a consequence,

$$
\max \{ (\|\Psi_\varepsilon\|_2^2 + \omega^2)^2, \|\Gamma_h\|_2^2, \omega^4 \} \le \frac{(\sigma_{\max}^2 + \omega^2)^2}{(1 - \vartheta)^4}
$$

*Proof.* By Lemma 1, we can write $\Psi_{\varepsilon}^2$ as a sum of Kronecker products (one for each block). Let $J_t$ be a matrix full of zeros, except for the subdiagonal of rank $t$, which is full of ones. Then we have:

$$
\Psi_{\varepsilon}^{2} = \mathrm{Cov}[X] = I \otimes \Gamma_{0}(\theta) + \sum_{t=1}^{T-1} [J_{t} \otimes \theta^{t} \Gamma_{0}(\theta) + J'_{t} \otimes \Gamma_{0}(\theta) \theta'^{t}]
$$

This gives us control over its spectral norm thanks to Lemma 24:

$$
\begin{align*}
\|\Psi_{\varepsilon}\|_{2}^{2} &= \|\Psi_{\varepsilon}^{2}\|_{2} \le \|I\|_{2} \times \|\Gamma_{0}(\theta)\|_{2} + \sum_{t=1}^{T-1} [\|J_{t}\|_{2} \times \|\theta^{t}\Gamma_{0}(\theta)\|_{2} + \|J'_{t}\|_{2} \times \|\Gamma_{0}(\theta)\theta^{t}\|_{2}] \\
&\le \|\Gamma_{0}(\theta)\|_{2} \left(1 + 2 \sum_{t=1}^{T-1} \|\theta\|_{2}^{t}\right)
\end{align*}
$$
**Remark.** This is where our proof corrects the one of Rao et al. [31]. Their supplementary material is available at https://web.stanford.edu/~milind/papers/system_id_icassp_proof.pdf. On page 2, in the middle of the right column, they use a Gaussian concentration inequality that is similar to unconditional Hanson-Wright (Lemma C.11). It is applied to a matrix $L_1 = B'_T A_T B_T$ (in their notations), where $A_T$ is a sum of the random projection coefficients $P(k)_{i,j}$. Thus, the matrix $L_1$ is itself random, and they should have applied conditional Hanson-Wright (Lemma C.12) instead. When we compare Lemma C.12 with Lemma C.11, we note the presence of an additional trace term, which needs to be analyzed separately. This term is missing from the proof in Rao et al. [31], which is why we deem it incomplete.

**Lemma A.11** (Applying Hanson-Wright). Let $\delta > 0$ and $u \in [0, 1]$. Assume that Equations (19) and (20) hold. Then the deviation probability for $(\hat{\Gamma}_h)_{d_1, d_2}$ satisfies

$$ \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \ge u) \le 4\delta + c_1 \exp \left( - \frac{c_2 u^2 T p^2}{\max \{ (\|\Psi_\epsilon\|_2^2 + \omega^2)^2, \|\Gamma_h\|_2^2, \omega^4 \}} \right). $$

*Proof.* The conclusion we had reached before our discrete interlude is given by Lemma A.5, and we can rewrite it as

$$ \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \ge u) \le \Delta_{\epsilon\epsilon} + \Delta_{\eta\epsilon} + \Delta_{\epsilon\eta} + \Delta_{\eta\eta}, $$

where each $\Delta_{ab}$ represents a deviation probability for the quadratic form $g'_a \Psi'_a L \Psi_b g_b$. Let us choose $\delta$ such that Equation (19) holds. By Lemmas A.7 and A.9, with probability at least $1 - \delta$, the following eight inequalities occur at the same time (we use Lemma C.3 to split the products):

$$ \|\Psi'_\epsilon L \Psi_\epsilon\|_F^2 \le \frac{c \|\Psi_\epsilon\|_2^4}{Tp^2} $$

$$ \|\Psi'_\eta L \Psi_\epsilon\|_F^2 \le \frac{c \|\Psi_\eta\|_2^2 \|\Psi_\epsilon\|_2^2}{Tp^2} $$

$$ \|\Psi'_\epsilon L \Psi_\eta\|_F^2 \le \frac{c \|\Psi_\epsilon\|_2^2 \|\Psi_\eta\|_2^2}{Tp^2} $$

$$ \|\Psi'_\eta L \Psi_\eta\|_F^2 \le \frac{c \|\Psi_\eta\|_2^4}{Tp^2} $$

$$ \|\Psi'_\epsilon L \Psi_\epsilon\|_2 \le \frac{c \|\Psi_\epsilon\|_2^2}{Tp^2} $$

$$ \|\Psi'_\eta L \Psi_\epsilon\|_2 \le \frac{c \|\Psi_\eta\|_2 \|\Psi_\epsilon\|_2}{Tp^2} $$

$$ \|\Psi'_\epsilon L \Psi_\eta\|_2 \le \frac{c \|\Psi_\epsilon\|_2 \|\Psi_\eta\|_2}{Tp^2} $$

$$ \|\Psi'_\eta L \Psi_\eta\|_2 \le \frac{c \|\Psi_\eta\|_2^2}{Tp^2}. $$

The spectral norm of $\Psi_\eta$ is easily seen to equal $\|\Psi_\eta\|_2 = \|\omega^2 I\|_2^{1/2} = \omega$, which allows us to lighten these expressions. From there, Lemma C.12 (applied with $X = g_a$, $Y = g_b$ and $A = \Psi'_a L \Psi_b$) provides the concentration bounds we need:

$$ \Delta_{\epsilon\epsilon} \le \delta + 2 \exp \left( -cTp^2 \min \left\{ \frac{(u/4)^2}{\|\Psi_\epsilon\|_2^4}, \frac{(u/4)}{\|\Psi_\epsilon\|_2^2} \right\} \right) + \mathbb{P} \left( |\mathrm{Tr}(\Psi'_\eta L \Psi_\epsilon) - \mathbb{E}[\Psi'_\eta L \Psi_\epsilon]| \ge u/8 \right) $$

$$ \Delta_{\eta\epsilon} \le \delta + 2 \exp \left( -cTp^2 \min \left\{ \frac{(u/4)^2}{\omega^2 \|\Psi_\epsilon\|_2^2}, \frac{(u/4)}{\omega \|\Psi_\epsilon\|_2} \right\} \right) $$
We now use Lemma 2:

$$ \|\Psi_{\varepsilon}\|_{2}^{2} \leq \frac{\sigma_{\max}^{2}}{1 - \vartheta^{2}} \left(1 + 2 \frac{\vartheta}{1 - \vartheta}\right) = \frac{\sigma_{\max}^{2}}{(1 - \vartheta)^{2}}. $$

We now turn to $\Gamma_h$ with Lemmas 1 and 2:

$$ \|\Gamma_h\|_2 = \|\theta^h \Gamma_0(\theta)\|_2 \leq \frac{\vartheta^h \sigma_{\max}^2}{1 - \vartheta^2}. $$

In particular, we have

$$ \begin{aligned} \max \{ (\|\Psi_{\varepsilon}\|_2^2 + \omega^2)^2, \|\Gamma_h\|_2^2, \omega^4 \} &\leq \max \left\{ \left( \frac{\sigma_{\max}^2}{(1-\vartheta)^2} + \omega^2 \right)^2, \left( \frac{\vartheta^h \sigma_{\max}^2}{1-\vartheta^2} \right)^2, \omega^4 \right\} \\ &\leq \frac{(\sigma_{\max}^2 + \omega^2)^2}{(1-\vartheta)^4} \end{aligned} $$

□

We can now control the error of the covariance estimator:

**Lemma 13** (Max norm convergence rate of the covariance estimator). Let $\delta > 0$ be small enough. Assume that Equations (20) and (23) hold. Then the covariance estimator $\hat{\Gamma}_h$ from Equation (9) satisfies

$$ \|\hat{\Gamma}_h - \Gamma_h\|_{\max} \leq c \frac{\sigma_{\max}^2 + \omega^2}{(1 - \vartheta)^2} \frac{\sqrt{\log(D/\delta)}}{\sqrt{Tpq_u}} = \text{err}(\delta) $$

with probability greater than $1 - \delta$.

*Proof.* Let us plug Lemma 12 into Lemma 11

$$ \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \geq u) \leq 4\delta + c_1 \exp \left( - \frac{c_2 (1 - \vartheta)^4 u^2 T pq_u}{(\sigma_{\max}^2 + \omega^2)^2} \right). $$

All that is left to do is choose $u$ such that

$$ \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \geq u) \leq 8\delta, $$

which will be true if

$$ c_1 \exp \left( - \frac{c_2 (1 - \vartheta)^4 T pq_u}{(\sigma_{\max}^2 + \omega^2)^2} u^2 \right) \leq 4\delta \iff u \geq \sqrt{\frac{\log(c_1/4\delta)(\sigma_{\max}^2 + \omega^2)^2}{c_2(1 - \vartheta)^4 T pq_u}}. $$

As long as $\delta$ is small enough, we can take

$$ u = c \frac{\sqrt{\log(1/\delta)} (\sigma_{\max}^2 + \omega^2)}{(1 - \vartheta)^2 \sqrt{Tpq_u}}. \qquad (22) $$

For Lemma 11 to apply, we must verify that $u \in [0, 1]$ and that Equation (21) is satisfied. In other words, we have to ensure that

$$ c \frac{\sqrt{\log(1/\delta)} (\sigma_{\max}^2 + \omega^2)}{(1 - \vartheta)^2 \sqrt{Tpq_u}} \leq \min\{1, 2(\|\Psi_{\varepsilon}\|_2^2 + \omega^2)\} $$

Using Lemma 12, this is implied by the condition

$$ \frac{\sqrt{\log(1/\delta)} \max\{1, (\sigma_{\max}^2 + \omega^2)^{-1}\}}{(1 - \vartheta)^2 \sqrt{Tpq_u}} \leq c \qquad (23) $$
$$
\begin{align*}
\Delta_{\varepsilon\eta} &\le \delta + 2 \exp \left( -cT p^2 \min \left\{ \frac{(u/4)^2}{\|\Psi_{\varepsilon}\|_2^2 \omega^2}, \frac{(u/4)}{\|\Psi_{\varepsilon}\|_2 \omega} \right\} \right) \\
\Delta_{\eta\eta} &\le \delta + 2 \exp \left( -cT p^2 \min \left\{ \frac{(u/4)^2}{\omega^4}, \frac{(u/4)}{\omega^2} \right\} \right) \\
&\quad + \mathbb{P} \left( |\operatorname{Tr}(\Psi'_{\eta} L \Psi_{\eta}) - \mathbb{E}[\Psi'_{\eta} L \Psi_{\eta}]| \ge u/8 \right).
\end{align*}
$$

The denominators inside the minima can be unified: for the left column,

$$
\max \{ \|\Psi_{\varepsilon}\|_{2}^{4}, \|\Psi_{\varepsilon}\|_{2}^{2} \omega^{2}, \omega^{4} \} \leq (\|\Psi_{\varepsilon}\|_{2}^{2} + \omega^{2})^{2},
$$

and for the right column,

$$
\max \{ \|\Psi_{\varepsilon}\|_{2}^{2}, \|\Psi_{\varepsilon}\|_{2} \omega, \omega^{2} \} \leq (\|\Psi_{\varepsilon}\|_{2} + \omega)^{2} \leq 2 (\|\Psi_{\varepsilon}\|_{2}^{2} + \omega^{2}).
$$

This means we can upper bound each of the four minima by

$$
\min \left\{ \left( \frac{u/4}{\|\Psi_{\varepsilon}\|_{2}^{2} + \omega^{2}} \right)^{2}, \frac{u/8}{\|\Psi_{\varepsilon}\|_{2}^{2} + \omega^{2}} \right\} .
$$

From now on, we additionally suppose that

$$
\frac{u/4}{\|\Psi_{\varepsilon}\|_{2}^{2} + \omega^{2}} \le \frac{1}{2} \qquad (20)
$$

This enables us to get rid of these minima by reducing them to the (smaller) quadratic term
on the left. We end up with

$$
\Delta_{\varepsilon\varepsilon} \le \delta + 2 \exp \left( - \frac{cu^2Tp^2}{\left( \|\Psi_\varepsilon\|_2^2 + \omega^2 \right)^2} \right) + \mathbb{P} \left( |\operatorname{Tr}(\Psi'_\varepsilon L \Psi_\varepsilon) - \mathbb{E}[\Psi'_\varepsilon L \Psi_\varepsilon]| \ge u/8 \right)
$$

$$
\Delta_{\eta\epsilon} \le \delta + 2 \exp \left( - \frac{cu^2Tp^2}{\left( \|\Psi_{\epsilon}\|_2^2 + \omega^2 \right)^2} \right)
$$

$$
\Delta_{\epsilon\eta} \le \delta + 2 \exp \left( - \frac{cu^2Tp^2}{\left( \|\Psi_{\epsilon}\|_{2}^{2} + \omega^{2} \right)^{2}} \right)
$$

$$
\Delta_{\eta\eta} \le \delta + 2 \exp \left( - \frac{cu^2Tp^2}{\left( \|\Psi_{\epsilon}\|_{2}^{2} + \omega^2 \right)^2} \right) + \mathbb{P} \left( |\operatorname{Tr}(\Psi'_{\eta} L \Psi_{\eta}) - \mathbb{E}[\Psi'_{\eta} L \Psi_{\eta}]| \ge u/8 \right).
$$

As for the trace terms, they are taken care of by Lemma A.10:

$$
\mathbb{P} \left( |\operatorname{Tr}(\Psi'_{\eta} L \Psi_{\varepsilon}) - \mathbb{E}[\Psi'_{\varepsilon} L \Psi_{\varepsilon}]| \ge u/8 \right) \le c_3 \exp \left( -c_4 \frac{(u/8)^2 T p^2}{\|\Gamma_h\|_2^2} \right)
$$

$$
\mathbb{P} \left( |\operatorname{Tr}(\Psi'_{\eta} L \Psi_{\eta}) - \mathbb{E}[\Psi'_{\eta} L \Psi_{\eta}]| \ge u/8 \right) \le c_3 \exp \left( -c_4 \frac{(u/8)^2 T p^2}{\omega^4} \right)
$$
<table><caption>Table 1: Notation correspondence between this paper and Han et al. [2015]</caption><thead><tr><th></th><th>This paper</th><th>Han et al. [2015]</th></tr></thead><tbody><tr><td>VAR def</td><td>X<sub>t</sub> = &theta;X<sub>t-1</sub> + &epsilon;<sub>t</sub></td><td>X<sub>t</sub> = A'<sub>1</sub>X<sub>t-1</sub> + Z<sub>t</sub></td></tr><tr><td>Covariance</td><td>&Gamma;<sub>h</sub> = Cov(X<sub>h</sub>, X<sub>0</sub>)</td><td>&Sigma;<sub>i</sub> = Cov(X<sub>0</sub>, X<sub>i</sub>)</td></tr><tr><td>Yule-Walker</td><td>&Gamma;<sub>h</sub> = &theta;<sup>h</sup>&Gamma;<sub>0</sub></td><td>&Sigma;<sub>i</sub> = &Sigma;<sub>0</sub>A<sub>1</sub><sup>i</sup></td></tr><tr><td>Covariance estimate</td><td>&Gamma;&#770;<sub>h</sub></td><td>S<sub>i</sub></td></tr><tr><td>Covariance error</td><td>err(&delta;)</td><td>&zeta;<sub>i</sub></td></tr><tr><td>Optimization constraint</td><td>||M&Gamma;&#770;<sub>0</sub> - &Gamma;&#770;<sub>1</sub>||<sub>max</sub> &le; &lambda;</td><td>||S<sub>0</sub>M - S<sub>1</sub>||<sub>max</sub> &le; &lambda;</td></tr><tr><td>Optimization objective</td><td>||vec(M)||<sub>1</sub></td><td>||vec(M)||<sub>1</sub></td></tr><tr><td>Threshold in proof</td><td>&nu;</td><td>&lambda;<sub>1</sub></td></tr></tbody></table>

Under these hypotheses, we just proved that with probability at least $1 - 8\delta$,

$$ |(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \le c \frac{\sigma_{\max}^2 + \omega^2}{(1 - \vartheta)^2} \frac{\sqrt{\log(1/\delta)}}{\sqrt{Tpq_u}}. $$

We finish with a union bound, applying the previous result to all pairs $(d_1, d_2) \in [D]^2$. With probability greater than $1 - 8D^2\delta$, we have:

$$ \max_{d_1, d_2} |(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| = \|\hat{\Gamma}_h - \Gamma_h\|_{\max} \le c \frac{\sigma_{\max}^2 + \omega^2}{(1 - \vartheta)^2} \frac{\sqrt{\log(1/\delta)}}{\sqrt{Tpq_u}}. $$

Replacing $\delta$ with $8D^2\delta$ gives us the result we wanted: with probability greater than $1 - \delta$,

$$ \|\hat{\Gamma}_h - \Gamma_h\|_{\max} \le c \frac{\sigma_{\max}^2 + \omega^2}{(1 - \vartheta)^2} \frac{\sqrt{\log(D/\delta)}}{\sqrt{Tpq_u}}. $$

□

## A.7 BEHAVIOR OF THE DANTZIG SELECTOR

We now walk the final steps from the error on $\hat{\Gamma}_h$ to the error on $\hat{\theta}$. In order to recover Theorem 1, we adapt the convergence proof from Han et al. [2015, Appendix A.1]. However, we use our own notations and our custom concentration results for $\hat{\Gamma}_h$. To make comparison between both papers easier, we provide a dictionary of the main notations in Table 1.

Our sparse transition estimator is defined as a solution to (8). The end goal is to control the error $\|\hat{\theta} - \theta\|_1$, where $\theta = \Gamma_1\Gamma_0^{-1}$ is the true transition matrix. We start by choosing a specific $\lambda$ such that $\theta$ is feasible with high probability.

**Lemma 14 (Feasibility of the real $\theta$).** If we select the penalization level

$$ \lambda = (\|\theta\|_{\infty} + 1) \text{err}(\delta), $$

then with probability at least $1 - \delta$, the real $\theta$ is a feasible solution to the optimization problem (8).

*Proof.*

$$ \begin{align*} \|\theta\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} &= \|\Gamma_1\Gamma_0^{-1}\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} \\ &= \|\Gamma_1\Gamma_0^{-1}\hat{\Gamma}_0 - \Gamma_1 + \Gamma_1 - \hat{\Gamma}_1\|_{\max} \\ &\le \|\Gamma_1\Gamma_0^{-1}\hat{\Gamma}_0 - \Gamma_1\Gamma_0^{-1}\Gamma_0\|_{\max} + \|\Gamma_1 - \hat{\Gamma}_1\|_{\max} \\ &= \|\theta(\hat{\Gamma}_0 - \Gamma_0)\|_{\max} + \|\Gamma_1 - \hat{\Gamma}_1\|_{\max} \end{align*} $$
We plug this in and rearrange to get:

$$ \Delta_{\varepsilon\varepsilon} + \Delta_{\eta\varepsilon} + \Delta_{\varepsilon\eta} + \Delta_{\eta\eta} \le 4\delta + c_1 \exp\left(-\frac{c_2 u^2 T p^2}{\max\{(\|\Psi_{\varepsilon}\|_2^2 + \omega^2)^2, \|\Gamma_h\|_2^2, \omega^4\}}\right). \quad \square $$

The following result will simplify the denominator inside the exponential.

**Lemma A.12 (Spectral norms of $\Psi_{\varepsilon}$ and $\Gamma_h$).** The matrices $\Psi_{\varepsilon}$ and $\Gamma_h$ satisfy $\|\Psi_{\varepsilon}\|_2^2 \le \frac{\sigma^2}{(1-\vartheta)^2}$ and $\|\Gamma_h\|_2 \le \frac{\vartheta^h \sigma^2}{1-\vartheta}$. As a consequence,

$$ \max \{ (\|\Psi_{\varepsilon}\|_2^2 + \omega^2)^2, \|\Gamma_h\|_2^2, \omega^4 \} \le \frac{(\sigma^2 + \omega^2)^2}{(1 - \vartheta)^4} $$

*Proof.* By Lemma A.1, we can write $\Psi_{\varepsilon}^2$ as a sum of Kronecker products $\otimes$ (one for each block). Let $J_t$ be a matrix full of zeros, except for the subdiagonal of rank $t$, which is full of ones. Then we have:

$$ \Psi_{\varepsilon}^2 = \mathrm{Cov}[X] = I \otimes \Gamma_0(\theta) + \sum_{t=1}^{T-1} [J_t \otimes \theta^t \Gamma_0(\theta) + J'_t \otimes \Gamma_0(\theta) \theta'''] $$

This gives us control over its spectral norm thanks to Lemma C.2:

$$ \begin{aligned} \|\Psi_{\varepsilon}\|_{2}^{2} &= \|\Psi_{\varepsilon}^{2}\|_{2} \le \|I\|_{2} \times \|\Gamma_{0}(\theta)\|_{2} + \sum_{t=1}^{T-1} \left[ \|J_{t}\|_{2} \times \|\theta^{t}\Gamma_{0}(\theta)\|_{2} + \|J'_{t}\|_{2} \times \|\Gamma_{0}(\theta)\theta''\|_{2} \right] \\ &\le \|\Gamma_{0}(\theta)\|_{2} \left( 1 + 2 \sum_{t=1}^{T-1} \|\theta\|_{2}^{t} \right) \end{aligned} $$

We now use Lemma A.2:

$$ \|\Psi_{\varepsilon}\|_2^2 \le \frac{\sigma^2}{1 - \vartheta^2} \left(1 + 2 \frac{\vartheta}{1 - \vartheta}\right) = \frac{\sigma^2}{(1 - \vartheta)^2}. $$

We now turn to $\Gamma_h$ with Lemmas A.1 and A.2:

$$ \|\Gamma_h\|_2 = \|\theta^h \Gamma_0(\theta)\|_2 \le \frac{\vartheta^h \sigma^2}{1 - \vartheta^2}. $$

In particular, we have

$$ \begin{aligned} \max \{ (\|\Psi_{\varepsilon}\|_2^2 + \omega^2)^2, \|\Gamma_h\|_2^2, \omega^4 \} &\le \max \left\{ \left( \frac{\sigma^2}{(1-\vartheta)^2} + \omega^2 \right)^2, \left( \frac{\vartheta^h \sigma^2}{1-\vartheta^2} \right)^2, \omega^4 \right\} \\ &\le \frac{(\sigma^2 + \omega^2)^2}{(1-\vartheta)^4} \quad \square \end{aligned} $$

We can now control the error of the covariance estimator:

**Lemma A.13 (Max norm convergence rate of the covariance estimator).** Let $\delta > 0$ be small enough. Assume that Equations (19) and (22) hold. Then the covariance estimator $\hat{\Gamma}_h$ from Equation (9) satisfies

$$ \|\hat{\Gamma}_h - \Gamma_h\|_{\max} \le c \frac{\sigma^2 + \omega^2 \sqrt{\log(D/\delta)}}{(1-\vartheta)^2 \sqrt{Tpq_u}} = \mathrm{err}(\delta) $$

with probability greater than $1 - \delta$.
By Lemma 26,

$$ \|\theta(\hat{\Gamma}_0 - \Gamma_0)\|_{\max} \le \|\theta\|_{\infty} \|\hat{\Gamma}_0 - \Gamma_0\|_{\max} $$

By Lemma 13, with probability greater than $1 - 2\delta$,

$$ \|\hat{\Gamma}_0 - \Gamma_0\|_{\max} \le \mathrm{err}(\delta) \quad \text{and} \quad \|\hat{\Gamma}_1 - \Gamma_1\|_{\max} \le \mathrm{err}(\delta) $$

which implies

$$ \|\theta\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} \le (\|\theta\|_{\infty} + 1) \mathrm{err}(\delta). $$

This is exactly the feasibility criterion for (8) if $\lambda = (\|\theta\|_{\infty} + 1) \mathrm{err}(\delta)$. $\square$

**Lemma 15** (Error on $\hat{\theta}$ in max norm). If we select $\lambda = (\|\theta\|_{\infty} + 1) \mathrm{err}(\delta)$, then with probability at least $1 - \delta$, the max norm error of $\hat{\theta}$ satisfies

$$ \|\hat{\theta} - \theta\|_{\max} \le 2\lambda \|\Gamma_0^{-1}\|_1. $$

*Proof.*

$$ \begin{aligned} \|\hat{\theta} - \theta\|_{\max} &= \|\hat{\theta} - \Gamma_1 \Gamma_0^{-1}\|_{\max} \\ &= \|(\hat{\theta}\Gamma_0 - \Gamma_1)\Gamma_0^{-1}\|_{\max} \\ &= \|(\hat{\theta}\Gamma_0 - \hat{\theta}\hat{\Gamma}_0 + \hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1 + \hat{\Gamma}_1 - \Gamma_1)\Gamma_0^{-1}\|_{\max} \\ &\le \|(\hat{\theta}\Gamma_0 - \hat{\theta}\hat{\Gamma}_0)\Gamma_0^{-1}\|_{\max} + \|(\hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1)\Gamma_0^{-1}\|_{\max} + \|(\hat{\Gamma}_1 - \Gamma_1)\Gamma_0^{-1}\|_{\max} \end{aligned} $$

By Lemma 26,

$$ \begin{aligned} \|\hat{\theta} - \theta\|_{\max} &\le (\|\hat{\theta}(\Gamma_0 - \hat{\Gamma}_0)\|_{\max} + \|\hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} + \|\hat{\Gamma}_1 - \Gamma_1\|_{\max}) \|\Gamma_0^{-1}\|_1 \\ &\le (\|\hat{\theta}\|_{\infty} \|\Gamma_0 - \hat{\Gamma}_0\|_{\max} + \|\hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} + \|\hat{\Gamma}_1 - \Gamma_1\|_{\max}) \|\Gamma_0^{-1}\|_1 \end{aligned} $$

We want to control $\|\hat{\theta}\|_{\infty}$ using $\|\theta\|_{\infty}$. Let us recall that the operator $\ell_{\infty}$ norm is equal to the maximum $\ell_1$ norm of the rows of a matrix. To control the rows of $\hat{\theta}$, we notice that the optimization problem defining $\hat{\theta}$, namely

$$ \min_{M \in \mathbb{R}^{D \times D}} \|\mathrm{vec}(M)\|_1 \quad \text{s.t.} \quad \|M\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} \le \lambda $$

is equivalent to the row-wise minimization

$$ \forall i, \quad \min_{M_{i, \cdot} \in \mathbb{R}^{1 \times D}} \|M_{i, \cdot}\|_1 \quad \text{s.t.} \quad \|M_{i, \cdot}\hat{\Gamma}_0 - (\hat{\Gamma}_1)_{i, \cdot}\|_{\max} \le \lambda $$

From this, we deduce that each row of the optimum $\hat{\theta}$ satisfies $\|\hat{\theta}_{i, \cdot}\|_1 \le \|\theta_{i, \cdot}\|_1$, which implies $\|\hat{\theta}\|_{\infty} \le \|\theta\|_{\infty}$. Going back to our error estimate, we get:

$$ \|\hat{\theta} - \theta\|_{\max} \le (\|\theta\|_{\infty} \|\Gamma_0 - \hat{\Gamma}_0\|_{\max} + \|\hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} + \|\hat{\Gamma}_1 - \Gamma_1\|_{\max}) \|\Gamma_0^{-1}\|_1 $$

Note that the middle term is smaller than $\lambda$ because the optimum $\hat{\theta}$ is a feasible solution. Meanwhile, the first and third term are smaller than $\mathrm{err}(\delta)$ with probability $1 - \delta$:

$$ \|\hat{\theta} - \theta\|_{\max} \le (\|\theta\|_{\infty} \mathrm{err}(\delta) + \lambda + \mathrm{err}(\delta)) \|\Gamma_0^{-1}\|_1 = 2\lambda \|\Gamma_0^{-1}\|_1 $$

To complete the proof of Theorem 1, we simply need to go from the max norm to the $\ell_{\infty}$ operator norm.
Proof. Let us plug Lemma A.12 into Lemma A.11

$$ \mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \ge u) \le 4\delta + c_1 \exp\left(-\frac{c_2(1-\vartheta)^4 u^2 T p^2}{(\sigma^2 + \omega^2)^2}\right). $$

All that is left to do is choose $u$ such that $\mathbb{P}(|(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \ge u) \le 8\delta$, which will be true if

$$ c_1 \exp\left(-\frac{c_2(1-\vartheta)^4 T p^2}{(\sigma^2 + \omega^2)^2} u^2\right) \le 4\delta \iff u \ge \sqrt{\frac{\log(c_1/4\delta)(\sigma^2 + \omega^2)^2}{c_2(1-\vartheta)^4 T p^2}}. $$

As long as $\delta$ is small enough, we can take

$$ u = c \frac{\sqrt{\log(1/\delta)(\sigma^2 + \omega^2)}}{(1 - \vartheta)^2 \sqrt{T p^2}}. \quad (21) $$

For Lemma A.11 to apply, we must verify that $u \in [0, 1]$ and that Equation (20) is satisfied. In other words, we have to ensure that

$$ c \frac{\sqrt{\log(1/\delta)(\sigma^2 + \omega^2)}}{(1 - \vartheta)^2 \sqrt{T p^2}} \le \min\{1, 2(\|\Psi_\epsilon\|_2^2 + \omega^2)\} $$

Using Lemma A.12, this is implied by the condition

$$ \frac{\sqrt{\log(1/\delta)} \max\{1, (\sigma^2 + \omega^2)^{-1}\}}{(1 - \vartheta)^2 \sqrt{T p^2}} \le c \quad (22) $$

Under these hypotheses, we just proved that with probability at least $1 - 8\delta$,

$$ |(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| \le c \frac{\sigma^2 + \omega^2}{(1 - \vartheta)^2} \frac{\sqrt{\log(1/\delta)}}{\sqrt{T p^2}}. $$

We finish with a union bound, applying the previous result to all pairs $(d_1, d_2) \in [D]^2$. With probability greater than $1 - 8D^2\delta$, we have:

$$ \max_{d_1, d_2} |(\hat{\Gamma}_h - \Gamma_h)_{d_1, d_2}| = \|\hat{\Gamma}_h - \Gamma_h\|_{\max} \le c \frac{\sigma^2 + \omega^2}{(1 - \vartheta)^2} \frac{\sqrt{\log(1/\delta)}}{\sqrt{T p^2}}. $$

Replacing $\delta$ with $8D^2\delta$ gives us the result we wanted: with probability greater than $1 - \delta$,

$$ \|\hat{\Gamma}_h - \Gamma_h\|_{\max} \le c \frac{\sigma^2 + \omega^2}{(1 - \vartheta)^2} \frac{\sqrt{\log(D/\delta)}}{\sqrt{T p^2}}. \quad \square $$

### A.7. Behavior of the Dantzig selector

We now walk the final steps from the error on $\hat{\Gamma}_h$ to the error on $\hat{\theta}$. In order to recover Theorem 3.1, we adapt the convergence proof from Han, Lu and Liu [17, Appendix A.1]. However, we use our own notations and our custom concentration results for $\hat{\Gamma}_h$. To make comparison between both papers easier, we provide a dictionary of the main notations in Table 1.

Our sparse transition estimator is defined as a solution to (8). The end goal is to control the error $\|\hat{\theta} - \theta\|_1$, where $\theta = \Gamma_1 \Gamma_0^{-1}$ is the true transition matrix. We start by choosing a specific $\lambda$ such that $\theta$ is feasible with high probability.
Proof. Let $\nu > 0$ be a threshold (to be chosen later). We define

$$s_1 = \max_i \sum_j \min \left\{ \frac{|\theta_{i,j}|}{\nu}, 1 \right\} \quad \text{and} \quad \mathcal{I}_i = \{j : |\theta_{i,j}| \ge \nu\}$$

With high probability, the following holds for any row $i$:

$$\begin{aligned}
\|\hat{\theta}_{i, \cdot} - \theta_{i, \cdot}\|_1 &\le \|\hat{\theta}_{i, \mathcal{I}_i^c} - \theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\
&\le \|\hat{\theta}_{i, \mathcal{I}_i^c}\|_1 + \|\theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\
&= (\|\hat{\theta}_{i, \cdot}\|_1 - \|\hat{\theta}_{i, \mathcal{I}_i}\|_1) + \|\theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\
&\le \|\theta_{i, \cdot}\|_1 - \|\hat{\theta}_{i, \mathcal{I}_i}\|_1 + \|\theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\
&= (\|\theta_{i, \mathcal{I}_i}\|_1 + \|\theta_{i, \mathcal{I}_i^c}\|_1) - \|\hat{\theta}_{i, \mathcal{I}_i}\|_1 + \|\theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\
&= 2\|\theta_{i, \mathcal{I}_i^c}\|_1 + (\|\theta_{i, \mathcal{I}_i}\|_1 - \|\hat{\theta}_{i, \mathcal{I}_i}\|_1) + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\
&\le 2\|\theta_{i, \mathcal{I}_i^c}\|_1 + 2\|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1
\end{aligned}$$

By definition of $\mathcal{I}_i$, for all $j \in \mathcal{I}_i^c$, $|\theta_{i,j}| \le \nu$, hence

$$\|\theta_{i, \mathcal{I}_i^c}\|_1 = \sum_{j \in \mathcal{I}_i^c} |\theta_{i,j}| = \sum_{j \in \mathcal{I}_i^c} \min\{|\theta_{i,j}|, \nu\} \le \sum_j \min\{|\theta_{i,j}|, \nu\} \le \nu s_1$$

Meanwhile, the second term satisfies

$$\|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \le |\mathcal{I}_i| \times \|\hat{\theta} - \theta\|_{\max}$$

And by definition of $\mathcal{I}_i$, for all $j \in \mathcal{I}_i$, $|\theta_{i,j}| \ge \nu$, hence

$$|\mathcal{I}_i| = \sum_{j \in \mathcal{I}_i} 1 = \sum_{j \in \mathcal{I}_i} \min \left\{ \frac{|\theta_{i,j}|}{\nu}, 1 \right\} \le \sum_j \min \left\{ \frac{|\theta_{i,j}|}{\nu}, 1 \right\} \le s_1$$

Combining all of this, we get that with high probability,

$$\|\hat{\theta}_{i, \cdot} - \theta_{i, \cdot}\|_1 \le 2(\nu + 2\lambda \|\Gamma_0^{-1}\|_1) s_1$$

Judging by the last Equation, it makes sense to choose $\nu = 2\lambda \|\Gamma_0^{-1}\|_1$. Furthermore, our sparsity hypothesis on $\theta$ implies that for all but $s$ of the coefficients of any row $i$, $\min\{|\theta_{i,j}|, \nu\} = |\theta_{i,j}| = 0$. We deduce that for every $i$,

$$\sum_j \min\{|\theta_{i,j}|, \nu\} \le s \max_j \min\{|\theta_{i,j}|, \nu\} \le \nu s$$

which directly implies

$$\nu s_1 = \max_i \sum_j \min\{|\theta_{i,j}|, \nu\} \le \nu s$$

We finally find that with high probability,

$$\|\hat{\theta}_{i, \cdot} - \theta_{i, \cdot}\|_1 \le 4\nu s_1 \le 4\nu s = 8\lambda \|\Gamma_0^{-1}\|_1 s$$

With the help of a union bound, again with high probability,

$$\|\hat{\theta} - \theta\|_{\infty} = \max_i \|\hat{\theta}_{i, \cdot} - \theta_{i, \cdot}\|_1 \le 8\lambda \|\Gamma_0^{-1}\|_1 s$$

We substitute the value of $\lambda$ and obtain

$$\|\hat{\theta} - \theta\|_{\infty} \le 8(\|\theta\|_{\infty} + 1) \text{err}(\delta) \|\Gamma_0^{-1}\|_1 s$$
<table><caption>Table 1
Notation correspondence between this paper and Han, Lu and Liu [17]</caption><thead><tr><th></th><th>This paper</th><th>Han, Lu and Liu [17]</th></tr></thead><tbody><tr><td>VAR def</td><td>X<sub>t</sub> = &theta;X<sub>t-1</sub> + &epsilon;<sub>t</sub></td><td>X<sub>t</sub> = A'<sub>1</sub>X<sub>t-1</sub> + Z<sub>t</sub></td></tr><tr><td>Covariance</td><td>&Gamma;<sub>h</sub> = Cov(X<sub>h</sub>, X<sub>0</sub>)</td><td>&Sigma;<sub>i</sub> = Cov(X<sub>0</sub>, X<sub>i</sub>)</td></tr><tr><td>Yule-Walker</td><td>&Gamma;<sub>h</sub> = &theta;<sup>h</sup>&Gamma;<sub>0</sub></td><td>&Sigma;<sub>i</sub> = &Sigma;<sub>0</sub>A<sub>1</sub><sup>i</sup></td></tr><tr><td>Covariance estimate</td><td>&Gamma;&#770;<sub>h</sub></td><td>S<sub>i</sub></td></tr><tr><td>Covariance error</td><td>err(&delta;)</td><td>&zeta;<sub>i</sub></td></tr><tr><td>Optimization constraint</td><td>||M&Gamma;&#770;<sub>0</sub> - &Gamma;&#770;<sub>1</sub>||<sub>max</sub> &le; &lambda;</td><td>||S<sub>0</sub>M - s<sub>1</sub>||<sub>max</sub> &le; &lambda;</td></tr><tr><td>Optimization objective</td><td>||vec(M)||<sub>1</sub></td><td>||vec(M)||<sub>1</sub></td></tr><tr><td>Threshold in proof</td><td>&nu;</td><td>&lambda;<sub>1</sub></td></tr></tbody></table>

**Lemma A.14** (Feasibility of the real $\theta$). If we select the penalization level $\lambda = (\|\theta\|_{\infty} + 1) \text{err}(\delta)$, then with probability at least $1-\delta$, the real $\theta$ is a feasible solution to the optimization problem (8).

*Proof.*

$$
\begin{aligned}
\|\theta\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} &= \|\Gamma_1\Gamma_0^{-1}\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} \\
&= \|\Gamma_1\Gamma_0^{-1}\hat{\Gamma}_0 - \Gamma_1 + \Gamma_1 - \hat{\Gamma}_1\|_{\max} \\
&\le \|\Gamma_1\Gamma_0^{-1}\hat{\Gamma}_0 - \Gamma_1\Gamma_0^{-1}\Gamma_0\|_{\max} + \|\Gamma_1 - \hat{\Gamma}_1\|_{\max} \\
&= \|\theta(\hat{\Gamma}_0 - \Gamma_0)\|_{\max} + \|\Gamma_1 - \hat{\Gamma}_1\|_{\max}
\end{aligned}
$$

By Lemma C.4, $\|\theta(\hat{\Gamma}_0 - \Gamma_0)\|_{\max} \le \|\theta\|_{\infty}\|\hat{\Gamma}_0 - \Gamma_0\|_{\max}$. By Lemma A.13, with probability greater than $1-2\delta$, $\|\hat{\Gamma}_0 - \Gamma_0\|_{\max} \le \text{err}(\delta)$ and $\|\hat{\Gamma}_1 - \Gamma_1\|_{\max} \le \text{err}(\delta)$, which implies $\|\theta\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} \le (\|\theta\|_{\infty} + 1) \text{err}(\delta)$. This is exactly the feasibility criterion for (8) if $\lambda = (\|\theta\|_{\infty} + 1) \text{err}(\delta)$. $\square$

**Lemma A.15** (Error on $\hat{\theta}$ in max norm). If we select $\lambda = (\|\theta\|_{\infty} + 1) \text{err}(\delta)$, then with probability at least $1-\delta$, the max norm error of $\hat{\theta}$ satisfies $\|\hat{\theta} - \theta\|_{\max} \le 2\lambda\|\Gamma_0^{-1}\|_1$.

*Proof.*

$$
\begin{aligned}
\|\hat{\theta} - \theta\|_{\max} &= \|\hat{\theta} - \Gamma_1 \Gamma_0^{-1}\|_{\max} \\
&= \|(\hat{\theta}\Gamma_0 - \Gamma_1)\Gamma_0^{-1}\|_{\max} \\
&= \|(\hat{\theta}\Gamma_0 - \hat{\theta}\hat{\Gamma}_0 + \hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1 + \hat{\Gamma}_1 - \Gamma_1)\Gamma_0^{-1}\|_{\max} \\
&\le \|(\hat{\theta}\Gamma_0 - \hat{\theta}\hat{\Gamma}_0)\Gamma_0^{-1}\|_{\max} + \|(\hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1)\Gamma_0^{-1}\|_{\max} + \|(\hat{\Gamma}_1 - \Gamma_1)\Gamma_0^{-1}\|_{\max}
\end{aligned}
$$

By Lemma C.4,

$$
\begin{aligned}
\|\hat{\theta} - \theta\|_{\max} &\le (\|\hat{\theta}(\Gamma_0 - \hat{\Gamma}_0)\|_{\max} + \|\hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} + \|\hat{\Gamma}_1 - \Gamma_1\|_{\max}) \|\Gamma_0^{-1}\|_1 \\
&\le (\|\hat{\theta}\|_{\infty} \|\Gamma_0 - \hat{\Gamma}_0\|_{\max} + \|\hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} + \|\hat{\Gamma}_1 - \Gamma_1\|_{\max}) \|\Gamma_0^{-1}\|_1
\end{aligned}
$$

We want to control $\|\hat{\theta}\|_{\infty}$ using $\|\theta\|_{\infty}$. Let us recall that the operator $\ell_{\infty}$ norm is equal to the maximum $\ell_1$ norm of the rows of a matrix. To control the rows of $\hat{\theta}$, we notice that the optimization problem defining $\hat{\theta}$, namely

$$ \min_{M \in \mathbb{R}^{D \times D}} \|\text{vec}(M)\|_1 \quad \text{s.t.} \quad \|M\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} \le \lambda $$
Once we plug in the value of err($\delta$), the resulting high-probability error bound reads

$$ \|\hat{\theta} - \theta\|_{\infty} \le c \frac{\|\theta\|_{\infty} + 1}{\|\Gamma_0^{-1}\|_1^{-1}} \frac{\sigma_{\max}^2 + \omega^2}{(1 - \vartheta)^2} \frac{s\sqrt{\log(D/\delta)}}{\sqrt{Tpq_u}} $$

Since $\vartheta$ only acted as an upper bound on $\|\theta\|_2$ in this proof, we can define

$$ \gamma_u(\theta) = \frac{\|\theta\|_{\infty} + 1}{(1 - \|\theta\|_2)^2} \frac{\sigma_{\max}^2 + \omega^2}{\|\Gamma_0^{-1}\|_1^{-1}} $$

to obtain the compressed expression

$$ \|\hat{\theta} - \theta\|_{\infty} \le c\gamma_u(\theta) \frac{s\sqrt{\log(D/\delta)}}{\sqrt{Tpq_u}}. $$

<figure><img src="image_10.png" /></figure>

# B PROOF OF THE MINIMAX LOWER BOUND

We now present the detailed proof of Theorem 2.

## B.1 OVERVIEW

Our argument is based on Fano's method, which we sum up in Lemma 27. For a detailed presentation, we refer the reader to Tsybakov [2008, Chapter 2]. Note that Wainwright [2019, Chapter 15] and Duchi [2019, Chapter 7] also offer good treatments of the subject.

Fano's method relies on choosing a set of parameters $\theta_0, \theta_1, ..., \theta_M$ satisfying two seemingly contradictory conditions: their induced distributions must be hard to distinguish, yet they must lie as far apart from one another as possible. In particular, the crucial requirement of Fano's method is a tight upper bound on the KL divergence between two distributions generated by different parameters $\theta_i$ and $\theta_0$. Taking the latter to be 0, we actually want to bound

$$ \frac{1}{M+1} \sum_{i=1}^{M} \mathrm{KL} \{ \mathbb{P}_{\theta_i}(\Pi, Y) \| \mathbb{P}_0(\Pi, Y) \} \le \max_i \mathrm{KL} \{ \mathbb{P}_{\theta_i}(\Pi, Y) \| \mathbb{P}_0(\Pi, Y) \} $$

By Lemma 28,

$$ \mathrm{KL} \{ \mathbb{P}_{\theta_i}(\Pi, Y) \| \mathbb{P}_0(\Pi, Y) \} = \mathrm{KL} \{ \mathbb{P}_{\theta_i}(\Pi) \| \mathbb{P}_0(\Pi) \} + \mathrm{E}_{\Pi} [ \mathrm{KL} \{ \mathbb{P}_{\theta_i}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi) \} ] $$

Since $\theta_i$ does not affect the distribution of the sampling process $\Pi$, the first term of the right-hand side is zero, and we will concentrate on the second term. First, we will upper-bound the random variable inside the expectation for a fixed realization of $\Pi$, and then we will average said bound over all possible projections.

We now give the structure of the argument in a coherent order, along with the most important intermediate results:

1. Compute the conditional covariance $\mathrm{Cov}_{\theta}[Y|\Pi]$ and decompose it into a constant term $Q_{\Pi}$ (corresponding to the independent case $\theta = 0$) plus a residual $R_{\Pi}(\theta)$ (Lemma 16).

2. Upper-bound the conditional KL divergence $\mathrm{KL}\{\mathbb{P}_{\theta}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi)\}$ using the “deviations from the identity” $\Delta_{\Pi}(\theta) = Q_{\Pi}^{-1/2} R_{\Pi}(\theta) Q_{\Pi}^{-1/2}$ (Lemma 17).

3. Control $\Delta_{\Pi}(\theta)$ using features of $R(\theta)$ scaled by sampling-related factors (Lemmas 18, 19 and 20).

4. Deduce an upper bound on the KL divergence $\mathrm{E}_{\Pi}[\mathrm{KL}\{\mathbb{P}_{\theta}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi)\}]$ (Lemma 21).

5. Apply Fano's method to a set of parameters $\theta_i$ constructed from a pruned binary hypercube of well-chosen radius.
is equivalent to the row-wise minimization

$$ \forall i, \quad \min_{M_{i, \cdot} \in \mathbb{R}^{1 \times D}} \|M_{i, \cdot}\|_1 \quad \text{s.t.} \quad \|M_{i, \cdot} \hat{\Gamma}_0 - (\hat{\Gamma}_1)_{i, \cdot}\|_{\max} \le \lambda $$

From this, we deduce that each row of the optimum $\hat{\theta}$ satisfies $\|\hat{\theta}_{i, \cdot}\|_1 \le \|\theta_{i, \cdot}\|_1$, which implies $\|\hat{\theta}\|_{\infty} \le \|\theta\|_{\infty}$. Going back to our error estimate, we get:

$$ \|\hat{\theta} - \theta\|_{\max} \le \left( \|\theta\|_{\infty} \|\Gamma_0 - \hat{\Gamma}_0\|_{\max} + \|\hat{\theta}\hat{\Gamma}_0 - \hat{\Gamma}_1\|_{\max} + \|\hat{\Gamma}_1 - \Gamma_1\|_{\max} \right) \|\Gamma_0^{-1}\|_1 $$

Note that the middle term is smaller than $\lambda$ because the optimum $\hat{\theta}$ is a feasible solution. Meanwhile, the first and third term are smaller than $\text{err}(\delta)$ with probability $1 - \delta$:

$$ \|\hat{\theta} - \theta\|_{\max} \le (\|\theta\|_{\infty} \text{err}(\delta) + \lambda + \text{err}(\delta)) \|\Gamma_0^{-1}\|_1 = 2\lambda \|\Gamma_0^{-1}\|_1 \quad \square $$

To complete the proof of Theorem 3.1, we simply need to go from the max norm to the $\ell_{\infty}$ operator norm.

*Proof.* Let $\nu > 0$ be a threshold (to be chosen later). We define

$$ s_1 = \max_i \sum_j \min \left\{ \frac{|\theta_{i,j}|}{\nu}, 1 \right\} \quad \text{and} \quad \mathcal{I}_i = \{j : |\theta_{i,j}| \ge \nu\} $$

With high probability, the following holds for any row $i$:

$$ \begin{aligned} \|\hat{\theta}_{i, \cdot} - \theta_{i, \cdot}\|_1 &\le \|\hat{\theta}_{i, \mathcal{I}_i^c} - \theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\ &\le \|\hat{\theta}_{i, \mathcal{I}_i^c}\|_1 + \|\theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\ &= (\|\hat{\theta}_{i, \cdot}\|_1 - \|\hat{\theta}_{i, \mathcal{I}_i}\|_1) + \|\theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\ &\le \|\theta_{i, \cdot}\|_1 - \|\hat{\theta}_{i, \mathcal{I}_i}\|_1 + \|\theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\ &= (\|\theta_{i, \mathcal{I}_i}\|_1 + \|\theta_{i, \mathcal{I}_i^c}\|_1) - \|\hat{\theta}_{i, \mathcal{I}_i}\|_1 + \|\theta_{i, \mathcal{I}_i^c}\|_1 + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\ &= 2\|\theta_{i, \mathcal{I}_i^c}\|_1 + (\|\theta_{i, \mathcal{I}_i}\|_1 - \|\hat{\theta}_{i, \mathcal{I}_i}\|_1) + \|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \\ &\le 2\|\theta_{i, \mathcal{I}_i^c}\|_1 + 2\|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \end{aligned} $$

By definition of $\mathcal{I}_i$, for all $j \in \mathcal{I}_i^c$, $|\theta_{i,j}| \le \nu$, hence

$$ \|\theta_{i, \mathcal{I}_i^c}\|_1 = \sum_{j \in \mathcal{I}_i^c} |\theta_{i,j}| = \sum_{j \in \mathcal{I}_i^c} \min\{|\theta_{i,j}|, \nu\} \le \sum_j \min\{|\theta_{i,j}|, \nu\} \le \nu s_1 $$

Meanwhile, the second term satisfies $\|\hat{\theta}_{i, \mathcal{I}_i} - \theta_{i, \mathcal{I}_i}\|_1 \le |\mathcal{I}_i| \times \|\hat{\theta} - \theta\|_{\max}$. And by definition of $\mathcal{I}_i$, for all $j \in \mathcal{I}_i$, $|\theta_{i,j}| \ge \nu$, hence

$$ |\mathcal{I}_i| = \sum_{j \in \mathcal{I}_i} 1 = \sum_{j \in \mathcal{I}_i} \min \left\{ \frac{|\theta_{i,j}|}{\nu}, 1 \right\} \le \sum_j \min \left\{ \frac{|\theta_{i,j}|}{\nu}, 1 \right\} \le s_1 $$

Combining all of this, we get that with high probability,

$$ \|\hat{\theta}_{i, \cdot} - \theta_{i, \cdot}\|_1 \le 2(\nu + 2\lambda \|\Gamma_0^{-1}\|_1)s_1 $$

Judging by the last Equation, it makes sense to choose $\nu = 2\lambda \|\Gamma_0^{-1}\|_1$. Furthermore, our sparsity hypothesis on $\theta$ implies that for all but $s$ of the coefficients of any row $i$, $\min\{|\theta_{i,j}|, \nu\} =$
## B.2 CHANGE OF NOTATIONS

For this part, we slightly modify the previous conventions: we now assume that all the rows of $\Pi_t$ that contain only zeros are removed. In other words, $\Pi_t$ is no longer the diagonal matrix diag($\pi_t$) but instead becomes a wide rectangular matrix with exactly one 1 per row and at most one 1 per column. We thus have $\Pi\Pi' = I$ unless all of the $\pi_{t,d}$ are zero, in which case the matrix $\Pi$ is empty, and so are the observations $Y$. Let us denote this very unlikely event by $E$, and its complement by $E^c$. If $\Pi$ is such that $E$ happens, we obviously have $\text{KL}\{\mathbb{P}_{\theta_i}(Y|\Pi) \parallel \mathbb{P}_0(Y|\Pi)\} = 0$, which means that

$$ \mathbb{E}_{\Pi}[\text{KL}\{\mathbb{P}_{\theta_i}(Y|\Pi) \parallel \mathbb{P}_0(Y|\Pi)\}] = \mathbb{E}_{\Pi}[\mathbf{1}_{E^c} \text{KL}\{\mathbb{P}_{\theta_i}(Y|\Pi) \parallel \mathbb{P}_0(Y|\Pi)\}] \quad (24) $$

For the beginning of the proof, we consider a fixed, non-empty realization of $\Pi$.

## B.3 COVARIANCE DECOMPOSITION

As we announced in the proof sketch, our reference parameter will be $\theta_0 = 0$, which is why it makes sense to express the conditional covariance of $Y$ as a deviation from the case without interactions. This is the aim of the following result.

**Lemma 16 (Conditional covariance decomposition).** The covariance matrix of $Y$ given $\Pi$ decomposes as

$$ \text{Cov}_{\theta}[Y|\Pi] = Q_{\Pi} + R_{\Pi}(\theta), $$

where $Q_{\Pi}$ is a constant term and $R_{\Pi}(\theta)$ is a residual which vanishes as $\theta \to 0$. They are defined as follows: the constant term is

$$ Q_{\Pi} = \Pi(\text{diag}_T \Sigma)\Pi' + \omega^2 I $$

whereas the residual equals

$$ R_{\Pi}(\theta) = \Pi R(\theta) \Pi' \quad \text{with} \quad R(\theta) = \begin{bmatrix} \theta \Gamma_0(\theta) \theta' & \Gamma_0(\theta) \theta'^1 & \Gamma_0(\theta) \theta'^2 & \dots \\ \theta^1 \Gamma_0(\theta) & \theta \Gamma_0(\theta) \theta' & \Gamma_0(\theta) \theta'^1 & \dots \\ \theta^2 \Gamma_0(\theta) & \theta^1 \Gamma_0(\theta) & \theta \Gamma_0(\theta) \theta' & \dots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}. $$

*Proof.* We use Equation (5) to see that the conditional distribution $\mathbb{P}_{\theta}(Y|\Pi)$ is a centered multivariate Gaussian with covariance

$$ \text{Cov}_{\theta}[Y|\Pi] = \omega^2 I + \Pi \text{Cov}_{\theta}[X] \Pi'. $$

We then use Lemma 1 to get an expression of $\text{Cov}_{\theta}[X]$ and deduce that its constant term (w.r.t to $\theta$) is a block-diagonal matrix filled with copies of $\Sigma$:

$$ \text{Cov}_{\theta}[Y|\Pi] = \omega^2 I + \Pi \text{diag}_T(\Sigma)\Pi' + \Pi (\text{Cov}_{\theta}[X] - \text{diag}_T(\Sigma))\Pi'. $$

Finally, we define $Q_{\Pi} = \omega^2 I + \Pi \text{diag}_T(\Sigma)\Pi'$, $R(\theta) = \text{Cov}_{\theta}[X] - \text{diag}_T(\Sigma)$ and $R_{\Pi}(\theta) = \Pi R(\theta) \Pi'$ to obtain the decomposition we announced. The diagonal blocks of $R(\theta)$ are easily computed by noticing that $\Gamma_0(\theta) - \Sigma = \theta \Gamma_0(\theta) \theta'$. $\square$

## B.4 FROM THE KL DIVERGENCE TO $\Delta_{\Pi}(\theta)$

Judging by Lemma 16, choosing a parameter $\theta$ close to 0 yields a conditional distribution for $Y$ whose covariance is close to $Q_{\Pi}$. In the next result, we translate this into a bound on the KL divergence between $\mathbb{P}_{\theta}(Y|\Pi)$ and $\mathbb{P}_0(Y|\Pi)$.

**Lemma 17.** Let us define the deviation from the identity:

$$ \Delta_{\Pi}(\theta) = Q_{\Pi}^{-1/2} R_{\Pi}(\theta) Q_{\Pi}^{-1/2}. $$

Then the conditional KL divergence is upper-bounded by:

$$ \text{KL}\{\mathbb{P}_{\theta}(Y|\Pi) \parallel \mathbb{P}_0(Y|\Pi)\} \le \frac{\|\Delta_{\Pi}(\theta)\|_F^2}{2(1 + \lambda_{\min}(\Delta_{\Pi}(\theta)))}. $$
$|\theta_{i,j}| = 0$. We deduce that for every $i$, $\sum_j \min\{|\theta_{i,j}|, \nu\} \le s \max_j \min\{|\theta_{i,j}|, \nu\} \le \nu s$ which directly implies $\nu s_1 = \max_i \sum_j \min\{|\theta_{i,j}|, \nu\} \le \nu s$.

We finally find that with high probability,

$$ \|\hat{\theta}_{i, \cdot} - \theta_{i, \cdot}\|_{1} \le 4\nu s_{1} \le 4\nu S = 8\lambda\|\Gamma_{0}^{-1}\|_{1}s $$

With the help of a union bound, again with high probability,

$$ \|\hat{\theta} - \theta\|_{\infty} = \max_{i} \|\hat{\theta}_{i, \cdot} - \theta_{i, \cdot}\|_{1} \le 8\lambda \|\Gamma_{0}^{-1}\|_{1}s $$

We substitute the value of $\lambda$ and obtain $\|\hat{\theta} - \theta\|_{\infty} \le 8(\|\theta\|_{\infty} + 1) \text{err}(\delta) \|\Gamma_{0}^{-1}\|_{1}s$. Once we plug in the value of $\text{err}(\delta)$, the resulting high-probability error bound reads

$$ \|\hat{\theta} - \theta\|_{\infty} \le c \frac{\|\theta\|_{\infty} + 1}{\|\Gamma_{0}^{-1}\|_{1}^{-1}} \frac{\sigma^2 + \omega^2}{(1 - \theta)^2} \frac{s\sqrt{\log(D/\delta)}}{\sqrt{Tp^2}} $$

Since $\theta$ only acted as an upper bound on $\|\theta\|_2$ in this proof, we can define

$$ \gamma_u(\theta) = \frac{\|\theta\|_{\infty} + 1}{(1 - \|\theta\|_2)^2} \frac{\sigma^2 + \omega^2}{\|\Gamma_0^{-1}\|_1^{-1}} $$

to obtain the compressed expression

$$ \|\hat{\theta} - \theta\|_{\infty} \le c \gamma_u(\theta) \frac{s\sqrt{\log(D/\delta)}}{\sqrt{Tp^2}}. $$

We finally leverage our additional assumption that $\theta$ commutes with its transpose. Then we have $\|\Gamma_0^{-1}(\theta)\|_1^{-1} = \|(\sigma^2(I - \theta\theta')^{-1})^{-1}\|_1^{-1} = \sigma^2\|I - \theta\theta'\|_1^{-1}$ and we can give a simpler expression of $\gamma_\ell$ and $\gamma_u(\theta)$:

$$ \gamma_u(\theta) = \frac{(\|\theta\|_{\infty} + 1)\|I - \theta\theta'\|_1}{(1 - \|\theta\|_2)^2} \frac{\sigma^2 + \omega^2}{\sigma^2}. \quad \square $$

## Appendix B: Proof of the minimax lower bound

We now present the detailed proof of Theorem 3.2.

### B.1. Overview

Our argument is based on Fano's method, which we sum up in Lemma C.5. For a detailed presentation, we refer the reader to Tsybakov [37, Chapter 2]. Note that Wainwright [39, Chapter 15] and Duchi [15, Chapter 7] also offer good treatments of the subject.

Fano's method relies on choosing a set of parameters $\theta_0, \theta_1, \dots, \theta_M$ satisfying two seemingly contradictory conditions: their induced distributions must be hard to distinguish, yet they must lie as far apart from one another as possible. In particular, the crucial requirement of Fano's method is a tight upper bound on the KL divergence between two distributions generated by different parameters $\theta_i$ and $\theta_0$. Taking the latter to be 0, we actually want to bound

$$ \frac{1}{M+1} \sum_{i=1}^{M} \mathrm{KL} \{ \mathbb{P}_{\theta_i}(\Pi, Y) \| \mathbb{P}_0(\Pi, Y) \} \le \max_i \mathrm{KL} \{ \mathbb{P}_{\theta_i}(\Pi, Y) \| \mathbb{P}_0(\Pi, Y) \} $$
Proof. The conditional KL divergence KL {$\mathbb{P}_{\theta}(Y|\Pi) \parallel \mathbb{P}_{0}(Y|\Pi)$} can be bounded using Lemma 30. Indeed, both conditional distributions are Gaussian and have the same expectation, and covariance matrices that are “close” in the following sense: by Lemma 16,

$$
\begin{align*}
\operatorname{Cov}_0(Y|\Pi) &= Q_{\Pi} = Q_{\Pi}^{1/2} (Q_{\Pi}^{1/2})' \\
\operatorname{Cov}_{\theta}(Y|\Pi) &= Q_{\Pi} + R_{\Pi}(\theta) = Q_{\Pi}^{1/2} \left( I + \underbrace{Q_{\Pi}^{-1/2} R_{\Pi}(\theta) Q_{\Pi}^{-1/2}}_{\Delta_{\Pi}(\theta)} \right) (Q_{\Pi}^{1/2})'.
\end{align*}
$$

By Lemma 23, there exists a real number $r_{\min} \ge s_{\min} (Q_{\Pi}^{1/2})^2 = s_{\min}(Q_{\Pi})$ such that

$$
\lambda_{\min}(\operatorname{Cov}_{\theta}(Y|\Pi)) = r_{\min}\lambda_{\min}(I + \Delta_{\Pi}(\theta)).
$$

Since $Q_{\Pi} \succeq \omega^2 I \succ 0$, its minimum singular value satisfies $s_{\min}(Q_{\Pi}) > 0$, so that $r_{\min} > 0$. In addition, $\operatorname{Cov}_{\theta}(Y|\Pi) \succeq \omega^2 I \succ 0$, so that $\lambda_{\min}(\operatorname{Cov}_{\theta}(Y|\Pi)) > 0$. Therefore,

$$
\lambda_{\min}(I + \Delta_{\Pi}(\theta)) = \frac{\lambda_{\min}(\operatorname{Cov}_{\theta}(Y|\Pi))}{r_{\min}} > 0 \quad \text{and} \quad \lambda_{\min}(\Delta_{\Pi}(\theta)) > -1,
$$

which means we can apply Lemma 30 with $\mathbb{P}_1 = \mathbb{P}_{\theta}(Y|\Pi)$ and $\mathbb{P}_0 = \mathbb{P}_0(Y|\Pi)$. $\square$

<figure><img src="image_11.png" /></figure>

B.5 FROM Δ<sub>Π</sub>(θ) TO R<sub>Π</sub>(θ)

Lemma 17 strongly suggests studying a certain fraction involving $\Delta_{\Pi}(\theta)$. In the following result, we boil it down to a function of the residual term $R_{\Pi}(\theta)$.

**Lemma 18.** *Assume* $\|R(\theta)\|_2 \le (\sigma_{\min}^2 + \omega^2)/2$. *We have the following upper bound:*

$$
\frac{\|\Delta_{\Pi}(\theta)\|_{F}^{2}}{2(1 + \lambda_{\min}(\Delta_{\Pi}(\theta)))} \le \frac{\|R_{\Pi}(\theta)\|_{F}^{2}}{(\sigma_{\min}^{2} + \omega^{2})^{2}}.
$$

*Proof.* Since the quantity $\lambda_{\min}(\Delta_{\Pi}(\theta))$ in the denominator is hard to control, we will work with the spectral norm instead. Indeed, whenever $\|\Delta_{\Pi}(\theta)\|_2 < 1$, we have the crude bound

$$
\frac{1}{1 - \lambda_{\min}(\Delta_{\Pi}(\theta))} \le \frac{1}{1 - \|\Delta_{\Pi}(\theta)\|_2}.
$$

Let us start by noticing that, thanks to Lemma 25,

$$
\begin{align*}
\|\Delta_{\Pi}(\theta)\|_{F}^{2} &= \|Q_{\Pi}^{-1/2} R_{\Pi}(\theta) Q_{\Pi}^{-1/2}\|_{F}^{2} \le \|Q_{\Pi}^{-1/2}\|_{2}^{4} \|R_{\Pi}(\theta)\|_{F}^{2} = \|Q_{\Pi}^{-1}\|_{2}^{2} \|R_{\Pi}(\theta)\|_{F}^{2} \\
\|\Delta_{\Pi}(\theta)\|_{2} &= \|Q_{\Pi}^{-1/2} \Pi R(\theta) \Pi' Q_{\Pi}^{-1/2}\|_{2} \le \|Q_{\Pi}^{-1/2} \Pi\|_{2}^{2} \|R(\theta)\|_{2}.
\end{align*}
$$

We will later see how the spectral and Frobenius norms of the residual $R(\theta)$ can be controlled as a function of $\theta$. For now, we must work to upper bound $\|Q_{\Pi}^{-1}\|_2$ and $\|Q_{\Pi}^{-1}\Pi\|_2^2$.

To simplify the following proof, we write $\Sigma_d = b\mathrm{diag}_T \Sigma$. Since $\Sigma_d$ is block-diagonal, its spectrum is the same as the spectrum of $\Sigma$ repeated $T$ times, hence $\lambda_{\min}(\Sigma_d) = \sigma_{\min}^2$. And since we assumed $E^c$ happens (non empty projection), we have $\Pi\Pi' = I$ and $\Pi'\Pi = \mathrm{diag}(\pi)$, which has at least one entry equal to 1.

We start with $\|Q_{\Pi}^{-1}\|_2$. Since $Q_{\Pi} \succeq \omega^2 I \succ 0$ is non-singular and symmetric,

$$
\|Q_{\Pi}^{-1}\|_2 = \lambda_{\max}(Q_{\Pi}^{-1}) = \frac{1}{\lambda_{\min}(Q_{\Pi})} = \frac{1}{\lambda_{\min}(\Pi\Sigma_d\Pi' + \omega^2 I)}.
$$

Remembering that $\Sigma_d \succeq \sigma_{\min}^2 I$, we get

$$
\Pi \Sigma_d \Pi' + \omega^2 I \succeq \sigma_{\min}^2 \Pi \Pi' + \omega^2 I = (\sigma_{\min}^2 + \omega^2) I
$$
By Lemma C.6,

$$KL \{\mathbb{P}_{\theta_i}(\Pi, Y) \| \mathbb{P}_0(\Pi, Y)\} = KL \{\mathbb{P}_{\theta_i}(\Pi) \| \mathbb{P}_0(\Pi)\} + \mathbb{E}_{\Pi} [KL \{\mathbb{P}_{\theta_i}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi)\}]$$

Since $\theta_i$ does not affect the distribution of the sampling process $\Pi$, the first term of the right-hand side is zero, and we will concentrate on the second term. First, we will upper-bound the random variable inside the expectation for a fixed realization of $\Pi$, and then we will average said bound over all possible projections.

We now give the structure of the argument in a coherent order, along with the most important intermediate results:

1. Compute the conditional covariance $\text{Cov}_{\theta}[Y|\Pi]$ and decompose it into a constant term $Q_{\Pi}$ (corresponding to the independent case $\theta = 0$) plus a residual $R_{\Pi}(\theta)$ (Lemma B.1).

2. Upper-bound the conditional KL divergence $KL \{\mathbb{P}_{\theta}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi)\}$ using the “derivations from the identity” $\Delta_{\Pi}(\theta) = Q_{\Pi}^{-1/2} R_{\Pi}(\theta) Q_{\Pi}^{-1/2}$ (Lemma B.2).

3. Control $\Delta_{\Pi}(\theta)$ using features of $R(\theta)$ scaled by sampling-related factors (Lemmas B.3, B.4 and B.5).

4. Deduce an upper bound on the KL divergence $\mathbb{E}_{\Pi}[KL \{\mathbb{P}_{\theta}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi)\}]$ (Lemma B.6).

5. Apply Fano's method to a set of parameters $\theta_i$ constructed from a pruned binary hypercube of well-chosen radius.

## B.2. Change of notations

For this part, we slightly modify the previous conventions: we now assume that all the rows of $\Pi_t$ that contain only zeros are removed. In other words, $\Pi_t$ is no longer the diagonal matrix diag($\pi_t$) but instead becomes a wide rectangular matrix with exactly one 1 per row and at most one 1 per column. We thus have $\Pi\Pi' = I$ (an identity matrix of random size) unless all the $\pi_{t,d}$ are zero, in which case the matrix $\Pi$ is empty, and so are the observations $Y$. Let us denote this very unlikely event by $E$, and its complement by $E^c$. If $\Pi$ is such that $E$ happens, we obviously have $KL \{\mathbb{P}_{\theta_i}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi)\} = 0$, which means that

$$\mathbb{E}_{\Pi} [KL \{\mathbb{P}_{\theta_i}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi)\}] = \mathbb{E}_{\Pi} [\mathbf{1}_{E^c} KL \{\mathbb{P}_{\theta_i}(Y|\Pi) \| \mathbb{P}_0(Y|\Pi)\}] \quad (23)$$

For the beginning of the proof, we consider a fixed, non-empty realization of $\Pi$.

## B.3. Covariance decomposition

As we announced in the proof sketch, our reference parameter will be $\theta_0 = 0$, which is why it makes sense to express the conditional covariance of $Y$ as a deviation from the case without interactions. This is the aim of the following result.

**Lemma B.1 (Conditional covariance decomposition).** The covariance matrix of $Y$ given $\Pi$ decomposes as

$$\text{Cov}_{\theta}[Y|\Pi] = Q_{\Pi} + R_{\Pi}(\theta),$$
and thus

$$ \|Q_{\Pi}^{-1}\|_2 \le \frac{1}{\sigma_{\min}^2 + \omega^2}. $$

We now continue with $\|Q_{\Pi}^{-1/2}\Pi\|_2^2$. By definition of the spectral norm,

$$ \|Q_{\Pi}^{-1/2}\Pi\|_2^2 = \lambda_{\max}(\Pi'Q_{\Pi}^{-1}\Pi) = \lambda_{\max}(\Pi'(\Pi\Sigma_d\Pi' + \omega^2 I)^{-1}\Pi). $$

Because matrix inversion is decreasing w.r.t. the Loewner order on positive semi-definite matrices,

$$ \begin{aligned} (\Pi\Sigma_d\Pi' + \omega^2 I)^{-1} &\le (\sigma_{\min}^2 + \omega^2)^{-1} I^{-1} \\ \Pi'(\Pi\Sigma_d\Pi' + \omega^2 I)^{-1}\Pi &\le \frac{1}{\sigma_{\min}^2 + \omega^2} \Pi'\Pi. \end{aligned} $$

It follows that

$$ \|Q_{\Pi}^{-1/2}\Pi\|_2^2 \le \frac{1}{\sigma_{\min}^2 + \omega^2} \lambda_{\max}(\Pi'\Pi) = \frac{1}{\sigma_{\min}^2 + \omega^2}. $$

The conclusion is within reach:

$$ \begin{aligned} \frac{\|\Delta_{\Pi}(\theta)\|_F^2}{1 + \lambda_{\min}(\Delta_{\Pi}(\theta))} &\le \frac{\|\Delta_{\Pi}(\theta)\|_F^2}{1 - \|\Delta_{\Pi}(\theta)\|_2} \\ &\le \frac{\|Q_{\Pi}^{-1}\|_2^2 \|R_{\Pi}(\theta)\|_F^2}{1 - \|Q_{\Pi}^{-1/2}\Pi\|_2^2 \|R(\theta)\|_2} \\ &\le \frac{\left(\frac{1}{\sigma_{\min}^2 + \omega^2}\right)^2 \|R_{\Pi}(\theta)\|_F^2}{1 - \frac{1}{\sigma_{\min}^2 + \omega^2} \|R(\theta)\|_2} \\ &\le \frac{2\|R_{\Pi}(\theta)\|_F^2}{(\sigma_{\min}^2 + \omega^2)^2} \end{aligned} $$

The last inequality is justified by our assumption $\|R(\theta)\|_2 \le (\sigma_{\min}^2 + \omega^2)/2$. Another consequence of this assumption is that

$$ \|\Delta_{\Pi}(\theta)\|_2 \le \|Q_{\Pi}^{-1/2}\Pi\|_2^2 \|R(\theta)\|_2 \le \frac{1}{\sigma_{\min}^2 + \omega^2} \frac{\sigma_{\min}^2 + \omega^2}{2} = \frac{1}{2} < 1 $$

which is sufficient for the first inequality to hold.

□

## B.6 FROM $R_{\Pi}(\theta)$ TO $R(\theta)$

As the previous Lemma underlines, the last step we need to get rid of the dependency in $\Pi$ is to study the average norm of $R_{\Pi}(\theta)$.

**Lemma 19.** Let $q_{\ell} = \max\{1-b, 2p-(1-b)\}$. Then

$$ \mathbb{E}_{\Pi} [\mathbf{1}_{E^c} \|R_{\Pi}(\theta)\|_F^2] \le p \operatorname{Tr}[R(\theta) \odot R(\theta)] + pq_{\ell} \|R(\theta)\|_F^2. $$

*Proof.* We first notice that for any matrix $A$,

$$ \begin{aligned} \mathbb{E}_{\Pi} [\mathbf{1}_{E^c} \|\Pi A \Pi'\|_F^2] &= \mathbb{E}_{\Pi} [\mathbf{1}_{E^c} \operatorname{Tr}[\Pi A \Pi' \Pi A' \Pi')] \\ &= \mathbb{E}_{\Pi} [\operatorname{Tr}[\operatorname{diag}(\pi)A \operatorname{diag}(\pi)A'] ] \\ &= \sum_{i,j} \mathbb{E}_{\Pi}[\pi_i \pi_j] A_{i,j}^2. \end{aligned} $$

We can apply this to $R_{\Pi}(\theta) = \Pi R(\theta) \Pi'$:

$$ \mathbb{E}_{\Pi} [\mathbf{1}_{E^c} \|R_{\Pi}(\theta)\|_F^2] = \sum_{i,j} \mathbb{E}_{\Pi}[\pi_i \pi_j] R(\theta)_{i,j}^2. $$
where $Q_Π$ is a constant term and $R_Π(θ)$ is a residual which vanishes as $θ \to 0$. They are defined as follows: the constant term is $Q_Π = (σ^2+ω^2)I$ whereas the residual equals $R_Π(θ) = ΠR(θ)Π'$ with

$$R(θ) = \begin{bmatrix} θΓ_0(θ)θ' & Γ_0(θ)θ'^1 & Γ_0(θ)θ'^2 & \dots \\ θ^1Γ_0(θ) & θΓ_0(θ)θ' & Γ_0(θ)θ'^1 & \dots \\ θ^2Γ_0(θ) & θ^1Γ_0(θ) & θΓ_0(θ)θ' & \dots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}$$

*Proof.* We use Equation (5) to see that the conditional distribution $P_θ(Y|Π)$ is a centered multivariate Gaussian with covariance $Cov_θ[Y|Π] = ω^2I + ΠCov_θ[X]Π'$. We then use Lemma A.1 to get an expression of $Cov_θ[X]$ and deduce that its constant term (w.r.t to $θ$) is a block-diagonal matrix filled with copies of $σ^2I$:

$$Cov_θ[Y|Π] = ω^2I + σ^2ΠΠ' + Π(Cov_θ[X] - σ^2I)Π'$$

Finally, we define $Q_Π = ω^2I + σ^2ΠΠ' = (ω^2 + σ^2)I$, as well as $R(θ) = Cov_θ[X] - σ^2I$ and $R_Π(θ) = ΠR(θ)Π'$ to obtain the decomposition we announced. The diagonal blocks of $R(θ)$ are easily computed by noticing that $Γ_0(θ) - σ^2I = θΓ_0(θ)θ'$. □

### B.4. From the KL divergence to $Δ_Π(θ)$

Judging by Lemma B.1, choosing a parameter $θ$ close to 0 yields a conditional distribution for $Y$ whose covariance is close to $Q_Π$. In the next result, we translate this into a bound on the KL divergence between $P_θ(Y|Π)$ and $P_0(Y|Π)$.

**Lemma B.2.** Let us define the deviation from the identity: $Δ_Π(θ) = Q_Π^{-1/2}R_Π(θ)Q_Π^{-1/2}$. Then the conditional KL divergence is upper-bounded by: KL{$P_θ(Y|Π) || P_0(Y|Π)} ≤ \frac{||Δ_Π(θ)||_F^2}{2(1+λ_{min}(Δ_Π(θ)))}$, where $λ_{min}$ denotes the minimum eigenvalue.

*Proof.* The conditional KL divergence $KL\{P_θ(Y|Π) || P_0(Y|Π)}$ can be bounded using Lemma C.8. Indeed, both conditional distributions are Gaussian and have the same expectation, and covariance matrices that are “close” in the following sense: by Lemma B.1,

$$ \begin{aligned} Cov_0(Y|Π) &= Q_Π = Q_Π^{1/2} (Q_Π^{1/2})' \\ Cov_θ(Y|Π) &= Q_Π + R_Π(θ) = Q_Π^{1/2} \underbrace{(I + Q_Π^{-1/2} R_Π(θ) Q_Π^{-1/2})}_{Δ_Π(θ)} (Q_Π^{1/2})' \end{aligned} $$

By Lemma C.1, there exists a real number $r_{min} ≥ s_{min} (Q_Π^{1/2})^2 = σ^2 + ω^2 > 0$ (where $s_{min}$ is the minimum singular value) such that $λ_{min}(Cov_θ(Y|Π)) = r_{min}λ_{min}(I + Δ_Π(θ))$. In addition, $Cov_θ(Y|Π) ≥ ω^2I > 0$, so that $λ_{min}(Cov_θ(Y|Π)) > 0$. Therefore, $λ_{min}(I+Δ_Π(θ)) = \frac{λ_{min}(Cov_θ(Y|Π))}{r_{min}} > 0$ and $λ_{min}(Δ_Π(θ)) > -1$, which means we can apply Lemma C.8 with $P_1 = P_θ(Y|Π)$ and $P_0 = P_0(Y|Π)$. □
The rest of the proof consists in plugging in the moments $E_{\Pi}[\pi_t \pi_j]$ from Lemma 4:

$$
\begin{align*}
\mathbb{E}_{\Pi} [\|R_{\Pi}(\theta)\|_{F}^{2}] &= \sum_{\substack{t_1, t_2, d_1, d_2 \\ (t_1, d_1) = (t_2, d_2)}} pR(\theta)_{(t_1, d_1), (t_2, d_2)}^{2} \\
&\quad + \sum_{\substack{t_1, t_2, d_1, d_2 \\ d_1 \neq d_2}} p^{2}R(\theta)_{(t_1, d_1), (t_2, d_2)}^{2} \\
&\quad + \sum_{\substack{t_1, t_2, d_1, d_2 \\ d_1 = d_2, t_1 \neq t_2}} (p^{2} + p(1-p)(1-a-b)^{|t_1-t_2|})R(\theta)_{(t_1, d_1), (t_2, d_2)}^{2}.
\end{align*}
$$

The sum in the last term can be crudely controlled as follows:

$$
\begin{align*}
\sum_{\substack{t_1, t_2, d \\ t_1 \neq t_2}} (1-a-b)^{|t_1-t_2|} R(\theta)_{(t_1,d), (t_2,d)}^2 &\le |1-a-b| \sum_{t_1 \neq t_2} \sum_d (R(\theta)_{[t_1,t_2]})_{d,d}^2 \\
&\le |1-a-b| \sum_{t_1 \neq t_2} \|R(\theta)_{[t_1,t_2]}\|_F^2 \\
&\le |1-a-b| \cdot \|R(\theta)\|_F^2
\end{align*}
$$

This yields a short, but probably suboptimal bound:

$$
\mathbb{E}_{\Pi}[\mathbf{1}_{E^c} \|R_{\Pi}(\theta)\|_{F}^{2}] \le p \operatorname{Tr}[R(\theta) \odot R(\theta)] + (p^2 + p(1-p)|1-a-b|) \|R(\theta)\|_{F}^{2}.
$$

In the previous part, we already saw that

$$
p + (1 - p)(1 - a - b) = 1 - b.
$$

Similarly, we obtain

$$
\begin{align*}
p + (1-p)(a+b-1) &= \frac{a}{a+b} + \frac{b}{a+b}(a+b-1) \\
&= \frac{a+ba+b^2-b}{a+b} = \frac{a(1+b)-b(1-b)}{a+b} \\
&= p(1+b) - (1-p)(1-b) = 2p - (1-b).
\end{align*}
$$

As a consequence,

$$
p + (1 - p) |1 - a - b| = \max\{1 - b, 2p - (1 - b)\} = q_{\ell},
$$

which yields the expected result.

<figure><img src="image_12.png" /></figure>

B.7 BOUNDING $R(\theta)$

Lemma 19 relates the bounds involving $R_{\Pi}(\theta)$ to features of the full residual $R(\theta)$, which we now study.

**Lemma 20.** The residual $R(\theta)$ satisfies the following inequalities:

$$
\begin{align*}
\|R(\theta)\|_2 &\le \frac{2\sigma_{\max}^2}{(1-\vartheta)^2} \|\theta\|_2 \\
\|R(\theta)\|_F^2 &\le \frac{2T\sigma_{\max}^4}{(1-\vartheta)^3} \|\theta\|_F^2 \\
\operatorname{Tr}[R(\theta) \odot R(\theta)] &\le \frac{T\sigma_{\max}^4}{(1-\vartheta)^2} \|\theta\|_2^2 \|\theta\|_F^2.
\end{align*}
$$