Submatrix localization via message passing

Introduction

The problem of submatrix detection and localization, also known as noisy biclustering , deals with finding a submatrix with an elevated mean in a large noisy matrix, which arises in many applications such as social network analysis and gene expression data analysis. A widely studied statistical model is the following:

where $\mu>0$, $\mathbf{1}_{C_{1}^{*}}$ and $\mathbf{1}_{C_{2}^{*}}$ are indicator vectors of the row and column support sets $C_{1}^{*}\subset[n_{1}]$ and $C_{2}^{*}\subset[n_{2}]$ of cardinality $K_{1}$ and $K_{2}$, respectively, and $Z$ is an $n_{1}\times n_{2}$ matrix consisting of independent standard normal entries. The objective is to accurately locate the submatrix by estimating the row and column support based on the large matrix $W$.

For simplicity we start by considering the symmetric version of this problem, namely, locating a principal submatrix, and later extend our theoretic and algorithmic findings to the asymmetric case. To this end, consider

where $C^{*}\subset[n]$ has cardinality $K$ and $Z$ is an $n\times n$ symmetric matrix with $\{Z_{ij}\}_{1\leq i\leq j\leq n}$ being independent standard normal. Given the data matrix $W$, the problem of interest is to recover $C^{*}$. This problem has been investigated in as a prototypical example of the hidden community problem,A slight variation of the model in is that the data matrix therein is assumed to have zero diagonal. As shown in , the absence of the diagonal has no impact on the statistical limit of the problem as long as $K\to\infty$, which is the case considered in the present paper. because the distribution of the entries exhibits a community structure, namely, $W_{i,j}\sim{\mathcal{N}}(\mu,1)$ if both $i$ and $j$ belong to $C^{*}$ and $W_{i,j}\sim{\mathcal{N}}(0,1)$ if otherwise.

Assuming that $C^{*}$ is drawn from all subsets of $[n]$ of cardinality $K$ uniformly at random, we focus on the following two types of recovery guarantees.Exact and weak recovery are called strong consistency and weak consistency in , respectively. Let $\xi=\mathbf{1}_{C^{*}}\in\{0,1\}^{n}$ denote the indicator of the community. Let $\widehat{\xi}=\widehat{\xi}(A)\in\{0,1\}^{n}$ be an estimator.

We say that $\widehat{\xi}$ weakly recovers $\xi$ if, as $n\to\infty$, $\ d(\xi,\widehat{\xi})/K\to 0$ in probability, where $d$ denotes the Hamming distance.

Intuitively, for a fixed matrix size $n$, as either the submatrix size $K$ or the signal strength $\mu$ decreases, it becomes more difficult to locate the submatrix. A key role is played by the parameter

which is the signal-to-noise ratio for classifying an index $i$ according to the statistic $\sum_{j}W_{i,j}$, which is distributed according to ${\mathcal{N}}(\mu K,n)$ if $i\in C^{*}$ and ${\mathcal{N}}(0,n)$ if $i\not\in C^{*}.$ As shown in Appendix A, it turns out that if the submatrix size $K$ grows linearly with $n$, the information-theoretic limits of both weak and exact recovery are easily attainable via thresholding. To see this, note that in the case of $K\asymp n$ simply thresholding the row sums can provide weak recovery in $O(n^{2})$ time provided that $\lambda\to\infty$, which coincides with the information-theoretic conditions of weak recovery as proved in . Moreover, in this case, one can show that this thresholding algorithm followed by a linear-time voting procedure achieves exact recovery whenever information-theoretically possible. Thus, this paper concentrates on weak and exact recovery in the sublinear regime of

We show that an optimized message passing algorithm provides weak recovery in nearly linear – $O(n^{2}\log n)$ – time if $\lambda>1/e$. This extends the sufficient conditions obtained in for the regime $K=\Theta(\sqrt{n}).$ Our algorithm is the same as the message passing algorithm proposed in , except that we find the polynomial that maximizes the signal-to-noise ratio via Hermite polynomials instead of using the truncated Taylor series as in . The proofs follow closely those in , with the most essential differences described in Remark 6. We observe that $\lambda>1/e$ is much more stringent than $\lambda>\frac{4K}{n}\log\frac{n}{K}$, the information-theoretic weak recovery threshold established in . It is an open problem whether any polynomial-time algorithm can provide weak recovery for $\lambda\leq 1/e$. In addition, we show that if $\lambda>1/e$, the message passing algorithm followed by a linear-time voting procedure can provide exact recovery whenever information theoretically possible. This procedure achieves the optimal exact recovery threshold determined in if $K\geq(\frac{1}{8e}+o(1))\frac{n}{\log n}$. See Section 1.2 for a detailed comparison with information-theoretic limits.

The message passing algorithm is simpler to formulate and analyze for the principal submatrix recovery problem; nevertheless, we show in Section 4 how to adapt the message passing algorithm and its analysis to the biclustering problem. Butucea et al. obtained sharp conditions for exact recovery for the bicluster problem. We show that calculations in with minor adjustments provide information theoretic conditions for weak recovery as well. The connection between weak and exact recovery via the voting procedure described in carries over to the biclustering problem.

1 Algorithms and main results

To avoid a plethora of factors $\frac{1}{\sqrt{n}}$ in the notation, we describe the message-passing algorithm using the scaled version $A=\frac{1}{\sqrt{n}}W.$ The entries of $A$ have variance $\frac{1}{n}$ and mean or $\frac{\mu}{\sqrt{n}}.$ This section presents algorithms and theoretical guarantees for the symmetric model (2). Section 4.2 gives adaptations to the asymmetric case for the biclustering problem (1).

with the initial conditions $\theta^{0}_{i\to j}\equiv 0.$ Moreover, let $\theta^{t+1}_{i}$ denote index $i$’s belief at iteration $t+1$, which is given by

