Conceptio › Archive › arXiv CS

arXiv CSopen access

When Does Model Collapse Occur in Structured Interactive Learning?

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

When Does Model Collapse Occur in Structured Interactive Learning? Yuchen Wu∗†

Kangjie Zhou∗ ‡

Weijie Su§

arXiv:2605.20151v1 [cs.LG] 19 May 2026

May 20, 2026

Abstract The proliferation of generative artificial intelligence has given rise to an interactive learning environment, where model parameters are continuously updated using not only data generated by natural processes, but also synthetic outputs produced by other models. This paradigm introduces two major challenges: (1) training data are no longer drawn exclusively from the target population, undermining a core assumption of classical statistical learning, and (2) model training processes become inherently correlated, as models interact with one another through repeated exposure to each other’s synthetic outputs in a potentially complex manner. Establishing reliable statistical inference in such structured interactive learning environments therefore remains an important open problem. In particular, there is growing concern about model collapse, a phenomenon in which the performance of generative models progressively degrades as they are trained on synthetic data produced by earlier model generations. Prior work on model collapse primarily focuses on a single model trained on its own output, failing to capture model performance in multi-model interactive settings. In this work, we fill this gap by investigating the performance of generative models in an interactive learning environment with general interaction patterns. In particular, we formalize model interactions using directed graphs and show that the occurrence of model collapse depends critically on the topology of the interaction graph. We further derive an explicit necessary and sufficient condition characterizing when model collapse occurs, and establish finite-sample results for linear regression and asymptotic guarantees for general M-estimators. We support our theoretical findings through extensive numerical experiments.

Contents 1 Introduction 1.1 Learning with self-generated content . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 A structured interactive learning environment . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Main contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 4 4 4 5 5

2 A directed graph framework for interactive learning 2.1 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 7

3 Results for linear regression 3.1 Interactive learning with linear regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 A matrix-form update equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Assumptions for the linear regression setting . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9 9 10 11 13

∗ The first two authors contribute equally to this work.

† School of Operations Research and Information Engineering, Cornell University;

Email: [email protected]. Email: [email protected]. § Department of Statistics and Data Science, University of Pennsylvania; Email: [email protected]. ‡ Department of Statistics, Columbia University;

4 Results for M-estimation 14 4.1 Interactive learning with general M-estimators . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.2 Assumptions for the M-estimation setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.3 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.4 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 5 Numerical experiments 19 5.1 Synthetic data experiment I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 5.2 Synthetic data experiment II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 5.3 Real data experiments with MNIST and CIFAR-10 . . . . . . . . . . . . . . . . . . . . . . . . 23 A Proofs of the main results 28 A.1 Proof of Theorem 3.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 A.2 Proof of Theorem 3.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 A.3 Proof of Theorem 4.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 A.4 Proof of Theorem 4.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 A.5 Proof of Theorem 4.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 B Technical lemmas for the linear regression setting 42 B.1 Proof of Lemma A.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 B.2 Proof of Lemma A.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 B.3 Proof of Lemma A.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 B.4 Proof of Lemma A.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 B.5 Proof of Lemma A.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 B.6 Proof of Lemma A.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 B.7 Proof of Lemma A.9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 C Technical lemmas for the M-estimation setting 52 ∗ . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 C.1 Consistency and asymptotic normality of βbt,µ C.2 Proof of Lemma A.10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 C.3 Proof of Lemma A.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 C.4 Auxiliary lemmas from empirical process theory . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Introduction

Generative modeling approaches are a powerful class of machine learning methods with broad applications (Yang et al., 2023; OpenAI, 2023; DeepMind, 2023; Abramson et al., 2024), driving numerous innovations across various fields such as image and text generation, drug discovery, and healthcare. Traditionally, generative models are trained on samples drawn from a target distribution. Given independent and identically distributed (i.i.d.) observations from this distribution, a probabilistic model is learned to approximate the underlying ground-truth distribution and subsequently generate new samples. The classical pipeline for training generative models has been extensively studied and is relatively well understood. However, its core assumptions may no longer hold in the modern data science era, where large-scale generative artificial intelligence (AI), such as GPT (OpenAI, 2023), Gemini (DeepMind, 2023), DALL-E (Betker et al., 2023) and Stable Diffusion (Podell et al., 2024), have revolutionized how data are generated and used. Generative models are increasingly being trained on not only data produced by human or natural processes, but also on synthetic data. We use several examples below to illustrate this trend. Evolution of online data ecosystem. Generative AI has become deeply integrated into our daily life. As a consequence, the internet is increasingly saturated with synthetic content that is indistinguishable from real data, and future training on data scraped from the web will inevitably be influenced by AI-generated data. Indeed, emerging evidence suggests that synthetic data has not only contaminated internet content (Sun et al., 2025; Spennemann, 2025), but has also infiltrated benchmark datasets (Alemohammad et al., 2024) and even human annotations produced via crowdsourcing (Veselovsky et al., 2023). 2

Accessing real data can be difficult. In certain applications, data generated by natural processes can be costly or hard to obtain. In such settings, synthetic data produced by generative AI offers a practical solution to mitigate data scarcity. For example, synthetic Electronic Health Record (EHR) data has been employed in scenarios where the original dataset is incomplete or inaccessible due to privacy concerns (Choi et al., 2017; Rankin et al., 2020; Vaid et al., 2023). Synthetic data also plays a central role in the development of Large Language Models (LLMs), enabling significant self-improvement through LLM-generated data (Huang et al., 2023; Lee et al., 2024). Beyond these settings, synthetic data has been widely adopted across diverse domains, including computer vision, finance, education, and agriculture. We refer the reader to Lu et al. (2023) for a comprehensive review. Knowledge distillation. Generative AI has demonstrated remarkable potential across a wide range of applications, yet its real-world deployment may be constrained by substantial computational demands. In particular, the high inference costs of these models can hinder adoption, especially in resource-constrained environments (e.g., mobile devices and edge computing platforms). Knowledge distillation offers a promising solution to address this challenge (Hinton et al., 2015; Xu et al., 2024). Specifically, the distillation process starts with a large “teacher” model, which typically has high capacity but can be expensive to deploy. During training, a smaller and more efficient “student” model is trained on the outputs of the teacher model. The goal is to build a model with fewer parameters and reduced inference-time costs, while maintaining performance that is comparable to the teacher model on selected tasks. Another feature of this synthetic learning regime is that the trained generative models are deployed to produce additional synthetic data, which may re-enter the training pool and be used to train subsequent generations of models. For example, a significant fraction of web-scraped data is now AI-generated, and such content may be incorporated into the training pipelines of future LLMs (Grattafiori et al., 2024; Achiam et al., 2023; Team et al., 2024). This process naturally induces interdependent training cycles whose dynamics differ significantly from those of classical regimes. We refer to Figure 1 for an illustration of this process. In particular, there is growing concern that continuously incorporating synthetic data into training may lead to model collapse (Shumailov et al., 2024), a phenomenon in which model performance degrades progressively with each training cycle until it becomes useless1 . Model collapse has been empirically observed in various applications (Hataya et al., 2023; Martı́nez et al., 2023a,b; Bohacek and Farid, 2023; Guo et al., 2024; Alemohammad et al., 2024; Bertrand et al., 2024), and has been rigorously characterized for a number of standard machine learning methods (Shumailov et al., 2024; Kazdan et al., 2025; Dohmatob et al., 2024a; Feng et al., 2025).

Figure 1: Illustration of the interdependent training cycle. In each iteration, models are updated using a mixture of synthetic and real data. Note that the synthetic data used for training may be generated by models from earlier generations, and not necessarily by the immediately preceding model. 1 The definition of model collapse varies across contexts, see Schaeffer et al. (2025) for a comprehensive discussion. In this work, we define model collapse as realized risk being significantly larger than the risk achieved when training on the same amount of naturally generated data. See Eq. (2) for a precise formulation.

1.1

Learning with self-generated content

Theoretical guarantees for learning from natural or human-generated data abound (van der Vaart, 2000), whereas comparatively fewer rigorous results exist for learning from synthetic data. Prior work in this direction has largely focused on settings in which a single model is recursively trained on natural data and its own outputs (Alemohammad et al., 2024; Shumailov et al., 2024; Martı́nez et al., 2023b,a; Dohmatob et al., 2024a,b; Gerstgrasser et al., 2024; Bertrand et al., 2024; Feng et al., 2025; Kazdan et al., 2025; He et al., 2025; Garg et al., 2025). Within this framework, one line of research examines the extreme setting where the training data at each iteration are entirely generated by the model from the previous iteration (Alemohammad et al., 2024; Shumailov et al., 2024, 2023). This regime is known to suffer from severe model collapse. Another line of work considers a more practically relevant setting, where the model is updated at each training cycle using a mixture of clean and synthetic data. This hybrid training regime has been theoretically shown to mitigate model collapse for a broad class of learning methods, covering linear models (Gerstgrasser et al., 2024; Dohmatob et al., 2025; Jain et al., 2024), Gaussian models (Kazdan et al., 2025; He et al., 2025), kernel methods (Jain et al., 2024; Kazdan et al., 2025), generalized linear models (Dey and Donoho, 2024; He et al., 2025), maximum likelihood estimation (Barzilai and Shamir, 2025), and softmax classifiers (Seddik et al., 2024). At a high level, model collapse tends to occur when the training process lacks sufficient access to natural data.

1.2

A structured interactive learning environment

The papers discussed in Section 1.1 focus primarily on the consequences of model training on a mixture of self-generated outputs and natural data, whereas relatively few works consider the more realistic setting in which models are trained on outputs generated also by other models. Such interactive learning regime captures a much broader class of practical applications. For instance, web-scale data generally consists of a mixture of real content and synthetic outputs produced by different generative AI models, as users may employ various AI tools to accomplish their tasks. In knowledge distillation, the outputs of a teacher model are used to train a separate student model, whose generated content may subsequently enter the web and be used to update the teacher model. In addition, in many real-world settings, interactions are structured and selective rather than all-to-all. For instance, in knowledge distillation, a student model learns exclusively from a preselected set of teacher models. In clinical applications, EHR foundation models trained on synthetic data typically rely on domainspecific synthetic EHR data rather than general web data. Moreover, language, regional, and licensing constraints may further restrict the data accessible to different AI developers. The complexity of real-world model interaction highlights the need for a generic theoretical framework to capture the learning dynamics of modern generative AI systems, thereby raising the following important question: How to assess model performance in a structured interactive learning environment, where models are repeatedly trained on each other’s synthetic outputs?

1.3

Main contributions

In this paper, we develop a general framework for analyzing interactive learning under arbitrary interaction patterns. We summarize our main contributions below: A structured interactive learning framework. We introduce a general theoretical framework that models interactive learning among generative models using directed graphs, thereby enabling the representation of arbitrary interaction patterns between models. As we show later, our framework is flexible enough to accommodate both synthetic and natural data sources. Within this framework, we assume that models are updated using architectures similar to those of previous generations. For example, models in the GPT1 through GPT-5 series descend from one another and partially inherit their predecessors’ architectures (Radford et al., 2018, 2019; Brown et al., 2020; Achiam et al., 2023). To the best of our knowledge, the only prior work on interactive learning is Vu et al. (2025), which studies linear regression under a dense interaction regime where every model learns from all others. Their results

suggest that interactive learning can improve model performance and homogenize model behaviors. Our framework generalizes their results by allowing for arbitrary model interaction patterns. Precise characterization of model collapse. Under mild regularity conditions, we establish explicit necessary and sufficient conditions for model collapse (as defined in Eq. (2)) in the interactive learning environment. In particular, we show that model collapse depends critically on the topology of the interaction graph, i.e., whether a model receives information from “unstable” models that lack a path to “stable” data sources, such as natural data. We present our findings formally in Section 2. Our results are consistent with prior work in single-model settings, where the occurrence of model collapse critically depends on the amount of natural data available to the model. Comprehensive theoretical guarantee. Our analysis covers a broad class of models used in modern machine learning. In particular, we provide a finite-sample characterization of model collapse under linear models, and an asymptotic characterization for general M-estimators learned via empirical risk minimization. In both settings, we show that under mild assumptions, model collapse is governed by the topology of the interaction graph. Our work provides theoretical insights into model behavior when multiple models co-evolve and the training datasets contain a mixture of natural data and synthetic data produced by different models. A key takeaway of our results is that model error remains controllable as long as the interaction structure is sufficiently rich to ensure that every model receives information (directly or indirectly) from reliable data sources, regardless of the number of training iterations. This finding aligns with the empirical evidence that generative models can still be trained successfully even when the training data is fully contaminated with synthetic samples (for instance, in knowledge distillation, models are trained on entirely synthetic labels).

1.4

Organization

The remainder of this paper is organized as follows. In Section 2, we introduce the interactive learning framework, using directed graphs to encode interactions between different models. In Section 3, we derive finite-sample results characterizing the occurrence of model collapse in the linear regression setting. We derive in Section 4 asymptotic results for general M-estimators and illustrate them with concrete examples, such as generalized linear models (GLMs). Finally, in Section 5, we perform extensive numerical experiments to validate our theoretical findings across a range of models and network structures. Proofs of all main theorems and supporting lemmas are deferred to the appendices.

1.5

Notation

For n ∈ N+ , we let [n] = {1, 2, · · · , n}. For two random objects X and Y , we say X ≡ Y if P(X = Y ) = 1, and write X ⊥ ⊥ Y if X is independent of Y . We use 1{·} to represent indicator functions. For r ∈ R, we denote by ⌊r⌋ the maximum integer that is no larger than r, and denote by ⌈r⌉ the minimum integer that is no smaller than r. For a, b ∈ R, we let a ∨ b = max{a, b} and a ∧ b = min{a, b}. For a matrix A, we denote by σmax (A) the largest eigenvalue of A, and denote by A† the pseudoinverse of A. For two symmetric matrices A and B, we write A ⪰ B if A − B is positive semidefinite, and A ≻ B if A − B is positive definite. For matrices A and B, we use A ⊗ B to represent the Kronecker product of A and B. We use ∥·∥2 to denote the Euclidean norm of a vector, and ∥·∥op and Tr(·) to denote the operator norm and trace of a matrix, respectively. For a set S in Euclidean space, we denote its interior by int S.

A directed graph framework for interactive learning

We consider an interactive learning environment with K models, where each model can learn from the outputs of other models. Let the collection of models be M = {µk : k ∈ [K]}. We represent their interactions using a directed graph G = (M, E), where the nodes correspond to the models, and the set of directed edges E encode their interactions. For µ, ν ∈ M, we use (ν, µ) to represent a directed edge connecting the origin ν to the destination µ. We say (ν, µ) ∈ E if and only if model µ learns from the outputs of model ν. Self-learning is captured by allowing self-loops (µ, µ) ∈ E whenever model µ learns from its own generated 5

data, thereby encompassing the setting of Section 1.1 (see Example 2.1 for more details). Our framework also accommodates natural data sources, which can be represented as nodes in G without incoming edges. Throughout this paper, we assume that all models share the common objective of learning a target population Ppop . Figure 2 presents an illustrative example of an interaction graph.

Figure 2: An example of an interaction graph. In this example, M = {µ1 , µ2 , µ3 , µ4 , µ5 , µ6 } and E = {(µ1 , µ3 ), (µ2 , µ3 ), (µ2 , µ5 ), (µ3 , µ6 ), (µ4 , µ4 ), (µ4 , µ5 )}. Nodes in Ml are colored red, and nodes in Mu are c nc colored blue. By definitions (3) and (4), we have M∞ l = {µ4 }, Ml = {µ4 , µ5 }, and Ml = {µ3 , µ6 }. For µ ∈ M, let Nµ,in = {ν : (ν, µ) ∈ E} denote the set of models whose outputs are used to train µ, and define Nmax = supµ∈M |Nµ,in |. We set Nµ,in = ∅ when µ does not update its parameters. We assume that interactions occur in rounds. Prior to any interactions, all models are initialized either using data drawn from the target distribution Ppop (corresponding to models trained on real data), or directly as Ppop itself, in which case the model represents a natural data source. For t ∈ N and µ ∈ M, we denote by Pbt,µ the learned distribution associated with model µ after the t-th interactive training cycle, where initialization is treated as the 0-th cycle. For each t ∈ N+ , during the t-th interactive training cycle and for every directed edge (ν, µ) ∈ E, model ν generates a dataset Dt,ν→µ according to its current distribution Pbt−1,ν , which is then passed to model µ. Model µ then aggregates the datasets {Dt,ν→µ : ν ∈ Nµ,in } to construct its updated distribution Pbt,µ . When Nµ,in = ∅, model µ receives no incoming information and therefore does not update, and Pbt,µ = Pb0,µ for all t ∈ N. Depending on whether a model learns from other models, we partition M into two disjoint subsets Ml and Mu , where Ml = {µ ∈ M : Nµ,in ̸= ∅} ,

Mu = {µ ∈ M : Nµ,in = ∅} .

(1)

For instance, for the interaction graph shown in Figure 2, we have Ml = {µ3 , µ4 , µ5 , µ6 } and Mu = {µ1 , µ2 }. Models in Mu do not update their parameters and can therefore be regarded as stable data sources. With these definitions, we summarize the interactive learning procedure as follows: - Initialization: For µ ∈ M, we either set Pb0,µ = Ppop (in this case, µ ∈ Mu represents a natural data source), or learn Pb0,µ by training on a dataset D0,µ drawn i.i.d. from Ppop . We denote by Mnature the collection of natural data sources, and assume that Mnature ⊆ Mu . Namely, Pb0,µ = Ppop for all µ ∈ Mnature ⊆ Mu , and Pb0,µ is learned from D0,µ for all µ ∈ M\Mnature . - Interactive training: Recursively for t ∈ N+ and each directed edge (ν, µ) ∈ E, model ν generates a dataset Dt,ν→µ consisting of i.i.d. samples from Pbt−1,ν . For each µ ∈ Ml , we construct Pbt,µ by training on the aggregated dataset ∪ν∈Nµ,in Dt,ν→µ . For all µ ∈ Mu , we set Pbt,µ = Pbt−1,µ . We study the evolution of Pbt,µ as t → ∞, focusing on characterizing its distance to the target distribution ∗ Ppop . When Pbt,µ is trained on nt,µ samples, a natural benchmark is the oracle estimator Pbt,µ , obtained by applying the same training algorithm to nt,µ i.i.d. samples drawn from Ppop . Specifically, we compare the

∗ ∗ risk rt,µ associated with Pbt,µ with the oracle risk rt,µ associated with Pbt,µ , where the precise definition of the risk depends on the specific context. We say that model collapse occurs for model µ if

lim sup t→∞

rt,µ ∗ = ∞. rt,µ

(2)

As discussed in Schaeffer et al. (2025), the definition of model collapse is context-dependent. In this work, we adopt the definition in Eq. (2). By construction, models in Mu retain fixed estimates and therefore do not evolve over time. For models ∗ in Ml , our results show that with reasonable risk choices, whether the risk ratio rt,µ /rt,µ remains bounded as t → ∞ depends critically on the structure of the interaction graph G. To present our results, we define M∞ l as the subset of nodes in Ml that are not reachable from any node in Mu : M∞ l = { µ ∈ Ml : ∀ν ∈ Mu , there is no directed path from ν to µ } .

(3)

Intuitively, M∞ l consists of the models in Ml that do not receive information from stable data sources. Taking the directed graph in Figure 2 as an example, M∞ l = {µ4 } since there is no directed path from {µ1 , µ2 } to µ4 . Our main results show that, for a model µ ∈ Ml , whether it collapses or not is completely determined by the existence of a path from some model in M∞ l to µ. Partitioning Ml into two disjoint : subsets, Mcl and Mnc l ∞ Mcl = M∞ l ∪ { µ ∈ Ml : ∃ν ∈ Ml , such that there is a directed path from ν to µ } ,

∞ ∞ Mnc l = { µ ∈ Ml \Ml : ∀ν ∈ Ml , there is no directed path from ν to µ } ,

(4)

we establish that across a wide range of scenarios, models in Mcl suffer from model collapse, whereas those c in Mnc l do not. For the network presented in Figure 2, this means that µ4 , µ5 ∈ Ml will collapse and nc µ3 , µ6 ∈ Ml will not collapse. We present a rigorous statement of these findings for linear regression in Section 3, followed by an extension to general M-estimators in Section 4. An important implication of our results is that, the set of models that will collapse in an interactive learning environment is uniquely determined by the graph structure, a property that is universal across many distinct learning scenarios. At a high level, one can think of M∞ l as the subset of models in Ml that do not receive information from stable data sources, causing their performance to vary widely. Consequently, models in Mcl are affected by these “unstable” data sources, leading to their collapse as the number of training cycles approaches infinity. Conversely, we demonstrate that models in Mnc l do not collapse because they exclusively receive information from “stable” sources. This dichotomy reveals the fundamental rationale behind model collapse in an interactive environment: it is not an inevitable consequence of learning from synthetic data, but rather hinges on contamination from unstable data sources. For example, as shown in Figure 2, model µ6 ∈ Mnc l does not collapse despite training exclusively on synthetic data generated by µ3 , since µ3 receives information only from stable data sources under our definition.

2.1

Examples

In this section, we present several illustrative examples of interaction graphs. Example 2.1 (Complete replacement training regime). In the complete replacement training regime, a single model is retrained at each cycle using data generated by itself with parameters from the previous round. The resulting training dynamics can thus be represented by a single-node interaction graph with M = {µ} and E = {(µ, µ)}. Example 2.2 (Accumulating training regime). In this setting, a single model is retrained at each cycle using all data generated in previous rounds. Specifically, we denote by Pbt the distribution estimate after training round t, with Pb0 representing the natural distribution. For t ∈ N+ , the estimate Pbt is obtained by b training on the aggregated dataset ∪t−1 i=0 Di , where each Di is sampled from Pi . We assume that training proceeds for a total of T rounds. We next show that this accumulating training regime can be represented by the interaction graph in Figure 3, where there is a directed edge from µa to µb for any a < b. To see this, note that µ1 in Figure 3 7

may represent a natural data source, so that Pbs,µ1 = Pb0 for all s ∈ N. If we further assume that µ1 generates the same dataset D0 and passes it to all µi with i ≥ 2 in every training cycle (indeed, D0 follows Pb0 = Pbs,µ1 ), then for all s ≥ 1, we have Pbs,µ2 = Pb1 . We may therefore assume that µ2 generates the dataset D1 and passes it to all µi with i ≥ 3 starting from round 2. Therefore, for all s ≥ 2 we have Pbs,µ3 = Pb2 . Continuing this construction inductively, we conclude that for all s ≥ t − 1 we have Pbs,µt = Pbt−1 . This holds for all t ∈ [T ], hence capturing Pbi for all i ∈ {0} ∪ [T − 1].

…… ……

Figure 3: Interaction graph that represents the accumulating training regime. Example 2.3 (Two patterns of model collapse). In many real-world scenarios, interactions are highly structured, leading to different patterns of model collapse. In this example, we consider a 5-node network G = (M, E), where M = {µi : 1 ≤ i ≤ 5} ,

E = {(µ1 , µ2 ), (µ5 , µ2 )} ∪ {(µi , µj ) : 3 ≤ i, j ≤ 5, i ̸= j} .

