An Unbiased Symmetric Matrix Estimator for Topology Inference under Partial Observability

Network topology inference is a fundamental problem in many applications of network science, such as locating the source of fake news, brain connectivity networks detection, etc. Many real-world situations suffer from a critical problem that only a limited part of observed nodes are available. This letter considers the problem of network topology inference under the framework of partial observability. Based on the vector autoregressive model, we propose a novel unbiased estimator for the symmetric network topology with the Gaussian noise and the Laplacian combination rule. Theoretically, we prove that it converges to the network combination matrix in probability. Furthermore, by utilizing the Gaussian mixture model algorithm, an effective algorithm called network inference Gauss algorithm is developed to infer the network structure. Finally, compared with the state-of-the-art methods, numerical experiments demonstrate the proposed algorithm enjoys better performance in the case of small sample sizes.


Introduction
In network topology inference, the goal is to identify the network structure from the collected signals, which plays a vital role in many applications [1][2][3]. In most existing research [4][5][6], the network topology was inferred based on a complete set of observed data for all the entities of interest.
However, many real-world applications suffer from the partial observability problem. For example, when investigating a social network with millions of members, only a limited number of observed nodes are available [7,8]. Thus, it is necessary to study the network topology inference based on partial observations, (see Fig. 1). Since the inference performance is affected by the missing observations [9][10][11][12], network topology inference is challenging when the signals of the hidden nodes are not observed.
Several works formulated partially observed network topology inference as sparse or convex optimization problems [13,14]. Recently, the problem of general network topology inference using Erdős-Rényi (ER) models [15] under partial observability was considered in [16]. Moreover, some attractive estimators, such as the Granger, one-lag, and residual estimators, were proposed in [16,17]. Furthermore, it was proven that these estimators can be used to recover the subgraph of observed network nodes when the network size grows.
Inspired by the estimators proposed in [17], we propose an efficient unbiased symmetric matrix estimator to infer the network topology under partial observability. Compared with the work in [17], we focus on the case in which the Notations: The total number of network nodes is denoted by N , and the number of samples on each node is n. For an N × N matrix Z, the submatrix spanning the rows and columns of Z indexed by set S ⊂ {1, . . . , N } is denoted by Z S . The fraction of the observed nodes is denoted as ξ = |S| N , where |S| represents the cardinality of the set S. Finally, I denotes the identity matrix.

Problem Formulation
We consider an undirected random network, such as ER and Barabási-Albert (BA) random graphs, with N nodes.
At time t, each node k collects the output signal y t (k) according to the following network diffusion process (i.e., running a VAR model) [7] where x t (k) is the input signal, µ ∈ (0, 1) is the corresponding weighting coefficient and is assumed to be known.
Here, x t (k) for t = 1, . . . , n are assumed to be independent and identically distributed (i.i.d.) standard normal random variables, i.e., x t (k) ∼ N (0, 1). The N × N matrix with entries A kj is denoted as the network combination matrix A that satisfies the following assumption. Assumption 1. The network combination matrix A is symmetric and obtained by applying the Laplacian combination rule [16]: where λ is a parameter with λ ≤ 1, G kj is the entry of matrix G, which is the adjacency matrix of the undirected random network, the degree of the k-th node is defined as d k = 1 + z =k G kz , and d max is the maximum degree of matrix G.
Since the network combination matrix A depends on the adjacency matrix G, it also represents the network topology structure. Specifically, A kj equals 0 if the k-th node disconnects from the j-th node, and A kj is a positive number (1 − µ) λ dmax , if the k-th node connects with the j-th node.
By stacking the input and output signals across the network at time t, equation (1) is rewritten as Assume that only a subset of the network nodes can be monitored. The monitored set is denoted as S ⊂ {1, . . . , N }. We only receive a partial output signal [y t ] S {y t (k)|k ∈ S}. However, according to (2), the observed signals are influenced by the hidden nodes.
In network topology inference, the goal is to estimate the partial combination matrix A S and recover the partial adjacency matrix G S by the data monitored [y t ] S .

Main Results
This section derives an unbiased estimator for A S and a new node pair clustering based on the GMM method. Finally, we present the NIG algorithm to infer the network topology.