The form of (4) is inspired by belief propagation algorithms, which have the natural non backtracking property: the message sent from $i$ to $j$ at time $t+1$ does not depend on the message sent from $j$ to $i$ at time $t,$ thereby reducing the effect of echoes of messages sent by $j.$

Suppose as $n\to\infty$, the messages $\theta^{t}_{i}$ (for fixed $t$) are such that the empirical distributions of $(\theta^{t}_{i}:i\in[n]\backslash C^{*})$ and $(\theta^{t}_{i}:i\in C^{*})$ converge to Gaussian distributions with a certain mean and variance. Specifically, $\theta^{t}_{i}$ is approximately ${\mathcal{N}}(\mu_{t},\tau_{t})$ for $i\in C^{\ast}$ and ${\mathcal{N}}(0,\tau_{t})$ for $i\notin C^{\ast}$. Then (4), (5), and the fact $\theta_{i\to j}^{t}\approx\theta_{i}^{t}$ for all $i,j$ suggest the following recursive equations for $t\geq 0$:

where $Z$ represents a standard normal random variable, and the initial conditions are $\mu_{0}=\tau_{0}=0.$ Following , we call (6) and (7) the state evolution equations, which are justified in Section 2. Thus, it is reasonable to estimate $C^{*}$ by selecting those indices $i$ such that $\theta_{i}^{t}$ exceeds a given threshold.

Suppose, for the time being, that message distributions are Gaussian with parameters accurately tracked by the state evolution equations. Then classifying an index $i$ based on $\theta_{i}^{t}$ boils down to testing two Gaussian hypotheses with signal-to-noise ratio $\frac{\mu_{t+1}^{2}}{\tau_{t+1}}.$ This gives guidance for selecting the functions $f(\cdot,t)$ based on $\mu_{t}$ and $\tau_{t}$ to maximize $\frac{\mu_{t+1}}{\sqrt{\tau_{t+1}}}$. For $t=0$ any choice of $f$ is equivalent, so long as $f(0,0)>0.$ Without loss of generality, for $t\geq 1,$ we can assume that the variances are normalized, namely, $\tau_{t}=1$ (e.g. we take $f(0,0)=1$ to make $\tau_{1}=1$) and choose $f(\cdot,t)$ to be the maximizer of

Clearly, the best $g$ aligns with $\rho$ and we obtain

With this optimized $f$, we have $\tau_{t}\equiv 1$ and the state evolution (6) reduces to

Therefore if $\lambda>1/e$, then (11) has no fixed point and hence $\mu_{t}\rightarrow\infty$ as $t\rightarrow\infty$.

Directly carrying out the above heuristic program, however, seems challenging. To rigorously justify the state evolution equations in Section 2 we rely on the the method of moments, requiring $f$ to be a polynomial, which prompts us to look for the best polynomial of a given degree that maximizes the signal-to-noise ratio. Denoting the corresponding state evolution by $(\widehat{\mu}_{t},\widehat{\tau}_{t})$, we aim to solve the following finite-degree version of (8):

As shown in Lemma 7, this problem can be easily solved via Hermite polynomials, which form an orthogonal basis with respect to the Gaussian measure, and the optimal choice, say, $f_{d}(\cdot,t)$ can be obtained by normalizing the first $d+1$ terms in the orthogonal expansion of relative density (9), i.e., the best degree-$d$ $L_{2}$-approximation. Compared to [11, Lemma 2.3] which shows the existence of a good choice of polynomial that approximates the ideal state evolution (11) based on Taylor expansions, our approach is to find the best message-passing rule of a given degree which results in the following state evolution that is optimal among all $f$ of degree $d$:

For any $\lambda>1/e$, there is an explicit choice of the degree $d$ depending only on $\lambda$,As $\lambda$ gets closer to the critical value $1/e$, we need a higher degree to ensure (13) diverges and in fact $d$ grows quite slowly as $\Theta(\log\frac{1}{\lambda e-1}/\log\log\frac{1}{\lambda e-1})$ See Remark 40. so that $\widehat{\mu}_{t}\to\infty$ as $t\to\infty$ and the state evolution (13) for fixed $t$ correctly predicts the asymptotic behavior of the messages when $n\to\infty$. As discussed above, $\widetilde{C}$ produced by thresholding messages $\theta_{i}^{t}$, is likely to contain a large portion of $C^{\ast}$, but since $K=o(n)$, it may (and most likely will) also contain a large number of indices not in $C^{\ast}$. Following [11, Lemma 2.4], we show that the power iterationNote that as far as statistical utility is concerned, we could replace $\widehat{u}$ produced by the power iteration by the leading singular vector of $A_{\widetilde{C}}$, but that would incur a higher time complexity because singular value decomposition in general takes $O(n^{3})$ time to compute. (a standard spectral method) in Algorithm 1 can remove a large portion of the outlier vertices in $\widetilde{C}.$

Combining message passing plus spectral cleanup yields the following algorithm for estimating $C^{\ast}$ based on the messages $\theta_{i}^{t}$.

The following theorem provides a performance guarantee for Algorithm 1 to approximately recover $C^{\ast}.$

Fix $\lambda>1/e.$ Let $K$ and $\mu$ depend on $n$ in such a way that $\frac{\mu^{2}K^{2}}{n}\to\lambda$ and $\Omega(\sqrt{n})\leq K\leq o(n)$ as $n\to\infty.$ Consider the model (2) with $|C^{*}|/K\to 1$ in probability as $n\to\infty.$ Define $d^{\ast}(\lambda)$ as in (40). For every $\eta\in(0,1),$ there exist explicit positive constants $t^{\ast},s^{\ast},c$ depending on $\lambda$ and $\eta$ such that Algorithm 1 returns $|\widehat{C}\Delta C^{\ast}|\leq\eta K$, with probability converging to one as $n\to\infty$, and the total time complexity is bounded by $c(\eta,\lambda)n^{2}\log n$, where $c(\eta,\lambda)\to\infty$ as either $\eta\to 0$ or $\lambda\to 1/e.$