Namely, µ3 , µ4 and µ5 distill from one another, µ1 represents a natural data source, and µ2 learns from both natural data and synthetic data produced by model µ5 . See Figure 4 for an illustration. In this example, we observe two types of collapsing models: • Models that do not receive information from µ1 : Note that models µ3 , µ4 and µ5 do not receive information (directly or indirectly) from the only stable data source µ1 , hence they belong to the set M∞ l and will collapse as the number of training cycles t goes to infinity. • Models that receive information from µ1 : Model µ2 is trained on a mixture of natural data produced by µ1 and synthetic data produced by µ5 , and the proportion of natural data in the training set can be large2 . Despite this, µ2 still collapses because it receives information from an unstable source µ5 . This phenomenon is not captured by results in the single-model setting, where sufficient access to natural data can typically prevent model collapse.

Figure 4: The 5-node interaction graph that appears in Example 2.3. 2 As we will see later, it can be any positive constant strictly less than 1.

Example 2.4 (Model collapse can be sensitive to learning patterns). In this example, we show that model collapse can be sensitive to small changes in learning patterns. Specifically, consider the two interaction graphs in the left and right panels of Figure 5, which differ by only a single edge (the left panel is obtained by removing the edge µ2 → µ1 in the right panel). The left panel presents a hierarchical distillation learning pattern, where the model µ1 does not update its parameters and serves as a stable data source, and for every 2 ≤ i ≤ 5, the model µi distills from {µj : 1 ≤ j ≤ i − 1}. In addition, we assume that apart from distilling from µ4 and µ5 , model µ3 is also trained on data produced by a natural data source µ6 . Since µ1 is a stable data source and µ2 , µ3 , µ4 and µ5 all receive information from it, no model collapses in the left panel. In the right panel, we add a single edge (µ2 , µ1 ), which moves µ1 from the stable data source set Mu to the unstable data source set M∞ l . Consequently, µ2 , µ3 , µ4 and µ5 all receive information from unstable sources, and all models collapse in the right panel (except for the natural data source µ6 ). This example also suggests that models should be selective about the datasets they learn from, as blindly incorporating synthetic data may be detrimental.

Figure 5: The interaction graph that appears in Example 2.4.

Results for linear regression

In this section, we present our main results for linear regression. Linear regression is a classical statistical model that relates a response vector y ∈ Rn to a covariate matrix X ∈ Rn×d through a linear relationship: y = Xβ∗ + ε.

(5)

Here, β∗ ∈ Rd is an unknown parameter vector, and ε ∈ Rn represents measurement noise. Given observations (y, X), the goal is to estimate β∗ that defines the conditional distribution of y given X. The most standard method for fitting model (5) is ordinary least squares (OLS). Specifically, when X has full column rank, the OLS estimate for β admits an explicit form βbOLS = (X ⊤ X)−1 X ⊤ y.

(6)

OLS forms the foundation of many classical and modern regression methods. In this section, we study interactive learning in the linear regression setting using the OLS estimator. As in other discriminative modeling frameworks, our objective is to model the conditional distribution of the response given the predictors, rather than the full joint distribution. For linear regression, this means that Ppop in Section 2 denotes the groundtruth conditional distribution parameterized by β∗ , and Pbt,µ denotes an estimated conditional distribution parameterized by βbt,µ . In what follows, we carry out a finite-sample analysis for linear regression. Asymptotic results for general M-estimators, including OLS, are deferred to Section 4.

3.1

Interactive learning with linear regression

We now describe the interactive learning regime for linear regression, under a given interaction graph G = (M, E). This serves as a concrete example of the general setup introduced in Section 2. 9

At initialization, for each model µ ∈ M, we either set βb0,µ = β∗ , in which case µ represents a natural data source (and hence µ ∈ Mnature ⊆ Mu ), or compute the OLS estimator βb0,µ using nature-generated data (y0,µ , X0,µ ). In particular, assuming X0,µ has full column rank, we compute y0,µ = X0,µ β∗ + ε0,µ ∈ Rn0,µ , ⊤ ⊤ βb0,µ = (X0,µ X0,µ )−1 X0,µ y0,µ .

(7)

At the t-th training cycle for t = 1, 2, · · · , models in Mu remain unchanged, whereas each model µ in Ml P ⊤ updates its parameters using data generated by the models in Nµ,in (assuming ν∈Nµ,in Xt,ν→µ Xt,ν→µ has full column rank): X −1 X ⊤ ⊤ Xt,ν→µ yt,ν→µ , βbt,µ = Xt,ν→µ Xt,ν→µ (8) ν∈Nµ,in ν∈Nµ,in n yt,ν→µ = Xt,ν→µ βbt−1,ν + εt,ν→µ ∈ R t,ν→µ . Here, Xt,ν→µ is an nt,ν→µ × d data matrix. The estimator βbt,µ corresponds to the OLS solution obtained by stacking the datasets (yt,ν→µ , Xt,ν→µ ) over all ν ∈ Nµ,in . Note that the training datasets (yt,ν→µ , Xt,ν→µ ) for different (ν, µ) are not necessarily independent. Example 3.1. To illustrate this updating rule for the parameters {βbt,µ }t≥1,µ∈M , we consider a simple twonode network: M = {µ1 , µ2 } and E = {(µ1 , µ2 ), (µ2 , µ2 )}. In this setup, µ1 acts as a static natural data source, while µ2 learns from both µ1 and its own outputs from the previous training cycle. The general update equations (8) for t ≥ 1 then reduce to: βbt,µ1 = βbt−1,µ1 = · · · = βb0,µ1 = β∗ , −1 ⊤ ⊤ ⊤ ⊤ βbt,µ2 = Xt,µ Xt,µ1 →µ2 + Xt,µ Xt,µ2 →µ2 Xt,µ y + Xt,µ y , 1 →µ2 2 →µ2 1 →µ2 t,µ1 →µ2 2 →µ2 t,µ2 →µ2 where for each k ∈ {1, 2}, the dataset (Xt,µk →µ2 , yt,µk →µ2 ) ⊆ (Xt,µk , yt,µk ) is generated via: yt,µk →µ2 = Xt,µk →µ2 βbt−1,µk + εt,µk →µ2 ,

3.2

(Xt,µk →µ2 , εt,µk →µ2 ) ⊆ (Xt,µk , εt,µk ).

A matrix-form update equation

In the linear regression setting of Section 3.1, the update equations for individual model parameters are given in Eq. (8). To establish our main results, we next assemble these equations to form a single matrix-form update equation. We begin by introducing some notation. For t ∈ N and µ, ν ∈ M, we define  P −1 ⊤ ⊤ Xt,m→µ Xt,ν→µ Xt,ν→µ if t ∈ N+ , µ ∈ Ml , and ν ∈ Nµ,in , Xt,m→µ  m∈N  µ,in  Id if t ∈ N+ , µ ∈ Mu , and µ = ν, (9) Tt,µ,ν =  I if t = 0 and µ = ν, d   0d×d otherwise. The d × d matrix Tt,µ,ν quantifies the information transmitted from model ν to model µ during training cycle t. The identity matrix represents full information transfer (e.g., copying parameter values from the source model), whereas the zero matrix indicates no information transfer. When the rows of each design matrix are i.i.d. according to a common distribution with invertible second-moment matrix Σ, and the sample size is sufficiently large, the law of large numbers implies con⊤ centration: Xt,ν→µ Xt,ν→µ ≈ nt,ν→µ Σ. Applying this approximation for all models in Nµ,in , we obtain P P ⊤ X X t,m→µ t,m→µ ≈ ( m∈Nµ,in m∈Nµ,in nt,m→µ )Σ. Consequently, for t ∈ N+ and (ν, µ) ∈ E, we have P ∗ Tt,µ,ν ≈ nt,ν→µ /( m∈Nµ,in nt,m→µ )Id . This observation motivates the definition of Tt,µ,ν as the “population” version of Tt,µ,ν :  P (nt,ν→µ /( m∈Nµ,in nt,m→µ ))Id if t ∈ N+ , µ ∈ Ml , and ν ∈ Nµ,in ,    Id if t ∈ N+ , µ ∈ Mu , and µ = ν, ∗ Tt,µ,ν = (10)  I if t = 0 and µ = ν, d   0d×d otherwise. 10

