# Learning performance in inverse Ising problems with sparse teacher couplings

Alia Abbara, Yoshiyuki Kabashima, Tomoyuki Obuchi, Yingying Xu

► To cite this version:

Alia Abbara, Yoshiyuki Kabashima, Tomoyuki Obuchi, Yingying Xu. Learning performance in inverse Ising problems with sparse teacher couplings. Journal of Statistical Mechanics: Theory and Experiment, 2020, 2020, pp.073402. 10.1088/1742-5468/ab8c3a. hal-04986028

HAL Id: hal-04986028

https://hal.science/hal-04986028v1

Submitted on 11 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
PAPER • OPEN ACCESS

Learning performance in inverse Ising problems
with sparse teacher couplings

To cite this article: Alia Abbara et al J. Stat. Mech. (2020) 073402

View the [article online](#) for updates and enhancements.

You may also like

- Statistical mechanics of the inverse Ising problem and the optimal objective function
  Johannes Berg

- A noniterative solution to the inverse Ising problem using a convex upper bound on the partition function
  Takashi Sano

- Ising models for neural activity inferred via selective cluster expansion: structural and coding properties
  John Barton and Simona Cocco
# Learning performance in inverse Ising problems with sparse teacher couplings

Alia Abbara¹, Yoshiyuki Kabashima²,³, Tomoyuki Obuchi⁴,<fnref n="6" /> and Yingying Xu²,⁵

¹ Laboratoire de Physique de l'Ecole Normale Supérieure, 24 rue Lhomond, 75005 Paris, France

² Department of Mathematical and Computing Science, Tokyo Institute of Technology, 2-12-1, Ookayama, Meguro-ku, Tokyo, Japan

³ Institute for Physics of Intelligence & Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

⁴ Department of Systems Science, Graduate School of Informatics, Kyoto University, Yoshida Hon-machi, Sakyo-ku, Kyoto-shi, Kyoto 606-8501, Japan

⁵ Mathematical Statistics Team, RIKEN Center for Advanced Intelligence Project, Japan

E-mail: obuchi@i.kyoto-u.ac.jp

Received 25 December 2019

Accepted for publication 15 April 2020

Published 3 July 2020

<figure><img src="image_2.png" /></figure>

CrossMark

Online at stacks.iop.org/JSTAT/2020/073402

https://doi.org/10.1088/1742-5468/ab8c3a

**Abstract.** We investigate the learning performance of the pseudolikelihood maximization method for inverse Ising problems. In the teacher-student scenario under the assumption that the teacher's couplings are sparse and the student does not know the graphical structure, the learning curve and order parameters are assessed in the typical case using the replica and cavity methods from statistical mechanics. Our formulation is also applicable to a certain class of cost functions having locality; the standard likelihood does not belong

<figure><img src="image_3.png" /></figure>

<footnote n="6">⁶ Author to whom any correspondence should be addressed.</footnote>
to that class. The derived analytical formulas indicate that the perfect inference of the presence/absence of the teacher's couplings is possible in the thermodynamic limit taking the number of spins $N$ as infinity while keeping the dataset size $M$ proportional to $N$, as long as $\alpha = M/N > 2$. Meanwhile, the formulas also show that the estimated coupling values corresponding to the truly existing ones in the teacher tend to be overestimated in the absolute value, manifesting the presence of estimation bias. These results are considered to be exact in the thermodynamic limit on locally tree-like networks, such as the regular random or Erdős–Rényi graphs. Numerical simulation results fully support the theoretical predictions. Additional biases in the estimators on loopy graphs are also discussed.

**Keywords:** learning theory, cavity and replica method, statistical inference, machine learning

## Contents

1. Introduction 3

2. Formulation 4

2.1. Inference framework 4

2.1.1. Teacher-student scenario 5

2.2. Statistical mechanical analysis 6

2.2.1. General framework 6

2.2.2. Revisit the fully-connected case 8

3. Details of the sparsely-connected case 10

3.1. Ansatz for the sparse case 10

3.2. Direct problem's properties 14

3.2.1. Marginal distribution of the teacher model 15

3.2.2. Inverse correlation function 15

3.3. Applicable range of the ansatz 15

4. Numerical experiments 18

4.1. RR graph case 20

4.2. ER graph case 23

4.3. Square lattice case for comparison 24

4.4. Comparison with interaction screening 25

5. Summary and discussion 26

Acknowledgment 28

Appendix A. Computations for L and S 28

Appendix B. Derivation of macroscopic parameters R and $\rho\rho$ 31
# Appendix C. Details of numerical experiments 31

## References 32

# 1. Introduction

Inference based on the classical Ising model is called the *inverse Ising problem* or *Boltzmann machine learning*, which is attracting more and more attention with the increasing interest in machine learning technologies. One recent application spurring this trend is for retinal neurons [1, 2], and subsequent applications to a wide range of systems have been conducted [3–10], showing the potential usefulness of the inverse Ising framework.

However, it is difficult to compute standard estimators such as the maximum likelihood estimator in this framework when the system size is large. Thus, certain approximations and/or algorithms must be tailored to ease this difficulty and meet the demands of advanced applications, which have been attempted in previous studies [11–23]. One of the most effective examples is the pseudolikelihood method [11, 20]. This method approximates the likelihood function by the product of conditional likelihood functions, each of which is for a single random variable conditioned by its neighboring variables. This is useful for a wide range of problems defined on graphical models, and enables us to treat large systems because local couplings directly connected to a focused single random variable are isolated from the other couplings, and thus can be estimated independently. This local nature is sometimes referred to as *local learning* [24].

Another benefit of local learning is its theoretical tractability in high-dimensional settings. Recent theoretical analyses based on the replica method revealed the tight limit of inference accuracy in the thermodynamic limit where the dataset size is comparable to model dimensionality [24–26]. This provides a firm theoretical basis for the inverse Ising framework.

Previous studies of [24–26] focused on fully-connected Ising models. In high-dimensional settings, however, sparsely-connected models are more interesting because the inference accuracy is expected to be much better than the dense case, and the inferred estimator is expected to have clearer interpretations owing to the sparsity. The present paper deals with this case. We investigate the so-called teacher-student scenario using the replica method by drawing on previous studies [24–26], but refine the theoretical treatment in [24] for dealing with the teacher with sparse connections. The cavity method is used but the cavity field is assumed to consist of separately treatable signal and noise; the functional form of the associated probability distribution is hypothesized. The network structure of couplings is assumed to be locally tree-like, and our theoretical result is expected to be exact on those networks in the thermodynamic limit. To check the accuracy, numerical experiments are also conducted.

The remainder of this paper is organized as follows. In section 2, we review the inverse Ising framework and the statistical mechanical formulation to analyze inference accuracy. We also briefly review the analysis of the fully-connected case given in [24] to demonstrate how the formulation is applied. In section 3, we extend the method in the statistical mechanical formulation to the sparsely-connected case. Some theoretical
treatments developed in [24] are refined and our ansatz to deal with the teacher model with sparse connections is stated. In section 4, our numerical analysis results are provided and compared with our theoretical results to check their accuracy. The final section presents a summary and discussion.

## 2. Formulation

In this section, we briefly review the inverse Ising framework and two associated inference methods: the maximum likelihood (ML) and pseudolikelihood (PL) methods. The class of local learning is also introduced and explained. In addition, the general framework of the statistical-mechanical analysis proposed in [24] is presented. We emphasize its technically important points and review the case of the fully-connected Ising model [24] to demonstrate its application.

### 2.1. Inference framework

Let us consider an Ising model consisting of $N$ spin variables $\mathbf{s} = (s_i = \pm 1)_{i=1}^N$ and obeying the following distribution:

$$P_{\text{Ising}}(\mathbf{s} | \mathbf{J}, \mathbf{H}) = \frac{1}{Z_{\text{Ising}}} e^{\sum_{i<j} J_{ij} s_i s_j + \sum_{i=1}^N H_i s_i}, \quad (1)$$

where $\mathbf{J} \in \mathbb{R}^{N \times N}$ and $\mathbf{H} \in \mathbb{R}^N$ are termed couplings and external fields, respectively. The setting of the inverse Ising problem is aimed at inferring the couplings and external fields from a given dataset of spin snapshots (samples) $D^M \equiv \{\mathbf{s}^{(\mu)}\}_{\mu=1}^M$, where $M$ denotes the dataset size. The representative statistical solution or estimator is the ML estimator defined by

$$\{\hat{J}^{\text{ML}}(D^M), \hat{H}^{\text{ML}}(D^M)\} = \arg \min_{\mathbf{J}, \mathbf{H}} \left\{ -\sum_{\mu=1}^{M} \log P_{\text{Ising}}(\mathbf{s}^{(\mu)} | \mathbf{J}, \mathbf{H}) \right\}. \quad (2)$$

Hereafter, the symbol $\hat{\cdot}$ is used to represent an estimator. This canonical estimator has some desirable properties such as consistency, unbiasedness, and efficiency when the dataset size $M$ is large enough while $N$ is kept finite (the asymptotic limit) [27, 28]. Consistency means that the estimator converges in probability to the true parameter in the large $M$ limit; unbiasedness means that the mean of the estimator (which is a random variable when the dataset is randomly generated) becomes identical to the true parameter; unbiased estimators yielding the lowest mean-squared error are termed efficient. Note that the ML estimator enjoys the unbiasedness and the efficiency only in the asymptotic sense, and for finite size datasets it can have bias and be not efficient. In this sense, the ML estimator is good, but is not always appropriate for the inverse Ising framework for some computational reasons. The PL method [11] overcomes this problem by replacing the likelihood with the conditional distribution $P(s_i | \mathbf{s}_{\setminus i}, \mathbf{J}_i, H_i)$ for each $s_i$, where $\mathbf{J}_i = (J_{ij})_j$ is the coupling vector connected to the
ith spin, and $\mathbf{s}_{\backslash i}$ denotes the spin vector without the ith component. The explicit form is

$$P(s_i | \mathbf{s}_{\backslash i}, \mathbf{J}_i, H_i) = \frac{1}{Z_i} e^{s_i (\sum_{j(\neq i)} J_{ij} s_j + H_i)}, \quad (3)$$

$$Z_i = 2 \cosh \left( \sum_{j(\neq i)} J_{ij} s_j + H_i \right). \qquad (4)$$

The PL estimator is obtained separately for each i by

$$\begin{aligned} \{\hat{\mathbf{J}}_i^{\text{PL}}(D^M), \hat{H}_i^{\text{PL}}(D^M)\} &= \arg \min_{\mathbf{J}_i, H_i} \left\{ -\sum_{\mu=1}^{M} \log P \left( s_i^{(\mu)} | \mathbf{s}_{\backslash i}^{(\mu)}, \mathbf{J}_i, H_i \right) \right\} \\ &= \arg \min_{\mathbf{J}_i, H_i} \left\{ \sum_{\mu=1}^{M} \ell^{\text{PL}} \left( s_i^{(\mu)} h_i(\mathbf{s}_{\backslash i}^{(\mu)}, \mathbf{J}_i, H_i) \right) \right\}, \end{aligned} \quad (5)$$

where

$$h_i(\mathbf{s}_{\backslash i}, \mathbf{J}_i, H_i) = \sum_{j(\neq i)} J_{ij} s_j + H_i, \qquad (6)$$

$$\ell^{\text{PL}}(x) = -x + \log 2 \cosh x. \qquad (7)$$

Two remarkable properties are held by the PL estimator. The first is its consistency. When the dataset size $M$ is sufficiently large, the PL estimator converges to the true values $\{\hat{\mathbf{J}}_i^{\text{PL}}, \hat{H}_i^{\text{PL}}\} \to \{\mathbf{J}_i, H_i\}$. The second is its locality. Owing to the factorized nature of PL, each coupling vector $\mathbf{J}_i$ can be assessed independently with low (polynomial) computational cost, which is in contrast to the ML estimator for which the exponentially large computational cost is necessary when assessing the partition function. However, coupling symmetry $J_{ij} = J_{ji}$ is inevitably lost ($\hat{J}_{ij}^{\text{PL}} \neq \hat{J}_{ji}^{\text{PL}}$ in general).

This local property of PL also simplifies the theoretical treatment, which inspired explorations of the ‘optimal cost function’ in a class of models with the same locality as the PL case [24, 26]. On the basis of these studies, we treat this local learning class and assume that the cost function argument is the only product of the local spin and effective field. The cost function and corresponding estimator are denoted as $\ell = \ell(s_i h(\mathbf{s}_{\backslash i}))$ and $\{\hat{\mathbf{J}}_i, \hat{H}_i\}$, respectively. When specifying a certain model in the class, the appropriate superscript will be attached as equations (5) and (7).

We have an additional remark on the lost symmetry in the PL estimator. This limitation of the PL estimator can be recovered by adding PL cost functions of different sites i and j and by performing the minimization with respect to a symmetric variable $J_{ij}$ shared in both terms [29]. Although our theoretical analysis presented in this paper can also be generalized to this case with a slight increase in the respective computational steps shown in section 3, we restrict ourselves on the simplest case without such a symmetrizing operation.
### 2.1.1. Teacher-student scenario.
Here, we consider the inverse Ising problem in the teacher-student scenario. Specifically, the dataset $D^M \equiv \{\mathbf{s}^{(\mu)}\}_{\mu=1}^M$ is assumed to be independently identically distributed (i.i.d.) samples from a teacher Ising model with the couplings $J^*$ and external fields $\mathbf{H}^*$, and a student Ising model attempts to infer the teacher couplings and fields from the dataset. The inference accuracy is quantified by the residual sum of squares (RSS) between the teacher couplings and student's estimator

$$ \mathcal{E} = \|\mathbf{J}_i^* - \hat{\mathbf{J}}_i\|_2^2 = \sum_{j(\neq i)} (\mathbf{J}_{ij}^* - \hat{\mathbf{J}}_{ij})^2 = R^* - 2\rho + R, \quad (8) $$

where we defined the following three macroscopic parameters:

$$ R^* = \sum_{j(\neq i)} (\mathbf{J}_{ij}^*)^2, \quad (9a) $$

$$ R = \sum_{j(\neq i)} (\hat{\mathbf{J}}_{ij})^2, \quad (9b) $$

$$ \rho = \sum_{j(\neq i)} \mathbf{J}_{ij}^* \hat{\mathbf{J}}_{ij}. \quad (9c) $$

These are later associated with the order parameters in the statistical mechanical analysis.

## 2.2. Statistical mechanical analysis

Here, we explain the statistical mechanical formulation developed in [24] to analyze the theoretical performance of local learning models. For simplicity of the theoretical analysis, we assume the absence of external fields $\mathbf{H} = \mathbf{0}$ both in the teacher and student models.

### 2.2.1. General framework.
The basic idea of statistical mechanical analysis is to introduce the following Hamiltonian and Boltzmann distribution induced by the cost function $\ell$:

$$ \mathcal{H}(\mathbf{J}|\mathbf{D}^M) = \sum_{\mu=1}^{M} \ell \left( s_{0}^{(\mu)} h(\mathbf{s}_{0}^{(\mu)}, \mathbf{J}) \right), \quad (10) $$

$$ P(\mathbf{J}|\mathbf{D}^M) = \frac{1}{Z} e^{-\beta \mathcal{H}(\mathbf{J}|\mathbf{D}^M)}, \quad (11) $$

where

$$ h(\mathbf{s}_{0}, \mathbf{J}) = \sum_{j=1}^{N-1} J_j s_j, \quad (12) $$

$$ Z = \mathrm{Tr} e^{-\beta \mathcal{H}(\mathbf{J}|\mathbf{D}^M)}. \quad (13) $$
Here, we reorder the spins and focus on the zeroth spin and its coupling vector **J**. The external field is set to zero as declared above; Tr<sub>**J**</sub> denotes the integration with respect to **J** with an appropriate measure, the detailed definition of which is given in appendix A. In the limit $\beta \to \infty$, the Boltzmann distribution converges to a pointwise measure on the estimator $\hat{\mathbf{J}} = \arg \min_{\mathbf{J}} \{\sum_{\mu=1}^{M} \ell(s_0^{(\mu)} h(\mathbf{s}^{(\mu)}, \mathbf{J}))\}$; thus, the Boltzmann distribution enables the capture of the estimator's properties. We are interested in the thermodynamic limit $N \to \infty$ while keeping $\alpha = M/N = O(1)$; hence, the free energy density $f = -(N\beta)^{-1}\log Z$ will show the self-averaging property. Thus, for typical datasets, the free energy density converges to its average in this limit. The averaged free energy density is denoted as

$$f^{\text{ave}} = -(N\beta)^{-1}[\log Z]_{D^M}, \quad (14)$$

where the square brackets $[\cdot]_{D^M}$ denote the average over the dataset, which is the average over the teacher Ising model:

$$[\cdot]_{D^M} = \sum_{\mathbf{s}^{(1)}, \dots, \mathbf{s}^{(M)}} (\cdot) \prod_{\mu=1}^{M} P_{\text{Ising}}(\mathbf{s}^{(\mu)} | \mathbf{J}^*). \quad (15)$$

The average of $\log Z$ over $D^M$ is, however, generally difficult to compute. The replica method is a prescription to overcome this difficulty and is symbolized by the following identity:

$$f^{\text{ave}} = -(N\beta)^{-1}[\log Z]_{D^M} = \lim_{n \to 0} - (nN\beta)^{-1} \log [Z^n]_{D^M}. \quad (16)$$

Assuming that $n$ is a positive integer, we can rewrite $[Z^n]_{D^M}$ as

$$[Z^n]_{D^M} = \text{Tr}_{\{\mathbf{J}^n\}_{a=1}^n} \left\{ \sum_s P_{\text{Ising}}(\mathbf{s} | \mathbf{J}^*) e^{-\beta \sum_{a=1}^n \ell(s_0 h(s_0, \mathbf{J}^a))} \right\}^M. \quad (17)$$

Rewriting this by introducing variables $(h^a = \sum_{j=1}^{N-1} J_j^a s_j)_a$ and $h^* = \sum_{j=1}^{N-1} J_j^* s_j$, which are hereafter called cavity fields, we get

$$\begin{aligned} [Z^n]_{D^M} &= \text{Tr}_{\{\mathbf{J}^n\}_{a=1}^n} \left\{ \sum_s \int dh^* \prod_{a=1}^n dh^a \delta\left(h^* - \sum_{j=1}^{N-1} J_j^* s_j\right) \right. \\ & \qquad \left. \times \prod_{a=1}^n \delta\left(h^a - \sum_{j=1}^{N-1} J_j^a s_j\right) P_{\text{Ising}}(\mathbf{s} | \mathbf{J}^*) e^{-\beta \sum_{a=1}^n \ell(s_0 h^a)} \right\}^M \\ &= \text{Tr}_{\{\mathbf{J}^n\}_{a=1}^n} \left\{ \sum_{s_0} \int dh^* \prod_{a=1}^n dh^a P_{\text{cav}}(h^*, \{h^a\}_{a=1}^n | \mathbf{J}^*, \{\mathbf{J}^a\}_{a=1}^n) \frac{1}{Z_0} e^{s_0 h^*} e^{-\beta \sum_{a=1}^n \ell(s_0 h^a)} \right\}^M, \end{aligned} \quad (18)$$
At the last equality, we took the summation over the spins except for $s_0$, yielding the joint distribution $P_{cav}(h^*, \{h^a\}_{a=1}^n | \mathbf{J}^*, \{\mathbf{J}^a\}_{a=1}^n)$ of the cavity fields. The normalization constant $Z_0$ is defined through its marginal as

$$Z_0 = \int dh^* P_{cav}(h^* | \mathbf{J}^*) 2 \cosh h^*. \qquad (19)$$

To proceed further with the computation, we need to specify the functional form of the cavity field distribution and take the average over it. This is possible when the teacher is a fully-connected model, as demonstrated in [24]. We review this result below as it contains some steps essential for the sparsely-connected case.

### 2.2.2. Revisit the fully-connected case.
When the teacher is a fully-connected model, we can use the central limit theorem and assume that the cavity fields are multivariate Gaussian variables with appropriate covariances and means. In [24], the authors assumed that the teacher system is in the paramagnetic phase and the replica symmetry (RS) holds in both the student and teacher systems. These assumptions imply that the following four order parameters are sufficient to describe the free energy density:

$$Q^* = \sum_{i,j} C_{ij}^{\wedge 0} J_i^* J_j^*, \qquad (20a)$$

$$Q = \sum_{i,j} C_{ij}^{\wedge 0} J_i^a J_j^a, \qquad (20b)$$

$$q = \sum_{i,j} C_{ij}^{\wedge 0} J_i^a J_j^b, (a \neq b), \qquad (20c)$$

$$m = \sum_{i,j} C_{ij}^{\wedge 0} J_i^* J_j^a, \qquad (20d)$$

where $C^{\wedge 0}$ is the correlation matrix between the spins:

$$C_{ij}^{\wedge 0} = \langle s_i s_j \rangle^{\wedge 0} - \langle s_i \rangle^{\wedge 0} \langle s_j \rangle^{\wedge 0} = \langle s_i s_j \rangle^{\wedge 0}, \qquad (21)$$

where $\langle \cdots \rangle^{\wedge 0}$ denotes the average over the teacher Ising model without the zeroth spin; the last equality is due to the paramagnetic assumption. The covariances of the cavity fields are described by equation (20) as $\langle h^a h^b \rangle^{\wedge 0} = Q\delta_{ab} + (1 - \delta_{ab})q$, $\langle h^* h^a \rangle^{\wedge 0} = m$, $\langle (h^*)^2 \rangle^{\wedge 0} = Q^*$. Upon assuming this, we can rewrite $[Z^n]_{DM}$ as

$$[Z^n]_{DM} = \int dQ dq dme^{NS(C^{\wedge 0}, \mathbf{J}^*, Q, q, m) + M \log L(Q^*, Q, q, m)}, \qquad (22)$$

where
$$e^{NS(C^{\setminus 0}, \mathbf{J}^*, Q, q, m)} \equiv \mathrm{Tr}_{\{\mathbf{J}^a\}_{a=1}^n} \prod_{a=1}^n \left\{ \delta\left(Q - \sum_{i,j} C_{ij}^{\setminus 0} J_i^a J_j^a\right) \delta\left(m - \sum_{i,j} C_{ij}^{\setminus 0} J_i^* J_j^a\right) \right\} \\ \times \prod_{a<b} \delta\left(q - \sum_{i,j} C_{ij}^{\setminus 0} J_i^a J_j^b\right), \qquad (23)$$

$$L(Q^*, Q, q, m) \equiv \sum_{s_0} \int dh^* \prod_{a=1}^{n} dh^a P_{\mathrm{cav}}(h^*, \{h^a\}_{a=1}^n | Q^*, Q, q, m) \frac{1}{Z_0} e^{s_0 h^*} e^{-\beta \sum_{a=1}^n \ell(s_0 h^a)}. \qquad (24)$$

Deferring the detailed computations to appendix A, we immediately have the result in the limit $n \to 0$:

$$\lim_{n \to 0} \frac{1}{n} S(C^{\setminus 0}, \mathbf{J}^*, Q, q, m) = \frac{1}{2} \left\{ \frac{Q - m^2/Q^*}{Q - q} + \log 2\pi + \log(Q - q) - \frac{1}{N} \mathrm{Tr} \log C^{\setminus 0} \right\}, \qquad (25)$$

$$\lim_{n \to 0} \frac{1}{n} \log L(Q^*, Q, q, m) = \int Dz e^{\sqrt{\frac{m^2}{q}z} - \frac{1}{2}\frac{m^2}{q}} \log \int Dv e^{-\beta \ell(\sqrt{Q-qv} + \sqrt{qz})}, \qquad (26)$$

where we introduce a Gaussian measure

$$\int Dx(\cdots) \equiv \int_{-\infty}^{\infty} \frac{dx}{\sqrt{2\pi}} e^{-\frac{1}{2}x^2(\cdots)}. \qquad (27)$$

Further, we take the limit $\beta \to \infty$, which requires the following relation:

$$\lim_{\beta \to \infty} \beta(Q - q) = \chi = O(1). \qquad (28)$$

After straightforward calculations, we get

$$f^{\mathrm{ave}}(\beta \to \infty) = \lim_{\beta \to \infty} -\frac{1}{\beta} \left( \lim_{n \to 0} \frac{1}{n} S(C^{\setminus 0}, \mathbf{J}^*, Q, q, m) + \alpha \lim_{n \to 0} \frac{1}{n} \log L(Q^*, Q, q, m) \right) \\ = -\mathrm{Extr}_{Q, \chi, m} \left\{ \frac{1}{2} \frac{Q - m^2/Q^*}{\chi} + \alpha \int Dz \max_y \left( -\frac{(y - \sqrt{Q}z - m)^2}{2\chi} - \ell(y) \right) \right\}, \qquad (29)$$

where $\alpha = M/N$, and $\mathrm{Extr}_x$ denotes the extremization condition with respect to $x$. The extremization condition yields the following equations of state (EOS):
$$0 = \frac{1}{\chi} - \frac{\alpha}{\sqrt{Q}} \int Dz \frac{\partial \ell}{\partial y} \Big|_{y=\hat{y}}, \qquad (30a)$$

$$0 = -\frac{m}{Q^*\chi} - \alpha \int Dz \frac{\partial \ell}{\partial y} \Big|_{y=\hat{y}}, \qquad (30b)$$

$$0 = -\frac{1}{\chi^2} \left( Q - \frac{m^2}{Q^*} \right) + \alpha \int Dz \left( \frac{\partial \ell}{\partial y} \Big|_{y=\hat{y}} \right)^2, \qquad (30c)$$

where

$$\hat{y}(z, Q, \chi, m) = \arg \max_y \left( -\frac{(y - \sqrt{Q}z - m)^2}{2\chi} - \ell(y) \right). \qquad (31)$$

Using the solution of equation (30), we can evaluate the macroscopic parameters (9) by

$$R = \left( Q - \frac{m^2}{Q^*} \right) \frac{1}{N} \mathrm{Tr} (C^{\wedge 0})^{-1} + R^* \left( \frac{m}{Q^*} \right)^2, \qquad (32a)$$

$$\rho = R^* \frac{m}{Q^*}. \qquad (32b)$$

These relations can be derived by a standard technique using auxiliary variables [24] and the derivation is given in appendix B. Once the values of $R$ and $\rho$ are obtained, the RSS is eventually computed by equation (8). Note that the values of $Q^*$ and $R^*$ are directly determined by the problem setting, as well as the inverse correlation function $(C^{\wedge 0})^{-1}$. To obtain these quantities, we need to separately solve the direct problem.

### 3. Details of the sparsely-connected case

This section provides the extension of the above result to the sparsely-connected case, which is the main contribution of this paper. To this end, we introduce an ansatz about the estimator's behavior as well as the functional form of the cavity field distribution. Under the ansatz, the cavity field is decomposed into a signal and a noise, and it is shown that the noise part obeys essentially the same EOS as the fully-connected case. To complete the computation under the ansatz, the tree-like structure of the coupling network of the teacher model is employed.

#### 3.1. Ansatz for the sparse case

In contrast to the fully-connected case, the cavity field distribution in the sparse case cannot be regarded as Gaussian; the distribution of $h^*$ actually becomes the sum of a few pointwise measures, which is far from Gaussian. Hence, we need a new ansatz to handle the cavity field distribution in the sparse case.
To obtain an idea of how to resolve this, let us consider an ideal situation where we know which couplings are nonzero. Let us suppose that the zeroth spin is connected to $c = O(1)$ neighboring spins, and introduce two sets of indices $\Omega = \{i|J_i^* \neq 0, i \in \{1, \dots, N-1\}\}$ and $\bar{\Omega} = \{i|J_i^* = 0, i \in \{1, \dots, N-1\}\}$, where $\Omega(\bar{\Omega})$ is called the active (inactive) set; $|\Omega| = c$ and $|\bar{\Omega}| = N-1-c$. If we know $\Omega$ and $\bar{\Omega}$ in advance, then the inference should be conducted only on $\{J_i|i \in \Omega\}$. Accordingly, the number of variables to be inferred is just $c = O(1)$; hence, the dataset size $M = O(N)$ is sufficiently large. Thus, we can apply the asymptotic theory of statistics, which implies that an estimator in this ideal case behaves as

$$ \hat{J}_i^{\text{oracle}} = \begin{cases} J_i^* + \Delta_i & (i \in \Omega) \\ 0 & (i \in \bar{\Omega}) \end{cases} \qquad (33) $$

We call this an *oracle* estimator. The ‘error’ from the true solution, $\Delta_i$, is a random variable. In the local learning class with appropriate (consistent) cost functions such as PL [30], $\Delta_i$ is considered to be zero mean with variance decreasing at the rate of $O(M^{-1}) = O(N^{-1})$. The RSS is written as $\mathcal{E} = \sum_{i \in \Omega} \Delta_i^2 = O(N^{-1})$, and vanishes in the thermodynamic limit.

Based on these observations about the oracle estimator, we assume that the (non-oracle) estimator obtained from consistent cost functions obeys the following form:

$$ \hat{J}_i \doteq \begin{cases} \bar{J}_i + \Delta_i & (i \in \Omega) \\ \Delta_i & (i \in \bar{\Omega}) \end{cases} \qquad (34) $$

where we again assume that $\Delta_i$ is a random variable which is asymptotically zero mean with variance scaled as $O(N^{-1})$; the correlations among $\{\Delta_i\}_i$ are also assumed to be sufficiently weak. The quantity $\bar{J}_i$ is interpreted as the mean value of the estimator and will deviate from the true value $J_i^*$ owing to the extensive number of noise terms $\{\Delta_i\}_i$. The values of $\{\bar{J}_i\}_{i \in \bar{\Omega}}$ are later computed by taking the minimization condition of the free energy as the order parameters. The applicable range of this ansatz is discussed in section 3.3.

Let us examine the consequence of the ansatz. The RSS can be written as

$$ \mathcal{E} \approx \sum_{i \in \Omega} (J_i^* - \bar{J}_i)^2 + \sum_{i \in \bar{\Omega}} \Delta_i^2 \qquad (35) $$

Now, there are two non-negligible contributions to the RSS coming from the bias in $\Omega$ and the noise in $\bar{\Omega}$; the RSS remains finite even in the limit $N \to \infty$ in contrast to the ideal case. The ansatz also allows us to decompose the cavity field as

$$ h^a = h_{\Omega} + h_{\Delta}^a, \qquad (36) $$

$$ h_{\Omega} = \sum_{j \in \Omega} \bar{J}_j s_j, \quad h_{\Delta}^a = \sum_j \Delta_j^a s_j \approx \sum_{j \in \bar{\Omega}} \Delta_j^a s_j, \qquad (37) $$

where $h_{\Delta}^a$ is termed as the ‘noise’ part. Furthermore, we can assume that $h_{\Omega}$ and $h_{\Delta}^a$ are asymptotically independent in the limit $N \to \infty$. This assumption is reasonable, and a schematic to explain this is given in figure 1. Owing to the tree-like structure, we
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic explaining the independence between ($s_0$, $\mathbf{s}_\Omega$) and $h_\Delta^a$. The majority in the sum in $h_\Delta^a$ is uncorrelated with spins in $\Omega$ and dominates $h_\Delta^a$ in the thermodynamic limit.</figcaption></figure>

can define the generation $g$ of a spin $s$ from $\Omega$ as the shortest path length between $s$ and any spin in $\Omega$ along the network. As $g$ grows, the correlation with $\{s_i | i \in \Omega\}$ decays exponentially fast, while the number of spins belonging to generation $g$ exponentially increases. If the correlation decay is sufficiently faster than the increase of the spins, then the majority of spins in the network can be regarded as uncorrelated with $\Omega$. Some terms in $h_\Delta^a$ are certainly correlated with $h^*$, but their contribution is $O(\sqrt{N}^{-1})$ because $\Delta_i = O(\sqrt{N}^{-1})$ and the number of correlating terms is $O(1)$ owing to the sufficiently fast decay of correlations. Hence, the contribution of correlating terms vanishes and the uncorrelated majority with $\Omega$ completely dominates $h_\Delta^a$ in the thermodynamic limit. These observations indicate that equation (17) can be now decomposed as follows:

$$ [Z^n]_{DM} = \mathrm{Tr}\{\mathcal{J}^a\}_{a=1}^n \left\{ \sum_s \int P_{\mathrm{Ising}}(\mathbf{s} | \mathbf{J}^*) e^{-\beta \sum_{a=1}^n \ell(s_0 h^a)} \right\}^M \\ \approx \mathrm{Tr}\{\Delta_j^a\}_{a=1}^n \left\{ \sum_s \int \prod_{a=1}^n dh_\Delta^a P_{\mathrm{Ising}}(\mathbf{s} | \mathbf{J}^*) \delta\left(h_\Delta^a - \sum_{j \in \Omega} \Delta_j^a s_j\right) e^{-\beta \sum_{a=1}^n \ell(s_0(\sum_{j \in \Omega} \bar{J}_j s_j + h_\Delta^a))} \right\}^M \\ = \mathrm{Tr}\{\Delta_j^a\}_{a=1}^n \left\{ \sum_{s_0, \mathbf{s}_\Omega} \int \prod_{a=1}^n dh_\Delta^a P(s_0, \mathbf{s}_\Omega, \{h_\Delta^a\}_a | \mathbf{J}^*, \{\Delta_j^a\}_a) e^{-\beta \sum_{a=1}^n \ell(s_0(\sum_{j \in \Omega} \bar{J}_j s_j + h_\Delta^a))} \right\}^M \\ \approx \mathrm{Tr}\{\Delta_j^a\}_{a=1}^n \left\{ \sum_{s_0, \mathbf{s}_\Omega} \int \prod_{a=1}^n dh_\Delta^a P_{\mathrm{Ising}}(s_0, \mathbf{s}_\Omega | \mathbf{J}^*) P_{\mathrm{Cav}}(\{h_\Delta^a\}_a | \{\Delta_j^a\}_a) e^{-\beta \sum_{a=1}^n \ell(s_0(\sum_{j \in \Omega} \bar{J}_j s_j + h_\Delta^a))} \right\}^M . \quad (38) $$

In the second line, we performed the variable transformation $\Delta^a = \mathbf{J}^u - \bar{\mathbf{J}}$ and neglected the contribution $\sum_{j \in \Omega} \Delta_j^a s_j$ in $h_\Delta^a$ as equation (37). In the third line, we denoted $\mathbf{s}_\Omega = \{s_i | i \in \Omega\}$ and $\bar{\mathbf{s}}_\Omega = \{s_i | i \in \bar{\Omega}\}$, and performed the sum over $\mathbf{s}_\Omega$, yielding the joint distribution $P(s_0, \mathbf{s}_\Omega, \{h_\Delta^a\}_a | \mathbf{J}^*, \{\Delta_j^a\}_a)$. In the fourth line, we used the asymptotic uncorrelatedness between $h_\Delta^a$ and $(s_0, \mathbf{s}_\Omega)$ discussed so far.
Now, the central limit theorem can be applied to the noise part {$h_{\Delta}^{a}$}$_{a}$, and they can be regarded as Gaussian variables. As the fully-connected case, two order parameters describing their covariances are introduced:

$$
Q \equiv \sum_{i,j} C_{ij}^{\setminus 0} \Delta_i^a \Delta_j^a, \qquad (39)
$$

$$
q \equiv \sum_{i,j} C_{ij}^{\setminus 0} \Delta_i^a \Delta_j^b, (a \neq b). \qquad (40)
$$

Counterparts of $m$ and $Q^*$ are unnecessary because the dependence on $(s_0, \mathbf{s}_\Omega)$ is
separately and explicitly treated in the present formulation. Then,

$$
[Z^n]_{DM} \approx \int dQ dq e^{NS(C^{\setminus 0}, Q, q) + M \log L(J^*, \bar{J}, Q, q)}, \quad (41)
$$

where

$$
e^{NS(C^{\setminus 0}, Q, q)} \equiv \mathrm{Tr}_{\{\Delta^a\}_{a=1}^n} \prod_{a=1}^n \delta \left( Q - \sum_{i,j} C_{ij}^{\setminus 0} \Delta_i^a \Delta_j^a \right) \prod_{a<b} \delta \left( q - \sum_{i,j} C_{ij}^{\setminus 0} \Delta_i^a \Delta_j^b \right), \quad (42)
$$

$$
L(\mathbf{J}^*, \bar{\mathbf{J}}, Q, q) = \sum_{s_0, s_\Omega} \int \prod_{a=1}^{n} dh_{\Delta}^{a} P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega} | \mathbf{J}^*) P_{\text{cav}}(\{h_{\Delta}^{a}\}_{a=1}^{n} | Q, q) e^{-\beta \sum_{a=1}^{n} \ell(s_0(\sum_{j \in \Omega} \bar{J}_j s_j + h_{\Delta}^{a}))}. \quad (43)
$$