An Unbiased Estimator of Matrix A S
First, by summing up x t in (2) from 0 to t with the assumption y 0 = µx 0 , x 0 (k) ∼ N (0, 1) for any k, we obtain Multiplying both sides of equation (2) by y T t and taking the expectations, we obtain Since where where the last equality comes from the i.i.d. property and unit variance of the random variables x t . Using (5), we obtain Remark 1. Inspired by the residual estimator in [17] and using (6)-(7), we can nullify the error arising from higherorder powers to obtain the combination matrix A by subtraction between R 1 (t) and A 3 R 0 (t − 1), which helps us to design an unbiased estimator for the partial combination matrix A S .
Since y t is an N × 1 vector for any t, we propose an N × N estimator of the matrix A at time t: and denote [A u (t)] S as the submatrix of A u (t), which is where the last equality uses the property of the block matrix product.
Equation (9) reveals that the estimator [A u (t)] S can be obtained by the observed signals. It is worth noting that the performance of the proposed estimator also depends on the hidden nodes, see Remark 2. Next, we show that the proposed estimator benefits from the following property.
Using the definition of R 1 (t) in (4), we obtain Additionally, using the recursion of (2), we obtain where the first equality uses the i.i.d. property. Using (6), we obtain Substituting (7) and (11) into (10), This proves that the estimator A u (t) is an unbiased estimator of A. Considering the partial observation with (9), we have Thus, [A u (t)] S is an unbiased matrix estimator of A S .