Recall that M = {µ1 , µ2 , · · · , µK }. For t ∈ N, we define Tt ∈ RdK×dK (resp. Tt∗ ∈ RdK×dK ) as a K × K ∗ block matrix with the (i, j)-th block being Tt,µi ,µj (resp. Tt,µ ). We also define the collection of estimated i ,µj coefficients at time t as: ⊤ ⊤ ⊤ )⊤ ∈ RdK . βbt = (βbt,µ , βbt,µ , · · · , βbt,µ 1 2 K

(11)

We make the convention that βb−1 = (β∗⊤ , β∗⊤ , · · · , β∗⊤ )⊤ ∈ RdK . For µ ∈ Ml and t ∈ N+ , define the “mixed noise vector” vt,µ =

⊤ Xt,ν→µ Xt,ν→µ

−1 X ν∈Nµ,in

ν∈Nµ,in

⊤ Xt,ν→µ εt,ν→µ ∈ Rd .

(12)

For µ ∈ Mu and t ∈ N+ , we simply set vt,µ = 0d . When t = 0, we let v0,µ = 0d for µ ∈ Mnature , and for µ ∈ M \ Mnature , we define: ⊤ v0,µ = X0,µ X0,µ

−1

⊤ X0,µ ε0,µ ∈ Rd .

(13)

⊤ ⊤ ⊤ Let vt = (vt,µ , vt,µ , · · · , vt,µ )⊤ ∈ RdK . With the definitions in (9), (11), and (12), we obtain the following 1 2 K recursive formula for all t ∈ N:

βbt = Tt βbt−1 + vt .

(14)

Note that for t ∈ N, the matrix Tt∗ can be expressed as the Kronecker product of a K × K matrix Pt and a d × d identity matrix. For i, j ∈ [K], let Pt,µi ,µj denote the (i, j)-th entry of Pt , defined as follows:  P nt,ν→µ /( m∈Nµ,in nt,m→µ ) if t ∈ N+ , µ ∈ Ml , and ν ∈ Nµ,in ,    1 if t ∈ N+ , µ ∈ Mu , and µ = ν, Pt,µ,ν =  1 if t = 0 and µ = ν,   0 otherwise.

(15)

It follows that Tt∗ = Pt ⊗ Id .

3.3

Assumptions for the linear regression setting

In this section, we state the assumptions required to establish our linear regression results, beginning with the assumption of temporal independence. Assumption 3.1 (Temporal independence). For t ∈ N+ , let Xt = (Xt,ν→µ )(ν,µ)∈E denote the collection of design matrices from the t-th training cycle, and let Et = (εt,ν→µ )(ν,µ)∈E denote the corresponding collection of noise vectors. When t = 0, we define X0 = (X0,µ )µ∈M\Mnature and E0 = (ε0,µ )µ∈M\Mnature . We assume that (Xt , Et ) are mutually independent for all t ∈ N. Remark 3.1. Assumption 3.1 does not preclude deterministic design matrices: the covariates may be fixed or random, and in the deterministic case, independence holds trivially. We next state the assumptions imposed on the data transmitting pattern. In words, Assumption 3.2 states that at each round, every model ν generates data using the parameters fitted in the previous round and distributes a subset of this data to every model µ for which (ν, µ) ∈ E. Assumption 3.2 (Data transmitting pattern). For each t ∈ N+ and µ ∈ M, we denote by (yt,µ , Xt,µ ) a dataset produced by model µ in the t-th training cycle, generated according to the model fitted in the previous round: yt,µ = Xt,µ βbt−1,µ + εt,µ ∈ Rnt,µ . We assume that for every (ν, µ) ∈ E, the dataset (yt,ν→µ , Xt,ν→µ ) sent from ν to µ is a subset of (yt,ν , Xt,ν ). 11

Remark 3.2. Note that for ν ∈ M, the datasets {(yt,ν→µ , Xt,ν→µ ) : (ν, µ) ∈ E} may overlap with each other, and are thus not necessarily independent. We then assume that for each (ν, µ) ∈ E, a nontrivial portion of the data used by model µ to update its parameters is generated by model ν. Additionally, we assume that the number of training samples in each round does not grow too quickly. Assumption 3.3. There exists α ∈ (0, 1), such that for any t ∈ N+ and (ν, µ) ∈ E, nt,ν→µ ≥ α. m∈Nµ,in nt,m→µ

pt,ν→µ = P

Recall nt,µ is from Assumption 3.2. For t ∈ N+ , we define nt,max = supµ∈M nt,µ . We also define n0,max = supµ∈M\Mnature n0,µ . We assume that sups∈{0}∪[t] ns,max /t → 0 as t → ∞. P By definition, for all t ∈ N+ and Nµ,in ̸= ∅, it holds that ν∈Nµ,in pt,ν→µ = 1. Recall that Pt is defined in Eq. (15). For t ≥ s ≥ 0, we define Jt,s+1 = Pt Pt−1 · · · Ps+1 , with the convention that Jt,t+1 = IK . One can verify that each row of Jt,s+1 sums to 1. By definition, we know that a model ν1 in Mnc l receives information from at least one model in Mu . In other words, there exists a directed path from some ν2 ∈ Mu to ν1 . Let S(ν1 ; 0) denote a model in Mu that achieves the shortest path to ν1 , and let ℓ(ν1 ) denote the length of this path. We further denote this path by S(ν1 ; 0) → S(ν1 ; 1) → · · · → S(ν1 ; ℓ(ν1 )) = ν1 . Let T0 = supν1 ∈Mnc ℓ(ν1 ). Leveraging Assumption 3.3, we know that for all ν1 ∈ Mnc l , t ∈ N+ and t0 ≥ T0 , l ℓ(ν1 )

X ν2 ∈Mu

Jt+t0 ,t+1,ν1 ,ν2 ≥

Y i=1

pt+i,S(ν1 ;i−1)→S(ν1 ;i) ≥ αT0 .

The above lower bound further implies the existence of a constant ω ∈ (0, 1), such that for any t0 ≥ T0 , X X sup Jt+t0 ,t+1,ν1 ,ν2 = 1 − inf nc Jt+t0 ,t+1,ν1 ,ν2 ≤ 1 − αT0 ≤ ω ∈ (0, 1). (16) t∈N+ , ν1 ∈Mnc l

ν2 ∈Ml

t∈N+ , ν1 ∈Ml

ν2 ∈Mu

We next state the assumptions on the noise vectors. Specifically, we assume that they are independent of the design matrices and that their variances are bounded both above and below. Assumption 3.4. We impose the following assumptions on the noise vectors: 1. Recall that Et , Xt are defined in Assumption 3.1. We assume that for all t ∈ N, Et is a mean-zero random vector that is independent of Xt . ⊤ ⊤ ⊤ ⊤ ⊤ 2. For t ∈ N+ , we let εt = (ε⊤ t,µ1 , εt,µ2 , · · · , εt,µK ) . For t = 0, we let ε0 = (ε0,µ )µ∈M\Mnature . We assume that there exists ρ1 > 0, such that for all t ∈ N, Cov[εt ] ⪯ ρ1 I.

3. We assume that there exists ρ2 > 0, such that for all t ∈ N, Cov[εt ] ⪰ ρ2 I. We next state assumptions on the design matrices. In particular, we assume that Tt,µ,ν defined in Eq. (9) ∗ concentrates around its population counterpart Tt,µ,ν defined in Eq. (10). Recall that Xt,µ ∈ Rnt,µ ×d is ⊤ defined in Assumption 3.2. We further assume that, the operator norm of Xt,µ Xt,µ /nt,µ is bounded with high probability. Assumption 3.5. There exist δ, κ > 0 and γ1 , γ2 ∈ (0, 1), such that for any t ∈ N+ , µ ∈ M, and (ν, µ) ∈ E, ∗ P Tt,µ,ν − Tt,µ,ν ≥ δ ≤ γ1 , op ⊤ P Xt,µ Xt,µ /nt,µ op ≥ κ ≤ γ2 . 12

When t = 0, for any µ ∈ M\Mnature , we assume ⊤ P X0,µ X0,µ /n0,µ op ≥ κ ≤ γ2 . In addition, we assume δ, γ1 and γ2 are small enough such that δ + 2γ1 K − αK ∈ (0, c(α, K)), α δ + 2γ1 c1 = 1 + (1 − αK )1/K ∈ (0, 1), α 9K 3 γ2 K(δ + 2γ1 ) , > c2 = α K − αK+1 α(1 − αK )1/K+1 (1 − c1 )2 0 −1 √ 2T −1 2TX p j 0 ω + K γ1 + 2γ1 + δNmax Nmax ∈ (0, 1),

j=1

where c(α, K) = 1 − αK /2 depends only on α and K, and we recall that Nmax = supµ∈M |Nµ,in |, α ∈ (0, 1) is from Assumption 3.3, and (ω, T0 ) are from Eq. (16). Remark 3.3. When the rows of the design matrices are i.i.d. samples from a common distribution and the sample size is sufficiently large, Assumption 3.5 holds by standard matrix concentration inequalities, with the exact bound determined by the tail behavior of the underlying distribution. Finally, we assume that the design matrices have full column rank and that the inverse of their Gram matrices has finite expectation. Assumption 3.6. We assume the following: 1. For all t ∈ N+ , (ν1 , ν2 ) ∈ E, and µ ∈ M\Mnature , the matrices Xt,ν1 →ν2 and X0,µ have full column rank with probability one. 2. There exists a constant C0 > 0, such that h i ⊤ −1 sup E Tr (Xt,ν X ) ≤ C0 , t,ν →ν →ν 1 2 1 2 t∈N+ , (ν1 ,ν2 )∈E

sup µ∈M\Mnature

h i ⊤ E Tr (X0,µ X0,µ )−1 ≤ C0 .

Remark 3.4. By Srivastava and Vershynin (2013), if the rows of X ∈ Rn×p are independent and identically distributed with a finite (2 + w)-th moment for some w > 0 and an invertible covariance matrix, and if n ≥ C ′ d for some constant C ′ > 0, then E[(X ⊤ X)−1 ] exists and is finite.

3.4

Main results

∗ In this section, we present our main results for the linear regression setting. For comparison, let βbt,µ denote the OLS estimator computed using the same covariate matrix as βbt,µ , but with responses generated ∗ according to the ground-truth linear model (5). Specifically, βbt,µ is trained on the dataset obtained by ∗ ∗ stacking (Xt,ν→µ , yt,ν→µ ), where, unlike yt,ν→µ (which is generated using βbt−1,ν ), the responses yt,ν→µ are generated from the true conditional model: ∗ yt,ν→µ = Xt,ν→µ β∗ + εt,ν→µ .

In this case, ∗ βbt,µ = β∗ +

⊤ Xt,ν→µ Xt,ν→µ

ν∈Nµ,in

−1 X ν∈Nµ,in

where we recall that (Xt,ν→µ , εt,ν→µ ) are introduced in Section 3.1. 13

⊤ Xt,ν→µ εt,ν→µ ,

As is standard in linear regression, we use the mean squared error to quantify the discrepancy between the fitted and target distributions. Specifically, for t ∈ N and µ ∈ M, we define rt,µ = E[∥βbt,µ − β ∗ ∥22 ] and ∗ ∗ rt,µ = E[∥βbt,µ − β ∗ ∥22 ]. Recall from Eq. (2) that, we say a model µ collapses if and only if the limit superior of the ratio between these two quantities diverges as t → ∞. Our findings imply that, under the assumptions in Section 3.3, whether a model collapses depends only on its relative position within the interaction graph G. In particular, our first main theorem states that every model in Mcl experiences model collapse. Theorem 3.1. Under Assumptions 3.1 – 3.6, for every µ ∈ Mcl , we have E[∥β∗ − βbt,µ ∥22 ] = ∞. t→∞ E[∥β − β b∗ ∥2 ] ∗ lim

t,µ 2

Proof of Theorem 3.1. We prove Theorem 3.1 in Appendix A.1. Our second main theorem states that all models in Mnc l do not collapse, in the sense that the associated risk ratios stay bounded as t → ∞. Theorem 3.2. Under Assumptions 3.1 – 3.6, and additionally assuming that the model performances are comparable, in the sense that ∗ ∥22 ] sups∈{0}∪[t] supν∈M E[∥β∗ − βbs,ν <∞ t≥1 E[∥β∗ − βb∗ ∥2 ]

sup

t,µ 2

for all µ ∈ Mnc l ,

then for every µ ∈ Mnc l , we have lim sup t→∞

E[∥β∗ − βbt,µ ∥22 ] < ∞. E[∥β∗ − βb∗ ∥2 ] t,µ 2

Proof of Theorem 3.2. We prove Theorem 3.2 in Appendix A.2. Remark 3.5. The additional assumption in Theorem 3.2 requires that the MSEs of different models, when trained solely on naturally generated data, do not differ substantially. Otherwise, if one model is intrinsically much worse than another when trained on nature-produced data (for example, due to having significantly fewer training samples), and its outputs are used by another model for parameter updates, then we would expect the second model to collapse regardless of the communication pattern.

Results for M-estimation

In this section, we study general M-estimators defined as minimizers of the empirical risk induced by a loss function L. Specifically, we consider generative models defined by a vector β ∈ Rd and a function φ : Rd × R 7→ Rp . We denote by Pβ the distribution associated with β ∈ Rd , and assume that z ∼ Pβ admits the representation d

z = φ(β, ε), where ε ∼ µε for some distribution µε on R. Given z1 , z2 , · · · , zn ∼i.i.d. Pβ , we estimate β via empirical risk minimization (ERM): n 1X b β = arg min L(β, zi ), β∈B n i=1

(17)

where B ⊆ Rd . In the case of unconstrained optimization, we simply take B = Rd . Estimators defined by Eq. (17) are quite general and cover a broad class of model training procedures.

4.1

Interactive learning with general M-estimators

We now describe the interactive learning procedure for M-estimation under the interaction graph G = (M, E), where M = {µ1 , µ2 , · · · , µK }. As in earlier sections, β∗ ∈ Rd denotes the ground-truth parameter, and βbt,µ ∈ Rd represents the estimate retained by model µ after the t-th training cycle. At initialization, we set βb0,µ = β∗ for µ ∈ Mnature ⊆ Mu , and for the remaining models, βb0,µ is set to be the M-estimator obtained via solving the ERM problem (17) over a dataset Z0,µ ∈ Rn0,µ ×p . Here, the rows of Z0,µ are i.i.d. generated from the ground-truth distribution Pβ∗ . To be specific, for µ ∈ M\Mnature , n0,µ 1 X b L(β, z0,µ,i ), β0,µ = arg min β∈B n0,µ i=1

where z0,µ,i ∈ Rp denotes the i-th row of Z0,µ , generated according to z0,µ,i = φ(β∗ , ε0,µ,i ) for ε0,µ,i ∼i.i.d. µε . In the t-th training cycle for t ∈ N+ , models in Mu remain unchanged and we set βbt,µ = βbt−1,µ for all µ ∈ Mu . As for a model µ in Ml , we update its parameters based on data generated by models in Nµ,in : βbt,µ = arg min P β∈B

nt,ν→µ

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

zt,ν→µ,i ∼i.i.d. Pβbt−1,ν

L(β, zt,ν→µ,i ),

i=1

for i = 1, 2, · · · , nt,ν→µ ,

where zt,ν→µ,i is the i-th row of the matrix Zt,ν→µ ∈ Rnt,ν→µ ×p , and Zt,ν→µ represents the data passed from model ν to model µ during the t-th training cycle. We write zt,ν→µ,i = φ(βbt−1,ν , εt,ν→µ,i ) for εt,ν→µ,i ∼ µε . Analogous to Assumption 3.2 in the linear setting, we assume that Zt,ν→µ is a subset of Zt,ν ∈ Rnt,ν ×p indexed by It,ν→µ ⊆ [nt,ν ], i.e., the rows of Zt,ν→µ consist of {zt,ν,i : i ∈ It,ν→µ }. Here, zt,ν,i ∈ Rp denotes the i-th row of Zt,ν . We further assume that the rows of Zt,ν are sampled i.i.d. from Pβbt−1,ν , and that the matrices Zt,ν are mutually independent. Note that {Zt,ν→µ : (ν, µ) ∈ E} are not necessarily mutually independent and may overlap with each other.

4.2

Assumptions for the M-estimation setting

In this section, we summarize the assumptions needed to establish our main results for M-estimation, beginning with our requirement on data independence. Assumption 4.1. For all (t, µ) ∈ N+ × M or (t, µ) ∈ {0} × (M\Mnature ), we denote by zt,µ,i ∈ Rp the i-th row of Zt,µ ∈ Rnt,µ ×p . For t ∈ N+ , we assume that (zt,µ,i )i∈[nt,µ ] are i.i.d. sampled from Pβbt−1,µ , and we write zt,µ,i = φ(βbt−1,µ , εt,µ,i ). For µ ∈ M\Mnature , we assume that (z0,µ,i )i∈[n0,µ ] are i.i.d. sampled from Pβ∗ , with z0,µ,i = φ(β∗ , ε0,µ,i ). We further assume that all relevant εt,µ,i are mutually independent and follow µε . We work in an asymptotic regime where the ratios of sample sizes in different training cycles converge to fixed constants as the sample sizes grow to infinity, as formalized in Assumption 4.2. P Assumption 4.2. For t ∈ N+ , let nt = P µ∈M nt,µ denote the total number of samples used for training in the t-th training cycle, and let n0 = µ∈M\Mnature n0,µ denote the total number of samples used for initialization. We assume that the sample sizes tend to infinity simultaneously, and there exist two positive constants α, α > 0, such that the following holds: 1. For all t ∈ N+ and (ν, µ) ∈ E, we assume nt,ν→µ /nt → pt,ν→µ ≥ α as the sample sizes tend to infinity. For all µ ∈ M\Mnature , we assume n0,µ /n0 → p0,µ ≥ α as the sample sizes tend to infinity. 2. For all t, s ∈ N, we assume that nt /ns → bt,s ∈ [α, α] as the sample sizes tend to infinity. 3. For all t ∈ N+ and ν1 , ν2 ∈ M, we assume that the proportion of common data points used to train ν1 and ν2 in the t-th round converges as the sample sizes tend to infinity, i.e., 1 |{(t, ν, i) : ν ∈ Nν1 ,in , i ∈ It,ν→ν1 } ∩ {(t, ν, i) : ν ∈ Nν2 ,in , i ∈ It,ν→ν2 }| → qt,ν1 ∩ν2 ∈ [0, 1]. nt 15

Under Assumption 4.2, we define a matrix P t ∈ RK×K as follows: for i, j ∈ [K], let P t,µi ,µj denote the (i, j)-th entry of P t , defined by:  P p /( m∈Nµ,in pt,m→µ ) if t ∈ N+ , µ ∈ Ml , and ν ∈ Nµ,in ,    t,ν→µ 1 if t ∈ N+ , µ ∈ Mu , and µ = ν, P t,µ,ν = (18)  1 if t = 0 and µ = ν,   0 otherwise. Next, we assume that the model is well specified and satisfies Fisher consistency, in the sense that the corresponding population risk is uniquely minimized at the true data-generating parameter. We further assume that the population Hessian at β∗ is positive definite. These assumptions are imposed to ensure the consistency and asymptotic normality of the estimated parameters. Assumption 4.3. We assume that the mapping β 7→ L(β, z) is twice differentiable for any z ∈ Rp , and that the mapping (β1 , β2 ) 7→ E[L(β1 , φ(β2 , ε))] is continuous. For any β ∈ B, we assume that β = arg min E [L(β ′ , φ(β, ε))] ′ β ∈B

is the unique minimizer. In the above equation, the expectation is taken over ε ∼ µε . Further, we assume that the population Hessian is strictly positive at β∗ : E ∇2β L(β, φ(β∗ , ε))

β=β∗

≻ 0.

We also assume that for all β1 , β2 ∈ B: E[L(β1 , φ(β2 , ε))] < ∞,

E[∇β1 L(β1 , φ(β2 , ε))] < ∞,

E[∇2β1 L(β1 , φ(β2 , ε))] < ∞,

and in addition, E[∥∇β L(β, φ(β∗ , ε))|β=β∗ ∥22 ] < ∞. We additionally assume that either L is a convex function, or B is a compact subset of Rd . This condition guarantees that the ERM problem (17) has a well-defined solution. Assumption 4.4. We assume β∗ ∈ int B. In addition, we assume at least one of the following conditions hold: (1) L(β, z) is convex in β; (2) B is compact. We next state our key assumptions on the loss function L. In particular, we require that L and its Hessian belong to Glivenko-Cantelli function classes, while the Jacobian of L forms a Donsker class. We begin by recalling the definitions of Glivenko-Cantelli and Donsker classes. For more details, see Chapter 2 of Van Der Vaart and Wellner (1996). Definition 1 (Glivenko–Cantelli class). We say that a function class F is P -Glivenko–Cantelli, if n

1X a.s. f (Xi ) − E[f (X1 )] → 0 n f ∈F i=1 sup

as n → ∞, for Xi ∼i.i.d. P . Definition 2 (Donsker class). We say that a function class F is P -Donsker, if d

in ℓ∞ (F), n 1 X f (Xi ) − E[f (X1 )] Gn (f ) = √ n i=1

Gn → G

for f ∈ F.

Here, Xi ∼i.i.d. P , and G is a mean-zero Gaussian process with covariance X ∼ P, ∀f, g ∈ F , Cov G(f ), G(g) = Cov f (X), g(X) , known as the P -Brownian Bridge. 16

We now state our assumptions on the loss function L: Assumption 4.5. For any compact set Ω ⊆ Rd , we assume the following: 1. The function classes

ε 7→ L (β1 , φ(β2 , ε)) (β1 , β2 ) ∈ Ω2 ,

ε 7→ ∇2β1 L (β1 , φ(β2 , ε)) (β1 , β2 ) ∈ Ω2

are µε -Glivenko-Cantelli. 2. For any fixed β1 ∈ Rd , the function class ε 7→ ∇β1 L (β1 , φ(β2 , ε)) β2 ∈ Ω is µε -Donsker. Remark 4.1. Assumption 4.5 ensures that tools from empirical process theory (such as the uniform law of large numbers and the uniform central limit theorem) can be applied to analyze the asymptotic behavior of M-estimators in this multi-round interactive learning setting. Notably, these conditions are weaker than many of the regularity assumptions commonly imposed in the literature, and hold for a broad class of models, including several canonical generalized linear models (GLMs). See Section 4.4 for further details. Finally, we assume that the population risk has sufficient regularity to allow for certain interchange of differentiation and integration. Assumption 4.6. The mapping (β1 , β2 ) 7→ Eε∼µε [∇β1 L(β1 , φ(β2 , ε))] is continuously differentiable, and ∇β1 Eε∼µε [∇β1 L(β1 , φ(β2 , ε))] = Eε∼µε [∇2β1 L(β1 , φ(β2 , ε))]. Later in Section 4.4, we provide examples demonstrating that the assumptions listed above are satisfied by many commonly used M-estimators.

4.3

Main results

⊤ ⊤ ⊤ Following the notation in Section 3.2, we define βbt = (βbt,µ , βbt,µ , · · · , βbt,µ )⊤ ∈ RdK as the collection of 1 2 K learned parameters after the t-th training cycle, with the convention that βb−1 = (β∗⊤ , β∗⊤ , · · · , β∗⊤ )⊤ ∈ RdK . ∗ denote the M-estimator trained using the same number of samples as βbt,µ , except that the samples Let βbt,µ are drawn directly from Pβ∗ . Using standard arguments in asymptotic statistics (cf. Chapter 3.2 of Van ∗ is consistent and asymptotically normal as the sample Der Vaart and Wellner (1996)), we know that βbt,µ ∗ sizes tend to infinity (see Appendix C.1 for a proof). Further, the asymptotic covariance of βbt,µ is given by

−1 −1 V∗ = E ∇2β L(β, φ(β∗ , ε)) E ∇β L(β, φ(β∗ , ε))∇β L(β, φ(β∗ , ε))⊤ E ∇2β L(β, φ(β∗ , ε))

(19)

β=β∗

Recall the limiting proportions p0,µ , pt,ν→µ , and qt,µ1 ∩µ2 from Assumption 4.2, we define the covariance matrices Σ0 ∈ RdK×dK and Vt ∈ RdK×dK for t ∈ N+ as follows: For i, j ∈ [K], the (i, j)-th block of Σ0 is a d × d matrix proportional to V∗ : Σ0,µi ,µj =

1{µi = µj , µi ̸∈ Mnature } V∗ , p0,µi

(20)

and the (i, j)-th block of Vt is also a constant multiple of V∗ : qt,µi ∩µj 1{µi ∈ Ml , µj ∈ Ml } P Vt,µi ,µj = P V∗ . ( ν∈Nµ ,in pt,ν→µi )( ν∈Nµ ,in pt,ν→µj ) i

(21)

With these definitions, our first main result establishes the large-sample properties of βbt as the sample sizes grow to infinity. 17

Theorem 4.1. Under Assumptions 4.1-4.6, as the sample sizes grow to infinity, each βbt,µ is consistent and asymptotically normal. In particular, √ b d nt βt − βb−1 → N (0, Σt ) , where the asymptotic covariance Σt ∈ RdK×dK is defined recursively as follows: Σt = bt,t−1 (P t ⊗ Id )Σt−1 (P t ⊗ Id )⊤ + Vt ,

(22)

where we recall that P t is defined in Eq. (18), bt,t−1 is from Assumption 4.2, and Σ0 and Vt are defined via Eqs. (20) and (21), respectively. Proof of Theorem 4.1. We prove Theorem 4.1 in Section A.3. Building on the recursive formula (22), our second main result characterizes the limiting behavior of Σt as ∗ t → ∞. To be specific, we define rt,µ = Tr(Σt,µ,µ ) and rt,µ = Tr(V∗ ) as proxies for the mean squared errors ∗ ∗ of βbt,µ and βbt,µ , respectively, where we recall that βbt,µ denotes the benchmark estimator trained exclusively on natural data using the same algorithm and the same number of training samples. We then compute the risk ratio Tr(Σt,µ,µ ) rt,µ = ∗ rt,µ Tr(V∗ ) ∗ in the large-sample limit. Similarly to quantify the asymptotic relative risk of βbt,µ compared to that of βbt,µ as before, we say that a model µ collapses if the limit superior of this ratio diverges as t → ∞. Since the transition matrices P t are completely determined by the network structure, we should expect the same qualitative behavior as in the linear regression setting: As t → ∞, every model in Mcl collapses, while every model in Mnc l remains non-collapsing. We establish this rigorously in Theorem 4.2 below.

Theorem 4.2. Under Assumptions 4.1-4.6, the following holds: (a) For all µ ∈ Mcl :

rt,µ Tr(Σt,µ,µ ) lim ∗ = lim = ∞. t→∞ rt,µ t→∞ Tr(V∗ )

(b) For all µ ∈ Mnc l : lim sup t→∞

rt,µ Tr(Σt,µ,µ ) = lim sup < ∞. ∗ rt,µ Tr(V∗ ) t→∞

Proof of Theorem 4.2. We prove Theorem 4.2 in Section A.4.

4.4

Examples

In this section, we present several examples of (L, φ) pairs that satisfy Assumptions 4.3-4.6. In particular, we focus on the widely used generalized linear models (GLMs). Specifically, for a GLM with true parameter β∗ ∈ Rd , the data z = (x, y) satisfies the following condition: Eβ∗ [y|x] = A′ (β∗⊤ x),

(23)

where A is a strictly convex and twice continuously differentiable cumulant function, x represents the vector of covariates, and y denotes the response variable. We present a few examples below. (i) Linear regression: A(ξ) = ξ 2 /2, and y = β∗⊤ x + ε,

ε ⊥⊥ x, E[ε] = 0, E[ε2 ] < ∞.

(ii) Logistic regression: A(ξ) = log(1 + exp(ξ)), and y | x ∼ Bernoulli 18

exp(β∗⊤ x) 1 + exp(β∗⊤ x)

(iii) Poisson regression: A(ξ) = exp(ξ), and ⊤ y | x ∼ Poisson eβ∗ x .

We consider the following standard loss function L for GLMs: L(β, z) = −yβ ⊤ x + A(β ⊤ x).

(24)

Below, we show that under mild moment conditions on the covariates, Assumptions 4.3–4.6 are satisfied for these specific GLMs. In Proposition A.11 in the appendix, we present some general conditions under which a GLM satisfies Assumptions 4.3–4.6. Theorem 4.3. The following is true: 4

(a) For linear regression, Assumptions 4.3–4.6 are satisfied, provided that E[∥x∥2 ] < ∞. (b) For logistic regression, Assumptions 4.3-4.6 are satisfied, provided that E[∥x∥32 ] < ∞. (c) For Poisson regression, Assumptions 4.3-4.6 are satisfied, provided that E[exp(R∥x∥2 )] < ∞ for any R > 0. Proof of Theorem 4.3. We prove Theorem 4.3 in Section A.5.

Numerical experiments

In this section, we present several numerical experiments to validate our theoretical findings, including synthetic experiments on linear regression, M-estimation, and learning with non-convex single-index models, as well as real data experiments with MNIST and CIFAR-10.

5.1

Synthetic data experiment I

Setting and algorithms. In our first experiment, we consider two interaction graphs presented in Figure 6. In the left panel, G = (M, E) contains 5 nodes, where M = {µ1 , µ2 , µ3 , µ4 , µ5 } and E = {(µ1 , µ2 ), (µ2 , µ5 ), (µ3 , µ4 ), (µ3 , µ5 ), (µ4 , µ3 )}. In this interaction graph, we have Ml = {µ2 , µ3 , µ4 , µ5 } and Mu = {µ1 }. Moreover, Mcl = {µ3 , µ4 , µ5 } and Mnc l = {µ2 }. Models in Mu are colored blue and those in Ml are colored red. In the right panel of Figure 6, we consider a slightly more complicated interaction graph that contains 8 nodes, where M = {µi : 1 ≤ i ≤ 8} and E = {(µ1 , µ3 ), (µ2 , µ4 ), (µ3 , µ4 ), (µ4 , µ8 ), (µ7 , µ8 )} ∪ {(µi , µj ) : 5 ≤ i, j ≤ 7, i ̸= j}. By definition, we have Ml = {µi : 3 ≤ i ≤ 8} and Mu = {µ1 , µ2 }. In addition, Mcl = {µ5 , µ6 , µ7 , µ8 } and Mnc l = {µ3 , µ4 }.

Figure 6: Interaction graphs considered in synthetic data experiment I. For edges in both examples, we use a common sample size by setting all nt,ν→µ and n0,µ equal to a fixed value nsample ∈ N+ . Here, we set nsample = 1000, d = 5, and train the models in M according to 19

the interaction graphs specified in Figure 6 for T = 50 rounds. We further assume that samples associated with different edges are independent of each other. We generate the ground-truth parameter β∗ randomly at the beginning of each experiment, and hold it fixed throughout the subsequent training cycles. We conduct numerical simulations for three models: (a) linear regression as described in Section 3; (b) logistic regression as described in Section 4.4; and (c) a single-index model with a quadratic link function, where z = (x, y),

y = (β∗⊤ x)2 + ε, and L(β, z) =

2 1 y − (β ⊤ x)2 . 2

Throughout, we assume that the noise in both the linear regression and single-index models follows a standard Gaussian distribution. The key difference between model (c) and the other two models is that its loss function is non-convex, making the global minimizer of the empirical risk computationally challenging to obtain. In our experiments, we use the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm to minimize the empirical risk and obtain the estimator βbt,µ for model (c). For models (a) and (b), we fit them using standard GLM implementations available in common Python libraries. Results of the experiments. compute the risk ratio

For both interaction graphs in Figure 6 and each model µ ∈ Ml , we rt,µ E[∥βbt,µ − β∗ ∥22 ] = . ∗ ∗ − β ∥2 ] rt,µ E[∥βbt,µ ∗ 2

∗ We plot the risk ratios as a function of t for t ∈ {1, · · · , 50}. The mean squared errors of βbt,µ and βbt,µ are estimated via Monte Carlo simulation using 1000 independent trials, with each trial based on 1000 independent samples. We present and discuss the results for the two interaction graphs in Figure 6 separately below. Experimental results for linear regression are presented in Figure 7, those for logistic regression in Figure 8, and those for the single-index model in Figure 9.

– Left panel of Figure 6: Our theoretical results in Sections 3 and 4 imply that models µ3 , µ4 and µ5 will collapse as t → ∞, while µ2 will not. We confirm this prediction in the left panels of Figures 7-9 for settings (a)-(c), all of which exhibit the same qualitative behavior. From this figure, we see that the risk ratio for µ2 remains bounded because it belongs to Mnc l and receives information directly from µ1 ∈ Mu . In contrast, the risk ratios for µ3 , µ4 and µ5 exhibit linear growth with respect to t, confirming that they will collapse as t → ∞. Notably, the risk ratios of µ3 and µ4 are nearly identical and they grow more rapidly than that of µ5 . This is because that µ3 and µ4 are confined to the loop µ3 ↔ µ4 , whereas µ5 benefits from information propagated along the directed path µ1 → µ2 → µ5 . – Right panel of Figure 6: We similarly plot the risk ratios for µ3 and µ4 (which learn only from stable data sources), µ8 (which learns from both stable and unstable data sources), and µ5 , µ6 , and µ7 (which are confined to an isolated synthetic loop hence learn only from unstable data sources), and observe the same qualitative behavior as in the first interaction graph. We additionally observe that, although neither µ3 nor µ4 collapses, the risk ratio of µ4 is slightly larger than that of µ3 . This is because that µ3 learns only from the natural data source µ1 , whereas µ4 also distills information from µ3 .

5.2

Synthetic data experiment II

In our second experiment, we compare the two interaction graphs presented in Example 2.4. Similar to the experiments in Section 5.1, we plot the evolution of risk ratios for models in the two interaction graphs under (a) linear regression, (b) logistic regression, and (c) a single-index model with a quadratic link function. The remaining settings are identical to those in Section 5.1. The numerical results are presented in Figure 10, from which we observe a sharp difference between the two graphs in terms of model stability: adding a single edge µ2 → µ1 to the first graph in Example 2.4 destroys its stable structure. Specifically, in the left panel of Figure 10, the models µi for 1 ≤ i ≤ 5 do not collapse, with risk ratios remaining bounded by 5. In contrast, in the right panel, all risk ratios grow approximately linearly with the number of training rounds, exceeding 10 after only a few iterations, and 20

Linear Regression 50

µ = µ3 µ = µ4 µ = µ5 µ = µ6 µ = µ7 µ = µ8

17.5 15.0 Risk Ratio

Risk Ratio

Linear Regression

20.0

µ = µ2 µ = µ3 µ = µ4 µ = µ5

30 20

12.5 10.0 7.5 5.0

2.5 0 0

20 30 Training Round

Figure 7: Plots of the risk ratios under the linear regression setting over the first 50 training cycles. The left panel shows results for the 5-node interaction graph in Figure 6, while the right panel shows results for the 8-node interaction graph in the same figure. Logistic Regression

60 50

µ = µ3 µ = µ4 µ = µ5 µ = µ6 µ = µ7 µ = µ8

20 Risk Ratio

40 Risk Ratio

Logistic Regression 25

µ = µ2 µ = µ3 µ = µ4 µ = µ5

30 20

10 0 0

20 30 Training Round

Figure 8: Plots of the risk ratios under the logistic regression setting over the first 50 training cycles. The left panel shows results for the 5-node interaction graph in Figure 6, while the right panel shows results for the 8-node interaction graph in the same figure. Non-convex Single-index Model 50

µ = µ3 µ = µ4 µ = µ5 µ = µ6 µ = µ7 µ = µ8

17.5 15.0 Risk Ratio

Risk Ratio

Non-convex Single-index Model

20.0

µ = µ2 µ = µ3 µ = µ4 µ = µ5

30 20

12.5 10.0 7.5 5.0

2.5 0 0

20 30 Training Round

Figure 9: Plots of the risk ratios under the non-convex single index model setting over the first 50 training cycles. The left panel shows results for the 5-node interaction graph in Figure 6, while the right panel shows results for the 8-node interaction graph in the same figure. with the largest surpassing 80 after 50 rounds. Once again, this confirms that model collapse can be highly sensitive to the underlying learning patterns.

Linear Regression (Stable Graph vs Unstable Graph) 5.0

4.5

70 60

Risk Ratio

4.0 µ = µ1 µ = µ2 µ = µ3 µ = µ4 µ = µ5

3.5 3.0 2.5

50 40 30

2.0

1.5

1.0

0 0

Training Round

Logistic Regression (Stable Graph vs Unstable Graph) 80 5

70 µ = µ1 µ = µ2 µ = µ3 µ = µ4 µ = µ5

Risk Ratio

60 50 40 30 20

10 1

0 0

Training Round

Non-convex Single-index Model (Stable Graph vs Unstable Graph) 5.0

4.5

Risk Ratio

4.0

60 µ = µ1 µ = µ2 µ = µ3 µ = µ4 µ = µ5

3.5 3.0 2.5

50 40 30

2.0

1.5

1.0

0 0

Training Round

Figure 10: Comparison of risk ratios for models in the two interaction graphs of Example 2.4. Top: linear regression. Middle: logistic regression. Bottom: Single-index model with quadratic link function. Left: the stable graph represented by the left panel of Figure 5. Right: the unstable graph represented by the right panel of Figure 5.

5.3

Real data experiments with MNIST and CIFAR-10

In this section, we present numerical experiments on some real-world datasets to validate our theoretical findings on model collapse. We first conduct a real data experiment using the MNIST dataset (LeCun et al., 1998). We again consider the interaction graphs in Example 2.4, and employ Generative Adversarial Networks (GANs) to train generative models for producing synthetic samples. As before, we assume that all nt,ν→µ and n0,µ are equal to a fixed value nsample ∈ N+ , and set nsample = 5000. We further assume that samples associated with different edges are mutually independent. To evaluate the quality of the generated samples, we compute the Fréchet Inception Distance (FID), a widely used metric for assessing the quality of generated images. Lower FID values indicate better image quality. Note that this setting does not fall within the framework proposed in Section 4, since the FID value is not the loss function used during training. In addition, to mimic the continuing-training procedure used in practice, we initialize each round of model training using the parameters obtained from the previous round. Although this example is not covered by our theoretical results, we qualitatively observe the same phenomena, suggesting that our theory may extend well beyond the current assumptions. We report the model FID ratios as a function of the number of rounds in Figure 11. Note that in the left panel of Figure 11, the FID ratios initially decrease and then stabilize as the training rounds proceed, potentially due to the use of previous-round parameters for model initialization. In contrast, in the right panel of the same figure the FID ratios increase with the training round. MNIST

MNIST

1.0

µ2 µ3 µ4 µ5

0.9

0.7

16 14 12 FID ratio

FID ratio

0.8

0.6

µ1 µ2 µ3 µ4 µ5

10 8

0.5

0.4

4 2

0.3 0

20 30 Training Round

Figure 11: FID ratios achieved by GANs trained on MNIST over 50 rounds in the interactive learning setting. The left panel shows results for the first interaction graph in Example 2.4, and the right panel corresponds to the second interaction graph. We also conduct the same experiment on the CIFAR-10 dataset (Krizhevsky et al., 2009) and observe similar behavior. The corresponding results are presented in Figure 12.

CIFAR − 10

CIFAR − 10 µ2 µ3 µ4 µ5

1.2

17.5 15.0 FID ratio

FID ratio

1.0 0.8 0.6

µ1 µ2 µ3 µ4 µ5

20.0

12.5 10.0 7.5

0.4

5.0 2.5

0.2

0.0 0

20 30 Training Round

Figure 12: FID ratios achieved by GANs trained on CIFAR-10 over 50 rounds in the interactive learning setting. Similar to Figure 11, the left panel shows results for the first interaction graph in Example 2.4, and the right panel corresponds to the second interaction graph.

References Josh Abramson, Jonas Adler, Jack Dunger, Richard Evans, Tim Green, Alexander Pritzel, Olaf Ronneberger, Lindsay Willmore, Andrew J Ballard, Joshua Bambrick, et al. Accurate structure prediction of biomolecular interactions with alphafold 3. Nature, 630(8016):493–500, 2024. Josh Achiam, Steven Adler, Sandhini Agarwal, Lama Ahmad, Ilge Akkaya, Florencia Leoni Aleman, Diogo Almeida, Janko Altenschmidt, Sam Altman, Shyamal Anadkat, et al. Gpt-4 technical report. arXiv preprint arXiv:2303.08774, 2023. Sina Alemohammad, Josue Casco-Rodriguez, Lorenzo Luzi, Ahmed Imtiaz Humayun, Hossein Babaei, Daniel LeJeune, Ali Siahkoohi, and Richard Baraniuk. Self-consuming generative models go MAD. In The Twelfth International Conference on Learning Representations, 2024. Daniel Barzilai and Ohad Shamir. When models don’t collapse: On the consistency of iterative MLE. arXiv preprint arXiv:2505.19046, 2025. URL https://arxiv.org/abs/2505.19046. Quentin Bertrand, Joey Bose, Alexandre Duplessis, Marco Jiralerspong, and Gauthier Gidel. On the stability of iterative retraining of generative models on their own data. In The Twelfth International Conference on Learning Representations, 2024. James Betker, Gabriel Goh, Li Jing, Tim Brooks, Jianfeng Wang, Linjie Li, Long Ouyang, Juntang Zhuang, Joyce Lee, Yufei Guo, et al. Improving image generation with better captions. Computer Science. https://cdn. openai. com/papers/dall-e-3. pdf, 2(3):8, 2023. Matyas Bohacek and Hany Farid. Nepotistically trained generative-ai models collapse. arXiv preprint arXiv:2311.12202, 2023. Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. Language models are few-shot learners. Advances in neural information processing systems, 33:1877–1901, 2020. Edward Choi, Siddharth Biswal, Bradley Malin, Jon Duke, Walter F. Stewart, and Jimeng Sun. Generating multi-label discrete patient records using generative adversarial networks. In Finale Doshi-Velez, Jim Fackler, David Kale, Rajesh Ranganath, Byron Wallace, and Jenna Wiens, editors, Proceedings of the 2nd Machine Learning for Healthcare Conference, volume 68 of Proceedings of Machine Learning Research, pages 286–305. PMLR, 18–19 Aug 2017. Google DeepMind. Gemini: Google’s next-generation ai models, 2023. technology/ai/introducing-gemini. 24

URL https://blog.google/

Apratim Dey and David Donoho. Universality of the π 2 /6 pathway in avoiding model collapse. arXiv preprint arXiv:2410.22812, 2024. Elvis Dohmatob, Yunzhen Feng, and Julia Kempe. Model collapse demystified: The case of regression. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024a. Elvis Dohmatob, Yunzhen Feng, Pu Yang, Francois Charton, and Julia Kempe. A tale of tails: Model collapse as a change of scaling laws. In Forty-first International Conference on Machine Learning, 2024b. URL https://openreview.net/forum?id=KVvku47shW. Elvis Dohmatob, Yunzhen Feng, Arjun Subramonian, and Julia Kempe. Strong model collapse. In The Thirteenth International Conference on Learning Representations, 2025. Yunzhen Feng, Elvis Dohmatob, Pu Yang, Francois Charton, and Julia Kempe. Beyond model collapse: Scaling up with synthesized data requires verification. In The Thirteenth International Conference on Learning Representations, 2025. URL https://openreview.net/forum?id=MQXrTMonT1. Anvit Garg, Sohom Bhattacharya, and Pragya Sur. Preventing model collapse under overparametrization: Optimal mixing ratios for interpolation learning and ridge regression. arXiv preprint arXiv:2509.22341, 2025. Matthias Gerstgrasser, Rylan Schaeffer, Apratim Dey, Rafael Rafailov, Henry Sleight, John Hughes, Tomasz Korbak, Rajashree Agrawal, Dhruv Pai, Andrey Gromov, Daniel A. Roberts, Diyi Yang, David L. Donoho, and Sanmi Koyejo. Is model collapse inevitable? Breaking the curse of recursion by accumulating real and synthetic data. In Proceedings of the First Conference on Language Modeling (COLM), 2024. Aaron Grattafiori, Abhimanyu Dubey, Abhinav Jauhri, Abhinav Pandey, Abhishek Kadian, Ahmad AlDahle, Aiesha Letman, Akhil Mathur, Alan Schelten, Alex Vaughan, et al. The llama 3 herd of models. In Neural Information Processing Systems. Curran Associates, 2024. Yanzhu Guo, Guokan Shang, Michalis Vazirgiannis, and Chloé Clavel. The curious decline of linguistic diversity: Training language models on synthetic text. Findings of the Association for Computational Linguistics: NAACL 2024, pages 3589–3604, 2024. Ryuichiro Hataya, Han Bao, and Hiromi Arai. Will large-scale generative models corrupt future datasets? 2023 IEEE/CVF International Conference on Computer Vision (ICCV), pages 20498–20508, 2023. Hengzhi He, Shirong Xu, and Guang Cheng. Golden ratio weighting prevents model collapse. arXiv preprint arXiv:2502.18049, 2025. Geoffrey Hinton, Oriol Vinyals, and Jeff Dean. Distilling the knowledge in a neural network. arXiv preprint arXiv:1503.02531, 2015. Jiaxin Huang, Shixiang Shane Gu, Le Hou, Yuexin Wu, Xuezhi Wang, Hongkun Yu, and Jiawei Han. Large language models can self-improve. In The 2023 Conference on Empirical Methods in Natural Language Processing, 2023. Ayush Jain, Andrea Montanari, and Eren Sasoglu. Scaling laws for learning with real and surrogate data. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024. Joshua Kazdan, Rylan Schaeffer, Apratim Dey, Matthias Gerstgrasser, Rafael Rafailov, David L Donoho, and Sanmi Koyejo. Collapse or thrive? perils and promises of synthetic data in a self-generating world. In Proceedings of the 42nd International Conference on Machine Learning (ICML), 2025. Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009. Yann LeCun, Corinna Cortes, and Christopher JC Burges. The mnist database of handwritten digits. http://yann.lecun.com/exdb/mnist/, 1998.

Harrison Lee, Samrat Phatale, Hassan Mansoor, Kellie Ren Lu, Thomas Mesnard, Johan Ferret, Colton Bishop, Ethan Hall, Victor Carbune, and Abhinav Rastogi. RLAIF: Scaling reinforcement learning from human feedback with AI feedback, 2024. Yingzhou Lu, Minjie Shen, Huazheng Wang, Xiao Wang, Capucine van Rechem, Tianfan Fu, and Wenqi Wei. Machine learning for synthetic data generation: a review. arXiv preprint arXiv:2302.04062, 2023. Gonzalo Martı́nez, Lauren Watson, Pedro Reviriego, José Alberto Hernández, Marc Juárez, and Rik Sarkar. Combining generative artificial intelligence (ai) and the internet: Heading towards evolution or degradation? arXiv preprint arXiv:2303.01255, 2023a. Gonzalo Martı́nez, Lauren Watson, Pedro Reviriego, José Alberto Hernández, Marc Juárez, and Rik Sarkar. Towards understanding the interplay of generative artificial intelligence and the internet. In Epi UAI, 2023b. OpenAI. Gpt-4 technical report, 2023. URL https://openai.com/research/gpt-4. Dustin Podell, Zion English, Kyle Lacey, Andreas Blattmann, Tim Dockhorn, Jonas Müller, Joe Penna, and Robin Rombach. Sdxl: Improving latent diffusion models for high-resolution image synthesis. In The Twelfth International Conference on Learning Representations, 2024. Alec Radford, Karthik Narasimhan, Tim Salimans, Ilya Sutskever, et al. Improving language understanding by generative pre-training. 2018. Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al. Language models are unsupervised multitask learners. OpenAI blog, 1(8):9, 2019. Debbie Rankin, Michaela Black, Raymond Bond, Jonathan Wallace, Maurice Mulvenna, and Gorka Epelde. Reliability of supervised machine learning using synthetic data in health care: model to preserve privacy for data sharing. JMIR medical informatics, 8(7):e18910, 2020. Rylan Schaeffer, Joshua Kazdan, Alvan Caleb Arulandu, and Sanmi Koyejo. Position: Model collapse does not mean what you think. arXiv preprint arXiv:2503.03150, 2025. Mohamed El Amine Seddik, Suei-Wen Chen, Soufiane Hayou, Pierre Youssef, and Merouane Abdelkader DEBBAH. How bad is training on synthetic data? a statistical analysis of language model collapse. In First Conference on Language Modeling, 2024. Ilia Shumailov, Zakhar Shumaylov, Yiren Zhao, Yarin Gal, Nicolas Papernot, and Ross Anderson. The curse of recursion: Training on generated data makes models forget. arXiv preprint arXiv:2305.17493, 2023. Ilia Shumailov, Zakhar Shumaylov, Yiren Zhao, Nicolas Papernot, Ross Anderson, and Yarin Gal. Ai models collapse when trained on recursively generated data. Nature, 631:755 – 759, 2024. Dirk HR Spennemann. Delving into: The quantification of ai-generated content on the internet (synthetic data). arXiv preprint arXiv:2504.08755, 2025. Nikhil Srivastava and Roman Vershynin. Covariance estimation for distributions with 2 + ϵ moments. The Annals of Probability, 41(5):3081–3111, 2013. Zhen Sun, Zongmin Zhang, Xinyue Shen, Ziyi Zhang, Yule Liu, Michael Backes, Yang Zhang, and Xinlei He. Are we in the ai-generated text world already? quantifying and monitoring aigt on social media. In Proceedings of the 63rd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 22975–23005, 2025. Gemini Team, Petko Georgiev, Ving Ian Lei, Ryan Burnell, Libin Bai, Anmol Gulati, Garrett Tanzer, Damien Vincent, Zhufeng Pan, Shibo Wang, et al. Gemini 1.5: Unlocking multimodal understanding across millions of tokens of context. arXiv preprint arXiv:2403.05530, 2024.

Akhil Vaid, Ashwin Sawant, Mayte Suarez-Farinas, Juhee Lee, Sanjeev Kaul, Patricia Kovatch, Robert Freeman, Joy Jiang, Pushkala Jayaraman, Zahi Fayad, et al. Implications of the use of artificial intelligence predictive models in health care settings: a simulation study. Annals of internal medicine, 176(10):1358– 1369, 2023. Aad Van Der Vaart and Jon A Wellner. Preservation theorems for glivenko-cantelli and uniform glivenkocantelli classes. In High dimensional probability II, pages 115–133. Springer, 2000. Aad W Van Der Vaart and Jon A Wellner. Weak convergence. In Weak convergence and empirical processes: with applications to statistics, pages 16–28. Springer, 1996. A.W. van der Vaart. Asymptotic Statistics. Asymptotic Statistics. Cambridge University Press, 2000. ISBN 9780521784504. Veniamin Veselovsky, Manoel Horta Ribeiro, and Robert West. Artificial artificial artificial intelligence: Crowd workers widely use large language models for text production tasks. arXiv preprint arXiv:2306.07899, 2023. Hung Ahn Vu, Galen Reeves, and Emily Wenger. What happens when generative ai models train recursively on each others’ generated outputs? arXiv preprint arXiv:2505.21677, 2025. Xiaohan Xu, Ming Li, Chongyang Tao, Tao Shen, Reynold Cheng, Jinyang Li, Can Xu, and Dacheng Tao. A survey on knowledge distillation of large language models. arXiv preprint arXiv:2402.13116, 2024. Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang. Diffusion models: A comprehensive survey of methods and applications. ACM computing surveys, 56(4):1–39, 2023.

Proofs of the main results

We present in this section proofs of our main results.

A.1

Proof of Theorem 3.1

Iterating Eq. (14), we see that βbt − βb0 = vt + Tt vt−1 + · · · + Tt · · · T3 v2 + Tt · · · T2 v1 =

t X

Ωt,s+1 vs ,

(25)

s=1

where for t ≥ s, we define Ωt,s = Tt · · · Ts ∈ RdK×dK , with the convention that Ωt,t+1 = IdK . The matrix Ωt,s can be viewed as a K × K block matrix, and we denote its (i, j)-th block by Ωt,s,µi ,µj ∈ Rd×d . The next lemma is a direct consequence of Assumption 3.2. Lemma A.1. Under Assumption 3.2, for all t ∈ N+ , there exist Nt ∈ N+ , a set of vectors {ζt,1 , ζt,2 , · · · , ζt,Nt } and a set of matrices {Zt,1 , Zt,2 , · · · , Zt,Nt }, such that the following holds: 1. For all s ∈ [t] and (ν, µ) ∈ E, εs,ν→µ and Xs,ν→µ admit the following decompositions: ⊤ ⊤ ⊤ εs,ν→µ = (ε⊤ s,ν→µ,1 , εs,ν→µ,2 , · · · , εs,ν→µ,hs,ν→µ ) ,

⊤ ⊤ ⊤ Xs,ν→µ = (Xs,ν→µ,1 , Xs,ν→µ,2 , · · · , Xs,ν→µ,h )⊤ , s,ν→µ

where hs,ν→µ ∈ N+ , and the number of rows of the matrix Xs,ν→µ,h is equal to the length of the vector εs,ν→µ,h for all h ∈ [hs,ν→µ ]. Additionally, there exists i ∈ [Nt ], such that εs,ν→µ,h ≡ ζt,i and Xs,ν→µ,h ≡ Zt,i . On the other hand, for any t ∈ N+ and i ∈ [Nt ], there exists (s, ν → µ, h) such that εs,ν→µ,h ≡ ζt,i and Xs,ν→µ,h ≡ Zt,i . 2. There exists C > 0 that does not depend on t, such that Nt ≤ Ct, and (recall ρ2 is from Assumption 3.4) ⊤ ⊤ ⊤ Cov[ζt ] = Cov (ζt,1 , ζt,2 , · · · , ζt,N )⊤ ⪰ ρ2 I. t Proof of Lemma A.1. We prove Lemma A.1 in Appendix B.1. For µ ∈ Ml , t ∈ N+ and i ∈ [Nt ], we define Pt,i = {(s, ν1 → ν2 , h) : εs,ν1 →ν2 ,h ≡ ζt,i }, X ⊤ ⊤ Ωt,s+1,µ,ν Ts,ν,m (Xs,m→ν Xs,m→ν )−1 Zt,i Zt,i . At,i,µ =

(26) (27)

(s,m→ν,h)∈Pt,i

The set Pt,i consists of indices for all pairs (εs,ν1 →ν2 ,h , Xs,ν1 →ν2 ,h ) such that (εs,ν1 →ν2 ,h , Xs,ν1 →ν2 ,h ) ≡ (ζt,i , Zt,i ). By the first point of Lemma A.1, for all t ∈ N+ and i ∈ [Nt ] we have Pt,i ̸= ∅. The next lemma establishes additional properties of Pt,i . Lemma A.2. Under Assumption 3.2, there does not exist s1 ̸= s2 with s1 , s2 ∈ [t], h1 , h2 ∈ N+ and ν1 , ν2 , µ1 , µ2 ∈ M, such that (s1 , ν1 → µ1 , h1 ), (s2 , ν2 → µ2 , h2 ) ∈ Pt,i for some t ∈ N+ and i ∈ [Nt ]. Consequently, the index s ∈ [t] associated with each Pt,i is unique, and we denote it by st,i . Proof of Lemma A.2. By Assumption 3.1 and 3.4, for any s1 ̸= s2 , h1 , h2 ∈ N+ and ν1 , ν2 , µ1 , µ2 ∈ M, the noise vectors εs1 ,ν1 →µ1 ,h1 and εs2 ,ν2 →µ2 ,h2 are independent and both nondegenerate. Therefore, they cannot be equal to the same ζt,i , completing the proof. With At,i,µ , we find a more convenient expression for βbt,µ , as stated in the next lemma. Lemma A.3. Under Assumptions 3.2 and 3.6, for all µ ∈ Ml and t ∈ N+ , it holds that βbt,µ = βb0,µ +

Nt X

⊤ ⊤ At,i,µ (Zt,i Zt,i )† Zt,i ζt,i .

i=1

Proof of Lemma A.3. We prove Lemma A.3 in Appendix B.2. Note that Cov[ζt ] ⪰ ρ2 I by the second point of Lemma A.1. Hence, we obtain the following lower bound: Nt Nt h h X X i i ⊤ ⊤ ⊤ † E ∥βb0,µ − βbt,µ ∥22 ≥ ρ2 E Tr At,i,µ (Zt,i = ρ . Zt,i )† A⊤ E Tr A A (Z Z ) 2 t,i,µ t,i t,i,µ t,i,µ t,i i=1

(28)

i=1

For any t ∈ N+ and i ∈ [Nt ], as Pt,i ̸= ∅ by the first point of Lemma A.1, there exists (st,i , νt,i,1 → νt,i,2 , ht,i ) ∈ Pt,i . Therefore, for any t ∈ N+ , i ∈ [Nt ] and µ ∈ Ml , h i ⊤ † E Tr A⊤ A (Z Z ) t,i,µ t,i t,i,µ t,i h (i) i ⊤ −1 ≥ E Tr A⊤ t,i,µ At,i,µ (Xst,i ,νt,i,1 Xst,i ,νt,i,1 ) h (ii) i −1 ⊤ ⊤ (29) ≥ n−1 κ E Tr A A 1{X X /n ⪯ κI } t,i,µ s ,ν s ,ν d st,i ,νt,i,1 t,i,µ st,i ,νt,i,1 t,i t,i,1 t,i t,i,1 h i −1 ⊤ ≥ n−1 E Tr A⊤ st,i ,max κ t,i,µ At,i,µ 1{Xst,i ,νt,i,1 Xst,i ,νt,i,1 /nst,i ,νt,i,1 ⪯ κId } h i −1 ⊤ ⊤ ≥ n−1 κ E Tr A A min 1{X X /n ⪯ κI } , t,i,µ s ,v s ,v d st,i ,max t,i,µ st,i ,v t,i t,i v∈M

where we recall that κ is from Assumption 3.5, and ns,max is from Assumption 3.3. In the above equations, (i) is because Xs⊤t,i ,νt,i,1 Xst,i ,νt,i,1 is invertible by Assumption 3.6, and (ii) is because ⊤ −1 Tr A⊤ t,i,µ At,i,µ (Xst,i ,νt,i,1 Xst,i ,νt,i,1 ) ⊤ −1 = Tr A⊤ 1{Xs⊤t,i ,νt,i,1 Xst,i ,νt,i,1 /nst,i ,νt,i,1 ⪯ κId } t,i,µ At,i,µ (Xst,i ,νt,i,1 Xst,i ,νt,i,1 ) ⊤ −1 + Tr A⊤ 1 − 1{Xs⊤t,i ,νt,i,1 Xst,i ,νt,i,1 /nst,i ,νt,i,1 ⪯ κId } t,i,µ At,i,µ (Xst,i ,νt,i,1 Xst,i ,νt,i,1 ) ⊤ ⊤ −1 ≥ n−1 1{Xs⊤t,i ,νt,i,1 Xst,i ,νt,i,1 /nst,i ,νt,i,1 ⪯ κId } st,i ,νt,i,1 Tr At,i,µ At,i,µ (Xst,i ,νt,i,1 Xst,i ,νt,i,1 /nst,i ,νt,i,1 ) −1 ⊤ ≥ n−1 Tr A⊤ st,i ,νt,i,1 κ t,i,µ At,i,µ 1{Xst,i ,νt,i,1 Xst,i ,νt,i,1 /nst,i ,νt,i,1 ⪯ κId } . The next lemma analyzes the sum

PNt

i=1 At,i,µ .

Lemma A.4. Under Assumptions 3.2 and 3.6, for all t ∈ N+ we have Nt X

At,i,µ =

t X X X s=1 ν∈Ml m∈M

i=1

Ωt,s+1,µ,ν Ts,ν,m .

Proof of Lemma A.4. We prove Lemma A.4 in Appendix B.3. We use v1 → v2 → · · · → vx to represent a path connecting nodes v1 , v2 , · · · , vx ∈ M, with (vi , vi+1 ) ∈ E for i = 1, 2, · · · , x − 1, and we define x − 1 as the length of this path. For ℓ ∈ N+ , we denote by Lℓ the collection of all paths of length ℓ in G. Throughout this paper, we make the convention that L0 = ∅. Then by Lemma A.4, Nt t hX i hX i X X E At,i,µ =E Ωt,s+1,µ,ν Ts,ν,m s=1 ν∈Ml m∈M

i=1

t X

x=1 v1 →v2 →···→vx →µ∈Lx t (i) X

x=1 v1 →v2 →···→vx →µ∈Lx

h i E Tt,µ,vx Tt−1,vx ,vx−1 · · · Tt−x+1,v2 ,v1 E[Tt,µ,vx ]E[Tt−1,vx ,vx−1 ] · · · E[Tt−x+1,v2 ,v1 ],

(30)

where (i) is by Assumption 3.1. For x ∈ N+ , ℓ ∈ {0} ∪ [x], and v ∈ M∞ l , we denote by Lx,ℓ,v the collection ∞ of paths of length x whose last visit to M∞ occurs at model v ∈ M l l at length x − ℓ: n o Lx,ℓ,v = v1 → v2 → · · · → vx → vx+1 : vx−ℓ+1 = v ∈ M∞ / M∞ . (31) l , vx−ℓ+2 , vx−ℓ+3 , · · · , vx+1 ∈ l For x ∈ N+ , we denote by L∗x the collection of paths of length x that do not hit models in M∞ l : n o L∗x = v1 → v2 → · · · → vx → vx+1 : vi ∈ / M∞ l for all i ∈ [x + 1] .

(32)

Similarly, we make the convention that L∗0 = L0,ℓ,v = ∅. With definitions (31) and (32), we then reformulate the last line of Eq. (30) as Nt t x hX i X X X E At,i,µ =

x=1 v∈M∞ l ℓ=0 v1 →v2 →···→vx →µ∈Lx,ℓ,v

i=1

t X

x=1 v1 →v2 →···→vx →µ∈L∗ x

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 (33)

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 .

We next analyze the above two sums on the right hand side of Eq. (33), starting from the second sum. Lemma A.5. Under Assumptions 3.1-3.3, for any µ ∈ Ml \M∞ l , t ∈ N+ and x ∈ [t], it holds that X v1 →v2 →···→vx →µ∈L∗ x

≤

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1

X v1 →v2 →···→vx →µ∈L∗ x

E Tt,µ,vx op E Tt−1,vx ,vx−1 op · · · E Tt−x+1,v2 ,v1 op

⌈x/K⌉ δ + 2γ1 K δ + 2γ1 K − αK . 1+ ≤ 1+ α α By Assumption 3.5, we further know that the last line above is no larger than 2c(α, K)⌈x/K⌉ . On the other hand, for any µ ∈ M∞ l , straightforwardly we have X E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 = 0d×d . v1 →v2 →···→vx →µ∈L∗ x

Proof of Lemma A.5. We prove Lemma A.5 in Appendix B.4. We then switch to analyze the first sum on the right hand side of Eq. (33). c Lemma A.6. Recall that Mcl and Mnc l are defined in Eq. (4). Under Assumptions 3.1-3.3, if µ ∈ Ml , then for t ≥ K, it holds that t x X X X

x=1 v∈M∞ l ℓ=0 v1 →v2 →···→vx →µ∈Lx,ℓ,v

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1

≥ t c2 − αK (K − 1),

where we recall c1 , c2 are from Assumption 3.5. On the other hand, if µ ∈ Mnc l , then straightforwardly, t x X X X

x=1 v∈M∞ l ℓ=0 v1 →v2 →···→vx →µ∈Lx,ℓ,v

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 = 0d×d .

Proof of Lemma A.6. We prove Lemma A.6 in Appendix B.5.

Substituting the bounds from Lemmas A.5 and A.6 into Eq. (33), we conclude that for t ≥ K and µ ∈ Mcl , Nt hX i E At,i,µ i=1

≥ t c2 − α (K − 1) − ≥ t c2 − αK (K − 1) −

t X x=1 t X

⌈x/K⌉ δ + 2γ1 K δ + 2γ1 K K 1+ 1+ −α α α

2c(α, K)⌈x/K⌉

(34)

x=1

≥ t c2 − αK (K − 1) − 2K

c(α, K) 1 − c(α, K)

≥ t c2 − αK (K − 1) − 4Kα−K ,

where the last inequality follows from taking c(α, K) = 1 − αK /2, as pointed out in Assumption 3.5. We next upper bound the operator norm of Nt hX i E At,i,µ max 1 − 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } .

(35)

v∈M

i=1

The following lemma provides a more tractable characterization of this sum. Lemma A.7. Under Assumptions 3.2 and 3.6, for all t ∈ N+ we have Nt X i=1

At,i,µ max 1 − 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } v∈M

t X

X X

s=1 ν∈Ml m∈M

Ωt,s+1,µ,ν Ts,ν,m max 1 − 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } . v∈M

Proof of Lemma A.7 . The proof is completely the same as that of Lemma A.4, by noting that st,i = s in the summation. We then leverage Lemma A.7 to upper bound the operator norm of (35). Lemma A.8. Under Assumptions 3.1-3.3, for any µ ∈ Ml and t ∈ N+ , it holds that Nt hX i E At,i,µ max 1 − 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } v∈M

i=1

≤

9 t K 3 γ2 . αK+1

Proof of Lemma A.8. We prove Lemma A.8 in Appendix B.6. We then combine Lemma A.8 and Eq. (34) to obtain the following operator norm lower bound for µ ∈ Mcl : Nt i hX E At,i,µ min 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } v∈M

i=1

Nt hX i ≥ E At,i,µ i=1

Nt hX i − E At,i,µ max 1 − 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } i=1

≥ t c2 − αK (K − 1) − 4Kα−K −

v∈M

9 t K 3 γ2 ≥ αK+1

