{'fraction_non_alphanumeric': 0.035220315924722245, 'fraction_numerical': 0.012999773259768725, 'mean_word_length': 4.348221503637833, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interaction Screening: Efficient and Sample-Optimal Learning of Ising Models

We consider the problem of learning the underlying graph of an unknown Ising model on p spins from a collection of i.i.d. samples generated from the model. We suggest a new estimator that is computationally efficient and requires a number of samples that is near-optimal with respect to previously established information-theoretic lower-bound. Our statistical estimator has a physical interpretation in terms of"interaction screening". The estimator is consistent and is efficiently implemented using convex optimization. We prove that with appropriate regularization, the estimator recovers the underlying graph using a number of samples that is logarithmic in the system size p and exponential in the maximum coupling-intensity and maximum node-degree.


Introduction
A Graphical Model (GM) describes a probability distribution over a set of random variables which factorizes over the edges of a graph. It is of interest to recover the structure of GMs from random samples. The graphical structure contains valuable information on the dependencies between the random variables. In fact, the neighborhood of a random variable is the minimal set that provides us maximum information about this variable. Unsurprisingly, GM reconstruction plays an important role in various fields such as the study of gene expression [1], protein interactions [2], neuroscience [3], image processing [4], sociology [5] and even grid science [6,7].
The origin of the GM reconstruction problem is traced back to the seminal 1968 paper by Chow and Liu [8], where the problem was posed and resolved for the special case of tree-structured GMs. In this special tree case the maximum likelihood estimator is tractable and is tantamount to finding a maximum weighted spanning-tree. However, it is also known that in the case of general graphs with cycles, maximum likelihood estimators are intractable as they require computation of the partition function of the underlying GM, with notable exceptions of the Gaussian GM, see for instance [9], and some other special cases, like planar Ising models without magnetic field [10].
A lot of efforts in this field has focused on learning Ising models, which are the most general GMs over binary variables with pairwise interaction/factorization. Early attempts to learn the Ising model structure efficiently were heuristic, based on various mean-field approximations, e.g. utilizing empirical correlation matrices [11,12,13,14]. These methods were satisfactory in cases when correlations decrease with graph distance. However it was also noticed that the mean-field methods perform poorly for the Ising models with long-range correlations. This observation is not surprising in light of recent results stating that learning the structure of Ising models using only their correlation matrix is, in general, computationally intractable [15,16].
Among methods that do not rely solely on correlation matrices but take advantage of higherorder correlations that can be estimated from samples, we mention the so-called regularized pseudolikelihood estimator [17]. This estimator, like the one we propose in this paper, is from the class of M-estimators i.e. estimators that are the minimum of a sum of functions over the sampled data. The regularized pseudo-likelihood estimator can be interpreted as a surrogate for the intractable likelihood estimator with an additive 1 -norm penalty to encourage sparsity of the reconstructed graph. The estimator offers guarantees for the structure reconstruction, but the result only applies to GMs that satisfy a certain condition that is rather restrictive and hard to verify. It was also proven that this estimator fails to reconstruct the structure of graphs with long-range correlations, even for simple test cases [18].
Principal tractability of structure reconstruction of an arbitrary Ising model from samples was proven only very recently. Bresler, Mossel and Sly in [19] suggested an algorithm which reconstructs the graph without errors in polynomial time. They showed that the algorithm requires number of samples that is logarithmic in the number of variables. Although this algorithm is of a polynomial complexity, it relies on an exhaustive neighborhood search, and the degree of the polynomial is equal to the maximal node degree.
Prior to the work reported in this manuscript the best known procedure for perfect reconstruction of an Ising model was through a greedy algorithm proposed by Bresler in [20]. Bresler's algorithm is based on the observation that the mutual information between neighboring nodes in an Ising model is lower bounded. This observation allows to reconstruct the Ising graph perfectly with only a logarithmic number of samples and in time quasi-quadratic in the number of variables. On the other hand, Bresler's algorithm suffers from two major practical limitations. First, the number of samples, hence the running time as well, scales double exponentially with respect to the largest node degree and with respect to the largest coupling intensity between pairs of variables. This scaling is rather far from the information-theoretic lower-bound reported in [21] predicting instead a single exponential dependency on the two aforementioned quantities. Second, Bresler's algorithm requires prior information on the maximum and minimum coupling intensities as well as on the maximum node degree, guarantees which, in reality, are not necessarily available.
In this paper we propose a novel estimator for the graph structure of an arbitrary Ising model which achieves perfect reconstruction in quasi-quartic time (although we believe it can be reduced to quasi-quadratic time) and with a number of samples logarithmic in the system size. The algorithm is near-optimal in the sense that the number of samples required to achieve perfect reconstruction, and the run time, scale exponentially with respect to the maximum node-degree and the maximum coupling intensity, thus matching parametrically the information-theoretic lower bound of [21]. Our statistical estimator has the structure of a consistent M-estimator [22] and is implemented via convex optimization. Moreover it allows intuitive interpretation in terms of what we coin the "interaction screening". We show that with a proper 1 -regularization our estimator reconstructs couplings of an Ising model from a number of samples that is near-optimal (in the sense of comparison with the information-theoretic bound). In addition, our estimator does not rely on prior information on the model characteristics, such as maximum coupling intensity and maximum degree.
The rest of the paper is organized as follows. In Section 2 we give a precise definition of the structure estimation problem for the Ising models and we describe in detail our method for structure reconstruction within the family of Ising models. The main results related to the reconstruction guarantees are provided by Theorem 1 and Theorem 2. In Section 3 we explain the strategy and the sequence of steps that we use to prove our main theorems. Proofs of Theorem 1 and Theorem 2 are summarized at the end of this Section. Section 4 illustrates performance of our reconstruction algorithm via simulations. Here we show on a number of test cases that the sample complexity of the suggested method scales logarithmically with the number of variables and exponentially with the maximum coupling intensity. In Section 5 we discuss possible generalizations of the algorithm and future work.

Main Results
Consider a graph G = (V, E) with p vertexes where V = {1, . . . , p} is the vertex set and E ⊂ V ×V is the undirected edge set. Vertexes i ∈ V are associated with binary random variables σ i ∈ {−1, +1} that are called spins. Edges (i, j) ∈ E are associated with non-zero real parameters θ * ij = 0 that are called couplings. An Ising model is a probability distribution µ over spin configurations σ = {σ 1 , . . . , σ p } that reads as follows: where Z is a normalization factor called the partition function.
Notice that even though the main innovation of this paper -the efficient "interaction screening" estimator -can be constructed for the most general Ising models, we restrict our attention in this paper to the special case of the Ising models with zero magnetic field. This simplification is not necessary and is done solely to simplify (generally rather bulky) algebra. Later in the text we will thus refer to the zero magnetic field model (2) simply as the Ising model.

Structure-Learning of Ising Models
Suppose that n sequences/samples of p spins σ (k) k=1,...,n are observed. Let us assume that each observed spin configuration (1). Based on these measurements/samples we aim to construct an estimator E of the edge set that reconstructs the structure exactly with high probability, i.e.
where ∈ 0, 1 2 is a prescribed reconstruction error. We are interested to learn structures of Ising models in the high-dimensional regime where the number of observations/samples is of the order n = O (ln p). A necessary condition on the number of samples is given in [21,Thm. 1]. This condition depends explicitly on the smallest and largest coupling intensity α := min and on the maximal node degree where the set of neighbors of a node i ∈ V is denoted by ∂i According to [21], in order to reconstruct the structure of the Ising model with minimum coupling intensity α, maximum coupling intensity β, and maximum degree d, the required number of samples should be at least We see from Eq. (6) that the exponential dependence on the degree and the maximum coupling intensity are both unavoidable. Moreover, when the minimal coupling is small, the number of samples should scale at least as α −2 . It remains unknown if the inequality (6) is achievable. It is shown in [21,Thm. 3] that there exists a reconstruction algorithm with error probability ∈ 0, 1 2 if the number of samples is greater than Unfortunately, the existence proof presented in [21] is non-constructive and thus it does not guarantee actual existence of an algorithm with low computational complexity. Notice also that the number of samples in (7) scales as exp (4βd) when d and β are asymptotically large and as α −4 when α is asymptotically small.

Regulrized Interaction Screening Estimator
The main contribution of this paper consists in presenting explicitly a structure-learning algorithm that is of low complexity and which is near-optimal with respect to bounds (6) and (7). Our algorithm reconstructs the structure of the Ising model exactly, as stated in Eq. (3), with an error probability ∈ 0, 1 2 , and with a number of samples which is proportional to exp (6βd) and α −2 . (See Theorem 1 and Theorem 2 below for mathematically accurate statements.) Our algorithm consists of two steps. First, we estimate couplings in the vicinity of every node. Then, on the second step, we threshold the estimated couplings that are sufficiently small to zero. Resulting zero coupling means that the corresponding edge is not present.
Denote the set of couplings around node u ∈ V by the vector θ * u ∈ R p−1 . In this, slightly abusive notation, we use the convention that if a coupling is equal to zero it reads as absence of the edge, i.e. θ * ui = 0 if and only if (u, i) / ∈ E. Note that if the node degree is bounded by d, it implies that the vector of couplings θ * u is non-zero in at most d entries. Our estimator for couplings around node u ∈ V is based on the following loss function coined the Interaction Screening Objective (ISO): The objective of the ISO is inversely proportional to the corresponding factor in the distribution (1). As a result, the effect of the objective is to reweight the distribution in a way that the spin σ u becomes independent of the rest; i.e., the dependencies of σ u on the rest of the spins are completely "screened", which explains our choice for the name of the loss function. This remarkable "screening" feature of the ISO suggests the following choice of the Regularized Interaction Screening Estimator (RISE) for the interaction vector around node u: where λ > 0 is a tunable parameter promoting sparsity through the additive 1 -penalty. Notice that the ISO is the empirical average of an exponential function of θ u which implies it is convex. Moreover, addition of the 1 -penalty preserves the convexity of the minimization objective in Eq. (9). As expected, the performance of RISE does depend on the choice of the penalty parameter λ. If λ is too small θ u (λ) is too sensitive to statistical fluctuations. On the other hand, if λ is too large θ u (λ) has too much of a bias towards zero. In general, the optimal value of λ is hard to guess. Luckily, the following theorem provides strong guarantees on the square error for the case when λ is chosen to be sufficiently large.
k=1,...,n be n realizations of p spins drawn i.i.d. from an Ising model with maximum degree d and maximum coupling intensity β. Then for any node u ∈ V and for any 1 > 0, the square error of the Regularized Interaction Screening Estimator (9) with penalty parameter λ ≥ 4 ln(3p/ 1) n is bounded with probability at least 1 − 1 by whenever n ≥ (9d) 5 e 6βd ln 3p 1 . Our structure estimator (for the second step of the algorithm), Structure-RISE, takes RISE output and thresholds couplings whose absolute value is less than α/2 to zero: Performance of the Structure-RISE is fully quantified by the following Theorem.
Proofs of Theorem 1 and Theorem 2 are given in Subsection 3.3. Theorem 1 states that RISE recovers not only the structure but also the correct value of the couplings up to an error based on the available samples. It is possible to improve the square-error bound (10) even further by first, running the Structure-RISE to recover edges, and then re-running RISE with λ = 0 for the remaining non-zero couplings.
The computational complexity of RISE is tantamount to minimizing the convex objective function and, as such, it scales at most as O np 3 . Therefore, computational complexity of the Structure-RISE scales at most as O np 4 simply because one has to call RISE at every node. We believe that this running time estimate can actually be reduced to a quasi-quadratic one by using first-order minimization techniques, in the spirit of [23].
Notice that in order to implement RISE there is no need for prior knowledge on the graph parameters. This is a considerable advantage in practical applications where the maximum degree or bounds on couplings is often unknown.

Analysis
The Regularized Interaction Screening Estimator (9) is from the class of the so-called regularized M-estimators. Naghaband et al. proposed in [22] a framework to analyze the square error of such estimators. As per [22], enforcing only two conditions on the loss function is sufficient to get a handle on the square error of an 1 -regularized M-estimator.
The first condition links the choice of the penalty parameter to the gradient of the objective function.
Condition 1. The 1 -penalty parameter strongly enforces regularization if it is greater than any partial derivatives of the objective function at θ u = θ * u , i.e.
Condition 1 guarantees that if the vector of couplings θ * u has at most d non-zero entries, then the estimation difference θ u (λ) − θ * u lies within the set The second condition ensure that the objective function is strongly convex in a restricted subset of R p−1 . Denote the reminder of the first-order Taylor expansion of the objective function by where ∆ u ∈ R p−1 is an arbitrary vector. Then the second condition reads as follows.

Condition 2.
The objective function is strongly convex at θ u = θ * u restricted to K on a neighborhood of radius R, if for all ∆ u ∈ K such that ∆ u 2 ≤ R, there exists a constant κ > 0 such that Strong regularization and restricted strong convexity enables us to control that the minimizer θ u of the full objective (9) lies in the vicinity of the sparse vector of parameters θ * u . The precise formulation is given in the proposition following from [22, Thm. 1]. Proposition 1. If the 1 -regularized M-estimator of the form (9) satisfies Condition 1 and Condition 2 with R > 3d λ κ then the square-error is bounded by

Gradient Concentration
Like the ISO (8), its gradient in any component l ∈ V \ u is an empirical average where the random variables X In order to prove that the ISO gradient concentrates we have to state few properties of the support, the mean and the variance of the random variables (19), expressed in the following three Lemmas. The first of the Lemmas states that at θ u = θ * u , the random variable X ul (θ * u ) has zero mean.

Lemma 1. For any Ising model with p spins and for all
Proof. By direct computation, we find that where in the last line we use the fact that the exponential terms involving σ u cancel, implying that the sum over σ u ∈ {−1, +1} is zero.
As a direct corollary of the Lemma 1, θ u = θ * u is always a minimum of the averaged ISO (8). The second Lemma proves that at θ u = θ * u , the random variable X ul (θ * u ) has a variance equal to one.

Lemma 2. For any Ising model with p spins and for all
Proof. As a result of direct evaluation one derives Notice that in the second line the first sum over edges (under the exponential) does not depend on σ u . Furthermore, the first sum is invariant under the change of variables, σ u → −σ u , while the second sum changes sign. This transformation results in appearance of the partition function in the numerator.
The next lemma states that at θ u = θ * u , the random variable X ul (θ * u ) has a bounded support.
Lemma 3. For any Ising model with p spins, with maximum degree d and maximum coupling intensity β, we guarantee that for all l = u ∈ V Proof. Observe that components of θ * u are smaller than β and at most d of them are non-zero. Recall that spins are binary, {−1, +1}, which results in the following estimate With Lemma 1, 2 and 3, and using Berstein's inequality we are now in position to prove that every partial derivative of the ISO concentrates uniformly around zero as the number of samples grows.
Proof. Let us first show that every term is individually bounded by the RHS of (26) with highprobability. We further use the union bound to prove that all components are uniformly bounded with high-probability. Utilizing Lemma 1, Lemma 2 and Lemma 3 we apply the Bernstein's Inequality Inverting the following relation and substituting the result in the Eq. (27) one derives where u = s n exp (βd) . When n ≥ s exp (2βd) Eq. (29) can be simplified to become independent of β and d Using s = ln 2p 3 and the union bound on every component of the gradient leads to the desired result.

Restricted Strong-Convexity
The remainder of the first-order Taylor-expansion of the ISO, defined in Eq. (15) reads where the random variables Y In order to prove that the reminder of Eq. (31) is lower bounded by a quadratic form we show a few properties of the relevant random variables in (32). The following Lemma provides a lower bound on the second moment of (32).

Lemma 5.
Consider an Ising model with p spins, with maximum degree d and maximum coupling intensity β. For all ∆ u ∈ R p−1 the following bound holds Proof. Our proof strategy here follows [16,Cor. 3.1]. Notice that the probability measure of the Ising model is symmetric with respect to the sign flip, i.e. µ (σ 1 , . . . , σ p ) = µ (−σ 1 , . . . , −σ p ). Thus any spin has zero mean, which implies that for every ∆ u ∈ R p−1 This allows to reinterpret (33) as a variance, using that σ 2 u = 1, Construct a subset A ⊂ V recursively as follows: (i) let i 0 = argmax j∈V \u ∆ 2 uj and define The set A possesses the following two main properties. First, every node i ∈ A does not have any neighbors in A and, second, We apply the law of total variance to (35) by conditioning on the set of spins σ A c whose indexes are from the complementary set A c .

Var
where in the last line one uses that the spins in A are conditionally independent given their neighbors σ A c . One concludes the proof by using relation (36) and the fact that the conditional variance of a spin given its neighbors are bounded from below: The following lemma proves a concentration property for the remainder of the Taylor expansion.
Lemma 6. Consider an Ising model with p spins, with maximum degree d and maximum coupling intensity β. For all ∆ u ∈ R p−1 and for all 4 > 0, the remainder of the Taylor expansion (31) satisfies with probability at least 1 − 4 , the following inequality whenever n ≥ exp (2βd) ln 1 4 . Proof. This concentration property is based on the Bernstein's inequality. First of all, observe that for all z ∈ R, the following bound holds This implies that the remainder of the Taylor expansion (31) is lower-bounded Using similar arguments than in Lemma 1 it is easy to show that for all ∆ ∈ R p−1 Then extending the technique of the type used in Lemma 2, one shows that for every integer l ∈ N Notice also that the support of Y u (∆ u ) is trivially upper bounded for all ∆ u ∈ R p−1 To finish the proof we apply the Bernstein's inequality to the right-hand side of (41), thus combining relations (42), (43), (44) and Lemma 5. Moreover when n ≥ e 2βd ln 1 4 further simplifications can be made in the way similar the one used in Lemma 4.
We are now in the position to prove that the ISO satisfies the the restricted strong convexity condition.

Proof of the main Theorems
Proof of Theorem 1 (Square Error of RISE). We seek to apply Proposition 1 to the Regularized Interaction Screening Estimator (9). Using 3 = 2 1 3 in Lemma 4 and letting λ ≥ 4 ln 3p/ 1 n , it follows that Condition 1 is satisfied with probability greater than 1 − 2 1 /3, whenever n ≥ (9d) 5 e 6βd ln 3p 1 . Using 4 = 1 /3 in Lemma 7, and observing that 3dλ e −3βd for n ≥ (9d) 5 e 6βd ln p 1 , we conclude that condition 2 is satisfied with probability greater than 1 − 1 3 . Theorem 1 then follows by using a union bound and then applying Proposition 1. The proof of Theorem 2 becomes an immediate application of Theorem 1.
Proof of Theorem 2 (Structure Learning of Ising Models). According to Theorem 1, one observes that, with probability 1 − 1 , the minimal amount of samples required to achieve an error of α/2 on every couplings around a single node is Let us choose 2 = 1 /p and use the union-bound to ensure that the couplings at every node (thresholded by α/2) is simultaneously recovered with probability greater than 1 − 2 .

Numerical Results
We test performance of the Struct-RISE, with the strength of the l 1 -regularization parametrized by λ = 4 ln(3p 2 / ) n , on Ising models over two-dimensional grid with periodic boundary conditions (thus degree of every node in the graph is 4).
We are interested to find the minimal number of samples, n min , such that the graph is perfectly reconstructed with probability 1 − ≥ 0.95. In our numerical experiments, we recover the value of n min as the minimal n for which Struct-RISE outputs the perfect structure 45 times from 45 different trials with n samples, thus guaranteeing that the probability of perfect reconstruction is greater than 0.95 with a statistical confidence of ≈ 90%.
We first verify the logarithmic scaling of n min with respect to the number of spins p. The couplings are chosen uniform and positive θ * ij = 0.7. This choice ensures that samples generated by Glauber dynamics are i.i.d. according to (1). Values of n min for p ∈ {9, 16, 25, 36, 49, 64} are shown on the left in Figure 1. Empirical scaling is, ≈ 1.1 × 10 5 ln p, which is orders of magnitude better than the rather conservative prediction of the theory, 4.9 × 10 15 ln p.
We also test the exponential scaling of n min with respect to the maximum coupling intensity β. The test is conducted over two different settings both with p = 16 spins: the ferromagnetic case where all couplings are uniform and positive, and the spin glass case where the sign of couplings is assigned uniformly at random. In both cases the absolute value of the couplings, θ * ij , is uniform and equal to β. To ensure that the samples are i.i.d, we sample directly from the exhaustive weighted list of the 2 16 possible spin configurations.
Experimental values of n min for different values of the maximum coupling intensity, β, are shown on the right in Fig. 1. Empirically observed exponential dependence on β is matched best by, Linear-exponential plot showing the observed relation between n min and β. The graph is the twodimensional 4 × 4 grid. In red the couplings are uniform and positive and in blue the couplings have uniform intensity but random sign.
exp (12.8β), in the ferromagnetic case and by, exp (5.6β), in the case of the spin glass. Theoretical bound for d = 4 predicts exp (24β). We observe that the difference in sample complexity depends significantly on the type of interaction. An interesting observation one can make based on these experiments is that the case which is harder from the sample-generating perspective is easier for learning and vice verse.

Conclusions and Path Forward
In this paper we construct and analyze the Regularized Interaction Screening Estimator (9). We show that the estimator is computationally efficient and needs an optimal number of samples for learning Ising models. The RISE estimator does not require any prior knowledge about the model parameters for implementation and it is based on the minimization of the loss function (8), that we call the Interaction Screening Objective. The ISO is an empirical average (over samples) of an objective designed to screen an individual spin/variable from its factor-graph neighbors. Even though we focus in this paper solely on learning pair-wise binary models, the "interaction screening" approach we introduce here is generic. The approach extends to learning other Graphical Models, including those over higher (discrete, continuous or mixed) alphabets and involving highorder (beyond pair-wise) interactions. These generalizations are built around the same basic idea pioneered in this paper -the interaction screening objective is (a) minimized over candidate GM parameters at the actual values of the parameters we aim to learn; and (b) it is an empirical average over samples. In the future, we plan to explore further theoretical and experimental power, characteristics and performance of the generalized screening estimator.