Ising Model Selection Using `1-Regularized Linear Regression: A Statistical Mechanics Analysis

Abstract:
We theoretically analyze the typical learning performance of `1-regularized linear regression (`1-LinR) for Ising model selection using the replica method from statistical mechanics. For typical random regular graphs in the paramagnetic phase, an accurate estimate of the typical sample complexity of `1-LinR is obtained. Remarkably, despite the model misspecification, `1-LinR is model selection consistent with the same order of sample complexity as `1-regularized logistic regression (`1-LogR), i.e., M = O (logN), where N is the number of variables of the Ising model. Moreover, we provide an efficient method to accurately predict the nonasymptotic behavior of `1-LinR for moderate M,N , such as precision and recall. Simulations show a fairly good agreement between theoretical predictions and experimental results, even for graphs with many loops, which supports our findings. Although this paper mainly focuses on `1-LinR, our method is readily applicable for precisely characterizing the typical learning performances of a wide class of `1-regularized M -estimators including `1-LogR and interaction screening.

heading: 1 Introduction
text: The advent of massive data across various scientific disciplines has led to the widespread use of undirected graphical models, also known as Markov random fields (MRFs), as a tool for discovering and visualizing dependencies among covariates in multivariate data [1]. The Ising model, originally proposed in statistical physics, is one special class of binary MRFs with pairwise potentials and has been widely used in different domains such as image analysis, social networking, gene network analysis [2, 3, 4, 5, 6, 7]. Among various applications, one fundamental problem of interest is called Ising model selection, which refers to recovering the underlying graph structure of the original Ising model from independent, identically distributed (i.i.d.) samples. A variety of methods have been proposed [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], demonstrating the possibility of successful Ising model selection even when the number of samples is smaller than that of the variables. Notably, it has been demonstrated that for the `1-regularized logistic regression (`1-LogR) [10, 16] and interaction screening (IS) [14, 15] estimators, M = O (logN) samples suffice for an Ising model with N spins under certain assumptions, which is consistent with respect to (w.r.t.) previously established ∗Corresponding author.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
information-theoretic lower-bound [11]. Both `1-LogR and IS are `1-regularized M -estimators [19] with logistic and IS objective (ISO) loss functions, respectively.
In this paper, we focus on one simpler linear estimator called `1-regularized linear regression (`1- LinR) and theoretically investigate its typical learning performance using the powerful replica method [20, 21, 22, 23] from statistical mechanics. The `1-LinR estimator, widely known as least absolute shrinkage and selection operator (LASSO) [24] in statistics and machine learning, is considered here mainly for two reasons. On the one hand, it is one representative example of model misspecification since the quadratic loss of `1-LinR does not match the true log-conditional-likelihood as `1-LogR, nor does it have the interaction screening property as IS. On the other hand, as one of the most popular linear estimator, `1-LinR is more computationally efficient than `1-LogR and IS, and thus it is of practical importance to investigate its learning performance for Ising model selection. Since it is difficult to obtain results for general graphs, as a first step we consider the random regular (RR) graphs GN,d,K0 in the paramagnetic phase [23], where GN,d,K0 denotes the ensemble of RR graphs with constant node degree d and uniform coupling strength K0 on the edges.

heading: 1.1 Contributions
text: The main contributions are summarized as follows. First, we obtain an accurate estimate of the typical sample complexity of `1-LinR for Ising model selection for typical RR graphs in the paramagnetic phase, which, remarkably, has the same order as `1-LogR. Specifically, for a typical RR graph G ∈ GN,d,K0 , using `1-LinR with a regularization parameter 0 < λ < tanh (K0), one can consistently reconstruct the structure with M > c(λ,K0) logNλ2 samples, where c (λ,K0) = 2(1−tanh2(K0)+dλ2) 1+(d−1) tanh2(K0) . The accuracy of our typical sample complexity prediction is verified by its excellent agreement with experimental results. To the best of our knowledge, this is the first result that provides an accurate typical sample complexity for Ising model selection. Interestingly, as λ→ tanh (K0), a lower bound M > 2 logN
tanh2(K0) of the typical sample complexity is obtained, which has the same scaling as the
information-theoretic lower bound M > c ′ logN K20 [11] for some constant c′ at high temperatures (i.e., small K0) since tanh (K0) = O (K0) as K0 → 0. Second, we provide a computationally efficient method to precisely predict the typical learning performance of `1-LinR in the non-asymptotic case with moderate M,N , such as precision, recall, and residual sum of square (RSS). Such precise non-asymptotic predictions of `1-LinR for Ising model selection have been unavailable even for `1-LogR [10, 16] and IS [14, 15], nor are they the same as previous asymptotic results of `1-LinR assuming fixed α ≡M/N [25, 26, 27, 28]. Moreover, although our theoretical analysis is based on a tree-like structure assumption, experimental results on two dimensional (2D) grid graphs also show a fairly good agreement, indicating that our theoretical result can be a good approximation even for graphs with many loops.
Third, while this paper mainly focuses on `1-LinR, our method is readily applicable to a wide class of `1-regularized M -estimators [19], including `1-LogR [10] and IS [14, 15]. Thus, an additional technical contribution is providing a generic approach for precisely characterizing the typical learning performances of various `1-regularized M -estimators for Ising model selection. Although the replica method from statistical mechanics is non-rigorous, our results are conjectured to be correct, which is supported by their excellent agreement with the experimental results. Additionally, several technical advances we propose in this paper, e.g., the entropy term computation by averaging over the Haar measure and the modification of EOS to address the finite-size effect, might be of general interest to those who use the replica method as a tool for performance analysis.

heading: 1.2 Related works
text: There has been some earlier works on the analysis of Ising model selection (also known as the inverse Ising problem) using the replica method [4, 5, 6, 7, 29] from statistical mechanics. For example, in [6], the performance of the pseudo-likelihood (PL) method [30] is studied. However, instead of graph structure learning, [6] focuses on the problem of parameter learning. Then, [7] extends the analysis to the Ising model with sparse couplings using logistic regression without regularization. The recent work [29] analyzes the performance of `2-regularized linear regression but the techniques invented there are not applicable to `1-LinR since the `1-norm breaks the rotational invariance property.
Regarding the study of `1-LinR (LASSO) under model misspecification, the past few years have seen a line of research in the field of signal processing with a specific focus on the single-index
model [31, 32, 27, 33, 34]. These studies are closely related to ours but there are several important differences. First, in our study, the covariates are generated from an Ising model rather than a Gaussian distribution. Second, we focus on model selection consistency of `1-LinR while most previous studies consider estimation consistency except [33]. However, [33] only considers the classical asymptotic regime while our analysis includes the high-dimensional setting where M N . As far as we have searched, there is no earlier study of `1-LinR estimator for Ising model selection, though some are found for Gaussian graphical models [35, 36]. One closely related work [15] states that at high temperatures when the coupling magnitude is approaching zero, both logistic and ISO losses can be approximated by a quadratic loss. However, their claim is only restricted to the very small magnitude near zero while our analysis extends the validity range to the whole paramagnetic phase. Moreover, they evaluate the minimum number of samples necessary for consistently reconstructing “arbitrary” Ising models, which, however, seems much larger than that actually needed. By contrast, we provide the first accurate assessment of typical sample complexity for consistently reconstructing typical samples of Ising models defined over the RR graphs. Furthermore, [15] does not provide precise predictions of the non-asymptotic learning performance as we do.

heading: 2 Background and Problem Setup

heading: 2.1 Ising Model
text: Ising model is one special class of MRFs with pairwise potentials and each variable takes binary values [22, 23], which is one classical model from statistical physics. The joint probability distribution of an Ising model with N variables (spins) s = (si) N−1 i=0 ∈ {−1,+1} N has the form
PIsing (s|J∗) = 1
ZIsing (J ∗)
exp ∑ i<j J∗ijsisj  , (1) where ZIsing (J∗) = ∑ s exp {∑ i<j J ∗ ijsisj } is the partition function and J∗ = ( J∗ij ) i,j
are the original couplings, respectively. In general, there are also external fields but here they are assumed to be zero for simplicity. The structure of Ising model can be described by an undirected graph G = (V, E), where V = {0, 1, ..., N − 1} is a collection of vertices at which the spins are assigned, and E = { (i, j) |J∗ij 6= 0 } is a collection of undirected edges, i.e., J∗ij = 0 for all pairs of (i, j) /∈ E. For each vertex i ∈ V, its neighborhood is defined as the subset N (i) ≡ {j ∈ V| (i, j) ∈ E}.
2.2 Neighborhood-based `1-regularized linear regression (`1-LinR)
The problem of Ising model selection refers to recovering the graph G (edge set E), given M i.i.d. samples DM = { s(1), ..., s(M) } from the Ising model. While the maximum likelihood method has nice properties of consistency and asymptotic efficiency, it suffers from high computational complexity. To tackle this difficulty, several local learning algorithms have been proposed, notably the `1-LogR estimator [10] and IS estimator [14]. Both `1-LogR and IS optimize a regularized local cost function ` (·) for each spin i.e., ∀i ∈ V,
Ĵ\i = arg min J\i
[ 1
M M∑ µ=1 ` ( s (µ) i h (µ) \i ) + λ ∥∥J\i∥∥1 ] , (2)
where h(µ)\i ≡ ∑ j 6=i Jijs (µ) j , J\i ≡ (Jij)j(6=i), and ‖·‖1 denotes the `1 norm. Specifically, ` (x) =
log ( 1 + e−2x ) for `1-LogR and ` (x) = e−x for IS, which correspond to the minus log conditional distribution [10] and the ISO [14], respectively. Consequently, the problem of recovering the edge set E is equivalently reduced to local neighborhood selection, i.e., recovering the neighborhood set N (i) for each vertex i ∈ V. In particular, given the estimates Ĵ\i in (2), the neighborhood set of vertex i can be estimated via the nonzero coefficients, i.e.,
N̂ (i) = { j|Ĵij 6= 0, j ∈ V \ i } , ∀i ∈ V. (3)
In this paper, we focus on one simple linear estimator, termed as the `1-LinR estimator, i.e., ∀i ∈ V,
Ĵ\i = arg min J\i
[ 1
2M M∑ µ=1 ( s (µ) i − h (µ) \i )2 + λ ∥∥J\i∥∥1 ] , (4)
which, recalling that s(µ)i ∈ {−1,+1}, corresponds to a quadratic loss ` (x) = 12 (x− 1) 2 in (2). The neighorbood set for each vertex i ∈ V is estimated in the same way as (3). Interestingly, the quadratic loss used in (4) implies that the postulated conditional distribution is Gaussian and thus inconsistent with the true one, which is one typical case of model misspecification. Furthermore, compared with nonlinear estimators `1-LogR and IS, the `1-LinR estimator is more efficient to implement.

heading: 3 Statistical Mechanics Analysis
text: In this section, a statistical mechanics analysis of the `1-LinR estimator is presented for typical RR graphs in the paramagnetic phase. Our analysis is applicable to any M-estimator of the form (2) and please refer to Appendix A for a unified analysis, including detailed results for the `1-LogR estimator.
To characterize the structure learning performance, the precision and recall are considered:
Precision = TP
TP + FP , Recall =
TP
TP + FN , (5)
where TP , FP , FN denote the number of true positives, false positives, and false negatives in the estimated couplings, respectively. The concept of model selection consistent for an estimator is defined in Definition 1, which is also known as the sparsistency property [10]. Definition 1. An estimator is called model selection consistent if both the associated precision and recall satisfy Precision→ 1 and Recall→ 1 as M →∞.
Additionally, if one is further interested in the specific values of the estimated couplings, our analysis can also yield the residual sum of squares (RSS) for the estimated couplings.
Our theoretical analysis of the learning performance builds on the statistical mechanics framework. Contrary to the probably almost correct (PAC) learning theory [37] in mathematical statistics, statistical mechanics tries to describe the typical (as defined in Definition 2) behavior exactly rather than to bound the worst case which is likely to be over-pessimistic [38]. Definition 2. “typical” means not just most probable but in addition the probability for situations different from the typical one can be made arbitrarily small as N →∞ [38].
Similarly, when referring to typical RR graphs, we mean tree-like RR graphs, i.e., when seen from a random node, they look like part of an infinite tree, which are typical realizations from the uniform probability distribution on the ensemble of RR graphs.

heading: 3.1 Problem Formulation
text: For simplicity and without loss of generality, we focus on spin s0. With a slight abuse of notation, we will drop certain subscripts in following descriptions, e.g., J\i will be denoted as J hereafter which represents a vector rather than a matrix. The basic idea of the statistical mechanical approach is to introduce the following Hamiltonian and Boltzmann distribution induced by the loss function ` (·)
H ( J |DM ) = M∑ µ=1 ` ( s (µ) 0 h (µ) ) + λM ‖J‖1 , (6)
P ( J |DM ) = 1
Z e−βH(J|D M), (7)
where Z = ∫ dJe−βH(J|D
M) is the partition function, and β (> 0) is the inverse temperature. In the zero-temperature limit β → +∞, the Boltzmann distribution (7) converges to a point-wise measure on the estimator (2). The macroscopic properties of (7) can be analyzed by assessing the free energy density f(DM ) = − 1Nβ logZ, from which, once obtained, we can evaluate averages of various quantities simply by taking its derivatives w.r.t. external fields [21]. In current case, f(DM ) depends on the predetermined randomness of DM , which plays the role of quenched disorder. As N,M →∞, f(DM ) is expected to show the self averaging property [21]: for typical datasets DM , f(DM ) converges to its average over the random data DM :
f = − 1 Nβ [logZ]DM , (8)
where [·]DM denotes expectation over DM , i.e. [·]DM = ∑ s(1),...,s(M) (·) ∏M µ=1 PIsing ( s(µ)|J∗ ) . Consequently, one can analyze the typical performance of any `1-regularized M-estimator of the form (2) via the assessment of (8), with `1-LinR in (4) being a special case with ` (x) = 12 (x− 1) 2.

heading: 3.2 Replica computation of the free energy density
text: Unfortunately, computing (8) rigorously is difficult. For practically overcoming this difficulty, we resort to the powerful replica method [20, 21, 22, 23] from statistical mechanics, which is symbolized using the following identity
f = − 1 Nβ [logZ]DM = − limn→0 1 Nβ
∂ log [Zn]DM
∂n . (9)
The basic idea is as follows. One replaces the average of logZ by that of the n-th power Zn which is analytically tractable for n ∈ N in the large N limit, and constructs an analytically continuable expression from N to R, then takes the limit n→ 0 by using the expression. Although the replica method is not rigorous, it has been empirically verified from extensive studies in disorder systems [22, 23] and also found useful in the study of high-dimensional models in machine learning [39, 40]. For more details of the replica method, please refer to [20, 21, 22, 23].
Specifically, with the HamiltonianH ( J |DM ) , assuming n ∈ N is a positive integer, the replicated partition function [Zn]DM in (9) can be written as
[Zn]DM = ∫ n∏ a=1 dJae−βλM ∑n a=1‖J a‖1 {∑ s PIsing (s|J∗) exp [ −β n∑ a=1 ` (s0h a) ]}M , (10)
where ha = ∑ j J a j sj will be termed as local field hereafter, and a (and b in the following) is index variable of the replicas. The analysis below essentially depends on the distribution of ha but it is nontrivial. To resolve it, we take a similar approach as [7, 29] and introduce the following ansatz.
Ansatz 1 (A1): Denote Ψ = {j|j ∈ N (0)} and Ψ̄ = {j|j = 1, ..., N − 1, j /∈ N (0)} as the active and inactive sets of spin s0, respectively, then for a typical RR graphG ∈ GN,d,K0 in the paramagnetic phase, i.e., (d− 1) tanh2 (K0) < 1, the `1-LinR estimator in (4) is a random vector determined by random realizations of DM and obeys the following form
Ĵj =
{ J̄j +
1√ N wj , j ∈ Ψ
1√ N wj , j ∈ Ψ̄
(11)
where J̄j is the mean value of the estimator and wj is a random variable which is asymptotically zero mean with variance scaled as O (1). The consistency of Ansatz 1 is checked in Appendix B. Under Ansatz 1, the local fields ha can be decomposed as ha = ∑ j∈Ψ J̄jsj + h a w where h a w ≡ ∑ j 1√ N waj sj is the “noise” part. According to the central limit theorem (CLT), haw can be approximated as multivariate Gaussian variables, which, under the replica symmetric (RS) ansatz [21], can be fully described by two order parameters
Q ≡ 1 N ∑ i,j wai C \0 ij w a j , q ≡ 1 N ∑ i,j wai C \0 ij w b j , (a 6= b), (12)
where C\0 ≡ {C\0ij } is the covariance matrix of the original Ising model without the spin s0. Since the difference between C\0 and that with s0 is not essential in the limit N → ∞, hereafter the superscript \0 will be discarded. As shown in Appendix A, for quadratic loss ` (x) = 12 (1− x)
2 of `1-LinR, the average free energy density (9) in the limit β →∞ can be computed as f (β →∞) = −Extr {−ξ + S} , (13) where Extr {·} denotes the extremum operation w.r.t. relevant variables and ξ, S denote the energy and entropy terms:
S = lim β→∞ lim n→0
1
Nβ
∂
∂n log I, (14)
I = ∫ n∏ a=1 dwa n∏ a=1 e−λβ‖w a‖1δ ∑ i,j wai Cijw a j −NQ ×∏ a<b δ ∑ i,j wai Cijw b j −Nq  , (15)
ξ = αEs,z
( s0 − ∑ j∈Ψ J̄jsj − √ Qz )2
2 (1 + χ) + αλ ∑ j∈Ψ ∣∣J̄j∣∣ , (16)
where α ≡ M/N,χ ≡ limβ→∞ β (Q− q), Es,z(·) denotes the expectation operation w.r.t. z ∼ N (0, 1) and (s0, sΨ) ∼ PIsing(s0, sΨ|J∗) ∝ es0 ∑ j∈Ψ J ∗ j sj [7]. For different losses ` (·), the free energy results (13) only differ in the energy term ξ, which in general is non-analytical (e.g., logistic loss for `1-LogR) but can be solved numerically. Please refer to Appendix A.3 for more details.
In contrast to the case of `2-norm in [29], the `1-norm in (15) breaks the rotational invariance property, i.e., ‖wa‖1 6= ‖Owa‖1 for general orthogonal matrix O, making it difficult to compute the entropy term S. To circumvent this difficulty, we employ an observation that, when considering the RR graph ensemble GN,d,K0 as the coupling network of the Ising model, the orthogonal matrix O diagonalizing the covariance matrix C appears to be distributed from the Haar orthogonal measure [41, 42]. Thus, it is assumed that I in (15) can be replaced by its average [I]O over the Haar-distributed O:
Ansatz 2 (A2): Denote C ≡ Es[ssT ], where Es[·] = ∑ s PIsing(s|J
∗)(·), as the covariance matrix of spin configurations s. Suppose that the eigendecomposition of C is C = OΛOT , where O is the orthogonal matrix, then O can be seen as a random sample generated from the Haar orthogonal measure and thus for typical graph realizations from GN,d,K0 , I in (15) is equal to the average [I]O. The consistency of Ansatz (A2) is numerically checked in Appendix C. Under Ansatz (A2), the entropy term S in (14) can be alternatively computed as lim
n→0 1 Nβ ∂ ∂n log [I]O, as shown in Appendix
A. Finally, under the RS ansatz, the average free energy density (9) in the limit β →∞ reads
f (β →∞) = −Extr Θ  − α2(1+χ)Es,z (( s0 − ∑ j∈Ψ J̄jsj − √ Qz )2) − λα ∑ j∈Ψ ∣∣J̄j∣∣ + (−ER+ Fη)G′ (−Eη) + 12EQ− 1 2Fχ+ 1 2KR− 1 2Hη
−Ezmin w { K 2 w 2 − √ Hzw + λM√ N |w| }
 , (17)
where z ∼ N (0, 1), and G (x) is a function defined as
G (x) = −1 2 log x− 1 2 + Extr Λ
{ −1
2
∫ log (Λ− γ) ρ (γ) dγ + Λ
2 x
} , (18)
and ρ (γ) is the eigenvalue distribution (EVD) of the covariance matrix C, and Θ is a collection of macroscopic parameters Θ = { χ,Q,E,R, F, η,K,H, { J̄j } j∈Ψ } . For details of these macroscopic parameters and ρ (γ), please refer to Appendix A and F, respectively. Note that in (17), apart from the ratio α ≡M/N , N and M also appear as λM/ √ N in the free energy result, which is different from previous results [7, 29, 39]. The reason is that, thanks to the `1-regularization term λM ‖J‖1 , the mean estimates { J̄j } j∈Ψ in the active set Ψ and the noise w in the inactive set Ψ̄ essentially give different scaling contributions to the free energy density.
Although there are no analytic solutions, these macroscopic parameters in (17) can be obtained by numerically solving the corresponding equations of state (EOS) employing the physics terminology. Specifically, for the `1-LinR estimator, the EOS can be obtained from the extremization condition in (17) as follows (for EOS of a general M-esimator and `1-LogR, please refer to Appendix A.3): E = α(1+χ) , F = α (1+χ)2 [ Es ( s0 − ∑ j∈Ψ sj J̄j )2 +Q ] , R = 1K2 [ ( H + λ 2M2 N ) erfc ( λM√ 2HN ) − 2λM √ H N 1√ 2π e− λ2M2 2HN ], Eη = − ∫ ρ(γ) Λ̃−γ dγ, Q = FE2 +RΛ̃− (−ER+Fη)η∫ ρ(γ) (Λ̃−γ)2 dγ , K = EΛ̃+ 1η , χ = 1E + ηΛ̃, H = Rη2 + FΛ̃+ (−ER+Fη)E∫ ρ(γ) (Λ̃−γ)2 dγ , η = 1K erfc ( λM√ 2HN ) ,
J̄j = soft(tanh(K0),λ(1+χ)) 1+(d−1) tanh2(K0) , j ∈ Ψ,
(19)
where Λ̃ satisfying Eη = −
∫ ρ(γ)
Λ̃−γ dγ is determined by the extremization condition in (18) and soft (z, τ) = sign (z) (|z| − τ)+ is the soft-thresholding function. Once the EOS is solved, the free energy density defined in (8) is readily obtained.

heading: 3.3 High-dimensional asymptotic result
text: One important result of our replica analysis is that, as derived (see Appendix A.3) from the free energy result (17), the original high dimensional `1-LinR estimator (4) is decoupled into a pair of scalar estimators, one for the active set and one for the inactive set, i.e.,
Ĵj =  soft (tanh (K0) , λ (1 + χ)) 1 + (d− 1) tanh2 (K0) ≡ J̄j , j ∈ Ψ (20) √ H
K √ N soft
( zj ,
λM√ HN
) , j ∈ Ψ̄ (21)
where zj ∼ N (0, 1) , j ∈ Ψ̄ are i.i.d. standard Gaussian random variables. The decoupling property asserts that, once the EOS (19) is solved, the asymptotic behavior of `1-LinR can be statistically described by a pair of simple scalar soft-thresholding estimators (see Figs. 1(a) and 1(b)).
In the high-dimensional setting where N is allowed to grow as a function of M , one important question is that what is the minimum number of samples M required to achieve model selection consistency as N → ∞. Though we obtain a pair of scalar estimators (20) and (21), there are no analytical solutions to EOS (19), making it difficult to derive an explicit condition. To overcome this difficulty, as shown in Appendix D, we perform a perturbation analysis of EOS (17) and obtain an asymptotic relation H ' F , Then, we obtain that for a RR graph G ∈ GN,d,K0 , given M i.i.d. samples DM , the `1-LinR estimator (4) can consistently recover the graph structure G as N →∞ if
M > c (λ,K0) logN
λ2 , 0 < λ < tanh (K0) , (22)
where c(λ,K0) is a constant value dependent on the regularization parameter λ and coupling strength K0 and a sharp prediction (as verified in Sec. 5) is obtained as
c (λ,K0) = 2 ( 1− tanh2 (K0) + dλ2 ) 1 + (d− 1) tanh2 (K0) . (23)
For details of the analysis, including the counterpart of `1-LogR, see Appendix D. Consequently, we obtain the typical sample complexity of `1-LinR for Ising model selection for typical RR graphs in the paramagnetic phase. The result in (22) is derived for `1-LinR with a fixed regularization parameter λ. Since the value of λ is upper bounded by tanh (K0) (otherwise false negatives occur as discussed in Appendix D), a lower bound of the typical sample complexity for `1-LinR is obtained as
M > 2 logN
tanh2 (K0) . (24)
Interestingly, the scaling in (24) is the same as the information-theoretic worst-case resultM > c logN K20 obtained in [11] at high temperatures (i.e., small K0) since tanh (K0) = O (K0) as K0 → 0.
3.4 Non-asymptotic result for moderate M,N
In practice, it is desirable to predict the non-asymptotic performance of the `1-LinR estimator for moderate M,N . However, the scalar estimator (20) for the active set (see Fig. 1(a)) fails to capture the fluctuations around the mean estimates. This is because in obtaining the energy term ξ (16) of the free energy density (17), the fluctuations around the mean estimates { J̄j } j∈Ψ are averaged out by
the expectation Es,z (·). To address this problem, we replace Es,z (·) in (17) with a sample average by accounting for the finite-size effect, thus obtaining a modified estimator for the active set as follows
{Ĵj} j∈Ψ = arg min Jj,j∈Ψ
 ∑M µ=1 ( sµ0 − ∑ j∈Ψ s µ j Jj − √ Qzµ )2 2 (1 + χ)M + λ ∑ j∈Ψ |Jj |  , (25) where sµ0 , s µ j,j∈Ψ ∼ P (s0, sΨ|J
∗) , zµ ∼ N (0, 1) , µ = 1...M . The modified d-dimensional estimator (25) (see Fig. 1(c) for a schematic) is equivalent to the scalar one (20) (Fig. 1(a)) asM →∞ but it enables us to capture the fluctuations of {Ĵj}j∈Ψ for moderate M . Note that due to the replacement of expectation with sample average in the free energy density (17), the EOS (19) also needs to be modified and it can be solved iteratively as sketched in Algorithm 1. The details are shown in Appendix E.1.
Algorithm 1: Method to solve EOS (19) together with (25) Input: M,N, λ,K0, ρ (γ) and TMC Output: χ,Q,E,R, F, η,K,H, {Ĵtj} j∈Ψ Initialization: χ,Q,E,R, F, η,K,H
1 MC sampling: For t = 1...TMC , draw random samples sµ,t0 , { sµ,tj } j∈Ψ ∼ P (s0, sΨ|J ∗) and
zµ,t ∼ N (0, 1), µ = 1...M 2 repeat 3 for t = 1 to TMC do
4 Solve {Ĵtj} j∈Ψ = arg min Jj,j∈Ψ
[∑M µ=1(s µ,t 0 − ∑ j∈Ψ s µ,t j Jj− √ Qzµ,t)2
2(1+χ)M + λ ∑ j∈Ψ |Jj | ] 5 Compute4t = 1
M ∑M µ=1 ( sµ,t0 − ∑ j∈Ψ s µ,t j Ĵ t j )2 6 Solve the EOS (19) with Es ( si − ∑ j∈Ψ sj J̄j )2 = 1
TMC
∑TMC t=1 4 t
7 Update values of χ,Q,E,R, F, η,K,H 8 until convergence
Consequently, for moderate M,N , the non-asymptotic statistical properties of the `1-LinR estimator can be characterized by the reduced d-dimensional `1-LinR estimator (25) (Fig. 1(c)) and scalar estimator (21) (Fig. 1(b)) using MC simulations. Denote {Ĵ tj}, t = 1, ..., TMC as the estimates in t-th MC simulation, where {Ĵ tj}j∈Ψ and {Ĵ tj}j∈Ψ̄ are solutions of (25) and (21), and TMC is the total number of MC simulations. Then, the Precision and Recall are computed as
Precision = 1
TMC TMC∑ t=1 ∥∥∥Ĵ tj,j∈Ψ∥∥∥ 0∥∥∥Ĵ tj,j∈Ψ∥∥∥ 0 + ∥∥∥Ĵ tj,j∈Ψ̄∥∥∥0 , Recall = 1 TMC TMC∑ t=1 ∥∥∥Ĵ tj,j∈Ψ∥∥∥ 0 d , (26)
where ‖·‖0 is the `0-norm indicating the number of nonzero elements. In addition, the RSS can be computed as RSS = ∑ j ∣∣∣Ĵj − J∗j ∣∣∣2 = 1TMC ∑TMCt=1 ∑j∈Ψ ∣∣∣Ĵ tj −K0∣∣∣2 +R.

heading: 4 Discussions
text: It might seem surprising that, despite apparent model misspecification due to the use of quadratic loss, the `1-LinR estimator can still correctly infer the structure with the same order of sample complexity as `1-LogR. Also, our theoretical analysis implies that the idea of using linear regression for binary data is not as outrageous as one might imagine. Here we provide an intuitive explanation of its success with a discussion of its limitations.
On the average, from (4), the condition for the `1-LinR estimator is given as 〈s0sk〉 − ∑ j 6=0 〈sjsk〉Jj = λ∂|Jk|, k = 1, . . . , N, (27)
where 〈·〉 and ∂|Jk| represent average w.r.t. the Boltzmann distribution (7) and the sub-gradient of |Jk|, respectively. In the paramagnetic phase, 〈sisj〉 decays in its magnitude exponentially w.r.t. the
distance of sites i and j. This guarantees that once the connections Jk of sites in the first nearest neighbor set Ψ are given so that
〈s0sk〉 = ∑ j∈Ψ 〈sjsk〉Jj + λsign(Jk), ∀k ∈ Ψ (28)
holds, the other conditions are automatically satisfied by setting all the other connections that are not from Ψ to zero. For appropriate choice of λ, (28) has solutions of sign(J∗k )Jk > 0,∀k ∈ Ψ. Namely ∀k ∈ Ψ, the estimate of Jk has the same sign as the true value J∗k . This implies that on average the `1-LinR estimator can successfully recover the network structure up to the connection signs if λ is chosen appropriately.
The key of the above argument is that 〈sisj〉 decays exponentially fast w.r.t. the distance of two sites, which does not hold after the phase transition. Thus, it is conjectured that the `1-LinR estimator will start to fail in the network recovery just at the phase transition point. However, it is worth noting that this is in fact not limited to `1-LinR: `1-LogR also exhibits similar behavior unless post-thresholding is used, as reported in [43].

heading: 5 Experimental results
text: In this section, we conduct numerical experiments to verify the accuracy of the theoretical analysis. The experimental procedures are as follows. First, a random graph G ∈ GN,d,K0 is generated and the Ising model is defined on it. Then, the spin snapshots are obtained using the Metropolis–Hastings algorithm [44, 45, 46] in the same way as [7], yielding the dataset DM . We randomly choose a center spin s0 and infer its neighborhood using the `1-LinR (4) and `1-LogR [10] estimators. To obtain standard error bars, we repeat the sequence of operations 1000 times. The RR graph G ∈ GN,d,K0 with node degree d = 3 and coupling strength K0 = 0.4 is considered, which satisfies the paramagnetic condition (d− 1) tanh2 (K0) < 1. The active couplings {Jij}(i,j)∈E have the same probability of taking both signs of +1 or −1 2. We first verify the precise non-asymptotic predictions of our method described in Sec.3.4. Fig. 2 (upper figure) shows the replica and experimental results of RSS, Precision, Recall for N = 200 with different values of α ≡M/N . It can be seen that for both `1-LinR and `1-LogR, there is a fairly
2Though this setting is different from the analysis where the nonzero couplings take a uniform sign, the result can be directly compared thanks to gauge symmetry [21].
good agreement between the theoretical predictions and experimental results, even for smallN = 200 and small α (equivalently small M ), verifying the correctness of the replica analysis. Interestingly, a quantitatively similar behavior between `1-LinR and `1-LogR is observed in terms of precision and recall. Regarding RSS, the two estimators actually behave differently, which can be clearly seen in
Fig. 7 in Appendix G: the RSS is much smaller for `1-LogR, which is natural since the estimates of `1-LogR are closer to the true ones due to the model match. As our theoretical analysis is based on the typical tree-like structure assumption, it is interesting to see if it is applicable to graphs with loops. To this end, we consider the 2D 4-nearest neighbor grid with periodic boundary condition, which is one common loopy graph. Fig. 2 (lower figure) shows the results for a 15× 15 2D grid with uniform constant coupling K0 = 0.2. The agreement between the theoretical and numerical results is fairly good, indicating that our theoretical result can be a good approximation even for loopy graphs. More results for different values of N and λ are shown in Fig. 7 and Fig. 8 in Appendix G.
Subsequently, the asymptotic result and sharpness of the critical scaling value c0 (λ,K0) ≡ c(λ,K0)λ2 in (22) are evaluated. First, Fig. 3 (left) shows comparison of c0 (λ,K0) between `1-LinR and `1-LogR for the RR graph G ∈ GN,d,K0 when d = 3,K0 = 0.4, indicating similar behavior of `1-LogR and `1-LinR. Then, we conduct experiments for M = c logN with different values of c around c0 (λ,K0), and investigate the trend of Precision andRecall asN increases. When λ = 0.3, Fig. 3 (middle and right) show the results of Precision and Recall, respectively. As expected, the Precision increases consistently with N when c > c0 (λ,K0) and decreases consistently with N when c < c0 (λ,K0) while the Recall increases consistently and approaches to 1 as N →∞, which verifies the sharpness of the critical scaling value prediction. The results for `1-LogR, including the case of λ = 0.1 for both `1-LinR and `1-LogR, are shown in Fig. 9 and Fig. 10 in Appendix G.

heading: 6 Conclusion
text: In this paper, we provide a unified statistical mechanics framework for the analysis of typical learning performances of `1-regularized M -estimators, `1-LinR in particular, for Ising model selection on typical paramagnetic RR graphs. Using the powerful replica method, the high-dimensional `1- regularized M-estimator is decoupled into a pair of scalar estimators, by which we obtain an accurate estimate of the typical sample complexity. It is revealed that, perhaps surprisingly, the misspecified `1-LinR estimator is model selection consistent using M = O (logN) samples, which is of the same order as `1-LogR. Moreover, with a slight modification of the scalar estimator for the active set to account for the finite-size effect, we further obtain sharp predictions of the non-asymptotic behavior of `1-LinR (also `1-LogR) for moderate M,N . There is an excellent agreement between theoretical predictions and experimental results, even for graphs with many loops, which supports our findings. Several key assumptions are made in our theoretical analysis, such as the paramagnetic assumption which implies that the coupling strength should not be too large. It is worth noting that the restrictive paramagnetic assumption is not only limited to `1-LinR, but also to other low-complexity estimators like `1-LogR unless post-thresholding is used [43]. These assumptions restrict the applicability of the presented result, and thus overcoming such limitations will be an important direction for future work.

heading: Acknowledgements
text: This work was supported by JSPS KAKENHI Nos. 17H00764, 18K11463, and 19H01812, and JST CREST Grant Number JPMJCR1912, Japan.

Review summaries:
The manuscript contains a theoretical investigation of the Ising model selection problem.  The task is to recover the support (the graph topology) of the interaction network from observations of Ising configurations.   An objective function in the M-estimator family, named L1-LinR by the authors, is analyzed.  It contains the typical quadratic data-dependent loss from linear regression and an L1 regularization term,  therefore it is a mismatched objective compared to L1-LogR (pseudolikelihood) and the Interaction Screening. 

The submission considers the problem of recovering the underlying graph structure of an Ising model defined on N nodes, given M i.i.d. samples. More specifically, the authors focus on Ising models associated with typical random regular graphs in the paramagnetic phase. Ideas from statistical physics are borrowed to study the performance of an estimator derived from an $\ell_1$-regularized linear regression problem. It is claimed that using $M = \mathcal{O}(\log(N))$ samples is enough to obtain model selection consistency. Non-asymptotic results are also given, and some numerical experiments are presented to support the findings.

The authors study the Ising model selection problem and compute analytically via the replica method the typical performances of the $\ell_1$-regularized linear regression. They confirm the theoretical result by running a good number of numerical simulations. They demonstrate the sample complexity is such that O(log N) samples are enough to correctly identify the right neighborhood of a generic variable. The theory presented is able to predict the precision and the recall rates also for finite values of the system size and a small number of samples. Moreover, the theory seems to work fine also for graphs with many loops. The methodology of the present work can be extended also to other estimators.

This paper studies the problem of estimating the edge couplings in a Markov random field using ell-one regularized least-squares estimation. It is assumed throughout that all nodes have the same degree and the nonzero coupling have the same value. Using the replica method, the authors derive a set of couple equations of state (EOS), which depend on the underlying parameters of the model (e.g., the degree and coupling strength),  the ratio between the number of samples and the number of nodes, and the spectral density of the covariance of the resulting MRF. By analyzing the solutions to the equations, the authors provide numerical predictions of the reconstruction performance (e.g, precision and recall) in the limit of high dimension. Further analysis of these expressions gives a sharp predicition for the number of samples needed for exact recovery as a function of the the degree and the coupling strength. Additionally, the authors suggest a modification of the EOS that uses the observed data to produce more accurate predictions. Because the analysis depends on the replica method the results are not rigorous. The accuracy of the theoretical analysis is supported by numerical simulations.