summaries: This paper studies the problem of implementing random walks in the Massively Parallel Computation (MPC) model using logarithmic number of rounds.  Very recently, [LMOS20] develop a MPC algorithm that using O(log \ell) rounds, computes (independent) random walks of length up to \ell from every node in an arbitrary undirected graph. Their main idea is based on stitching: In order to compute walks of length r, we can stick together two walks of length r/2 and we can do this stitching procedure recursively using O(log \ell) rounds. One issue with this algorithm is that if we want to do a few walks starting from a few nodes, this algorithm still implements all walks for all nodes. In particular, for say B walks of length \ell, the algorithm of [LMOS20] uses O(mB) total memory (for all machines).  The current paper claims to solve this issue by proposing a MPC algorithm that uses a total memory of O(mn^{eps}\ell^4+Bn^{eps}\ell) and computes B walks of a given vertex r. The number of communication rounds of this new algorithm is slightly more and is O(log(\ell)/eps).

summaries: This submission gives an algorithm in the MPC model (Massively Parallel Computation) for computing independent random walks from a single vertex or a subset of vertices. A previous paper, [LMOS], on the topic gave an algorithm that allowed for generating random walks from all vertices simultaneously (more specifically, the starting points of generated walks were from the stationary distribution). Suppose that you want to generate many random walks from a single vertex. [LMOS] would require generating lots of random walks from all vertices and because of that would require a lot of total space available in the system, and therefore, it would require a lot of machines. The algorithm introduced in this submission can be significantly more space efficient in this type of setting.

summaries: The paper introduces an algorithm for computing multiple random walks in parallel. Performance is measured using the MPC model. The authors also show an interesting application to local graph clustering. The proposed algorithm improves upon memory requirements of SOTA methods at the cost a logarithmically many more iterations.

summaries: This paper studies the problem of parallelizing single-source and subset random walk computations. The primary technical contribution is an algorithm with the number of parallel rounds that depends on the length of the random walk in a logarithmic way. In addition, two applications are discussed in the article: personalized PageRank estimation and local clustering.

abstractText: Random walks are a fundamental primitive used in many machine learning algorithms with several applications in clustering and semi-supervised learning. Despite their relevance, the first efficient parallel algorithm to compute random walks has been introduced very recently (Łącki et al.) Unfortunately their method has a fundamental shortcoming: their algorithm is non-local in that it heavily relies on computing random walks out of all nodes in the input graph, even though in many practical applications one is interested in computing random walks only from a small subset of nodes in the graph. In this paper, we present a new algorithm that overcomes this limitation by building random walk efficiently and locally at the same time. We show that our technique is both memory and round efficient, and in particular yields an efficient parallel local clustering algorithm. Finally, we complement our theoretical analysis with experimental results showing that our algorithm is significantly more scalable than previous approaches.

authors:
Alekh Agarwal
Soumen Chakrabarti

title: Learning random walks to rank nodes in graphs

venue: In Proceedings of the 24th international conference on Machine learning,

authors:
Lorenzo Alvisi
Allen Clement
Alessandro Epasto
Silvio Lattanzi
Alessandro Panconesi

title: Communities, random walks, and social sybil defense

venue: Internet Math.,

authors:
R. Andersen
F. Chung
K. Lang

title: Local graph partitioning using pagerank vectors

venue: 47th Annual IEEE Symposium on Foundations of Computer Science

authors:
Sepehr Assadi
Xiaorui Sun
Omri Weinstein

title: Massively parallel algorithms for finding well-connected components in sparse graphs

venue: In Proceedings of the 2019 ACM Symposium on Principles of Distributed Computing,

authors:
Bahman Bahmani
Kaushik Chakrabarti
Dong Xin

title: Fast personalized pagerank on mapreduce

venue: In Proceedings of the 2011 ACM SIGMOD International Conference on Management of data,

title: Bookmark-coloring approach to personalized pagerank computing

venue: Technical report, Internet Mathematics,

authors:
Paul Beame
Paraschos Koutris
Dan Suciu

title: Communication steps for parallel query processing

venue: In Proceedings of the 32nd ACM SIGMOD-SIGACT-SIGAI symposium on Principles of database systems,

authors:
Sergey Brin
Lawrence Page

title: The anatomy of a large-scale hypertextual web search

venue: engine. Comput. Networks,

authors:
Ashish Chiplunkar
Michael Kapralov
Sanjeev Khanna
Aida Mousavifar
Yuval Peres

title: Testing graph clusterability: Algorithms and lower bounds

venue: IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS),

authors:
Artur Czumaj
Pan Peng
Christian Sohler

title: Testing cluster structure of graphs

venue: In Proceedings of the forty-seventh annual ACM symposium on Theory of Computing,

authors:
Fan Chung
Olivia Simpson

title: Distributed algorithms for finding local clusters using heat kernel pagerank

venue: In International Workshop on Algorithms and Models for the Web-Graph,

authors:
Jeffrey Dean
Sanjay Ghemawat

title: Mapreduce: simplified data processing on large clusters

venue: Communications of the ACM,

authors:
Atish Das Sarma
Sreenivas Gollapudi
Rina Panigrahy

title: Estimating pagerank on graph streams

venue: Journal of the ACM (JACM),

authors:
Atish Das Sarma
Anisur Rahaman Molla
Gopal Pandurangan
Eli Upfal

title: Fast distributed pagerank computation

venue: Theoretical Computer Science,

authors:
Atish Das Sarma
Danupon Nanongkai
Gopal Pandurangan
Prasad Tetali

title: Distributed random walks

authors:
Dániel Fogaras
Balázs Rácz

title: Towards scaling fully personalized pagerank

venue: Third International Workshop,

authors:
Grzegorz Gluch
Michael Kapralov
Silvio Lattanzi
Aida Mousavifar
Christian Sohler

title: Spectral clustering oracles in sublinear time

venue: In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA),

authors:
Ullas Gargi
Wenjun Lu
Vahab Mirrokni
Sangho Yoon

title: Large-scale community detection on youtube for topic discovery and exploration

