summaries: The paper considers the problem of robustly learning fixed structure Bayesian networks in nearly-linear time. Previous work by Cheng et al. gives a runtime of O(Nd^2/eps). The paper improves this to O(Nd). The algorithm works by directly relating the problem to robust mean estimation, and then leveraging the algorithm of Dong et al. for robust mean estimation which works in nearly-linear time. The authors have to modify the runtime analysis of the algorithm of Dong et al. to work in time linear in the sparsity, rather than dimension. The paper studies the problem of robust learning of fixed-structure Bayesian networks under the eps-adversarial corruptions model. Fixed-structure means a known structure of the underlying Bayesian network. Robust learning is an important area of research and this particular question has been studied in prior work. The main contribution of this work is in improving the running time of the algorithm. On a d-node Bayes net, let m denote the total number of parental configurations possible. Prior work of Cheng et al showed a robust learning algorithm using O(m/eps^2) samples and runs in time O(md^2/eps^2). This paper studies the problem of learning Bayes nets using adversarially corrupted data. The model is that $N$ samples are made from a Bayes net on d nodes, out of which an unknown $\varepsilon$ fraction are changed arbitrarily. The structure of the Bayes net is already given, but it remains to learn the probability distribution.
abstractText: We study the problem of learning Bayesian networks where an -fraction of the samples are adversarially corrupted. We focus on the fully-observable case where the underlying graph structure is known. In this work, we present the first nearlylinear time algorithm for this problem with a dimension-independent error guarantee. Previous robust algorithms with comparable error guarantees are slower by at least a factor of (d/ ), where d is the number of variables in the Bayesian network and is the fraction of corrupted samples. Our algorithm and analysis are considerably simpler than those in previous work. We achieve this by establishing a direct connection between robust learning of Bayesian networks and robust mean estimation. As a subroutine in our algorithm, we develop a robust mean estimation algorithm whose runtime is nearly-linear in the number of nonzeros in the input samples, which may be of independent interest.
authors.name: NEARLY-LINEAR TIME
authors.name: Yu Cheng
authors.name: Honghao Lin
references.authors: P. Abbeel D. Koller A.Y. Ng
references.title: Learning factor graphs in polynomial time and sample complexity
references.venue: J. Mach. Learn. Res.,
references.authors: E. Amaldi V. Kann
references.title: The complexity and approximability of finding maximum feasible subsystems of linear relations
references.venue: Theoretical Computer Science,
references.authors: A. Anandkumar D.J. Hsu F. Huang S. Kakade
references.title: Learning mixtures of tree graphical models
references.venue: In Proc. 26th Advances in Neural Information Processing Systems (NeurIPS),
references.authors: S. Balakrishnan S.S. Du J. Li A. Singh
references.title: Computationally efficient robust sparse estimation in high dimensions
references.venue: In Proc. 30th Conference on Learning Theory (COLT),
references.authors: G. Bresler
references.title: Efficiently learning Ising models on arbitrary graphs
references.venue: In Proc. 47th Annual ACM Symposium on Theory of Computing (STOC),
references.authors: G. Bresler E. Mossel A. Sly
references.title: Reconstruction of Markov random fields from samples: Some observations and algorithms
references.venue: SIAM J. Comput.,
references.authors: G. Bresler D. Gamarnik D. Shah
references.title: Structure learning of antiferromagnetic Ising models
references.venue: In Proc. 28th Advances in Neural Information Processing Systems (NeurIPS),
references.authors: C.L. Canonne I. Diakonikolas D.M. Kane A. Stewart
references.title: Testing Bayesian networks
references.venue: In Proc. 30th Conference on Learning Theory (COLT),
references.authors: T.M. Chan
references.title: An optimal randomized algorithm for maximum Tukey depth
references.venue: In Proc. 15th ACM-SIAM Symposium on Discrete Algorithms (SODA),
references.authors: M. Charikar J. Steinhardt G. Valiant
references.title: Learning from untrusted data
references.venue: In Proc. 49th Annual ACM Symposium on Theory of Computing (STOC),
references.authors: M. Chen C. Gao Z. Ren
references.title: Robust covariance and scatter matrix estimation under Huber’s contamination model
references.venue: The Annals of Statistics,
references.authors: Y. Cheng I. Diakonikolas D.M. Kane A. Stewart
references.title: Robust learning of fixed-structure Bayesian networks
references.venue: In Proc. 32nd Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Y. Cheng I. Diakonikolas R. Ge
references.title: High-dimensional robust mean estimation in nearly-linear time
references.venue: In Proc. 30th ACM-SIAM Symposium on Discrete Algorithms (SODA),
references.authors: Y. Cheng I. Diakonikolas R. Ge D.P. Woodruff
references.title: Faster algorithms for high-dimensional robust covariance estimation
references.venue: In Proc. 32nd Conference on Learning Theory (COLT),
references.authors: Y. Cheng I. Diakonikolas R. Ge M. Soltanolkotabi
references.title: High-dimensional robust mean estimation via gradient descent
references.venue: In Proc. 37th International Conference on Machine Learning (ICML),
references.authors: Y. Cherapanamjeri E. Aras N. Tripuraneni M.I. Jordan N. Flammarion P.L. Bartlett
references.title: Optimal robust linear regression in nearly linear time
references.venue: arXiv preprint arXiv:2007.08137,
references.authors: Y. Cherapanamjeri S. Mohanty M. Yau
references.title: List decodable mean estimation in nearly linear time
references.venue: In Proc. 61st IEEE Symposium on Foundations of Computer Science (FOCS),
references.authors: C. Chow C. Liu
references.title: Approximating discrete probability distributions with dependence trees
references.venue: IEEE Trans. Inf. Theor.,
references.authors: K.L. Clarkson D. Eppstein G.L. Miller C. Sturtivant S.-H. Teng
references.title: Approximating center points with iterated Radon points
references.venue: In Proc. 9th Annual Symposium on Computational Geometry (SoCG), pp
references.authors: R. Daly Q. Shen S. Aitken
references.title: Learning Bayesian networks: approaches and issues
references.venue: The Knowledge Engineering Review,
references.authors: S. Dasgupta
references.title: The sample complexity of learning fixed-structure Bayesian networks
references.venue: Machine Learning,
references.authors: J. Depersin G. Lecué
references.title: Robust subgaussian estimation of a mean vector in nearly linear time
references.venue: arXiv preprint arXiv:1906.03058,
references.authors: L. Devroye L. Györfi
references.title: Nonparametric Density Estimation: The L1 View
references.authors: I. Diakonikolas D.M. Kane
references.title: Recent advances in algorithmic high-dimensional robust statistics
references.venue: arXiv preprint arXiv:1911.05911,
references.authors: I. Diakonikolas G. Kamath D.M. Kane J. Li A. Moitra A. Stewart
references.title: Being robust (in high dimensions) can be practical
references.venue: In Proc. 34th International Conference on Machine Learning (ICML),
references.authors: I. Diakonikolas D.M. Kane A. Stewart
references.title: Statistical query lower bounds for robust estimation of high-dimensional Gaussians and Gaussian mixtures
references.venue: In Proc. 58th IEEE Symposium on Foundations of Computer Science (FOCS),
references.authors: I. Diakonikolas D.M. Kane A. Stewart
references.title: List-decodable robust mean estimation and learning mixtures of spherical Gaussians
references.venue: In Proc. 50th Annual ACM Symposium on Theory of Computing (STOC),
references.authors: I. Diakonikolas G. Kamath D.M. Kane J. Li A. Moitra A. Stewart
references.title: Robust estimators in high dimensions without the computational intractability
references.venue: SIAM Journal on Computing,
references.authors: I. Diakonikolas G. Kamath D.M. Kane J. Li J. Steinhardt A. Stewart
references.title: SEVER: A robust meta-algorithm for stochastic optimization
references.venue: In Proc. 36th International Conference on Machine Learning (ICML),
references.authors: I. Diakonikolas W. Kong A. Stewart
references.title: Efficient algorithms and lower bounds for robust linear regression
references.venue: In Proc. 30th ACM-SIAM Symposium on Discrete Algorithms (SODA),
references.authors: I. Diakonikolas D.M. Kane D. Kongsgaard J. Li K. Tian
references.title: List-decodable mean estimation in nearly-PCA time
references.venue: arXiv preprint arXiv:2011.09973,
references.authors: Y. Dong S.B. Hopkins J. Li
references.title: Quantum entropy scoring for fast robust mean estimation and improved outlier detection
references.venue: In Proc. 33rd Advances in Neural Information Processing Systems (NeurIPS),
references.authors: S.B. Hopkins J. Li
references.title: Mixture models, robustness, and sum of squares proofs
references.venue: In Proc. 50th Annual ACM Symposium on Theory of Computing (STOC),
references.authors: P.J. Huber
references.title: Robust estimation of a location parameter
references.venue: Ann. Math. Statist., 35(1):73–101,
references.authors: F.V. Jensen T.D. Nielsen
references.title: Bayesian Networks and Decision Graphs
references.authors: D.S. Johnson F.P. Preparata
references.title: The densest hemisphere problem
references.venue: Theoretical Computer Science,
references.authors: A. Klivans P. Kothari R. Meka
references.title: Efficient algorithms for outlier-robust regression
references.venue: In Proc. 31st Conference on Learning Theory (COLT),
references.authors: D. Koller N. Friedman
references.title: Probabilistic Graphical Models: Principles and Techniques - Adaptive Computation and Machine Learning
references.authors: P.K. Kothari J. Steinhardt D. Steurer
references.title: Robust moment estimation and improved clustering via sum of squares
references.venue: In Proc. 50th Annual ACM Symposium on Theory of Computing (STOC),
references.authors: K.A. Lai A.B. Rao S. Vempala
references.title: Agnostic estimation of mean and covariance
references.venue: In Proc. 57th IEEE Symposium on Foundations of Computer Science (FOCS),
references.authors: J. Li G. Ye
references.title: Robust Gaussian covariance estimation in nearly-matrix multiplication time
references.venue: In Proc. 34th Advances in Neural Information Processing Systems (NeurIPS),
references.authors: L. Liu Y. Shen T. Li C. Caramanis
references.title: High dimensional robust sparse regression
references.venue: In Proc. 23rd International Conference on Artificial Intelligence and Statistics (AISTATS),
references.authors: P.L. Loh M.J. Wainwright
references.title: Structure estimation for discrete graphical models: Generalized covariance matrices and their inverses
references.venue: In Proc. 26th Advances in Neural Information Processing Systems (NeurIPS),
references.authors: G.L. Miller D. Sheehy
references.title: Approximate centerpoints with proofs
references.venue: Comput. Geom.,
references.authors: R.E. Neapolitan
references.title: Learning Bayesian Networks
references.venue: Prentice-Hall, Inc.,
references.authors: A. Prasad A.S. Suggala S. Balakrishnan P. Ravikumar
references.title: Robust estimation via robust gradient estimation
references.venue: Journal of the Royal Statistical Society: Series B (Statistical Methodology),
references.authors: N.P. Santhanam M.J. Wainwright
references.title: Information-theoretic limits of selecting binary graphical models in high dimensions
references.venue: IEEE Trans. Information Theory,
references.authors: J. Steinhardt M. Charikar G. Valiant
references.title: Resilience: A criterion for learning in the presence of arbitrary outliers
references.venue: In In Proc. 9th Innovations in Theoretical Computer Science Conference (ITCS),
references.authors: J.W. Tukey
references.title: Mathematics and picturing of data
references.venue: In Proceedings of ICM,
references.authors: R. Vershynin
references.title: Introduction to the non-asymptotic analysis of random matrices
references.venue: arXiv preprint arXiv:1011.3027,
references.authors: M.J. Wainwright M.I. Jordan
references.title: Graphical models, exponential families, and variational inference
references.venue: Found. Trends Mach. Learn.,
references.authors: M.J. Wainwright P. Ravikumar J.D. Lafferty
references.title: High-dimensional graphical model selection using `1-regularized logistic regression
references.venue: In Proc. 20th Advances in Neural Information Processing Systems (NeurIPS),
references.authors: B. Zhu J. Jiao J. Steinhardt
references.title: Robust estimation via generalized quasi-gradients
references.venue: arXiv preprint arXiv:2005.14073,
references.authors: Dong
references.title: section, we give a robust mean estimation algorithm that runs in nearly input-sparsity time
references.venue: We build on the following lemma,
references.authors: Dong
references.title: 2019)). Given an -corrupted version of an ( , β, γ)-stable set S of N samples w.r.t. a d-dimensional distribution with mean μX
references.authors: Dong
references.title: This is because we use a more concise stability condition than the one in Dong et al. (2019). We will show that Lemma C.1 is equivalent to the version in Dong et al. (2019) in Appendix C.1
references.venue: The computational bottleneck of the algorithm in Dong et al
references.authors: Dong
references.title: We say a set of points S is ( ,γ1,γ2,β1,β2)-good with respect to a distribution D with true mean μ if the following two properties hold
references.authors: Then Dong
references.title: 2019) showed the following result. Lemma C.6. Let D be a distribution on R with unknown mean μ. Let 0 < < 0 for some universal constant 0. Let S be a set of N samples with S = Sg ∪ Sb\Sr
sections.heading: 1 INTRODUCTION
sections.text: Probabilistic graphical models (Koller & Friedman, 2009) offer an elegant and succinct way to represent structured high-dimensional distributions. The problem of inference and learning in probabilistic graphical models is an important problem that arises in many disciplines (see Wainwright & Jordan (2008) and the references therein), which has been studied extensively during the past decades (see, e.g., Chow & Liu (1968); Dasgupta (1997); Abbeel et al. (2006); Wainwright et al. (2006); Anandkumar et al. (2012); Santhanam & Wainwright (2012); Loh & Wainwright (2012); Bresler et al. (2013; 2014); Bresler (2015)).
  Bayesian networks (Jensen & Nielsen, 2007) are an important family of probabilistic graphical models that represent conditional dependence by a directed graph (see Section 2 for a formal definition). In this paper, we study the problem of learning Bayesian networks where an -fraction of the samples are adversarially corrupted. We focus on the simplest setting: all variables are binary and observable, and the structure of the Bayesian network is given to the algorithm.
  Formally, we work with the following corruption model:
  Definition 1.1 ( -Corrupted Set of Samples). Given 0 < < 1/2 and a distribution family P on Rd, the algorithm first specifies the number of samples N , and N samples X1, X2, . . . , XN are drawn from some unknown P ∈ P . The adversary inspects the samples, the ground-truth distribution P , and the algorithm, and then replaces N samples with arbitrary points. The set of N points is given to the algorithm as input. We say that a set of samples is -corrupted if it is generated by this process.
  This is a strong corruption model which generalizes many existing models. In particular, it is stronger than Huber’s contamination model (Huber, 1964), because we allow the adversary to add bad samples and remove good samples, and he can do so adaptively.
  Our goal is to design robust algorithms for learning Bayesian networks with dimension-independent error. More specifically, given as input an -corrupted set of samples drawn from some groundtruth Bayesian network P and the graph structure of P , we want the algorithm to output a Bayesian network Q, such that the total variation distance between P and Q is upper bounded by a function that depends only on (the fraction of corruption) but not d (the number of variables in P ).
  In the fully-observable fixed-structure setting, the problem is straightforward when there is no corruption. We know that the empirical estimator (which computes the empirical conditional probabilities) is sample efficient and runs in linear time (Dasgupta, 1997).
  It turns out that the problem becomes much more challenging when there is corruption. Even for robust learning of binary product distributions (i.e., a Bayesian network with an empty dependency graph), the first computational efficient algorithms with dimension-independent error was only discovered in (Diakonikolas et al., 2019a). Subsequently, (Cheng et al., 2018) gave the first polynomialtime algorithms for robust learning of fixed-structured Bayesian networks. The main drawback of the algorithm in (Cheng et al., 2018) is that it runs in time Ω(Nd2/ ), which is slower by at least a factor of (d/ ) compared to the fastest non-robust estimator.
  Motivated by this gap in the running time, in this work we want to resolve the following question:
  Can we design a robust algorithm for learning Bayesian networks in the fixedstructure fully-observable setting that runs in nearly-linear time?
sections.heading: 1.1 OUR RESULTS AND CONTRIBUTIONS
sections.text: We resolve this question affirmatively by proving Theorem 1.2. We say a Bayesian network is cbalanced if all its conditional probabilities are between c and 1 − c. For the ground-truth Bayesian network P , let m be the size of its conditional probability table and α be its minimum parental configuration probability (see Section 2 for formal definitions).
  Theorem 1.2 (informal statement). Consider an -corrupted set of N = Ω̃(m/ 2) samples drawn from a d-dimensional Bayesian network P . Suppose P is c-balanced and has minimum parental configuration probability α, where both c and α are universal constants. We can compute a Bayesian network Q in time Õ(Nd) such that dTV (P,Q) ≤ √ ln(1/ ). 1
  For simplicity, we stated our result in the very special case where both c and α are Ω(1). Our approach works for general values of α and c, where our error guarantee degrades gracefully as α and c gets smaller. A formal version of Theorem 1.2 is given as Theorem 4.1 in Section 4.
  Our algorithm has optimal error guarantee, sample complexity, and running time (up to logarithmic factors). There is an information-theoretic lower bound of Ω( ) on the error guarantee, which holds even for Bayesian networks with only one variable. A sample complexity lower bound of Ω(m/ 2) holds even without corruption (see, e.g., (Canonne et al., 2017)).
  Our Contributions. We establish a novel connection between robust learning of Bayesian networks and robust mean estimation. At a high level, we show that one can essentially reduce the former to the latter. This allows us to take advantage of the recent (and future) advances in robust mean estimation and apply the algorithms almost directly to obtain new algorithms for learning Bayesian networks.
  Our algorithm and analysis are considerably simpler than those in previous work. For simplicity, consider learning binary product distributions as an example. Cheng et al. (2018) tried to remove samples to make the empirical covariance matrix closer to a diagonal matrix (since the true covariance matrix is diagonal because each coordinate is independent). They used a “filtering” approach which requires proving specific tail bounds on the samples. In contrast, we show that it suffices to use any robust mean estimation algorithms which minimize the spectral norm of the empirical covariance matrix (regardless of whether it is close to being diagonal or not).
  As a subroutine in our approach, we develop the first robust mean estimation algorithm that runs in nearly input-sparsity time (i.e., in time nearly linear in the total number of nonzero entries in the input), which may be of independent interest. The main computation bottleneck of current nearlylinear time robust mean estimation algorithms (Cheng et al., 2019a; Depersin & Lecué, 2019; Dong et al., 2019) is running matrix multiplication weight update with the Johnson-Lindenstrauss lemma, which we show can be done in nearly input-sparsity time.
  1Throughout the paper, we use Õ(f) to denote O(f polylog(f)).
sections.heading: 1.2 RELATED WORK
sections.text: Bayesian Networks. Probabilistic graphical models (Koller & Friedman, 2009) provide an appealing and unifying formalism to succinctly represent structured high-dimensional distributions. The general problem of inference in graphical models is of fundamental importance and arises in many applications across several scientific disciplines (see Wainwright & Jordan (2008) and references therein). The problem of learning graphical models from data (Neapolitan, 2003; Daly et al., 2011) has many variants: (i) the family of graphical models (e.g., directed, undirected), (ii) whether the data is fully or partially observable, and (iii) whether the graph structure is known or not. This learning problem has been studied extensively (see, e.g., Chow & Liu (1968); Dasgupta (1997); Abbeel et al. (2006); Wainwright et al. (2006); Anandkumar et al. (2012); Santhanam & Wainwright (2012); Loh & Wainwright (2012); Bresler et al. (2013; 2014); Bresler (2015)), resulting in a beautiful theory and a collection of algorithms in various settings.
  Robust Statistics. Learning in the presence of outliers has been studied since the 1960s (Huber, 1964). For the most basic problem of robust mean estimation, it is well-known that the empirical median works in one dimension. However, most natural generalizations of the median to high dimensions (e.g., coordinate-wise median, geometric median) would incur an error of Ω( √ d), even in the infinite sample regime (see, e.g., Diakonikolas et al. (2019a); Lai et al. (2016)). After decades of work, sample-efficient robust estimators have been discovered (e.g., the Tukey median (Tukey, 1975; Devroye & Györfi, 1985; Chen et al., 2018)). However, the Tukey median is NP-Hard to compute in the worse case (Johnson & Preparata, 1978; Amaldi & Kann, 1995) and many heuristics for approximating it perform poorly as the dimension scales (Clarkson et al., 1993; Chan, 2004; Miller & Sheehy, 2010).
  Computational Efficient Robust Estimators. Recent work (Diakonikolas et al., 2019a; Lai et al., 2016) gave the first polynomial-time algorithms several high-dimensional unsupervised learning tasks (e.g., mean and covariance estimation) with dimension-independent error guarantees. After the dissemination of (Diakonikolas et al., 2019a; Lai et al., 2016), algorithmic high-dimensional robust statistics has attracted a lot of recent attention and there has been a flurry of research that obtained polynomial-time robust algorithms for a wide range of machine learning and statistical tasks (see, e.g., Balakrishnan et al. (2017); Charikar et al. (2017); Diakonikolas et al. (2017a;b); Steinhardt et al. (2018); Diakonikolas et al. (2018); Hopkins & Li (2018); Kothari et al. (2018); Prasad et al. (2020); Diakonikolas et al. (2019b); Klivans et al. (2018); Diakonikolas et al. (2019c); Liu et al. (2020); Cheng et al. (2020); Zhu et al. (2020)). In particular, the most relevant prior work is (Cheng et al., 2018), which gave the first polynomial-time algorithms for robust learning of fixed-structure Bayesian networks.
  Faster Robust Estimators. While recent work gave polynomial-time robust algorithms for many tasks, these algorithms are often significantly slower than the fastest non-robust ones (e.g., sample average for mean estimation). Cheng et al. (2019a) gave the first nearly-linear time algorithm for robust mean estimation and initiated the research direction of designing robust estimators that are as efficient as their non-robust counterparts. Since then, there have been several works that develop faster robust algorithms for various learning and statistical tasks, including robust mean estimation for heavy-tailed distributions Dong et al. (2019); Depersin & Lecué (2019), robust covariance estimation Cheng et al. (2019b); Li & Ye (2020), robust linear regression Cherapanamjeri et al. (2020a), and list-decodable mean estimation Cherapanamjeri et al. (2020b); Diakonikolas et al. (2020).
  Organization. In Section 2, we define our notations and provide some background on robust learning of Bayesian networks and robust mean estimation. In Section 3, we give an overview of our approach and highlight some of our key technical results. In Section 4, we present our algorithm for robust learning of Bayesian networks and prove our main result.
sections.heading: 2 PRELIMINARIES
sections.text: Bayesian Networks. Fix a d-node directed acyclic graph H whose nodes are labelled [d] = {1, 2, . . . , d} in topological order (every edge goes from a node with smaller index to one with larger index). Let Parents(i) be the parents of node i in H . A probability distribution P on
  {0, 1}d is a Bayesian network (or Bayes net) with graph H if, for each i ∈ [d], we have that PrX∼P [Xi = 1 | X1, . . . , Xi−1] depends only on the values Xj where j ∈ Parents(i).
  Conditional Probability Table. Let P be a Bayesian network with graph H . Let Γ = {(i, a) : i ∈ [d], a ∈ {0, 1}|Parents(i)|} be the set of all possible parental configurations. Let m = |Γ|. For (i, a) ∈ Γ, the parental configuration Πi,a is defined to be the event that X(Parents(i)) = a. The conditional probability table p ∈ [0, 1]m of P is given by pi,a = PrX∼P [X(i) = 1 | Πi,a] . In this paper, we often index p as anm-dimensional vector. We use the notation pk and the associated events Πk, where each k ∈ [m] stands for an (i, a) ∈ Γ lexicographically ordered.
  Notations. For a vector v, let ‖v‖2 and ‖v‖∞ be the `2 and `∞ norm of v respectively. We write√ v and 1/v for the entrywise square root and entrywise inverse of a vector v respectively. For two vectors x and y, we write x>y for their inner product, and x ◦ y for their entrywise product. We use I to denote the identity matrix. For a matrix M , let Mi be the i-th column of M , and let ‖M‖2 be the spectral norm of M . For a vector v ∈ Rn, let diag(v) ∈ Rn×n denote a diagonal matrix with v on the diagonal.
  Throughout this paper, we use P to denote the ground-truth Bayesian network. We use d for the dimension (i.e., the number of nodes) of P , N for the number of samples, for the fraction of corrupted samples, and m = ∑d i=1 2
  |Parents(i)| for the size of the conditional probability table of P . We use p ∈ Rm to denote the (unknown) ground-truth conditional probabilities of P , and q ∈ Rm for our current guess of p.
  Let G? be the original set of N uncorrupted samples drawn from P . After the adversary corrupts an -fraction of G?, let G ⊆ G? be the remaining set of good samples, and B be the set of bad samples added by the adversary. The set of samples S = G ∪ B is given to the algorithm as input. Let X ∈ Rd×N denote the sample matrix whose i-th column Xi ∈ Rd is the i-th input sample. Abusing notation, we sometimes also use X as a random variable (e.g., a sample drawn from P ).
  We use πP ∈ Rm to denote the parental configuration probabilities of P . That is, πPk = PrX∼P [X ∈ Πk]. For a set S of samples, we use πS ∈ Rm to denote the empirical parental configuration probabilities over S: πSk = PrX [X ∈ Πk] where X is uniformly drawn from S.
  Balance and Minimum Configuration Probability. We say a Bayesian network P is c-balanced if all conditional probabilities of P are between c and 1− c. We use α for the minimum probability of parental configuration of P : α = mink πPk .
  In this paper, we assume that the ground-truth Bayesian network is c-balanced, and its minimum parental configuration probability α satisfies that α = Ω ( ( √ ln(1/ ))2/3c−1/3) ) . Without loss of generality, we further assume that both c and α are given to the algorithm.
sections.heading: 2.1 TOTAL VARIATION DISTANCE BETWEEN BAYESIAN NETWORKS
sections.text: Let P and Q be two distributions supported on a finite domain D. For a set of outcomes A, let P (A) = PrX∼P [X ∈ A]. The total variation distance between P and Q is defined as
  dTV (P,Q) = max A⊆D
  |P (A)−Q(A)| .
  For two balanced Bayesian networks that share the same structure, it is well-known that the closeness in their conditional probabilities implies their closeness in total variation distance. Formally, we use the following lemma from Cheng et al. (2018), which upper bounds the total variation distance between two Bayesian networks in terms of their conditional probabilities.
  Lemma 2.1 (Cheng et al. (2018)). Let P and Q be two Bayesian networks that share the same structure. Let p and q denote the conditional probability tables of P and Q respectively. We have
  (dTV (P,Q)) 2 ≤ 2 ∑ k √ πPk π Q k
  (pk − qk)2
  (pk + qk)(2− pk − qk) .
sections.heading: 2.2 EXPANDING THE DISTRIBUTION TO MATCH CONDITIONAL PROBABILITY TABLE
sections.text: Lemma 2.1 states that to learn a known-structure Bayesian network P , it is sufficient to learn its conditional probabilities p. However, a given coordinate of X ∼ P may contain information about multiple conditional probabilities (depending on which parental configuration happens).
  To address this issue, we use a similar approach as in Cheng et al. (2018). We expand each sampleX into an m-dimensional vector f(X, q), such that each coordinate of f(X, q) corresponds to an entry in the conditional probability table. Intuitively, q ∈ Rm is our current guess for p, and initially we set q to be the empirical conditional probabilities. We use q to fill in the missing entries in f(X, q) for which the parental configurations fail to happen. Definition 2.2. Let f(X, q) for {0, 1}d × Rm → Rm be defined as follows:
  f(X, q)i,a = { Xi − qi,a X ∈ Πi,a 0 otherwise
  When X ∼ P and q = p, the distribution of f(X, p) has many good properties. Using the conditional independence of Bayesian networks, we can compute the first and second moment of f(X, p) and show that f(X, p) has subgaussian tails. Lemma 2.3. For X ∼ P and f(X, p) as defined in Definition 2.2, we have (i) E(f(X, p)) = 0. (ii) Cov[f(X, p)] = diag(πP ◦ p ◦ (1− p)).
  (iii) For any unit vector v ∈ Rm, we have PrX∼P [ |v>f(X, p)| ≥ T ] ≤ 2 exp(−T 2/2).
  We defer the proof of Lemma 2.3 to Appendix A. A slightly stronger version of Lemma 2.3 was proved in Cheng et al. (2018), which discusses tail bounds for f(X, q). For our analysis, Lemma 2.3 is sufficient.
  For general values of q, we can similarly compute the mean of f(X, q): Lemma 2.4. Let πP denote the parental configuration of P . For X ∼ P and f(X, q) as defined in Definition 2.2, we have E[f(X, q)] = πP ◦ (p− q).
sections.heading: 2.3 DETERMINISTIC CONDITIONS ON GOOD SAMPLES
sections.text: To avoid dealing with the randomness of the good samples, we require the following deterministic conditions to hold for the original set G? of N good samples (before the adversary’s corruption).
  We prove in Appendix A that these three conditions hold simultaneously with probability at least 1− τ if we draw N = Ω(m log(m/τ)/ 2) samples from P .
  The first condition states that we can obtain a good estimation of p from G?. Let pG ?
  denote the empirical conditional probabilities over G?. We have∥∥∥√πP ◦ (p− pG?)∥∥∥
  2 ≤ O( ) . (1)
  The second condition says that we can estimate the parental configuration probabilities πP from any (1− 2 )-fraction of G?. Formally, for any subset T ⊂ G? with |T | ≥ (1− 2 )N , we have∥∥πT − πP∥∥∞ ≤ O( ) . (2) The third condition is that the empirical mean and covariance of any (1 − 2 )-fraction of G? are very close to the true mean and covariance of f(X, p). Formally, for any subset T ⊂ G? with |T | ≥ (1− 2 )N , we require the following to hold for δ1 =
  √ ln 1/ and δ2 = ln(1/ ):∥∥∥∥∥ 1|T |∑
  i∈T f(Xi, p) ∥∥∥∥∥ 2 ≤ O(δ1) , ∥∥∥∥∥ 1|T |∑ i∈T f(Xi, p)f(Xi, p) > − Σ ∥∥∥∥∥ 2 ≤ O(δ2) , (3)
  where Σ = Cov[f(X, p)] = diag(πP ◦ p ◦ (1− p)).
sections.heading: 2.4 ROBUST MEAN ESTIMATION AND STABILITY CONDITIONS
sections.text: Robust mean estimation is the problem of learning the mean of a d-dimensional distribution from an -corrupted set of samples. As we will see in later sections, to robustly learn Bayesian networks, we repeatedly use robust mean estimation algorithms as a subroutine.
  Recent work (Diakonikolas et al., 2019a; Lai et al., 2016) gave the first polynomial-time algorithms for robust mean estimation with dimension-independent error guarantees. The key observation in Diakonikolas et al. (2019a) is the following: if the empirical mean is inaccurate, then many samples must be far from the true mean in roughly the same direction. Consequently, these samples must alter the variance in this direction more than they distort the mean. Therefore, if the empirical covariance behaves as we expect it to be, then the empirical mean provides a good estimate to the true mean.
  Many robust mean estimation algorithms follow the above intuition, and they require the following stability condition to work (Definition 2.5). Roughly speaking, the stability condition states that the mean and covariance of the good samples are close to that of the true distribution, and more importantly, this continues to hold if we remove any 2 -fraction of the samples.
  Definition 2.5 (Stability Condition (see, e.g., Diakonikolas & Kane (2019))). Fix 0 < < 12 . Fix a d-dimensional distribution X with mean µX . We say a set S of samples is ( , β, γ)-stable with respect to X , if for every subset T ⊂ S with |T | ≥ (1− 2 )|S|, the following conditions hold:
  (i) ∥∥∥ 1|T |∑X∈T (X − µX)∥∥∥
  2 ≤ β , (ii) ∥∥∥ 1|T |∑X∈T (X − µX) (X − µX)> − I∥∥∥ 2 ≤ γ .
  Subsequent work (Cheng et al., 2019a; Dong et al., 2019; Depersin & Lecué, 2019) improved the runtime of robust mean estimation to nearly-linear time. Formally, we use the following result from Dong et al. (2019). A set S is an -corrupted version of a set T if |S| = |T | and |S \T | ≤ |S|. Lemma 2.6 (Robust Mean Estimation in Nearly-Linear Time (Dong et al., 2019)). Fix a set of N samples G? in Rd. Suppose G? is ( , β, γ)-stable with respect to a d-dimensional distribution X with mean µX ∈ Rd. Let S be an -corrupted version of G?. Given as input S, , β, γ, there exists an algorithm that can output an estimator µ̂ ∈ Rd in time Õ(Nd), such that with high probability, ‖µ̂− µX‖2 ≤ O( √ γ + β + √ log 1/ ) .
  As we will see later, a black-box use of Lemma 2.6 does not give the desired runtime in our setting. Instead, we extend Lemma 2.6 to handle sparse input such that it runs in time nearly-linear in the number of non-zeros in the input (see Lemma 3.3).
sections.heading: 3 OVERVIEW OF OUR APPROACH
sections.text: In this section, we give an overview of our approach and highlight some of our key technical results.
  To robustly learn the ground-truth Bayesian network P , it is sufficient to learn its conditional probabilities p ∈ Rm. At a high level, we start with a guess q ∈ Rm for p and then iteratively improve our guess to get closer to p. For any q ∈ Rm, we can expand the input samples into m-dimensional vectors f(X, q) as in Definition 2.2. We first show that the expectation of f(X, q) gives us useful information about (p− q). Recall that πP is the parental configuration probabilities of P . By Lemma 2.4, we have
  EX∼P [f(X, q)] = πP ◦ (p− q) .
  Note that if we had access to this expectation and the vector πP , we could recover p immediately: we can set q′ = E[f(X, q)] ◦ (1/πP ) + q which simplifies to q′ = p. Note that since S is an -corrupted set of samples ofP , we know that {f(Xi, q)}i∈S is an -corrupted set of samples of the distribution f(X, q) (with X ∼ P ). Therefore, we can run robust mean estimation algorithms on {f(Xi, q)}i∈S to learn E[f(X, q)]. It turns out a good approximation of E[f(X, q)] can help us improve our current guess q.
  There are two main difficulties in getting this approach to work.
  The first difficulty is that, to use robust mean estimation algorithms, we need to show that f(X, q) satisfies the stability condition in Definition 2.5. This requires us to analyze the first two moments and tail bounds of f(X, q). Consider the second moment for example. Ideally, we would like to have a statement of the form Cov[f(X, q)] ≈ Cov[f(X, p)] + (p − q)(p − q)>, but this is false because we only have f(X, p)k − f(X, q)k = (p − q)k if the k-th parental configuration happens for X . Intuitively, the “error” (p − q) is shattered into all samples where each sample only gives d out of m coordinates of (p− q), and there is no succinct representation for Cov[f(X, q)]. The second difficulty is that f(X, q) ism-dimensional. We cannot explicitly write down all the samples {f(Xi, q)}Ni=1, because this takes time Ω(Nm), which could be much slower than our desired running time of Õ(Nd). Similarly, a black-box use of nearly-linear time robust mean estimation algorithms (e.g., Lemma 2.6) runs in time Ω(Nm), which is too slow.
  In the rest of this section, we explain how we handle these two issues.
  Stability Condition of f(X, q). Because the second-order stability condition in Lemma 2.3 is defined with respect to I , we first scale the samples so that the covariance of f(X, p) becomes I . Lemma 2.3 shows that Cov[f(X, p)] = diag(πP ◦p◦(1−p)). To make it close to I , we can multiply the k-th coordinate of f(X, p) by (πPk pk(1 − pk))−1/2. However, we do not know the exact value of πP or p, instead we use the corresponding empirical estimates πS and qS (see Algorithm 1). Definition 3.1. Let πS and qS denote the parental configuration probabilities and conditional means estimated over S. Let s = 1/ √ (πS ◦ qS ◦ (1− qS)). Throughout this paper, for a vector v ∈ Rm,
  we use v̂ ∈ Rm to denote v ◦ s. In particular, we have X̂i = Xi ◦ s (and similarly p̂, q̂, f̂(x, q)).
  Now we analyze the concentration bounds for f̂(X, q). Formally, we prove the following lemma. Lemma 3.2. Assume the conditions in Section 2.3 hold for the original set of good samples G?. Then, for δ1 = √ log 1/ and δ2 = log(1/ ), the set of samples { f̂(Xi, q) } i∈G? is(
  , O ( δ1√
  αc + ‖p̂− q̂‖2 ) , O ( δ2 αc +B + √ B )) -stable,
  where B = ‖ √ πP ◦ (p̂− q̂)‖22.
  We provide some intuition for Lemma 3.2 and defer its proof to Appendix B.
  For the first moment, the difference between E[f̂(X, q)] and the empirical mean of f̂(X, q) comes from several places. Even if q = p, we would incur an error of δ1 from the concentration bound in Equation equation 3, which is at most δ1(αc)−1/2 after the scaling by s. Moreover, on average πPk fraction of the samples gives us information about (p̂ − q̂)k. Since an -fraction of the samples are removed when proving stability, we may only have (πPk − )-fraction instead, which introduces an error of ‖p̂− q̂‖2. This is why the first-moment parameter is ( δ1(αc) −1/2 + ‖p̂− q̂‖2 ) .
  For the second moment, after the scaling, we have Cov[f̂(X, p)] ≈ I . Ideally, we would like to prove Cov[f̂(X, q)] ≈ I + (πP ◦ (p̂− q̂))(πP ◦ (p̂− q̂))>, but this is too good to be true. For two coordinates k 6= `, whether a sample gives information about (p̂− q̂)k or (p̂− q̂)` is not independent. We can upper bound the probability that both parental configurations happen by min(πPk , π P ` ). If they were independent we would have a bound of πPk π P ` . The difference in these two upper bounds
  is intuitively why √ πP appears in the second-moment parameter. See Appendix B for more details.
  Robust Mean Estimation with Sparse Input. To overcome the second difficulty, we exploit the sparsity of the expanded vectors. Observe that each vector f(X, q) is guaranteed to be d-sparse because exactly d parental configuration can happen (see Definition 2.2). The same is true for f̂(X, q) because scaling does not change the number of nonzeros. Therefore, there are in total O(Nd) nonzero entries in the set of samples {f̂(X, q)}i∈S . We develop a robust mean estimation algorithm that runs in time nearly-linear in the number of nonzeros in the input. Combined with the above argument, if we only invoke this mean estimation algorithm polylogarithmic times, we can get the desired running time of Õ(Nd).
  Lemma 3.3. Consider the same setting as in Lemma 2.6. Suppose ‖X‖2 ≤ R for all X ∈ S. There is an algorithm Amean with the same error guarantee that runs in time Õ(logR ·(nnz(S)+N+d)) where nnz(S) is the total number of nonzeros in S. That is, given an -corrupted version of an ( , β, γ)-stable set of N samples w.r.t. a d-dimensional distribution with mean µX , the algorithm Amean outputs an estimator µ̂ ∈ Rd in time Õ(logR · (nnz(S) + N + d)) such that with high probability, ‖µ̂− µX‖2 ≤ O( √ γ + β + √ log 1/ ) .
  We prove Lemma 3.3 by extending the algorithm in Dong et al. (2019) to handle sparse input. The main computation bottleneck of recent nearly-linear time robust mean estimation algorithms (Cheng et al., 2018; Dong et al., 2019) is in using the matrix multiplicative weight update (MMWU) method. In each iteration of MMWU, a score is computed for each sample. Roughly speaking, this score indicates whether one should continue to increase the weight on the corresponding sample. Previous algorithms use the Johnson-Lindenstrauss lemma to approximate the scores for all N samples simultaneously. We show that the sparsity of the input vectors allows for faster application of the Johnson-Lindenstrauss lemma, and all N scores can be computed in time nearly-linear in nnz(S).
  We defer the proof of Lemma 3.3 to Appendix C.
sections.heading: 4 ROBUST LEARNING OF BAYESIAN NETWORKS IN NEARLY-LINEAR TIME
sections.text: In this section, we prove our main result. We present our algorithm (Algorithm 1) and prove its correctness and analyze its running time (Theorem 4.1). Theorem 4.1. Fix 0 < < 0 where 0 is a sufficiently small universal constant. Let P be a c-balanced Bayesian network on {0, 1}d with known structure H . Let α be the minimum parental configuration probability of P . Assume α = Ω̃( 2/3c−1/3).
  Let S be an -corrupted set of N = Ω̃(m/ 2) samples drawn from P . Given H , S, , c, and α, Algorithm 1 outputs a Bayesian network Q in time Õ(Nd) such that, with probability at least 9/10, dTV (P,Q) ≤ O( √ log(1/ )/ √ αc).
  The c and α terms in the error guarantee also appear in prior work (Cheng et al., 2018). Removing this dependence is an important technical question that is beyond the scope of this paper.
  Theorem 4.1 follows from three key technical lemmas. At the beginning of Algorithm 1, we first scale all the input vectors as in Definition 3.1. We maintain a guess q for p and gradually move it closer to p. In our analysis, we track our progress by the `2-norm of πP ◦ (p̂− q̂). Initially, we set q to be the empirical conditional mean over S. Lemma 4.2 proves that∥∥πP ◦ (p̂− q̂′)∥∥
  2 is not too large for our first guess. Lemma 4.3 shows that, as long as q is still relatively far from p, we can compute a new guess such that ∥∥πP ◦ (p̂− q̂)∥∥
  2 decreases by a constant factor. Lemma 4.4 states that, when the algorithm terminates and ∥∥πP ◦ (p̂− q̂)∥∥
  2 is small,
  we can conclude that the output Q is close to the ground-truth P in total variation distance.
  In the following three lemmas, we consider the same setting as in Theorem 4.1 and assume the conditions in Section 2.3 hold. Lemma 4.2 (Initialization). In Algorithm 1, we have∥∥πP ◦ (p̂− q̂0)∥∥
  2 ≤ O(
  √ d/ √ αc) .
  Lemma 4.3 (Iterative Refinement). Fix an iteration t in Algorithm 1. Assume the robust mean estimation algorithm Amean succeeds. If ∥∥πP ◦ (p̂− q̂t)∥∥ 2 ≤ ρt and ρt = Ω( √ log(1/ )/ √ αc),
  then we have ∥∥πP ◦ (p̂− q̂t+1)∥∥
  2 ≤ c1ρt for some universal constant c1 < 1.
  Lemma 4.4. Let Q be a Bayesian network that has the same structure as P . Suppose that (1) P is c-balanced, (2) α = Ω(r + /c), and (3) ∥∥πP ◦ (p̂− q̂)∥∥ 2 ≤ r/2. Then we have dTV (P,Q) ≤ r.
  We defer the proofs of Lemmas 4.2, 4.3, and 4.4 to Appendix D and we first prove Theorem 4.1.
  Proof of Theorem 4.1. We first prove the correctness of Algorithm 1.
  Algorithm 1: Robustly Learning Bayesian Networks Input : The dependency graph H of a c-balanced Bayesian network P with minimum
  parental configuration α, an -corrupted set S of N = Ω̃(m/ 2) samples {Xi}Ni=1 drawn from P , and the values of , c and α.
  Output: A Bayesian network Q such that, with probability at least 9/10, dTV (P,Q) ≤ O( √ log(1/ )/ √ αc). Compute the empirical probabilities πS where πS(i, a) = PrX∈S [Πi,a]; Compute the empirical conditional probabilities qS where qS(i, a) = PrX∈S [X(i) = 1 | Πi,a]; Compute the scaling vector s = 1/ √ (πS ◦ qS ◦ (1− qS));
  Let T = O(log d) and q0 = qS ; Let ρ0 = O( √ d/ √ αc). (We maintain upper bounds ρt s.t. ∥∥πP ◦ (p̂− q̂t)∥∥ 2 ≤ ρt for all t); for t = 0 to T − 1 do βt = O( ρt/α), γt = O((ρt)2/α+ ρt/ √ α) ;
  Solve a robust mean estimation problem. Let ν = Amean ( {f̂(Xi, qt)i∈S}, , βt, γt);
  qt+1 = ν ◦ (1/s) ◦ (1/πS) + qt; ρt+1 = c1ρt; return the Bayesian network Q with graph H and conditional probabilities qT ;
  The original set of N = Ω(m log(m/ )/ 2) good samples drawn from P satisfies the conditions in Section 2.3 with probability at least 1− 120 . With high probability, the robust mean estimation oracle Amean succeeds in all iterations. For the rest of this proof, we assume the above conditions hold, which by a union bound happens with probability at least 9/10.
  From Lemma 4.2, we have the following condition on the initial estimate q0.∥∥πP ◦ (p̂− q̂0)∥∥ 2 = O( √ d/ √ αc) . We start with an upper bound ρ0 of ∥∥πP ◦ (p̂− q̂0)∥∥
  2 = O(
  √ d/ √ αc). By Lemma 4.3, in each iteration, if ρt = Ω( √ log(1/ )/ √ αc), we can obtain a new estimate qt+1 and an upper bound ρt+1
  on ∥∥πP ◦ (p̂− q̂t)∥∥
  2 such that ρt+1 is smaller than ρt by a constant factor. Hence after O(log(d))
  iterations, we can get a vector qt such that∥∥πP ◦ (p̂− q̂t)∥∥ 2 = O( √ log(1/ )/ √ αc) .
  Let Q be the Bayesian network with conditional probability table qt. The assumption that α = Ω̃( 2/3c−1/3) allows us to apply Lemma 4.4 with r = O( √ log(1/ )/ √ αc), which gives the claimed upper bound on dTV (P,Q).
  Now we analyze the runtime of Algorithm 1. First, qS and πS can be computed in time Õ(Nd) because each sample only affects d entries of q. We do not explicitly write down f(X, q). In each iteration, we solve a robust mean estimation problem with input { f̂(Xi, q t) } i∈S , which takes time Õ(Nd). This is because there are N input vectors, each vector is d-sparse, and the robust mean estimation algorithm runs in time nearly-linear in the number of nonzeros in the input (Lemma 3.3). We can compute qt+1 = ν ◦ (1/s) ◦ (1/πS) + qt in time in time O(m). Since there are O(log d) iterations, the overall running time is
  Õ(Nd) +O(log d) ( Õ(Nd) +O(m) ) = Õ(Nd) .