c2 −

9K 3 γ2 αK+1

Leveraging Eqs. (28), (29) and the above lower bound, we see that E ∥βb0,µ − βbt,µ ∥2 31

t−

5K . αK

≥ ρ2 κ−1

Nt X i=1

≥ ρ2 κ−1 −1

= ρ2 κ

h i ⊤ ⊤ n−1 E Tr A A min 1{X X /n ⪯ κI } t,i,µ s ,v s ,v d st,i ,max t,i,µ st,i ,v t,i t,i v∈M

sup s∈{0}∪[t]

≥ ρ2 κ−1 ≥ ρ2 κ−1 ≥ ρ2 κ−1

ns,max

sup

Nt h i −1 X ⊤ E Tr A⊤ t,i,µ At,i,µ min 1{Xst,i ,v Xst,i ,v /nst,i ,v ⪯ κId }

Nt h i −1 X 2 ns,max E At,i,µ min 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } F

ns,max

−1

Nt hX i 2 Nt−1 E At,i,µ min 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } v∈M

i=1

ns,max

−1

s∈{0}∪[t]

sup

v∈M

i=1

s∈{0}∪[t]

sup

v∈M

i=1

Nt hX i 2 Nt−1 E At,i,µ min 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } v∈M

i=1

ns,max

−1

Nt−1

s∈{0}∪[t]

c2 −

5K 9K 3 γ2 t− K αK+1 α

2 .

Recall Nt = O(t) by Lemma A.1, and sups∈{0}∪[t] ns,max /t → 0 as t → ∞ by Assumption 3.3. This implies that the last line above tends to infinity as t → ∞ under the current set of assumptions. By the temporal independence assumption (Assumption 3.1), we have lim inf E[∥β∗ − βbt,µ ∥22 ] = lim inf E ∥βb0,µ − βbt,µ ∥22 + E ∥β∗ − βb0,µ ∥22 → ∞ as t → ∞. t→∞

t→∞

(36)

Under Assumptions 3.1 – 3.6, we have h h X −1 i i ∗ ⊤ lim sup E[∥β∗ − βbt,µ ∥22 ] ≤ lim sup ρ1 E Tr Xt,ν→µ Xt,ν→µ < ∞. t→∞

t→∞

(37)

ν∈Nµ,in

Eqs. (36) and (37) together complete the proof of the theorem.

A.2

Proof of Theorem 3.2

Recall that for t ≥ s, we define Ωt,s = Tt Tt−1 · · · Ts ∈ RdK×dK , with the convention that Ωt,t+1 = IdK . The matrix Ωt,s has a K × K block structure, where the (i, j)-th block is denoted by Ωt,s,µi ,µj ∈ Rd×d . According nc c to the definition of Mnc l in Eq. (4), for every µ ∈ Ml , there is no path originating in Ml and terminating , at µ. Invoking this observation together with Eq. (25), we conclude that for all µ ∈ Mnc l βbt,µ − βb0,µ =

t X

Ωt,s+1,µ,ν vs,ν ,

(38)

s=1 ν∈M\Mcl

where we recall that vs,ν is defined in Eq. (12). Under Assumption 3.1, vs is independent of Ωt,s+1 . Combining this observation with Eq. (38) yields t h h X 2i E βbt,µ − βb0,µ = E 2

 t X ≤ 

 t X = 

(i)

s=1 ν∈M\Mcl

Ωt,s+1,µ,ν vs,ν

2i 2

2 h i1/2 ⊤  E vs,ν Ω⊤ t,s+1,µ,ν Ωt,s+1,µ,ν vs,ν

s=1 ν∈M\Mcl

(ii)

2 h i1/2 ⊤  , E vs,ν E[ Ω⊤ t,s+1,µ,ν Ωt,s+1,µ,ν ] vs,ν

s=1 ν∈M\Mcl

(39)

where (i) follows from the Cauchy–Schwarz inequality, and (ii) follows from the fact that vs is independent of Ωt,s+1 . The remainder of this section is devoted to establishing an upper bound for h i1/2 X ⊤ E vs,ν E[ Ω⊤ Ω ] v Et,s,µ = s,ν t,s+1,µ,ν t,s+1,µ,ν ν∈M\Mcl

for all t ∈ N+ , s ∈ [t] and µ ∈ Mnc l . Next, we employ an induction argument to show that this quantity decays exponentially in t. Specifically, we prove that Et+t0 ,s,µ ≤ θEt,s,µ for suitably chosen constants θ ∈ (0, 1) and t0 ∈ N+ . Recall that Ωt,s+1 = Tt Tt−1 · · · Ts+1 . Therefore, for any t0 ∈ N+ , we have Ωt+t0 ,s+1 = Ωt+t0 ,t+1 Ωt,s+1 . c Therefore, for all µ ∈ Mnc l and ν ∈ M\Ml , X X Ωt+t0 ,t+1,µ,m Ωt,s+1,m,ν vs,ν . (40) Ωt+t0 ,s+1,µ,ν vs,ν = Ωt+t0 ,t+1,µ,m Ωt,s+1,m,ν vs,ν = m∈Mnc l

m∈M

Substituting Eq. (40) into the definition of Et,s,µ , we obtain i 1/2 X h X ⊤ ⊤ Et+t0 ,s,µ = E vs,ν Ω⊤ Ω Ω Ω v t,s+1,m1 ,ν t+t0 ,t+1,µ,m1 t+t0 ,t+1,µ,m2 t,s+1,m2 ,ν s,ν m1 ,m2 ∈Mnc l

ν∈M\Mcl

X ν∈M\Mcl

h i 1/2 ⊤ ⊤ E vs,ν Ω⊤ , t,s+1,m1 ,ν E Ωt+t0 ,t+1,µ,m1 Ωt+t0 ,t+1,µ,m2 Ωt,s+1,m2 ,ν vs,ν

X m1 ,m2 ∈Mnc l

where the last equality above follows from Assumption 3.1 and the observation that (Ωt,s+1 , vs ) is a function of (Xi , εi )s≤i≤t , and Ωt+t0 ,t+1 is a function of (Xi )t+1≤i≤t+t0 . From the last line above, we further conclude that X X ⊤ ⊤ 1/2 Et+t0 ,s,µ ≤ E ∥vs,ν Ωt,s+1,m1 ,ν ∥2 ∥Ωt,s+1,m2 ,ν vs,ν ∥2 E Ω⊤ Ω t+t ,t+1,µ,m t+t0 ,t+1,µ,m1 0 2 op ν∈M\Mcl m1 ,m2 ∈Mnc l

(i)

≤ =

h i1/2 1/2 ⊤ E vs,ν Ω⊤ E Ω⊤ t,s+1,m,ν Ωt,s+1,m,ν vs,ν t+t0 ,t+1,µ,m Ωt+t0 ,t+1,µ,m op

ν∈M\Mcl m∈Mnc l

1/2 Et,s,m E Ω⊤ , t+t0 ,t+1,µ,m Ωt+t0 ,t+1,µ,m

(41)

m∈Mnc l

where (i) is obtained by applying the matrix Cauchy-Schwarz inequality as follows: ⊤ ⊤ ⊤ ⊤ 1/2 ⊤ ⊤ 1/2 E ∥vs,ν Ωt,s+1,m1 ,ν ∥2 ∥Ωt,s+1,m2 ,ν vs,ν ∥2 ≤ E vs,ν Ωt,s+1,m1 ,ν Ωt,s+1,m1 ,ν vs,ν E vs,ν Ωt,s+1,m2 ,ν Ωt,s+1,m2 ,ν vs,ν , ⊤ ⊤ 1/2 ⊤ 1/2 E Ωt+t0 ,t+1,µ,m1 Ωt+t0 ,t+1,µ,m2 op ≤ E Ωt+t0 ,t+1,µ,m1 Ωt+t0 ,t+1,µ,m1 op E Ωt+t0 ,t+1,µ,m2 Ωt+t0 ,t+1,µ,m2 op . Using Eq. (41), we have ( sup Et+t0 ,s,µ ≤ sup Et,s,µ × sup

µ∈Mnc l

) ⊤ 1/2 E Ωt+t0 ,t+1,µ,m Ωt+t0 ,t+1,µ,m .

(42)

m∈Mnc l

To establish that supµ∈Mnc Et,s,µ decays exponentially in t, we leverage Assumption 3.5, which characterizes l ∗ ∗ the concentration behavior of Tt . Recall Tt∗ is defined in Eq. (10). Define Ω∗t,s+1 = Tt∗ Tt−1 · · · Ts+1 for ∗ 0 ≤ s ≤ t, with the convention that Ωt,t+1 = IdK . The following lemma shows that, under Assumption 3.5, Ω∗t,s+1 provides a close approximation to Ωt,s+1 with high probability. Lemma A.9. Recall Nmax = supµ∈M |Nµ,in |. Under Assumption 3.5, for any 0 ≤ s ≤ t and ν1 , ν2 ∈ M, Ω∗t,s+1,ν1 ,ν2 − Ωt,s+1,ν1 ,ν2 op ≤ γ1 with probability at least 1 − δ

Pt−s

j j=1 Nmax .

t−s X j=1

j Nmax

(43)

Proof of Lemma A.9. We prove Lemma A.9 in Appendix B.7. For all t ≥ s ≥ 0 and µ, ν ∈ M, define the event: ( Ω∗t,s+1,µ,ν − Ωt,s+1,µ,ν op ≤ γ1

Et,s+1,µ,ν =

t−s X

) j Nmax .

j=1

By Lemma A.9, we have P (Et,s+1,µ,ν ) ≥ 1 − δ

t−s X

j Nmax .

j=1

On the event Et,s+1,µ,ν , we have ∗⊤ ∗ Ω⊤ t,s+1,µ,ν Ωt,s+1,µ,ν − Ωt,s+1,µ,ν Ωt,s+1,µ,ν op (i)

≤ 2 Ωt,s+1,µ,ν − Ω∗t,s+1,µ,ν op + Ωt,s+1,µ,ν − Ω∗t,s+1,µ,ν op

≤ 2γ1

t−s X

j Nmax + γ12

j=1

2 j Nmax ,

j=1

where (i) follows from the fact that ∥Ω∗t,s+1,µ,ν ∥op ≤ 1. By induction, it follows that t−s ∥Ωt,s+1,µ,ν ∥op ≤ Nmax .

Combining the above estimates, we get E Ω⊤ t,s+1,µ,ν Ωt,s+1,µ,ν op h i ⊤ c ≤ E Ωt,s+1,µ,ν Ωt,s+1,µ,ν 1Et,s+1,µ,ν op + E Ω⊤ t,s+1,µ,ν Ωt,s+1,µ,ν 1Et,s+1,µ,ν ≤ Ω∗t,s+1,µ,ν Ω∗t,s+1,µ,ν op + 2γ1 ≤ Ω∗t,s+1,µ,ν Ω∗t,s+1,µ,ν

+ 2γ1

t−s X

j Nmax + γ12

j=1 t−s X

t−s X

j Nmax

2(t−s) c + Nmax P(Et,s+1,µ,ν )

j=1 j Nmax + γ12

t−s X

+δ

j+2(t−s) Nmax ,

j=1

t−s X

which implies that 1/2 E Ω⊤ ≤ t,s+1,µ,ν Ωt,s+1,µ,ν op

Ω∗t,s+1,µ,ν

t−s √ t−s X p j + γ1 + 2γ1 + δNmax Nmax .

(44)

j=1

Leveraging Eq. (44), we conclude that for all µ ∈ Mnc l : X

1/2 E Ω⊤ t+t0 ,t+1,µ,m Ωt+t0 ,t+1,µ,m op

m∈Mnc l

≤

Ω∗t+t0 ,t+1,µ,m op + K

γ1 +

2γ1 +

m∈Mnc l

By Eq. (45), to upper bound

√

t0 δNmax

t0 X

(45) j Nmax .

j=1 1/2 nc ∥E[Ω⊤ t+t0 ,t+1,µ,m Ωt+t0 ,t+1,µ,m ]∥op for µ ∈ Ml , it suffices to bound m∈Mnc l

Ω∗t+t0 ,t+1,µ,m op

m∈Mnc l

for µ ∈ Mnc Ω∗ and the Markov chain P defined in l . To this end, we exploit the connection between P Eq. (15). By definition, Pt is row-stochastic matrix, since ν∈M Pt,µ,ν = 1 for all µ ∈ M. Consequently, the sequence {Pt }t≥0 defines a (possibly time-inhomogeneous) Markov chain. Recall for t ≥ s ≥ 0 we have Jt,s+1 = Pt Pt−1 · · · Ps+1 , and Jt,t+1 = IK . For any µ ∈ Mnc l , we have X X Jt+t0 ,t+1,µ,m ≤ 1. (46) Ω∗t+t0 ,t+1,µ,m op = m∈Mnc l

m∈Mnc l

From Eq. (16), we know that for all t0 ≥ T0 , we have X Jt+t0 ,t+1,ν1 ,ν2 ≤ ω ∈ (0, 1). sup t∈N+ , ν1 ∈Mnc l

ν2 ∈Mnc l

Combining the above upper bound with Eqs. (45) and (46), we conclude that under Assumption 3.5, for all 2T0 − 1 ≥ t0 ≥ T0 ,    X ⊤ 1/2  E Ωt+t0 ,t+1,µ,m Ωt+t0 ,t+1,µ,m op sup   µ∈Mnc nc l m∈Ml

t0 √ t X p j 0 ≤ ω + K γ1 + 2γ1 + δNmax Nmax ≤ α0 ∈ (0, 1). j=1

As a consequence of the above upper bound and Eq. (42), we know that for all T0 ≤ t0 ≤ 2T0 − 1, we have sup Et+t0 ,s,µ ≤ α0 × sup Et,s,µ ,

µ∈Mnc l

which implies that for all t ≥ s + T0 : ⌊(t−s)/T0 ⌋

sup Et,s,µ ≤ α0

µ∈Mnc l

⌊(t−s)/T0 ⌋

× sup Es,s,µ = α0 µ∈Ml

X ν∈M\Mcl

⊤ 1/2 E vs,ν vs,ν (47)

⊤ 1/2 ⌊(t−s)/T0 ⌋ ≤ Kα0 · sup E vs,ν vs,ν , ν∈M

where we recall that vs,ν is from Eq. (12). On the other hand, for s ≤ t ≤ s + T0 − 1, there exists C ′ > 0 that depends only on (G, α), such that ⊤ 1/2 sup Et,s,µ ≤ C ′ · sup E vs,ν vs,ν . (48) µ∈Mnc l

ν∈M

Combining Eqs. (39), (47) and (48), we conclude that for all µ ∈ Mnc l , there exists CG > 0 that depends only on (G, α), such that h i h i h i E ∥βbt,µ − β∗ ∥22 ≤2E ∥βbt,µ − βb0,µ ∥22 + 2E ∥βb0,µ − β∗ ∥22 ≤

t X

Et,s,µ

⊤ + 2E v0,µ v0,µ

s=1

≤CG ·

sup

⊤ sup E vs,ν vs,ν .

s∈{0}∪[t] ν∈M

⊤ T By assumption supt≥1 (sups∈{0}∪[t] supν∈M E vs,ν vs,ν )/(E[vt,µ vt,µ ]) < ∞, hence completing the proof.

A.3

Proof of Theorem 4.1

We prove our conclusion by induction. For t = 0, the consistency and asymptotic normality of βb0 follows from standard arguments in empirical process theory and our assumptions. See Appendix C.1 for a more √ d rigorous treatment. Furthermore, we have n0 (βb0 − βb−1 ) → N(0, Σ0 ), where Σ0 is defined in Eq. (20). 35

Now we proceed with the induction step. Assume that the conclusion already holds for rounds 0, 1, · · · , t− 1, we first show that the βbt,µ ’s are all consistent estimators of β∗ . For µ ∈ Mu , this is automatically true since βbt,µ = βbt−1,µ . For µ ∈ Ml , note that by definition: 1

βbt,µ = arg min P

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

β∈B

where we use the shorthand

L (β, Zt,ν→µ ) ,

nt,ν→µ

L (β, Zt,ν→µ ) =

L(β, zt,ν→µ,i ).

i=1

Recall that zt,ν→µ,i = φ(βbt−1,ν , εt,ν→µ,i ) for εt,ν→µ,i ∼i.i.d. µε . To establish consistency of βbt,µ , we invoke Theorem 3.2.3 (ii) in Van Der Vaart and Wellner (1996). Define 1

ct,µ (β) = P M

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

ft,µ (β) = P M

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

L (β, Zt,ν→µ ) , h i E L (β, Zt,ν→µ ) βbt−1,ν ,

M (β) = E [L(β, φ(β∗ , ε))] . By Assumption 4.3, M (β) is continuous and uniquely minimized at β∗ . It then suffices to show that (i) P ct,µ (β) − M (β)| → supβ∈Ω |M 0 for any compact set Ω; (ii) The sequence βbt,µ is uniformly tight. We first prove claim (i). For any δ > 0, define the event o n βbt−1,ν − β∗ ≤ δ, ∀ν ∈ M . Ξt−1,δ = 2

Then, our induction hypothesis implies that P(Ξt−1,δ ) → 1 as the sample sizes tend to infinity. Conditioned on Ξt−1,δ , we have ct,µ (β) − M ft,µ (β) sup M

β∈Ω

≤ P

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

sup β∈Ω bt−1,ν −β∗ ∥2 ≤δ ∥β

h i L (β, Zt,ν→µ ) − E L (β, Zt,ν→µ ) βbt−1,ν .

By point 1 of Assumption 4.5, the function class ε 7→ L (β1 , φ(β2 , ε)) β1 ∈ Ω, ∥β2 − β∗ ∥2 ≤ δ is µε -Glivenko-Cantelli. Therefore, we can apply the uniform law of large numbers to deduce that h i 1 P sup L (β, Zt,ν→µ ) − E L (β, Zt,ν→µ ) βbt−1,ν → 0 nt,ν→µ β∈Ω bt−1,ν −β∗ ∥2 ≤δ ∥β

for all ν ∈ Nµ,in . This implies that

P ct,µ (β) − M ft,µ (β) → sup M 0

β∈Ω

ft,µ (β) − as the sample sizes tend to infinity. To finish the proof of (i), it remains to show that supβ∈Ω |M P

M (β)| → 0. Again, conditioning on Ξt−1,δ , we obtain that ft,µ (β) − M (β)| ≤ sup |M

β∈Ω

sup β∈Ω ∥β ′ −β∗ ∥2 ≤δ

|E[L(β, φ(β ′ , ε))] − E[L(β, φ(β∗ , ε))]| .

By Assumption 4.3, we know that the mapping (β1 , β2 ) 7→ E[L(β1 , φ(β2 , ε))] is continuous, and therefore uniformly continuous on any compact set. Consequently, for any δ ′ > 0, we may choose δ sufficiently small such that sup |E[L(β, φ(β ′ , ε))] − E[L(β, φ(β∗ , ε))]| ≤ δ ′ , β∈Ω ∥β ′ −β∗ ∥2 ≤δ

which further implies that ! ft,µ (β) − M (β)| ≥ δ ′ sup |M

β∈Ω

→ 0.