Again, using the techniques in appendix A we get

$$
\lim_{n \to 0} \frac{1}{n} S(C^{\setminus 0}, Q, q) = \frac{1}{2} \left\{ \frac{Q}{Q-q} + \log 2\pi + \log(Q-q) - \frac{1}{N} \mathrm{Tr} \log C^{\setminus 0} \right\}, \quad (44)
$$

$$
\lim_{n \to 0} \frac{1}{n} \log L(\mathbf{J}^*, \bar{\mathbf{J}}, Q, q) = \sum_{s_0, s_\Omega} P_{\text{Ising}}(s_0, \mathbf{s}_\Omega | \mathbf{J}^*) \int Dz \log \int Dv e^{-\beta \ell(s_0(\sum_{j \in \Omega} \bar{J}_j s_j + \sqrt{Q-q}v + \sqrt{q}z))}. \quad (45)
$$

Employing the relation (28) and taking the $\beta \to \infty$ limit, we get
$$f^{\text{ave}} = -\text{Extr}_{Q,\chi} \left\{ \frac{1}{2\chi} + \alpha \sum_{s_0, \mathbf{s}_{\Omega}} P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega} | \mathbf{J}^*) \int Dz \max_y \right. \\ \left. \times \left( -\frac{(y - s_0(\sqrt{Q}z + \sum_{j \in \Omega} \bar{J}_j s_j))^2}{2\chi} - \ell(y) \right) \right\}. \quad (46)$$