venue: In Proceedings of the International AAAI Conference on Web and Social Media,

authors:
David F Gleich
C Seshadhri

title: Vertex neighborhoods, low conductance cuts, and good seeds for local community methods

venue: In Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining,

authors:
Michael T Goodrich
Nodari Sitchinava
Qin Zhang

title: Sorting, searching, and simulation in the mapreduce framework

venue: In International Symposium on Algorithms and Computation,

title: Topic-sensitive pagerank: a context-sensitive ranking algorithm for web search

venue: IEEE Transactions on Knowledge and Data Engineering,

authors:
Jeremy G Hoskins
Cameron Musco
Christopher Musco
Charalampos E Tsourakakis

title: Inferring networks from random walk-based node similarities

venue: In Proceedings of the 32nd International Conference on Neural Information Processing Systems,

authors:
Shay Halperin
Uri Zwick

title: An optimal randomised logarithmic time connectivity algorithm for the erew pram

venue: Journal of Computer and System Sciences,

authors:
Michael Isard
Mihai Budiu
Yuan Yu
Andrew Birrell
Dennis Fetterly

title: Dryad: distributed data-parallel programs from sequential building blocks

venue: In ACM SIGOPS operating systems review,

title: Simulating random walks on graphs in the streaming model

venue: In 10th Innovations in Theoretical Computer Science Conference,

authors:
Glen Jeh
Jennifer Widom

title: Scaling personalized web search

venue: In Proceedings of the 12th International Conference on World Wide Web,

authors:
David R Karger
Noam Nisan
Michal Parnas

title: Fast connected components algorithms for the erew pram

venue: In Proceedings of the fourth annual ACM symposium on Parallel algorithms and architectures,

authors:
Howard Karloff
Siddharth Suri
Sergei Vassilvitskii

title: A model of computation for mapreduce

venue: In Proceedings of the twenty-first annual ACM-SIAM symposium on Discrete Algorithms,

authors:
Bryan Klimt
Yiming Yang

title: Introducing the enron corpus

venue: In CEAS,

authors:
Jure Leskovec
Andrej

title: Krevl. SNAP Datasets: Stanford large network dataset collection

authors:
Jure Leskovec
Jon M. Kleinberg
Christos Faloutsos

title: Graph evolution: Densification and shrinking diameters

venue: ACM Trans. Knowl. Discov. Data,

authors:
Jure Leskovec
Kevin J. Lang
Anirban Dasgupta
Michael W. Mahoney

title: Community structure in large networks: Natural cluster sizes and the absence of large well-defined clusters

authors:
Jakub Łącki
Slobodan Mitrović
Krzysztof Onak
Piotr Sankowski

title: Walking randomly, massively, and efficiently

venue: In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing,

authors:
Marina Meila
Jianbo Shi

title: Learning segmentation by random walks

venue: Advances in neural information processing systems,

authors:
Lawrence Page
Sergey Brin
Rajeev Motwani
Terry Winograd

title: The pagerank citation ranking: Bringing order to the web

venue: Technical Report 1999-66, Stanford InfoLab, November

authors:
Julian Shun
Farbod Roosta-Khorasani
Kimon Fountoulakis
Michael W. Mahoney

title: Parallel local graph clustering

venue: Proc. VLDB Endow.,

authors:
Daniel A Spielman
Shang-Hua Teng

title: Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems

venue: In Proceedings of the thirty-sixth annual ACM symposium on Theory of computing,

authors:
Joyce Jiyoung Whang
David F Gleich
Inderjit S Dhillon

title: Overlapping community detection using seed set expansion

venue: In Proceedings of the 22nd ACM international conference on information & knowledge management,

title: Hadoop: The definitive guide

authors:
Jaewon Yang
Jure Leskovec

title: Defining and evaluating network communities based on ground-truth

authors:
Matei Zaharia
Mosharaf Chowdhury
Michael J Franklin
Scott Shenker
Ion Stoica

title: Spark: Cluster computing with working

venue: sets. HotCloud,

text: Random walks are a fundamental primitive used in many machine learning algorithms with several applications in clustering and semi-supervised learning. Despite their relevance, the first efficient parallel algorithm to compute random walks has been introduced very recently (Łącki et al.) Unfortunately their method has a fundamental shortcoming: their algorithm is non-local in that it heavily relies on computing random walks out of all nodes in the input graph, even though in many practical applications one is interested in computing random walks only from a small subset of nodes in the graph. In this paper, we present a new algorithm that overcomes this limitation by building random walk efficiently and locally at the same time. We show that our technique is both memory and round efficient, and in particular yields an efficient parallel local clustering algorithm. Finally, we complement our theoretical analysis with experimental results showing that our algorithm is significantly more scalable than previous approaches.

heading: 1 Introduction