After the message passing algorithm and spectral cleanup are applied in Algorithm 1, a final linear-time voting procedure is deployed to obtain weak or exact recovery, leading to Algorithm 2 next. As in , we consider a threshold estimator for each vertex $i$ based on a sum over $\widehat{C}$ given by $r_{i}=\sum_{j\in\widehat{C}}A_{ij}$. Intuitively, $r_{i}$ can be viewed as the aggregated “votes” received by the index $i$ in $\widehat{C}$, and the algorithm picks the set of $K$ indices with the most significant “votes”. To show that this voting procedure succeeds in weak recovery, a key step is to prove that $r_{i}$ is close to $\sum_{j\in C^{\ast}}A_{ij}.$ If $\mu=\Theta(1)$ as in , given that $|\widehat{C}\triangle C^{*}|=o(K),$ the error incurred by summing over $\widehat{C}$ instead of over $C^{*}$ could be bounded by truncating $A_{ij}$ to a large magnitude. However, for $\mu\to 0$ that approach fails (see Remark 6 for more details). Our approach is to introduce the clean-up procedure in Algorithm 2 based on the successive withholding method described in (see also for variants of this method). In particular, we randomly partition the set of vertices into $1/\delta$ subsets. One at a time, one subset, say $S$, is withheld to produce a reduced set of vertices $S^{c}$, on which we apply Algorithm 1. The estimate obtained from $S^{c}$ is then used by the voting procedure to classify the vertices in $S$. The analysis of the two stages is decoupled because conditioned on $C^{*}$, the outcome of Algorithm 1 depends only on $A_{S^{c}}$, which is independent of $A_{SS^{c}}$ used in the voting.

The following theorem provides a sufficient condition for the message passing plus voting cleanup procedure (Algorithm 2) to achieve weak recovery, and, if the information-theoretic sufficient condition is also satisfied, exact recovery.

Let $K$ and $\mu$ depend on $n$ in such a way that $\frac{\mu^{2}K^{2}}{n}\to\lambda$ for some fixed $\lambda>1/e$. and $\Omega(\sqrt{n})\leq K\leq o(n)$ as $n\to\infty.$ Consider the model (2) with $|C^{*}|\equiv K$. Let $\delta>0$ be such that $\lambda e(1-\delta)>1$. Define $d^{*}=d^{\ast}(\lambda(1-\delta))$ as per (40). Then there exist positive constants $t^{\ast},s^{\ast},c$ determined explicitly by $\delta$ and $\lambda$, such that

(Weak recovery) Algorithm 2 returns $C^{\prime}$ with $|C^{\prime}\Delta C^{\ast}|/K\to 0$ in probability as $n\to\infty$.

(Exact recovery) Furthermore, assume that

Let $\delta>0$ be chosen such that for all sufficiently large $n,$

The total time complexity is bounded by $c(\delta,\lambda)n^{2}\log n$, where $c(\delta,\lambda)\to\infty$ as $\delta\to 0$ or $\lambda\to 1/e$.

As shown in [13, Theorem 7], if there is an algorithm that can approximately recover $|C^{*}|$ even if $|C^{*}|$ is random and only approximately equal to $K,$ then that algorithm can be combined with a linear-time voting procedure to achieve exact recovery. By Theorem 1, Algorithm 1 indeed works for such random $|C^{*}|$ and so the second part of Theorem 2 follows from Theorem 1 and the general results of .

Theorem 2 ensures Algorithm 2 achieves exact recovery if both (14) and $\lambda>1/e$ hold; it is of interest to compare these two conditions. Note that

Hence, if $\liminf_{n\to\infty}K\log n/n\geq\frac{1}{8e}$, (14) implies $\lambda>1/e$ and thus (14) alone is sufficient for Algorithm 2 to succeed; if $\limsup_{n\to\infty}K\log n/n\leq\frac{1}{8e}$, then $\lambda>1/e$ implies (14) and thus $\lambda>1/e$ alone is sufficient for Algorithm 2 to succeed. The asymptotic regime considered in entails $K=\Theta(\sqrt{n}),$ in which case the condition $\lambda>1/e$ is sufficient for exact recovery.

2 Comparison with information theoretic limits

As noted in the introduction, in the regime $K=\Theta(n)$, a thresholding algorithm based on row sums provides weak and, if a voting procedure is also used, exact recovery whenever it is informationally possible. In this subsection, we compare the performance of the message passing algorithms to the information-theoretic limits on the recovery problem in the regime (3). Notice that the comparison here takes into account the sharp constant factors. Information-theoretic limits for the biclustering problem are discussed in Section 4.1.

The information-theoretic threshold for weak recovery has been determined in [13, Theorem 2], which, in the regime of (3), boils down to the following: Weak recovery is possible if

This implies that the minimal signal-to-noise ratio for weak recovery is

for any $\epsilon>0$, which vanishes in the sublinear regime of $K=o(n)$. In contrast, in the regime (3), to achieve weak recovery message passing (Algorithm 1) demands a non-vanishing signal-to-noise ratio, namely, $\lambda>1/e$. No polynomial-time algorithm is known to succeed if $\lambda\leq 1/e$, suggesting that computational complexity might incur a severe penalty on the statistical optimality when $K=o(n)$.

Exact recovery

In the regime of (3), the information limits of exact recovery (see [13, Theorem 4 and Remark 7]) are as follows: Exact recovery is possible if (14) holds, and impossible if

