Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 8: Top 50% of accepted NeurIPS papers, clear accept
review_content: 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.
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.
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.
I 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.

I 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.
Contrary to logistic regression and interaction screening, where the loss functions go to zero at large values of the argument.
The 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.
Is 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.

In the numerical simulations, the authors study only the spin-glass model, that is couplings with a random sign.
Why 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.

The 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.

Minor point: in several places the authors write "intersection screening" instead of "interaction screening"
limitations_and_societal_impact: I strongly suggest the authors discuss more the limitations of their method.
In particular, what happens approaching the critical point.

I don't see any potential negative societal impact.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.
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.
main_review: ### Originality

#### Are the tasks or methods new?

Using an $\ell_1$-regularized linear regression estimator to recover the underlying graph structure of Ising models

#### Is the work a novel combination of well-known techniques? (This can be valuable!)

The authors borrow classical ideas from statistical physics (replica method, free energy) and optimization (leave one out) to study the performance of their estimator

#### Is it clear how this work differs from previous contributions?

The difference with previous work is explained in the introduction section.

#### Is related work adequately cited?

Related 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.
Please give references to the sampling algorithms used to generate the datasets $\mathcal{D}_M$

---

### Quality

#### Is the submission technically sound?

I didn't carefully check the mathematical details.

#### Are claims well supported (e.g., by theoretical analysis or experimental results)?

The supplementary material contains the proofs of the claims.

#### Are the methods used appropriate?

The methods used seem appropriate.

#### Is this a complete piece of work or work in progress?

The submission is work in progress regarding the scope of applicability of the results.

#### Are the authors careful and honest about evaluating both the strengths and weaknesses of their work?

It 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.

---

### Clarity

#### Is the submission clearly written?

The paper is well written.

#### Does it adequately inform the reader? (Note that a superbly written paper provides enough information for an expert reader to reproduce its results.)

Please

* define what a random regular graph with coupling strength $K_0$ is
* explain the choice for the specific class of random regular graphs in the paramagnetic phase
* define the notions of paramagnetic phase and model consistency earlier
* detail the misspecification of the l1-LinR model
* define what a 'typical' random regular graph is. It seems a critical insight to call for the replica method (l148)
* explain why the overall analysis is more difficult with a larger class of graphs than 'typical' random regular graphs
* explain why Ansatz 2 can only be partly checked, and if it is a
* give more details on the free energy and explain why this is a key quantity
* 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.
* 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".

* define the notations
  + $J_j$ l164
  + $\| \cdot \|_1$ equation (2)
  + $Extr$ l182
* clarify l193 "$I$ in (15) is equal to $[I]_O$"
* consider adopting the notation $m$ in place of $\mu$
* make the attributes of Figure 2 bigger

---

### Significance

#### Are the results important?

It is an interesting study of a natural estimator derived from an l1-regularized linear regression problem for the Ising model selection problem.

#### Are others (researchers or practitioners) likely to use the ideas or build on them?

The presentation of the ideas borrowed from statistical physics lacks of clarity for the ML community to build on it.
Besides, 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.
limitations_and_societal_impact: NTR
review_title: null

final_score: rating: 5: Marginally below the acceptance threshold
review_content: 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. 
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. 

One 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. 

I 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. 

Some technical comments and questions are the following:

- 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?

- It seems that the discussion of the CLT  in lines 172-174 is meaningful only if the degree is large.

- In the context of Ansantz 1, does the expression (21) also describe the noise term for the estimator in the active set?

- 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. 

- How the regular graphs sampled for the numerical experiments?

- What are the values of alpha used in Figure 1?
limitations_and_societal_impact: These are adequately addressed.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: summary: 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.

main_review: The paper is clear and well written. 

The analysis is carried out using the replica technique,
with a quenched disorder given by the random
regular graph ensemble. Couplings are uniform and ferromagnetic.
The paper mostly builds on the analyses of Ref. [7] and Ref. [25]

- The authors show the mismatched L1-linR is a consistent 
  estimator with performance similar to L1-logR

- The authors provide the “typical” sample complexity, as opposed to
  worst-case bounds given in the literature 

- They provide nonasymptotic results replacing an expectation
  with sampling in the state evolution equations.
  
- The authors manage to carry out a quite impressive replica calculation, 
  relying on a few ansatz and approximations that could be of general
  interest.  

- They numerically validate their theoretical predictions for recall and precision
  finding excellent agreement with experiments. 
  Moreover, also results on 2d grids, which are beyond the scope of the 
  technique employed, are surprisingly accurate.
  

Other comments:


- Results are obtained through a saddle point equation, but some non-asymptotic effects are 
are retained using the finite sampling of eq. 25 instead of the full expectation. It is not
clear how these non-asymptotic effects are more relevant than e.g. gaussian fluctuations around 
the saddle point that are not taken into account.

- Related to the above question: Does the extremization condition of the free energy functional come from the
large M or N limit or from the zero-temperature limit?
From the appendix I assume the latter, but this could be made more explicit.

- 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
the equation of state. It is not clear at this point what is expected in the infinite system size limit, what are 
the possible sensible scaling for M and N.

- Computation of eq. 25 is the only one in which not O(1) terms are involved in the EOS. 
 Would it be possible to express it analytically through a 1/M expansion instead?

Minor comments:
- Sec 3: clarify that J is a vector (not the full matrix)

- Ansatz 1: specify that the \hat{J} is a random variable over the realization of the dataset 
(if I understand correctly)

- Eq. 13, lim \beta \to \infty also on the r.h.s.? 

- The authors should comment on the additional complexity given by using different losses, such as the logistic one. 
It seems to me that this change would only affect the energetic term in equation 16, and correspondingly the computation
of the F order parameter (and possibly others), without much complexity increase.

- 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
same line this marginal distribution is assumed to be proportional to exp(s0\sum_j∈Ψ J^* j sj) which seems wrong at first sight. 
I assume cavity fields coming from the RRG are involved, those fields being zero in the paramagnetic phase.
How to sample from this marginal distribution should be explicit (e.g. after line 447 of the appendix), it's hard for the
reader to dig the relevant piece of information from Ref. [7].

- Also, since J^* is a matrix, the notation needs to be made more consistent.   

limitations_and_societal_impact: A major critic is that L1-LinR objective has little relevance in 
the Ising model selection context, where practitioners use 
matched M-estimators. 
Nonetheless, given the quality of the results obtained and the complexity
of these kinds of analysis, I think the paper would represent a valid contribution
to NeurIPS literature.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.
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.
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_title: null

final_score: rating: 5: Marginally below the acceptance threshold
review_content: 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.
main_review: ### Originality
Using 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. 

### Quality
- 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.

### Clarity
I 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.

### Significance
While 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.

Additional comments below:

Section 1 and 2:
- IS is interaction screening, and not intersection screening as the authors refer multiple times in the introduction.
- 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.
- 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.
- The authors address to certain things as "typical" in Section 1.1. What is "typical"?
- Line 40: assessing learning performance is not very novel. refer Wu et al (2019) for a comparison.
- Typo Line 92: IOS losses --> ISO losses

Section 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.

- 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.

- 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$.

- 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?

- 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).

- 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.

- 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?

- 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?

- 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?

- 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.

Section 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.
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.

There are no direct societal impacts of this work.
review_title: null

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)