Ising Model Selection Using $\ell_{1}$-Regularized Linear Regression: A Statistical Mechanics Analysis
Authors: Xiangming Meng
Tomoyuki Obuchi
Yoshiyuki Kabashima
Venue: NeurIPS 2021 Poster (2021)
Keywords: Ising model
Lasso
replica method
statistical physics

Abstract: We theoretically analyze the typical learning performance of $\ell_{1}$-regularized linear regression ($\ell_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 $\ell_1$-LinR is obtained.   Remarkably, despite the model misspecification, $\ell_1$-LinR is model selection consistent with the same order of sample complexity as $\ell_{1}$-regularized logistic regression ($\ell_1$-LogR), i.e., $M=\mathcal{O}\left(\log N\right)$,  where $N$ is the number of variables of the Ising model. Moreover, we provide an efficient method to accurately predict the non-asymptotic behavior of $\ell_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 $\ell_1$-LinR, our method is readily applicable for precisely characterizing the typical learning performances of a wide class of  $\ell_{1}$-regularized $M$-estimators including $\ell_1$-LogR and interaction screening. 

Ising Model Selection Using l--Reuularized Linear Regression: A Statistical Mechanics Analysis
Xiangming Meng"
Tomoyuki Obuchi
Institute for Physics of Intelligence
Department of Systems Science
The University of Tokyo
Kyoto University Kyoto 606-8501, Japan obuchi@i wootoiiikootuuuaa ac.
7-3-1, Hongo, Tokyo 113-0033, Japan meng0g. ecc. w.c..u.tooo..at ac. jp
Yoshiyuki Kabashima
Institute for Physies of Intelligence
The University of Tokyo
7-3-1, Hongo, Tokyo 113-0033, Japan kabaüphys. S.ei:ppyy....uooky... w-wtokyo.acc..... ac.
Abstract
We theoretically analyze the typical learning performance e-rregularized linear regression (& (LinR) for Ising model selection using the replica method from statis- tical mechanics. For typical random regular graphs in the paramagnetic phase, an accurate estimate of the typical sample complexity of LLiiR is obtained. Remark- ably, despite the model misspecification, Li-LinR is model selection consistent with the same order of sample complexity & reguularized logistic regression (L1-LogR), i.e., M = (log where is the number variables the Ising model. Moreover, provide efficient method to accurately predict the asymptotic behavior of Li-LinR for moderate M, N, such as precision and recall. Simulations show fairly good agreement between theoretical predictions and experimental results, even for graphs with many loops, which supports our findings. Although this paper mainly focuses e LLLiR,, our method is readily applicable for precisely characterizing the typical learning performances of a wide class of rereeuularzzed M-estimators including fj-LogR and interaction screening.
Introduction
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 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
6 (7). Among various applications, one fundamental problem of interest is called Ising model selection, which refers to recovering the underlying graph structure the original Ising model from independent, identically distributed (i.i.d.) samples. variety of methods have been [5] 16, 17 [8]. demonstrating the possibility successful Ising model selection even when the number of samples is smaller than that of the variables. Notably, has been demonstrated that for the & r-reeullaiized logistic regression (& (oggR) [10] [6] and interaction screening (IS) [14 [ estimators, M O (log N) samples suffice for Ising model with N spins under certain assumptions, which is consistent with respect to (w.r.t)) previously established
(5,
proposed (8] 10 12 [3] 14]
"Coreesoodiig author.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

information -hhooretic lower- [ Both & Lo-LogR and IS are <1-regularized -estimattors 19 with logistic and IS objective (ISO) loss functions, respectively.
In this paper, we focus on one simpler linear estimator called & f1-regularized linear regression (61- Link) and theoretically investigate its rypical learning performance using the powerful replica method [20 22 [3] from statistical mechanics. The E Li-LinR estimator, widely known as least absolute shrinkage and selection operator (LASSO) [24] in statistics and machine learning, considered here mainly for two reasons. On the one hand, it is one representative example of model misspecification since the quadratic loss of f LinLi does not match the true log- segg--odttiooaa-ikkeeee likelliollooo as LooLo, nor does it have the interaction screening property as IS. On the other hand, as one of the most popular linear estimator, 4 LinR is more computationally efficient than Logg and IS, and thus is of practical importance to investigate its learning performance for Ising model selection. Since it difficult to obtain results for general graphs, as a first step we consider the random regular (RR) graphs Gv,d,K, Ko in the paramagnetic phase [23], where GN.d.K, denotes the ensemble of RR graphs with constant node degree d and uniform coupling strength Ko on the edges.
1.1 Contributions
The main contributions are summarized as follows, First, we obtain an accurate estimate of the typical sample complexity of RLnin for Ising model selection for typical RR graphs in the paramagnetic phase, which, remarkably. has the same order as Lo-Log Specifically, for typical RR graph G e GN.dKK. using l Linn with regularization parameter 0 < À < tanh (Ko), one can consistently (aali (do)) I+(d-1) tanb"
reconstruct the structure with M e(). log N samples, where e (A, Kn) 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 rypical sample complexity for Ising model selection. Interestingly, as À - tanh (Ko), a lower bound of the typical sample complexity is obtained, which has the same scaling as the M tanh [Ka]
X²
2log
information-theoreti lower bound M > d log
N for some constant c' at high temperatures (i.e.,
Kg
small Ko) since tanh (Ko)= O (Ko) as Ko O.
Second, we provide a computationally efficient method to precisely predict the typical learning performance of RLiRR 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 LLnR for Ising model selection have been unavailable even for LogL [10] [6] and IS 14 15). nor are they the
same as previous asymptotic results of -LinR assuming fixed a = 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 LinLi,,, our method is readily applicable to wide class of Errreuularized -esstmaators [19, including LogL (10) and IS [4[15). Thus, an additional technical contribution is providing generic approach for precisely characterizing the typical learning performances of various l rerggllarized M Meestimators 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 the replica method tool for performance analysis.
1.2 Related works
There has been some earlier works the analysis of Ising model selection (also known as the inverse Ising problem) using the replica method 4
from statistical mechanics. For example, in (PL) method [30 is studied. However, instead of graph [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 fzrregularized linear regression but the techniques invented there are not applicable to -LinR since the E -norm breaks the rotational invariance property.
[6]. the performance of the pseudo-likelihen
structure learning,
Regarding the study of LinR (LASSSO) 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
2

---

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 L1-LinR while most previous studies consider estimation consistency except
3 However, [ 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 Li-Ln estimator for Ising model selection, though some are found for Gaussian graphical models [35] 36.. One closely related work [TS] 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 rypical sample complexity for consistently reconstructing typical samples of Ising models defined over the RR graphs. Furthermore, does not provide precise predictions of the non-asymptotic learning performance we do.
[15]
2
Background and Problem Setup
2.1 Ising Model
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 has the form
of Ising model with N variables (spins) s (s:) N-1 E +1,
1
Prsing (s|J")
Zising (J") exp
where Zising (J)) - -. exp ii Jijs18} } is the partition function and J = -J on 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 undirected graph G = (V, E), where V = {0, 1,., - }} is a collection of vertices at which the spins are assigned, and E (i [J # 0} a collection of undirected edges, i.e., Jú = 0 for all pairs of (i,j) & For each vertex E V, its neighborhood is defined as the subset N (i) = {j E v (i.j) E).
2.2 Neighborhood-base é reggularreed linear regression (e, -LinR)
The problem of Ising model selection refers to recovering the graph G' (edge set E), given M i.i.d. samples DMM = { ..., 8 { 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 l LogR estimator [10] and IS estimator [14]. Both LogL and IS optimize regularized local cost function é (-) for each spin i.e., Vi E V,
i arg min
Σ
(1()) (()) + X +11J
(2)
where n(4)) = Ejji 5M) J. = (J))) and I|- denotes the E, norm. Specifically, E (xx
log (1 11+e-20 -2x) for E LLogR and ll (xx =- for IS, which correspond to the minus log conditional distribution I0) 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 E V. In partieular, given the estimates J in (2), the neighborhood set of vertex i can be estimated via the nonzero coefficients, i.e.,
(i) J + O.JEV ViEV.
(3)
In this paper, we focus on one simple linear estimator, termed as the e -LinR estimator, i.e., Vi E V,
J = arg min
Σ
1 - m(/1)² + + +AllJool
(4)
2M
J x
==1
3

---

which, recalling that w/!)) E (11,+1) +1), corresponds quadratic loss 8 (x) (z 1)) in 2 The neighorbood set for each vertex E 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 & LogR and IS, the l LinR estimator is more efficient to implement.
3 Statistical Mechanics Analysis
In this section, a statistical mechanics analysis of the LinLi 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 unified analysis, including detailed results for the & -Log estimator.
To characterize the structure learning performance, the precision and recall are considered:
TP
TP
Precision TP + FP" Recall T
(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 estimator is defined in Definition[] 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 + and Recall
as M + DC.
Additionally, if one is further interested 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 the learning performance builds on the statistical mechanics framework. Contrary to the probably almost correct (PAC) learning theory (37 in mathematical statistics, statisti- cal mechanies tries to describe the rypical (as defined Definition behavior exactly rather than to bound the worst case which is likely to be over wer-pessimistic (38..
Definition 2. "rypical" means just most probable but in addition the probability for situations different the rypical one can be made arbitrarily small as N oo [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 infinite tree, which are rypical realizations from the uniform probability distribution on the ensemble of RR graphs.
3.1 Problem Formulation
For simplicity and without loss of generality, we focus spin sp. With slight abuse of notation, will drop certain subscripts in following descriptions, e.g., J will be denoted as hereafter which represents vector rather than matrix. The basic idea of the statistical mechanical approach is to introduce the following Hamiltonian and Boltzmann distribution induced by the loss function E (-)
M
H. (J = Σ
(00)))
hw)) + AM ||
(6)
1=1 (HIDDM z
P ((J(DM)
(7)
where 3 f (H( is the partition function, and 3 (> 0) is the inverse temperature. In the zero-temperature limit 3 to0, the Boltzmann distribution ? converges
wisn
measure on the estimator (2). The macroscopic propertics of (7) can be analyzed by assessing the
free energy density f(DMM )
No log Z, 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, ((DM) depends the predetermined randomness of DM which plays the role of quenched disorder. As
N. M + oo, J(DM) is expected to show the self averaging property (1] for typical datasets D f(DM) converges to its average over the random data DM:
1
[log Zl
NB
where p denotes expectation over DM, i.e. HP = Extt An OI - Pking (s(m) Consequently, one can analyze the typical performance of any fi-regularized M-estimator of the form via the assessment of (8) with LiiLi in (4) being a special case with E (x) = 1 (x 1)².

---

3.2 Replica computation of the free energy density
Unfortunately, computing (8) rigorously is difficult. For practically overcoming this difficulty, we resort to the powerful replica method 20 using the following identity
21
22
23 from statistical mechanics, which is symbolized
1 alog [2"]
llog 2I = = lim
(9)
N3
N3 Ôn
The basic idea is as follows. One replaces the average of log Z by that of the n-th power 2" which analytically tractable for E N in the large N limit, and constructs an analytically continuable expression from N to IR, then takes the limit n 0 by using the expression. Although the replica method is not rigorous, has been empirically verified from extensive studies in disorder systems [22 [33 and also found useful in the study of high sigh-dimensional models in machine learning [39] 40)
For more details of the replica method, please refer to 120 21
23
Specifically, with the Hamiltonian H. (J(DM), assuming n E N is a positive integer, the replicated partition function [Z"]
DAr in 9 can be written as
n 1 1 / II dJ BAM PAMES-LIJ
Σ Eising (s|J") exp 3 > e (
1} (10)
where Ata x==>x E, JJ will be termed as local field hereafter, and a (and in the following) is index variable of the replicas. The analysis below essentially depends on the distribution of ht but is and introduce the following ansatz.
nontrivial. To resolve it, we take similar approach as 7 [9]
- 1. & N (0)) as the active and inactive sets spin So, respectively, then for typical RR graph G E &N.d,Ko Ko in the paramagnetic phase, i.e., (d - 1) tanh? 2 (Ko) < 1, the & Li-LinR estimator in (4) is random vector determined by random realizations of D and obeys the following form
Ansatz (AI): Denote = -ilj E EN (0)) and v = -نj==,
J + VN
Wys
jE V
J
(11)
VN ")
jE
where Jj is the mean value of the estimator and wy is a random variable which is asymptotically zero mean with variance scaled as O (1).
The consistency of Ansatz. is checked in Appendix B Under Ansatz 1, the local fields ht" can be decomposed ha Ejew JJs; + ha where ha the central limit theorem (CLT), hü can be approximated as multivariate Gaussian variables, which, under the replica symmetric (RS) ansatz [21], can be fully described by two order parameters
E3 vN 53 is the "noise" part. According to
1 Q = N ==NZ
Σ w wg: q=
w C. ug, (a b),
(12)
N
where C (C)]] is the covariance matrix of the original Ising model without the spin so. Since the difference between C and that with so is not essential in the limit N - oc, hereafter the superscript 10 will be discarded. As shown in Appendix A for quadratic loss l (z) ]1 x)² of
LinR, the average free energy density (9) in the limit 3 - +
can be computed as
[B + oo) = Extr (-E + S}
(13)
where Extr {.} denotes the extremum operation w.r.t. relevant variables and &, S denote the energy and entropy terms:
S = lim lim 8- N3 On
log 1,
(14)
dw a Il
e Aal
Σ in
w Cifuu - NQ
x Il 8 axch
Σ f
w
Ngg
a=1 1=1
(15)
alEg,s 80 - -jeE J19j Jys)=VQx)²²
vQz)²
g {=
+ +ad Σ |J,]
(16)
2 + x)
jEW
5

---

where = MM/Nx = lo B(Q 9). (.) denotes the expectation operation w.r.t. 2 ~ NO 1) and (5a, s4) ( sq [") & e 2 di [7]. For different losses é (-), the free energy results (13) only differ in the energy term E, which in general is non-analytical (e.g., logistic loss for LLoogR) but can be solved numerically. Please refer to Appendix |A.3 for more details.
In contrast to the case of no-m in [29]. the norm in (1)) breaks the rotational invariance property, Le., ||0 II * ||Ow# I|, 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. .d Ko, 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 / in (1)) can be replaced by its average (To over the Haar-distributed O:
Ansatz 2 (A2): Denote C' CE= Es(ss"), where Ex[:] = E. Pising (s)"nn as the covariance matrix of spin configurations 8. Suppose that the eigendecomposition of C is C = ONOT, 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 &d 1 in 15) is equal to the average (4o.-
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 1 0 log [ o. as shown in Appendix n-0 ND an
Finally, under the RS ansatz, the average free energy density (9) in the limit 3 - co reads
x1+x) B.s (000 -sev - ) da JE |
J,s;
(B + 00) = Extrr
+ (C-ER- + FT) G' (-En) + EE _ + 1 1KR 1
min K w= v VHaw + AM
VN |oa|
(17)
where 2 ~~ (0,1), and G (z) is a function defined as
1 2 log + Extr 1 2 f log (A - (A-7) dy +
{
A 2
G (x)
(18)
A
and p (7) is the eigenvalue distribution (EVD) of the covariance matrix C, and O is collection of E, R, parameters and p (1). please refer to Appendix A and F respectively. Note that in 17 7. apart from the ratio a = M/N, N and M also appear as XM/ /M/VVV the free energy result, which different
macroscopic parameters θ = fx. E,
F, n,, K, H, (J,} see) } For details of these macroscopic
from previous results 7 [2 (39. The reason is that, thanks to the E f-rreuularization term XM ||]|| the mean estimates } jew in the active set V and the noise w in the inactive set V essentially give different scaling contributions to the free energy density.
Although there are no analytic solutions, these macroscopic parameters in (7) can be obtained by numerically solving the corresponding equations of state (EOS) employing the physics terminology. Specifically, for the RiiRi estimator, the EOS can be obtained from the extremization condition in as follows (for EOS of a general M-esimator and l LooL,,, please refer to Appendix A.3).
E (1+x)*
Eg 50 --Ejje oJ)) +
F' (1+x)²
a
R = 1
H +
X²112) N
erfe
1 c - 1, V2HN) ) 22M V N /2m
AM
Ka
EΙj-
pl9)
d7;
A-3
Q EN
F +
RÁ-
(-ER+Fmn
[a-1)2 dy
(19)
K EA +
1 + T
X = E
R H T
+ FF+ +
(- EERRFAE dy
I erffe
AM
" = K
V2HN
J, I+(d-1) 1) tanh-
soft(tanh(Ka, Ka), A(1+x))
jj E V

---

soft (2.
soft 2.
Reduced (.Link
2f (000)
scalar estimator active
(+) scalar eximutor
reduced didimensional F. Link
Figure 1: Equivalent low-dimensional estimators for high-dimensional LLin obtained from the statistical mechanics analysis. (a) and (b) are diagrams the pair of scalar estimators Eqs. (20) and (2T). (c) is
schematic description the modified estimator in Eq. (25 which takes into account the finite-size effect.
where Á satisfying En J p(r) dr is determined by the extremization condition in (18 and - A-7 soft (z, r) = sign (2) (|2) - + is the soft setethhreshlldin function. Once the EOS is solved, the free energy density defined in (8) is readily obtained.
3.3 High-dimensional asymptotic result
One important result of our replica analysis that, as derived (see Appendix A.3) from the free energy result (17). the original high dimensional E RinLinR estimator (4) is decoupled into a pair of scalar estimators, one for the active set and one for the inactive set, i.e.,
soft (tanh (Ko) A x)) = Jj, jEV
(20)
1+ (d
1) tanh (Ko)
J,= =
VH
AM
soft
Zj11
jeW v
(21)
KVN
VHN
where 2j ~ (0.1). v are i.i.d standard Gaussian random variables. The decoupling property asserts that, once the EOS (19) solved, the asymptotic behavior of L can be statistically described by pair of simple scalar soft-thresholding estimators (see Figs. a and (b)).
In the high-dimensional setting where N allowed to grow a functioon of M, one important question is that what is the minimum number of samples M required to achieve model selection consistency as N - Dc. Though we obtain 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 E GN.d, given M i.i.d. samples D the & -LnR estimator (4) can consistently recover the graph structure G as N + co if
M > c cA,, Ko) log
0 X < tanh (Ku):
(22)
X²
where c(A, Ko) is constant value dependent on the regularization parameter and coupling strength Ko and sharp prediction (as verified in Sec. 5) is obtained as
c (A) = 2 (1 tanh? 2 (Ko) + dd2)
(23)
+ (d - 1) tanli (Ko)
For details of the analysis, including the counterpart of L LogR, see Appendix D Consequently, we obtain the typical sample complexity of -LinR for Ising model selection for typical RR graphs the paramagnetic phase. The result in (22) is derived for Li-LinR with a fixed regularization parameter X. Since the value of À is upper bounded by tanh (Kg) (otherwise false negatives occur as discussed in Appendix D). lower bound of the typical sample complexity for & LinR is obtained as
2 log N tanh? (Ko)
M
(24)
Interestingly, the scaling in 24 is the same as the information-ther theoretic worst-case result M > clog N
obtained []] at high temperatures (i.e., small Ko) since tanh (Ko) = O (Ko) as Ko + 0.
3.4 Non-asymptotic result for moderate M, N
In practice, it is desirable to predict the non-asee con-asymptotic performance of the & -LiRR estimator for moderate M, .N. However, the scalar estimator (20) for the active set (see Fig. [(a) fails to capture the fluctuations around the mean estimates. This because in obtaining the energy term (16) of the free energy density (17), the fluctuations around the mean estimates {J,} jev are averaged out by
7

---

the expectation E,,, (+). To address this problem, we replace 5,, (-) in (1 with sample average by accounting for the finite-size effect, thus obtaining a modified estimator for the active set follows
Su=i (s6 -je of vQ2*
(J)) jEv arg min J,
[3 JEW
(25)
2(1 + x) M
where og 8 J.jev P (S0, 8w J*), z# = (0,1) ;M L The modified dddimensional esti- I...M mator (25) (see Fig. ()) for schematic) is equivalent to the scalar one (20) (Fig. a)) as M + oo but it enables us to capture the fluctuations of {J,} hew 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 modi-
fied and it can be solved iteratively as sketched in Algorithm The details are shown in Appendix E.
Algorithm 1: Method to solve EOS (19) together with (25)
Input: M,, N, A, Ko, p (7) and Tuc
Output: x. Q, E, R. F, K, H, (Jf)
Jev
Initialization: X. Q, E, R. F, T), K,
MC sampling: For = Mc, draw random samples sb* [s sew ~ P (600 (ulssll))) (000,801J1)
and
g/t N (0, 1), 1...M
repeat
for = to TMc do
- Dise ( -Eeo a -,
Solve (J jEN = arg min
E ESE |را|
2(1+xiM
Compute A' M Σ Ey=1
s E Ejen s4t 13))
Solve the EOS 19 with B. (s. 2a-Da -iee SyJ,)
LByc A
TMe
Update values
X. Q, F., R, P. K,H
until convergence
Consequently, for moderate M, N, the non-asymptotic statistical properties of the LinLi estimator can be characterized by the reduced d-dimensional LinR estimator (25) (Fig. [c) and scalar
estimator (21) (Fig. (b)) using MC simulations. Denote {J}}.t t 1, TMC as the estimates in
1-th MC simulation, where (Ji) and {J)) are solutions of (25) and 21 and TMC the total number of MC simulations. Then, the Precision and Recall are computed as
jt
jt JEE
1 TMC
1 Tse
Precision
Recall
(26)
TMC
Jijew
jt
+
jt j
TMC
d
lo
where Illllo the fo-norm indicating the number of nonzero elements. In addition, the RSS can be computed RSS E3 J J7
TMC LYmO Ejev Jj Ko + R.
Discussions
might seem surprising that, despite apparent model misspecification due the use of quadratic loss. the f LinR estimator can still correctly infer the structure with the same order of sample complexity as LogRo 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 L,-Li estimator is given as
(SDSk)
Σ (sssa) Jxl, k 140
(27)
where (.) and J&l represent average w.r.t. the Boltzmann distribution 7 and the sub-gradient of [Je], respectively. In the paramagnetic phase, (sts)) decays in its magnitude exponentially w.r.t. the

---

distance of sites i and j. This guarantees that once the connections JR of sites in the first nearest neighbor set are given so that
(SoSk)
> (sySk)] + Asign(Ja), Ak E v jEV
(28)
holds, the other conditions are automatically satisfied by setting all the other connections that not from to zero. For appropriate choice of X, (28) has solutions of sign( (Jz) > O, Ak E V. Namely yk E V the estimate of J& has the same sign as the true value J This implies that average the L--in estimator can successfully recover the network structure up to the connection signs if is chosen appropriately.
The key of the above argument that (s,8,) 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 & Li-Li 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 & LLin::: LooL also exhibits similar behavior unless post-thresholding
is used, as reported in [43]
5 Experimental results
RSS, N=200, À =0.3
Precision, N=200, =0.3
Recall, N=200, A =0.3
0.4
of throoryet
0.38
0.9
ingurimenntt -LogR
0.8
RSS
0.98
theorysf -log
Precision
Recall
0.8
0.6
0.7
0.34
expermmett( ( theoryst
experimectsf, Linßt
0.6
shouryyly Linlt
experiment:0 -L.gg
experimer tody dboory.f-LLoft
thoory.d-LogR
RSS, 2D grid, N=225, =0.15
Precision, 2D grid, N=225, 0.15
Recall, 2D grid, N2225, 0015
0.35
esperimentzff-Llle
theory:é, LinR
0.8
o.8
experineettt LogR
Precision
theorysl -Loght
0.6
Recall
0.8
0.4
experiment:f (, -Lnk
0.4
experimnnntt Linlt
0.15
theorryty-Lin
throey- fi -LinR expestment. ft Loooo theoey f. LogR
6.1
0.2
experiment:/ LLogt
Log
0.2
0.05
2
Figure 2: Theoretical and experimental results of RSS, Precision and Recall for RR graph and 2D grid using LinL and é LogR with different values of a M/N. The standard error bars obtained from 1000 random runs. good agreement between theory and experiment achieved, even small-size 2D grid graph with many loops. For more results, please refer to Appendix
In this section, we conduct numerical experiments verify the accuracy of the theoretical analysis. The experimental procedures are as follows. First, random graph G E 9N.d. K. is generated and the Ising model is defined on it. Then, the spin snapshots are obtained using the Metropolis Metropolis-Hastings Hastings algorithm
44 45| 46) in the same way as (7]. yielding the dataset DM We randomly choose a center spin So and infer its neighborhood using the &i-Li (4) and LogL [10] estimators. To obtain standard error bars, we repeat the sequence of operations 1000 times. The RR graph G E GN, with node degree d = 3 and coupling strength Ko = 0.4 considered, which satisfies
the paramagnetic condition (d 1) tanh? (Kn) < 1. The active couplings (Jig) have the same probability of taking both signs of +1 or
We first verify the precise non-asymptotic predictions of our method described in Sec 3.4 Fig. (upper figure) shows the replica and experimental results of RSS. Precision, Recall for N 200 with different values of a = M/N. It can be seen that for both LinLnR and E, LogR, there is a fairly
TToough this setting is different from the analysis where the nonzero couplings take uniform sign, the result be directly compared thanks to gauge symmetry [2]].

---

12
1p
18.00
10.0
13005
0231
ord iquareet
12
C:18 Pho0)
Eel Shrme
12
070
177
03
04
Figure 3: Left: critical scaling value co (A, Ko) Ka).Ko) of Li-Ri and L for the RR graph G E &x with 3, Ko = 0.4. Middle and Right: Precision and Recall for RR graph using Li-Li with À = 0.3. Experimental results are shown for N = 200, 400, 800, 1600, 3200. When e > cu (A, Ko) (co (A = 0.3, Ko) = 19.41 in this case), the Precision increases consistently with N and approaches I as N while it decreases consistently with N when < co (A, Ko). The Recall inereases consistently and approach to as + c0. For more results, please refer to Appendix G
good agreement between the theoretical predictions and experimental results, even for small N 200 and small a (equivalently small M), verifying the correctness of the replica analysis. Interestingly, quantitatively similar behavior between & LLinR and LLgR 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 [ the RSS much smaller for fooLog,, which is natural since the estimates of Logk 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 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 x 15 2D grid with uniform constant coupling Ko- = 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. & in Appendix G
Subsequently, the asymptotic result and sharpness of the critical scaling value co (A, Ko) Ko) in (22 are evaluated. First, Fig. (left) shows comparison of cu (A. Ko) between Lg-LinR and e LogR for the RR graph G IN,d, when d = 3, Ko 0.4, indicating similar behavior of LogL and L,LLinR. Then, we conduct experiments for M = clog N with different values of e around co (A, Kn), and investigate the trend of Precision and Recall N increases. When X 0.3, Fig. [3 (middle and right) show the results of Precision and Recall, respectively. As expected, the Precision increases consistently with N when e > co (A, Ko) and decreases consistently with N when e < co (A, Ko) while the Recall increases consistently and approaches to I as N - co, which verifies the sharpness of the critical scaling value prediction. The results for LoLLo,,, including the
case of x = 0.1 for both & LinR and E LogR, are shown in Fig. 9 and Fig. 10 in Appendix G
6 Conclusion
In this paper, we provide a unified statistical mechanics framework for the analysis of typical learning performances of E-rregularized M-estimators, LinL 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 accurate estimate of the typical sample complexity. It is revealed that, perhaps surprisingly, the misspecified é Li-LinR estimator is model selection consistent using M O (log N) samples, which is of the same order as LogR, Moreover, with 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 Liini (also & Lo-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 LLiR,, but also to other low-complexity estimators like LogL unless post-thresholding used 43 These assumptions restrict the applicability of the presented result, and thus overcoming such limitations will be important direction for future work.
10

---

Acknowledgements
This work was supported by JSPS KAKENHI Nos. 17H00764, 18K1 1463, and 19H01812, and JST CREST Grant Number JPMJCR1912, Japan.
References
[1] Martin J Wainwright and Michael Irwin Jordan. Graphical models, exponential families, and variational inference. Now Publishers Inc, 2008.
[2] Chau Nguyen, Riccardo Zecchina, and Johannes Berg. Inverse statistical problems: from the inverse ising problem to data science. Advances in Physics, 66(3):197-261, 2017.
[3] Erik Aurell and Magnus Ekeberg. Inverse ising inference using all the data. Physical review letters, 10899),0 2012.
[4] Ludovica Bachschmid Baanchhmid-Romnnnoo and Manfred Opper. Learning of couplings for random asym- metric kinetic ising models revisited: random correlation matrices and learning curves. Journal of Statistical Mechanics: Theory and Experiment, 2015(9):P09016, 2015.
[5] Johannes Berg. Statistical mechanics of the inverse ising problem and the optimal objective function. Journal of Statistical Mechanics: Theory and Experiment, 2017(8)0083402, 2017.
[6] Ludovica Bachschmid- Rooomn and Manfred Opper. A statistical physics approach to learning curves for the inverse ising problem. Journal of Statistical Mechanics: Theory and Experiment, 2017(6)003306, 2017.
[7] Alia Abbara, Yoshiyuki Kabashima, Tomoyuki Obuchi, and Yingying Xu. Learning performance inverse ising problems with sparse teacher couplings. Journal of Statistical Mechanics: Theory and Experiment, 2020(7):073402, 2020.
[8] Martin Wainwright, John D Lafferty, and Pradeep K Ravikumar, High-dimensional graphical model selection using f I-regularized logistic regression. In Advances neural information processing systems, pages 1465-1472, 2007.
[9] Holger Höfling and Robert Tibshirani. Estimation of sparse binary pairwise markov networks using seeudo-likelhoood Journal of Machine Learning Research, 10(4):883-906, 2009.
[10] Pradeep Ravikumar, Martin J Wainwright, John D Lafferty, et al. High-dimensional ising model selection using reguulrrie logistic regression. The Annals of Statistics, 38(3):127-1339, 2010.
[11] Narayana P Santhanam and Martin J Wainwright. Information-theoretic limits of selecting binary graphical models high dimensions. IEEE Transactions on Information Theory, 58(7):41174113, 2012.
[12] Aurélien Decelle and Federico Ricci-Tersenghi. Pseudolikelihood decimation algorithm im- proving the inference of the interaction network in general class of ising models. Physical review letters. 12((7):076033, 2014.
[13] Guy Bresler. Efficiently learning ising models on arbitrary graphs. In Proceedings of the forty-seventh annual ACM symposium on Theory of computing, pages 771-782, 2015.
[14] Marc Vuffray, Sidhant Misra, Andrey Lokhov, and Michael Chertkov. Interaction screening: Efficient and sample-optimal learning of ising models. In Advances in Neural Information Processing Systems, pages
2016.
[15] Andrey Y Lokhov, Mare Vuffray, Sidhant Misra, and Michael Chertkov. Optimal structure and parameter learning of ising models. Science advances, 1707917 2018.
[16] Shanshan Wu, Sujay Sanghavi, and Alexandros G Dimakis. Sparse logistic regression learns all discrete pairwise graphical models. Advances in Neural Information Processing Systems, 32:8071-80$3 2019.
11

---

[17] Guy Bresler, Elchanan Mossel, and Allan Sly. Reconstruction of markov random fields from samples: Some observations and algorithms. In Approximation, Randomization and Combina- torial Optimization. Algorithms and Techniques, pages 343-356. Springer, 20088
[18] Adarsh Prasad, Vishwak Srinivasan, Sivaraman Balakrishnan, and Pradeep Ravikumar. On learning ising models under huber's contamination model. Advances Neural Information Processing Systems, 33, 2020.
[19] Sahand N Negahban, Pradeep Ravikumar, Martin J Wainwright, and Bin Yu. A unified framework for digeesinnnii analysis of m-estimators with decomposable regularizers. Statistical science, 27(4):538-55) 557, 2012.
[20] Mare Mézard, Giorgio Parisi, and Miguel Angel Virasoro. Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, volume 9. World Scientific Publishing
Company, 1987.
[21] Hidetoshi Nishimori. Statistical physics of spin glasses and information processing: introduction. Number 111. Clarendon Press, 2001.
[22] Manfred Opper and David Saad. Advanced mean field methods: Theory and practice. MIT press, 2001.
[23] Marc Mezard and Andrea Montanari. Information, physics, and computation. Oxford University Press, 2009.
[24] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267-288, 1996.
[25] Mohsen Bayati and Andrea Montanari. The lasso risk for gaussian matrices. IEEE Transactions on Information Theory, 58(4):1997-2017, 2011.
[26] Sundeep Rangan, Alyson K Fletcher, and Vivek K Goyal. Asymptotic analysis of map esti- mation via the replica method and applications to compressed sensing. IEEE Transactions on Information Theory,
2012.
[27] Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi. Lasso with non-linear measurements is equivalent to one with linear measurements. Advances in Neural Information Processing Systems, 28:3420-3423 2015.
[28] Cédric Gerbelot, Alia Abbara, and Florent Krzakala. Asymptotic errors for convex penalized linear regression beyond gaussian matrices. arXiv preprint arXiv:2002.04337, 2020.
[29] Xiangming Meng, Tomoyuki Obuchi, and Yoshiyuki Kabashima. Structure learning in inverse ising problems using fz-regularized linear estimator. Journal of Statistical Mechanics: Theory and Experiment, 2021 2021.
[30] Julian Besag. Statistical analysis of non-lattice data. Journal of the Royal Statistical Society: Series D (The Statistician), 24 1975.
[31] David R Brillinger. A generalized linear model with gaussian regressor variables. Festschrift for Erich I. Lehmann, page 97-114. 1982.
[32] Yaniv Plan and Roman Vershynin. The generalized lasso with non-linear observations. IEEE Trunsuctions on information theory, 62(3): 1228-11528 2016.
[33] Yue Zhang, Weihong Guo, and Soumya Ray. On the consistency of feature selection with lasso for non-linear targets. In International Conference Machine Learning, pages 183-191. PMLR, 2016.
[34] Martin Genzel. High-dimensional estimation of structured signals from non-linear observations with general convex loss functions. IEEE Transactions on Information Theory, 630));1601-1619 2016.
[35] Nicolai Meinshausen, Peter Bühlmann, et al. High-dimensional graphs and variable selection with the lasso. The annals of statistics, 34(3): 1ee11144344422 2006.
12

---

[36] Peng Zhao and Bin Yu. On model selection consistency of lasso. The Journal of Machine Learning Research, 7:2541-2563, 2006.
[37] Leslie Valiant. theory of the learnable. Communications of the ACM, 273111:113--114. 1984.
[38] Andreas Engel and Christian Van den Broeck. Statistical mechanics of learning. Cambridge University Press, 2001.
[39] Federica Gerace, Bruno Loureiro, Florent Krzakala, Mare Mézard, and Lenka Zdeborová. Generalisation error in learning with random features and the hidden manifold model. In International Conference on Machine Learning, pages 3452-3462. PMLR, 2020.
[40] Madhu Advani and Surya Ganguli. Statistical mechanics of optimal convex inference in high dimensions. Physical Review X. 6(3):01033, 2016.
[41] Persi Diaconis and Mehrdad Shahshahani. On the eigenvalues random matrices. Journal of Applied Probubility, (1()):99662, 1994.
[42] Kurt Johansson. On random matrices from the compact classical groups. Annals of mathematics, pages 519-545, 1997.
[43] José Bento and Andrea Montanari. Which graphical models are difficult to learn? arXiv preprint 10.51:0910.5761, 2009.
[44] Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta Teller, and Edward Teller Equation of state calculations by fast computing machines. The journal of chemical physics, 1186):0077-109, 1953.
[45] W Keith Hastings. Monte carlo sampling methods using markov chains and their applications.
1970.
[46] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions on pattern analysis and machine intelligence, (6):721- 741, 1984.
[47] Federico ficci-Tersengh.. The bethe approximation for solving the inverse ising problem: a comparison with other inference methods. Journal of Statistical Mechanics: Theory and Experiment, 2012(05):P08015, 2012.
[48] Chau Nguyen and Johannes Berg. Bethe-peierls approximation and the inverse ising problem. Journal of Statistical Mechanics: Theory and Experiment, 2012 2012(03),P03004, 2012.
[49] Brendan D McKay. The expected eigenvalue distribution of large regular graph. Linear Algebra and Applications, 40:203-216, 1981.
Checklist
For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] See Section 6
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of theoretical results? [Yes] See Section 3.2
(b) Did you include complete proofs of all theoretical results? [Yes] See Section and Appendices A
3. If you experiments...
13