In view of Remark 2, we conclude that Algorithm 2 achieves the sharp threshold of exact recovery if

We note that a counterpart of this conclusion for the biclustering problem is obtained in Remark 7 in terms of the submatrix sizes.

To further the discussion on weak and exact recovery, consider the regime

where $s\geq 1$, $\rho\in(0,1)$, and $\mu_{0}>0$ are fixed constants. Throughout this regime, weak recovery is information theoretically possible because the left-hand side of (15) is $\Omega(\frac{\log n}{\log\log n})\to\infty$. On one hand, in view of (14) and (17), exact recovery is possible if $\frac{\rho\mu_{0}^{2}}{8}>1$ and impossible if $\frac{\rho\mu_{0}^{2}}{8}<1$. On the other hand, $\lambda=\rho^{2}\mu_{0}^{2}(\log n)^{2-s},$ yielding:

When $1\leq s<2$, then $\lambda=\Omega(\log^{2-s}n)\to\infty.$ Thus weak recovery is achievable in polynomial-time by the message passing algorithm, spectral methods, or even row-wise thresholding. If $\frac{\rho\mu_{0}^{2}}{8}>1$, exact recovery is attainable in polynomial-time by combining the weak recovery algorithm with a linear time voting procedure as shown in .

When $s=2$, then $\lambda=\rho^{2}\mu_{0}^{2},$ and weak recovery by the message passing algorithm is possible if $\rho^{2}\mu_{0}^{2}e>1.$ Fig. 1 shows the curve $\{(\mu_{0},\rho):\rho^{2}\mu_{0}^{2}e=1\}$ corresponding to the weak recovery condition by the message passing algorithm, and the curve $\{(\mu_{0},\rho):\rho\mu_{0}^{2}/8=1\}$ corresponding to the information-theoretic exact recovery condition. When $\rho\geq\frac{1}{8e}$, the latter curve dominates the former and Algorithm 2 achieves optimal exact recovery.

When $s>2$, $\lambda\to 0$, and no polynomial-time procedure is known to provide weak, let alone exact, recovery.

3 Comparison with the spectral limit

It is reasonable to conjecture that $\lambda>1$ is the spectral limit for recoverability by spectral estimation methods. This conjecture is rather vague, because it is difficult to define what constitutes spectral methods. Nevertheless, some evidence for this conjecture is provided by [11, Proposition 1.1], which, in turn, is based on results on the spectrum of a random matrix perturbed by adding a rank-one deterministic matrix [17, Theorem 2.7].

For the submatrix detection problem, namely, testing $\mu=0$ (pure noise) versus $\mu>0$, as opposed to support recovery, if $\lambda$ is fixed with $\lambda>1,$ a simple thresholding test based on the largest eigenvalue of the matrix $A$ provides detection error probability converging to zero , while if $\lambda<1$ no test based solely on the eigenvalues of $A$ can achieve vanishing probability of error . It remains, however, to establish a solid connection between the detection and estimation problem for submatrix localization for spectral methods.

4 Computational barriers

A recent line of work has uncovered a fascinating interplay between statistical optimality and computational efficiency for the recovery problem and the related detection and estimation problem.The papers considered the biclustering version of the submatrix localization problem (1). Assuming the hardness of the planted clique problem, rigorous computational lower bounds have been obtained in through reduction arguments. In particular, it is shown in that when $K=n^{\alpha}$ for $0<\alpha<2/3$, merely achieving the information-theoretic limits of detection within any constant factor (let alone sharp constants) is as hard as detecting the planted clique; the same hardness also carries over to exact recovery in the same regime. Furthermore, it is shown that the hardness of estimating this type of matrix, which is both low-rank and sparse, highly depends on the loss function [19, Section 5.2]. For example, for $K=\Theta(\sqrt{n})$, entry-wise thresholding attains an $O(\log n)$ factor of the minimax mean-square error; however, if the error is gauged in squared operator norm instead of Frobenius norm, attaining an $O(\sqrt{n}/\log n)$ factor of the minimax risk is as hard as solving planted clique. Similar reductions have been shown in for exact recovering of the submatrix of size $K=n^{\alpha}$ and the planted clique recovery problem for any $0<\alpha<1$.

The results in revealed that the difficulty of submatrix localization crucially depends on the size and planted clique hardness kicks in if $K=n^{1-\Theta(1)}$. In search of the exact phase transition point where statistical and computational limits depart, we further zoom into the regime of $K=n^{1-o(1)}$. We showed in no computational gap exists in the regime $K=\omega(n/\log n)$, since a semi-definite programming relaxation of the maximum likelihood estimator can achieve the information limit for exact recovery with sharp constants. The current paper further pushes the boundary to $K\geq\frac{n}{\log n}(\frac{1}{8e}+o(1))$, in which case the sharp information limits can be attained in nearly linear-time via message passing plus clean-up. However, as soon as $\limsup_{n\to\infty}K\log n/n<\frac{1}{8e}$, there is a gap between the information limits and the sufficient condition of message passing plus clean-up, given by $\lambda>1/e$. For weak recovery, a similar departure emerges whenever $K=o(n)$.

Justification of state evolution equations

In this section, we justify the state evolution equations by establishing the following key lemma. The method of moments is used, closely following . Remark 6 describes the main differences between the analysis here and in .

where $\mu_{t}$ and $\tau_{t}$ are defined in (6) and (7), respectively.

Let ${\mathcal{T}}^{t}$ denote the family of labeled trees $T$ with exactly $t$ generations satisfying the conditions:

For a vertex $u$ that is not the root or a leaf, the mark $r(u)$ is set to the number of children of $v$.

Note that $t=\max_{v\in L(T)}|v|$. All leaves $u$ with $|u|\leq t-1$ have mark .