P ft,µ (β) − M (β)| → Since δ ′ is arbitrary, we deduce that supβ∈Ω |M 0, which completes the proof of claim (i). As for claim (ii), note that Assumption 4.4 implies that either L is convex or B is compact, which ensures the uniform tightness of βbt,µ . It follows that βbt,µ is a consistent estimator of β∗ . We next establish asymptotic normality of βbt and derive the recursive formula (22) for computing its √ asymptotic variance. The first step is to obtain an asymptotic equivalent expression for nt (βbt − βb−1 ). Using Assumption 4.4 and the consistency of βbt,µ , we know that βbt,µ ∈ int B with probability 1 − o(1). Therefore, βbt,µ must satisfy the following first-order condition (∇β1 denotes the gradient with respect to the first argument): X 1 P ∇β1 L βbt,µ , Zt,ν→µ = 0, ν∈Nµ,in nt,ν→µ ν∈Nµ,in

which further implies that 1

−P

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

= P

∇β1 L (β∗ , Zt,ν→µ )

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

 = P

∇β1 L βbt,µ , Zt,ν→µ − ∇β1 L (β∗ , Zt,ν→µ )

X Z 1

ν∈Nµ,in nt,ν→µ ν∈Nµ,in



∇2β1 L αβbt,µ + (1 − α)β∗ , Zt,ν→µ dα βbt,µ − β∗ .

√ d As our goal is to show nt (βbt − βb−1 ) → N(0, Σt ), we can assume without loss of generality that ∥βbt,µ −β∗ ∥2 ≤ 1, which happens with high probability by consistency. Using point 1 of Assumption 4.5 and the uniform law of large numbers, it follows that X Z 1 1 P ∇2β1 L αβbt,µ + (1 − α)β∗ , Zt,ν→µ dα ν∈Nµ,in nt,ν→µ ν∈Nµ,in 0 Z 1 h i X 1 P P nt,ν→µ E ∇2β1 L αβbt,µ + (1 − α)β∗ , φ(βbt−1,ν , ε) dα → 0 − n 0 ν∈Nµ,in t,ν→µ ν∈Nµ,in

as the sample sizes tend to infinity, where the expectation above is taken with respect to ε ∼ µε . Since βbt,µ and the βbt−1,ν ’s are consistent estimators of β∗ , we deduce from Assumption 4.6 that for all ν ∈ Nµ,in , Z 1 h i P E ∇2β1 L αβbt,µ + (1 − α)β∗ , φ(βbt−1,ν , ε) dα → E ∇2β1 L (β∗ , φ(β∗ , ε)) = E ∇2β L(β, φ(β∗ , ε)) β=β∗

as the sample sizes tend to infinity. Further, Assumption 4.3 implies that E[∇2β1 L(β∗ , φ(β∗ , ε))] is invertible, we thus obtain that √ b nt βt,µ − β∗ √ X −1 nt P = − E ∇2β1 L (β∗ , φ(β∗ , ε)) ∇β1 L (β∗ , Zt,ν→µ ) + oP (1). ν∈Nµ,in nt,ν→µ ν∈Nµ,in

Now for each ν ∈ Nµ,in , we analyze the term ∇β1 L (β∗ , Zt,ν→µ ) in more details. Note that √ =√ =√

1 nt,ν→µ

∇β1 L (β∗ , Zt,ν→µ ) = √ X

nt,ν→µ

i=1

∇β1 L(β∗ , zt,ν→µ,i )

∇β1 L(β∗ , φ(βbt−1,ν , εt,ν→µ,i ))

nt,ν→µ

i=1

nt,ν→µ

i=1

nt,ν→µ

√

nt,ν→µ

h i ∇β1 L β∗ , φ(βbt−1,ν , εt,ν→µ,i ) − E ∇β1 L β∗ , φ(βbt−1,ν , ε)

h i nt,ν→µ E ∇β1 L β∗ , φ(βbt−1,ν , ε) := (I) + (II).

For (I), we note that point 2 of Assumption 4.5 implies uniform CLT for βbt−1,ν in a neighborhood of β∗ . Since the limiting Gaussian process is continuous and βbt−1,ν − β∗ = oP (1), we know that (I) = √ =√

1 nt,ν→µ 1 nt,ν→µ

nt,ν→µ

X i=1

(∇β1 L (β∗ , φ(β∗ , εt,ν→µ,i )) − E [∇β1 L (β∗ , φ(β∗ , ε))]) + oP (1)

nt,ν→µ

X i=1

∇β1 L (β∗ , φ(β∗ , εt,ν→µ,i )) + oP (1),

where the last equality follows from Assumption 4.3: E[∇β1 L(β∗ , φ(β∗ , ε))] = 0. For (II), we invoke Assumption 4.6 together with the induction hypothesis that βbt−1,ν is asymptotically normal to deduce that ⊤√ (II) = ∇β2 E [∇β1 L (β∗ , φ(β∗ , ε))] nt,ν→µ βbt−1,ν − β∗ + oP (1) (i) ⊤√ = − ∇β1 E [∇β1 L (β∗ , φ(β∗ , ε))] nt,ν→µ βbt−1,ν − β∗ + oP (1) √ (ii) = − E ∇2β1 L (β∗ , φ(β∗ , ε)) nt,ν→µ βbt−1,ν − β∗ + oP (1), where (i) follows from Assumption 4.3: for any β ∈ Rd , 0=

d E [∇β1 L (β, φ(β, ε))] = ∇β1 E [∇β1 L (β, φ(β, ε))] + ∇β2 E [∇β1 L (β, φ(β, ε))] , dβ

and (ii) follows from Assumption 4.6. Combining the above estimates, and denoting H∗ = E ∇2β1 L (β∗ , φ(β∗ , ε)) , we obtain that √

√

X nt H∗−1 ∇β1 L (β∗ , Zt,ν→µ ) + oP (1) ν∈Nµ,in nt,ν→µ

nt βbt,µ − β∗ = − P

ν∈Nµ,in

√

X nt,ν→µ X nt H∗−1 ∇β1 L (β∗ , φ(β∗ , εt,ν→µ,i )) n t,ν→µ ν∈Nµ,in ν∈Nµ,in i=1 √ X nt nt,ν→µ βbt−1,ν − β∗ + oP (1). +P ν∈Nµ,in nt,ν→µ

=−P

ν∈Nµ,in

For µ ∈ Ml , define ut,µ = − P

nt,ν→µ

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

X i=1

H∗−1 ∇β1 L (β∗ , φ(β∗ , εt,ν→µ,i ))

and recall the definition of P t from Eq. (18). Using Assumption 4.2, the above expression for can be further simplified to √ √ X √ b nt βt,µ − β∗ = nt P t,µ,ν βbt−1,ν − β∗ + nt ut,µ + oP (1).

√

nt (βbt,µ − β∗ ) (49)

ν∈M

Since βbt,µ = βbt−1,µ for µ ∈ Mu , as long as we define ∀µ ∈ Mu ,

ut,µ = 0,

Eq. (49) will still be valid for any µ ∈ Mu . It finally follows that √ √ b √ nt βt − βb−1 = (P t ⊗ Ip ) nt βbt−1 − βb−1 + nt ut + oP (1). According to our data generating process, ut is asymptotically normal (see Appendix C.1) and independent √ d of βbt−1 . In particular, we have that nt ut → N(0, Vt ) as the sample sizes tend to infinity, where we recall that Vt is defined according to Eq. (21). Further, since nt /nt−1 → bt,t−1 by Assumption 4.2, we deduce √ d from our induction hypothesis that nt (βbt−1 − βb−1 ) → bt,t−1 Σt−1 . We finally conclude that βbt is also asymptotically normal, in the sense that √ b d nt βt − βb−1 → N(0, Σt ), Σt = bt,t−1 (P t ⊗ Id )Σt−1 (P t ⊗ Id )⊤ + Vt . This closes the induction argument and completes the proof of Theorem 4.1.

A.4

Proof of Theorem 4.2

To facilitate analysis, we first establish some useful lower and upper bounds on ΣT,µ,µ in the lemma below. Lemma A.10. For t ≥ s ≥ 1, define Jt,s = P t P t−1 · · · P s , with the convention that Jt,t+1 = IK . Recall that V∗ is defined in Eq. (19), and that p0,µ , pt,ν→µ and bt,s are defined in Assumption 4.2. Then, the following holds: (a) For any µ ∈ Ml and T ∈ N, we have ΣT,µ,µ ⪰

! P T X bT,t ( ν∈Ml JT,t+1,µ,ν )2 bT,0 P P +P V∗ . ν∈Ml ν∈M p0,ν m∈Nν,in pt,m→ν t=1

(b) For any µ ∈ Ml and T ∈ N, we have ΣT,µ,µ ⪯

T X bT,t ( t=1

2 ν∈Ml JT,t+1,µ,ν )

inf ν∈Ml

m∈Nν,in pt,m→ν

bT,0

inf ν∈M p0,ν

V∗ .

Proof of Lemma A.10. We prove Lemma A.10 in Section C.2. We now continue the proof of Theorem 4.2. Under Assumption 4.2, our lower and upper bounds in Lemma A.10 implies that, for all µ ∈ Ml and T ∈ N:   !2 !2 T T X X X X Tr(Σ ) T,µ,µ ΣT,µ,µ ≃  JT,t+1,µ,ν + 1 V∗ =⇒ ≃ JT,t+1,µ,ν + 1, Tr(V∗ ) t=1 t=1 ν∈Ml

ν∈Ml

where ≃ hides constants that are independent of T . Therefore, a model µ ∈ Ml collapses as T → ∞ if and only if !2 T X X lim sup JT,t+1,µ,ν = ∞. T →∞ t=1

ν∈Ml

∞ We first show that ∀µ ∈ M∞ l will collapse. By definition of Ml , we know that for all t ≥ s and ν ∈ Mu , Jt,s,µ,ν = 0. Therefore, for all t ≥ s we have X Jt,s,µ,ν = 1, ν∈Ml

PT P and hence t=1 ( ν∈Ml JT,t+1,µ,ν )2 = T → ∞ as T → ∞. c ′ ∞ We next consider µ ∈ Ml \M∞ l . If µ ∈ Ml , i.e., there is a directed path from some ν ∈ Ml to µ, then Assumption 4.2 implies that there exists ε > 0 and T0 ∈ N, such that Jt+T0 ,t,µ,ν ′ ≥ ε, ∀t ∈ N. This further implies that as long as t − s ≥ T0 , we have X X Jt,s,µ,ν ≥ Jt,t−T0 ,µ,ν ′ Jt−T0 −1,s,ν ′ ,ν = Jt,t−T0 ,µ,ν ′ ≥ ε. ν∈Ml

ν∈Ml

Therefore, as T → ∞,

T X

t=1

ν∈Ml

≥ (T − T0 )ε2 → ∞,

JT,t+1,µ,ν

which means that model µ will collapse. This proves part (a) of Theorem 4.2. ∞ Finally, we consider µ ∈ Mnc / M∞ l , i.e., µ ∈ l and there is no directed path from Ml to µ. By definition, there is a directed path from Mu to µ, meaning that there exists ε(µ) > 0 and T0 (µ) ∈ N, such that X X Jt+T0 (µ),t,µ,ν ≥ ε(µ) ⇐⇒ Jt+T0 (µ),t,µ,ν ≤ 1 − ε(µ) ν∈Mu

ν∈Ml

for all t ∈ N. Since Jt,s,ν,ν = 1 for all ν ∈ Mu , the above inequality remains true if T0 (µ) is replaced by any T ≥ T0 (µ). Of course, the above inequality also holds for all µ ∈ Mu , so it actually holds for all µ ∈ / Mcl . Let us now define ε = minc ε(µ), T0 = maxc T0 (µ). µ∈M / l

µ∈M / l

We deduce that for all µ ∈ / Mcl and t − s ≥ T0 : X ν∈Ml

Jt,s,µ,ν ≤ 1 − ε.

Using this inequality, we know that for all µ ∈ / Mcl and t − s ≥ 2T0 + 1, X X X Jt,s,µ,ν = Jt,t−T0 ,µ,ν ′ Jt−T0 −1,s,ν ′ ,ν ν∈Ml

ν ′ ∈M / cl ∪Mu

≤ (1 − ε)

ν∈Ml

X ν ′ ∈M / cl ∪Mu

Jt,t−T0 ,µ,ν ′ ≤ (1 − ε)2 .

Iterating further, we can use induction to show that for all µ ∈ / Mcl and t − s ≥ kT0 + k − 1, X Jt,s,µ,ν ≤ (1 − ε)k . ν∈Ml

This immediately implies that !2

T X

T →∞ t=1

ν∈Ml

lim sup

JT,t+1,µ,ν

< ∞,

so that model µ will not collapse. This proves part (b) of Theorem 4.2. 40

A.5

Proof of Theorem 4.3

To begin with, we state some general conditions under which Assumptions 4.3–4.6 are satisfied. In this proof, we assume that the response variable y can be represented as follows: y = φ(β∗⊤ x, ε). Note that, with a slight abuse of notation, we still use φ to denote the link function and µε to denote the joint distribution of (x, ε). The above representation of y yields the following expression for the loss function: L β1 , φ(β2⊤ x, ε) = −φ(β2⊤ x, ε) · β1⊤ x + A(β1⊤ x). Proposition A.11. For a GLM as specified in Eq. (23), assume that (i) E[∥x∥2 ] < ∞. Further, for any compact set Ω, h i E sup ∥x∥2 |φ(β ⊤ x, ε)| < ∞, β∈Ω

h i E sup |A(β ⊤ x)| < ∞. β∈Ω

(ii) The function classes {(x, ε) 7→ φ(β ⊤ x, ε) : β ∈ Ω}, {x 7→ A(β ⊤ x) : β ∈ Ω} and {x 7→ A′′ (β ⊤ x)xx⊤ : β ∈ Ω} are all µε -Glivenko-Cantelli for any compact set Ω. (iii) The function class {(x, ε) 7→ φ(β ⊤ x, ε)x : β ∈ Ω} is µε -Donsker for any compact set Ω. (iv) For all β ∈ Rd , both E[A′ (β ⊤ x)2 xx⊤ ] and E[A′′ (β ⊤ x)xx⊤ ] are finite. Further, the mapping β 7→ E[A′′ (β ⊤ x)xx⊤ ] is continuous. Then, Assumptions 4.3-4.6 are satisfied. Proof. We prove Lemma A.11 in Section C.3. We next use Lemma A.11 to show that the three GLMs (linear regression, logistic regression and Poisson regression) all satisfy Assumptions 4.3-4.6. Throughout, for a compact set Ω ∈ Rd , let rΩ > 0 be such that Ω ⊂ B(0, rΩ ). Linear regression. In this setting, we have A(ξ) = ξ 2 /2 and φ(β ⊤ x, ε) = β ⊤ x + ε with E[ε] = 0, 4 E[ε2 ] < ∞. We next show that Assumptions 4.3-4.6 are satisfied, provided that E[∥x∥2 ] < ∞. To this end, it suffices to verify conditions (i)-(iv) in Lemma A.11. Condition (i) can be verified by noting that 2

sup ∥x∥2 |φ(β ⊤ x, ε)| ≤ rΩ ∥x∥2 + ∥x∥2 |ε|,

β∈Ω

sup |A(β ⊤ x)| ≤

β∈Ω

1 2 2 r ∥x∥2 . 2 Ω

Condition (ii) can be verified by applying Lemma C.1 to the function classes {(x, ε) 7→ β ⊤ x|β ∈ Ω} and 4 {(x, ε) 7→ ε}, and noting that E[∥x∥2 ] is finite. Similarly, (iii) follows from Lemma C.2, the compactness of 4 Ω, and our assumption E[∥x∥2 ] < ∞. Finally, (iv) can be verified by direct calculation. Logistic regression.

In this setting, we have A(ξ) = log(1 + eξ ), and ⊤

⊤

φ(β ⊤ x, ε) = 1{ε ≤ A′ (β ⊤ x)} = 1{ε ≤ eβ x /(1 + eβ x )}, where ε ∼ Unif[0, 1]. We show that Assumptions 4.3-4.6 are satisfied, provided that E[∥x∥32 ] < ∞. To this end, it suffices to verify conditions (i)-(iv) in the statement of Lemma A.11. Condition (i) can be easily verified, since φ is uniformly bounded by 1 and A(ξ) ≤ 1 + |ξ| for all ξ ∈ R. To verify condition (ii), we first note that the function class {(x, ε) 7→ β ⊤ x : β ∈ Ω} is µε -Glivenko-Cantelli since E[∥x∥32 ] < ∞, and A′′ (ξ) ∈ [0, 1/4] for all ξ ∈ R. We can then invoke Lemma C.1 to show that the function classes {(x, ε) 7→ A(β ⊤ x) : β ∈ Ω},

{(x, ε) 7→ A′′ (β ⊤ x)xx⊤ : β ∈ Ω}, 41

{(x, ε) 7→ φ(β ⊤ x, ε) : β ∈ Ω}

are µε -Glivenko-Cantelli. To verify condition (iii), we use Lemma C.3 to show that the function class {(x, ε) 7→ φ(β ⊤ x, ε)x : β ∈ Ω} is µε -Donsker. To apply Lemma C.3, we set f1 (t) = A′ (t) and fk (t) = 0 for k ≥ 2. We also set L(t) = 1, f (t) = f1 (t) and g(x) = x. Since A′ (z) ∈ [0, 1] and A′′ (z) ∈ [0, 1/4] for all z ∈ R, conditions (i) and (ii) in Lemma C.3 are automatically satisfied. Further, condition (iii) in Lemma C.3 is verified by our assumption E[∥x∥32 ] < ∞. Finally, condition (iv) can be established using the boundedness of A′ and A′′ and the dominated convergence theorem. Poisson regression. Unif[0, 1], and

In this setting, we have A(ξ) = eξ , and φ(β ⊤ x, ε) = fk (z) :=

∞ X exp(mz − ez ) m=k

P∞

k=1 1ε≤fk (β ⊤ x) , where ε ∼

= P (Poisson(ez ) ≥ k) .

We show that Assumptions 4.3-4.6 are satisfied, provided that E[exp(R∥x∥2 )] < ∞ for any R > 0. To this end, it suffices to verify conditions (i)-(iv) in Lemma A.11. For condition (i), note that our assumption already implies E[∥x∥2 ] < ∞. Further, by the definition of φ and the monotonicity of the fk ’s, we have ∞ h i h h X i i E sup ∥x∥2 |φ(β ⊤ x, ε)| ≤ E ∥x∥2 fk rΩ ∥x∥2 = E ∥x∥2 exp rΩ ∥x∥2 < ∞, β∈Ω

k=1

h i h i E sup |A(β ⊤ x)| ≤ E exp rΩ ∥x∥2 < ∞, β∈Ω

where rΩ = sup{∥β∥2 : β ∈ Ω}. Condition (iv) follows immediately from our moment assumption. To verify conditions (ii) and (iii), we first use the same argument as that in logistic regression to establish that the function classes {(x, ε) 7→ A(β ⊤ x) : β ∈ Ω}, {(x, ε) 7→ A′′ (β ⊤ x)xx⊤ : β ∈ Ω} are µε -Glivenko-Cantelli. It now remains to show that the function classes {(x, ε) 7→ φ(β ⊤ x, ε) : β ∈ Ω},

{(x, ε) 7→ φ(β ⊤ x, ε)x : β ∈ Ω}

are µε -Donsker, since any µε -Donsker function class is also µε -Glivenko-Cantelli. To this end, we invoke Lemma C.3, with g(x) = 1 or x. We next show that conditions (i)-(iii) in the statement of Lemma C.3 are satisfied. To this end, we set L(t) = L and f (t) = e3t + 3e2t + et for an absolute constant L > 0 (to be determined later). This already gives E[g(x)2 ∥x∥2 L(R∥x∥2 )] < ∞ and E[g(x)2 f (R∥x∥2 )] < ∞ for any R > 0, completing the proof of (iii) in Lemma C.3. To establish upper bounds on fk′ , note that by direct calculation: k (ii) √ 1 exp(kz − ez ) (i) (k/e) ≤ sup λk e−λ = ≤ L k, fk′ (z) = (k − 1)! (k − 1)! λ>0 (k − 1)! where (i) follows from direct computation, and (ii) is due to Stirling’s formula. Taking c1 = −1/2 then verifies condition (i) in Lemma C.3. To establish upper bounds on fk , we resort to Markov inequality: fk (z) = P (Poisson(ez ) ≥ k) ≤ k −3 E Poisson(ez )3 = k −3 (e3z + 3e2z + ez ) = k −3 f (z), we can then take c2 = 3. This concludes the proof.

Technical lemmas for the linear regression setting

In this section, we provide proofs of the technical lemmas supporting our linear regression results.

B.1

Proof of Lemma A.1

By Assumption 3.2, we know that (Xs,ν→µ , εs,µ→ν ) is a subset of (Xs,ν , εs,ν ). Therefore, for any fixed s ∈ [t] and ν ∈ M, there exists a partition of (Xs,ν , εs,ν ): N

s,ν (Xs,ν , εs,ν ) = ∪j=1 (Xs,ν,j , εs,ν,j ),

such that Ns,ν ≤ 2|M| = 2K , and for each µ ∈ M, there exists a subset As,ν→µ ⊂ [Ns,ν ], such that (Xs,ν→µ , εs,µ→ν ) = ∪j∈As,ν→µ (Xs,ν,j , εs,ν,j ). (If (ν, µ) ∈ / E, one can simply take As,ν→µ = ∅.) We can therefore define {(Zt,i , ζt,i ) : 1 ≤ i ≤ Nt } := {(Xs,ν,j , εs,ν,j ) : s ∈ [t], ν ∈ M, j ∈ As,ν→µ for some (ν, µ) ∈ E}. This proves part 1. To show part 2, note that for any fixed pair (s, ν), we have Ns,ν ≤ 2K . Hence, Nt ≤ Ct, where the constant C = K2K does not depend on t. Finally, by our construction, we know that ζt = ⊤ ⊤ ⊤ (ζt,1 , ζt,2 , · · · , ζt,N )⊤ is a subvector of εt , which naturally implies that Cov(ζt ) ⪰ ρ2 I by Assumption 3.4. t This completes the proof of this lemma.

B.2

Proof of Lemma A.3

By Eq. (25), we have for all µ ∈ Ml : βbt,µ − βb0,µ = t (i) X

t X X

Ωt,s+1,µ,ν vs,ν

s=1 ν∈M

Ωt,s+1,µ,ν

s=1 ν∈Ml

t X X

m∈Nν,in

s=1 ν∈Ml m∈Nν,in

t X X

(ii)

⊤ Xs,m→ν εs,m→ν

m∈Nν,in

⊤ ⊤ Ωt,s+1,µ,ν Ts,ν,m (Xs,m→ν Xs,m→ν )−1 Xs,m→ν εs,m→ν

⊤ Ωt,s+1,µ,ν Ts,ν,m (Xs,m→ν Xs,m→ν )−1

i=1 (s,m→ν,h)∈Pt,i Nt X

−1 X

hs,m→ν

s=1 ν∈Ml m∈Nν,in Nt X

⊤ Xs,m→ν Xs,m→ν

⊤ Xs,m→ν,h εs,m→ν,h

h=1

⊤ ⊤ Ωt,s+1,µ,ν Ts,ν,m (Xs,m→ν Xs,m→ν )−1 Zt,i ζt,i

⊤ ⊤ At,i,µ (Zt,i Zt,i )† Zt,i ζt,i ,

i=1