The extremization condition with respect to $Q$ and $\chi$ gives

$$0 = \frac{1}{\chi} - \frac{\alpha}{\sqrt{Q}} \sum_{s_0, \mathbf{s}_{\Omega}} P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega} | \mathbf{J}^*) s_0 \int Dz \frac{\partial \ell}{\partial y} \Big|_{y=\hat{y}}, \quad (47a)$$

$$0 = -\frac{Q}{\chi^2} + \alpha \sum_{s_0, \mathbf{s}_{\Omega}} P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega} | \mathbf{J}^*) \int Dz \left( \left. \frac{\partial \ell}{\partial y} \right|_{y=\hat{y}} \right)^2, \quad (47b)$$

where

$$\hat{y}(z, s_0, \mathbf{s}_{\Omega} | Q, \chi, \{\bar{J}_j\}_{j \in \Omega}) = \arg \max_y \left( -\frac{(y - s_0 (\sqrt{Q}z + \sum_{j \in \Omega} \bar{J}_j s_j))^2}{2\chi} - \ell(y) \right). \quad (48)$$

Further, the mean estimates $\{\bar{J}_j\}_{j \in \Omega}$ are also evaluated by the extremization condition. The result for $\bar{J}_j$ is given by

$$0 = \sum_{s_0, \mathbf{s}_{\Omega}} P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega} | \mathbf{J}^*) \int Dz \left. \frac{\partial \ell}{\partial y} \right|_{y=\hat{y}} s_0 s_j. \quad (49)$$