Let ${\mathcal{T}}^{t}_{i\to j}\subset{\mathcal{T}}^{t}$ be the subfamily satisfying the following additional conditions:

The root has a single child with type distinct from $i$ and $j$.

Similarly, let ${\mathcal{T}}^{t}_{i}\subset{\mathcal{T}}^{t}$ be the subfamily satisfying the following:

The root has a single child with type distinct from $i$.

We point out that under the above definition, a vertex of a tree in ${\mathcal{T}}^{t}$ can have siblings of the same type and mark. Also two trees in ${\mathcal{T}}^{t}$ are considered to be the same if and only if the labels of all nodes are the same, with the understanding that the order of the children of any given node matters. In addition, the mark of a leaf $u$ with $|u|=t$ is not specified and can possibly take any value in $\{0,\ldots,d\}$. The following lemma is proved by induction on $t$ and the proof can be found in [11, Lemma 3.3].

Since the initial messages are zero, $f(\theta_{i\to j}^{0},0)=q_{0}^{0}.$ Thus, for notational convenience in what follows, we can assume without loss of generality that $f(x,0)\equiv q_{0}^{0}$, i.e., $f(x,0)$ is a degree zero polynomial. With this assumption, it follows that for a labeled tree $T\in{\cal T}^{t},$ $\Gamma(T,\mathbf{q},t)=0$ unless the mark of every leaf of $T$ is zero. If the mark of every leaf is zero, then $\theta(T)=1,$ because in this case $\theta(T)$ is a product of terms of the form $0^{0},$ which are all one, by convention. Therefore, $\Gamma(T,\mathbf{q},t)\theta(T)=\Gamma(T,\mathbf{q},t)$ for all $T\in{\cal T}_{t}.$ Consequently, the factor $\theta(T)$ can be dropped from the representations of $\theta_{i\to j}^{t},$ $\theta_{i}^{t},$ $\xi_{i\to j}^{t},$ and $\xi_{i}^{t}$ given in Lemma 2. Applying Lemma 2, we can prove that all finite moments of $\theta_{i}^{t}$ and $\xi_{i}^{t}$ are asymptotically the same.

For any $t\geq 1$, there exists a constant $c$ independent of $n$ and dependent on $m,t,d,C$ such that for any $i\in[n]$:

As explained just before the lemma, the assumption that $f(x,0)\equiv q_{0}^{0}$ implies that the factor $\theta(T)$ can be dropped in the representations given in Lemma 2. Therefore, it follows from Lemma 2 that for $t\geq 1,$

$Q$ consists of $m$-tuples of trees $(T_{1},\ldots,T_{m})$ such that there exists an edge $(r,s)$ in $E(G_{2})\cup E_{J}$ which is covered exactly once.

$R_{1}$ consists of $m$-tuples of trees $(T_{1},\ldots,T_{m})$ such that all edges in $E(G_{2})\cup E_{J}$ are covered at least twice and at least one of them is covered at least $3$ times.

$R_{2}$ consists of $m$-tuples of trees $(T_{1},\ldots,T_{m})$ such that each edge in $E(G_{2})\cup E_{J}$ is covered exactly twice and the graph $G$ contains a cycle.

$R_{3}$ consists of $m$-tuples of trees $(T_{1},\ldots,T_{m})$ such that each edge in $E(G_{2})\cup E_{J}$ is covered exactly twice and the graph $G$ is a tree.

First consider $R_{1}$. Further, divide $R_{1}$ according to the total number of edges in $T_{1},\ldots,T_{m}$ and the number of edges in $E(G_{1})$ which are covered exactly once. In particular, for $\alpha=1,\ldots,m(d+1)^{t}$ and $k=0,1,\ldots,\alpha$, let $R_{1,\alpha,k}$ denote the subset of $R_{1}$ consisting of $m$-tuples of trees $T_{1},\ldots,T_{m}$ such that there are $\alpha$ edges in $T_{1},\ldots,T_{m}$ and there are $k$ edges in $E(G_{1})$ which are covered exactly once. It suffices to show that

Fix $\alpha,k$ and an $m$-tuple of trees $(T_{1},\ldots,T_{m})\in R_{1,\alpha,k}$. Then

We consider breaking $R_{1,\alpha,k}$ down into a large number of smaller sets. While large, the number of these smaller sets depends on $m,t,d$, but not on $n.$ One way to describe these sets is that they are equivalence classes for the following equivalence relation over $R_{1,\alpha,k}:$ Two $m$-tuples in $R_{1,\alpha,k}$ are equivalent if there is a permutation of the set of types $[n]$ such that $i$ maps to $i$, $C^{*}$ maps $C^{*}$, and the second $m$-tuple is obtained by applying the permutation to the types of the vertices of the first $m$-tuple. In particular the marks of the two $m$-tuples must be the same.

Another way to think about these equivalence classes is the following. Given an $m$-tuple $(T_{1},\ldots,T_{m})$ in $R_{1,\alpha,k}$, form the graph $G$ as described above. Let the type of each vertex in $G$ be the common type of the vertices it represents in the $m$-tuple. For convenience, refer to the vertex of $G$ with type $i$ as vertex $i.$ Let $V_{1}$ be the set of vertices in $G$ with types in $C^{*}\backslash\{i\}$ and $V_{2}$ be the set of vertices in $G$ with types in $([n]-C^{*})\backslash\{i\}.$ Record $V_{1}$ and $V_{2},$ and then erase the types of the vertices in $G\backslash\{i\}.$ Then the class of $m$-tuples equivalent to $(T_{1},\ldots,T_{m})$ is the set of $m$-tuples in $R_{1,\alpha,k}$ that can be obtained by assigning distinct types to the vertices of $G$ (which are inherited by the corresponding vertices in the $m$-tuple of trees) consistent with the specified vertex of type $i$ and sets $V_{1}$ and $V_{2}.$ Note that the marks (as opposed to the types) of all $m$-tuples in the equivalence class are the same as the marks on the representative $m$-tuple.