sections.heading: ACKNOWLEDGMENTS
sections.text: Part of this work was done while Yu Cheng was visiting the Institute of Advanced Study. Part of this work was done while Honghao Lin was an undergraduate student at Shanghai Jiao Tong University.
sections.heading: A DETERMINISTIC CONDITIONS ON GOOD SAMPLES
sections.text: In this section, we will first prove Lemma 2.3, then we prove that the deterministic conditions in Section 2.3 hold with high probability if we take enough samples. Lemma A.1. For X ∼ P and f(X, p) as defined in Definition 2.2, we have (i) E(f(X, p)) = 0.
  (ii) Cov[f(X, p)] = diag(πP ◦ p ◦ (1− p)). (iii) For any unit vector v ∈ Rm, we have PrX∼P [ |v>f(X, p)| ≥ T ] ≤ 2 exp(−T 2/2).
  Proof. We first claim that EX∼P [f(X, p)k|f(X, p)1, ..., f(X, p)k−1] = 0 for all k ∈ [m]. Let k = (i, a), conditioned on f(X, p)1, ..., f(X, p)k−1, the event πi,a may or may not happen. A simple calculation shows that in both cases, we have EX∼P [f(X, p)k|f(X, p)1, ..., f(X, p)k−1] = 0. For (i), we have E[f(X, p)] = πP ◦ p+ (1− πP ) ◦ p− p = 0. For (ii), we first show that for any (i, a) 6= (j, b), we have E[f(X, p)i,af(X, p)j,b] = 0. For the case i = j, we can see at least one of Πi,a and Πj,b does not occur, so f(X, p)i,af(X, p)j,b is always 0. For the case i 6= j, we assume without loss of generality that i > j, then we have E[f(X, p)i,a|f(X, p)j,b] = 0.
  For all (i, a) ∈ [m], we have E[f(X, p)2i,a] = πPi,aE[(X − pi,a)2|Πi,a] = πPi,api,a(1 − pi,a). Combining these two, we get Cov[f(X, p)] = diag(πP ◦ p ◦ (1− p)). For (iii), we recall that EX∼P [f(X, p)k|f(X, p)1, ..., f(X, p)k−1] = 0, thus the sequence∑` k=1 vkf(X, q)k for 1 ≤ ` ≤ m is a martingale, and we can apply Azuma’s inequality. Note
  that |vk| ≥ |vkf(X, p)k|, hence we have PrX∼P [ |v>f(X, p)| ≥ T ] ≤ 2 exp(−T 2/2 ‖v‖22) = 2 exp(−T 2/2).
  The conditions in Equations equation 1 and equation 2 are proved in Lemma A.2, and the conditions in Equation equation 3 are proved in Corollary A.4. Lemma A.2. Let P be a Bayesian network. LetG? be a set of Ω((m log(m/τ))/ 2) samples drawn from P . Let πG ? and pG ?
  be the empirical parental configuration probabilities and conditional probabilities of P given by G?. Then with probability 1− τ , the following conditions hold: (i) For any subset T ⊂ G? with |T | ≥ (1− 2 )N , we have
  ∥∥πT − πP∥∥∞ ≤ O( ) . (ii) ∥∥∥√πP ◦ (p− pG?)∥∥∥
  2 ≤ O( ) .
  Proof. For (i), first consider the case of T = G? and fix an entry 1 ≤ k ≤ m in the conditional probability table. Because each sample is drawn independently from P , by the Chernoff bound, we have that when N = Ω(log(m/τ)/ 2), |πPk − πG ?
  k | ≤ holds with probability at least 1 − τ/m. Hence, after taking an union bound over k, we have that ∥∥πT − πP∥∥∞ ≤ holds with probability at least 1− τ . Now for a general subset T ⊂ G?, notice that removing O( )-fraction of samples can change πT by at most O( ). Thus, condition (i) holds with probability at least 1− τ .
  For (ii), for any k = (i, a), note that pG ? k is estimated from π G?
  k N samples. In these samples, the parental configuration Πk happens and the value of Xi is decided independently. By the Chernoff bound and the union bound, we get that when N = Ω((m log(m/τ))/ 2), |pG?k −pk| ≤ / √ mπG ? k holds for every k with probability at least 1− τ , which implies∥∥∥√πG? ◦ (p− pG?)∥∥∥ 2 ≤ O( ).
  Combining this with ∥∥πP − πG?∥∥∞ ≤ , we get that condition (ii) holds.
  To prove Equation equation 3, we use the following concentration bounds for subgaussian distributions. Recall that a distribution D on Rd with mean µ is subgaussian if for any unit vector v ∈ Rd we have Prx∼D[|〈v, x− µ〉| ≥ t] ≤ exp(−ct2), where c is a universal constant. Lemma A.3. Let G? be a set of N = Ω(( √ log 1/ )−2(d + log(1/τ))) samples drawn from a d-dimensional subgaussian distribution with mean µ and covariance matrix Σ I . Here A B means that B − A is a positive semi-definite matrix. Then, with probability 1 − τ , the following conditions hold:
  For δ1 = c1( √
  log 1/ ) and δ2 = c1( log 1/ ) where c1 is an universal constant, we have that for any subset T ⊂ G? with |T | ≥ (1− 2 )N ,∥∥∥∥∥ 1|T |∑
  i∈T (Xi − µ) ∥∥∥∥∥ 2 ≤ δ1 , ∥∥∥∥∥ 1|T |∑ i∈T (Xi − µ)(Xi − µ)> − Σ ∥∥∥∥∥ 2 ≤ δ2 (4)
  A special case of Lemma A.3 where Σ = I is proved in Diakonikolas et al. (2019a). The proof for the general case where Σ I is almost identical. In particular, the concentration inequalities used in Diakonikolas et al. (2019a) for subgaussian distributions still hold when Σ I (see, e.g., Vershynin (2010)).
  From Lemma A.3 and 2.3, we have the following corollary: Corollary A.4. Let G? be a set of N = Ω(( √
  log 1/ )−2(m + log(1/τ))) samples drawn P . Then, with probability 1 − τ , the following conditions to hold: For δ1 = c1( √ log 1/ ) and δ2 = c1( log 1/ ), where c1 is an universal constant, we have that for any subset T ⊂ G? with |T | ≥ (1− 2 )N ,∥∥∥∥∥ 1|T |∑
  i∈T f(Xi, p) ∥∥∥∥∥ 2 ≤ O(δ1) , ∥∥∥∥∥ 1|T |∑ i∈T f(Xi, p)f(Xi, p) > − Σ ∥∥∥∥∥ 2 ≤ O(δ2) , (5)
  where Σ = Cov[f(X, p)] = diag(πP ◦ p ◦ (1− p)).
  B STABILITY CONDITION OF f̂(X, q)
  In this section, we prove the stability condition for the samples f̂(X, q) (Lemma 3.2). Recall the definition of f̂(X, q) from Definitions 2.2 and 3.1. We first restate Lemma 3.2.
  Lemma 3.2. Assume the conditions in Section 2.3 hold for the original set of good samples G?. Then, for δ1 = √ log 1/ and δ2 = log(1/ ), the set of samples { f̂(Xi, q) } i∈G? is(
  , O ( δ1√
  αc + ‖p̂− q̂‖2 ) , O ( δ2 αc +B + √ B )) -stable,
  where B = ‖ √ πP ◦ (p̂− q̂)‖22.
  We will prove the stability of f(X, q). The stability of f̂(X, q) follows directly. We introduce a matrix CD,q which is crucial in proving the stability of f(X, q). Intuitively, the matrix CD,q is related to the difference in the covariance of f(X, p) and that of f(X, q) on the sample set D. Definition B.1. For any set D of samples {Xi}i∈D, we define the following m×m matrix
  CD,q = 1 |D| ∑ i∈D (f(Xi, p)− f(Xi, q))(f(Xi, p)− f(Xi, q))> .
  Observe that for x ∈ {0, 1}d with x /∈ Πk, we have f(x, p)k = f(x, q)k = 0. On the other hand, if x ∈ Πk for some k = (i, a), we have f(x, p)k − f(x, q)k = (xi − pk)− (xi − qk) = qk − pk. In the very special case where all parental configurations happen (i.e., a binary product distribution), we would have CD,q = (p− q)(p− q)>. However, in general the information related to (p− q) is spread among the samples. We show that even though CD,q does not have a succinct representation, we can prove the following upper bound on the spectral norm of CD,q.
  Lemma B.2. ‖CD,q‖2 ≤ ∑ k π D k (pk − qk)2 .
  Proof. For notational convenience, let C = CD,q. For every 1 ≤ k, ` ≤ m, we have
  |Ck,`| = ∣∣∣(Pr D [Πk ∧Π`])(pk − qk)(p` − q`) ∣∣∣ ≤ min{πDk , πD` } · |(pk − qk)(p` − q`)| ≤ (√ πDk |(pk − qk)| ) · (√ πD` |(p` − q`)| )
  We can upper bound the spectral norm of C in term of its Frobenius norm:
  ‖C‖22 ≤ ‖C‖ 2 F = ∑ k,` C2k,` ≤ ∑ k,` ( πDk (pk − qk)2 ) ( πD` (p` − q`)2 ) ≤ (∑ k πDk (pk − qk)2 )2 .
  The following lemma essentially proves the stability of f(X, q), except that in the second-order condition, we should have Cov(f(X, q)) instead of Σ. We will bridge this gap in Lemma B.4. Lemma B.3. Assume the conditions in Section 2.3 hold. For δ1 = √ log 1/ and δ2 = log(1/ ), we have that for any subset T ⊂ G? with |T | ≥ (1− )|G?|,∥∥∥∥∥ 1|T |∑ i∈T (f(Xi, q)− πP ◦ (p− q)) ∥∥∥∥∥ 2
  ≤ O(δ1 + ‖p− q‖2) , and∥∥∥∥∥ 1|T |∑ i∈T (f(Xi, q)− πP ◦ (p− q))(f(Xi, q)− πP ◦ (p− q))> − Σ ∥∥∥∥∥ 2 ≤ O(δ2 +B + √ B)
  where B = ∥∥∥√πP ◦ (p− q)∥∥∥2
  2 ≤ 1α
  ∥∥πP ◦ (p− q)∥∥2 2 , and Σ = diag(πP ◦ p ◦ (1 − p)) is the true
  covariance of f(X, p) .
  Proof. For the first moment, we have∥∥∥∥∥ 1|T |∑ i∈T (f(Xi, q)− πP ◦ (p− q)) ∥∥∥∥∥ 2
  = ∥∥∥∥∥ 1|T |∑ i∈T (f(Xi, p) + f(Xi, q)− f(Xi, p)− πP ◦ (p− q)) ∥∥∥∥∥ 2
  ≤ ∥∥∥∥∥ 1|T |∑ i∈T f(Xi, p) ∥∥∥∥∥ 2 + ∥∥∥∥∥ 1|T |∑ i∈T (f(Xi, q)− f(Xi, p)− πP (p− q)) ∥∥∥∥∥ 2
  = ∥∥∥∥∥ 1|T |∑ i∈T f(Xi, p) ∥∥∥∥∥ 2 + ∥∥(πT − πP ) ◦ (p− q)∥∥ 2 = O(δ1 + ‖p− q‖2) .
  For the second moment, consider any unit vector v ∈ Rm. We have
  v>
  ( 1
  |T | ∑ i∈T f(Xi, q)f(Xi, q) >
  ) v = 1
  |T | ∑ i∈T 〈f(Xi, q), v〉2
  = 1 |T | ∑ i∈T ( 〈f(Xi, p), v〉2 + 〈f(Xi, p)− f(Xi, q), v〉2 + 2〈f(Xi, p), v〉〈f(Xi, p)− f(Xi, q), v〉 ) ≤ 1 |T | ∑ i∈T ( 〈f(Xi, p), v〉2 + 〈f(Xi, p)− f(Xi, q), v〉2 )
  + 2 √√√√ 1 |T | (∑ i∈T 〈f(Xi, p), v〉2 )( 1 |T | ∑ i∈T 〈f(Xi, p)− f(Xi, q), v〉2 )
  where the last inequality follows from the Cauchy-Schwarz inequality. Therefore, we have∥∥∥∥∥ 1|T |∑ i∈T f(Xi, q)f(Xi, q) > − Σ ∥∥∥∥∥ 2
  ≤ ∥∥∥∥∥ 1|T |∑ i∈T f(Xi, p)f(Xi, p) > − Σ ∥∥∥∥∥ 2
  + ∥∥∥∥∥ 1|T |∑ i∈T (f(Xi, p)− f(Xi, q))(f(Xi, p)− f(Xi, q))> ∥∥∥∥∥ 2
  + 2 √√√√∥∥∥∥∥ 1|T |∑ i∈T f(Xi, p)f(Xi, p)> ∥∥∥∥∥ 2 ∥∥∥∥∥ 1|T |∑ i∈T (f(Xi, p)− f(Xi, q))(f(Xi, p)− f(Xi, q))> ∥∥∥∥∥ 2
  ≤ δ2 + ‖CT,q‖2 + 2 √ 1 + δ2 √ ‖CT,q‖2 = O ( δ2 + ‖CT,q‖2 + √ ‖CT,q‖2 ) Finally, we show that the second moment matrix 1|T | ∑ i∈T f(Xi, q)f(Xi, q)
  > is not too far from the empirical covariance matrix of f(X, q).∥∥∥∥∥ 1|T |∑
  i∈T f(Xi, q)f(Xi, q) > − 1 |T | ∑ i∈T (f(Xi, q)− πP ◦ (p− q))(f(Xi, q)− πP ◦ (p− q))> ∥∥∥∥∥ 2
  ≤ 2 ∥∥∥∥∥ 1|T |∑ i∈T (f(Xi, q)− πP ◦ (p− q)) ∥∥∥∥∥ 2 ‖πP ◦ (p− q)‖2 + ‖πP ◦ (p− q)‖ 2 2
  ≤ O ((δ1 + ‖p− q‖2) ‖πP ◦ (p− q)‖2) ≤ O ( ‖πP ◦ (p− q)‖22 ) .
  Putting everything together and using Lemma B.2, we conclude this proof.
  The stability of f̂(X, q) follows from the stability of f(X, q) (Lemma B.3), scaling all samples by s, and replacing Σ̂ with I in the second-order condition using Lemma B.4.
  Lemma B.4. Assume the conditions in Section 2.3 hold. Then after scaling, we have∥∥∥Σ̂− I∥∥∥ 2 ≤ O ( αc ) .
  where Σ̂ is the covariance matrix of f̂(X, p).
  Proof. We recall that πPk pk(1 − pk) ≥ 12π P k min(pk, 1 − pk) = Ω(αc). Because ‖s‖∞ =
  O(1/ √ αc), it suffices to show that
  ‖Cov(s ◦ f(X, p))− I‖2 ≤ O( ) .
  In other words, we need to show∥∥πP ◦ p ◦ (1− p)− πS ◦ qS ◦ (1− qS)∥∥∞ = O( ) . Let πG ? and pG ?
  be the empirical parental configuration probabilities and conditional probabilities of P given by G?. We first prove that∥∥∥πG? ◦ pG? ◦ (1− pG?)− πS ◦ qS ◦ (1− qS)∥∥∥
  ∞ = O( ) .
  Note that ∥∥πG? − πS∥∥∞ = O( ) and qSk (1− qSk ) < 1, so it is sufficient to show that∥∥∥πG? ◦ pG? ◦ (1− pG?)− πG? ◦ qS ◦ (1− qS)∥∥∥
  ∞ = O( ) .
  We have∥∥∥πG? ◦ pG? ◦ (1− pG?)− πG? ◦ qS ◦ (1− qS)∥∥∥ ∞
  = ∥∥∥πG? ◦ pG? − πG? ◦ qS − πG? ◦ (pG? + qS) ◦ (pG? − qS)∥∥∥
  ∞ ≤ 3 ∥∥∥πG? ◦ (pG? − qS)∥∥∥ ∞ .
  Let nk denote the number of times that the event Πk happens over G?, and let tk be the number of times that X(i) = 1 when Πk = Πi,a happens. Because S is obtained by changing at most N samples in G?, we can get
  |πG ? k ◦ (pG ? k − qSk )| ≤ nk N · ( tk + N nk − N − tk nk ) = nk(nk + tk) N Nnk(nk − N) ≤ 2nk N 0.5Nnk = 4 .
  The last inequality follows from tk ≤ nk and nk − N ≥ 0.5nk (because we assume the minimum parental configuration probability is Ω( )).
  This concludes ∥∥πG? ◦ pG? ◦ (1− pG?)− πS ◦ qS ◦ (1− qS)∥∥∞ = O( ). Similarly, in order to prove that∥∥∥πP ◦ p ◦ (1− p)− πG? ◦ pG? ◦ (1− pG?)∥∥∥ ∞ = O( ) . We just need to show ∥∥πP ◦ (p− pG?)∥∥
  2 = O( ), which is follows from (ii) in Lemma A.2.
  An application of triangle inequality finishes this proof.
  We are now ready to prove Lemma 3.2.
  Proof of Lemma 3.2. By Lemma B.3 and the fact that ‖s‖∞ = O(1/ √ αc), we know that for any subset T ⊂ G? with |T | ≥ (1− )|G?|, we have∥∥∥∥∥ 1|T |∑ i∈T (f̂(Xi, q)− πP ◦ (p̂− q̂)) ∥∥∥∥∥ 2 ≤ O( δ1√ αc
  + ‖p̂− q̂‖2) , and∥∥∥∥∥ 1|T |∑ i∈T (f̂(Xi, q)− πP ◦ (p̂− q̂))(f̂(Xi, q)− πP ◦ (p̂− q̂))> − Σ̂ ∥∥∥∥∥ 2 ≤ O( δ2 αc +B + √ B)
  where B = ∥∥∥√πP ◦ (p̂− q̂)∥∥∥2
  2 , and Σ̂ is the true covariance of f̂(X, p). This is because the scaling
  is applied to all vectors on both sides of the inequalities, so we only need to scale the scalars δ1 and δ2 appropriately.
  We conclude the proof by replacing Σ̂ in the second-order condition with I using Lemma B.4.
sections.heading: C ROBUST MEAN ESTIMATION WITH SPARSE INPUT
sections.text: In this section, we give a robust mean estimation algorithm that runs in nearly input-sparsity time. We build on the following lemma, which is essentially the main result of Dong et al. (2019).
  Lemma C.1 (essentially Dong et al. (2019)). Given an -corrupted version of an ( , β, γ)-stable set S of N samples w.r.t. a d-dimensional distribution with mean µX . Suppose further that ‖X‖2 ≤ R for all X ∈ S, there is an algorithm outputs an estimator µ̂ ∈ Rd such that with high probability,
  ‖µ̂− µX‖2 ≤ O( √ γ + β + √ log 1/ ) .
  Moreover, this algorithm runs in time Õ((nnz(S)+N +d+T (Oapx)) · logR), where nnz(S) is the total number of nonzeros in the samples in S and T (Oapx)) is the runtime of an approximate score oracle defined in Definition C.2.
  Lemma C.1 is different from the way it is stated in Dong et al. (2019). This is because we use a more concise stability condition than the one in Dong et al. (2019). We will show that Lemma C.1 is equivalent to the version in Dong et al. (2019) in Appendix C.1.
  The computational bottleneck of the algorithm in Dong et al. (2019) is logarithmic uses of matrix multiplicative weight update (MMWU). In each iteration of every MMWU, they need to compute a score for each sample. Intuitively, these scores help the algorithm decide whether it should continue to increase the weight on each sample or not.
  We define some notations before we formally define the approximate score oracle. Let ∆N = {w ∈ RN : 0 ≤ wi ≤ 1, ∑ wi = 1} be theN -dimensional simplex. Given a set ofN samplesX1, ...,XN
  and a weight vector w ∈ ∆N , let µ(w) = 1|w| ∑ wiXi and Σ(w) = 1|w| ∑ wi(Xi − µ(w))(Xi − µ(w))> denote the empirical mean and covariance weighted by w. Definition C.2 (Approximate Score Oracle). Given as input a set of N samples X1, . . . , XN ∈ Rd, a sequence of t + 1 = O(log(d)) weight vectors w0, . . . , wt ∈ ∆N , and a parameter α > 0, an approximate score oracle Oapx outputs (1± 0.1)-approximations (τ̃i)Ni=1 to each of the N scores
  τi = (Xi − µ(wt))>U(Xi − µ(wt))
  for
  U = exp(α
  ∑t−1 i=0 Σ(w i))
  tr exp(α ∑t−1 i=0 Σ(w i)) .
  In addition, Oapx outputs a scalar q̃ such that |q̃ − q| ≤ 0.1q + 0.05 ∥∥Σ(wt)− I∥∥
  2 ,where q = 〈Σ(wt)− I, U〉 .
  These scores are computed using the Johnson-Lindenstrauss lemma. Our algorithm for computing these scores are given in Algorithm 2.
  Let r = O(logN log(1/δ)), ` = O(log d), and Q ∈ Rr×d be a matrix with i.i.d entries drawn from N (0, 1/r). Algorithm 2 computes an r × d matrix
  A = Q · P`
  ( α
  2 t−1∑ i=0 Σ(wi)
  ) . (6)
  where P`(Y ) = ∑` j=0 1 j!Y j is a degree-` Taylor approximation to exp(Y ).
  The estimates for the individual scores are then given by
  τ̃i = 1
  tr(AA>) ∥∥A(Xi − µ(wt))∥∥2 (7) and the estimate for q is given by
  q̃ = N∑ i=1 (τ̃i − 1) . (8)
  Algorithm 2: Nearly-linear time approximate score computation Input: A set S of N samples X1, . . . , XN ∈ Rd, a sequence of weight vectors w0, ..., wt, a parameter α, and a failure probability δ > 0. Let r = O(logN log(1/δ)) and ` = O(log d); Let Q ∈ Rr×d have entries drawn i.i.d. from N (0, 1/r); Compute the matrix A ∈ Rr×d as in Equation equation 6; return (τ̃i)Ni=1 given by Equation equation 7 and q̃ given by Equation equation 8;
  The correctness of Algorithm 2 was proved in Dong et al. (2019). Lemma C.3 (Dong et al. (2019)). With probability at least 1− δ, the output of Algorithm 2 satisfies |q̃ − q| ≤ 0.1q + 0.05 ‖Σ(wt)− I‖2 and |τ̃i − τi| ≤ 0.1τi for all 1 ≤ i ≤ N .
  Consequently, we only need to analyze the runtime of Algorithm 2.
  Lemma C.4. Algorithm 2 runs in time Õ(t · (N + d+ nnz(S)) · log(1/δ)).
  Proof. We first show that matrix A ∈ Rr×d can be computed in time
  Õ(t · (N + d+ nnz(S)) · log 1/δ) . We will multiply each row of Q (from the left) through the matrix polynomial to obtain A. Let v> ∈ R1×d be one of the rows of Q and let w ∈ RN be any weight vector. Observe that we can compute all ( v> (Xi − µ(w)) )N i=1 in time
  O ( N∑ i=1 nnz(Xi) +N + d ) = O(nnz(S) +N + d) .
  This is because we can compute µ(w) and v>µ(w) just once, and then compute v>Xi for every i and subtract v>µ(w) from it.
  Then, we can compute
  v>Σ(w) = v> ( N∑ i=1 wi(Xi − µ(w))(Xi − µ(w))> )
  = N∑ i=1 wi ( v>(Xi − µ(w)) ) X>i − ( N∑ i=1 wi ( v>(Xi − µ(w) )) µ(w)>
  as the sum of N sparse vectors subtracting a dense vector in time O(nnz(S) +N + d). Therefore, for any v ∈ Rm, we can evaluate v> ∑t−1 i=0 Σ(w i) in time O(t · (nnz(S) + d+N)).
  Because P` is a degree-` matrix polynomial of ∑t−1 i=0 Σ(w i), we can use Horner’s method for poly-
  nomial evaluation to compute v>P` ( −α2 ∑t−1 i=0 Σ(w i) )
  in timeO(`·t·(nnz(S)+d+N)). We need to multiply each of r rows of A through, we can compute A in time O(r · ` · t · (nnz(S) + d+N)). It remains to show that (τ̃i)Ni=1 and q̃ as defined in Equations 7 and 8 can be computed quickly. Note that tr(AA>) is the entrywise inner product of A with itself, so it can be computed in time O(rd). The vectors (A(Xi − µ(wt))) N i=1 can be computed in time O (r · ( ∑ i nnz(Xi) + d)) = O(r · (nnz(S) + d)), because each AXi can be compute in time O(r ·nnz(Xi)) and Aµ(wt) can be computed only once in time O(rd). Because r = O(logN log(1/δ)), we can compute all τ̃i in time O(r · (nnz(S) + d)). Given the τ̃i’s, q̃ can be computed in O(N) time. Recall that r = O(logN log(1/δ)) and ` = O(log d). Putting everything together, the overall runtime of the oracle is
  O(r ·` · t ·(nnz(S)+d+N))+O(r ·(nnz(S)+d)+N) = Õ(t ·(nnz(S)+d+N) · log(1/δ)) .
  By Lemma C.4 and the fact that t = O(log d) and δ = 1/ poly(d), we can implement an approximate score oracle that succeeds with high probability and runs in time Õ(nnz(S) +N + d).
  Lemma 3.3 follows from Lemma C.1 and the correctness and the runtime of the approximate score oracle (Lemmas C.3, and C.4). (Note that we have nnz(S) = Nd, N = N , and d = m when invoking these lemmas.)
  C.1 PROOF OF LEMMA C.1
  In this section, we will show the equivalence between Lemma C.1 and Lemma C.6. Lemma C.6 is a restatement of the result of Dong et al. (2019) using their stability notations.
  We first state the stability condition used throughout Dong et al. (2019). Definition C.5 (Dong et al. (2019)). We say a set of points S is ( ,γ1,γ2,β1,β2)-good with respect to a distribution D with true mean µ if the following two properties hold:
  • ‖µ(S)− µ‖2 ≤ γ1, ∥∥∥ 1|S|∑i∈S (Xi − µ(S))(Xi − µ(S))> − I∥∥∥
  2 ≤ γ2 .
  • For any subset T ⊂ S so that |T | = 2 |S|, we have∥∥∥∥∥ 1|T |∑ i∈T Xi − µ ∥∥∥∥∥ 2 ≤ β1, ∥∥∥∥∥ 1|T |∑ i∈T (Xi − µ(S))(Xi − µ(S))> − I ∥∥∥∥∥ 2 ≤ β2 .
  Then, Dong et al. (2019) showed the following result.
  Lemma C.6. Let D be a distribution on Rd with unknown mean µ. Let 0 < < 0 for some universal constant 0. Let S be a set of N samples with S = Sg ∪ Sb\Sr where |Sb|, |Sr| ≤ |S|, and Sg is ( , γ1, γ2, β1, β2)-good with respect to D. Let Oapx be an approximate score oracle for S. Suppose ‖X‖2 ≤ R for all X ∈ S. Then, there is an algorithm QUEScoreFilter(S, Oapx, δ) that outputs a µ̂ such that with high probability,
  ‖µ̂− µ‖2 ≤ O( √ log 1/ + √ ξ + γ1) ,
  where ξ = γ2 + 2γ 2 1 + 4 2β21 + 2 β2 +O( log 1/ ) .
  Moreover, QUEScoreFilter makes O(logR log d) calls to the score oracle Oapx, and the rest of the algorithm runs in time Õ(N log(R)).
  We first show the connection between our stability notion (Definition 2.5) and theirs (Definition C.5).
  Lemma C.7. Fix a d-dimensional distribution D with mean µ, if a set S of N samples is ( , β, γ)stable with respect to D, then S is ( , β, γ, β/ , γ/ + 3β2/ 2)-good with respect to D.
  Proof. For any subset T ⊂ S with |T | = 2 |S|, we have∥∥∥∥∥ 1|T |∑ i∈T Xi − µ ∥∥∥∥∥ 2 = ∥∥∥∥∥∥ ( 1 |T | ∑ i∈S Xi − 1 2 µ ) −  1 |T | ∑ i∈S\T Xi − ( 1 2 − 1 ) µ ∥∥∥∥∥∥ 2
  = ∥∥∥∥∥∥ 12 ( 1 N ∑ i∈S Xi − µ ) − (1− 2 ) 2  1 (1− 2 )N ∑ i∈S\T Xi − µ ∥∥∥∥∥∥ 2
  ≤ 1 2 ∥∥∥∥∥ 1N ∑ i∈S Xi − µ ∥∥∥∥∥ 2 + 1− 2 2 ∥∥∥∥∥∥ 1(1− 2 )N ∑ i∈S\T Xi − µ ∥∥∥∥∥∥ 2
  ≤ 1 2 β + 1 2 β = β .
  The last line follows from the assumption that S is ( , β, γ)-stable with respect to D.
  Similarly, we have∥∥∥∥∥ 1|T |∑ i∈T (Xi − µ)(Xi − µ)> − I ∥∥∥∥∥ 2
  = ∥∥∥∥∥∥ 12 ( 1 N ∑ i∈S (Xi − µ)(Xi − µ)> − I ) − (1−2 )2 ( 1 (1−2 )N ∑ i∈S\T (Xi − µ)(Xi − µ)> − I )∥∥∥∥∥∥
  2
  ≤ 12 ∥∥∥∥∥ 1N ∑ i∈S (Xi − µ)(Xi − µ)> − I ∥∥∥∥∥ 2 + 1−2 2 ∥∥∥∥∥∥ 1(1−2 )N ∑ i∈S\T (Xi − µ)(Xi − µ)> − I ∥∥∥∥∥∥ 2 ≤ 1 2 γ + 1 2 γ = γ .
  Notice that 1 |T | ∑ i∈T (Xi − µ(S))(Xi − µ(S))>
  = 1 |T | ∑ i∈T (Xi − µ)(Xi − µ)> + (µ− µ(S))(µ− µ(S))>
  + (µ− µ(S))
  ( 1
  |T | ∑ i∈T Xi − µ
  )T + ( 1
  |T | ∑ i∈T Xi − µ
  ) (µ− µ(S))T .
  Combining the above two inequalities, by the triangle inequality, we get∥∥∥∥∥ 1|T |∑ i∈T (Xi − µ(S))(Xi − µ(S))> − I ∥∥∥∥∥ 2 ≤ γ + ‖µ− µ(S)‖22 + 2β ‖µ− µ(S)‖2
  = γ +
  3β2
  2 .
  When S is ( , β, γ)-stable, by Lemma C.6, we know that S is ( , γ1 = β, γ2 = γ, β1 = β/ , β2 = γ/ + 3β2/ 2)-good. The parameter ξ in Lemma C.7 is
  ξ = 3γ + 6β2 + 6β2/ ,
  and error guarantee in Lemma C.7 translates to √ ξ + γ1 = O( √ γ + β), which is exactly what is needed in Lemma C.6.
sections.heading: D OMITTED PROOFS FROM SECTION 4
sections.text: In this section, we prove the technical lemmas in Section 4. We restate each lemma before proving it.
  Lemma 4.2 states that the (scaled) initial estimation is not too far from the true conditional probabilities p.
  Lemma 4.2. Consider the same setting as in Theorem 4.1. Assume the conditions in Section 2.3 hold. In Algorithm 1, we have ∥∥πP ◦ (p̂− q̂0)∥∥
  2 ≤ O(
  √ d/ √ αc) .
  Proof. Recall that q0 = qS is the empirical conditional probabilities over S, and v̂ = v ◦ s where s is the scaling vector with ‖s‖∞ ≤ O(1/ √ αc).
  Let πG ? and pG ?
  be the empirical parental configuration probabilities and conditional probabilities given by G?. We first show that ∥∥∥πG? − πS∥∥∥ 2 ≤ √ 2d .
  Let nG ?
  k and n S k denote the number of times that Πk happens in G ? and S. Note that changing one sample in G? can increase or decrease nG ?
  k by at most 1. Moreover, in a single sample, exactly d parental configuration events happen, so changing a sample can affect at most 2d nG ?
  k ’s. Since S is obtained from G? by changing N samples, we have |nG?k − nSk | ≤ N for all k, and ∑ k |nG ?
  k − nSk | ≤ 2 dN . Together they imply ∥∥πG? − πS∥∥ 2 ≤ √ 2d. By a similar argument, we can show that∥∥∥πG? ◦ pG? − πS ◦ qS∥∥∥ 2 ≤ √ 2d ,
  because πG ? k p G? k is the probability that Πk happens and X(k) = 1 over G ?.
  By the triangle inequality, we have∥∥∥πG? ◦ (pG? − qS)∥∥∥ 2 ≤ ∥∥∥πG? ◦ pG? − πS ◦ qS∥∥∥ 2 + ∥∥∥πG? − πS∥∥∥ 2 ≤ 3 √ d . Using the condition in Equation equation 1 from Section 2.3, i.e., ∥∥πG? ◦ (pG? − p)∥∥
  2 ≤ O( ), we
  get ∥∥∥πG? ◦ (p− qS)∥∥∥ 2 ≤ O( √ d) .
  Now by Equation equation 2 from Section 2.3 and the assumption that the minimum parental configuration probability mink πPk = α = Ω( ), we have π P k ≤ πG ? k +O( ) ≤ O(πG ?
  k ), and hence∥∥πP ◦ (p− qS)∥∥ 2 ≤ O( √ d) .
  After scaling by s, we have ∥∥∥πG? ◦ (p̂− q̂S)∥∥∥ 2 ≤ O( √ d/ √ αc) . Lemma 4.3 shows that, when q is relatively far from p, the algorithm can find a new q such that∥∥πP ◦ (p̂− q̂)∥∥ 2 decreases by a constant factor.
  Lemma 4.3. Consider the same setting as in Theorem 4.1. Assume the conditions in Section 2.3 hold. Fix an iteration t in Algorithm 1. Assume the robust mean estimation algorithm Amean succeeds. If∥∥πP ◦ (p̂− q̂t)∥∥
  2 ≤ ρt and ρt = Ω(
  √ log(1/ )/
  √ αc), then we have∥∥πP ◦ (p̂− q̂t+1)∥∥
  2 ≤ c1ρt
  for some universal constant c1 < 1. Proof. We assume ρt > c4( √ log(1/ )/ √ αc) and α > c5 for some sufficiently large universal constants c4 and c5.
  Because ∥∥πP ◦ (p̂− q̂t)∥∥
  2 ≤ ρt, Lemma 3.2 shows that
  { f̂(Xi, q t) } i∈G?
  is( , O ( √
  log 1/ √ αc + α ρt
  ) , O ( log 1/
  αc +
  (ρt)2
  α + ρt√ α
  )) -stable.
  By Lemma 3.3, the robust mean estimation oracle Amean, which we assume to succeed, outputs a ν ∈ Rm such that, for some universal constant c3,∥∥ν − πP ◦ (p̂− q̂t)∥∥
  2 ≤ c3
  (√
  α ρt + √ √ α ρt + α ρt + √ log(1/ )√ αc )
  < ( c3√ c5 + c3√ c4 + c3 c5 + c3 c4 ) ρt .
  From Section 2.3, we have ∥∥πS − πP∥∥∞ = O( ), which implies∥∥(πS − πP ) ◦ (p̂− q̂t)∥∥
  2 ≤ α ∥∥πP ◦ (p̂− q̂t)∥∥ 2 ≤ α ρt .
  By the triangle inequality, we have∥∥ν − πS ◦ (p̂− q̂t)∥∥ 2 ≤ ( c3√ c5 + c3√ c4 + c3 + 1 c5 + c3 c4 ) ρt .
  Algorithm 1 sets q̂t+1 = ν ◦ (1/πS) + q̂t, which is equivalent to
  πS ◦ (p̂− q̂t+1) = πS ◦ (p̂− q̂t)− ν .
  Since ∥∥πS − πP∥∥∞ = O( ) and α = Ω( ), we have
  πPi ≤ 1.1πSi ∀1 ≤ i ≤ m .
  Putting everything together, letting c1 = 1.1 ( c3√ c5 + c3√c4 + c3+1 c5 + c3c4 ) , we have ∥∥πP ◦ (p̂− q̂t+1)∥∥ 2 ≤ 1.1 ∥∥πS ◦ (p̂− q̂t+1)∥∥ 2 < 1.1 ( c3√ c5 + c3√ c4 + c3 + 1 c5 + c3 c4 ) ρt = c1ρ t .
  Because c4 and c5 can be sufficiently large, we have c1 < 1 as needed.
  Lemma 4.4 shows that when the algorithm terminates, we can conclude that the output Q is close to the ground-truth P in total variation distance.
  Lemma 4.4. Consider the same setting as in Theorem 4.1. Assume the conditions in Section 2.3 hold. Let Q be a Bayesian network that shares the same structure with P . Suppose that (1) P is c-balanced, (2) α = Ω(r + /c), and (3) ∥∥πP ◦ (p̂− q̂)∥∥ 2 ≤ r/2. Then we have
  dTV (P,Q) ≤ r .
  Proof of Lemma 4.4. We have (pk + qk)(2− pk − qk) ≥ pk(1− pk). Hence,∑ k √ πPk π Q k (pk − qk)2 (pk + qk)(2− pk − qk) ≤ ∑ k √ πPk π Q k π P k (pk − qk)2 πPk pk(1− pk) ,
  From the proof of Lemma B.4, we know |πPk pk(1−pk)− 1s2k | = O( ) and π P k pk(1−pk) ≥ πPk pk 2 = Ω(r), so we have∑ k √ πPk π Q k π P k (pk − qk)2 πPk pk(1− pk) ≤ 1.1 ∑ k √ πPk π Q k π P k (pk − qk)2s2k
  = 1.1 ∑ k √ πPk π Q k π P k (p̂k − q̂k)2 .
  It suffices to show that |πPk − π Q k | ≤ r, which implies π Q k ≤ 1.1πPk and further implies dTV (P,Q) ≤ 2 ∥∥πP ◦ (p̂− q̂)∥∥
  2 .
  Let P≤i and Q≤i be the distributions of the first i coordinates of P and Q respectively. We prove |πPk − π Q k | ≤ r by induction on i. Suppose that for 1 ≤ j < i and all a′ ∈ {0, 1}|parents(j)|, |πPj,a′ − π Q j,a′ | ≤ r, then we have dTV (P≤(i−1), Q≤(i−1)) ≤ r. Because that events Πi,a only depends on j < i, |πPi,a − π Q i,a| ≤ dTV (P≤(i−1), Q≤(i−1)) ≤ r for all a. Consequently, we have dTV (P,Q) = dTV (P≤d, Q≤d) ≤ r.