Using the parameters $Q, \chi, \{\bar{J}_i\}_{i \in \Omega}$ computed by equations (47)–(49), we can evaluate the RSS which is expressed in the present setting as

$$\mathcal{E} \approx \sum_{i \in \Omega} (J_i^* - \bar{J}_i)^2 + \sum_{i \in \bar{\Omega}} \Delta_i^2 = \sum_{i \in \Omega} (J_i^* - \bar{J}_i)^2 + \frac{Q}{N} \text{Tr}(C^{\langle 0 \rangle})^{-1}. \quad (50)$$

The quantity $(C^{\langle 0 \rangle})^{-1}$ will be computed in another discussion on the direct problem as the fully-connected case. In addition, $P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega} | \mathbf{J}^*)$ will also be assessed separately in the sparse case. These points are addressed in the next subsection.

## 3.2. Direct problem's properties

The inverse problem essentially requires certain information from its direct problem counterpart. Necessary information to compute the quantities of interest depends on the system's properties. In the fully-connected case, two-body quantities such as $(C^{\langle 0 \rangle})^{-1}$ and $\sum_{i,j} C_{ij}^{\langle 0 \rangle} J_i^* J_j^*$ are sufficient. However, in the sparse case, higher-order information is needed because the central limit theorem does not fully dominate the system. Hence, the functional form of $P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega})$ becomes necessary, as seen in equation (46). Techniques for computing such quantities in the sparse case largely advanced in the '90-'00s. Here, we quote a portion of the results to compute the necessary quantities. For readers
interested in the detailed techniques, please refer to [31, 32]. Although these techniques are applied in general situations, to obtain compact analytic forms of the quantities of interest, we rely on the assumptions that the teacher model is in the paramagnetic phase and the external fields are absent.

### 3.2.1. Marginal distribution of the teacher model.
The marginal distribution $P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega}|\mathbf{J}^*)$ is computed by marginalizing the whole distribution $P_{\text{Ising}}(\mathbf{s}|\mathbf{J}^*)$ with respect to $\mathbf{s}_{\bar{\Omega}}$. In general, this operation requires nontrivial computations and the resultant distribution becomes dependent on parameters among the marginalized spins. However, under the present assumptions, such dependencies do not exist and the expression becomes rather simple:

$$P_{\text{Ising}}(s_0, \mathbf{s}_{\Omega}|\mathbf{J}^*) = \frac{1}{Z_{\Omega_c}} e^{s_0 \sum_{j \in \Omega} J_j^* s_j}, Z_{\Omega_c} = \sum_{s_0, \mathbf{s}_{\Omega}} e^{s_0 \sum_{j \in \Omega} J_j^* s_j}, \quad (51)$$

where $\Omega_c$ denotes the union index set of 0 and $\Omega$. This form is applied to equations (47) and (49) to obtain the order parameters.

### 3.2.2. Inverse correlation function.
Next, we compute the inverse correlation function $C^{-1}$; hereafter, we treat the whole system and discard the superscript $\backslash 0$ as it is not essential.

On the tree-like networks, the so-called *Bethe approximation* gives an exact result as the system size tends to infinity, and it provides a compact analytic form of the inverse correlation function. The derivation is given in [33, 34] and we here just quote the result. Under the present assumptions without magnetizations and external fields, it becomes

$$(C^{-1})_{ij} = \left( \sum_{k \in \partial i} \frac{1}{1 - \tanh^2(J_{ik})} - (c_i - 1) \right) \delta_{ij} - \frac{\tanh(J_{ij})}{1 - \tanh^2(J_{ij})} (1 - \delta_{ij}), \quad (52)$$

where $\partial i$ denotes the index set of nodes connected to $i$ and $c_i$ denotes the connectivity or the number of edges connecting to node $i$.

### 3.3. Applicable range of the ansatz

Here we consider the applicable range of the ansatz (34). This ansatz is a strong statement since it allows us to not think of any possible biases of the estimators outside the active set $\Omega$. When is this valid? How does it relate to the tree-like network structure?

To obtain answers to these questions, we rethink equation (49). An important observation is that this equation is merely the zero-gradient condition of $\ell$ with respect to $J_j$, ($j \in \Omega$) averaged over the datasets, as shown below. Denoting the empirical average on the dataset $D^M$ by $\overline{\cdots}^{D^M}$, and using the statistical mechanical analysis explained so far, we can write the zero-gradient condition with respect to $J_j$ for $j \in \Omega$ as

$$0 = \frac{\overline{\partial \ell(s_0 h(\mathbf{s}_{\backslash 0}, \mathbf{J})))}}{\partial J_j} \Big|_{J=\tilde{J}(D^M)} = \frac{\overline{\partial \ell(y)}}{\partial y} s_0 s_j, \quad (53)$$
where we put $y = s_0h(s_{\backslash 0}, \hat{\mathbf{J}}) = s_0(\sum_i \hat{J}_i s_i)$. With this expression, we replace the estimator $\hat{\mathbf{J}}$ with the average over equation (11), take the average over the dataset, and use the replica method. The result is

$$0 = \lim_{\beta \to \infty} \lim_{n \to 0} \mathrm{Tr}_{\{\mathbf{J}^a\}_{a=1}^n} \left( \sum_s P_{\mathrm{Ising}}(\mathbf{s} | \mathbf{J}^*) e^{-\beta \sum_{a=1}^n \ell(y^a)} \right)^{M-1} \\ \times \left( \sum_s P_{\mathrm{Ising}}(\mathbf{s} | \mathbf{J}^*) e^{-\beta \sum_{a=1}^n \ell(y^a)} \frac{\partial \ell(y^1)}{\partial y^1} s_0 s_j \right). \quad (54)$$

where $y^a \equiv s_0(\sum_i J_i^a s_j)$. The ansatz (34) and RS used in section 3.1, in short, say

$$y^a \stackrel{\text{ansatz}}{=} s_0 \left\{ \sum_{i \in \Omega} \bar{J}_i s_i + h_{\Delta}^a \right\} \stackrel{\text{RS}}{=} s_0 \left\{ \sum_{i \in \Omega} \bar{J}_i s_i + \sqrt{Q - q} v^a + \sqrt{q} z \right\}. \quad (55)$$

where $v^a, z \sim \mathcal{N}(0, 1)$. Applying this form and following the same line of computations as in section 3.1, we get

$$\left( \sum_s P_{\mathrm{Ising}}(\mathbf{s} | \mathbf{J}^*) e^{-\beta \sum_{a=1}^n \ell(y^a)} \frac{\partial \ell(y^1)}{\partial y^1} s_0 s_j \right) \\ \xrightarrow{N \to \infty} \sum_{s_0, s_{\Omega}} P_{\mathrm{Ising}}(s_0, s_{\Omega} | \mathbf{J}^*) \int Dz \left( \int Dv e^{-\beta \ell(y(z,v))} \right)^n \frac{\int Dv^1 e^{-\beta \ell(y(z,v^1))} \frac{\partial \ell(y^1)}{\partial y^1} s_0 s_j}{\int Dv^1 e^{-\beta \ell(y(z,v^1))}} \\ \xrightarrow{n \to 0} \sum_{s_0, s_{\Omega}} P_{\mathrm{Ising}}(s_0, s_{\Omega} | \mathbf{J}^*) \int Dz \frac{\int Dv e^{-\beta \ell(y(z,v))} \frac{\partial \ell(y^1)}{\partial y^1} s_0 s_j}{\int Dv e^{-\beta \ell(y(z,v))}} \\ \xrightarrow{\beta \to \infty} \sum_{s_0, s_{\Omega}} P_{\mathrm{Ising}}(s_0, s_{\Omega} | \mathbf{J}^*) \int Dz \frac{\partial \ell(y^1)}{\partial y^1} \bigg|_{y=\hat{y}} s_0 s_j, \quad (56)$$

and the factor $(\sum_s P_{\mathrm{Ising}}(\mathbf{s} | \mathbf{J}^*) e^{-\beta \sum_{a=1}^n \ell(y^a)})^{M-1}$ in equation (54) becomes unity when taking the $n \to 0$ limit, yielding the identical result to equation (49).

This computation naturally leads to the following question: should we compute all the zero-gradient conditions not only for $\Omega$ but also for $\bar{\Omega}$? This point is important because if this is the case, then the ansatz (34) is insufficient as it only suffices those for the active set $j \in \Omega$. To be consistent, the answer is considered to be yes in general; hence, we need to take into account the zero-gradient conditions for $k \in \bar{\Omega}$. This implies that the ansatz (34) should be modified and we need to introduce mean estimates $\bar{J}_k$ for $k \in \bar{\Omega}$ in general situations.

Fortunately, if the network is tree-like, we can show that all the zero-gradient conditions are automatically satisfied once those for $\forall j \in \Omega$ are met. Hence, the ansatz (34) is consistent on such networks; we show proof of this below. For technical reasons, we recover the external field $\mathbf{H}^*$ in the remaining of this subsection. When the external field exists, the student model should also have the external field variable, and hence
the replica result is slightly modified. That modification is accomplished by replacing $(\sum_{j \in \Omega} \bar{J}_j s_j)$ with $(\sum_{j \in \Omega} \bar{J}_j s_j + \bar{H}_0)$ in equations (46), (48) and (55). Here, $\bar{H}_0$ denotes the mean estimate of the external field variable acting on the focused spin $s_0$ of the student model, and is determined by the extremization condition of the free energy, yielding

$$0 = \sum_{s_0, s_\Omega} P_{\text{Ising}}(s_0, \mathbf{s}_\Omega | \mathbf{J}^*, \mathbf{H}^*) \int Dz \frac{\partial \ell}{\partial y} \Big|_{y=\hat{y}} s_0. \quad (57)$$

Under the above setup, we show the consistency of equation (34) on tree-like networks. The first step is to write down the zero-gradient condition for $k \in \bar{\Omega}$. The result of applying the averages and replica method is the replacement of $s_j$ with $s_k$ in equation (54). With this expression, we perform the following transformation:

$$\begin{aligned} & \sum_s P_{\text{Ising}}(\mathbf{s} | \mathbf{J}^*, \mathbf{H}^*) e^{-\beta \sum_{a=1}^n \ell(y^a)} \frac{\partial \ell(y^1)}{\partial y^1} s_0 s_k \\ & = \frac{\partial}{\partial H_k^*} \left( \sum_s P_{\text{Ising}}(\mathbf{s} | \mathbf{J}^*, \mathbf{H}^*) e^{-\beta \sum_{a=1}^n \ell(y^a)} \frac{\partial \ell(y^1)}{\partial y^1} s_0 \right) \\ & \quad + \left( \sum_s P_{\text{Ising}}(\mathbf{s} | \mathbf{J}^*, \mathbf{H}^*) e^{-\beta \sum_{a=1}^n \ell(y^a)} \frac{\partial \ell(y^1)}{\partial y^1} s_0 \right) \langle s_k \rangle, \end{aligned} \quad (58)$$

where $\langle \cdots \rangle$ denotes the average over $P_{\text{Ising}}(\mathbf{s} | \mathbf{J}^*, \mathbf{H}^*)$. By following the same computations so far, the second term vanishes in the limit $\lim_{\beta \to \infty} \lim_{n \to 0} \lim_{N \to \infty}$ because the coefficient converges to the right-hand side of (57) giving zero. Meanwhile, in the same limit, the first term can be transformed as

$$\begin{aligned} & \frac{\partial}{\partial H_k^*} \left( \sum_s P_{\text{Ising}}(\mathbf{s} | \mathbf{J}^*, \mathbf{H}^*) e^{-\beta \sum_{a=1}^n \ell(y^a)} \frac{\partial \ell(y^1)}{\partial y^1} s_0 \right) \to \frac{\partial}{\partial H_k^*} \left( \sum_{s_0, s_\Omega} P_{\text{Ising}}(s_0, \mathbf{s}_\Omega | \mathbf{J}^*, \mathbf{H}^*) \right. \\ & \qquad \left. \times \int Dz \frac{\partial \ell(y)}{\partial y} \Big|_{y=\hat{y}} s_0 \right), \end{aligned} \quad (59)$$

and the dependence on $H_k^*$ appears only in the marginal distribution $P_{\text{Ising}}(s_0, \mathbf{s}_\Omega | \mathbf{J}^*, \mathbf{H}^*)$. On tree-like networks, the marginal distribution necessarily takes the following form:

$$P_{\text{Ising}}(s_0, \mathbf{s}_\Omega | \mathbf{J}^*, \mathbf{H}^*) = \frac{1}{Z} e^{s_0 (\sum_{i \in \Omega} J_i^* s_i + H_0^*) + \sum_{i \in \Omega} h_i^{\Omega 0} s_i}, \quad (60)$$

where $h_i^{\Omega 0}$ is the effective field obtained by marginalizing the descendant spins of $i$, and is usually termed as cavity field. An important point of equation (60) is the absence of higher-order interactions among active set spins because of the tree-like structure. Hence, the differentiation of $H_k^*$ appears only through that of the effective fields. Furthermore, owing to the tree-like structure, only one of the effective fields is dependent on $H_k^*$. Specifying the corresponding index as $j \in \Omega$), we get
$$ 
\begin{aligned}
& \frac{\partial}{\partial H_k^*} \left( \sum_{s_0, s_\Omega} P_{\text{Ising}}(s_0, \mathbf{s}_\Omega | \mathbf{J}^*, \mathbf{H}^*) \int Dz \frac{\partial \ell(y)}{\partial y} \bigg|_{y=\hat{y}} s_0 \right) \\
& = \frac{\partial h_j^{\perp 0}}{\partial H_k^*} \frac{\partial}{\partial h_j^{\perp 0}} \left( \sum_{s_0, s_\Omega} P_{\text{Ising}}(s_0, \mathbf{s}_\Omega | \mathbf{J}^*, \mathbf{H}^*) \int Dz \frac{\partial \ell(y)}{\partial y} \bigg|_{y=\hat{y}} s_0 \right) \\
& = \frac{\partial h_j^{\perp 0}}{\partial H_k^*} \{(\text{rhs of equation (49)}) - (\text{rhs of equation (57)}) \times \langle s_j \rangle\} = 0. \quad (61)
\end{aligned}
$$

Hence, the zero-gradient conditions on the inactive set $\bar{\Omega}$ are satisfied once those of the active set $\Omega$ hold, proving our statement.

The above proof also provides a perspective for loopy graphs. If loops exist, then higher-order interactions emerge in $P_{\text{Ising}}(s_0, \mathbf{s}_\Omega)$; they generally depend on $H_k^*$ in a complex manner and yield some additional terms as a result of differentiation. In such situations, additional mean estimates $\bar{J}_k$ for $k \in \bar{\Omega}$ will be necessary to satisfy the corresponding zero-gradient conditions; however, treating all variables in $\bar{\Omega}$ is clearly infeasible. Tailoring good approximations in such cases may be interesting in future work, although in section 4.3 we show an example in which our present theoretical treatment becomes a good approximation even for loopy graphs.

## 4. Numerical experiments

In this section, we conduct numerical experiments to check the accuracy of the theoretical computations. The actual behavior of the order parameters and related quantities depends on the details of the coupling ensembles. Hence, we treat the regular random (RR) graph and Erdős–Rényi (ER) graph as representative examples of sparse tree-like graphs. The RR graph is characterized by one connectivity parameter $c$, while the ER graph is characterized by the connection probability $p$. To keep the generated graph sparse enough in the ER case, we assume the probability is scaled as $p = d/N$, yielding the mean degree $d$. Furthermore, we also assume that the couplings of the teacher model have the same probability of taking both signs and the strength is constant: $|J_i^*| = K > 0$. The coupling strength $K$ is assumed to be small enough to satisfy the paramagnet assumption of the teacher model. In particular, for the RR graph, the paramagnetic condition is

$$ (c - 1) \tanh^2 K < 1, \quad (62) $$

while that of the ER one is

$$ d \tanh^2 K < 1. \quad (63) $$

For readers interested in the derivation, please refer to [31, 32]. The cost function is fixed to that of PL in the following, as the simplest and commonly used case. The result of the RR graph case is shown below in section 4.1, and that of the ER graph is in section 4.2. For comparison, some numerical results on the square lattice are shown in section 4.3.
focusing on the approximation nature of the present theoretical results. Furthermore, as another common cost function, the so-called interaction screening (IS) method [21] belonging to the local learning class is examined and compared with PL in section 4.4.

Owing to the uniformity of the coupling strength, the strength of mean estimates $\{\bar{J}_i\}_{i \in \Omega}$ can also be set to a uniform value $|\bar{J}_i| = \bar{K} = \hat{b}K$, where the bias factor

$$ \hat{b} \equiv \bar{K}/K, \qquad (64) $$

is introduced. By the same reason, the marginal distribution can be simplified by again introducing the cavity field $h^* = \sum_{j \in \Omega} J_j^* s_j$ as

$$ \sum_{s_0, \boldsymbol{s}_{\Omega}} P(s_0, \boldsymbol{s}_{\Omega} | \boldsymbol{J}^*)(\cdots) = \sum_{s_0} \int dh^* P_{\text{cav}}(h^* | \boldsymbol{J}^*) \frac{e^{s_0 h^*}}{Z_0}(\cdots), \qquad (65) $$

where $Z_0 = \int dh^* P_{\text{cav}}(h^* | \boldsymbol{J}^*) 2 \cosh h^*$. If the focused spin's connectivity is $c$, then the cavity field distribution becomes

$$ P_{\text{cav}}(h^* | \boldsymbol{J}^*) = P_{\text{cav}}(h^* | K, c) \equiv \frac{1}{2^c} \sum_{k=0}^{c} \binom{c}{k} \delta(h^* - K(c - 2k)). \qquad (66) $$

Applying the reduction (65) in equations (47) and (49) with replacement $\sum_{j \in \Omega} \bar{J}_j s_j \to \hat{b}h^*$ in equation (48) reduces the computation of mean estimates to that of the bias factor $\hat{b}$. The theoretically evaluated $\hat{b}$ was compared with that obtained by numerical experiments to check the validity of our theoretical treatment.

The numerical computation of the order parameter $Q$ will be conducted below, but it has some delicate points. In our actual computations, the following procedure was adopted: from the generated teacher model we first compute the inverse correlation function $(C^{\Omega})^{-1}$ by the cavity formula (52) and numerically invert it to obtain $C^{\Omega}$; then we introduce $\{\hat{\Delta}_i = \hat{J}_i - \bar{J}_i\}_{i=1}^{N-1}$ from the learning result $\hat{\boldsymbol{J}}$ and the mean estimate $\bar{\boldsymbol{J}}$ which is obtained as $\bar{\boldsymbol{J}} = \hat{b}\boldsymbol{J}^*$ into which the theoretically evaluated value of $\hat{b}$ is inserted; finally we get a numerical value of $Q$ by $Q = \sum_{i,j} C_{ij}^{\Omega} \hat{\Delta}_i \hat{\Delta}_j$. Although it is also possible to evaluate $C^{\Omega}$ by the Monte Carlo (MC) sampling instead of using the formula (52), this method is better for controlling fluctuations and reducing computational cost.

The outline of our numerical experiment is as follows. We first generate teacher model, next perform MC sampling, and finally choose a center spin and conduct learning by numerically minimizing the PL cost function (5). The obtained estimator is used to compute relevant quantities such as RSS. As a result, there are some distinctive sources of fluctuation in the estimate, but we do not discriminate them below. The error bar is accordingly defined by the standard error coming from those fluctuations. The number of datasets used to compute the error bar is hereafter denoted as $N_{\text{set}}$. More details of the numerical experiment are summarized in appendix C.
Learning performance in inverse Ising problems with sparse teacher couplings

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Plots of $\mathcal{E}$ (left), $Q$ (middle), and $\hat{b}$ (right) against $\alpha$ for $K = 0.2$ and $0.4$ at $(N, c) = (200, 3)$. Dotted lines and color markers are the theoretical and numerical values, respectively. The agreement between them is fairly good. The left and middle panels are plotted in the double log scale because $\mathcal{E}$ and $Q$ drastically diverge in the limit $\alpha \to 2$. The error bars obtained from $N_{\text{set}} = 100$ datasets are shown, although they tend to be comparable with the size of markers.</figcaption></figure>

### 4.1. RR graph case

In the case of the RR graph with connectivity $c$, using equation (52) the trace of the inverse correlation function becomes

$$ \frac{1}{N} \text{Tr} C^{-1} = \frac{c}{1 - \tanh^2 K} - c + 1. \qquad (67) $$

Substituting this in conjunction with the parameters obtained by equations (47)–(49) into equation (50), we obtain the RSS. Below, we compare these theoretical values with the numerically evaluated ones.

We start by comparing the theoretical and numerical values of $\mathcal{E}$, $Q$, and $\hat{b}$. In figure 2, these quantities are plotted against $\alpha$ for $K = 0.2$ and $0.4$ at $N = 200$ and $c = 3$. In all the plots, the agreement between the theoretical (dotted lines) and numerical (color markers) results is fairly good, supporting the validity of our analytical treatment.

Next, we consider the distributions of the estimators in figure 3, which were normalized as probability distribution functions. The left panel is the distribution of the estimators on the active set $\Omega$. We can observe that two peaks are located around the theoretical prediction $\pm \hat{b}K$. In the middle panel, the estimator distribution on the inactive set $\bar{\Omega}$ is shown, yielding a Gaussian-like distribution with zero mean. Similar behavior is observed for the noise part on the active set, $\{\hat{\Delta}_i = \hat{J}_i - \bar{J}_i\}_{i \in \Omega}$, the distribution of which is given in the right panel. Here, the mean estimates $\{\bar{J}_i\}_{i \in \Omega}$ are computed by multiplying the theoretically evaluated bias $\hat{b}$ by the true coupling $\{J_i^*\}_{i \in \Omega}$. These observations are again consistent with our theoretical analysis.

Thirdly, we check the finite size effect. In figure 4, against the system size $N$, the RSS and rescaled variance (multiplied by $N$) of the noise parts $\hat{\Delta} = \hat{J} - \bar{J}$ are plotted in the
Learning performance in inverse Ising problems with sparse teacher couplings

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Distribution of the estimators $\hat{J}$ on the active and inactive sets are given in the left and middle panels, respectively. The right panel is the distribution of the noise part on the active set, $\{\hat{\Delta}_i = \hat{J}_i - \bar{J}_i\}_{i \in \Omega}$. The system parameters are $(N, K, \alpha, c) = (200, 0.4, 5, 3)$. The middle and right panels imply that the noise parts obey the zero-mean Gaussian distribution and have no discriminative difference between the active and inactive sets. Here, the histograms are generated from $N_{\text{set}} = 500$ datasets; from each dataset, the number of obtained estimators is $c = 3$ for $\Omega$ while that for $\bar{\Omega}$ is $N - c - 1 = 196$.</figcaption></figure>