The number of equivalence classes is bounded by a function of $m,t,d$ alone, because the total number of vertices of an $m$-tuple $(T_{1},\ldots,T_{m})$ is bounded independently of $n$, therefore so are the number of ways to partition these vertices to be identified with each other to form vertices in a graph $G,$ along with binary designations on the subsets of the partitions of whether the types of the vertices in the subset are in $C^{*}$ or not (i.e. determining $V_{1}$ and $V_{2}$) and the number of ways to assign marks to the vertices of the trees. Not all partitions with binary designations on the partition subsets correspond to valid equivalence classes because valid partitions must respect the non-backtracking rule and they should have all the root vertices in the same partition set. Also, whether the type of the subset of the partition containing the root vertices corresponds to a type in $C^{*}$ or not is already determined by $i.$ The purpose here is only to verify that the number of such equivalence classes is bounded above by a function of $m,t,d,$ independently of $n.$

Hence, fix such an equivalence class $S\subset R_{1,\alpha,k}.$ It follows from (22)

Note that $|S|\leq K^{n_{1}}n^{n_{2}},$ where $n_{i}=|V_{i}|$ for $i=1,2,$ because there are at most $K$ choices of type for each vertex in $V_{1}$ and fewer than $n$ choices of type for each vertex in $V_{2}.$ The graph $G$ is connected (because all the trees have a root of type $i$), so $n_{1}+n_{2}$ (the number of vertices of $G$ minus one) is less than or equal to the number of edges in $G.$ The number of edges in $G$ is at most $k+\frac{\alpha-k-1}{2}$ because there are $k$ edges in $G$ covered once, and the rest are covered at least twice, with one edge covered at least three times. So $n_{1}+n_{2}\leq k+\frac{\alpha-k-1}{2}.$ Also, since $k$ of the edges in $G$ have both endpoints in $C^{*}$, and the vertices of $V_{2}$ have types in $[n]-C^{*},$ there are at most $\frac{\alpha-k-1}{2}$ edges in $G$ with at least one endpoint in $V_{2}.$ Therefore, since $G$ is connected, $n_{2}\leq\frac{\alpha-k-1}{2}$; otherwise, there must exist a node in $V_{2}$ which has no neighbors in $G$, contradicting the connectedness of $G$. The bound $K^{n_{1}}n^{n_{2}}$ is maximized subject to $n_{1}+n_{2}\leq k+\frac{\alpha-k-1}{2}$ and $n_{2}\leq\frac{\alpha-k-1}{2}$ by letting equality hold in both constraints, yielding $|S|\leq(K)^{k}n^{\frac{\alpha-k-1}{2}}.$ Combining with (23) shows that

where we’ve used the fact that $\frac{\mu K}{\sqrt{n}}$ is bounded independently of $n.$ In a similar way, it can be shown that

Since the number of equivalence classes $S$ does not depend on $n,$ (21) follows.

Next consider $R_{2}$. The previous argument carries over with a minor adjustment. In particular, define $R_{2,\alpha,k}$ accordingly as $R_{1,\alpha,k}$ and then consider an equivalence class $S\subset R_{2,\alpha,k}$ corresponding to some representative $m$-tuple in $R_{2,\alpha,k}.$ Let $G$ and the partition of its vertices into $\{i\}$, $V_{1},$ and $V_{2}$ be determined by the $m$-tuple as before. The number of edges in $G$ is at most $k+\frac{\alpha-k}{2}$ because there are $k$ edges in $G$ covered once, and the rest are covered at least twice. Since $G$ has $n_{1}+n_{2}+1$ vertices, is connected, and has a cycle, $n_{1}+n_{2}$ is less than or equal to the number of edges of $G$ minus one, so $n_{1}+n_{2}\leq k+\frac{\alpha-k-2}{2}.$ Also, since $k$ of the edges in $G$ have both endpoints with types in $C^{*}$, and $V_{2}$ has types in $[n]-C^{*},$ there are at most $\frac{\alpha-k}{2}$ edges in $G$ with at least one endpoint in $V_{2}.$ Therefore, since $G$ is connected, $n_{2}\leq\frac{\alpha-k}{2}$. The bound $K^{n_{1}}n^{n_{2}}$ is maximized subject to these constraints by letting equality hold in both constraints, yielding $|S|\leq K^{k-1}n^{\frac{\alpha-k}{2}}.$ So $|S|\mu^{k}n^{-\alpha/2}\leq\left(\frac{\mu K}{\sqrt{n}}\right)^{k}/K\leq c/K\leq cn^{-1/2},$ and the reminder of the proof for bounding the contribution of $R_{2}$ is the same as for $R_{1}$ above.

The second step is to compute the moments of $\xi^{t}$ in the asymptotic limit $n\to\infty$. We need the following lemma to ensure that all moments of $\xi^{t}$ are bounded by a constant independent of $n$.

For any $t\geq 1$, there exists a constant $c$ independent of $n$ and dependent on $m,t,d,C$ such that for any $i,j\in[n]$

We prove the claim for $\xi_{i}^{t}$; the claim for $\xi_{i\to j}^{t}$ follows by the similar argument. Since $\xi^{0}_{i}=\theta^{0}_{i}=0$ for all $i\in[n]$, it follows from Lemma 2 that

Hence, we only need to check $R_{3}$. Again divide $R_{3}$ according to the total number of edges in $T_{1},\ldots,T_{m}$ and the number of edges in $E(G_{1})$ which are covered exactly once. In particular, $R_{3}=\cup_{1\leq\alpha\leq m(d+1)^{t},0\leq k\leq\alpha}R_{3,\alpha,k}$, where $R_{3,\alpha,k}$ is defined in the similar way as $R_{1,\alpha,k}$. Furthermore, consider dividing $R_{3,\alpha,k}$ into a number of equivalence classes, the number of which depends only on $m,t,d,$ as in the proof of Lemma 3. To prove the lemma, it suffices to show that for any such equivalence class $S,$