text:
Random walks are key components of many machine learning algorithms with applications in computing graph partitioning [ST04, GKL+21], spectral embeddings [CPS15, CKK+18], or network inference [HMMT18], as well as learning image segmentation [MS00], ranking nodes in a graph [AC07] and many other applications. With the increasing availability and importance of large scale datasets it is important to design efficient algorithms to compute random walks in large networks.
Several algorithms for computing random walks in parallel and streaming models have been proposed in the literature. In the streaming setting, Sarma, Gollapudi and Panigrahy [DSGP11] introduced multipass streaming algorithms for simulating random walks, and recently Jin [Jin19] gave algorithms for generating a single random walk from a prespecified vertex in one pass. The first efficient parallel algorithms for this problem have been introduced in the PRAM model [KNP92, HZ96].
In a more recent line of work, Bahmani, Chakrabarti, and Xin [BCX11] designed a parallel algorithm that constructs a single random walk of length ` from every node in O(log `) rounds in the massively parallel computation model (MPC), with the important caveat that these walks are not independent (an important property in many applications). This was followed by the work of Assadi, Sun and Weinstein [ASW19], which gave an MPC algorithm for generating random walks in an undirected regular graph. Finally, Łącki et al. [ŁMOS20] presented a new algorithm to compute random walks of length ` from every node in an arbitrary undirected graph. The algorithm of [ŁMOS20] still uses only O(log `) parallel rounds, and walks computed are now independent.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
From a high level perspective, the main idea behind all the MPC algorithms presented in [BCX11, ASW19, ŁMOS20] is to compute random walks of length ` by stitching together random walks of length /̀2 in a single parallel round. The walks of length /̀2 are computed by stitching together random walks of length /̀4 and so on. It is possible to prove that such strategy leads to algorithms that run in O(log `) parallel rounds as shown in previous work (this is also optimal under the 1-vs-2 cycle conjecture, as shown in [ŁMOS20]). Note that this technique in order to succeed computes in round i several random walks of length 2i for all the nodes in the network in parallel. This technique is very effective if we are interested in computing random walks from all the nodes in the graph, or, more precisely, when the number of walks computed out of a node is proportional to its stationary distribution. However, this approach leads to significant inefficiencies when we are interested in computing random walks only out of a subset of nodes or for a single node in the graph. This is even more important when we consider that in many applications as in clustering [GLMY11, GS12, WGD13] we are interested in running random walks only from a small subset of seed nodes. This leads to the natural question: Is it possible to compute efficiently and in parallel random walks only from a subset of nodes in a graph?
In this paper we answer this question in the affirmative, and we show an application of such a result in local clustering. Before describing our results in detail, we discuss the precise model of parallelism that we use in this work.
The MPC model. We design algorithms in the massively parallel computation (MPC) model, which is a theoretical abstraction of real-world system, such as MapReduce [DG08], Hadoop [Whi12], Spark [ZCF+10] and Dryad [IBY+07]. The MPC model [KSV10, GSZ11, BKS13] is the de-facto standard for analyzing algorithms for large-scale parallel computing.
Computation in MPC is divided into synchronous rounds over multiple machines. Each machine has memory S and at the beginning data is partitioned arbitrarily across machines. During each round, machines process data locally and then exchange data with the restriction that no machine receives more than S bits of data. The efficiency of an algorithm in this model is measured by the number of rounds it takes for the algorithm to terminate, by the size of the memory of every machine and by the total memory used in the computation. In this paper we focus on designing algorithm in the most restrictive and realistic regime where S ∈ O(nδ) for a small constant δ ∈ (0, 1) – these algorithms are called fully scalable.
Our contributions. Our first contribution is an efficient algorithm for computing multiple random walks from a single node in a graph efficiently. Theorem 1. There exists a fully scalable MPC algorithm that, given a graph G = (V,E) with n vertices and m edges, a root vertex r, and parameters B∗, ` and λ, can simulate B∗ independent random walks on G from r of length ` with an arbitrarily low error, in O(log ` logλB
∗) rounds and Õ(mλ`4 +B∗λ`) total space.
Our algorithm also applies to the more general problem of generating independent random walks from a subset of nodes in the graph: Theorem 2. There exists a fully scalable MPC algorithm that, given a graph G = (V,E) with n vertices and m edges and a collection of non-negative integer budgets (bu)u∈V for vertices in G such that ∑ u∈V bu = B
∗, parameters ` and λ, can simulate, for every u ∈ V , bu independent random walks on G of length ` from u with an arbitrarily low error, in O(log ` logλB
∗) rounds and Õ(mλ`4 +B∗λ`) total space. The generated walks are independent across starting vertices u ∈ V .
The following remark clarifies the effect of parameter λ on the number of machines. Remark 1. The parameter λ has nothing to do with the input of the algorithm, but is a trade-off parameter between space and round complexity. It is useful to think of it as λ = n for some (not necessarily a constant), in which case we get a round complexity of O(log ` logB∗/( log n)) ≤ O(log `/ ) and a total memory of Õ(mn `4 +B∗n `). We can set to, for example, 1/ log log n, to get nearly optimal total space and Õ(log `) round complexity.
If we compare our results with previous works, our algorithm computes truly independent random walks as [ŁMOS20] does. This is in contrast with the algorithm of [BCX11], which introduces dependent constructs not independent walks. Our algorithm has significantly better total memory than [ŁMOS20], which would result in memory Ω(m · B∗) for generating B∗ walks out of a root
node r. This comes at the cost of a slightly higher number of rounds (logλB ∗, a factor that in many applications can be considered constant).
The main idea is to preform multiple cycles of stitching algorithms, changing the initial distribution of the random walks adaptively. More precisely, in an initial cycle we construct only a few random walks, distributed according to the stationary distribution – this is known to be doable from previous work. Then, in each cycle we increase the number of walks that we build for node r by a factor of λ and we construct the walks only by activating other nodes in the graph that contribute actively in the construction of the random walks for r. In this way we obtain an algorithm that is significantly more work efficient in terms of total memory, compared with previous work.
Our second contribution is to present an application of our algorithm to estimating Personalized PageRank and to local graph clustering. To the best of our knowledge, our algorithm is the first local clustering algorithm that uses a number of parallel rounds that only have a logarithmic dependency on the length of the random walk used by the local clustering algorithm.
Theorem 3. For any λ > 1, let α ∈ (0, 1] be a constant and let C be a set satisfying that the conductance of C, Φ(C), is at most α/10 and Vol(C) ≤ 23Vol(G). Then there is an MPC algorithm for local clustering that usesO(log ` · logλB∗) = O(log log n · logλ(Vol(C))) rounds of communication and total memory Õ(mλ`4 + B∗λ`) = Õ(mλ+ λVol(C)2), where B∗ := 10
6 log3 n η2α2 , ` := 10 logn α and η = 110Vol(C) , and outputs a cluster with conductance O( √ α log(Vol(C))).
Finally we present an experimental analysis of our results where we show that our algorithm to compute random walk is significantly more efficient than previous work [ŁMOS20], and that our technique scale to very large graphs.
Additional related works. Efficient parallel random walks algorithm have also been presented in distributed computing [DSNPT13] and using multicores [SRFM16]. Although the algorithms in [DSNPT13] require a number of rounds linear in the diameter of the graph. The results in [SRFM16] are closer in spirit to our work here but from an algorithmic perspective the challenges in developing algorithms in multicore and MPC settings are quite different. In our setting, most of the difficulty is in the fact that there is no shared memory and coordination between machines. As a result, bounding communication between machines and number of rounds is the main focus of this line of research. From an experimental perspective an advantage of the MPC environment is that it can scale to larger instances, in contrast an advantage of the multicore approach is that it is usually faster in practice for medium size instances. In our case study, where one is interested in computing several local clusters from multiple nodes(for example to detect sybil attack(look at [ACE+14] and following work for applications) the MPC approach is often more suitable. This is due to the fact that the computation of multiple Personalized PageRank vectors at the same time often requires a large amount of space.
Several parallel algorithm have also been presented for estimating PageRank or PersonalizedPageRank [DSGP11, DSMPU15, BCX11], but those algorithms either have higher round complexity or introduce dependencies between the PersonalizedPageRank vectors computed for different nodes. Finally there has been also some work in parallelize local clustering algorithm [CS15], although all previously known methods have complexity linear in the number of steps executed by the random walk/process used in the algorithm(in fact, our method could potentially be used to speed-up this work as well).
Notation. We work on undirected unweighted graphs, which we usually denote by G = (V,E), where V and E are the set of vertices and the set of edges, respectively. We also have n := |V | and m := |E|, unless otherwise stated. We define matrix D as the diagonal matrix of degrees, i.e., Di,i = d(vi). Also, we let A be the adjacency matrix, where Ai,j = 1 if and only if there is an edge joining vi and vj , and Ai,j = 0, otherwise.

heading: 2 MPC random walks

text:
In this section, we present our main result to compute B∗ random walks from a single root vertex r up to a length of `, then we generalize it to multiple sources. As mentioned before, our main idea is to carefully stitch random walks adaptively to activate nodes locally. In the rest of the section, we start by presenting our main theorem and giving a brief overview of our algorithm. We then describe
our stitching algorithm, and analyze the number of random walks we must start from each node so that our algorithms work. Finally, we present the extension of our result to the setting with multiple sources. Theorem 1. There exists a fully scalable MPC algorithm that, given a graph G = (V,E) with n vertices and m edges, a root vertex r, and parameters B∗, ` and λ, can simulate B∗ independent random walks on G from r of length ` with an arbitrarily low error, in O(log ` logλB
∗) rounds and Õ(mλ`4 +B∗λ`) total space.

heading: 2.1 Overview of our Algorithm

text:
Here, we explain the frameworks of stitching and budgeting, which are the two key tools making up our algorithm. For simplicity and without loss of generality, we assume that each vertex v has its own machine that stores its neighborhood, its budgets, and its corresponding random walks.1
Remark 2. For ease of notation we assume that ` = 2j for some integer j. One can see that this assumption is without loss of generality, because otherwise one can always round ` to the smallest power of 2 greater than `, and solve the problem using the rounded `. This affects the complexity bounds by at most a constant factor.
Stitching.
Here, we explain the framework of stitching, which is a key tool for our algorithm. At each point in time the machine corresponding to v stores sets of random walks of certain lengths, each starting in v. Each edge of each walk is labeled by a number from 1 to `, denoting the position it will hold in the completed walk. Thus, a walk of length s could be labeled (k + 1, . . . , k + s) for some k. Initially each vertex generates a pre-determined number of random edges (or walks of length one) with each label from 1 to `. Thus at this point, we would find walks labeled 1, 2, 3, . . . on each machine. After the first round of stitching, these will be paired up into walks of length two, and so we will see walks labeled by (1, 2), (3, 4), (5, 6), . . . on each machine. After the second round of stitching we will see walks of length 4, such as (1, 2, 3, 4), and so on. Finally, after the last round of stitching, each machine will contain some number of walks of length ` (labeled from 1 to `), as desired.
At any given time let Wk(v) denote the set of walks stored by v whose first label is k and B(v, k) denotes their cardinality – in the future, we will refer to the function B as the budget. After the initial round of edge generation, Wk(v) consists of B(v, k) individual edges adjacent to v, for each v and k.
The rounds of communication proceed as follows: in the first round of stitching, for each edge (or length one walk) e in Wk(v), for any odd k, v sends a request for the continuation of the walk to z, where z is the other endpoint of e. That is, v sends a request to z for an element of Wk+1(z). Each vertex sends out all such requests simultaneously in a single MPC round. Following this each vertex replies to each request by sending a walk from the appropriate set. Crucially, each request must be answered with a different, independent walk. If the number of requests for Wk+1(z) exceeds |Wk+1(z)| = B(z, k + 1), the vertex z declares failure and the algorithm terminates. Otherwise, all such requests are satisfied simultaneously in a single MPC round. Finally, each vertex v increases the length of each of its walks in Wk(v) to two when k is odd, and deletes all remaining walks in Wk(v) when k is even (see Figure 1). For a more formal definition see Algorithm 1.
Budgeting.
A crucial aspect of stitching is that the budgets B(v, k) need to be carefully prescribed. If at any point in time a vertex receives more requests than it can serve, the entire algorithm fails. In the case where B(·, 1) follows the stationary distribution, this can be done (see for instance [ŁMOS20]), since the number of requests – at least in expectation – will also follow the stationary distribution. In our setting however, when B(·, 1) follows the indicator distribution of r, this is much more difficult. We should assign higher budgets to vertices nearer r; however, we have no knowledge of which vertices these are.
In other words, the main challenge in making stitching work with low space and few rounds is to set the vertex budgets (B(v, k)) accurately enough for stitching to succeed – this is the main technical contribution of this paper.
1In reality multiple vertices may share a machine, if they have low degree; or if a vertex has a high degree it may be accommodated by multiple machines in a constant-depth tree structure.
Our technique is to run multiple cycles of stitching sequentially. In the first cycle, we simply start from the stationary distribution. Then, with each cycle, we shift closer and closer to the desired distribution, in which the budget of r is much greater than the budgets of other vertices. We do this by augmenting B(r, 1) by some parameter λ each cycle. This forces us to augment other budgets as well: For example, for u in the neighborhood of r we expect to have a significantly increased budget B(u, 2). In order to estimate the demand on u (and all other vertices) we use data from the previous cycle.
Though initially only a few walks simulated by our algorithm start in r (we call these rooted walks), we are still able to derive some information from them. For instance, we can count the number of walks starting in r and reaching u as their second step. If κ rooted walks visited u as their second step in the previous cycle, we expect this number to increase to λ · κ in the following cycle. Hence, we can preemptively increase B(u, 2) to approximately λ · κ. More precisely, we set the initial budget of each vertex to ∼ B0 · deg(v) – an appropriately scaled version of the stationary distribution. This guarantees that the first round of stitching succeeds. Afterwards we set each budget B(v, k) individually based on the information gathered from the previous cycle. We first count the number of rooted walks that ended up at v as their kth step (Line 9). If this number is sufficiently large to be statistically significant (above some carefully chosen threshold θ in our case, Line 10), then we estimate the budget B(v, k) to be approximately λ ·κ in the following cycle (Line 11). On the other hand, if κ is deemed too small, it means that rooted random walks rarely reach v as their kth step, and the budget B(v, k) remains what it was before.
Algorithm 1 Stitching algorithm 1: procedure STITCH(G,B) 2: for v ∈ V in parallel do 3: for k ∈ [`] do 4: Wk(v)← a set of B(v, k) independent uniformly random edges adjacent on v 5: for j = 1 . . . log2 ` do . See Remark 2 6: for v ∈ V in parallel do 7: for k ≡ 1 (mod 2j) do 8: for walk p ∈Wk(v) do 9: send (v, k) to z, where z ← end vertex of p
10: for z ∈ V in parallel do 11: for each message (v, k) do 12: if Wk+2j−1(z) = ∅ then 13: return Fail 14: send (q, k, z) to v, where the walk q ← a randomly chosen walk in Wk+2j−1(v) 15: Wk+2j−1(z)←Wk+2j−1(z)\{q} 16: for v ∈ V in parallel do 17: for each message (q, k, z) do 18: p← any walk of length 2j−1 in Wk(v) with end vertex z 19: Wk(v)←Wk(v)\{p} ∪ {p+ q} . p+ q is the concatenated walk 20: for k ≡ 2j−1 (mod 2j) do 21: Wk(v)← ∅ 22: return W1(v) for all v ∈ V
Algorithm 2 Main Algorithm (Budgeting) 1: procedure MAIN(G, r, `, B∗, λ) 2: θ ← 10C`2 log n, B0 ← 30Cλ`3 log n, τ ← 1 + √ 20C logn
θ . Parameter settings 3: ∀v ∈ V , B0(v)← B0 · deg(v) 4: ∀v ∈ V, ∀k ∈ [`] : B(v, k)← B0 · deg(v) · τ3k−3 5: for i = 1 . . . blogλB∗c do 6: W1 ← STITCH(G,B) 7: W ←W1(r) 8: for v ∈ V, k ∈ [`] do 9: κ← |{w ∈W |wk = v}|
10: if κ ≥ θ then 11: B(v, k)← (B0(v) + λi · κ|W | ) · τ 3k−3 12: else 13: B(v, k)← B0(v) · τ3k−3 14: W1 ← STITCH(G,B) 15: W ←W1(r) 16: return W
Remark 3. In the above pseudocode (Algorithm 2) τ is a scaling parameter slightly greater then one. We augment all budgets B(·, k) by a factor τ3k−3, to insure that there are always slightly more walks with higher labels, and ensure that stitching succeeds with high probability.
Analysis. We are now ready to present the main properties of our algorithm: Lemma 1 (Correctness and complexity). Algorithm 2 takes O(log ` · logλB∗) rounds of MPC communication, STITCH terminates without failures with high probability, and the total amount of memory used for walks is ∑ v∈V ∑` k=1B(v, k) = O(mλ` 4 log n+B∗λ`).
Proof. Let P k be distribution of random walks of length k starting from r. That is, the probability that such a walk ends up in v ∈ V after k steps is P k(v).
To bound the round complexity we note that each call of STITCH takes only O(log `) rounds of communication, and it is called blogλB∗c+ 1 times; this dominates the round complexity. Further rounds of communication are needed to update the budgets. However this can be done in parallel for each vertex, and thus takes only one round per iteration of the outer for-loop (Line 5).
To prove that the algorithm fails with low probability, we must show the following crucial claim about the budgets B(V,K). Recall that the ideal budget in the ith iteration would be B(v, k) ≈ B0(v) + λ
i · P k(v). We show that in reality, the budgets do not deviate too much from this. Claim 1. After iteration i of the outer for-loop (Line 5) in Algorithm 2, with high probability B is set such that
∀v ∈ V, k ∈ [`] : B(v, k) ∈ [ (B0(v) + λ i · P k(v)) · τ3k−4, (B0(v) + λi · P k(v)) · τ3k−2 ] .
Proof. We first note how B(r, 1) – the budget of walks starting at the root vertex – evolves. For v = r and k = 1, κ is always equal to |W | – since r is the root vertex – and greater than θ. Therefore, B(r, 1) is set in Line 11 of Algorithm 2 to (B0(r) + λi). This is important, because it means that when setting other budgets in iteration i > 1, |W | is always (B0(r) + λi−1), the number of walks rooted from r in the previous round. The exception is the first iteration, when |W | is simply |B0(r)|. In both cases we may say that |W | ≥ λi−1. There are two options we have to consider: If after the ith round of STITCH κ exceeded θ, in which case our empirical estimator κ/|W | for P k(v) is deemed reliable. We then use this estimater to set the budget for the next round (see Line 11). Alternately, if κ did not exceed θ, the imperical estimator is deemed too unreliable; we then simply set B(v, k) proportionally to B0(v) (see Line 13).
Case I (κ < θ) then intuitively, κ is too small to provide an accurate estimator of P k(v). In this case we are forced to argue that the (predictable) term B0(v) dominates the (unknown) term P k(v). Since κ < θ, E(κ) = P k(v) · |W | ≤ 2θ. (The opposite happens with low probability2 by Chernoff bounds, since θ ≥ 10C log n.) Therefore,
B(v, k) = B0(v) · τ3k−3 ≤ (B0(v) + λi · P k(v)) · τ3k−2, and
B(v, k) = B0(v) · τ3k−3 = (B0(v) + λi · P k(v)) · ( 1− λ i · P k(v)
B0(v) + λi · P k(v)
) · τ3k−3.
So, we need to prove that 1− λ i·Pk(v) B0(v)+λi·Pk(v) ≥ τ −1. Now, by the above bound on E(κ) as well as the fact that |W | ≥ λi−1, we have 2θ ≥ P k(v) · |W | ≥ P k(v) ·λi−1, which results in λi ·P k(v) ≤ 2λθ. Consequently,(
1− λ i · P k(v)
B0(v) + λi · P k(v)
) ≥ (
1− 2λθ B0(v) + 2λθ
) ≥ τ−1.
Here we used B0(v) ≥ B0 ≥ 3λθ · ( √ θ/(20C log n)), which holds by definition of B0(v) and our setting of parameters B0 and θ.
Case II (κ ≥ θ), then intuitively κ is robust enough to provide a reliable estimator for P k(v). More precisely, κ/|W | ∈ [ E(κ/|W |) · τ−1,E(κ/|W |) · τ ] with high probability – indeed τ is defined in terms of θ deliberately in exactly such a way that this is guaranteed by Chernoff bounds. E(κ/|W |) = P k(v), therefore
λi · κ |W |
∈ [ λi · P k(v) · τ3k−4, λi · P k(v) · τ3k−2 ] ,
and
B(v, k) = (B0(v) + λ i · κ |W | ) · τ3k−3
∈ [( B0(v) + λ i · P k(v) ) · τ3k−4, (B0(v) + λi · P k(v)) · τ3k−2 ] .
2Throughout the proof, we say ’low probability’ to mean probability of n−Ω(C) where C can be set arbitrarily high.
STITCH only reports failure if for some v ∈ V and k ∈ [2, `], vertex v receives more requests for walks in Wk(v) than |Wk(v)| = B(v, k). Number of such request is upper bounded by the number of edges ending in v generated by neighbors of v, say w at level k−1. That is, the number of requests for Wk(v) is in expectation at most∑
w∈Γ(v)
1
d(w) B(w, k − 1) ≤ ∑ w∈Γ(v) 1 d(w) (B0(w) + λ i · P k−1(w)) · τ3k−5
=  ∑ w∈Γ(v) 1 d(w) B0(w) + λ i · P k(v)  · τ3k−5 = ( d(v)B0 + λ i · P k(v) ) · τ3k−5
= ( B0(v) + λ i · P k(v) ) · τ3k−5.
Since this is greater than θ, the actual number of requests is at most (B0(v) + λi · P k(v)) · τ3k−4 ≤ B(v, k), with high probability by Chernoff. Therefore, STITCH indeed does not fail.
Finally, we prove the memory bound. By setting of parameter θ, τ3k−2 is at most a constant. Also by the setting of parameters we have,
B(v, k) ≤ (B0(v) + λblogλ B ∗c+1 · P k(v)) · τ3k−2 = O(B0(v) +B∗λ · P k(v)),
and ∑ v∈V ∑` k=1B(v, k) = O(1) · ∑ v∈V ∑` k=1(B0(v)+B ∗λ ·P k(v)) = O(mλ`4 log n+B∗λ`).
This gives us the proof of Theorem 1. Now, one can easily extend this result to the case when multiple sources for the starting vertex is considered. Theorem 2 is proven in Appendix A.

heading: 3 From random walks to local clustering

text:
We now present two applications for our algorithm to compute random walks. In particular we show how to use it to compute PageRank vectors and how to use it to compute local clustering. In interest of space, we only state here our main results and we defer all the technical definition and proofs to the Appendix.
Approximating PageRank using MPC random walks Interestingly, we can show that we can use our algorithm as a primitive to compute PersonalizedPageRank for a node of for any input vector3
Theorem 4 (Approximating PersonalizedPageRank using MPC random walks). For any starting single vertex vector s (indicator vector), any α ∈ (0, 1) and any η, there is a MPC algorithm that usingO(log `·logλB∗) rounds of communication and the total amount of memory ofO(mλ`4 log n+ B∗λ`), outputs a vector q̃, such that q̃ is a η-additive approximation to prα(s), whereB ∗ := 10 6 log3 n η2α2 and ` := 10 lognα .
The proof is deferred to Appendix C.
Using approximate PersonalizedPageRank vectors to find sparse cuts Now we can use the previous result on PersonalizedPageRank to find sets with relatively sparse cuts. Roughly speaking, we argue that for any set C of conductance O(α), for many vertices v ∈ C, if we calculate an approximate prα(v) using our algorithms and perform a sweep cut over it, we can find a set of conductance O( √ α log(Vol(C))). This result is stated in Theorem 3. The proof of this result is very similar to the proofs of Section 5 of [ACL06], however since our approximation guarantees are slightly different, we need to modify some parts of the proof for completeness. The full proof is presented in Appendix Our main result of this subsection is stated below.
3We note that [BCX11] also propose an algorithm to compute PersonalizedPageRank vector but with the limitation that this could be computed only for a single node and not for a vector.

heading: 4 Empirical Evaluation

text:
In this Section we present empirical evaluations of our algorithms for random walk generation, as well as clustering. As our datasets, we use several real-world graphs form the Stanford Network Analysis Project [LK14, LKF07, LLDM08, KY04, YL12]. The graphs are undirected and come mostly (though not in all cases) from ’Networks with ground truth communities’, where the clustering application is most relevant. In order to demonstrate the scalability of our main Algorithm we use graphs of varying sizes, as demonstrated in the table below.
The experiments were performed on Amazon’s Elastic Map-Reduce system using the Apache Hadoop library. The clusters consisted of 30 machines, each of modest memory and computing power (Amazon’s m4.large instance) so as to best adhere to the MPC setting. Each experiment described in this section was repeated 3 times to minimize the variance in performance inherent in distributed systems like this.
Practical considerations. In Section 2 we worked with the guarantee that no walks "fail" in the stitching phase. This assumption can be fulfilled at nearly no expense to the (asymptotic) guarantees in space and round complexity of Theorem 1, and make the proof much cleaner. In practice however, allowing some small fraction of the walks to fail allows for a more relaxed setting of the parameters, and thus better performance.
Each experiment is performed with 15 roots, selected uniformly at random. The main parameters defining the algorithm are as follows: ` — the length of a target random walk (16 and 32 in various experiments), C — The number of cycles (iterations of the for-loop in Line 5 of Algorithm 2) performed. , B0 — the initial budget-per-degree of each vertex, λ — the approximate scaling of the budgets of the root vertices each cycle, τ — a parameter defining the amount of excess budget used in stitching. This is used somewhat differently here than in Algorithm 2. For more details see Appendix D.

heading: 4.1 Scalability

text:
In this section we present the results of our experiments locally generating random walks simultaneously from multiple root vertices. We use the graphs COM-DBLP, COM-YOUTUBE, COMLIVEJOURNAL, and COM-ORKUT, in order to observe how the runtime of Algorithm 2 scales with the size of the input graph. In each of these graphs, 15 root vertices have been randomly chosen. We ran three experiments with various settings of the parameters, of which one is presented below. For additional experiments see Appendix D. B0 is set to be proportional to n/m – that is inverse proportional to the average degree – since the initial budget of each of each vertex is set to B0 times the degree of the vertex.
We report the execution time in the Amazon Elastic Map-Reduce cluster, as well as the number of rooted walks generated. Finally, under ’Walk failure rate’, we report the percentage of rooted walks that failed in the last cycle of stitching. This is the crucial quantity; earlier cycles are used only to calibrate the vertex budgets for the final cycle.4
4For completeness, we report the empirical standard deviation everywhere. Note, however, that due to the high resource requirement of these experiments, each one was only repeated three times.
We observe that Algorithm 2 successfully generates a large number of rooted walks – far more than the initial budgets of the root vertices. As predicted, execution time scales highly sublinearly with the size of the input (recall for example, that COM-ORKUT is more than a hundred times larger than COM-DBLP). The failure rate of walks decreases with the size of the graph in this dataset, with that of COM-ORKUT reaching as low as 3.4% on average; this may be due to the the higher average degree of our larger graphs, leading to the random walks spreading out more.
In Appendix D we report the results of two more experiments, including one with longer walks.

heading: 4.2 Comparison

text:
In this section we compare to the previous work of [ŁMOS20] for generating random walks in the MPC model. This work heavily relies upon generating random walks from all vertices simultaneously, with the number of walks starting from a given vertex v being proportional d(v). In many applications, however, we are interested in computing random walks from a small number of root vertices. The only way to implement this using the methods of [ŁMOS20] is to start with an initial budget large enough to guarantee the desired number of walks from each vertex.
We perform similar experiments to those in the previous section – albeit on much smaller graphs. Each graph has 15 root vertices chosen randomly, from which we wish to sample random walks. In Table 3 we set B0 to 1 and perform C = 3 cycles of Algorithm 2 with λ = 10, effectively augmenting the budget of root vertices by a factor 100 by the last cycle. Correspondingly, we implement the algorithm of [ŁMOS20] – which we call UNIFORM STITCHING – by simply setting the initial budget 100 times higher, and performing only a single cycle of stitching, ie.: B0 = 100, C = 1.
We observe that the running time of our Algorithm 2 barely differs across the three graphs, despite the nearly 100-factor difference between the sizes of CA-GRQC and COM-DBLP. At this small size, the execution time is dominated by setting up the 12 Map-Reduce rounds required to execute Algorithm 2 with these parameters. As expected, the baseline far outperforms our algorithm on the smallest graph; as the size of the input graph grows, however, its running time deteriorates quickly. For larger setting of the parameters UNIFORMSTITCHING can no longer complete on the cluster, due to the higher memory requirement, as we can see in Table 6 in Appendix D.
Acknowledgments and Disclosure of Funding
This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 759471).