where (i) follows from the definition of vs,ν , (ii) follows from Lemma A.1 and A.2, and the last line is due ⊤ ⊤ ⊤ ⊤ to the definition of At,i,µ , and the identity Zt,i Zt,i (Zt,i Zt,i )† Zt,i = Zt,i . This completes the proof.

B.3

Proof of Lemma A.4

By definition, we have Nt X

At,i,µ =

i=1

Nt X

i=1 (s,m→ν,h)∈Pt,i

t X X

⊤ ⊤ Ωt,s+1,µ,ν Ts,ν,m (Xs,m→ν Xs,m→ν )−1 Zt,i Zt,i hs,m→ν

s=1 ν∈Ml m∈Nν,in

⊤ Ωt,s+1,µ,ν Ts,ν,m (Xs,m→ν Xs,m→ν )−1

X h=1

⊤ Xs,m→ν,h Xs,m→ν,h

t X X

s=1 ν∈Ml m∈Nν,in

t X X X s=1 ν∈Ml m∈M

⊤ ⊤ Ωt,s+1,µ,ν Ts,ν,m (Xs,m→ν Xs,m→ν )−1 (Xs,m→ν Xs,m→ν )

Ωt,s+1,µ,ν Ts,ν,m ,

where the last line follows from the fact that Ts,ν,m = 0d×d if m ∈ / Nν,in . This completes the proof.

B.4

Proof of Lemma A.5

In this section we prove Lemma A.5. The first inequality follows immediately from triangle inequality: X E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 op

v1 →v2 →···→vx →µ∈L∗ x

≤

X v1 →v2 →···→vx →µ∈L∗ x

E[Tt,µ,vx ] op E[Tt−1,vx ,vx−1 ] op · · · E[Tt−x+1,v2 ,v1 ] op .

(50)

By Assumption 3.5, for all s ∈ N+ and (ν1 , ν2 ) ∈ E, the two matrices E[Ts,ν2 ,ν1 ] and ps,ν1 →ν2 Id are close, and we provide an upper bound for their difference: E Ts,ν2 ,ν1 − ps,ν1 →ν2 Id op (i)

≤ E [∥Ts,ν2 ,ν1 − ps,ν1 →ν2 Id ∥op ] h i ≤ E Ts,ν2 ,ν1 − ps,ν1 →ν2 Id 1{∥Ts,ν2 ,ν1 − ps,ν1 →ν2 Id ∥op ≤ δ} op h i +E Ts,ν2 ,ν1 − ps,ν1 →ν2 Id 1{∥Ts,ν2 ,ν1 − ps,ν1 →ν2 Id ∥op > δ} op

(ii)

(51)

(iii)

≤ δ + 2P (∥Ts,ν2 ,ν1 − ps,ν1 →ν2 Id ∥op > δ)

≤ δ + 2γ1 ,

where we recall ps,ν1 →ν2 is from Assumption 3.3. In the above display, (i) is by Jensen’s inequality, (ii) is by triangle inequality, and (iii) is because ∥Ts,ν2 ,ν1 ∥op ≤ 1. Leveraging Eq. (51), we derive the following upper bound for the quantity displayed in the last line of Eq. (50): X E[Tt,µ,vx op E[Tt−1,vx ,vx−1 ] op · · · E[Tt−x+1,v2 ,v1 ] op v1 →v2 →···→vx →µ∈L∗ x

≤

X v1 →v2 →···→vx →µ∈L∗ x

(pt,vx →µ + δ + 2γ1 )(pt−1,vx−1 →vx + δ + 2γ1 ) · · · (pt−x+1,v1 →v2 + δ + 2γ1 ).

(52)

x Recall that K is the number of models. We let ax = K⌈ K ⌉, and note that X (pt,vx →µ + δ + 2γ1 )(pt−1,vx−1 →vx + δ + 2γ1 ) · · · (pt−x+1,v1 →v2 + δ + 2γ1 ) v1 →v2 →···→vx →µ∈L∗ x

⌈x/K⌉

 ≤

sup

vK+1 ∈Ml \M∞ l , s≥K v1 →v2 ···→vK →vK+1 ∈L∗ K

(ps,vK →vK+1 + δ + 2γ1 ) · · · (ps−K+1,v1 →v2 + δ + 2γ1 )

" ×

sup vx−ax +1 ∈Ml \M∞ l , s≥x−ax v →v ···→v 1

∗ x−ax →vx−ax +1 ∈Lx−ax

(ps,vx−ax →vx−ax +1 + δ + 2γ1 ) · · ·

# · · · (ps−x+ax +1,v1 →v2 + δ + 2γ1 ) ∨ 1  ≤

sup

vK+1 ∈Ml \M∞ l , s≥K v1 →v2 ···→vK →vK+1 ∈L∗ K

ps,vK →vK+1 · · · ps−K+1,v1 →v2 +

K X K δ + 2γ1 i

i=1

⌈x/K⌉ 

" ×

sup vx−ax +1 ∈Ml \M∞ l , s≥x−ax v →v ···→v 1

K X K δ + 2γ1 i i=1

∗ x−ax →vx−ax +1 ∈Lx−ax

ps,vx−ax →vx−ax +1 · · · ps−x+ax +1,v1 →v2

∨ 1. (53)

where the last inequality is because under Assumption 3.3, for both k ∈ {K, x − ax } (note that x − ax ∈ {0, 1, · · · , K − 1}) and a fixed vk+1 (that is to say, we do not sum over vk+1 ), it holds that X (ps,vk →vk+1 + δ + 2γ1 ) · · · (ps−k+1,v1 →v2 + δ + 2γ1 ) v1 →v2 ···→vk →vk+1 ∈L∗ k

X v1 →v2 ···→vk →vk+1 ∈L∗ k

≤

X v1 →v2 ···→vk →vk+1 ∈L∗ k

ps,vk →vk+1 · · · ps−k+1,v1 →v2 ps,vk →vk+1 · · · ps−k+1,v1 →v2

k X

i Y

i=1 1≤j1 <···<ji ≤k

ps,vk →vk+1 · · · ps−k+1,v1 →v2 +

δ + 2γ1 p r=1 s−k+jr ,vjr →vjr +1

K X K δ + 2γ1 i i=1

ν∈Nµ,in pt,ν→µ = 1 for all µ ∈ Ml , hence

where the last inequality is because X v1 →v2 ···→vk →vk+1 ∈L∗ k

ps,vk →vk+1 · · · ps−k+1,v1 →v2 ≤ 1

for k ∈ {K, x − ax }.

(54)

In addition, under Assumption 3.3, i Y

δ + 2γ1

p r=1 s−k+jr ,vjr →vjr +1

≤

δ + 2γ i 1

For any fixed vK+1 and a path v1 → v2 · · · → vK → vK+1 ∈ L∗K , by definition of L∗K we know vK+1 ∈ Ml \M∞ l . Therefore, there exists a path u1 → u2 → · · · → uJ → vK+1 with J ≤ K − 1, such that u1 ∈ Mu . Observe that (the sum below is taken over v1 , v2 , · · · , vK ) X ps,vK →vK+1 · · · ps−K+1,v1 →v2 + ps,uJ →vK+1 ps−1,uJ−1 →uJ · · · ps−J+1,u1 →u2 ≤ 1, v1 →v2 ···→vK →vK+1 ∈L∗ K

ps,uJ →vK+1 ps−1,uJ−1 →uJ · · · ps−J+1,u1 →u2 ≥ αK , We then conclude that fixing vK+1 , X v1 →v2 ···→vK →vK+1 ∈L∗ K

by Assumption 3.3.

ps,vK →vK+1 · · · ps−K+1,v1 →v2 ≤ 1 − αK .

(55)

Putting together Eqs. (53) to (55), we arrive at the following upper bound: X (pt,vx →µ + δ + 2γ1 )(pt−1,vx−1 →vx + δ + 2γ1 ) · · · (pt−x+1,v1 →v2 + δ + 2γ1 ) v1 →v2 →···→vx →µ∈L∗ x

(56)

⌈x/K⌉ δ + 2γ1 K δ + 2γ1 K K ≤ 1+ 1+ −α , α α

which under Assumption 3.5 is no larger than 2c(α, K)⌈x/K⌉ . The proof then follows via putting together Eqs. (50), (52) and (56). Finally, the “on the other hand” part follows directly from the definition of L∗x .

B.5

Proof of Lemma A.6

The proof proceeds by analyzing the first sum in Eq. (33) separately for different values of ℓ. We first introduce a decomposition of this sum: First sum in Eq. (33) =

t t X X X

ℓ=1 x=ℓ v∈M∞ l v1 →v2 →···→vx →µ∈Lx,ℓ,v

+ 1{µ ∈ M∞ l }

t X

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1

x=1 v1 →v2 →···→vx →µ∈Lx,0,µ

(57)

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 .

∞ We then separately consider two cases µ ∈ Mcl \M∞ l and µ ∈ Ml .

Case I: µ ∈ Mcl \M∞ In this case, we only need to consider terms associated with a strictly positive ℓ l . in Eq. (57), as the second sum is zero. Note that for a fixed ℓ ∈ N+ , t X X

x=ℓ v∈M∞ l v1 →v2 →···→vx →µ∈Lx,ℓ,v

= (t − ℓ + 1)

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 X

E[Tt,µ,vℓ−1 ]E[Tt−1,vℓ−1 ,vℓ−2 ]

(58)

∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

· · · E[Tt−ℓ+2,v2 ,v1 ]E[Tt−ℓ+1,v1 ,v ]. We then prove that the matrix given in Eq. (58) can be well approximated by a scalar times an identity matrix. Observe that X X E[Tt,µ,vℓ−1 ]E[Tt−1,vℓ−1 ,vℓ−2 ] · · · E[Tt−ℓ+2,v2 ,v1 ]E[Tt−ℓ+1,v1 ,v ] ∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

− pt,vℓ−1 →µ pt−1,vℓ−2 →vℓ−1 · · · pt−ℓ+2,v1 →v2 pt−ℓ+1,v→v1 Id op X X E[Tt,µ,vℓ−1 ]E[Tt−1,vℓ−1 ,vℓ−2 ] · · · E[Tt−ℓ+2,v2 ,v1 ]E[Tt−ℓ+1,v1 ,v ]. (59) ≤ ∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

− pt,vℓ−1 →µ pt−1,vℓ−2 →vℓ−1 · · · pt−ℓ+2,v1 →v2 pt−ℓ+1,v→v1 Id

. op

For i ∈ {0, 1, · · · , ℓ − 1}, we express E[Tt−i,νℓ−i ,νℓ−i−1 ] (treating µ as vℓ and v as v0 ) as the sum of pt−i,νℓ−i−1 →νℓ−i Id and E[Tt−i,νℓ−i ,νℓ−i−1 ] − pt−i,νℓ−i−1 →νℓ−i Id . We then apply the operator norm bound ∥E[Tt−i,νℓ−i ,νℓ−i−1 ] − pt−i,νℓ−i−1 →νℓ−i Id ∥op ≤ δ + 2γ1 from Eq. (51), and obtain an upper bound for the final line of Eq. (59). Specifically, let Di = E[Tt−i,νℓ−i ,νℓ−i−1 ] − pt−i,νℓ−i−1 →νℓ−i Id for i ∈ {0, 1, · · · , ℓ − 1}, then The last line of Eq. (59) X X =

(D0 + pt,vℓ−1 →vℓ Id ) · · · (Dℓ−1 + pt−ℓ+1,v0 →v1 Id )

∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

− pt,vℓ−1 →vℓ · · · pt−ℓ+1,v0 →v1 Id ≤ ≤

op ℓ−1 Y

pt−i,vℓ−i−1 →vℓ−i

∗ i=0 v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

ℓ−1 Y

ℓ X

j Y

j=1 0≤h1 <···<hj ≤ℓ−1 b=1

pt−i,vℓ−i−1 →vℓ−i

∗ i=0 v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

ℓ X ℓ δ + 2γ1 j j=1

δ + 2γ1 pt−hb ,νℓ−hb −1 →νℓ−hb

(60)

Let aℓ−1 = K⌈(ℓ − 1)/K⌉. Note that when ℓ − 1 ≥ K, similar to the derivation of Eq. (55), it holds that X

ℓ−1 Y

pt−i,vℓ−i−1 →vℓ−i

∗ i=0 v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

≤K

ℓ−2 Y

pt−i,vℓ−i−1 →vℓ−i

i=0 v1 →v2 →···→vℓ−1 →µ∈L∗ ℓ−1

⌈(ℓ−1)/K⌉

 ≤K

sup

vK+1 ∈Ml \M∞ l , s≥K v1 →v2 ···→vK →vK+1 ∈L∗ K

ps,vK →vK+1 ps−1,vK−1 →vK · · · ps−K+1,v1 →v2 

" ×

# X

sup vℓ−aℓ−1 ∈Ml \M∞ l , s≥ℓ−1−aℓ−1 v →···→v 1

∗ ℓ−aℓ−1 ∈Lℓ−1−a ℓ−1

ps,vℓ−1−aℓ−1 →vℓ−aℓ−1 · · · ps−ℓ+aℓ−1 +2,v1 →v2 ∨ 1

≤ K(1 − αK )⌈(ℓ−1)/K⌉ .

(61)

Substituting Eq. (61) into Eq. (60), we conclude that δ + 2γ1 ℓ −1 the last line of Eq. (59) ≤K(1 − α ) 1+ α K δ + 2γ1 ℓ K ℓ/K K ℓ/K ≤ (1 − α ) − (1 − α ) . 1+ α (1 − αK )1/K+1 K ⌈(ℓ−1)/K⌉

(62)

For the sake of simplicity, for ℓ ∈ N+ , we define the following d × d matrices: X X Qℓ = E[Tt,µ,vℓ−1 ]E[Tt−1,vℓ−1 ,vℓ−2 ] · · · E[Tt−ℓ+2,v2 ,v1 ]E[Tt−ℓ+1,v1 ,v ], ∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

Q∗ℓ =

∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

pt,vℓ−1 →µ pt−1,vℓ−2 →vℓ−1 · · · pt−ℓ+2,v1 →v2 pt−ℓ+1,v→v1 Id .

Under Assumption 3.5, Q∗ℓ can be regarded as a population version of Qℓ . Note that t X t X X

ℓ=1 x=ℓ v∈M∞ l v1 →v2 →···→vx →µ∈Lx,ℓ,v

t X E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 = (t − ℓ + 1)Qℓ . ℓ=1

As µ ∈ Mcl \M∞ length ℓ ∈ [K] from l , by definition there exists a path v → v1 → v2 → · · · → vℓ−1 → µ of P K ∞ ∗ K some v ∈ M∞ to µ, where the intermediate nodes v , v , · · · , v ∈ / M . Therefore, 1 2 ℓ−1 l l ℓ=1 ∥Qℓ ∥op ≥ α . Together, these results yield a lower bound for the matrix largest eigenvalue (recall σmax (A) refers to the largest eigenvalue of A): σmax

t X (t − ℓ + 1)Qℓ ℓ=1

t t X (i) X ≥ (t − ℓ + 1)∥Q∗ℓ ∥op − (t − ℓ + 1)∥Qℓ − Q∗ℓ ∥op ℓ=1

ℓ=1

(ii)

≥ (t − K + 1)αK − t

∞ X

δ + 2γ1 ℓ K K ℓ/K K ℓ/K 1 + (1 − α ) − (1 − α ) α (1 − αK )1/K+1 ℓ=1

(iii)

tK(δ + 2γ1 ) α(1 − αK )1/K+1 (1 − c1 )2 K(δ + 2γ1 ) − αK (K − 1) = t c2 − αK (K − 1), = t αK − α(1 − αK )1/K+1 (1 − c1 )2 ≥ (t − K + 1)αK −

where (i) uses the fact that Q∗ℓ is a positive scalar times an identity matrix, (ii) is by Eq. (62), and (iii) is by Assumption 3.5. The proof is done. Case II: µ ∈ M∞ l . it holds that

In this case, the second sum in Eq. (57) is zero. On the other hand, when µ ∈ M∞ l ,

t X

x=1 v1 →v2 →···→vx →µ∈Lx,0,µ

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1 = tId .

(63)

Therefore, t x X X X

x=1 v∈M∞ l ℓ=0 v1 →v2 →···→vx →µ∈Lx,ℓ,v

E Tt,µ,vx E Tt−1,vx ,vx−1 · · · E Tt−x+1,v2 ,v1

≥ t > t c2 − αK (K − 1). The proof is done. Finally, the “on the other hand” part follows directly from the definition of Mnc l .

B.6

Proof of Lemma A.8

In this section we prove Lemma A.8. According to Lemma A.7, Nt hX i E At,i,µ max 1 − 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } i=1

v∈M

t hX X X i ⊤ =E Ωt,s+1,µ,ν Ts,ν,m max 1 − 1{Xs,v Xs,v /ns,v ⪯ κId } v∈M

s=1 ν∈Ml m∈M

t X

h i ⊤ E Tt,µ,vx Tt−1,vx ,vx−1 · · · Tt−x+1,v2 ,v1 max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } v∈M

x=1 v1 →···→vx →µ∈Lx

t X

⊤ E[Tt,µ,vx ]E[Tt−1,vx ,vx−1 ] · · · E[Tt−x+1,v2 ,v1 max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } ]. v∈M

x=1 v1 →···→vx →µ∈Lx

Recall Lx,ℓ,v and L∗x are defined in Eqs. (31) and (32). With these definitions, we have Nt hX i E At,i,µ max 1 − 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId }

i=1 t X

v∈M

x X X

x=1 v∈M∞ l ℓ=0 v1 →···→vx →µ∈Lx,ℓ,v

t X

x=1 v1 →···→vx →µ∈L∗ x

⊤ E Tt,µ,vx · · · E Tt−x+1,v2 ,v1 max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } v∈M

⊤ E Tt,µ,vx · · · E Tt−x+1,v2 ,v1 max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } . v∈M

(64) We denote by S1 and S2 the first and second sums in Eq. (64), respectively. We separately upper bound these two terms below, starting from S2 . ∞ Upper bounding S2 . When µ ∈ M∞ l we have S2 = 0d×d . In what follows, we assume µ ∈ Ml \Ml . Note that for all s ∈ N+ and (ν1 , ν2 ) ∈ E, by Eq. (51) we have ∥E[Ts,ν2 ,ν1 − ps,ν1 →ν2 Id ]∥op ≤ δ + 2γ1 . In addition,

⊤ E Ts,ν2 ,ν1 max 1 − 1{Xs,v Xs,v /ns,v ⪯ κId } v∈M

(i) op

⊤ ≤ E ∥Ts,ν2 ,ν1 ∥op max 1 − 1{Xs,v Xs,v /ns,v ⪯ κId } v∈M

(ii)

⊤ ≤ 1 − P ∩v∈M {Xs,v Xs,v /ns,v ⪯ κId }

(iii)

≤ Kγ2 , where (i) is by Jensen’s inequality, (ii) is because ∥Ts,ν2 ,ν1 ∥op ≤ 1, and (iii) is by Assumption 3.5. As a consequence, ∥S2 ∥op ≤ ≤ ≤

t X

x=1 v1 →···→vx →µ∈L∗ x t X

x=1 v1 →···→vx →µ∈L∗ x t Kγ2 X α x=1

⊤ ∥E Tt,µ,vx ∥op · · · ∥E Tt−x+1,v2 ,v1 max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } ∥op v∈M

Kγ2 (pt,vx →µ + δ + 2γ1 )(pt−1,vx−1 →vx + δ + 2γ1 ) · · · (pt−x+2,v2 →v3 + δ + 2γ1 )

X v1 →···→vx →µ∈L∗ x

(pt,vx →µ + δ + 2γ1 )(pt−1,vx−1 →vx + δ + 2γ1 ) · · · (pt−x+1,v1 →v2 + δ + 2γ1 ).

Similar to the derivation of Eq. (56), the above quantity is no larger than ⌈x/K⌉ t Kγ2 X δ + 2γ1 K δ + 2γ1 K 1+ 1+ − αK , α x=1 α α 2

γ2 which is further no larger than 4K under Assumption 3.5. In summary, for any µ ∈ Ml : αK+1

∥S2 ∥op ≤

4K 2 γ2 . αK+1

(65)

Upper bounding S1 . We then upper bound S1 . When µ ∈ Mnc l , we immediately see that S1 = 0d×d . c Below we consider µ ∈ Ml . Note that S1 = M1 + M2 , where M1 =

t x X X X

x=1 v∈M∞ l ℓ=1 v1 →···→vx →µ∈Lx,ℓ,v

⊤ E Tt,µ,vx · · · E Tt−x+1,v2 ,v1 max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } , v∈M

M2 =

t X

x=1 v1 →v2 →···→vx →µ∈Lx,0,µ

⊤ E Tt,µ,vx · · · E Tt−x+1,v2 ,v1 max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } . v∈M

When µ ∈ M∞ l , M1 = 0d×d . As for M2 , observe that M2 =

t h X i ⊤ E max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } Id , v∈M

x=1

t X

=⇒ ∥M2 ∥op ≤

x=1

⊤ 1 − P ∩x∈M {Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } ≤ t Kγ2 .

On the other hand, when µ ∈ Mcl \M∞ l , then M2 = 0d×d . We then upper bound ∥M1 ∥op . Note that M1 =

t X

M1,ℓ ,

ℓ=1 t X

M1,ℓ =

⊤ E Tt,µ,vx · · · E Tt−x+1,v2 ,v1 max 1 − 1{Xt−x+1,v Xt−x+1,v /nt−x+1,v ⪯ κId } . v∈M

x=ℓ v∈M∞ l v1 →···→vx →µ∈Lx,ℓ,v

In addition, for any ℓ ∈ [t − 1], X M1,ℓ = (t − ℓ + 1)

E[Tt,µ,vℓ−1 ]E[Tt−1,vℓ−1 ,vℓ−2 ]

∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

h i ⊤ · · · E[Tt−ℓ+2,v2 ,v1 ]E[Tt−ℓ+1,v1 ,v ]E max 1 − 1{Xt−x+1,η Xt−x+1,η /nt−x+1,η ⪯ κId } . η∈M

As a consequence, ∥M1,ℓ ∥op Kγ2 (t − ℓ + 1) X ≤

(66) X

∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

Leveraging Lemma A.5, we see that X X ∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

≤

E[Tt,µ,vℓ−1 ]E[Tt−1,vℓ−1 ,vℓ−2 ] · · · E[Tt−ℓ+2,v2 ,v1 ]E[Tt−ℓ+1,v1 ,v ]

∗ v∈M∞ l v1 →v2 →···→vℓ−1 →µ∈Lℓ−1 , v∈Nv1 ,in

(i)

≤K

X v1 →v2 →···→vℓ−1 →µ∈L∗ ℓ−1

E[Tt,µ,vℓ−1 ]E[Tt−1,vℓ−1 ,vℓ−2 ] · · · E[Tt−ℓ+2,v2 ,v1 ]E[Tt−ℓ+1,v1 ,v ]

. op

E[Tt,µ,vℓ−1 ] op · · · E[Tt−ℓ+2,v2 ,v1 ] op E[Tt−ℓ+1,v1 ,v ] op

E[Tt,µ,vℓ−1 ] op E[Tt−1,vℓ−1 ,vℓ−2 ] op · · · E[Tt−ℓ+2,v2 ,v1 ] op

⌈(ℓ−1)/K⌉ δ + 2γ1 K δ + 2γ1 K 1+ − αK , ≤K 1+ α α

(67)