In the proof of Lemma 3, we have shown that

We can bound $|S|$ in the similar way as we did for $|R_{1,\alpha,k}|$, with the only adjustment being we cannot use the assumption that there exists at least one edge which is covered at least three times. Fix a representative $m$-tuple $(T_{1},\ldots,T_{m})$ for $S$ and let $G$ and the partition of the vertices of $G$: $\{i\},$ $V_{1}$, $V_{2}$ , be as in the proof of Lemma 3. Let $n_{i}=|V_{i}|$ as before. There are $n_{1}+n_{2}+1$ vertices in the connected graph $G$ and, since the $m$-tuple is in $R_{3,\alpha,k}$, there are at most $k+\frac{\alpha-k}{2}$ edges in $G$, so $n_{1}+n_{2}\leq k+\frac{\alpha-k}{2}.$ Also, at most $\frac{\alpha-k}{2}$ edges of $G$ have at least one endpoint in $V_{2}$ so $n_{2}\leq\frac{\alpha-k}{2}.$ Therefore, $|S|\leq K^{n_{1}}n^{n_{2}}\leq K^{k}n^{\frac{\alpha-k}{2}}.$ It follows that

We also need the following lemma to show the convergence of $\frac{1}{|C^{\ast}|}\sum_{i\in C^{\ast}}(\xi_{i}^{t})^{m}$ in probability using the Chebyshev inequality.

For any $t\geq 1$, $m\geq 1$ and $i,j\in[n]$,

where the same also holds when replacing $\xi^{t}$ by $\theta^{t}$.

We prove the first claim; the other claim follows by a similar argument. Notice that

$Q$ consists of $2m$-tuples of trees such that $G$ is connected and there exists an edge $(r,s)$ in $E(G_{2})\cup E_{J}$ which is covered exactly once.

$R$ consists of $2m$-tuples of trees such that $G$ is connected and all edges in $E(G_{2})\cup E_{J}$ are covered at least twice.

so that for any of the equivalence classes $S\subset R_{\alpha,k}:$

Given a representative $2m$-tuple $(T_{1},\ldots,T_{m},T^{\prime}_{1},\ldots,T^{\prime}_{m})\in R_{\alpha,k},$ the corresponding equivalence class is defined as in Lemma 3. However, in this case there are two distinguished vertices, $i$ and $j$, in the graph $G$, corresponding to the type of the root vertices of the first $m$ trees and the second $m$ trees, respectively. We then let $V_{1}$ be the set of vertices in $G\backslash\{i,j\}$ with types in $C^{*}$ and $V_{2}$ be the set of vertices in $G\backslash\{i,j\}$ with types in $[n]-C^{*}.$ As before, let $n_{1}=|V_{1}|$ and $n_{2}=|V_{2}|.$ There are $n_{1}+n_{2}+2$ vertices in the connected graph $G$ and at most $k+\frac{\alpha-k}{2}$ edges, so $n_{1}+n_{2}\leq k-1+\frac{\alpha-k}{2}.$ At most $\frac{\alpha-k}{2}$ edges have at least one endpoint in $V_{2}$ and $G$ is connected, so $n_{2}\leq\frac{\alpha-k}{2}.$ Thus, $|S|\leq K^{n_{1}}n^{n_{2}}\leq K^{k-1}n^{\frac{\alpha-k}{2}}.$ Hence,

In conclusion, $\mathsf{var}\left(\frac{1}{K}\sum_{i\in C^{\ast}}(\xi_{i}^{t})^{m}\right)\leq c/K$ and the first claim follows. ∎

With Lemma 4 and Lemma 5 in hand, we are ready to compute the moments of $\xi^{t}$ in the asymptotic limit $n\to\infty$.

We prove the first two claims; the last two follows by the similar argument. We prove by induction over $t$. Suppose the following identities hold for $t$ and all $m\geq 1$:

where $Z_{t}\sim{\mathcal{N}}(0,1)$. We aim to show they also hold for $t+1$. Notice that the above identities hold for $t=0$, because $\xi^{0}_{i\to j}=0$ for all $i\neq j$ and $\mu_{0}=\tau_{0}=0$. Let ${\mathcal{F}}_{t}$ denote the $\sigma$-algebra generated by $A^{0},\ldots,A^{t-1}$.

In view of Lemma 5 and Chebyshev’s inequality,

We now fix $i\notin C^{\ast}$. Following the previous argument, one can easily check that

Using the central limit theorem and Chebyshev’s inequality, one can further show that

Fix a subsequence $n_{k}$. In view of Lemmas 3 and 6, for any fixed integer $m$,

Combining Lemma 5 with Chebyshev’s inequality,

Since Gaussian distribution are determined by its moments, by the method of moments (see, for example, [7, Theorem 4.5.5]), applied for each outcome $\omega$ in the underlying probability space (excluding some subset of probability zero), it follows that the sequence of empirical distribution of $\theta_{i}^{t}$ for $i\in C^{\ast}$ weakly converges to ${\mathcal{N}}(\mu_{t},\tau_{t}^{2})$, which, since Gaussian density is bounded, is equivalent to convergence in the Kolmogorov distance,This follows from the fact that when one of the distributions has bounded density the Lévy distance, which metrizes weak convergence, is equivalent to the Kolmogorov distance (see, e.g. [23, 1.8.32]). proving the desired (32). ∎

Proofs of algorithm correctness