title: Efficient and Local Parallel Random Walks


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

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

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

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

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

authors:
Martin J Wainwright
Michael Irwin Jordan

title: Graphical models, exponential families, and variational inference

venue: Now Publishers Inc,

authors:
H Chau Nguyen
Riccardo Zecchina
Johannes Berg

title: Inverse statistical problems: from the inverse ising problem to data science

venue: Advances in Physics,

authors:
Erik Aurell
Magnus Ekeberg

title: Inverse ising inference using all the data

venue: Physical review letters,

authors:
Ludovica Bachschmid-Romano
Manfred Opper

title: Learning of couplings for random asymmetric kinetic ising models revisited: random correlation matrices and learning curves

venue: Journal of Statistical Mechanics: Theory and Experiment,

title: Statistical mechanics of the inverse ising problem and the optimal objective function

title: A statistical physics approach to learning curves for the inverse ising problem

authors:
Alia Abbara
Yoshiyuki Kabashima
Tomoyuki Obuchi
Yingying Xu

title: Learning performance in inverse ising problems with sparse teacher couplings

authors:
Martin J Wainwright
John D Lafferty
Pradeep K Ravikumar

title: High-dimensional graphical model selection using `_1-regularized logistic regression

venue: In Advances in neural information processing systems,

authors:
Holger Höfling
Robert Tibshirani

title: Estimation of sparse binary pairwise markov networks using pseudo-likelihoods

venue: Journal of Machine Learning Research,

authors:
Pradeep Ravikumar
Martin J Wainwright
John D Lafferty

title: High-dimensional ising model selection using `1-regularized logistic regression