upper and lower panels, respectively. Although the finite size effect behaves in different ways depending on the parameters and quantities, we can see that the numerical results (markers) fairly matched the theoretical values (black dotted lines) as the system size is large. Here, the rescaled variance corresponds to the quantity $Q\text{Tr}(C^{\top 0})^{-1}/N$ in our theoretical computation, which is consistent with equation (50). These results again confirm the validity of our computations.

Finally, we have some noteworthy remarks. The results shown in figures 3 and 4 imply the possibility of an efficient method of debiasing and hypothesis testing. The bias factor $\hat{b}$ can be computed from our analytical result, and hence we can debias our estimator in an efficient manner. The residual after debiasing $\hat{\Delta}$ is considered to obey a Gaussian distribution, as shown in figure 3, and is supported by our analytical computations in appendix A. Thus, we can efficiently compute the P-value according to the standard hypothesis testing method, enabling us to judge the relevance of the estimated couplings. Moreover, in the thermodynamic limit $N \to \infty$, we can show that the perfect reconstruction of the teacher's network is possible for any $\alpha > 2$. To do so, we need to evaluate the probability of getting false positives in the estimator. To control false positives, we introduce a constant threshold value $K_{\text{th}}(> 0)$, and consider estimated couplings with absolute values less than $K_{\text{th}}$ as negligible and set to zero; we independently repeat this procedure for all $i = 1, \dots, N$. Let us evaluate the probability of successfully screening out false positives using this method. The observations so far imply, on the inactive set $\bar{\Omega}$, that the estimator behaves as

$$ \hat{J}_i \sim \mathcal{N}\left(0, \frac{\sigma_i^2}{N}\right), (\forall i \in \bar{\Omega}), \qquad (68) $$

where $\sigma_i^2(= O(1))$ is the rescaled variance of the estimate, with relation
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (Upper): plot of $\mathcal{E}$ against the system size $N$ for $\alpha = 5$ (left) and 50 (right) at $(K, c) = (0.4, 3)$. The black dotted lines denote the theoretical result and the markers are the numerical ones. The numerical results tend to converge with the theoretical results as the system size grows, although the finite size effect seems to be different between the left and right panels. The error bar is obtained from $N_{set} = 500$ datasets for $N = 50-200$, $N_{set} = 400$ for $N = 400$, and $N_{set} = 50$ for $N = 800$. (Lower): the rescaled variance (multiplied by $N$) of the noise part $\hat{\Delta} = \hat{J} - \bar{J}$ is plotted against the system size $N$. The parameters are the same as those of their counterparts in the upper panels. Although in this closeup scale there is a small gap between the numerical and theoretical results within the one standard error, this gap can be eliminated by taking a larger number of samples. Here, the error bar was obtained using the bootstrap method by considering each realization and component of $\hat{\Delta}$ as i.i.d.</figcaption></figure>

$(1/N) \sum_{i \in \bar{\Omega}} \sigma_i^2 \approx Q \text{Tr}(C^{\backslash 0})^{-1}/N$. Hence, the probability of successfully screening out these estimators on $\bar{\Omega}$ is

$$ \prod_{i \in \bar{\Omega}} \text{Prob}(|\hat{J}_i| < K_{\text{th}}) = \prod_{i \in \bar{\Omega}} \left(1 - 2 \int_{\sqrt{\sigma_i^2}}^{\infty} dz \frac{e^{-\frac{1}{2}z^2}}{\sqrt{2\pi}}\right) \\ \approx \prod_{i \in \bar{\Omega}} \left(1 - \frac{2}{\sqrt{2\pi}} \frac{e^{-\frac{1}{2}\frac{N}{\sigma_i^2}K_{\text{th}}^2}}{\sqrt{\frac{N}{\sigma_i^2}K_{\text{th}}}}\right) \rightarrow 1, (N \rightarrow \infty). \quad (69) $$
The second approximate equality comes from the asymptotic formula of the integral, which can be justified for large $N$. The last limiting behavior holds as long as $\sigma_i$ is bounded from above, because the exponential factor $e^{-\frac{1}{2} \frac{N}{\sigma_i^2} K_{\text{th}}^2}$ decays fast enough compared with the number of products $|\bar{\Omega}| = N - c - 1$. Hence, we can completely suppress the false positives in the limit $N \to \infty$. Meanwhile, we also desire to accurately reproduce the presence of couplings on $\Omega$. This can be done by tuning the threshold value $K_{\text{th}}$ as smaller than the true coupling strength $K$ (the mean estimates $\bar{J}$ are larger than $K$ in the absolute value). In practical situations, we do not know the true coupling strength $K$ in advance, and thus it is nontrivial to correctly tune $K_{\text{th}}$. In such cases, it may be better to tune $K_{\text{th}}$ by monitoring the distribution of estimators such as figure 3, and to find a value that effectively separates the modes of distribution. The present theoretical analysis supports this process by manifesting the behavior of estimators in the limit $N \to \infty$.

## 4.2. ER graph case

For the ER graph with connection probability $p = d/N$, the evaluation of the order parameters and related quantities is slightly more complex than the RR case because of the distributed nature of the connectivity. In the thermodynamic limit, the distribution of connectivity $c$ in the ER graph obeys the Poisson distribution:

$$P_{\text{po}}(c|d) = e^{-d} \frac{d^c}{c!} \qquad (70)$$

The trace of the inverse correlation function fortunately becomes simple in the limit:

$$\frac{1}{N} \text{Tr} C^{-1} \xrightarrow{N \to \infty} \sum_{c=0}^{\infty} \left( \frac{c}{1 - \tanh^2 K} - c + 1 \right) P_{\text{po}}(c|d) = \left( \frac{d}{1 - \tanh^2 K} - d + 1 \right) . \quad (71)$$

When focusing on spin $i$ with connectivity $c_i$ in the ER graph, its associated order parameters are computed by equation (47) with $P(h^*|c_i)$ defined in equation (66), and the RSS is given by

$$\mathcal{E}_i(c_i) = (1 - \hat{b}(c_i))^2 K^2 + Q(c_i) \left( \frac{d}{1 - \tanh^2 K} - d + 1 \right) . \qquad (72)$$

This explicit dependence of the order parameter on $c_i$ is the complex point of the ER case. Upon realizing this property, we can easily compute the mean RSS for the whole network, which is written as

$$\mathcal{E}_{\text{mean}} = \frac{1}{N} \sum_{i=1}^{N} \mathcal{E}_i \xrightarrow{N \to \infty} \sum_{c=0}^{\infty} \left\{ (1 - \hat{b}(c))^2 K^2 + \left( \frac{d}{1 - \tanh^2 K} - d + 1 \right) Q(c) \right\} P_{\text{po}}(c|d) . \qquad (73)$$

These provide explicit formulas of the RSSs for the ER case.
Learning performance in inverse Ising problems with sparse teacher couplings

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Plots of $\mathcal{E}$ (left), $Q$ (middle), and $\hat{b}$ (right) against $c$ at $(N, \alpha, d) = (400, 10, 4)$. Black dotted lines and markers are the theoretical and numerical values, respectively; the different colors correspond to different $K$. The agreement between them is fairly good.</figcaption></figure>

As an interesting departure from the RR case, we here examine the connectivity dependence of our quantities of interest. The plots of $\mathcal{E}(c)$, $Q(c)$, and $\hat{b}(c)$ at $(N, \alpha, d, K) = (400, 10, 4, 0.4)$ are given in figure 5. In this experiment, we generated ten different ER networks, performed two independent MC samplings, and conducted learning for all $i = 1, \dots, N$. The error bars were placed using the obtained datasets, and thus $N_{\text{set}}$ varied among different connectivity $c$. The agreement between the theoretical and numerical results is fairly good, supporting our theoretical result. Although a slight deviation at large $c$ in $\mathcal{E}(c)$ and $Q(c)$ was observed, this was presumably attributed to the finite size effect, which increased at large $c$ because of insufficient system size for generating nodes with large $c$. We have tried to control this deviation but found it is difficult to conduct experiments of sufficiently large systems in reasonable time: the generation probability of node with, say, $c = 13$ can be estimated as $P_{\text{po}}(13|4) \approx 2 \times 10^{-4}$, and hence for stably generating networks with such large degree nodes we need at least $N \approx 5000$, which is too much in our experiment. We thus leave this problem untouched.

We also computed the mean RSS (73) for the whole network. The theoretical value is $\mathcal{E}_{\text{mean}} = 0.4780$, while the present experimental value is $\mathcal{E}_{\text{mean}} = 0.4907 \pm 0.0041$. The slight difference between these is again attributed to the finite size effect. Here, the theoretical value was obtained by taking the sum of equation (73) up to $c = 20$; the effect of this truncation was found to be small.

### 4.3. Square lattice case for comparison

The cavity method in direct problems is known to yield good approximations even for loopy graphs, when correlations among spins are weak; it is sometimes referred to as Bethe approximation. Here, we examined this approximation nature of the present theoretical computation of inverse problems. To this end, we compared our theoretical result for $c = 4$ with the simulation result on the square lattice with periodic boundary condition. To avoid possible complexity due to frustration, the present teacher couplings were assumed to be all positive and constant, $J_i^* = K > 0$, $(i \in \Omega)$.

In figure 6, we plotted $\mathcal{E}$ and $\hat{b}$ against $\alpha$ for $K = 0.2$ on the square lattice of size $20 \times 20$, in comparison with our theoretical result (dotted line) computed with the
Learning performance in inverse Ising problems with sparse teacher couplings

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Plots of $\mathcal{E}$ (left) and $\hat{b}$ (right) against $\alpha$ for $K = 0.2$ on the square lattice of size $20 \times 20$. For comparison, the theoretical results derived by assuming the tree-like structure of the coupling network are plotted as the dotted lines. The agreement between the markers (numerical results) and lines is fairly good. The error bars obtained from $N_{set} = 400$ datasets are shown.</figcaption></figure>

assumption of the tree-like network structure. The agreement between the theoretical and numerical results is fairly good. This indicates that our theoretical result can be a good approximation even for loopy graphs.

Another interesting phenomenon for loopy graphs is the possible presence of bias in the estimated couplings for spins in $\bar{\Omega}$, as discussed in section 3.3. To examine this, in the upper panels of figure 7 we show the distributions of the coupling estimates corresponding to the next nearest neighbors (NNN) from the center spin $s_0$ in the teacher model for the square lattice (left) and for the RR graph with $c = 4$ (right). To make a fair comparison, the present teacher couplings for the RR graph case are all positive and constant as the square lattice case. These two distributions are very similar, implying that the bias in coupling estimates for remote spins is, even if it exists in loopy graphs, very weak for the present situation. For further quantitative information, the means of those distributions were plotted against the system size in the lower panels. Again, we observed no clear deviation from zero and no significant difference between the two cases of the square lattice and RR graph. These suggest the practicality of the theoretical results for wider situations than tree-like networks.

### 4.4. Comparison with interaction screening

We here examine the IS cost function [21, 22], which is another common method for the inverse Ising problem. The IS cost function is given by

$$ \ell^{\text{IS}}(x) = e^{-x}. \qquad (74) $$

The quantitative comparison with the PL method is of our interest.

In figure 8, we plot $\mathcal{E}$, $Q$, and $\hat{b}$ against $\alpha$ for the RR graph at $(c, K) = (3, 0.4)$, with two theoretical curves of IS (dashed line) and PL (dotted line). Numerical results (asterisks) are also shown to validate the theoretical result of the IS case. The important observation is that the IS result consistently gives larger values of $\mathcal{E}$, $Q$, $\hat{b}$ than those of PL. This implies that the IS method gives larger error, variance, and bias than PL.
Learning performance in inverse Ising problems with sparse teacher couplings