Remark 2. The performance of the proposed estimator [A u (t)]
S not only depends on the observed signals but also on the hidden nodes since the hidden nodes influence the observed signals. Moreover, we also showed the hidden nodes do impact the variance of the proposed unbiased estimator in the supplementary material section.
Proposition 1 shows that we obtain an unbiased estimator for A S at each time t. However, this estimator may have a large variance. To overcome this problem, we propose another unbiased matrix estimator that converges in probability to A S . Based on the unbiased estimator and the cumulative samples, we denote Then, an efficient matrix estimator using cumulative samples is proposed as n ] S is also unbiased. Besides, it has the following property. Theorem 1. Under Assumption 1, if the observed signals are collected by the VAR model with Gaussian noise in (1), Proof. For brevity, the detailed proof is available in the supplementary material.
Theorem 1 demonstrates that the unbiased estimator [ A u n ] S converges to A S in probability. However, for a fixed network size N , the estimators proposed in [17] cannot guarantee convergence to A S as the number of samples grows. In comparison, the advantage of the Granger estimator is that it does not rely on the symmetry assumption. Therefore, it can be used in non-symmetric scenarios.
The detailed comparisons are shown in Table 1. [ ( * ) means that for a fixed network size N, the estimator converges in probability to A S when n → ∞.

Node Pair Clustering
We have developed an unbiased estimator [ A u n ] S , as denoted in (13), that quantifies the strength of the connections among the network nodes. Considering the binary property of the adjacency matrix, we require the standard clustering algorithms, which allow grouping the entries of the unbiased estimator into two clusters, indicating the 'disconnection' and 'connection' pairs of the observed network.
To select a suitable clustering method, we complete experiments on ER models to visualize the distribution of the entries of the unbiased matrix estimator [ A u n ] S in Fig. 2. For convenience, the 'disconnection' and 'connection' clusters are denoted as '0' and '1', respectively. 2 Fig. 2 (a) shows that entries are largely mixed under the small size situation. In addition, their distributions also partly overlap. Thus, it is difficult to infer the network topology when the sample size is small. However, when the sample size n is large, as shown in Fig. 2 (b), the entries are perfectly separated into two clusters, which is consistent with Theorem 1. The histograms in Fig. 2 reveal that their distribution seems to be a Gaussian mixture distribution.
Thus, we adopt the GMM algorithm [19] to group the entries into two clusters. It is known that the GMM algorithm is implemented by the expectation-maximum algorithm, which is sensitive to the initial conditions [20]. Since the value for the disconnected pairs under the Laplacian rule is 0, the expectation of the unbiased estimator for the disconnected pairs is also 0 (see Fig. 2). Thus, we set 0 as the initial mean for the disconnected pairs when using the GMM algorithm.

Proposed Method
The proposed NIG method is summarized in Algorithm 1, which is used to infer network topology under partial observability. Steps 1-3 are completed to generate the unbiased estimator [ A u n ] S . Step 4 is completed using the GMM clustering method for the entries of the unbiased matrix estimator.

Numerical Experiments
In this section, we examine the numerical performance of the proposed NIG algorithm and present comparison results with three existing methods, namely, one-lag (One), residual (Res), and Granger (Gra) estimators followed by   the K-means clustering algorithm in [17]. If not otherwise specified, the NIG algorithm uses the GMM clustering algorithm. The error rates are adopted for performance evaluation, which are the fractions of network node pairs that are incorrectly recovered.
Consider a network with N nodes, where, according to the ER model, edges are determined by independent Bernoulli distributions with link probability p. If not otherwise specified, we set N = 400, ξ = 0.2, µ = 0.1, p = 0.1. The parameter λ of the Laplacian rule is set as 0.99. To assess the performance of the proposed unbiased estimator in the NIG algorithm, we consider the following scenarios, where the results are averaged over 50 experiments.
• The proposed algorithm is compared with other methods in terms of error rate, bias, and variance in Fig. 3.
• An ablation study is provided to clarify where the gain of the proposed algorithm arises. Specifically, in Fig.  3 (d), we show the results from the K-means and GMM algorithms for the proposed unbiased estimator and Granger estimator, respectively.
• We examine the error rates and 95% confidence intervals of the NIG algorithm when varying N , p, and ξ in Fig. 4. Moreover, in Fig. 4 (d), we compare these algorithms in terms of a popular metric, the false-negative score [21], which is defined as the fraction of disconnected node pairs declared as connected. Fig. 3 (a) shows that the proposed estimator has the lowest error rate compared to the other estimators in the case of small sample sizes. However, when the sample size becomes larger than 0.8 * 10 5 , the proposed algorithm has higher error rates than the Granger and one-lag estimators. The reason may be that the proposed estimator has the smallest bias, which mainly affects the performance and makes the proposed unbiased estimator perform best in the small sample size situation. However, with the gradual increase in sample sizes, the variance becomes the main factor affecting the performance. Thus, the larger variance leads to the higher error rates of the proposed NIG algorithm than baselines (see Fig. 3 (b)-(c)). In addition, one can observe that its variance converges to 0 as n increases, which is consistent with Theorem 1. Fig. 3 (d) shows that both the proposed and Granger estimators with the GMM algorithm outperform those with the Kmeans algorithm. Thus, the GMM algorithm provides some gains for both estimators. When the sample size is small, one can also observe that the proposed estimator has a better performance than the Granger estimator with the same GMM algorithm. The reason is that the unbiased estimator provides a better initialization for the GMM clustering algorithm. Thus, the gains of the proposed method arise from both the new estimator and the GMM algorithm.
Comparing the proposed method with the baselines in Fig. 4 (a) and Fig. 3 (a), we see that the performance gap grows with larger network sizes. Fig. 4 (b) shows that all the estimators have higher error rates when the graph becomes denser, but the proposed estimator still performs best with small sample sizes. Fig.4 (c) shows that, with a higher fraction of observed nodes, the proposed NIG algorithm has a narrower confidence interval and a lower average error rate. The reason may be that the impact of the hidden nodes on the observed signals decreases when the fraction of the observed nodes increases since there are fewer hidden nodes. In addition, Fig. 4 (d) shows the false-negative score. A satisfactory estimator should have small false-negative scores. Fig. 4 (d) shows that the proposed method has the lowest false-negative score and demonstrates that the proposed method has the highest probability of recovering the disconnected pairs correctly. On the other hand, the proposed NIG algorithm has contrary performance under the false-positive score [21].

Conclusion
In this letter, the novel NIG algorithm, including an unbiased estimator and a GMM algorithm for network topology inference under partial observability, is proposed. Theoretically, we have proven that the unbiased estimator converges to the network combination matrix in probability and helps us obtain a better initialization for the GMM algorithm.
The simulations have shown that the proposed NIG algorithm obtains a better performance with small sample sizes, sparse networks, and high fractions of observed nodes. Future work aims to improve the algorithms by combining the advantages of the biased estimators and the proposed method.

Proof of Theorem 1
Let us define the combination matrix estimator A u (S1) Then, applying (S1) on A u n , we have where Tr is the trace of a matrix. From the VAR model in (2), we have y t+j = A j y t + µ( j k=1 A j−k x t+k ), which is multiplied on both sides by y T t , yielding y t+j y T t = A j y t y T t + µ j k=1 A j−k x t+k y T t , ∀j ≥ 1. Considering M 1 + M 2 + M 3 + M 4 , using (S1), we have Based on the definitions of R 1 (n), R 3 (n − 1), (S1), and (S3) we have n ] S is an unbiased estimator of A S , using the multidimensional Chebyshev's inequality [24], we have for any > 0. Then, substituting Therefore, we derive [ A u n ] S P → A S when n → ∞. Cov(x k+j y T k , x t+j y T t )], we first consider the trace of the following covariance, where the first equality comes from the definition of covariance and uses the fact that where the last equation comes from y t ∼ N (0, R 0 (t)) in equation (3). Finally, combining (S8) with (S6), we obtain Tr Cov 1 n n t=1 x t+j y T t = N n 2 Tr n t=1 R 0 (t) .
Using the Stolz Theorem [25], we have lim n→∞ 1 n n t=1 R 0 (t) = lim n→∞ R 0 (n). Recall that R 0 (t) = µ 2 (I + A 2 + · · · + A 2t ), yielding that lim According to the VAR model (2), we have y j y T where the first equality uses (S6) and the last equality uses property where ⊗ means the Kronecker product. Since y t ∼ N (0, R 0 (t)), the streaming data y t y T t follows a Wishart distribution. Then we have E[y t y T t ] = R 0 (t) and Cov(y t y T t ) = 2R 0 (t) ⊗ R 0 (t) in [26,Proposition 8.3]. Substituting this into (S10), we obtain Cov y j y T j , y t y T . Thus, we have Then, we apply the corollary of the Stolz theorem [25, Example 3.75], which says that two convergent sequences a n , b n satisfy that lim n→∞ a n = a, lim n→∞ b n = b, then, lim n→∞ a 1 b n + a 2 b n−1 + · · · + a n b 1 n = ab. (S12) One may notice that a i+1 , b n−i ∈ R. However, (S12) also holds when (a n ) n∈N and (b n ) n∈N are matrix sequences by just using a matrix norm (for example, the Frobenius norm) to replace the absolute value in the proof of [25, Example 3.75]. We provide the proof for the matrix version here: Let two matrix sequences (A n ) n∈N and (B n ) n∈N be given and let us define the sequence (C n ) n∈N by C n = A 1 B n + A 2 B n−1 + · · · + A n B 1 n , n = 1, 2, . . . Since the sequence (B n ) n∈N is bounded, i.e., B n F ≤ b for every n ∈ N, using the sub-additive and submultiplicative properties of matrix norm, i.e., M 1 + M 2 F ≤ M 1 F + M 2 F , M 1 M 2 F ≤ M 1 F M 2 F , we have C n F ≤ A1 F + A2 F +···+ An F n b and lim n→∞ C n F ≤ lim n→∞ b n n k=1 A k F = 0. This means that lim n→∞ C n = 0.
Proof (b): Let us put lim n→∞ A n = a and X n = A n − A for every n ∈ N. Then it holds lim n→∞ X n = 0 and C n = (X 1 + A) B n + · · · + (X n + A) B 1 n = X 1 B n + · · · + X n B 1 n + A · B n + · · · + B 1 n = F n + G n .
Since lim n→∞ B n = B, the sequence (B n ) n∈N is bounded, according to (a), it follows that lim n→∞ F n = lim n→∞ X 1 B n + · · · + X n B 1 n = 0 while lim n→∞ G n = lim n→∞ A · B n + · · · + B 1 n = AB.
Therefore lim n→∞ C n = AB. Then we have completed the proof of the corollary of matrix version.
Back to the proof of 6.3, we can set that B n−t = 1 + 2 n j>t A j−t ⊗ A j−t n = 1 + 2A ⊗ A + · · · + 2A n−t ⊗ A n−t n .
Then, the matrix in trace operation on the right side of (S11) has the shape of (S12). Moreover, using the Stolz Theorem and ρ(A) < 1, we have lim n→∞ B n = lim