venue: The Annals of Statistics,

authors:
Narayana P Santhanam
Martin J Wainwright

title: Information-theoretic limits of selecting binary graphical models in high dimensions

venue: IEEE Transactions on Information Theory,

authors:
Aurélien Decelle
Federico Ricci-Tersenghi

title: Pseudolikelihood decimation algorithm improving the inference of the interaction network in a general class of ising models

title: Efficiently learning ising models on arbitrary graphs

venue: In Proceedings of the forty-seventh annual ACM symposium on Theory of computing,

authors:
Marc Vuffray
Sidhant Misra
Andrey Lokhov
Michael Chertkov

title: Interaction screening: Efficient and sample-optimal learning of ising models

venue: In Advances in Neural Information Processing Systems,

authors:
Andrey Y Lokhov
Marc Vuffray
Sidhant Misra
Michael Chertkov

title: Optimal structure and parameter learning of ising models

venue: Science advances,

authors:
Shanshan Wu
Sujay Sanghavi
Alexandros G Dimakis

title: Sparse logistic regression learns all discrete pairwise graphical models

venue: Advances in Neural Information Processing Systems,

authors:
Guy Bresler
Elchanan Mossel
Allan Sly

title: Reconstruction of markov random fields from samples: Some observations and algorithms. In Approximation, Randomization and Combinatorial Optimization