<figure><img src="image_10.png" /><figcaption>**Figure 7.** (Upper): distributions of the NNN estimators for the 20 × 20 square lattice (left) and for the RR graph with (N, c) = (400, 4). In both cases, other parameters are set to be (K, α) = (0.2, 5) and Nset = 400. No clear positive/negative tendency is observed in both cases. (Lower): plots of the mean of the NNN estimate distribution against the system size for the square lattice (left) and RR graph (right). The other parameters are similar to those of the corresponding upper panels. The means are quite small, and no clear deviation from zero is observed. The dataset sizes are Nset = 600, 600, 400, 200, 40 for N = 100, 225, 400, 900, 1600, respectively. The error bars are obtained using the bootstrap method.</figcaption></figure>

Although it is known that IS achieves the optimal scaling of the dataset size *M* [21, 22, 35] necessary for perfect reconstruction of the coupling network structure, it does not necessarily mean the better accuracy in terms of variance/bias. The present analysis demonstrates this in a concise way.

**5. Summary and discussion**

We proposed a theory to evaluate the reconstruction performance in inverse Ising problems with sparse couplings by maximizing of the pseudolikelihood in the thermodynamic limit. A large part of the theory relies on the statistical mechanical formulation in [24], but we refined the theoretical treatment in the cavity method to handle the teacher model with sparse couplings. The resultant expression requires a full functional form of the cavity field distribution, which is far from Gaussian but was obtained by appropriate
Learning performance in inverse Ising problems with sparse teacher couplings

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of the PL and IS results. Plots of $\mathcal{E}$ (left), $Q$ (middle), and $\hat{b}$ (right) against $\alpha$ at $(c, K) = (3, 0.4)$ for the RR graph are given. The theoretical curves for PL and IS are described by dotted and dashed lines, respectively. The color markers are the numerical results for the IS case at $N = 200$, whose error bars are obtained from $N_{\text{set}} = 100$ datasets. The consistency of the numerical and theoretical results is fairly good. The IS result provides larger values in all these three quantities than those of PL, implying that IS gives larger error, variance, and bias.</figcaption></figure>

consideration of the direct problem counterpart. The theoretical result shows fairly good agreement with numerical experiments conducted on the RR and ER graphs, justifying our theoretical treatment. This agreement holds even for the case of the square lattice, suggesting the practicality of the present result as an approximation for loopy graphs.

The crucial assumptions of our treatment are the asymptotic behavior of the estimator (34) and the paramagnet assumption of the teacher model, leading to the decoupled distributions of the cavity fields. The former assumption implies that the teacher's couplings can be reconstructed by the student almost perfectly, as discussed at the end of section 4.1 according to the hypothesis testing framework, providing a theoretical reasoning to use the inverse Ising framework. The latter assumption requires the smallness of the coupling strength, implying that strongly-correlated datasets cannot be treated by the proposed theory. It will be a challenge to overcome this applicability limitation.

Concerning to the perfect reconstruction of the sparse network in the inverse Ising framework, earlier studies reported similar results in empirical and theoretical ways [29, 35–37]. Especially, a series of analyses by Wainwright and the collaborators [35–37] derived the necessary and sufficient conditions for the perfect reconstruction in the high-dimensional limit, clarifying that the necessary size of the dataset scales as $M = O(\log N)$ when the maximum degree of the network is bounded from above. Compared to this scaling, our result on the scaling $M = O(N)$ is rather conservative. Our formulation, however, has some nontrivial advantages by deriving more detailed information about the system. For example, our formulation can treat the ER graph, the maximum degree of which is not bounded and the proofs established in [35–37] are not applicable. Our result also clarifies by directly assessing the estimator's fluctuation that the hypothesis testing can actually achieve the perfect reconstruction, which provides another efficient way of reconstruction than the $\ell_1$ regularization used in the above earlier studies. The explicit computation of the bias on the estimator is also another profit. In this way, the present formulation can provide finer information in concrete cases and will be beneficial in wider situations.
For handling real-world datasets, finite magnetizations as well as possible loop structures in the network should be taken into account. For such realistic situations, the computation of $(1/N)\text{Tr}C^{-1}$ and $P_{\text{cav}}(h^*|J^*)$ will be more complicated. For evaluating those quantities, advanced techniques such as Bethe approximation [33, 34], high-temperature expansion, and MC samplings will be useful. The ansatz (34) should be also modified for the case of loopy graphs, as discussed in section 3.3. The presented result can be still practical as an approximation for treating such situations, as demonstrated in section 4.3. Certain data analysis utilizing these theoretical results will be interesting and useful.

A clear drawback of the estimator treated in this paper is that it is not informative in the region $\alpha \le 2$, as indicated by the divergent RSS in the limit $\alpha \to 2+0$ shown in section 4 and [24]. To overcome this, the use of regularizations will be promising. The $l_1$ regularization will be particularly useful to control false positives in the estimated couplings. It is also possible to employ hypothesis testing in conjunction with $l_2$ regularization [38]. An extension of the present analysis to these cases is interesting and will be our focus in future work.

Another interesting extension of the present analysis might be the analysis of the model-mismatched cases where the student model cannot be equal to the teacher one. Even in such cases, some limited information in the teacher, such as the coupling network structure, might be recovered in some conditions [10]. Pursuing this possibility will be interesting because it can provide a better justification for applications of the inverse Ising framework to the analysis of real-world datasets. The analysis of model-mismatched cases may also have a connection to ‘criticality’ observed in earlier studies [39].

The inverse Ising problem or Boltzmann machine has been treated in this paper. Although this model is much simpler than the current models of machine learning communities, we believe that it is important to enhance theoretical knowledge on such simple models to maintain the reliability and interpretability of the results given by machine learning technologies. We hope that our study contributes to this direction, which will lead to a better understanding of and relationship with more complex models.

## Acknowledgment

This work was supported by JSPS KAKENHI Nos. 25120013, 17H00764, 18K11463 and 19H01812, JST CREST Grant Number JPMJCR1912, Japan, and RIKEN Special Postdoctoral Researcher (SPDR) Grant No. 8232.

## Appendix A. Computations for L and S

For computing $L$, the following decomposition of the cavity fields becomes useful:

$$h^* = \sqrt{Q^* - \frac{m^2}{q}} v_* + \sqrt{\frac{m^2}{q}} z, \qquad (75)$$
$$h^a = \sqrt{Q - q} v_a + \sqrt{q} z, \qquad (76)$$

where $v_a, v_*, z$ are i.i.d. Gaussian variables with zero mean and unit variance. It is easy to confirm that this decomposition reproduces the covariances among {$h^*, h^1, \dots, h^n$}. Using this and performing the integration with respect to $v_*$, we get

$$L(Q^*, Q, q, m) = \int Dz e^{\sqrt{\frac{m^2}{q}}z - \frac{1}{2}\frac{m^2}{q}} \left( \int Dv e^{-\beta\ell(\sqrt{Q-q}v + \sqrt{q}z)} \right)^n, \qquad (77)$$

where we use the relation $Z_0 = 2e^{\frac{1}{2}Q^*}$, which was canceled with a factor appearing by the integration of $v_*$. Equation (26) is easily derived from this.

For computing the entropic term $S(C^{\lambda^0}, J^*, Q, q, m)$, we use the rescaled variable $W = \sqrt{N}J$ and set the integration measure as $\text{Tr}_J = \int dW$. Here, we use the uniform measure because in the present setting the student has no prior information about the teacher couplings. If certain prior knowledge is available such as the teacher coupling sparseness, it can be suitable to introduce another measure such as Laplace prior. Further, we represent the delta functions by the Fourier expressions as follows:

$$\delta \left( Q - \frac{1}{N} \sum_{i,j} C_{ij}^{\lambda^0} W_i^a W_j^a \right) = C_1 \int d\tilde{Q} e^{\frac{1}{2} N \tilde{Q} Q - \frac{1}{2} \tilde{Q} \sum_{i,j} C_{ij}^{\lambda^0} W_i^a W_j^a}, \qquad (78a)$$

$$\delta \left( q - \frac{1}{N} \sum_{i,j} C_{ij}^{\lambda^0} W_i^a W_j^b \right) = C_2 \int d\tilde{q} e^{-N\tilde{q}q + \tilde{q} \sum_{i,j} C_{ij}^{\lambda^0} W_i^a W_j^b}, \qquad (78b)$$

$$\delta \left( m - \frac{1}{N} \sum_{i} C_{ij}^{\lambda^0} W_{i}^{*} W_{j}^{a} \right) = C_2 \int d\tilde{m} e^{-N\tilde{m}m + \tilde{m} \sum_{i} C_{ij}^{\lambda^0} W_{i}^{*} W_{j}^{a}}, \qquad (78c)$$

where the integration contour is the imaginary axis and $C_1, C_2$ are appropriate normalization constants; however, these points are irrelevant and ignored hereafter. Inserting this into equation (23), we get

$$e^{NS} = \int d\tilde{Q} d\tilde{q} d\tilde{m} e^{S_X} \int \prod_a d\mathbf{W}^a e^U, \qquad (79)$$

where

$$S_X = N \left( \frac{1}{2} n \hat{Q} Q - \frac{1}{2} n(n-1) \tilde{q} q - n \tilde{m} m \right), \qquad (80)$$

$$\begin{aligned} U &= -\frac{1}{2} \tilde{Q} \sum_a (\mathbf{W}^a)^\top C^{\lambda^0} \mathbf{W}^a + \tilde{q} \sum_{a<b} (\mathbf{W}^a)^\top C^{\lambda^0} \mathbf{W}^b + \tilde{m} \sum_a (\mathbf{J}^*)^\top C^{\lambda^0} \mathbf{W}^a \\ &= -\frac{1}{2} (\tilde{Q} + \tilde{q}) \sum_a (\mathbf{W}^a)^\top C^{\lambda^0} \mathbf{W}^a + \frac{1}{2} \tilde{q} \sum_{a,b} (\mathbf{W}^a)^\top C^{\lambda^0} \mathbf{W}^b + \tilde{m} \sum_a (\mathbf{J}^*)^\top C^{\lambda^0} \mathbf{W}^a. \end{aligned} \qquad (81)$$

To decouple different replicas and components of {$\mathbf{W}^a$}$_a$, we use the expression $C^{\lambda^0} = O^\top \Lambda O$, where $\Lambda$ is the diagonal matrix consisting of the eigenvalues {$\lambda_i$}$_i$ and $O$ is the
appropriate orthogonal matrix. Performing the variable transformation $\tilde{W} = O W$ and applying the Hubbard-Stratonovich transformation, we get

$$ 
\begin{aligned}
\int \prod_a d\mathbf{W}^a e^U &= \int \prod_i Dz_i \int \prod_a d\tilde{\mathbf{W}}^a e^{-\frac{1}{2}(\tilde{Q}+\tilde{q})\sum_a \sum_i \lambda_i (\tilde{W}_i^a)^2 + \sum_a \sum_i (\lambda_i \tilde{W}_i^* \tilde{m}_i + \sqrt{\lambda_i} \tilde{q} z_i) \tilde{W}_i^a} \\
&= \int \prod_i Dz_i e^{\sum_i \left\{ \frac{1}{2} \frac{(\sqrt{\lambda_i} \tilde{W}_i^* \tilde{m}_i + \sqrt{\tilde{q}} z_i)^2}{(\tilde{Q}+\tilde{q})} + \frac{1}{2} (\log 2\pi - \log \lambda_i - \log(\tilde{Q}+\tilde{q})) \right\}} \\
&= e^{-\frac{N}{2} \log\left(1-\frac{n\tilde{q}}{\tilde{Q}+\tilde{q}}\right) + \frac{1}{2} n \sum_i \frac{\lambda_i (\tilde{W}_i^*)^2 \tilde{m}_i^2}{\tilde{Q}+\tilde{q}(1-n)} + \frac{N}{2} \sum_i (\log 2\pi - \log \lambda_i - \log(\tilde{Q}+\tilde{q}))} \equiv e^{S_J}. \quad (82)
\end{aligned}
$$

Note that this quadratic form with respect to $\tilde{W}$ implies that $\tilde{W}$ essentially obeys Gaussian, and thus the estimator $\hat{J}$ also does. This knowledge of the distribution can be used for hypothesis testing as addressed in the main text.

In the thermodynamic limit $N \to \infty$, we can use the saddle-point (or Laplace) method to avoid the explicit integrations with respect to $\tilde{Q}, \tilde{q}, \tilde{m}$. This yields

$$ S = \text{Extr}_{\tilde{Q}, \tilde{q}, \tilde{m}} \left\{ \frac{S_X + S_J}{N} \right\} = \text{Extr}_{\tilde{Q}, \tilde{q}, \tilde{m}} \left\{ \frac{1}{2} n \tilde{Q} Q - \frac{1}{2} n(n-1) \tilde{q} q - n \tilde{m} m \right. \\ \left. - \frac{1}{2} \log \left(1 - \frac{n \tilde{q}}{\tilde{Q} + \tilde{q}}\right) + \frac{1}{2} n \frac{Q^* \tilde{m}^2}{\tilde{Q} + \tilde{q}(1-n)} + \frac{n}{2} \left(\log 2\pi - \log(\tilde{Q} + \tilde{q})\right) - \frac{n}{2N} \text{Tr} \log C^{\setminus 0} \right\}. \quad (83) $$