---

(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either the supplemental material or URL)? [Yes] See Section [
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [N/A]
(c) Did you report error bars (e.g., with respect the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [N/A]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [N/A]
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]
(d) Did you discuss whether and how consent obtained from people whose data you's re using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsoureing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
14


Reviews:
[{"review_id": "v9NfNWCAq1", "reviewer": "Reviewer_Y3wG", "summary": "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.\nThey demonstrate the sample complexity is such that O(log N) samples are enough to correctly identify the right neighborhood of a generic variable.\nThe 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.\nThe methodology of the present work can be extended also to other estimators.", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The present work looks scientifically valid and robust. The authors seem to master the replica computation. The results obtained are very original and are presented in a very clear way.\nI consider these results very important as they show how the replica method can be used to compute the relevant estimators and the sample complexity. The numerical simulations run for not too large values of the system size and number of samples convincingly support the validity of the theory.\n\nI have only one main concern. The quadratic loss function does not look like the appropriate function to minimize in order to correctly infer the interaction graph. Indeed on variables where the field is very large, even aligning the spin to the field, one gets a very large penalty.\nContrary to logistic regression and interaction screening, where the loss functions go to zero at large values of the argument.\nThe authors should explain how and why the quadratic loss function that penalizes a spin that aligns with a very strong field can provide the right answer to the Ising model selection problem.\nIs it maybe because the authors are considering only the paramagnetic phase and so the fields are not very large? If it is so, it would be worth discussing it, also to clarify the limits of this approach.\n\nIn the numerical simulations, the authors study only the spin-glass model, that is couplings with a random sign.\nWhy not studying also the ferromagnetic model? In the latter case probably the coupling constant $K_0$ could be made larger and the theory tested also in the case where the model has a ferromagnetic phase.\n\nThe results for LinR and LogR are always so close maybe because the authors are studying the problem in a regime where the 2 losses are very similar. It would be very interesting to see what happens in the other regimes. How different would be the two estimators and which one is closer to the truth.\n\nMinor point: in several places the authors write \"intersection screening\" instead of \"interaction screening\"", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I strongly suggest the authors discuss more the limitations of their method.\nIn particular, what happens approaching the critical point.\n\nI don't see any potential negative societal impact.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "uNqjHklrrRI", "reviewer": "Reviewer_8TQ3", "summary": "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.\nMore specifically, the authors focus on Ising models associated with typical random regular graphs in the paramagnetic phase.\nIdeas from statistical physics are borrowed to study the performance of an estimator derived from an $\\ell_1$-regularized linear regression problem.\nIt is claimed that using $M = \\mathcal{O}(\\log(N))$ samples is enough to obtain model selection consistency.\nNon-asymptotic results are also given, and some numerical experiments are presented to support the findings.", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### Originality\n\n#### Are the tasks or methods new?\n\nUsing an $\\ell_1$-regularized linear regression estimator to recover the underlying graph structure of Ising models\n\n#### Is the work a novel combination of well-known techniques? (This can be valuable!)\n\nThe authors borrow classical ideas from statistical physics (replica method, free energy) and optimization (leave one out) to study the performance of their estimator\n\n#### Is it clear how this work differs from previous contributions?\n\nThe difference with previous work is explained in the introduction section.\n\n#### Is related work adequately cited?\n\nRelated work is adequately cited, but the ideas and notions from statistical physics must be explained more carefully for the ML community to integrate it and build on it.\nPlease give references to the sampling algorithms used to generate the datasets $\\mathcal{D}_M$\n\n---\n\n### Quality\n\n#### Is the submission technically sound?\n\nI didn't carefully check the mathematical details.\n\n#### Are claims well supported (e.g., by theoretical analysis or experimental results)?\n\nThe supplementary material contains the proofs of the claims.\n\n#### Are the methods used appropriate?\n\nThe methods used seem appropriate.\n\n#### Is this a complete piece of work or work in progress?\n\nThe submission is work in progress regarding the scope of applicability of the results.\n\n#### Are the authors careful and honest about evaluating both the strengths and weaknesses of their work?\n\nIt is not clear whether the assumptions made (random regular graphs, paramagnetic phase) are very restrictive and how these assumptions differ from other works like l1-LogR estimator for example.\n\n---\n\n### Clarity\n\n#### Is the submission clearly written?\n\nThe paper is well written.\n\n#### Does it adequately inform the reader? (Note that a superbly written paper provides enough information for an expert reader to reproduce its results.)\n\nPlease\n\n* define what a random regular graph with coupling strength $K_0$ is\n* explain the choice for the specific class of random regular graphs in the paramagnetic phase\n* define the notions of paramagnetic phase and model consistency earlier\n* detail the misspecification of the l1-LinR model\n* define what a 'typical' random regular graph is. It seems a critical insight to call for the replica method (l148)\n* explain why the overall analysis is more difficult with a larger class of graphs than 'typical' random regular graphs\n* explain why Ansatz 2 can only be partly checked, and if it is a\n* give more details on the free energy and explain why this is a key quantity\n* clarify what l264 \"spin snapshots are obtained with MC sampling\". Are the samples exactly distributed according to the target Ising model or are they obtained using an MCMC techniques.\n* clarify the mismatch between the claim l267 \"to obtain standard error bars, we repeat the sequence of operations many times\" and Figure 1 were the \"error bars were obtained with only 5 runs\".\n\n* define the notations\n  + $J_j$ l164\n  + $\\| \\cdot \\|_1$ equation (2)\n  + $Extr$ l182\n* clarify l193 \"$I$ in (15) is equal to $[I]_O$\"\n* consider adopting the notation $m$ in place of $\\mu$\n* make the attributes of Figure 2 bigger\n\n---\n\n### Significance\n\n#### Are the results important?\n\nIt is an interesting study of a natural estimator derived from an l1-regularized linear regression problem for the Ising model selection problem.\n\n#### Are others (researchers or practitioners) likely to use the ideas or build on them?\n\nThe presentation of the ideas borrowed from statistical physics lacks of clarity for the ML community to build on it.\nBesides, the emphasis l70 that this work \"might be of particular interest for the statistical mechanics community\" may seem surprising for a submission to a ML conference.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "NTR", "review_text": "", "extra_scores": {"ethical_concerns": "NTR", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "rtUIEG_Fkd5", "reviewer": "Reviewer_QMbt", "summary": "The manuscript contains a theoretical investigation of the Ising model selection problem. \nThe task is to recover the support (the graph topology) of the interaction network from observations of Ising configurations. \n\nAn objective function in the M-estimator family, named L1-LinR by the authors, is analyzed. \nIt contains the typical quadratic data-dependent loss from linear regression and an L1 regularization term, \ntherefore it is a mismatched objective compared to L1-LogR (pseudolikelihood) and the Interaction Screening.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is clear and well written. \n\nThe analysis is carried out using the replica technique,\nwith a quenched disorder given by the random\nregular graph ensemble. Couplings are uniform and ferromagnetic.\nThe paper mostly builds on the analyses of Ref. [7] and Ref. [25]\n\n- The authors show the mismatched L1-linR is a consistent \n  estimator with performance similar to L1-logR\n\n- The authors provide the “typical” sample complexity, as opposed to\n  worst-case bounds given in the literature \n\n- They provide nonasymptotic results replacing an expectation\n  with sampling in the state evolution equations.\n  \n- The authors manage to carry out a quite impressive replica calculation, \n  relying on a few ansatz and approximations that could be of general\n  interest.  \n\n- They numerically validate their theoretical predictions for recall and precision\n  finding excellent agreement with experiments. \n  Moreover, also results on 2d grids, which are beyond the scope of the \n  technique employed, are surprisingly accurate.\n  \n\nOther comments:\n\n\n- Results are obtained through a saddle point equation, but some non-asymptotic effects are \nare retained using the finite sampling of eq. 25 instead of the full expectation. It is not\nclear how these non-asymptotic effects are more relevant than e.g. gaussian fluctuations around \nthe saddle point that are not taken into account.\n\n- Related to the above question: Does the extremization condition of the free energy functional come from the\nlarge M or N limit or from the zero-temperature limit?\nFrom the appendix I assume the latter, but this could be made more explicit.\n\n- Comment on the large M and N limit at e.g. fixed alpha=M/N. Below eq. 19 the authors argue that 2 different scalings are involved in\nthe equation of state. It is not clear at this point what is expected in the infinite system size limit, what are \nthe possible sensible scaling for M and N.\n\n- Computation of eq. 25 is the only one in which not O(1) terms are involved in the EOS. \n Would it be possible to express it analytically through a 1/M expansion instead?\n\nMinor comments:\n- Sec 3: clarify that J is a vector (not the full matrix)\n\n- Ansatz 1: specify that the \\hat{J} is a random variable over the realization of the dataset \n(if I understand correctly)\n\n- Eq. 13, lim \\beta \\to \\infty also on the r.h.s.? \n\n- The authors should comment on the additional complexity given by using different losses, such as the logistic one. \nIt seems to me that this change would only affect the energetic term in equation 16, and correspondingly the computation\nof the F order parameter (and possibly others), without much complexity increase.\n\n- On line 182: (s0, s\\psi) seem to be distributed as P_{ising} once the non-active set is marginalized out. It is not clear though why on the\nsame line this marginal distribution is assumed to be proportional to exp(s0\\sum_j∈Ψ J^* j sj) which seems wrong at first sight. \nI assume cavity fields coming from the RRG are involved, those fields being zero in the paramagnetic phase.\nHow to sample from this marginal distribution should be explicit (e.g. after line 447 of the appendix), it's hard for the\nreader to dig the relevant piece of information from Ref. [7].\n\n- Also, since J^* is a matrix, the notation needs to be made more consistent.   \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "A major critic is that L1-LinR objective has little relevance in \nthe Ising model selection context, where practitioners use \nmatched M-estimators. \nNonetheless, given the quality of the results obtained and the complexity\nof these kinds of analysis, I think the paper would represent a valid contribution\nto NeurIPS literature. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "omNXDMYo0wz", "reviewer": "Reviewer_gune", "summary": "This paper investigates the L1-regularized linear regression (L1-LinR) for learning Ising models. The authors apply the replica methods from statistical physics and obtain meaningful learning performance results for random regular graphs in the paramagnetic phase.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The Ising model is one of the most popular statistical models that are used in various fields of computer science. It belongs to the class of undirected graphical models (also called Markov random fields) for modeling high dimensional binary data. This paper studies the problem of Ising model selection, which tries to learn the edge set of the Ising model from iid samples. The authors focus on the case of random regular graphs in the paramagnetic phase. As a main result, the authors propose the L1-regularized linear regression (L1-LinR) and show that it has an optimal logarithm sample complexity for typical random regular graphs. This estimator is simpler and computationally more efficient compared with other estimators like L1-LogR and IS. It is also representative for model misspecification. The authors provide detailed analysis for the learning performance of L1-LinR in the non-asymptotic case with strong experimental supports. The methods in this paper are potentially applicable to other estimators and models.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "It would be helpful to include some discussions about why the L1-LinR estimators perform well (only) in the paramagnetic phase. How good is L1-LinR beyond the paramagnetic phase? Would it become inaccurate since the replica methods fail? Some intuition from statistical physics is sufficient.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1"}, "extra_text": {}}, {"review_id": "aoxm8TcyqXA", "reviewer": "Reviewer_gfUX", "summary": "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. ", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper provides a sharp analysis of an interesting problem. On the positive side, I think there is value is obtaining working through the details of the replica method for this problem. Also, the modifications used for the moderate dimensional setting are interesting and could be useful for other problem settings as well. \n\nOne weakness of the paper is that the results are not rigorous because they depend on the replica method. In the context of NeurIPS, I think there is a higher bar (in terms of the novelty of the insight and level of interest in the particular problem) for acceptance of a non-rigorous theoretical analysis. While I think the results in this paper and worth publishing, it is not clear to me if they meet the very high bar given by NeurIPS. \n\nI am also a bit unclear about whether the spectral density used in (18) is something that is known in practice. In Appendix F it seems the authors derive this density for the problem of interest, but I am not sure if there are any more assumptions that are needed here. In any event, I think it would help to say a bit more what is and is not known about this density in the main text since it plays a prominent role in the analysis. \n\nSome technical comments and questions are the following:\n\n- In Ansantz 1 the term \\bar{J} is defined as the the mean-valued of the estimator. I am not sure what this means. Is this the conditional mean given that j is in the neighborhood of index 0?\n\n- It seems that the discussion of the CLT  in lines 172-174 is meaningful only if the degree is large.\n\n- In the context of Ansantz 1, does the expression (21) also describe the noise term for the estimator in the active set?\n\n- With respect to the upper bound on lambda discussed in lines 235-237, would it be possible to consider a two-stage approach where one the final estimate is based on adding certain ``residuals'' to the estimates obtained from the ell-one regularized problem and then thresholding a second time. Such approaches have been shown to work for support recovery in sparse linear regression problems. \n\n- How the regular graphs sampled for the numerical experiments?\n\n- What are the values of alpha used in Figure 1?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "These are adequately addressed. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "RyqUPMjuF5Y", "reviewer": "Reviewer_hC3k", "summary": "This paper analyses the behaviour of $\\ell_{1}$ regularized linear regression for estimating Ising models, as compared to the canonical approach of $\\ell_{1}$ regularized logistic regression. The authors use tools from statistical mechanics to arrive at their result, which is in the paramagnetic regime, and corroborate some of their claims using experiments.", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### Originality\nUsing techniques from statistical mechanics to show results relevant in the statistical learning theory literature is a nice unification. Some of the claims in the work seem to be unsubstantiated as I describe in later comments. \n\n### Quality\n- I found it interesting that $c(\\lambda, K_0)$ decreases with $d$, standard results have $M > max(d^{2}, \\lambda^{-2}) \\log(N)$ to estimate graph structure. Furthermore, the precision, recall can be easily computed once we've learned the graph, and it appears that the authors do the same from Algorithm 1 by estimating $J$ which is sufficient, so it is not clear where the benefit in their approach arises from. Non-asymptotic rates for logistic regression have been derived in Wu et al (2019) more recently.\n\n### Clarity\nI think the writing can be significantly improved to be more accessible to members of the theoretical machine learning community. A lot of the notation and techniques used are very esoteric, and specialized to the physics community (even remarked by the authors). As someone familiar with statistical learning theory literature more than statistical mechanics, my difficulty in reading the manuscript is reflected in my detailed comments below.\n\n### Significance\nWhile it is interesting to see such a result, it is unclear about the utility of the result w.r.t. predicting the recall and precision. Often times, we intend on finding the edges of the graph rather than the recall and precision alone, and estimating the graph gives us this automatically.\n\nAdditional comments below:\n\nSection 1 and 2:\n- IS is interaction screening, and not intersection screening as the authors refer multiple times in the introduction.\n- The model misspecification appears to be deliberate, because if we know that the datapoints are binary valued, then a linear regression approach is probably too much.\n- lasso vs. $\\ell_{1}$ regularized logistic regression comparison: why is the former \"computationally efficient\" than the latter? Are there any computational complexity analyses to validate this statement. Practically, software like glmnet have bridged any gap between regularized linear and logistic regression.\n- The authors address to certain things as \"typical\" in Section 1.1. What is \"typical\"?\n- Line 40: assessing learning performance is not very novel. refer Wu et al (2019) for a comparison.\n- Typo Line 92: IOS losses --> ISO losses\n\nSection 3: for someone with a statistical learning theory background, a lot of the notation and use of techniques are alien, and my review of this section will reflect this inexperience.\n\n- It would have been nice if the authors provided a primer for this. For instance, it is unclear to me what the exact purpose of the Hamiltonian and Boltzmann distribution derived in eqns. (6) and (7) are, and subsequent propositions about convergence to a point wise measure when the temperature tends to 0. Additionally the notation in line 143 defining the partition function is also conflicting because it is unclear if that definition and the one in eqn (1) are equivalent for Ising models.\n\n- It appears that the self-averaging property is similar to the law of large numbers. However, it is not clear why a result like this would hold only when $N \\to \\infty$.\n\n- The expectation under a dataset is defined in a strange way: are we considering multiple draws of M samples from the Ising model (as the summation indicates)? If so, a) how many draws?\n\n- I feel that the replica method could have been briefly described, or at least the essence of it. Eqn (9) has undefined terminology as well: $Z^{n}$, derivatives of the log partition function w.r.t. n. The authors describe the RHS of the equation in the paragraph, but it is not immediately apparent why this is relevant. On further reading, it seems as though the authors have made use of this technique due to its prior use in disorder systems and in a generalisation analysis. The citation to the paper by Bresler on learning Ising models is irrelevant (there is no mention of \"replica\" in that paper).\n\n- A similar issue appears with Eqn (10) i.e., undefined notation: $||J^{a}||\\_{1}$, $J^{a}\\_{j}$. It would have been useful to include a citation for this equation, because a derivation isn't available and $Z^{n}$ has not been defined prior.\n\n- The Ansatz and subsequent corollaries are very hard to read due to excessive notation and terminology. Line 173 has $w_{i}$ where the index of summation is $j$. The purpose of $Q$ and $q$ are unclear, the introduction of $a$ and $b$ on Equation 12 is also unprecedented. Eqn (13) introduces $\\texttt{Extr}$ which is previously undefined too. Eqn (14) -- (16) are unhelpful, because it is not clear what the purpose of those equations are. It is also unclear as to what the estimator is in this case: is it the L1-LinR estimator that is being referred to here?\n\n- Line 189: why is this true? Is this a valid assumption? The Ansatz appears to be partly checked, does this mean that the claim is verified but not proven?\n\n- Unfortunately, despite the authors' effort, I find the majority of page 6 inaccessible, as many of the definitions and concepts are unknown in the standard statistical learning theory. Are the EOS equations directly relevant to the results in this work? Additionally, why are the authors focusing on the temperature tending to 0 case for a result in the paramagnetic regime?\n\n- Eqn (20) is derived using the Ansatz stated earlier and Eqn (17). However, Eqn (17) talks about the expression in the limit $\\beta \\to \\infty$. Additionally, how expensive is solving the EOS, because they appear to be a set of highly nonlinear equations.\n\nSection 4 contain experimental results, and here it would have been interesting to note the variation of their results with $d$, which the authors seem to have not considered.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors describe the only limitation of the work being that the paramagnetic regime is too constrained. I agree with this assessment.\n\nThere are no direct societal impacts of this work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors in this paper study the learning of a Boltzmann machine, aka the Inverse Ising model, a classical problem in graphical models and Markov Random Fields. They propose an analysis of the performance of the l1-regularised linear regression estimator in finding the underlying non-zero coefficients. The analysis is performed using a non-rigourous but powerful heuristic approach from statistical mechanics, the replica method, that has been used in many other machine learning problems, in idealistic situations. The author confirmed the theoretical result by running a large number of numerical simulations. They show the sample complexity is such that O(log N) samples are enough to correctly identify the right neighbourhood of a generic variable. The theory presented seems to be powerful enough to predict both  the precision and the recall rates for finite dimension, and seems to give fairly good prediction for graphs with many loops. Interestingly, the methodology of the present work can be apparently generalised and extended also to other estimators.

The review process was intense, with six reviewer and a large number of forum replies. All reviewers testified of the quality of presented results. From a technical point of view, computation contained in this manuscript was found to  involved and complex, yet it was found provides a non-trivial analytical result. The comparison with the numerics however strongly supports the claim that the replica computation is giving the correct answer, even though it is not rigorous. This paper was thus judged to represent a non-trivial contribution from statistical physics to machine learning. However, the use of the non-rigorous replica method and its acrobatic non-rigorous mathematics was judged to be sometimes dazzling.

The rebuttal saw a number of discussions, and most reviewers agreed that the paper was impressive, and score were increased during the process. A criticism that remained, though, from a minority of reviewers, was that the derivation is likely inaccessible for non-experts. While this is true, it can be said of many theoretical at Neurips, including rigorous ones, and there is a long standing tradition of welcoming such non-rigorous paper at Neurips. 

Given the agreement on the quality of the results, the very good ranking & grading of the paper, the extensive numerical simulations that confirm the validity of theory, we believe the paper to be largely worthy of a publication at a venue like Neurips, and recommend acceptance.


Decision: Accept (Poster)