venue: Algorithms and Techniques,

authors:
Adarsh Prasad
Vishwak Srinivasan
Sivaraman Balakrishnan
Pradeep Ravikumar

title: On learning ising models under huber’s contamination model

authors:
Sahand N Negahban
Pradeep Ravikumar
Martin J Wainwright
Bin Yu

title: A unified framework for high-dimensional analysis of m-estimators with decomposable regularizers

venue: Statistical science,

authors:
Marc Mézard
Giorgio Parisi
Miguel Angel Virasoro

title: Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, volume 9

venue: World Scientific Publishing Company,

title: Statistical physics of spin glasses and information processing: an introduction

authors:
Manfred Opper
David Saad

title: Advanced mean field methods: Theory and practice

venue: MIT press,

authors:
Marc Mezard
Andrea Montanari

title: Information, physics, and computation

title: Regression shrinkage and selection via the lasso

venue: Journal of the Royal Statistical Society: Series B (Methodological),

authors:
Mohsen Bayati
Andrea Montanari

title: The lasso risk for gaussian matrices

authors:
Sundeep Rangan
Alyson K Fletcher
Vivek K Goyal

title: Asymptotic analysis of map estimation via the replica method and applications to compressed sensing

authors:
Christos Thrampoulidis
Ehsan Abbasi
Babak Hassibi