where we used the relations $\sum_i \lambda_i (\tilde{W}_i^*)^2 = NQ^*$, $\sum_i \log \lambda_i = \text{Tr} \log C^{\setminus 0}$. The limit $n \to 0$ leads to

$$ \lim_{n \to 0} \frac{S}{n} = \text{Extr}_{\tilde{Q}, \tilde{q}, \tilde{m}} \left\{ \frac{1}{2} \tilde{Q} Q + \frac{1}{2} \tilde{q} q - \tilde{m} m + \frac{1}{2} \frac{\tilde{q} + Q^* \tilde{m}^2}{\tilde{Q} + \tilde{q}} + \frac{1}{2} \left( \log 2\pi - \log(\tilde{Q} + \tilde{q}) \right) - \frac{1}{2N} \text{Tr} \log C^{\setminus 0} \right\}, \quad (84) $$

and the extremization condition gives

$$ \tilde{Q} = \frac{Q - 2q + m^2/Q^*}{(Q - q)^2}, \quad \tilde{q} = \frac{q - m^2/Q^*}{(Q - q)^2}, \quad \tilde{m} = \frac{m/Q^*}{Q - q}. \quad (85) $$

Substituting these relations into equation (84), we obtain equation (25). If we ignore the terms related to $m$ and $\tilde{m}$, we have equation (44).
## Appendix B. Derivation of macroscopic parameters $R$ and $\rho$

To derive the expressions of $R$ and $\rho$, we can employ the technique of auxiliary variables. We introduce two terms $h_R \sum_a (W^a)^\top W^a$ and $h_\rho \sum_a (W^*)^\top W^a$ in equation (81), and perform the same line of computations as in appendix A. As a result, the entropic term is modified to the following expression:

$$ \lim_{n \to 0} \frac{S}{n} = \mathrm{Extr}_{\tilde{Q}, \tilde{q}, \tilde{m}} \left\{ \frac{1}{2} \tilde{Q} Q + \frac{1}{2} \tilde{q} q - \tilde{m} m + \frac{1}{2N} \sum_i \left( \frac{(\tilde{m} \lambda_i + h_\rho)^2 (\tilde{W}_i^*)^2 + \lambda_i \tilde{q}}{(\tilde{Q} + \tilde{q}) \lambda_i - 2h_R} + \log 2\pi - \log ((\tilde{Q} + \tilde{q})(\lambda_i - 2h_R)) \right) \right\}. \quad (86) $$

Taking the differentiation with respect to $h_\rho$ and taking the limit $h_\rho, h_R \to 0$, we get

$$ \rho = \lim_{h_\rho, h_R \to 0} \frac{\partial}{\partial h_\rho} \lim_{n \to 0} \frac{S}{n} = \frac{1}{N} \sum_i \frac{\tilde{m} (\tilde{W}_i^*)^2}{\tilde{Q} + \tilde{q}} = \frac{m}{Q^*} R^*. \quad (87) $$

The last expression is obtained by using equation (85). Similarly,

$$ \begin{aligned} R &= \lim_{h_\rho, h_R \to 0} \frac{\partial}{\partial h_R} \lim_{n \to 0} \frac{S}{n} = \frac{1}{N} \sum_i \left( \frac{\tilde{m}^2 (\tilde{W}_i^*)^2}{(\tilde{Q} + \tilde{q})^2} + \frac{\tilde{Q} + 2\tilde{q}}{(\tilde{Q} + \tilde{q})^2 \lambda_i} \right) \\ &= \left( \frac{m}{Q^*} \right)^2 R^* + \left( Q - \frac{m^2}{Q^*} \right) \frac{1}{N} \mathrm{Tr}(C^{\otimes 0})^{-1}. \end{aligned} \quad (88) $$

These give equation (32).

In the sparse case, we need to compute $\sum_{i \in \bar{\Omega}} \Delta_i^2$ for computing the RSS. By construction, this is equivalent with $R$ when $m$ is absent. Hence $\sum_{i \in \bar{\Omega}} \Delta_i^2$ is given by putting $m=0$ in equation (88), leading to equation (50).

## Appendix C. Details of numerical experiments

The actual experimental procedures are summarized as follows. We first generated a random graph and the teacher couplings on it, and obtained spin snapshots using MC sampling. Then, we randomly chose a center spin $s_0$ from the whole spins and learned the couplings connected to $s_0$ by minimizing the PL cost function defined with a dataset obtained from the sampled spin configurations. This single sequence of operations provided single values of the quantities of interest, such as $\mathcal{E}$ and $Q$. To obtain the error
bars of those quantities, we repeated this sequence many times. Here, the experiment
had three different sources of fluctuations: the generated teacher model (graph shape
and couplings), the choice of the center spin, and the MC sampling. We did not dis-
criminate between these three fluctuations unless explicitly mentioned, and we defined
the error bar as the standard error among the obtained values according to their recur-
rence; the number of datasets obtained this way is hereafter denoted as N\_set. In the MC
sampling, we started from a random initial configuration and updated the state by the
standard Metropolis method; one MC step (MCS) is defined by N trial flips of spins,
where N is the total number of spins. We discarded the first 10^5 MCSs as burn-in to
avoid systematic errors from the initialization. Furthermore, to avoid possible correla-
tions in samples, each dataset for learning was generated by subsampling from a much
larger dataset, which consists of all the configurations recorded after every few numbers
of MCS. The size of the subsampled dataset was chosen to be at least five times smaller
than that of the larger dataset. The optimization algorithm is a standard trust-region
method using the second-order expansion of the objective function.

References

[1] Schneidman E, Berry M J, Ronen S and Bialek W 2006 Weak pairwise correlations imply strongly correlated network states in a neural population *Nature* **440** 1007–12

[2] Shlens J, Field G D, Gauthier J L, Grivich M I, Dumitru P, Alexander Sher A, Litke M and Chichilnisky E J 2006 The structure of multi-neuron firing patterns in primate retina *J. Neurosci.* **26** 8254–66

[3] Tang A *et al* 2008 A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro *J. Neurosci.* **28** 505–18

[4] Hamilton L S, Sohl-Dickstein J, Huth A G, Carels V M, Deisseroth K and Bao S 2013 Optogenetic activation of an inhibitory network enhances feedforward functional connectivity in auditory cortex *Neuron* **80** 1066–76

[5] Watanabe T, Hirose S, Wada H, Imai Y, Machida T, Shirouzu I, Konishi S, Miyashita Y and Masuda N 2012 A pairwise maximum entropy model accurately describes resting-state human brain networks *Nat. Commun.* **4** 1370

[6] Watanabe T, Masuda N, Fukuda M, Kanai R and Rees G 2014 Energy landscape and dynamics of brain activity during human bistable perception *Nat. Commun.* **5** 4765

[7] Tkačik G, Marre O, Amodei D, Schneidman E, Bialek W and Berry M J 2014 Searching for collective behavior in a large network of sensory neurons *PLoS Comput. Biol.* **10** e1003408

[8] Tavoni G, Ferrari U, Battaglia F P, Cocco S and Monasson R 2015 Inferred model of the prefrontal cortex activity unveils cell assemblies and memory replay (bioRxiv)

[9] Terada Y, Obuchi T, Isomura T and Kabashima Y 2018 Objective procedure for reconstructing couplings in complex systems (arXiv:1803.04738)

[10] Terada Y, Obuchi T, Isomura T and Kabashima Y 2018 Objective and efficient inference for couplings in neuronal networks *Advances in Neural Information Processing Systems* vol 31 ed S Bengio, H Wallach, H Larochelle, K Grauman, N Cesa-Bianchi and R Garnett (Red Hook, NY: Curran Associates Inc.) pp 4971–80

[11] Julian B 1974 Spatial interaction and the statistical analysis of lattice systems *J. R. Stat. Soc. B* **36** 192–225

[12] Kappen H J and Rodríguez F 1998 Efficient learning in Boltzmann machines using linear response theory *Neural Comput.* **10** 1137–56

[13] Tanaka T 1998 Mean-field theory of Boltzmann machine learning *Phys. Rev. E* **58** 2302

[14] Hinton G E 2002 Training products of experts by minimizing contrastive divergence *Neural Comput.* **14** 1771–800

[15] Broderick T, Dudik M, Tkacik G, Schapire R E and Bialek W 2007 Faster solutions of the inverse pairwise Ising problem (arXiv:0712.2437)

[16] Sessak V and Monasson R 2009 Small-correlation expansions for the inverse Ising problem *J. Phys. A: Math. Theor.* **42** 055001

[17] Roudi Y and Hertz J Jan 2011 Mean field theory for nonequilibrium network reconstruction *Phys. Rev. Lett.* **106** 048702
[18] Mézard M and Sakellariou J 2011 Exact mean-field inference in asymmetric kinetic Ising systems *J. Stat. Mech.* L07001

[19] Cocco S and Monasson R 2011 Adaptive cluster expansion for inferring Boltzmann machines with noisy data *Phys. Rev. Lett.* **106** 090601

[20] Aurell E and Ekeberg M 2012 Inverse Ising inference using all the data *Phys. Rev. Lett.* **108** 090201

[21] Vuffray M, Misra S, Lokhov A and Chertkov M 2016 Interaction screening: efficient and sample-optimal learning of Ising models *Advances in Neural Information Processing Systems* (Red Hook, NY: Curran Associates Inc.) pp 2595–603

[22] Lokhov A Y, Vuffray M, Misra S and Chertkov M 2018 Optimal structure and parameter learning of Ising models *Sci. Adv.* **4** e1700791

[23] Vuffray M, Misra S, Andrey Y and Lokhov 2019 Efficient learning of discrete graphical models (arXiv: 1902.00600)

[24] Bachschmid-Romano L and Opper M 2017 A statistical physics approach to learning curves for the inverse Ising problem *J. Stat. Mech.* 063406

[25] Bachschmid-Romano L and Opper M 2015 Learning of couplings for random asymmetric kinetic Ising models revisited: random correlation matrices and learning curves *J. Stat. Mech.* P09016

[26] Berg J 2017 Statistical mechanics of the inverse Ising problem and the optimal objective function *J. Stat. Mech.* 083402

[27] Cover T M and Thomas J A 2012 *Elements of Information Theory* (New York: Wiley)

[28] Nguyen H C, Zecchina R and Berg J 2017 Inverse statistical problems: from the inverse Ising problem to data science *Adv. Phys.* **66** 197–261

[29] Decelle A and Ricci-Tersenghi F 2014 Pseudolikelihood decimation algorithm improving the inference of the interaction network in a general class of Ising models *Phys. Rev. Lett.* **112** 070603

[30] Hyvärinen A 2006 Consistency of pseudolikelihood estimation of fully visible Boltzmann machines *Neural Comput.* **18** 2283–92

[31] Opper M and Saad D 2001 *Advanced Mean Field Methods: Theory and Practice* (Neural Information Processing Series) (Cambridge, MA: MIT Press)

[32] Mezard M and Montanari A 2009 *Information, Physics, and Computation* (Oxford: Oxford University Press)

[33] Ricci-Tersenghi F 2012 The Bethe approximation for solving the inverse Ising problem: a comparison with other inference methods *J. Stat. Mech.* P08015

[34] Nguyen H C and Berg J 2012 Bethe–Peierls approximation and the inverse ising problem *J. Stat. Mech.* P03004

[35] Santhanam N P and Wainwright M J 2012 Information-theoretic limits of selecting binary graphical models in high dimensions *IEEE Trans. Inf. Theory* **58** 4117–34

[36] Wainwright M J, Jordan M I et al 2008 Graphical models, exponential families, and variational inference *Found. Trends Mach. Learn.* **1** 1–305

[37] Ravikumar P, Wainwright M J and Lafferty J D 2010 High-dimensional Ising model selection using $\ell_1$-regularized logistic regression *Ann. Stat.* **38** 1287–319

[38] Xu Y, Puranen S, Corander J and Kabashima Y 2018 Inverse finite-size scaling for high-dimensional significance analysis *Phys. Rev. E* **97** 062112

[39] Mastromatteo I and Marsili M 2011 On the criticality of inferred models *J. Stat. Mech.* P10012