Guillaume Dalle
email: guillaume.dalle@enpc.fr

Yohann De Castro
email: yohann.de-castro@ec-lyon.fr

MINIMAX ESTIMATION OF PARTIALLY-OBSERVED VECTOR AUTOREGRESSIONS

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 statespace 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 nearoptimality 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.

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?

CONTEXT OF THE STUDY

To answer this question, we study a state-space model where the state X t ∈ 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 θ ∈ R D×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 π t,d = 1, where π 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 = θX t-1 + N (0, σ 2 I) π t,d ∼ B(p) Y t = diag(π t )X t + N (0, ω 2 I).

(1)

When we write π t,d ∼ 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 π 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 θ, whose nonasymptotic 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 θ satisfies

θ -θ ∞ 1 + ω 2 σ 2 s p √ T (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 ω 2 /σ 2 , which means that precise recovery of θ 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.

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 ∈ E. This congestion evolves according to a VAR process, whose transition matrix θ represents pairwise interactions between edges. The sparsity structure of θ 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 π 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 π 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 θ, which is essential to help railway operatives dimension their data sets or evaluate prediction uncertainty.

RELATED WORKS

The theory of VAR processes has been known for a long time: the book of [START_REF] Lütkepohl | New Introduction to Multiple Time Series Analysis[END_REF] provides a detailed account. If we have full and noiseless observations of the process X t , we can use conditional Least Squares to estimate θ by minimizing the quadratic error t X t -θX t-1 2 2 . This is equivalent to solving the Yule-Walker equation Γ h = θΓ h-1 , where we replace the autocovariance matrix Γ h = Cov[X t+h , X t ] with its empirical counterpart Γ 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 [START_REF] Buonaccorsi | Measurement Error: Models, Methods, and Applications[END_REF]. The same goes for incomplete data sets ; an extensive survey was recently published by [START_REF] Little | Statistical Analysis with Missing Data[END_REF]. According to their terminology, our work deals with data that is missing completely at random (MCAR), which means that the projection π t is independent from the underlying process X t . We also assume to know the distribution of the missingness indicators 1 -π 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 [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. Most of the time, exact or approximate inference is achievable using some version of the Kalman filter [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF] or particle methods [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF], whereas parameter estimation typically involves the Expectation-Maximization (EM) algorithm [START_REF] Shumway | An Approach to Time Series Smoothing and Forecasting Using the EM Algorithm[END_REF]. 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 [START_REF] Loh | High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity[END_REF], and numerous other papers followed [see for example [START_REF] Basu | Regularized estimation in sparse high-dimensional time series models[END_REF][START_REF] Kock | Oracle inequalities for high dimensional vector autoregressions[END_REF][START_REF] Melnyk | Estimating Structured Vector Autoregressive Models[END_REF][START_REF] Jalali | Missing Data in Sparse Transition Matrix Estimation for Sub-Gaussian Vector Autoregressive Processes[END_REF].

As an alternative to the LASSO, the Dantzig selector [START_REF] Candes | The Dantzig selector: Statistical estimation when p is much larger than n[END_REF] 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. [START_REF] Han | A Direct Estimation of High Dimensional Stationary Vector Autoregressions[END_REF] 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 best1 .

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. [START_REF] Rao | Fundamental estimation limits in autoregressive processes with compressive measurements[END_REF] 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 [START_REF] Rao | Fundamental estimation limits in autoregressive processes with compressive measurements[END_REF], 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.

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. Appendix D contains a summary of the main notations and symbols.

THE PARTIALLY-OBSERVED VAR PROCESS AND ITS SPARSE ESTIMATOR

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

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∈[T ] follows a stationary VAR process of order 1. This process has dimension D and the following recursive definition:

X t = θX t-1 + ε t with ε t ∼ N (0, Σ). (3)
Here θ ∈ R D×D is the transition matrix and Σ ∈ R D×D is the covariance matrix of the innovations (in the introduction, we assumed Σ = σ 2 I).

To ensure stationarity of the VAR process, we must constrain the spectral radius of θ to satisfy ρ(θ) < 1. Throughout the paper, we actually make the following (slightly stronger) assumption on the spectral norm of θ: there exists ϑ ∈ (0, 1) such that for all the values of θ we consider, θ 2 ≤ ϑ < 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 θ ∈ Θ s where

Θ s = {θ ∈ R D×D : θ 2 ≤ ϑ < 1 and ∀i, θ i,• 0 ≤ s}. (4)
We denote by σ 2 min = λ min (Σ) and σ 2 max = λ max (Σ) the minimum and maximum eigenvalues of the covariance matrix Σ.

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 π t . Then, independent Gaussian noise with variance ω 2 is added to these selected components, and we observe the result. If we denote by Π t = diag(π t ) the diagonal projection matrix, we have

Y t = Π t X t + η t with η t ∼ N (0, ω 2 I). (5)
An essential hypothesis we make is the mutual independence between our three sources of randomness: the innovations ε t , the projections π t and the observation noise η t .

A major feature of the present work is the non-deterministic selection of observed state components, that is, the fact that π 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 ∈ (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 = P(π 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 = 1-a a b 1-b . Its coefficients are chosen to make the chain stationary with invariant measure ( b a+b , 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 χ such that 0 < χ ≤ a, b ≤ 1 -χ < 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 N T in the resulting error bounds.

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):

Γ h (θ) = Cov θ [X t+h , X t ] = θΓ h-1 (θ) = θ h Γ 0 (θ) (6)
We can use it to define a simple two-step procedure:

1. For a given h 0 , build estimators Γ h0 and Γ h0+1 of the covariances Γ h0 and Γ h0+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 Γ h0 is singular):

θ dense = Γ h0+1 Γ † h0 . (7)
The problem with this procedure is that is does not guarantee sparsity of θ. To obtain a sparse result, we follow [START_REF] Han | A Direct Estimation of High Dimensional Stationary Vector Autoregressions[END_REF] and cast Equation ( 6) as a soft constraint enforcing proximity between Γ h0+1 and θ Γ h0 . This amounts to solving the following constrained optimization problem:

θ ∈ argmin M ∈R D×D vec(M ) 1 subject to M Γ h0 -Γ h0+1 max ≤ λ. (8)
Here vec(•) 1 denotes the sum of the absolute values of all the coefficients of a matrix, while • 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 Γ h .

THE COVARIANCE ESTIMATOR we use is a variant of the empirical covariance. Since Y t = Π t X t + η t where η t is zero-mean, a natural proxy for X t is obtained by inverting the sampling operator: X t = Π † t Y t . It would therefore seem logical to build an estimator of Γ 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) = E[π t+h π t ]. Intuitively, since X t+h 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 -ω 2 I. This correction becomes unnecessary for h ≥ 1 since the observation noise η t is independent across time.

In conclusion, we obtain the following covariance estimator:

Γ h = 1 S(h) 1 T -h T -h t=1 Π † t+h Y t+h Π † t Y t -1 {h=0} ω 2 I. (9)
The coefficients of the scaling matrix S(h) are computed in Lemma 4.

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"...

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 θ 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

γ u (θ) = θ ∞ + 1 (1 -θ 2 ) 2 σ 2 max + ω 2 Γ 0 (θ) -1 -1 1 and q u = min{p, 1 -b} ≤ p. ( 10
)
Then there is a value of λ such that the following upper bound holds with probability at least 1 -δ:

θ -θ ∞ ≤ c γ u (θ)s √ T pq u log(D/δ). ( 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 ζ, when we pick the best possible estimator θ and nature replies by choosing the worst possible parameter θ:

P(ζ) = inf θ sup θ∈Θs P θ θ -θ ∞ ≥ ζ . (12)
More precisely, we want to find a threshold ζ such that the probability of exceeding it is non-negligible, for instance P(ζ) ≥ 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

γ = (1 -ϑ) 3/2 σ 2 min + ω 2 σ 2 max and q = max{1 -b, 2p -(1 -b)} ≥ p. ( 13
)
Then the following minimax lower bound holds:

inf θ sup θ∈Θs P θ θ -θ ∞ ≥ c γ s √ T pq ≥ 1 2 . ( 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 P θ0 (Π, Y ) and P θ1 (Π, Y ), where θ 0 and θ 1 are sufficiently far apart. See Appendix B for more details.

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/ √ T is optimal.

THE SAMPLING PARAMETERS appear as 1/ √ pq u in the upper bound, whereas the lower bound scales as 1/ √ pq 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 ≥ p ≥ 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 = p). So in the case of independent sampling, 1/p is indeed the optimal rate. THE 2 NORM OF THE TRANSITION MATRIX plays opposite roles on each side. In the lower bound, 1 -ϑ = 1 -max θ∈Θs θ 2 appears in the numerator, whereas in the upper bound, 1 -θ 2 appears in the denominator. It is likely that these dependencies are suboptimal, but at least they are compatible with one another: as θ 2 → 1, that is, as the VAR process becomes unstable, the lower bound tends to 0 and the upper bound to +∞. This is a reflection of the fact that our proofs make heavy use of the distance between θ and the unit sphere, which means they become meaningless when θ gets too large.

THE VARIANCES Σ AND ω 2 are involved in γ for the lower bound, and in γ u (θ) for the upper bound. In both cases, the ratio γ 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 Σ, while the magnitude of Y t is related to the spectrum of Σ + ω 2 I.

If the latter is significantly larger than the former, recovering X t (and thus θ) is a hopeless endeavor.

To simplify the comparison, let us assume in this discussion that Σ = σ 2 I, and that θ commutes with its transpose. Then we have

Γ -1 0 (θ) -1 1 = σ 2 (I -θθ ) -1 -1 -1 1 = σ 2 I -θθ -1
1 , and we can give a simpler expression of γ and γ u (θ):

γ u (θ) = ( θ 1 + 1) I -θθ 1 (1 -ϑ) 2 σ 2 + ω 2 σ 2 γ = (1 -ϑ) 3/2 σ 2 + ω 2 σ 2 .
We recognize the same dependency in both bounds, namely γ ∝ 1 + σ 2 ω 2 . Lemma 38 gives a heuristic argument linking this functional form to the asymptotic behavior of the MLE.

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 = θ 1 X t-1 + θ 2 X t-2 + ... + θ K X t-K + ε t .
Then we can represent this as a VAR(1) process using augmented variables [START_REF] Lütkepohl | New Introduction to Multiple Time Series Analysis[END_REF]. Indeed, observe that defining X

t = X t X t-1 • • • X t-K+1 and ε t = ε t 0 • • • 0 yields X t = θ X t-1 + ε t with θ =        θ 1 θ 2 • • • θ K-1 θ K I D 0 • • • 0 0 0 I D 0 0 . . . . . . . . . . . . 0 0 • • • I D 0        .
Unfortunately, by this reasoning, the Markov sampling mechanism that generates Π t gives rise to a new distribution for Π t which is no longer part of the same family. Indeed, the augmented sampling process Π 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.

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 [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF], linear optimization problems were modeled using JuMP [START_REF] Dunning | JuMP: A Modeling Language for Mathematical Optimization[END_REF] and solved with the COIN-OR Clp solver [START_REF] Forrest | Coin-or/Clp: Release releases/1.17.7[END_REF]. 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.

DATA GENERATION

Simulating a partially-observed VAR process with known transition matrix θ allows us to compute the estimation error θ -θ ∞ and study the influence of parameters such as T , D, s, p, ω, etc. Real values for θ were drawn using independent standard Gaussian distributions for each coefficient, and then normalized to satisfy θ 2 = ϑ = 1 2 . To simplify comparison with the theoretical bounds, we used a diagonal innovation covariance Σ = σ 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 ω is known):

T = 10000 D = 5 σ = 1.0 ω = 0.1 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 λ is paramount: indeed, Theorem 1 is only valid for a specific value of λ (which is not known in practice, but we can hope to approximate this near-optimal choice).

A standard way to tune λ would be cross-validation. However, evaluating a choice of λ (and the resulting estimate θ) requires inferring the hidden state sequence X t from the observations Y t . If the projection matrices Π t were deterministic, the inference could be performed with Kalman filtering [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF], 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 λ, we suppose that the sparsity level of the real transition matrix θ is known. We then use this target sparsity s to guide a dichotomy search on λ, until we find a transition matrix estimate θ whose row sparsity level is sufficiently close to s.

RESULTS

The main results are presented on Figure 1. With the exception of 1f, all plots have the estimation error θ -θ ∞ 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 θ. When a straight line is added to a scatter plot, it is the result of a Theil-Sen regression [START_REF] Sen | Estimates of the Regression Coefficient Based on Kendall's Tau[END_REF] applied to the points of the same color: its slope is denoted by α in the legend.

Figure 1a confirms that the error decreases as 1/

√

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 ∝ 1/T , in which case the error would increase as

√

T instead of decreasing. Figure 1b exhibits three clearly identifiable regimes with respect to the noise variance. In the first one, corresponding to ω/σ 1, the error remains small and constant. Then, the error increases when ω/σ 1. In the third phase, corresponding to ω/σ 1, the error remains high and volatile. This is consistent with the theoretical dependency in 1 + ω 2 /σ 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 θ dense scales linearly with D, while its sparse counterpart θ 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 ∞ 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 γ u (θ) also depends on the sparsity level in complicated ways through θ, especially since the real values are renormalized to satisfy θ 2 = 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 ω to perform estimation. If this parameter is unknown, we can choose h 0 ≥ 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 ω = 0, effectively trading lower variance for a higher bias. 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 from our Theorems could be refined and brought together with a more careful theoretical analysis.

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.

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 Γ h by showing that it defines an unbiased estimator of Γ h (Lemmas 3 and 4).

3. Fixing two indices d 1 and d 2 , rewrite ( Γ h -Γ h ) d1,d2 using quadratic forms g a Ψ a LΨ 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 Ψ a LΨ b g b (Lemma 11).

6. Obtain a high-probability control on Γ h -Γ h max with a union bound (Lemma 13).

7. Deduce the error of θ from the error of Γ h0 and Γ h0+1 by drawing inspiration from [START_REF] Han | A Direct Estimation of High Dimensional Stationary Vector Autoregressions[END_REF] (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:

Γ 0 (θ) = Cov θ [X t ] = ∞ k=0 θ k Σθ k Γ h (θ) = Cov θ [X t+h , X t ] = θ h Γ 0 (θ).
Proof. We start by noting that according to Equation (3), the stacked vector X = (X t ) t∈[T ] follows a T D-dimensional centered multivariate Gaussian distribution. The covariance matrix of X t can be deduced from the recursion:

Γ 0 (θ) = Cov θ [X t ] = Cov θ [θX t-1 + ε t ] = θ Cov θ [X t-1 ]θ + Σ = θΓ 0 (θ)θ + Σ.
There is a unique stationary solution:

Γ 0 (θ) = ∞ k=0 θ k Σθ k .
The covariance matrix between X t+h and X t is obtained similarly:

Γ h (θ) = Cov θ [X t+h , X t ] = E[X t+h X t ] = E[(θX t+h-1 + ε t+h )X t ] = θ Cov θ [X t+h-1 , X t ] = θ h Cov θ [X t , X t ] = θ h Γ 0 (θ). And Cov θ [X t , X t+h ] = Cov θ [X t+h , X t ] .
In other words, we just proved that

Cov θ [X] =  θ T -1 Γ 0 (θ) Γ 0 (θ)       
The following result will come in handy later.

Lemma 2 (Norm of Γ 0 (θ)). The covariance matrix Γ 0 (θ) satisfies

Γ 0 (θ) 2 ≤ σ 2 max 1 -ϑ 2 Proof. By Lemma 1, Γ 0 (θ) 2 ≤ ∞ k=0 θ k Σθ k 2 ≤ ∞ k=0 θ k 2 Σ 2 θ k 2 = Σ 2 1 -θ 2 2 ≤ σ 2 max 1 -ϑ 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 ∈ {h 0 , h 0 + 1} = {0, 1}.

Lemma 3 (Bias of the covariance estimator). The estimator Γ h given by Equation (9) for the covariance matrix Γ h is unbiased.

Proof. First, let us remember that since Π t = diag(π t ) is diagonal and binary, we also have

Π † t = Π t = Π t . By Equation (5), (Π † t+h Y t+h )(Π † t Y t ) = Π † t+h (Π t+h X t+h + η t+h )(X t Π t + η t )Π † t = diag(π t+h ) (X t+h X t + X t+h η t + η t+h X t + η t+h η t ) diag(π t ). (15)
Taking the conditional expectation and removing the cross-product terms (by independence of X and Π), we get:

E[(Π † t+h Y t+h )(Π † t Y t ) |Π] = diag(π t+h ) (E[X t+h X t ] + E[η t+h η t ]) diag(π t ).
Since E [X t+h X t ] = Γ h and E[η t+h η t ] = 1 {h=0} ω 2 I, we are left with:

E[(Π † t+h Y t+h )(Π † t Y t ) |Π] = (π t+h π t ) Γ h + 1 {h=0} ω 2 diag(π t ).
where is the elementwise Hadamard product. We now take the expectation w.r.t. Π:

E[(Π † t+h Y t+h )(Π † t Y t ) ] = E[π t+h π t ] Γ h + 1 {h=0} ω 2 E[diag(π t )].
Dividing elementwise by the scaling matrix

S(h) = E[π t+h π t ], we get E 1 E[π t+h π t ] (Π † t+h Y t+h )(Π † t Y t ) = Γ h + 1 {h=0} ω 2 E[diag(π t )] 1 E[π t π t ] = Γ h + 1 {h=0} ω 2 diag E[π t ] E[π 2 t ] = Γ h + 1 {h=0} ω 2 I
which shows that our estimator

Γ h = 1 T -h T -h t=1 1 E[π t+h π t ] (Π † t+h Y t+h )(Π † t Y t ) -1 {h=0} ω 2 I is unbiased.
Note that since the process (Π 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 π are given by

S(h) d1,d2 = E [π t+h,d1 π t,d2 ] =      p 2 if d 1 = d 2 p if d 1 = d 2 and h = 0 p 2 + p(1 -p)(1 -a -b) h if d 1 = d 2 and h ≥ 1
In particular, every coefficient of the scaling matrix S(h) is lower-bounded by min d1,d2,h S(h) d1,d2 = min{p 2 , p(1 -b)} = pq u where q u = min{p, 1 -b}.

Proof. Let i = (t+h, d 1 ) and j = (t, d 2 ) be two indices in [T ]×[D]. We have E[π i ] = E[π 2 i ] = p. If d 1 = d 2 ,
then the variables π i and π j belong to independent Markov chains, and thus E[π i π 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

E[π i π j ] = P(π t,d = 1) × P(π t+h,d = 1|π t,d = 1) = p(T h ) 11 .
When diagonalizing the transition matrix T , we see that the bottom-right coefficient of T h is

(T h ) 11 = a + b(1 -a -b) h a + b = p + (1 -p)(1 -a -b) h .
Plugging this in, we get

E[π i π j ] = p 2 + p(1 -p)(1 -a -b) h .
Among all the possible values of S(h) d1,d2 , the smallest one is

p 2 if 1 -a -b ≥ 0, and p 2 + p(1 -p)(1 -a -b) otherwise. But since p + (1 -p)(1 -a -b) = a a + b + b a + b (1 -a -b) = a + b -ab -b 2 a + b = a(1 -b) + b(1 -b) a + b = 1 -b, we conclude min d1,d2,h S(h) d1,d2 = 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 Γ h , coefficient by coefficient. Let us fix two indices d 1 and d 2 : our goal is to control the deviation of ( Γ h ) d1,d2 around its mean.

Lemma 5 (Deviation of ( Γ h ) d1,d2 ). The deviation probability for ( Γ h ) d1,d2 can be decomposed as follows:

P(|( Γ h -Γ h ) d1,d2 | ≥ u) ≤ P (|g ε Ψ ε LΨ ε g ε -E [g ε Ψ ε LΨ ε g ε ] | ≥ u/4) + P |g η Ψ η LΨ ε g ε -E g η Ψ η LΨ ε g ε | ≥ u/4 + P (|g ε Ψ ε LΨ η g η -E [g ε Ψ ε LΨ η g η ] | ≥ u/4) + P |g η Ψ η LΨ η g η -E g η Ψ η LΨ η g η | ≥ u/4
where the random matrix L is defined in Equation ( 16), Ψ ε and Ψ η are defined in Equation (18), and g ε and g η 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),

( Γ h + 1 {h=0} ω 2 I) d1,d2 = 1 T -h T -h t=1 1 S(h) (Π † t+h Y t+h )(Π † t Y t ) d1,d2 = 1 T -h T -h t=1 1 S(h) d1,d2 e d1 (Π † t+h Y t+h )(Π † t Y t ) e d2 = 1 T -h T -h t=1 Tr e d2 e d1 S(h) d1,d2 (Π † t+h Y t+h )(Π † t Y t )
Equation ( 15) allows us to rewrite

(Π † t+h Y t+h )(Π † t Y t ) : ( Γ h + 1 {h=0} ω 2 I) d1,d2 = 1 T -h T -h t=1 X t diag(π t ) e d2 e d1 S(h) d1,d2 diag(π t+h )X t+h + 1 T -h T -h t=1 η t diag(π t ) e d2 e d1 S(h) d1,d2 diag(π t+h )X t+h + 1 T -h T -h t=1 X t diag(π t ) e d2 e d1 S(h) d1,d2 diag(π t+h )η t+h + 1 T -h T -h t=1 η t diag(π t ) e d2 e d1 S(h) d1,d2 diag(π t+h )η t+h
Let us denote by P t the projection of R T D keeping only the components associated with time t, i.e. such that X t = P t X and η t = P t η. We recognize the following matrix L in all four lines of the expression above:

L = 1 T -h T -h t=1 P t diag(π t )
e d2 e d1 S(h)

This leads to:

( Γ h + 1 {h=0} ω 2 I) d1,d2 = X LX + η LX + X Lη + η Lη
Since X and η both follow centered multivariate Gaussian distributions, we can express them as linear combinations of standard Gaussian vectors g ε and g η of dimension T D (indexed by the source of randomness):

X = Ψ ε g ε and η = Ψ η g η (17)
where Ψ ε and Ψ η are the square roots of the respective covariance matrices

Ψ ε = Cov[X] 1/2 and Ψ η = Cov[η] 1/2 = ωI. ( 18
)
We substitute X and η to get:

( Γ h + 1 {h=0} ω 2 I) d2,d1 = g ε Ψ ε LΨ ε g ε + g η Ψ η LΨ ε g ε + g ε Ψ ε LΨ η g η + g η Ψ η LΨ η g η ,
which implies

( Γ h -Γ h ) d1,d2 = g ε Ψ ε LΨ ε g ε -E[g ε Ψ ε LΨ ε g ε ] + g η Ψ η LΨ ε g ε -E[g η Ψ η LΨ ε g ε ] + g ε Ψ ε LΨ η g η -E[g ε Ψ ε LΨ η g η ] + g η Ψ η LΨ η g η -E[g η Ψ η LΨ η g η ].
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 π, it is random, and so are the products Ψ a LΨ b (where a, b ∈ {ε, η}). 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 Ψ a LΨ b 2 • The Frobenius norm Ψ a LΨ b 2 F • The shifted trace Tr(Ψ a LΨ b -E[Ψ a LΨ 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 = max t∈[T -h] π t+h,d1 π t,d2 (T -h)S(h) d1,d2 and
L 2 F = 1 (T -h) 2 S(h) d1,d2 T -h t=1 π t+h,d1 π t,d2 .
Proof. We first notice that L has a block-superdiagonal structure of rank h:

L = 1 T -h T -h t=1 P t L [t,t+h] P t+h with L [t,t+h] = π t+h,d1 π t,d2 S(h) d1,d2 e d2 e d1 . (19)
The spectral and Frobenius norms of such a matrix can easily be deduced from those of its blocks. Since e d2 e d1 2 = e d2 e d1 F = 1 and the π t are binary-valued, this leads to the following formulas:

L 2 = 1 T -h max t∈[T -h] L [t,t+h] 2 = 1 (T -h)S(h) d1,d2 max t∈[T -h] π t+h,d1 π t,d2 L 2 F = 1 (T -h) 2 T -h t=1 L [t,t+h] 2 F = 1 (T -h) 2 S(h) 2 d1,d2 T -h t=1 π t+h,d1 π t,d2 .
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 ≤ c T pq u
Proof. Note that S(h) d1,d2 ≥ pq u , and since h ∈ {0, 1}, we can state that T -h ≥ cT . By Lemma 6, we deduce

L 2 = max t∈[T -h] π t+h,d1 π t,d2 (T -h)S(h) d1,d2 ≤ 1 (T -h)S(h) d1,d2 ≤ 1 (T -h)pq u ≤ 1 cT pq u .
The Frobenius norm requires a little more work because of the sum it contains.

Lemma 8 (Concentration of the sampling Bernoullis). For all u ∈ [0, 1],

P 1 T -h T -h t=1 π t+h,d1 π t,d2 -S(h) d1,d2 ≥ uS(h) d1,d2 ≤ c 1 exp(-c 2 u 2 T S(h) d1,d2 ).
Proof. We distinguish three cases:

• When d 1 = d 2 and h = 0, we have π t+h,d1 = π t,d2 , which is a 2-state Markov chain with transition matrix T ⊗ I, depicted on Figure 2a.

• When d 1 = d 2 , the couple (π t,d2 , π t+h,d1 ) is a 4-state Markov chain with transition matrix T ⊗ T since the chains π t+h,d1 and π t,d2 evolve along independent dimensions. It is shown on Figure 2b.

• When d 1 = d 2 and h ≥ 1, we must study the (h + 1)-tuple (π t,d1 , π t+1,d1 , ..., π t+h,d1 ). 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.

(0, 0)

(1, 1) In all of these cases, our variable of interest π t+h,d1 π t,d2 is a function of the underlying Markov chain. The relevant functions are:

(a) d1 = d2 & h = 0: transition T (0, 0) (0, 1) (1, 0) (1, 1) (b) d1 = d2: transition T ⊗ T (0, 0) (0, 1) (1, 0) (1, 1) (c) d1 = d2 & h = 1: transition S(1).
f 1 : x → x f 2 : (x, y) → yx f 3 : (x 0 , ..., x h ) → x h x 0 .
We note that since χ ≤ a, b ≤ 1 -χ, all the coefficients of T are greater than χ. Furthermore, all the coefficients of T ⊗ T are greater than χ 2 . Finally, all the coefficients of S(h) h+1 are greater than χ h+1 , because all pairs of states are connected after h + 1 steps. Let us illustrate this with h = 1:

S(1) =     1 -a a 0 0 0 0 b 1 -b 1 -a a 0 0 0 0 b 1 -b     S(1) 2 =     (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     .
Subsequently, all the transition matrices R we are interested in, namely R ∈ {T , T ⊗T , S(h) h+1 }, satisfy the Doeblin condition with r = h + 1 and δ = χ h+1 :

R h+1 ≥ χ h+1    1 • • • 1 . . . . . . . . . 1 • • • 1    .
Since we only consider h ∈ {0, 1} and since χ is fixed for our purposes, these chains fulfill the assumptions of Lemma 34. We thus conclude:

P 1 T -h T -h t=1 π t+h,d1 π t,d2 -S(h) d1,d2 ≥ uS(h) d1,d2 ≤ c 1 exp -c 2 u 2 (T -h)S(h) d1,d2 .
We finally replace T -h with cT 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 9 (Frobenious norm bound for L). For any δ such that Equation (20) holds, with probability at least 1 -δ, the Frobenius norm L 2 F satisfies L 2 F ≤ c T pq u .

Proof. By Lemma 8: for all u ∈ [0, 1],

P 1 T -h T -h t=1 π t+h,d1 π t,d2 ≥ (1 + u)S(h) d1,d2 ≤ c 1 exp(-c 2 u 2 T S(h) d1,d2 ).
We remember the expression of Lemma 6 for L 2 F and notice that:

P L 2 F ≥ 1 + u (T -h)S(h) d1,d2 = P 1 (T -h)S(h) d1,d2 1 T -h T -h t=1 π t+h,d1 π t,d2 S(h) d1,d2 ≥ 1 (T -h)S(h) d1,d2
(1 + u)

≤ c 1 exp(-c 2 u 2 T S(h) d1,d2 )
We finally recall that S(h) d1,d2 ≥ pq u and T -h ≥ cT , so that

P L 2 F ≥ 1 + u cT pq u ≤ P L 2 F ≥ 1 + u (T -h)pq u ≤ P L 2 F ≥ 1 + u (T -h)S(h) d1,d2 ≤ c 1 exp -c 2 u 2 T S(h) d1,d2 ≤ c 1 exp -c 2 u 2 T pq u .
All we need to make sure that P L 2 F ≥ 1+u cT pqu ≤ δ is to choose u such that

c 1 exp -c 2 u 2 T pq u ≤ δ ⇐⇒ u ≥ log(c 1 /δ) c 2 T pq u
Note that we can replace log(c 1 /δ) by a constant times log(1/δ) to simplify expressions: this is possible as long as δ 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 ≤ 1 is implied by

log(1/δ) T pq u ≤ c. ( 20
)
If this holds, then we have

P L 2 F ≥ 2 cT pq u ≤ P L 2 F ≥ 1 + u cT pq u ≤ δ.
This yields the result we wanted.

We now move on to studying the shifted trace of Ψ a LΨ b , which is the last ingredient we need for our application of Lemma 37.

Lemma 10 (Trace bound for the L matrices). For all u ∈ [0, 1],

P(| Tr(Ψ ε LΨ ε -E[Ψ ε LΨ ε ])| ≥ u) ≤ c 1 exp - c 2 u 2 T pq u Γ h 2 2 P(| Tr(Ψ η LΨ η -E[Ψ η LΨ η ])| ≥ u) ≤ c 1 exp - c 2 u 2 T pq u ω 4 .
Proof. We can compute an explicit formula thanks to Equation ( 19): if a ∈ {ε, η} then

Tr(Ψ a LΨ a ) = Tr 1 T -h T -h t=1 Ψ a P t π t+h,d1 π t,d2 S(h) d1,d2 e d2 e d1 P t+h Ψ a = 1 T -h T -h t=1 π t+h,d1 π t,d2 S(h) d1,d2
Tr Ψ a P t e d2 e d1 P t+h Ψ a

= 1 T -h T -h t=1 π t+h,d1 π t,d2 S(h) d1,d2 e d1 P t+h Ψ a Ψ a P t e d2 = 1 T -h T -h t=1 π t+h,d1 π t,d2 S(h) d1,d2 (Ψ a Ψ a ) [t+h,t] d1,d2
where (Ψ a Ψ a ) [t+h,t] d1,d2 denotes the (d 1 , d 2 ) coefficient of the (t + h, t) block of Ψ a Ψ a . Now is the time to look back on Equation ( 18), which tells us that both Ψ ε Ψ ε and Ψ η Ψ η are constant along their superdiagonal of rank h. We thus find that

Tr(Ψ ε LΨ ε -E[Ψ ε LΨ ε ]) = (Γ h ) d1,d2 1 T -h T -h t=1 π t+h,d1 π t,d2 S(h) d1,d2 -1 Tr(Ψ η LΨ η -E[Ψ η LΨ η ]) = (1 {h=0} ω 2 I) d1,d2 1 T -h T -h t=1 π t+h,d1 π t,d2 S(h) d1,d2 -1
Like before, we can apply Lemma 8: for all u ∈ [0, 1],

P 1 T -h T -h t=1 π t+h,d1 π t,d2 S(h) d1,d2 -1 ≥ u ≤ c 1 exp(-c 2 u 2 T S(h) d1,d2 ) ≤ c 1 exp(-c 2 u 2 T pq u ). Since |(Γ h ) d1,d2 | ≤ Γ h 2 and (1 {h=0} ω 2 I) d1,d2
≤ ω 2 , we can deduce

P (| Tr(Ψ ε LΨ ε -E[Ψ ε LΨ ε ])| ≥ u Γ h 2 ) ≤ c 1 exp(-c 2 u 2 T pq u ) P | Tr(Ψ η LΨ η -E[Ψ η LΨ η ])| ≥ uω 2 ≤ c 1 exp(-c 2 u 2 T pq 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 δ > 0 and u ∈ [0, 1]. Assume that Equations (20) and (21) hold. Then the deviation probability for ( Γ h ) d1,d2 satisfies

P(|( Γ h -Γ h ) d1,d2 | ≥ u) ≤ 4δ + c 1 exp - c 2 u 2 T pq u max {( Ψ ε 2 2 + ω 2 ) 2 , Γ h 2 2 , ω 4 } .
Proof. The conclusion we had reached before our discrete interlude is given by Lemma 5, and we can rewrite it as

P(|( Γ h -Γ h ) d1,d2 | ≥ u) ≤ p εε + p ηε + p εη + p ηη ,
where each p ab represents a deviation probability for a specific quadratic form g a Ψ a LΨ b g b . Let us choose δ such that Equation (20) holds. By Lemmas 7 and 9, with probability at least 1 -δ, the following eight inequalities occur at the same time (we use Lemma 25 to split the products):

Ψ ε LΨ ε 2 F ≤ c Ψ ε 4 2 T pq u Ψ ε LΨ ε 2 ≤ c Ψ ε 2 2 T pq u Ψ η LΨ ε 2 F ≤ c Ψ η 2 2 Ψ ε 2 2 T pq u Ψ η LΨ ε 2 ≤ c Ψ η 2 Ψ ε 2 T pq u Ψ ε LΨ η 2 F ≤ c Ψ ε 2 2 Ψ η 2 2 T pq u Ψ ε LΨ η 2 ≤ c Ψ ε 2 Ψ η 2 T pq u Ψ η LΨ η 2 F ≤ c Ψ η 4 2 T pq u Ψ η LΨ η 2 ≤ c Ψ η 2 2
T pq u .

The spectral norm of Ψ η is easily seen to equal

Ψ η 2 = ω 2 I 1/2 2
= ω, which allows us to lighten these expressions. From there, Lemma 37 (applied with X = g a , Y = g b and A = Ψ a LΨ b ) provides the concentration bounds we need2 :

p εε ≤ δ + 2 exp -cT pq u min (u/4) 2 Ψ ε 4 2 , (u/4) Ψ ε 2 2 + P | Tr(Ψ η LΨ ε ) -E[Ψ ε LΨ ε ])| ≥ u/8 p ηε ≤ δ + 2 exp -cT pq u min (u/4) 2 ω 2 Ψ ε 2 2 , (u/4) ω Ψ ε 2 p εη ≤ δ + 2 exp -cT pq u min (u/4) 2 Ψ ε 2 2 ω 2 , (u/4) Ψ ε 2 ω p ηη ≤ δ + 2 exp -cT pq u min (u/4) 2 ω 4 , (u/4) ω 2 + P | Tr(Ψ η LΨ η ) -E[Ψ η LΨ η ])| ≥ u/8 .
The denominators inside the minima can be unified: for the left column,

max Ψ ε 4 2 , Ψ ε 2 2 ω 2 , ω 4 ≤ Ψ ε 2 2 + ω 2 2 ,
and for the right column,

max Ψ ε 2 2 , Ψ ε 2 ω, ω 2 ≤ ( Ψ ε 2 + ω) 2 ≤ 2 Ψ ε 2 2 + ω 2 .
This means we can upper bound each of the four minima by

min u/4 Ψ ε 2 2 + ω 2 2 , u/8 Ψ ε 2 2 + ω 2 .
From now on, we additionally suppose that u/4

Ψ ε 2 2 + ω 2 ≤ 1 2 (21)
This enables us to get rid of these minima by reducing them to the (smaller) quadratic term on the left. We end up with

p εε ≤ δ + 2 exp - cu 2 T pq u ( Ψ ε 2 2 + ω 2 ) 2 + P (| Tr(Ψ ε LΨ ε ) -E[Ψ ε LΨ ε ])| ≥ u/8) p ηε ≤ δ + 2 exp - cu 2 T pq u ( Ψ ε 2 2 + ω 2 ) 2 p εη ≤ δ + 2 exp - cu 2 T pq u ( Ψ ε 2 2 + ω 2 ) 2 p ηη ≤ δ + 2 exp - cu 2 T pq u ( Ψ ε 2 2 + ω 2 ) 2 + P | Tr(Ψ η LΨ η ) -E[Ψ η LΨ η ])| ≥ u/8 .
As for the trace terms, they are taken care of by Lemma 10:

P | Tr(Ψ η LΨ ε ) -E[Ψ ε LΨ ε ])| ≥ u/8 ≤ c 3 exp -c 4 (u/8) 2 T pq u Γ h 2 2 P | Tr(Ψ η LΨ η ) -E[Ψ η LΨ η ])| ≥ u/8 ≤ c 3 exp -c 4 (u/8) 2 T pq u ω 4
We plug this in and rearrange to get:

p εε + p ηε + p εη + p ηη ≤ 4δ + c 1 exp - c 2 u 2 T pq u max {( Ψ ε 2 2 + ω 2 ) 2 , Γ h 2 2 , ω 4 } .
The following result will simplify the denominator inside the exponential.

Lemma 12 (Spectral norms of Ψ ε and Γ h ). The matrices Ψ ε and Γ h satisfy:

Ψ ε 2 2 ≤ σ 2 max (1 -ϑ) 2 and Γ h 2 ≤ ϑ h σ 2 max 1 -ϑ .
As a consequence,

max ( Ψ ε 2 2 + ω 2 ) 2 , Γ h 2 2 , ω 4 ≤ (σ 2 max + ω 2 ) 2 (1 -ϑ) 4
Proof. By Lemma 1, we can write Ψ 2 ε as a sum of Kronecker products (one for each block). Let J t he a matrix full of zeros, except for the subdiagonal of rank t, which is full of ones. Then we have:

Ψ 2 ε = Cov[X] = I ⊗ Γ 0 (θ) + T -1 t=1 J t ⊗ θ t Γ 0 (θ) + J t ⊗ Γ 0 (θ)θ t
This gives us control over its spectral norm thanks to Lemma 24:

Ψ ε 2 2 = Ψ 2 ε 2 ≤ I 2 × Γ 0 (θ) 2 + T -1 t=1 J t 2 × θ t Γ 0 (θ) 2 + J t 2 × Γ 0 (θ)θ t 2 ≤ Γ 0 (θ) 2 1 + 2 T -1 t=1 θ t 2 22
We now use Lemma 2:

Ψ ε 2 2 ≤ σ 2 max 1 -ϑ 2 1 + 2 ϑ 1 -ϑ = σ 2 max (1 -ϑ) 2 .
We now turn to Γ h with Lemmas 1 and 2:

Γ h 2 = θ h Γ 0 (θ) 2 ≤ ϑ h σ 2 max 1 -ϑ 2 .
In particular, we have

max ( Ψ ε 2 2 + ω 2 ) 2 , Γ h 2 2 , ω 4 ≤ max σ 2 max (1 -ϑ) 2 + ω 2 2 , ϑ h σ 2 max 1 -ϑ 2 2 , ω 4 ≤ (σ 2 max + ω 2 ) 2 (1 -ϑ) 4
We can now control the error of the covariance estimator:

Lemma 13 (Max norm convergence rate of the covariance estimator). Let δ > 0 be small enough. Assume that Equations ( 20) and ( 23) hold. Then the covariance estimator Γ h from Equation (9) satisfies

Γ h -Γ h max ≤ c σ 2 max + ω 2 (1 -ϑ) 2 log(D/δ) √ T pq u = err(δ)
with probability greater than 1 -δ.

Proof. Let us plug Lemma 12 into Lemma 11

P(|( Γ h -Γ h ) d1,d2 | ≥ u) ≤ 4δ + c 1 exp - c 2 (1 -ϑ) 4 u 2 T pq u (σ 2 max + ω 2 ) 2 .
All that is left to do is choose u such that

P(|( Γ h -Γ h ) d1,d2 | ≥ u) ≤ 8δ,
which will be true if

c 1 exp - c 2 (1 -ϑ) 4 T pq u (σ 2 max + ω 2 ) 2 u 2 ≤ 4δ ⇐⇒ u ≥ log(c 1 /4δ)(σ 2 max + ω 2 ) 2 c 2 (1 -ϑ) 4 T pq u .
As long as δ is small enough, we can take

u = c log(1/δ)(σ 2 max + ω 2 ) (1 -ϑ) 2 √ T pq u . ( 22
)
For Lemma 11 to apply, we must verify that u ∈ [0, 1] and that Equation ( 21) is satisfied. In other words, we have to ensure that

c log(1/δ)(σ 2 max + ω 2 ) (1 -ϑ) 2 √ T pq u ≤ min{1, 2( Ψ ε 2 2 + ω 2 )}
Using Lemma 12, this is implied by the condition

log(1/δ) max{1, (σ 2 max + ω 2 ) -1 } (1 -ϑ) 2 √ T pq u ≤ c (23) This paper Han et al. [2015] VAR def X t = θX t-1 + ε t X t = A 1 X t-1 + Z t Covariance Γ h = Cov(X h , X 0 ) Σ i = Cov(X 0 , X i ) Yule-Walker Γ h = θ h Γ 0 Σ i = Σ 0 A i 1 Covariance estimate Γ h S i Covariance error err(δ) ζ i Optimization constraint M Γ 0 -Γ 1 max ≤ λ S 0 M -S 1 max ≤ λ Optimization objective vec(M ) 1 vec(M ) 1 Threshold in proof ν λ 1
Table 1: Notation correspondence between this paper and [START_REF] Han | A Direct Estimation of High Dimensional Stationary Vector Autoregressions[END_REF] Under these hypotheses, we just proved that with probability at least 1 -8δ,

|( Γ h -Γ h ) d1,d2 | ≤ c σ 2 max + ω 2 (1 -ϑ) 2 log(1/δ) √ T pq u .
We finish with a union bound, applying the previous result to all pairs (d 1 , d 2 ) ∈ [D] 2 . With probability greater than 1 -8D 2 δ, we have:

max d1,d2 |( Γ h -Γ h ) d1,d2 | = Γ h -Γ h max ≤ c σ 2 max + ω 2 (1 -ϑ) 2 log(1/δ) √ T pq u .
Replacing δ with 8D 2 δ gives us the result we wanted: with probability greater than 1 -δ,

Γ h -Γ h max ≤ c σ 2 max + ω 2 (1 -ϑ) 2 log(D/δ) √ T pq u .

A.7 BEHAVIOR OF THE DANTZIG SELECTOR

We now walk the final steps from the error on Γ h to the error on θ. 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 Γ 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 θ -θ 1 , where θ = Γ 1 Γ -1 0 is the true transition matrix. We start by choosing a specific λ such that θ is feasible with high probability.

Lemma 14 (Feasibility of the real θ). If we select the penalization level λ = ( θ ∞ + 1) err(δ), then with probability at least 1 -δ, the real θ is a feasible solution to the optimization problem (8).

Proof.

θ Γ 0 -Γ 1 max = Γ 1 Γ -1 0 Γ 0 -Γ 1 max = Γ 1 Γ -1 0 Γ 0 -Γ 1 + Γ 1 -Γ 1 max ≤ Γ 1 Γ -1 0 Γ 0 -Γ 1 Γ -1 0 Γ 0 max + Γ 1 -Γ 1 max = θ( Γ 0 -Γ 0 ) max + Γ 1 -Γ 1 max By Lemma 26, θ( Γ 0 -Γ 0 ) max ≤ θ ∞ Γ 0 -Γ 0 max
By Lemma 13, with probability greater than 1 -2δ, Γ 0 -Γ 0 max ≤ err(δ) and Γ 1 -Γ 1 max ≤ err(δ)

which implies θ Γ 0 -Γ 1 max ≤ ( θ ∞ + 1) err(δ).
This is exactly the feasibility criterion for (8) if λ = ( θ ∞ + 1) err(δ).

Lemma 15 (Error on θ in max norm). If we select λ = ( θ ∞ + 1) err(δ), then with probability at least 1 -δ, the max norm error of θ satisfies θ -

θ max ≤ 2λ Γ -1 0 1 .
Proof.

θ -θ max = θ -Γ 1 Γ -1 0 max = ( θΓ 0 -Γ 1 )Γ -1 0 max = ( θΓ 0 -θ Γ 0 + θ Γ 0 -Γ 1 + Γ 1 -Γ 1 )Γ -1 0 max ≤ ( θΓ 0 -θ Γ 0 )Γ -1 0 max + ( θ Γ 0 -Γ 1 )Γ -1 0 max + ( Γ 1 -Γ 1 )Γ -1 0 max By Lemma 26, θ -θ max ≤ θ(Γ 0 -Γ 0 ) max + θ Γ 0 -Γ 1 max + Γ 1 -Γ 1 max Γ -1 0 1 ≤ θ ∞ Γ 0 -Γ 0 max + θ Γ 0 -Γ 1 max + Γ 1 -Γ 1 max Γ -1 0 1
We want to control θ ∞ using θ ∞ . Let us recall that the operator ∞ norm is equal to the maximum 1 norm of the rows of a matrix. To control the rows of θ, we notice that the optimization problem defining θ, namely min

M ∈R D×D vec(M ) 1 s.t. M Γ 0 -Γ 1 max ≤ λ is equivalent to the row-wise minimization ∀i, min Mi,•∈R 1×D M i,• 1 s.t. M i,• Γ 0 -( Γ 1 ) i,• max ≤ λ
From this, we deduce that each row of the optimum θ satisfies

θ i,• 1 ≤ θ i,• 1 , which implies θ ∞ ≤ θ ∞ .
Going back to our error estimate, we get:

θ -θ max ≤ θ ∞ Γ 0 -Γ 0 max + θ Γ 0 -Γ 1 max + Γ 1 -Γ 1 max Γ -1 0 1
Note that the middle term is smaller than λ because the optimum θ is a feasible solution. Meanwhile, the first and third term are smaller than err(δ) with probability 1 -δ:

θ -θ max ≤ ( θ ∞ err(δ) + λ + err(δ)) Γ -1 0 1 = 2λ Γ -1 0 1
To complete the proof of Theorem 1, we simply need to go from the max norm to the ∞ operator norm.

Proof. Let ν > 0 be a threshold (to be chosen later). We define

s 1 = max i j min |θ i,j | ν , 1 and I i = {j : |θ i,j | ≥ ν}
With high probability, the following holds for any row i:

θ i,• -θ i,• 1 ≤ θ i,I c i -θ i,I c i 1 + θ i,Ii -θ i,Ii 1 ≤ θ i,I c i 1 + θ i,I c i 1 + θ i,Ii -θ i,Ii 1 = ( θ i,• 1 -θ i,Ii 1 ) + θ i,I c i 1 + θ i,Ii -θ i,Ii 1 ≤ θ i,• 1 -θ i,Ii 1 + θ i,I c i 1 + θ i,Ii -θ i,Ii 1 = ( θ i,Ii 1 + θ i,I c i 1 ) -θ i,Ii 1 + θ i,I c i 1 + θ i,Ii -θ i,Ii 1 = 2 θ i,I c i 1 + ( θ i,Ii 1 -θ i,Ii 1 ) + θ i,Ii -θ i,Ii 1 ≤ 2 θ i,I c i 1 + 2 θ i,Ii -θ i,Ii 1 By definition of I i , for all j ∈ I c i , |θ i,j | ≤ ν, hence θ i,I c i 1 = j∈I c i |θ i,j | = j∈I c i min{|θ i,j |, ν} ≤ j min{|θ i,j |, ν} ≤ νs 1
Meanwhile, the second term satisfies

θ i,Ii -θ i,Ii 1 ≤ |I i | × θ -θ max
And by definition of I i , for all j ∈ I i , |θ i,j | ≥ ν, hence

|I i | = j∈Ii 1 = j∈Ii min |θ i,j | ν , 1 ≤ j min |θ i,j | ν , 1 ≤ s 1
Combining all of this, we get that with high probability,

θ i,• -θ i,• 1 ≤ 2(ν + 2λ Γ -1 0 1 )s 1
Judging by the last Equation, it makes sense to choose ν = 2λ Γ -1 0 1 . Furthermore, our sparsity hypothesis on θ implies that for all but s of the coefficients of any row i, min{|θ i,j |, ν} = |θ i,j | = 0. We deduce that for every i,

j min {|θ i,j |, ν} ≤ s max j min {|θ i,j |, ν} ≤ νs which directly implies νs 1 = max i j min {|θ i,j |, ν} ≤ νs
We finally find that with high probability,

θ i,• -θ i,• 1 ≤ 4νs 1 ≤ 4νs = 8λ Γ -1 0 1 s
With the help of a union bound, again with high probability,

θ -θ ∞ = max i θ i,• -θ i,• 1 ≤ 8λ Γ -1 0 1 s
We substitute the value of λ and obtain

θ -θ ∞ ≤ 8( θ ∞ + 1) err(δ) Γ -1 0 1 s
Once we plug in the value of err(δ), the resulting high-probability error bound reads

θ -θ ∞ ≤ c θ ∞ + 1 Γ -1 0 -1 1 σ 2 max + ω 2 (1 -ϑ) 2 s log(D/δ) √ T pq u
Since ϑ only acted as an upper bound on θ 2 in this proof, we can define

γ u (θ) = θ ∞ + 1 (1 -θ 2 ) 2 σ 2 max + ω 2 Γ -1 0 -1 1 to obtain the compressed expression θ -θ ∞ ≤ cγ u (θ) s log(D/δ) √ T pq u .

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 θ 0 , θ 1 , ..., θ M satisfying two seemingly contradictory conditions: their induced distributions must be hard to distinguish, yet they must lie as fart 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 θ i and θ 0 . Taking the latter to be 0, we actually want to bound

1 M + 1 M i=1 KL {P θi (Π, Y ) P 0 (Π, Y )} ≤ max i KL {P θi (Π, Y ) P 0 (Π, Y )} By Lemma 28, KL {P θi (Π, Y ) P 0 (Π, Y )} = KL {P θi (Π) P 0 (Π)} + E Π [KL {P θi (Y |Π) P 0 (Y |Π)}]
Since θ i does not affect the distribution of the sampling process Π, 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 Π, 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 Cov θ [Y |Π] and decompose it into a constant term Q Π (corresponding to the independent case θ = 0) plus a residual R Π (θ) (Lemma 16).

2. Upper-bound the conditional KL divergence KL {P θ (Y |Π) P 0 (Y |Π)} using the "deviations from the identity"

∆ Π (θ) = Q -1/2 Π R Π (θ)Q -1/2 Π (Lemma 17).
3. Control ∆ Π (θ) using features of R(θ) scaled by sampling-related factors (Lemmas 18,19 and 20). 5. Apply Fano's method to a set of parameters θ 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 Π t that contain only zeros are removed. In other words, Π t is no longer the diagonal matrix diag(π t ) but instead becomes a wide rectangular matrix with exactly one 1 per row and at most one 1 per column. We thus have ΠΠ = I unless all of the π t,d are zero, in which case the matrix Π is empty, and so are the observations Y . Let us denote this very unlikely event by E, and its complement by E c . If Π is such that E happens, we obviously have

KL {P θi (Y |Π) P 0 (Y |Π)} = 0, which means that E Π [KL {P θi (Y |Π) P 0 (Y |Π)}] = E Π [1 E c KL {P θi (Y |Π) P 0 (Y |Π)}] (24)
For the beginning of the proof, we consider a fixed, non-empty realization of Π.

B.3 COVARIANCE DECOMPOSITION

As we announced in the proof sketch, our reference parameter will be θ 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 Π decomposes as

Cov θ [Y |Π] = Q Π + R Π (θ),
where Q Π is a constant term and R Π (θ) is a residual which vanishes as θ → 0. They are defined as follows: the constant term is

Q Π = Π(bdiag T Σ)Π + ω 2 I whereas the residual equals R Π (θ) = ΠR(θ)Π with R(θ) =      θΓ 0 (θ)θ Γ 0 (θ)θ 1 Γ 0 (θ)θ 2 • • • θ 1 Γ 0 (θ) θΓ 0 (θ)θ Γ 0 (θ)θ 1 θ 2 Γ 0 (θ) θ 1 Γ 0 (θ) θΓ 0 (θ)θ . . . . . .      .
Proof. We use Equation ( 5) to see that the conditional distribution

P θ (Y |Π) is a centered multivariate Gaussian with covariance Cov θ [Y |Π] = ω 2 I + Π Cov θ [X]Π .
We then use Lemma 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 Σ:

Cov θ [Y |Π] = ω 2 I + Π bdiag T (Σ)Π + Π (Cov θ [X] -bdiag T (Σ)) Π .
Finally, we define Lemma 17. Let us define the deviation from the identity:

Q Π = ω 2 I + Π bdiag T (Σ)Π , R(θ) = Cov θ [X] -bdiag T (Σ)
∆ Π (θ) = Q -1/2 Π R Π (θ)Q -1/2 Π .
Then the conditional KL divergence is upper-bounded by:

KL {P θ (Y |Π) P 0 (Y |Π)} ≤ ∆ Π (θ) 2 F 2(1 + λ min (∆ Π (θ))) .
Proof. The conditional KL divergence KL {P θ (Y |Π) P 0 (Y |Π)} 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,

Cov 0 (Y |Π) = Q Π = Q 1/2 Π (Q 1/2 Π ) Cov θ (Y |Π) = Q Π + R Π (θ) = Q 1/2 Π I + Q -1/2 Π R Π (θ)Q -1/2 Π ∆Π(θ) (Q 1/2 Π ) .
By Lemma 23, there exists a real number r min ≥ s min Q

1/2 Π 2 = s min (Q Π ) such that λ min (Cov θ (Y |Π)) = r min λ min (I + ∆ Π (θ)).
Since Q Π ω 2 I 0, its minimum singular value satisfies s min (Q Π ) > 0, so that r min > 0. In addition, Cov θ (Y |Π) ω 2 I 0, so that λ min (Cov θ (Y |Π)) > 0. Therefore,

λ min (I + ∆ Π (θ)) = λ min (Cov θ (Y |Π)) r min > 0 and λ min (∆ Π (θ)) > -1,
which means we can apply Lemma 30 with P 1 = P θ (Y |Π) and P 0 = P 0 (Y |Π).

B.5 FROM ∆ Π (θ) TO R Π (θ)

Lemma 17 strongly suggests studying a certain fraction involving ∆ Π (θ). In the following result, we boil it down to a function of the residual term R Π (θ).

Lemma 18. Assume R(θ) 2 ≤ (σ 2 min + ω 2 )/2. We have the following upper bound:

∆ Π (θ) 2 F 2(1 + λ min (∆ Π (θ))) ≤ R Π (θ) 2 F (σ 2 min + ω 2 ) 2 .
Proof. Since the quantity λ min (∆ Π (θ)) in the denominator is hard to control, we will work with the spectral norm instead. Indeed, whenever ∆ Π (θ) 2 < 1, we have the crude bound

1 1 -λ min (∆ Π (θ)) ≤ 1 1 -∆ Π (θ) 2 .
Let us start by noticing that, thanks to Lemma 25,

∆ Π (θ) 2 F = Q -1/2 Π R Π (θ)Q -1/2 Π 2 F ≤ Q -1/2 Π 4 2 R Π (θ) 2 F = Q -1 Π 2 2 R Π (θ) 2 F ∆ Π (θ) 2 = Q -1/2 Π ΠR(θ)Π Q -1/2 Π 2 ≤ Q -1/2 Π Π 2 2 R(θ) 2 .
We will later see how the spectral and Frobenius norms of the residual R(θ) can be controlled as a function of θ. For now, we must work to upper bound Q -1 Π 2 and Q -1 Π Π 2 2 . To simplify the following proof, we write Σ d = bdiag T Σ. Since Σ d is block-diagonal, its spectrum is the same as the spectrum of Σ repeated T times, hence λ min (Σ d ) = σ 2 min . And since we assumed E c happens (non empty projection), we have ΠΠ = I and Π Π = diag(π), which has at least one entry equal to 1.

We start with Q -1 Π 2 . Since Q Π ω 2 I 0 is non-singular and symmetric,

Q -1 Π 2 = λ max (Q -1 Π ) = 1 λ min (Q Π ) = 1 λ min (ΠΣ d Π + ω 2 I) .
Remembering that Σ d σ 2 min I, we get

ΠΣ d Π + ω 2 I σ 2 min ΠΠ + ω 2 I = (σ 2 min + ω 2 )I
and thus

Q -1 Π 2 ≤ 1 σ 2 min + ω 2 .
We now continue with

Q -1/2 Π Π 2 2
. By definition of the spectral norm,

Q -1/2 Π Π 2 2 = λ max Π Q -1 Π Π = λ max Π (ΠΣ d Π + ω 2 I) -1 Π .
Because matrix inversion is decreasing w.r.t. the Loewner order on positive semi-definite matrices,

(ΠΣ d Π + ω 2 I) -1 (σ 2 min + ω 2 ) -1 I -1 Π (ΠΣ d Π + ω 2 I) -1 Π 1 σ 2 min + ω 2 Π Π. It follows that Q -1/2 Π Π 2 2 ≤ 1 σ 2 min + ω 2 λ max (Π Π) = 1 σ 2 min + ω 2 .
The conclusion is within reach:

∆ Π (θ) 2 F 1 + λ min (∆ Π (θ)) ≤ ∆ Π (θ) 2 F 1 -∆ Π (θ) 2 ≤ Q -1 Π 2 2 R Π (θ) 2 F 1 -Q -1/2 Π Π 2 2 R(θ) 2 ≤ 1 σ 2 min +ω 2 2 R Π (θ) 2 F 1 - 1 σ 2 min +ω 2 R(θ) 2 ≤ 2 R Π (θ) 2 F (σ 2 min + ω 2 ) 2
The last inequality is justified by our assumption R(θ) 2 ≤ (σ 2 min + ω 2 )/2. Another consequence of this assumption is that

∆ Π (θ) 2 ≤ Q -1/2 Π Π 2 2 R(θ) 2 ≤ 1 σ 2 min + ω 2 σ 2 min + ω 2 2 = 1 2 < 1
which is sufficient for the first inequality to hold.

B.6 FROM R Π (θ) TO R(θ)

As the previous Lemma underlines, the last step we need to get rid of the dependency in Π is to study the average norm of R Π (θ).

Lemma 19. Let q = max{1 -b, 2p -(1 -b)}. Then E Π 1 E c R Π (θ) 2 F ≤ p Tr[R(θ) R(θ)] + pq R(θ) 2 F .
Proof. We first notice that for any matrix A,

E Π 1 E c ΠAΠ 2 F = E Π 1 E c Tr ΠAΠ ΠA Π = E Π Tr diag(π)A diag(π)A = i,j E Π [π i π j ]A 2 i,j .
We can apply this to R Π (θ) = ΠR(θ)Π :

E Π 1 E c R Π (θ) 2 F = i,j E Π [π i π j ]R(θ) 2 i,j .
The rest of the proof consists in plugging in the moments E Π [π i π j ] from Lemma 4:

E Π R Π (θ) 2 F = t1,t2,d1,d2 (t1,d1)=(t2,d2) pR(θ) 2 (t1,d1),(t2,d2) + t1,t2,d1,d2 d1 =d2 p 2 R(θ) 2 (t1,d1),(t2,d2) + t1,t2,d1,d2 d1=d2,t1 =t2 (p 2 + p(1 -p)(1 -a -b) |t1-t2| )R(θ) 2 (t1,d1),(t2,d2) .
The sum in the last term can be crudely controlled as follows:

t1,t2,d t1 =t2 (1 -a -b) |t1-t2| R(θ) 2 (t1,d),(t2,d) ≤ |1 -a -b| t1,t2 t1 =t2 d (R(θ) [t1,t2] ) 2 d,d ≤ |1 -a -b| t1 =t2 R(θ) [t1,t2] 2 F ≤ |1 -a -b| • R(θ) 2 F
This yields a short, but probably suboptimal bound:

E Π 1 E c R Π (θ) 2 F ≤ p Tr[R(θ) R(θ)] + (p 2 + p(1 -p)|1 -a -R(θ) 2 F .
In the previous part, we already saw that

p + (1 -p)(1 -a -b) = 1 -b.
Similarly, we obtain

p + (1 -p)(a + b -1) = a a + b + b a + b (a + b -1) = a + ba + b 2 -b a + b = a(1 + b) -b(1 -b) a + b = p(1 + b) -(1 -p)(1 -b) = 2p -(1 -b). As a consequence, p + (1 -p)|1 -a -b| = max{1 -b, 2p -(1 -b)} = q ,
which yields the expected result.

B.7 BOUNDING R(θ)

Lemma 19 relates the bounds involving R Π (θ) to features of the full residual R(θ), which we now study.

Lemma 20. The residual R(θ) satisfies the following inequalities:

R(θ) 2 ≤ 2σ 2 max (1 -ϑ) 2 θ 2 R(θ) 2 F ≤ 2T σ 4 max (1 -ϑ) 3 θ 2 F Tr[R(θ) R(θ)] ≤ T σ 4 max (1 -ϑ) 2 θ 2 2 θ 2 F .
Proof. We start by giving a formula for the blocks of R(θ): by Lemma 16,

R(θ) [t,s] =      θ t-s Γ 0 (θ) if s ∈ [1, t -1] θΓ 0 (θ)θ if s = t Γ 0 (θ)θ t-s if s ∈ [t + 1, T ].
These individual blocks can be bounded using Lemmas 25 and 2: if r ≥ 1, then

θ r Γ 0 (θ) 2 F ≤ Γ 0 (θ) 2 2 θ r 2 F ≤ Γ 0 (θ) 2 2 θ 2 F θ r-1 2 2 ≤ σ 4 max (1 -ϑ) 2 θ 2 F θ 2(r-1) 2 Γ 0 (θ)θ r 2 F ≤ Γ 0 (θ) 2 2 θ r 2 F ≤ Γ 0 (θ) 2 2 θ 2 F θ r-1 2 2 ≤ σ 4 max (1 -ϑ) 2 θ 2 F θ 2(r-1) 2 θΓ 0 (θ)θ 2 F ≤ θ 2 2 Γ 0 (θ) 2 2 θ 2 F ≤ σ 4 max (1 -ϑ) 2 θ 2 F θ 2 2 .
Since we control the norm of each block of R(θ), we control the norm of the whole matrix:

R(θ) 2 F = T t=1 t-1 s=1 θ t-s Γ 0 (θ) 2 F + θΓ 0 (θ)θ 2 F + T s=t+1 Γ 0 (θ)θ s-t 2 F ≤ σ 4 max θ 2 F (1 -ϑ 2 ) 2 T t=1 t-1 s=1 θ 2(t-s-1) 2 + θ 2 2 + T s=t+1 θ 2(s-t-1) 2 ≤ σ 4 max θ 2 F (1 -ϑ 2 ) 2 T t=1 t-1 s=-∞ θ 2(t-1-s) 2 + θ 2 2 + +∞ s=t+1 θ 2(s-1-t) 2 = σ 4 max θ 2 F (1 -ϑ 2 ) 2 T 1 1 -θ 2 2 + θ 2 2 + 1 1 -θ 2 2
We now remember our hypothesis θ 2 ≤ ϑ < 1:

R(θ) 2 F ≤ σ 4 max θ 2 F (1 -ϑ 2 ) 2 T 1 1 -ϑ 2 + ϑ 2 + 1 1 -ϑ 2 = σ 4 max θ 2 F (1 -ϑ 2 ) 2 T 2 + ϑ 2 (1 -ϑ 2 ) 1 -ϑ 2 ≤ σ 4 max θ 2 F (1 -ϑ 2 ) 2 T 2 + 2ϑ 1 -ϑ 2 = σ 4 max θ 2 F (1 -ϑ 2 ) 2 T 2 1 -ϑ = 2T σ 4 max θ 2 F (1 -ϑ) 3 .
Now that we have a handle on the Frobenius norm of R(θ), we move on to its spectral norm. Notice that R(θ) can be written as a sum of Kronecker products with the subdiagonal matrices J t :

R(θ) = I ⊗ θΓ 0 (θ)θ + T -1 t=1 J t ⊗ θ t Γ 0 (θ) + J t ⊗ Γ 0 (θ)θ t .
We can use Lemma 24 and write:

R(θ) 2 ≤ I 2 × θΓ 0 (θ)θ 2 + T -1 t=1 J t 2 × θ t Γ 0 (θ) 2 + J t 2 × Γ 0 (θ)θ t 2 ≤ Γ 0 (θ) 2 θ 2 2 + 2 T -1 t=1 θ t 2 ≤ σ 2 max 1 -ϑ 2 θ 2 2 + 2 θ 2 1 -θ 2 ≤ σ 2 max θ 2 1 -ϑ 2 ϑ + 2 1 -ϑ ≤ σ 2 max θ 2 1 -ϑ 2 + 2ϑ 1 -ϑ 2 = 2 σ 2 max θ 2 (1 -ϑ) 2 .
We finish with the trace of the Hadamard product R(θ) R(θ).

Tr[R(θ) R(θ)] = T Tr[(θΓ 0 (θ)θ ) (θΓ 0 (θ)θ )] ≤ T θΓ 0 (θ)θ 2 F ≤ T σ 4 max θ 2 2 θ 2 F (1 -ϑ) 2 .

B.8 UPPER BOUND ON THE KL DIVERGENCE

We now have all the tools in hand to extract a KL divergence bound.

Lemma 21 (Final KL bound). Assume θ ∈ Θ s satisfies

θ 2 ≤ (1 -ϑ) 2 (σ 2 min + ω 2 ) 4σ 2 max
then the expected conditional KL divergence is upper-bounded as follows:

E Π [KL {P θ (Y |Π) P 0 (Y |Π)}] ≤ KL avg ( θ 2 , θ F )
where we defined

γ = (1 -ϑ) 3/2 σ 2 min + ω 2 σ 2 max and KL avg ( θ 2 , θ F ) = 2T p( θ 2 2 + q ) θ 2 F γ .
Proof. Let us start with Lemma 17 on the conditional KL divergence between P θ (Y |Π) and P 0 (Y |Π): for any nonempty Π,

KL {P θ (Y |Π) P 0 (Y |Π)} ≤ ∆ Π (θ) 2 F 2(1 + λ min (∆ Π (θ)))
We continue with Lemma 18 linking ∆ Π (θ) to R Π (θ). As long as R(θ) 2 ≤ (σ 2 min + ω 2 )/2 (we will see to that at the end), we have

KL {P θ (Y |Π) P 0 (Y |Π)} ≤ R Π (θ) 2 F (σ 2 min + ω 2 ) 2 .
Taking the expectation on the event E c yields:

E Π [1 E c KL {P θ (Y |Π) P 0 (Y |Π)}] ≤ E Π 1 E c R Π (θ) 2 F (σ 2 min + ω 2 ) 2
We can now apply Lemma 19:

E Π [1 E c KL {P θ (Y |Π) P 0 (Y |Π)}] ≤ p Tr[R(θ) R(θ)] + pq R(θ) 2 F (σ 2 min + ω 2 ) 2
We substitute the residual bounds from Lemma 20:

E Π [1 E c KL {P θ (Y |Π) P 0 (Y |Π)}] ≤ p × T σ 4 max (1-ϑ) 2 θ 2 2 θ 2 F + pq × 2T σ 4 max (1-ϑ) 3 θ 2 F (σ 2 min + ω 2 ) 2 ≤ σ 2 max σ 2 min + ω 2 2 2T p( θ 2 2 + q ) θ 2 F (1 -ϑ) 3 = 2T p( θ 2 2 + q ) θ 2 F γ .
By Equation ( 24), this is equivalent to bounding the expected KL divergence regardless of the event E, hence the result. Note that our assumption on θ, combined with Lemma 20, implies

R(θ) 2 ≤ 2σ 2 max (1 -ϑ) 2 θ 2 ≤ 2σ 2 max (1 -ϑ) 2 (1 -ϑ) 2 (σ 2 min + ω 2 ) 4σ 2 max ≤ σ 2 min + ω 2 2
B.9 APPLICATION OF FANO'S METHOD

Given the KL bound we just obtained, we are finally able to prove Theorem 2.

Proof. Fano's method requires finding M + 1 parameters θ i such that θ 0 = 0 and θ i -θ j F ≥ 2τ for i = j (with τ to be specified), while keeping control upon the average KL divergence between the probability distributions P θi and P 0 . Judging by Lemma 21, one way to achieve this control on the KL divergence is to bound the θ i F uniformly in i (in other words, to choose them all inside a ball of fixed radius). We will then have to see how many 2τ -separated matrices we can fit in such a ball. Let us consider the set H(r) of all block-diagonal D × D matrices with coefficients in {0, r} such that each block has size s × s (we assume s divides D). In particular, these matrices are all row-and column-sparse, with no more than s non-zero coefficients per row or column. In terms of dimensionality, we are dealing with the (scaled) matrix equivalent of a Ds-dimensional hypercube, hence the notation H(r). It has cardinality 2 Ds and for every θ ∈ H, we have the following norm bounds: θ 2 ≤ rs and θ F ≤ r √ Ds.

The spectral norm bound on θ is obtained as the maximum spectral norm of each block, which we in turn control using the Frobenius norm of each block.

Unfortunately, in this hypercube, not all pairs of vertices are well-separated. That is why we need the Gilbert-Varshamov bound of Lemma 35: according to this result, there exists a pruned subset K(r) ⊂ H(r) containing 0 and such that |K(r)| ≥ |H(r)| 1/8 = 2 Ds/8 and vec(θ i ) -vec(θ j ) 1 ≥ rDs 8 for all pairs of distinct vertices θ i and θ j in K(r). We choose our set of parameters θ 0 , θ 1 , ..., θ M to be exactly this pruned subset K(r), in particular M + 1 = |K(r)|.

The missing ingredient is an upper bound on the maximum average KL divergence between P θi and P 0 : we can obtain it using Lemma 21. We only need to assume

θ i F ≤ r √ Ds ≤ min ϑ, (1 -ϑ) 2 (σ 2 min + ω 2 ) 4σ 2 max to get the upper bound max i E Π [KL {P θi (Y |Π) P θ0 (Y |Π)}] ≤ max i KL avg ( θ i 2 , θ i F ) ≤ KL avg (rs, r √ Ds).
Since we must satisfy the constraint from Equation ( 28) in Fano's method, we will choose r so that:

KL avg (rs, r √ Ds) ≤ α log(M ) = α log 2 Ds/8 -1 with α = log 3-log 2 2 log 2 . We want to solve the previous inequality for r, and for that we start by replacing KL avg (rs, r √ Ds) with its value from Lemma 21, replacing γ with γ to lighten notations:

KL avg (rs, r √ Ds) ≤ α log 2 Ds/8 -1 ⇐⇒ 2 γ T p (rs) 2 + q (r √ Ds) 2 ≤ cDs ⇐⇒ Ds 3 r 4 + q Dsr 2 -c γ 2 Ds T q ≤ 0.
If we consider this as a degree two polynomial in the variable r 2 , its determinant is

∆ = q 2 D 2 s 2 + 4Ds 3 c γ 2 Ds T p .
For β to be small enough, r 2 must remain below the only positive root of the polynomial, namely

r 2 ≤ -q Ds + q 2 D 2 s 2 + c γ 2 D 2 s 4 T p 2Ds 3 = q 2s 2 1 + c γ 2 s 2 T pq 2 -1 .
If we assume the quantity c

γ 2 s 2 T pq 2 inside the square root is smaller than 1, i.e. γ s √ pq √ T ≤ c, (25)
then we can lower-bound √ 1 + x by its chord ( √ 2 -1)x. In other words, a sufficient condition for r 2 to remain small enough is given by

r 2 ≤ q 2s 2 × ( √ 2 -1)c γ 2 s 2 T pq 2 = c γ 2
T pq .

To sum up, we have three constraints on r:

rs ≤ ϑ rs ≤ (1 -ϑ) 2 (σ 2 min + ω 2 ) 4σ 2 max = √ 1 -ϑ 4 γ r ≤ c γ 2 T pq .
We can therefore choose r as the largest value satisfying all three of them:

r = 1 s min ϑ, γ √ 1 -ϑ 4 , c γ s √ T pq (26)
To reach our conclusion, we simply need to remark that the vectorized 1 distance between any two matrices in K(r)

gives us a lower bound on the operator ∞ distance that separates them:

θ i -θ j ∞ = max k∈[D] l∈[D] |(θ i -θ j )| k,l ≥ 1 D 1≤k,l≤D |(θ i -θ j )| k,l = 1 D vec(θ i ) -vec(θ j ) 1 ≥ rDs 8D = rs 8
Subsequently, our parameters θ i are 2τ -separated (in ∞ operator distance) with τ = rs/8. As soon as the minimum in Equation ( 26) is reached by the third value, i.e. whenever

γ s √ T pq ≤ c min{ϑ, γ √ 1 -ϑ} (27)
Lemma 30 (KL divergence between close Gaussians). Let ∆ be a symmetric matrix of size n such that λ min (∆) > -1, and let M be a rectangular matrix such that M M 0. Then the KL divergence between Proof. From Lemma 29 (beware of the switch between P 0 and P 1 ) we get:

KL {P 1 P 0 } = 1 2 Tr(Σ 1 Σ -1 0 ) + (µ 0 -µ 1 ) Σ -1 0 (µ 0 -µ 1 ) -n + logdet(Σ 0 Σ -1 1 ) = 1 2 Tr(M (I + ∆)M -1 ) -n -logdet(M (I + ∆)M -1 ) = 1 2 (Tr(∆) -logdet(I + ∆)) .

As it happens, for small deviations from the identity, the log-determinant is almost equal to the trace. Indeed, since Lemma 32 (Doeblin condition and mixing time). Let (X t ) be an irreducible aperiodic Markov chain with state space X , transition matrix P and stationary distribution µ. Suppose that (X t ) satisfies the Doeblin condition:

∃r ∈ N, ∃δ > 0, ∀(x, y) ∈ X 2 , P r (x, y) ≥ δµ(y).

Then the mixing time of X t , defined as We only have to compute the Fisher information matrix I T (θ). The covariance matrix of Y is given by Lemma 16, but in our case the sampling matrix is constant, and we obtain the simpler (unconditional) result

Cov θ [Y ] = (σ 2 + ω 2 )I T + R(θ),

where the residual R(θ) is of order 1 in θ. Indeed, our simplifying assumptions imply Γ 0 (θ) = σ 2 1-θ 2 and therefore

R(θ) = σ 2 1 -θ 2      θ 2 θ 1 θ 2 • • • θ 1 θ 2 θ 1 θ 2 θ 1 θ 2 . . . . . .      ∂ θ R(θ) = σ 2      0 1 0 • • • 1 0 1 0 1 0 . . . . . .      + O(θ).
The Fisher information of Y with respect to θ has an explicit formula [Malagò and Pistone, 2015, Section 3.5]:

I T (θ) = 1 2 Tr Cov θ [Y ] -1 ∂ θ Cov θ [Y ] Cov θ [Y ] -1 ∂ θ Cov θ [Y ] = 1 2 Tr I + R(θ) σ 2 + ω 2 -1 ∂ θ R(θ) σ 2 + ω 2 I + R(θ) σ 2 + ω 2 -1 ∂ θ R(θ) σ 2 + ω 2 .
If assume θ is small and perform a Taylor expansion, we get:

I T (θ) ≈ 1 2(σ 2 + ω 2 ) 2 Tr (∂ θ R(θ)) 2 .
Incidentally, we also note that at the lowest order in θ,

Tr[(∂ θ R(θ)) 2 ] = ∂ θ R(θ) 2 F ≈ 2σ 4 (T -1).
Which gives us an approximate information matrix for T steps:

I T (θ) ≈ Tr[(∂ θ R(θ)) 2 ] 2(σ 2 + ω 2 ) 2 ≈ T 2 σ 2 σ 2 + ω 2 2 .
Taking the temporal limit yields:

I ∞ (θ) = lim T →∞ I T (θ) T ≈ 1 2 σ 2 σ 2 + ω 2 2 .
In conclusion, this informal analysis reveals an asymptotic error equivalent to

1 √ T I ∞ (θ) -1 ≈ √ 2 √ T 1 + ω 2 σ 2 ,
which is coherent with the dependency we identified in Theorem 2. Given a real number a, we denote by |a| its absolute value. Given a vector x, we denote by x 2 (resp. x 1 , x ∞ , x 0 ) its Euclidean norm (resp. 1 norm, ∞ norm, number of nonzero entries). The notation e i stands for a vector with a single non-zero coordinate at position i.

A matrix can be defined by its coefficients M = (M i,j ) i,j or by its blocks M = (M [b1,b2] ) b1,b2 . We write I for the identity matrix, and J r for the square matrix entirely filled with zeros, except for the subdiagonal of rank r which is filled with ones. The notation diag(λ) stands for the diagonal matrix with coefficients λ 1 , ..., λ n , while bdiag T (M ) stands for a block-diagonal matrix with T copies of M on the diagonal and zeros elsewhere. We write vec(M ) for the column-wise flattening of matrix M into a vector. When we want to apply a function elementwise, we often use notation that is standard for real numbers but not for matrices: for instance, √ M = ( M i,j ) i,j and 1/M = (1/M i,j ) i,j . Given a real matrix M , we denote by • M its transposition, M † its Moore-Penrose pseudo-inverse and M -1 its inverse;

• Tr(M ) its trace and det(M ) its determinant;

• λ max (M ) (resp. λ min (M ), λ i (M )) its maximum (resp. minimum, i-th largest) eigenvalue, so that

λ max (M ) = λ 1 (M ) ≥ λ 2 (M ) • • • ≥ λ n (M ) = λ min (M )

ppFigure 1 :

 1 Figure 1: Impact of model parameters on the estimation error

Figure 2 :

 2 Figure 2: State space and transitions for the Markov chains used in the discrete concentration result

4.

  Deduce an upper bound on the KL divergence E Π [KL {P θ (Y |Π) P 0 (Y |Π)}] (Lemma 21).

  and R Π (θ) = ΠR(θ)Π to obtain the decomposition we announced. The diagonal blocks of R(θ) are easily computed by noticing that Γ 0 (θ) -Σ = θΓ 0 (θ)θ .B.4 FROM THE KL DIVERGENCE TO ∆ Π (θ)Judging by Lemma 16, 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 |Π).

P 1 = 2 F2

 12 N (µ, M (I + ∆)M ) and P 0 = N (µ, M M ) satisfies KL {P 1 P 0 } ≤ ∆ (1 + λ min (∆)).

  min (∆) .Lemma 31 (Chernoff inequality for Bernoulli variables). Let (X t ) be sequence of independent B(p) variables. Their average satisfies∀u ∈ [0, 1], P 1 T T t=1 X t -p ≥ up ≤ c 1 exp -c 2 u 2 T p .Proof. SeeDubhashi and Panconesi [2009, Theorem 1.1].

t

  mix ( ) = min t ∈ N : max x∈X P t (x, •) -µ TV ≤ , satisfies:t mix ( ) ≥ r 1 + log 1 log 1 1-δ .• The true parameter θ is identifiable and does not lie on the boundary of Θ s . Identifiability is easily deduced from Lemma 1 by observing that θ = Γ 1 (θ)Γ 0 (θ) -1 can be entirely deduced from distribution moments.Since all of these prerequisites hold here,Douc et al. [2014, Proposition 2.14] gives us a Central Limit Theorem for the MLE of linear Gaussian models: √ T ( θ -θ) I ∞ (θ) -1 ) where I ∞ (θ) = limT →∞ I T (θ) T .

  n, let [n] = {1, ..., n}. The symbol 1 {...} stands for an indicator function. When dealing with random variables, we write P(X = x) for a probability density, E[X] for an expectation, Var[X] for a variance (scalar of vector) and Cov[X, Y ] for a covariance (scalar or matrix). The symbols B(p) and N (µ, Σ) denote a Bernoulli distribution and a (possibly multivariate) Gaussian distribution. When we write log(x), we mean the natural (base-e) logarithm.

•

  T : transition matrix for Markov sampling • a, b: transition probabilities for Markov sampling • χ: minimum distance between a or b and {0, 1} (considered constant) • Y t : observations • η t : noise • ω 2 : variance of η t Estimation: • h: covariance time lag • h 0 : minimum covariance time lag for transition estimation • S(h): scaling matrix for covariance estimation • pq u : smallest coefficient of the scaling matrix Other: • g: standard Gaussian vector • Ψ ε (resp. Ψ η ): link between X (resp. η) and a standard Gaussian vector • L: random bilinear form • u: threshold in concentration inequalities • δ: small probability • Q Π : constant term in the conditional variance of Y • R Π (θ): varying term in the conditional variance of Y • ∆ Π (θ): deviation from the identity • γ (resp. γ u (θ)): signal-to-noise ratio in the lower bound (resp. the upper bound)

Indeed, the combination of discrete and Gaussian concentration inequalities as performed on page

(middle of right column) of the supplementary material forRao 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.

      Γ 0 (θ) Γ 0 (θ)θ 1 Γ 0 (θ)θ 2 • • • Γ 0 (θ)θ T -1 θ 1 Γ 0 (θ) Γ 0 (θ) Γ 0 (θ)θ 1 θ 2 Γ 0 (θ) θ 1 Γ 0 (θ) Γ 0 (θ) . . . . . .

The additional trace terms that appear when applying Lemma

(as opposed to the non-conditional version of Lemma 36) are absent from the papers byRao et al. [2017a,b], which is why we think their upper bound proofs are incomplete.

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.

we can simplify the expression of τ :

In this case, by Lemma 27, we can conclude:

C USEFUL LEMMAS

C.1 LINEAR ALGEBRA

The following set of results will sometimes be used in matrix calculations without explicit justifications.

Lemma 22 (Weyl's inequality). Let A and B be two n × n symmetric matrices. Then for all i we have:

In particular, λ min (A) + λ min (B) ≤ λ min (A + B).

Proof. See Horn and Johnson [2012, Theorem 4.3.1].

Lemma 23 (Ostrowski). Let S and A be two n × n matrices with S symmetric. For all i, there is a real number r i ∈ [s min (A) 2 , s max (A) 2 ] such that λ i (ASA ) = r i λ i (S), where s min (resp. s max ) denotes the minimum (resp. maximum) singular value.

Proof. See Horn and Johnson [2012, Theorem 4.5.9 and Corollary 4.5.11] Lemma 24 (Singular values of the Kronecker product). Let A and B be two matrices. Then

Proof. See Horn and Johnson [1994, Theorem 4.2.15].

Lemma 25. For any two matrices A and B, we have:

Proof. The Loewner order on symmetric matrices satisfies the following properties:

The first inequality is true because if x is a vector, x R (Q -P )Rx = (Rx) (Q -P )(Rx) ≥ 0 due to the Loewner positivity of Q -P . The second inequality can be directly deduced from the relation between the spectra of P and Q. Therefore, since A A is symmetric,

The proof for the other inequality is identical.

Lemma 26. Let A and B be two matrices with compatible sizes: then

We easily deduce:

Lemma 27 (Fano's method). Let θ 0 , ..., θ M be M + 1 parameters that are 2τ -separated w.r.t. a distance

and such that the average KL divergence between P θi and P θ0 is small enough

Then the minimax probability of an error at threshold τ satisfies:

≥ 1/2 is enough to obtain a minimax risk greater than α, as soon as M ≥ 3.

Lemma 28 (Chain rule for KL divergence). If P 0 and P 1 are probability densities on a product space X × Y with X discrete, then:

Proof. See Cover and Thomas [2012, Theorem 2.5.3].

Lemma 29 (KL divergence between Gaussians). The KL divergence between two multivariate Gaussian distributions

Proof. See Duchi [2007, page 13].

Proof. The proof of Levin and Peres [2017, Theorem 5.4] shows that with our assumptions,

From which we can deduce a sufficient condition for -mixing:

The result follows easily.

Lemma 33 (Chernoff inequality for Markov chains). Let (X t ) be an ergodic stationary Markov chain with finite state space X . We consider a function f : Chung et al. [2012, Theorem 3] Lemma 34 (Chernoff inequality for Markov chains under Doeblin condition). Under the hypotheses of the previous two Lemmas (32 and 33), if the parameters r and δ in the Doeblin condition are constants, then we have:

Proof. By Lemma 32, since r and δ are constants, the 1 8 -mixing time of (X t ) can be bounded by a constant

which we merge with the c 2 inside the exponential of Lemma 33. Lemma 36 (Hanson-Wright inequality: Gaussian case). Let A be a square matrix. If X and Y are two independent standard Gaussian vectors, we have:

Proof. See Vershynin [2018, Theorem 6.2.1] for the first inequality. We will see that it implies the second one. Let us define

We note that A F = A F and A 2 = A 2 . Applying the first inequality to X A X = X AY yields the expected result.

Lemma 37 (Conditional Hanson-Wright inequality). Let A be a random square matrix such that with probability 1-δ,

If X and Y are two independent standard Gaussian vectors independent of A, we have:

Proof. We start with the first case. Since A is a discrete random matrix with a finite set A of possible values,

Using independence between X and A gives us

We now split the set of feasible values A into

Since we assumed P(A ∈ A > ) = a∈A> P(A = a) ≤ δ, we get:

Unfortunately, Lemma 36 only lets us bound

(notice the change inside the expectation), which means we need an additional step. For a fixed a ∈ A ≤ , we use independence and normality to obtain

We are now ready to decompose, with the help of the union bound:

This implies:

By definition of A ≤ ,

The proof for X AY follows the same lines, except that we replace E[XX ] = I by E[XY ] = 0, which removes the trace term in the final expression.

Lemma 38 (Heuristic optimality of the signal-to-noise ratio). In the one-dimensional setting with full observations, the dependency of the error in 1 + σ 2 ω 2 is "coherent" with the asymptotic behavior of the MLE. Proof. Let us consider the case where D = 1 and p = 1, since we are mainly interested in the role of the parameters σ 2 and ω 2 . In this case, Theorem 2 argues that the error of any estimator should grow at least like γ = 1 + ω 2 σ 2 . We also note that in this simple scenario, Theorem 1 states that γ u ∝ γ .

We will compare this to the asymptotic error of the Maximum Likelihood Estimator (MLE) θ, which (for wellbehaved models) is given by the inverse of the Fisher information matrix. To make this statement more precise, we will invoke Douc et al. [2014, Proposition 2.14]. Let us verify the conditions:

• The process is stable, i.e. ρ(θ) < 1. We made sure of that by assuming θ 2 ≤ ϑ < 1.

• The sampling matrix Π t is constant across time. Although this assumption is not essential, it is true here since p = 1 and D = 1 hence Π t = I 1 .

• The model has the smallest possible dimension.

• s max (M ) (resp. s min (M ), s i (M )) its maximum (resp. minimum, i-th largest) singular value;

• M 1 = sup M x 1 x 1 = max j i |M i,j | its operator 1 norm, which is the maximum 1 norm of a column of M ;

• M 2 = sup M x 2

x 2 = |s max (M )| = λ max (M M ) its operator 2 norm, also known as the spectral norm;

| its operator ∞ norm, which is the maximum 1 norm of a row of M ;

• M F = vec(M ) 2 = Tr(M M ) its Frobenius norm;

• M max = vec(M ) ∞ = max i,j |M i,j | the maximum absolute value of its entries;

• ρ(M ) its spectral radius.

See [START_REF] Petersen | The matrix cookbook[END_REF] for a collection of inequalities relating all of these quantities. Given two real matrices A and B, we denote by

• A ⊗ B their Kronecker product;

• A B Hadamard (elementwise) product;

• A B or A B the (partial) Loewner order on symmetric matrices.

D.2 FREQUENT SYMBOLS

Here is a list of the most frequent symbols and their meaning.

Dimensions:

•