title: Lasso with non-linear measurements is equivalent to one with linear measurements

authors:
Cédric Gerbelot
Alia Abbara
Florent Krzakala

title: Asymptotic errors for convex penalized linear regression beyond gaussian matrices

venue: arXiv preprint arXiv:2002.04372,

authors:
Xiangming Meng
Tomoyuki Obuchi
Yoshiyuki Kabashima

title: Structure learning in inverse ising problems using `2-regularized linear estimator

title: Statistical analysis of non-lattice data

venue: Journal of the Royal Statistical Society: Series D (The Statistician),

title: A generalized linear model with gaussian regressor variables. In A Festschrift for Erich L

authors:
Yaniv Plan
Roman Vershynin

title: The generalized lasso with non-linear observations

venue: IEEE Transactions on information theory,

authors:
Yue Zhang
Weihong Guo
Soumya Ray

title: On the consistency of feature selection with lasso for non-linear targets

venue: In International Conference on Machine Learning,

title: High-dimensional estimation of structured signals from non-linear observations with general convex loss functions

authors:
Nicolai Meinshausen
Peter Bühlmann

title: High-dimensional graphs and variable selection with the lasso

venue: The annals of statistics,

authors:
Peng Zhao
Bin Yu

title: On model selection consistency of lasso

venue: The Journal of Machine Learning Research,

title: A theory of the learnable

venue: Communications of the ACM,

authors:
Andreas Engel
Christian Van den Broeck

title: Statistical mechanics of learning

authors:
Federica Gerace
Bruno Loureiro
Florent Krzakala
Marc Mézard
Lenka Zdeborová

title: Generalisation error in learning with random features and the hidden manifold model

authors:
Madhu Advani
Surya Ganguli

title: Statistical mechanics of optimal convex inference in high dimensions

venue: Physical Review X,

authors:
Persi Diaconis
Mehrdad Shahshahani

title: On the eigenvalues of random matrices

venue: Journal of Applied Probability,

title: On random matrices from the compact classical groups

venue: Annals of mathematics,

authors:
José Bento
Andrea Montanari

title: Which graphical models are difficult to learn

venue: arXiv preprint arXiv:0910.5761,

authors:
Nicholas Metropolis
Arianna W Rosenbluth
Marshall N Rosenbluth
Augusta H Teller
Edward Teller

title: Equation of state calculations by fast computing machines

venue: The journal of chemical physics,

title: Monte carlo sampling methods using markov chains and their applications

authors:
Stuart Geman
Donald Geman

title: Stochastic relaxation, gibbs distributions, and the bayesian restoration of images

venue: IEEE Transactions on pattern analysis and machine intelligence,

title: The bethe approximation for solving the inverse ising problem: a comparison with other inference methods

authors:
H Chau Nguyen
Johannes Berg

title: Bethe–peierls approximation and the inverse ising problem

heading: 1 Introduction

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

heading: 1.1 Contributions

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

heading: 1.2 Related works

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

heading: 2 Background and Problem Setup

heading: 2.1 Ising Model

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

heading: 3 Statistical Mechanics Analysis

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

heading: 3.1 Problem Formulation

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

heading: 3.2 Replica computation of the free energy density

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

heading: 3.3 High-dimensional asymptotic result

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

heading: 4 Discussions

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

heading: 5 Experimental results

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

heading: 6 Conclusion

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

heading: Acknowledgements

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

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