where (i) is because for any s ∈ N+ and (µ1 , µ2 ) ∈ E, we have ∥E[Ts,µ2 ,µ1 ]∥op ≤ 1. Combining Eqs. (66) and (67), we conclude that for all ℓ ∈ [t − 1], ⌈(ℓ−1)/K⌉ δ + 2γ1 K δ + 2γ1 K ∥M1,ℓ ∥op ≤ K 2 γ2 (t − ℓ + 1) 1 + 1+ − αK . (68) α α When ℓ = t, note that M1,t =

E[Tt,µ,vt−1 ]E[Tt−1,vt−1 ,vt−2 ]

∗ v∈M∞ l v1 →v2 →···→vt−1 →µ∈Lt−1 , v∈Nv1 ,in

h i ⊤ · · · E[T2,v2 ,v1 ]E T1,v1 ,v max 1 − 1{X1,η X1,η /n1,η ⪯ κId } . η∈M

Using Lemma A.5 again, we see that X X ∥M1,t ∥op ≤

E[Tt,µ,vt−1 ] op E[Tt−1,vt−1 ,vt−2 ] op

∗ v∈M∞ l v1 →v2 →···→vt−1 →µ∈Lt−1 , v∈Nv1 ,in

h i ⊤ · · · × ∥E[T2,v2 ,v1 ] op E T1,v1 ,v max 1 − 1{X1,η X1,η /n1,η ⪯ κId } η∈M

(i)

≤ K 2 γ2

X v1 →v2 →···→vt−1 →µ∈L∗ t−1

E[Tt,µ,vt−1 ] op E[Tt−1,vt−1 ,vt−2 ] op · · · E[T2,v2 ,v1 ] op

(69)

⌈(ℓ−1)/K⌉ δ + 2γ1 K δ + 2γ1 K K ≤ K γ2 1 + 1+ −α . α α 2

Putting together Assumption 3.5, Eqs. (68) and (69), we conclude that for any µ ∈ Ml , ∥M1 ∥op ≤

t X ℓ=1

∥M1,ℓ ∥op ≤

4 t K 3 γ2 . αK

(70)

Finally, by Eqs. (65) and (70) we see that Nt hX i E At,i,µ max 1 − 1{Xs⊤t,i ,v Xst,i ,v /nst,i ,v ⪯ κId } v∈M

i=1

≤

4 t K 3 γ2 4K 2 γ2 9 t K 3 γ2 + K+1 + tKγ2 ≤ , K α α αK+1

since K ≥ 1 and α < 1. The proof is done.

B.7

Proof of Lemma A.9

If ν1 ∈ Mu , then one can verify that Ω∗t,s+1,ν1 ,ν2 = Ωt,s+1,ν1 ,ν2 = 1ν1 =ν2 Id . In the rest parts of the proof we consider ν1 ∈ Ml . The proof proceeds by induction on t − s. If t − s = 0, then by definition Ω∗s,t+1,µ,ν = Ωs,t+1,µ,ν = 1µ=ν Id . Now suppose Eq. (43) holds for t − s = k. When t − s = k + 1, by definition X X ∗ ∗ Ω∗t,s+1,ν1 ,ν2 = Tt,ν Ω∗t−1,s+1,m,ν2 = Tt,ν Ω∗t−1,s+1,m,ν2 , 1 ,m 1 ,m m∈M

Ωt,s+1,ν1 ,ν2 =

X m∈M

m∈Nν1 ,in

Tt,ν1 ,m Ωt−1,s+1,m,ν2 =

X m∈Nν1 ,in

Tt,ν1 ,m Ωt−1,s+1,m,ν2 .

Therefore, Ω∗t,s+1,ν1 ,ν2 − Ωt,s+1,ν1 ,ν2 op = ≤

X m∈Nν1 ,in

∗ Tt,ν Ω∗t−1,s+1,m,ν2 − 1 ,m

X m∈Nν1 ,in

Tt,ν1 ,m Ωt−1,s+1,m,ν2

∗ Tt,ν − Tt,ν1 ,m op Ω∗t−1,s+1,m,ν2 op + Ω∗t−1,s+1,m,ν2 − Ωt−1,s+1,m,ν2 op Tt,ν1 ,m op . 1 ,m

By Assumption 3.5, we know that with probability at least 1 − ∗ Tt,ν − Tt,ν1 ,m op ≤ γ1 , 1 ,m

By the induction hypothesis, with probability at least 1 − δ Ω∗t−1,s+1,m,ν2 − Ωt−1,s+1,m,ν2 op ≤ γ1 51

µ∈Nν1 ,in δ ≥ 1 − Nmax δ,

∀m ∈ Nν1 ,in . Pk+1

j j=2 Nmax ,

k X j=1

j Nmax ,

∀m ∈ Nν1 ,in .

Pk+1 j Combining the above arguments, we see that with probability at least 1 − δ j=1 Nmax , we have   k X X j γ1 Ω∗t−1,s+1,m,ν + γ1 Nmax ∥Tt,ν1 ,m ∥op  Ω∗t,s+1,ν1 ,ν2 − Ωt,s+1,ν1 ,ν2 op ≤ 2 op j=1

m∈Nν1 ,in

 (i)

≤ |Nν1 ,in | γ1 + γ1

k X j=1

 j  ε ≤ γ1 Nmax

k+1 X

Nmax .

j=1

In the above display, (i) follows from the fact that Ω∗t−1,s+1,m,ν2 op ≤ 1,

∥Tt,ν1 ,m ∥op ≤ 1,

which can be verified by straightforward calculation. This completes the induction step and the proof of Eq. (43).

Technical lemmas for the M-estimation setting

C.1

∗ Consistency and asymptotic normality of βbt,µ

∗ is consistent and asymptotically normal under Assumptions 4.1-4.6. By In this section, we prove that βbt,µ definition, nt,ν→µ

∗ βbt,µ = arg min P

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

β∈B

L(β, φ(β∗ , εt,ν→µ,i )).

i=1

We first prove consistency. If B is compact, then under Assumption 4.5, as the sample sizes tend to infinity we have

Consistency. sup

β1 ,β2 ∈B

nt,ν→µ

1 P

ν∈Nµ,in nt,ν→µ

ν∈Nµ,in

i=1

L(β1 , φ(β2 , εt,ν→µ,i )) − E[L(β1 , φ(β2 , ε))] → 0.

Therefore, nt,ν→µ

∗ E[L(βbt,µ , φ(β∗ , ε))] ≤ P

ν∈Nµ,in nt,ν→µ

+ sup β1 ,β2 ∈B

≤P

β1 ,β2 ∈B

ν∈Nµ,in

i=1

∗ L(βbt,µ , φ(β∗ , εt,ν→µ,i ))

nt,ν→µ

L(β1 , φ(β2 , εt,ν→µ,i )) − E[L(β1 , φ(β2 , ε))]

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

i=1

nt,ν→µ

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

+ sup

L(β∗ , φ(β∗ , εt,ν→µ,i ))

i=1 nt,ν→µ

L(β1 , φ(β2 , εt,ν→µ,i )) − E[L(β1 , φ(β2 , ε))]

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

≤ E[L(β∗ , φ(β∗ , ε)] + 2 sup

β1 ,β2 ∈B

i=1

nt,ν→µ

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

i=1

L(β1 , φ(β2 , εt,ν→µ,i )) − E[L(β1 , φ(β2 , ε))]

= E[L(β∗ , φ(β∗ , ε)] + oP (1). By Assumption 4.3, β∗ is the unique minimizer of β 7→ E[L(β, φ(β∗ , ε))], and the mapping β 7→ E[L(β, φ(β∗ , ε))] P ∗ is continuous. Therefore, we conclude that βbt,µ → β∗ . 52

We next prove consistency when L is convex and B is not necessarily compact. Since β∗ ∈ int B, there exists r > 0, such that {x : ∥x − β∗ ∥2 ≤ r} ⊆ B. Under Assumption 4.5, we have nt,ν→µ

sup β1 ,β2 ∈{x:∥x−β∗ ∥2 ≤r}

ν∈Nµ,in

i=1

ν∈Nµ,in nt,ν→µ

L(β1 , φ(β2 , εt,ν→µ,i )) − E[L(β1 , φ(β2 , ε))] → 0.

∗ We next show that for all k ∈ N+ satisfying 1/k < r, P(∥βbt,µ − β∗ ∥2 > 1/k) → 0 as the sample sizes tend to ∗ infinity. observe that if ∥βbt,µ − β∗ ∥2 > 1/k, then by convexity, there exists βe that satisfies ∥βe − β∗ ∥2 = 1/k, such that nt,ν→µ

e φ(β∗ , εt,ν→µ,i )) ≤ P L(β,

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

i=1

nt,ν→µ

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

L(β∗ , φ(β∗ , εt,ν→µ,i )).

i=1

By triangle inequality, e φ(β∗ , ε))] − E[L(β∗ , φ(β∗ , ε))] E[L(β, ≤2

nt,ν→µ

sup β1 ,β2 ∈{x:∥x−β∗ ∥2 ≤r}

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

i=1

L(β1 , φ(β2 , εt,ν→µ,i )) − E[L(β1 , φ(β2 , ε))] = oP (1).

By Assumption 4.3 we see that inf ∥β−β∗ ∥2 =1/k E[L(β, φ(β∗ , ε))] > E[L(β∗ , φ(β∗ , ε))]. As a consequence, the ∗ − β∗ ∥2 > 1/k) → 0 as the sample sizes above event occurs with vanishingly small probability, thus P(∥βbt,µ tend to infinity. This holds for all k ∈ N+ , hence completing the proof of consistency. ∗ Asymptotic normality. We then prove asymptotic normality. By consistency, we know that P(∥βbt,µ − β∗ ∥2 < r) → 1 as the sample sizes tend to infinity. Therefore, with probability 1 − o(1) the following first-order condition holds: nt,ν→µ

ν∈Nµ,in

i=1

∇β L(β, φ(β∗ , εt,ν→µ,i ))

b∗ β=β t,µ

= 0.

Rewriting the expression in integral form, we obtain 1 P

nt,ν→µ Z 1

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

∇2β L(β, φ(β∗ , εt,ν→µ,i ))

b∗ β=(1−t)β∗ +tβ t,µ

∗ (βbt,µ − β∗ )dt

nt,ν→µ

= −P

i=1

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

X i=1

∇β L(β, φ(β∗ , εt,ν→µ,i ))

. β=β∗

By Assumption 4.5, 1

nt,ν→µ Z 1

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

i=1

∇2β L(β, φ(β∗ , εt,ν→µ,i ))

b∗ β=(1−t)β∗ +tβ t,µ

P dt → E ∇2β L(β, φ(β∗ , ε)) β=β . ∗

By the central limit theorem, − qP

nt,ν→µ

ν∈Nµ,in nt,ν→µ ν∈Nµ,in

X i=1

∇β L(β, φ(β∗ , εt,ν→µ,i ))

d → N 0, E[∇β L(β, φ(β∗ , ε))∇β L(β, φ(β∗ , ε))⊤ ] β=β . ∗

The proof is done. 53

β=β∗

C.2

Proof of Lemma A.10

Proof of (a). Using the recursive formula of Theorem 4.1, we obtain that ΣT =

T X t=1

bT,t (JT,t+1 ⊗ Ip )Vt (JT,t+1 ⊗ Ip )⊤ + bT,0 (JT,1 ⊗ Ip )Σ0 (JT,1 ⊗ Ip )⊤ .

Therefore, for any µ ∈ Ml , ΣT,µ,µ =

T X

bT,t

JT,t+1,µ,ν1 JT,t+1,µ,ν2 Vt,ν1 ,ν2 + bT,0

ν1 ,ν2 ∈Ml

t=1

JT,1,µ,ν1 JT,1,µ,ν2 Σ0,ν1 ,ν2 ,

ν1 ,ν2 ∈M

where we use the fact that Vt,ν1 ,ν2 = 0 if ν1 ∈ Mu or ν2 ∈ Mu . Since by definition, Vt,ν1 ,ν2 ⪰ 0 for all ν1 , ν2 ∈ M, and Σ0,ν1 ,ν2 = 0 if ν1 ̸= ν2 , we deduce that ΣT,µ,µ ⪰ =

T X

bT,t

t=1 T X

ν∈Ml

bT,t

t=1 (i)

⪰ (ii)

2 JT,t+1,µ,ν Vt,ν,ν + bT,0

2 JT,1,µ,ν Σ0,ν,ν

ν∈M

X P ν∈Ml

m∈Nν,in pt,m→ν

2 JT,t+1,µ,ν + bT,0

ν∈M

p0,ν

! 2 JT,1,µ,ν

V∗

! P 2 ν∈M JT,1,µ,ν ) P P + bT,0 V∗ ν∈Ml ν∈M p0,ν m∈Nν,in pt,m→ν t=1 ! P T X bT,t ( ν∈Ml JT,t+1,µ,ν )2 bT,0 P P +P V∗ , ν∈Ml ν∈M p0,ν m∈Nν,in pt,m→ν t=1 T X

( bT,t P

2 ν∈Ml JT,t+1,µ,ν )

(

P where in (i) we use Cauchy-Schwarz inequality, and (ii) follows from the fact that ν∈M JT,1,µ,ν = 1. This completes the proof of part (a). Proof of (b). Following our calculations in the proof of part (a), we have for all µ ∈ Ml : P X ( ν∈Ml JT,t+1,µ,ν )2 P V∗ . JT,t+1,µ,ν1 JT,t+1,µ,ν2 Vt,ν1 ,ν2 ⪯ inf ν∈Ml m∈Nν,in pt,m→ν ν1 ,ν2 ∈Ml

Similarly, X

JT,1,µ,ν1 JT,1,µ,ν2 Σ0,ν1 ,ν2 =

ν1 ,ν2 ∈M

X ν∈M

2 JT,1,µ,ν Σ0,ν,ν ⪯

P ( ν∈M JT,1,µ,ν )2 1 V∗ = V∗ . inf ν∈M p0,ν inf ν∈M p0,ν

It finally follows that ΣT,µ,µ ⪯

! P T X bT,t ( ν∈Ml JT,t+1,µ,ν )2 bT,0 P + V∗ , inf ν∈Ml m∈Nν,in pt,m→ν inf ν∈M p0,ν t=1

completing the proof of part (b).

C.3

Proof of Lemma A.11

We verify Assumptions 4.3-4.6 one by one below. Verifying Assumption 4.3. By definition of L, we know that the twice differentiability of the mapping β 7→ L(β, z) follows from the twice differentiability of A. Further, we have E[L(β1 , φ(β2 , ε))] = E −A′ (β2⊤ x)β1⊤ x + A(β1⊤ x) , 54

whose continuity follows from the continuous differentiability of A′ and dominated convergence theorem. Next, we have for any fixed β: E ∇β1 L(β, φ(β ⊤ x, ε)) = E −φ(β ⊤ x, ε)x + A′ (β ⊤ x)x = E x(A′ (β ⊤ x) − E[y|x]) = 0, since E[y|x] = A′ (β ⊤ x). Combining this with the convexity of L verifies that β is the unique minimizer of β ′ 7→ E[L(β ′ , φ(β, ε))]. Further, the population Hessian is positive definite since A is strictly convex. Finally, the integrability conditions in Assumption 4.3 can be directly verified. Verifying Assumption 4.4. from the convexity of A.

This is easily verified because B = Rd , and the convexity of L follows directly

Verifying point 1 of Assumption 4.5. We first show that the function class (x, ε) 7→ −φ(β2⊤ x, ε) · β1⊤ x + A(β1⊤ x) (β1 , β2 ) ∈ Ω2

(71)

is µε -Glivenko-Cantelli. Since E[∥x∥2 ] < ∞, the function class {x 7→ β1⊤ x|β1 ∈ Ω} is µε -Glivenko-Cantelli. By our assumption, {x 7→ −φ(β2⊤ x, ε)|β2 ∈ Ω} and {x 7→ A(β1⊤ x)|β1 ∈ Ω} are also µε -Glivenko-Cantelli. Let rΩ > 0 be such that Ω ⊂ B(0, rΩ ), then h i h i h i E sup |L (β1 , φ(β2 , ε))| ≤ rΩ E sup ∥x∥2 |φ(β ⊤ x, ε)| + E sup |A(β ⊤ x)| < ∞, (β1 ,β2 )∈Ω2

β∈Ω

meaning that the function class (71) has an L1 -integrable envelope. Applying Lemma C.1 shows that the function class (71) is µε -Glivenko-Cantelli. Finally, note that ∇2β1 L β1 , φ(β2⊤ x, ε) = A′′ (β1⊤ x)xx⊤ , which is assumed to be µε -Glivenko-Cantelli over compact sets. This completes the verification of point 1 of Assumption 4.5. Verifying point 2 of Assumption 4.5.

Note that by definition: ∇β1 L β1 , φ(β2⊤ x, ε) = x A′ (β1⊤ x) − φ(β2⊤ x, ε) .

Since E[A′ (β1⊤ x)2 xx⊤ ] < ∞, we know that the function class {x 7→ A′ (β1⊤ x)x|β2 ∈ Ω} is µε -Donsker. Further, we have assumed that {(x, ε) 7→ φ(β2⊤ x, ε)x|β2 ∈ Ω} is µε -Donsker, applying Lemma C.2 completes the verification of point 2 of Assumption 4.5. Verifying Assumption 4.6.

By direct calculation, we know that ∇2β1 L(β1 , φ(β2 , ε)) = A′′ (β1⊤ x)xx⊤ .

Therefore, Assumption 4.6 directly follows from point (iv) in the statement of Lemma A.11 and the dominated convergence theorem.

C.4

Auxiliary lemmas from empirical process theory

We need the following auxiliary results for the preservation of Glivenko-Cantelli and Donsker classes under composition of functions. Lemma C.1 (Preservation of Glivenko-Cantelli classes, Theorem 3 in Van Der Vaart and Wellner (2000)). Let F1 , . . . , Fk be P -Glivenko-Cantelli classes of measurable functions, and let φ : Rk → R be a continuous function. Define the class H as the image of these classes under φ: H = {x 7→ φ(f1 (x), . . . , fk (x)) : fj ∈ Fj for j = 1, . . . , k} . Then, the class H is P -Glivenko-Cantelli, provided that H has an integrable envelope H (i.e., H ∈ L1 (P ) such that |h(x)| ≤ H(x) for all h ∈ H). 55

Lemma C.2 (Preservation of Donsker classes, Theorem 2.10.6 in Van Der Vaart and Wellner (1996)). Let F1 , . . . , Fk be P -Donsker classes of measurable functions. Let φ : Rk → R be an L-Lipschitz function for some constant L > 0. Define the class H as: H = {x 7→ φ(f1 (x), . . . , fk (x)) : fj ∈ Fj } . If there exists at least one tuple (f1 , . . . , fk ) such that φ(f1 , . . . , fk ) is square-integrable (i.e., belongs to L2 (P )), then the class H is a P -Donsker class. The next lemma establishes the P -Donsker property for a wide range of function classes: Lemma C.3. Let K ⊆ Rp be a compact set, ε ∼ Unif[0, 1] be a random variable that is independent of the random vector x ∈ Rp , and g : Rp → R be a measurable function. Let {fk }∞ k=1 be a sequence of non-decreasing functions taking values in [0, 1]. We assume the following statements hold: (i) There exist a non-decreasing function L(t) > 0 and a constant c1 < 2, such that for all k ≥ 1, |fk′ (t)| ≤ k −c1 L(t) for all t ∈ R; (ii) There exist a constant c2 > 2 and a non-negative non-decreasing function f , such that fk (t) ≤ k −c2 f (t) for all k ≥ 1; (iii) For any R > 0: E g(x)2 ∥x∥2 L R∥x∥2 < ∞, Then, the function class HK =

E g(x)2 f R∥x∥2 < ∞.

∞ n o X (x, ε) 7→ g(x)1ε≤fk (β ⊤ x) : β ∈ K k=1

is P -Donsker, where P represents the joint distribution of (x, ε). Proof. We invoke the Donsker theorem (cf. Section 2.5.2 of Van Der Vaart and Wellner (1996)) to show that HK is a P -Donsker class. Without loss of generality, we may assume that g is nonnegative. Otherwise, we decompose g as g = g+ − g− , and note that the function classes obtained by replacing g with either g+ or g− in the definition of HK are both P -Donsker. It then follows from Lemma C.2 that HK is also P -Donsker. To prove the lemma, it suffices to show Z 1q log N[ ] (γ, HK ∪ {0}, L2 (P )) dγ < ∞, (72) 0

where N[ ] denotes the bracketing number. We prove Eq. (72) by upper bounding the bracketing number N[ ] (γ, HK ∪ {0}, L2 (P )) for all γ > 0. For δ > 0 (to be determined later), let N (δ, K) = {β1 , · · · , βM } be a δ-net of K (under ℓ2 -norm). Then, we know that p CK M = |N (δ, K)| ≤ , δ where CK > 0 is a constant that depends only on K. We next upper bound the bracketing number by constructing brackets. For j = 1, · · · , M , define lj (x, ε) =

∞ X

g(x)1ε≤fk (βj⊤ x−δ∥x∥2 ) ,

uj (x, ε) =

k=1

∞ X

g(x)1ε≤fk (βj⊤ x+δ∥x∥2 ) .

k=1

We next show that HK ⊂ ∪M j=1 [lj , uj ]. To this end, note that for any β ∈ K, there exists j ∈ [M ] such that ∥β − βj ∥2 ≤ δ, which implies βj⊤ x − δ∥x∥2 ≤ β ⊤ x ≤ βj⊤ x + δ∥x∥2 . The claim then follows from our assumption that fk is non-decreasing and g is non-negative. Next, we upper bound ∥uj − lj ∥L2 (P ) . Observe that for any fixed positive integer B, we have ∥uj − lj ∥L2 (P ) 56

≤

B X k=1

g(x) 1ε≤fk (βj⊤ x+δ∥x∥2 ) − 1ε≤fk (βj⊤ x−δ∥x∥2 )

L2 (P )

∞ X

g(x)1ε≤fk (βj⊤ x+δ∥x∥2 )

k=B+1

L2 (P )

X B ∞ X 1/2 1/2 2 ⊤ ⊤ 2 ⊤ ≤ E g(x) fk (βj x + δ∥x∥2 ) − fk (βj x − δ∥x∥2 + 2 kE g(x) fk (βj x + δ∥x∥2 ) k=1

≤

k=B+1

B √ X k=1

≤ CB

X ∞ q 1/2 −c2 +1 2 ⊤ ⊤ −c 2 1 2k δ E g(x) ∥x∥2 L(βj x + δ∥x∥2 ) + 2 k E g(x) f (βj x + δ∥x∥2 ) k=B+1

1−c1 /2

q δE g(x)2 ∥x∥2 L (rK + δ)∥x∥2 + CB 1−c2 /2 E g(x)2 f (rK + δ)∥x∥2 ,

where C > 0 is a constant that depends only on c1 and c2 , and rK > 0 is a constant that depends only on ′ ′ K. Note that setting B ≍ γ −c1 and δ ≍ γ c2 for some c′1 , c′2 > 0 gives ∥uj − lj ∥L2 (P ) ≤ γ. We finally deduce that 1 1 log N[ ] γ, HK ∪ {0}, L2 (P ) ≲ p log ≲ p log . δ γ This establishes Eq. (72), and hence completes the proof of the lemma.

Record · ID 204770 · SHA-256 b481b1f68274ce56
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.