Theorems 1-2 are proved in this section. Lemma 1 implies that if $i\in C^{\ast}$, then $\theta_{i}^{t}\sim{\mathcal{N}}(\mu_{t},\tau_{t}^{2})$; if $i\notin C^{\ast}$, then $\theta_{i}^{t}\sim{\mathcal{N}}(0,\tau_{t}^{2})$. Ideally, one would pick the optimal $f(x,t)=e^{\mu_{t}(x-\mu_{t})}$ which result in the optimal state evolution $\mu_{t+1}=\sqrt{\lambda}e^{\mu_{t}^{2}/2}$ and $\tau_{t}=1$ for all $t\geq 1$. Furthermore, if $\lambda>1/e$, then $\mu_{t}\to\infty$ as $t\to\infty$, and thus we can hope to estimate $C^{\ast}$ by selecting the indices $i$ such that $\theta_{i}^{t}$ exceeds a certain threshold. The caveat is that Lemma 1 needs $f$ to be a polynomial of finite degree. Next we proceed to find the best degree-$d$ polynomial for iteration $t$, denoted by $f_{d}(\cdot,t)$, which maximizes the signal to noise ratio.

Recall that the Hermite polynomials $\{H_{k}:k\geq 0\}$ are the orthogonal polynomials with respect to the standard normal distribution (cf. [25, Section 5.5]), given by

Truncating and normalizing the series at the first $d+1$ terms immediately yields the solution to (12) as the best degree-$d$ $L_{2}$-approximation to the relative density, described as follows:

where $G_{d}(\mu)=\sum_{k=0}^{d}\frac{\mu^{k}}{k!}$. Define

where $a_{k}\triangleq\frac{\widehat{\mu}_{t}^{k}}{k!}(\sum_{k=0}^{d}\frac{\widehat{\mu}_{t}^{2k}}{k!})^{-1/2}$. Then $f_{d}(\cdot,t)$ is the unique maximizer of (12) and the state evolution (6) and (7) with $f=f_{d}$ coincides with $\tau_{t}=1$ and $\mu_{t}=\widehat{\mu}_{t}$. Furthermore, for any $d\geq 2$ the equation

has a unique positive solution, denoted by $a^{*}_{d}$. Let $\lambda^{*}_{d}=\frac{1}{G_{d-1}(a^{*}_{d})}$ and define $\lambda^{*}_{1}=1$. Then

$\lambda^{*}_{d}\downarrow 1/e$ monotonically as $d\to\infty$ according to $\lambda^{*}_{d}=1/e-\frac{1/e^{2}+o(1)}{(d+1)!}$.

The best affine update gives $\lambda^{*}_{1}=1$; for the best quadratic update, $a_{2}^{*}=\sqrt{2}$ and hence $\lambda^{*}_{2}=\frac{1}{1+\sqrt{2}}\approx 0.414$. More values of the threshold are given below, which converges to $1/e\approx 0.368$ rapidly.

which is finite for any $\lambda>1/e$. Then for any $d\geq d^{*}$, $\widehat{\mu}_{t}\to\infty$ as $t\to\infty$. We note that as $\lambda$ approaches the critical value $1/e$, the degree $d^{*}(\lambda)$ blows up according to $d^{*}(\lambda)=\Theta(\log\frac{1}{\lambda e-1}/\log\log\frac{1}{\lambda e-1})$, as a consequence of the last part of Lemma 7.

For $d=1$, the best state evolution is given by

and the corresponding optimal update rule is

This is strictly better than $f(x,t)=x$ described in Section 1.3 which gives $\widehat{\mu}_{t+1}^{2}=\lambda\widehat{\mu}_{t}^{2}$; nevertheless, in order to have $\widehat{\mu}_{t}\to\infty$ we still need to assume the spectral limit $\lambda\geq 1$.

Next we analyze the behavior of the iteration (36). The case of $d=1$ follows from the obvious fact that $\widehat{\mu}_{t+1}^{2}=\lambda(\widehat{\mu}_{t}^{2}+1)$ diverges if and only if $\lambda\geq 1$. For $d\geq 2$, note that $G_{d}$ is a strictly convex function with $G_{d}(0)=1$ and $G_{d}^{\prime}=G_{d-1}$. Also, $(G_{d}(a)-aG_{d-1}(a))^{\prime}=-aG_{d}^{\prime\prime}(a)<0$. Thus, $G_{d}(a)-aG_{d-1}(a)$ is strictly decreasing on $a>0$ with value $\frac{1}{d!}$ at $a=1$ and limit $-\infty$ as $a\to\infty,$ so (38) has a unique positive solution $a^{*}_{d}$ and it satisfies $a^{*}_{d}>1.$ Furthermore, (G_{d}(a)-aG_{d-1}(a))^{\prime}\big{|}_{a=1}=\sum_{k=0}^{d-2}\frac{1}{k!}, so by Taylor’s theorem,

Consider next the values of $\lambda$ such that $\widehat{\mu}_{t}$ diverges. For very large $\lambda$, $G_{d}(a)$ dominates $a/\lambda$ pointwise and $\widehat{\mu}_{t}$ diverges. The critical value of $\lambda$ is when $G_{d}(a)$ and $a/\lambda$ meet tangentially, namely,

whose solution is given by $a=a_{d}^{*}$ and $\lambda=\lambda^{*}_{d},$ where

Thus, $\lambda_{d}^{*}$ is the minimum value such that for all $\lambda>\lambda^{*}_{d}$, $\lambda G_{d}(a)>a$ for all $a>0,$ so that starting from any $\widehat{\mu}_{t}\geq 0$ we have $\widehat{\mu}_{t}\to\infty$ monotonically. The fact $\lambda_{d}^{*}$ is decreasing in $d$ follows from the fact $G_{d}$ is pointwise increasing in $d.$ ∎