ADDITIVE SEPARABLE GRAPHON MODELS

• ``Theoretically, it raises the question: why is an additive separable scheme considered reasonable? Why not a multiplicative or another type of schema?''

Answer: Thank you for your question. The additive separable form is a natural choice, grounded in the insight provided by the Hilbert–Schmidt theorem (see for example [1]), which implies that any bounded graphon function $f(u,v)$ can be decomposed as: $f(u,v) = \sum_{j=1}^\infty \lambda_j G_j(u) G_j(v),$ where $G_j$ are orthonormal eigenfunctions and $\sum_{j} \lambda_j^2 < \infty$ [2].

A practical approach to estimating this infinite series is to truncate it at $r$, retaining the $r$ eigenfunctions corresponding to the $r$ largest absolute eigenvalues. This results in a low-rank structure that aligns naturally with our model and facilitates efficient estimation. An additional advantage of this approach is that the resulting connection probability matrix $P$—formed by evaluating the graphon function at the nodes—naturally inherits the same low-rank property as the truncated graphon. We believe this model provides a robust and insightful framework for investigation. That said, exploring multiplicative or other alternative schemas remains an important direction for future research.

[1] Szegedy, B. (2011). Limits of kernel operators and the spectral regularity lemma. European Journal of Combinatorics, 32(7), 1156–1167.

[2] Bickel, P. J., Chen, A., & Levina, E. (2011). The method of moments and degree distributions for network models. The Annals of Statistics, 39(5), 2280–2301.

Thank you for your interest in our paper and for your thoughtful review. We appreciate your positive feedback on the clarity of our method, as well as our theoretical results and simulation results. We summarize the main changes according to the valuable suggestions from reviewers as follows.

1. We provide a complete theoretical framework for general $r$, updated in Lines 203–338 on pages 4–7.
2. Computationally, we introduce a new approximation approach for computations, updated in Lines 339–381 on pages 7–8. None of the theoretical results are affected, while it achieves a time complexity matching that of matrix multiplication ($O(n^{2.373})$), eliminating the need for subsampling techniques.
3. In Appendix A.3, we provide a method for selecting $r$ when it is unknown. We also point out this in Remark 2 of the main article.
4. We update the simulation results and included an additional simulation for $r = 3$.
5. In Appendix A.1, we include a real data analysis, demonstrating the use of our method for selecting $r$, as well as our algorithms for estimating the connection matrix and the graphon function.
6. We update the mathematical proofs for all theoretical results.

In response to your questions, we have made revisions and improvements to the paper. Below, we address your points one by one.

• ``Empirically, as shown in the comparison tables, the proposed method yields only marginal improvement over existing methods.''

Answer: We would like to clarify that the improvement in performance is influenced by the sparsity of the graphon. For relatively dense graphons (see Table 2), our method performs similarly to USVT, as both achieve minimax optimal rates for estimating the connection probability matrix. However, for sparse graphons (see Table 3), our method significantly outperforms USVT, achieving MSE reductions of over 50%. This underscores the versatility of our approach in effectively addressing both dense and sparse networks.

Importantly, this superior performance is accompanied by significantly lower computational requirements, with a complexity of $O(n^{2.373})$ (see Table 2). A further advantage of our method is that it is tuning-parameter-free, unlike the SAS method (Chan and Airoldi, 2014) and the network histogram method (Olhede and Wolfe, 2014), both of which require careful selection of tuning parameters. Consequently, these methods can be sensitive to the choice of parameters, potentially impacting their performance (see tables below).

Additionally, our method directly leverages the low-rank structure of the graphon to estimate the graphon function itself, not just on discrete grids (as represented by connection probability matrices), but across the entire domain of $[0,1] \times [0,1]$. In contrast, neither the USVT method nor the neighborhood smoothing method can estimate the full graphon function.

Thus, our method offers a computationally efficient (tuning-free with optimal complexity), theoretically optimal (achieving minimax rates), and comprehensive framework for graphon estimation, outperforming its competitors in these key aspects.

The MSE of the SAS (Chan and Airoldi, 2014) under different parameter selections.

The MSE of the Nethist (Olhede and Wolfe, 2014) under different parameter selections.

We have made a second revision to our article, and the main changes are listed below for your reference. In the newly revised version, we have revised the wording to reduce the emphasis on the novelty of the model itself. Instead, we frame our estimation approach as a new attempt to estimate a low-rank representation. Additionally, we have added a new remark (Remark 1) that discusses the power iteration method and included it in our simulation studies (see Tables 2, 3, and 5). Furthermore, we have added Remark 3 (line 289) in the revised manuscript to further highlight the key differences between our method and spectral methods, including differences in motivation, estimation procedure, assumptions, and empirical performance for sparse graphons. It is worth noting that our approach allows for the simultaneous estimation of both the connection probability matrix and the graphon function. For estimating the connection probability matrix, our method shows clear advantages over spectral methods in the sparse region (see Table 3).

Thank you for your interest in our paper and for your thoughtful review. We appreciate your acknowledgment of the logical structure of our writing and the positive feedback on the simulation results. We summarize the main changes according to the valuable suggestions from reviewers as follows.

1. We provide a complete theoretical framework for general $r$, updated in Lines 203–338 on pages 4–7.
2. Computationally, we introduce a new approximation approach for computations, updated in Lines 339–381 on pages 7–8. None of the theoretical results are affected, while it achieves a time complexity matching that of matrix multiplication ($O(n^{2.373})$), eliminating the need for subsampling techniques.
3. In Appendix A.3, we provide a method for selecting $r$ when it is unknown. We also point out this in Remark 2 of the main article.
4. We update the simulation results and included an additional simulation for $r = 3$.
5. In Appendix A.1, we include a real data analysis, demonstrating the use of our method for selecting $r$, as well as our algorithms for estimating the connection matrix and the graphon function.
6. We update the mathematical proofs for all theoretical results.

In response to your questions, we have made revisions and enhancements to the manuscript. Below, we address your questions point by point.

• ``Overall, I would regard the proposed method as a method for estimating the network moment or motifs. That said, I doubt that it would scale easily to even moderate values of $r$, especially due to computational demands. This constraint could make it less practical for most commonly studied stochastic block models, where the connection probability matrix rank is often higher than 2.''

Answer: Thank you for your suggestion, which motivated us to update the theory and implementation for the general case of $r$. In the revised version of our paper, we have included a comprehensive theoretical study for the case $r > 2$, which is now detailed in Lines 203–338 on page 4-7 of the revised version of the paper. In the more general setting, we obtain the same consistent result as before. Computationally, we find a new approximation approach for computing and update it in Lines 339-381 on pages 7-8. Specifically, we approximate simple paths with paths allowing node repetition and prove their asymptotic equivalence, ensuring that this approximation does not affect any theoretical results. This approach enables our method to achieve a time complexity matching that of matrix multiplication ($O(n^{2.373})$), eliminating the need for subsampling techniques. We update all the simulation results using the new methods, and consider a new simulation scenario with rank $r=3$, see ID 7 in Tables 2. Additionally, we consider the problem of selecting $r$ in practice, see discussion on Algorithm 4 in Appendix A.3. and Table 6 for the corresponding simulation results.

• ``Despite what’s suggested in the discussion, I think it would help if the authors could explicitly address the estimation process for general ASG(r) as this is particularly important for practitioners looking to apply it more broadly.''

Answer: Similar to the previous answer, we have now incorporated the theory as well as the methods for choosing $r$ in practice. Thank you for your suggestion.

• ``As a follow-up question, deciding the appropriate rank $r$ could be a challenge in practice. A brief discussion on rank selection—whether there’s a heuristic, or a data-driven way to guide practitioners would make the approach more usable.''

Answer: Thanks a lot! We have proposed a computational feasible methods to sequentially determine $r$, which works well in our numerical analysis, see discussion on Algorithm 4 in Appendix A.3. and Table 6 for the corresponding simulation results.

• ``Given the wide range of applications inspired by network literature, I suggest the authors analyze at least one real-world dataset to demonstrate the practical utility of their method.''

Answer: Thank you for your suggestion. We have added a real data example in section A.1 of the appendix (page 12-13), which comes from contact records in a primary school. Using this data, we selected $r = 4$ and estimated the corresponding connection matrix and graphon function.

Thank you for reviewing our manuscript and providing valuable feedback. It is encouraging that you could acknowledge the clarity, theoretical foundation, and simulation results of our paper. We summarize the main changes according to the valuable suggestions from reviewers as follows.

1. We provide a complete theoretical framework for general $r$, updated in Lines 203–338 on pages 4–7.
2. Computationally, we introduce a new approximation approach for computations, updated in Lines 339–381 on pages 7–8. None of the theoretical results are affected, while it achieves a time complexity matching that of matrix multiplication ($O(n^{2.373})$), eliminating the need for subsampling techniques.
3. In Appendix A.3, we provide a method for selecting $r$ when it is unknown. We also point out this in Remark 2 of the main article.
4. We update the simulation results and included an additional simulation for $r = 3$.
5. In Appendix A.1, we include a real data analysis, demonstrating the use of our method for selecting $r$, as well as our algorithms for estimating the connection matrix and the graphon function.
6. We update the mathematical proofs for all theoretical results.

In response to the weaknesses and questions you pointed out, we have made some additions and clarifications. Below, we address your questions point by point.

• ``Regarding the proposed model Additive Separable Graphons ASG(r), the authors state in line 076 that this model is new. However, it appears to me as a generic low-rank graphon, which is a well-recognized particular case of a general graphon already discussed in the literature (e.g., Chan and Airoldi, 2014; Xu, 2018). Indeed, for any rank-r graphon, the decomposition given in Eq. (1) follows directly from the spectral theorem. I struggle to see this as a new model, and in my opinion, it would be beneficial for the authors to justify why ASG(r) is indeed novel.''

Answer: Thank you for your question. As you noted, our model has appeared in other contexts and represents a natural formulation of low-rank models. However, this work is likely the first to leverage this model specifically for estimating both the connection probability matrix $P$ and, more importantly, the graphon function $f(u,v)$.

While USVT can partially utilize low-rank information to estimate the connection probability matrix, it cannot estimate the graphon function over the full domain $[0,1] \times [0,1]$. It is worth noting that the connection probability matrix is derived by evaluating the graphon function at a finite set of discrete points. In contrast, our model allows direct estimation of the graphon function, leveraging its low-rank structure in a way that distinguishes our approach.

Moreover, our method is tuning-parameter-free, making it computationally faster in practice compared to existing methods that require tuning parameters for graphon estimation. For further discussion, please refer to the third point below.

• ``The power iteration is known to converge in one iteration for matrices of rank 1 and fast for low-rank matrices. Comparing with this approach may be beneficial, as the direct competitors here are spectral methods, and it’s possible the methods are analogous.''

Answer: Thank you for your insightful comments! The connection between our method and spectral methods is indeed an interesting topic. When $r = 1$, the results from the spectral method closely align with those of our method, as shown in Table 5 on page 14.

Regarding power iteration, we would like to point out that although the connection matrix $P$ is rank-1, the corresponding adjacency matrix $E$ may still be full-rank. Additionally, both the spectral method and the power iteration method face challenges when extended to estimating graphon functions over $[0,1]$.

For $r > 1$, the relationship between these methods and ours remains unclear. We have not identified any definitive connections, either through theoretical analysis or simulations. Moreover, our motivations differ fundamentally, as our goal is to estimate the graphon function itself, rather than just the connection probability matrix. It is important to note that the graphon function captures the core of the statistical model, while the connection matrix is merely a realization of the graphon function evaluated at a finite set of discrete points. Furthermore, the connection probability matrix is not a proper limit object in random graph models, whereas the key advantages of graphon models arise precisely in this limiting sense.

From a practical perspective, our method tends to be slightly faster than USVT, as implemented in the R package provided by the authors of USVT.

• ``The notation $O_p$ used throughout the paper is not defined. In the finite sample results, such as in Theorem 1, the results should hold with high probability—this is currently unstated.''

Answer: Thank you for pointing this out. We have now included the definition of $O_p$ in Section 1 in Line 88 on page 2. After properly defining $O_p$, the statement of Theorem 1 inherently includes the meaning of ``hold with high probability''.

• ``In line 434, the authors mention that the absence of a tuning parameter "enhances robustness." Could they provide a clearer justification for this claim?''

Answer: For nonparametric methods (e.g., the SAS method by Chan and Airoldi, 2014, and the network histogram method by Olhede and Wolfe, 2014), a bandwidth parameter $h$ (or $K=n/h$) is required, which affects estimation accuracy and can be challenging to select in practice. In our simulations, we used their suggested default value, which achieves the optimal theoretical convergence rate. Additionally, we conducted extra simulations and found that choosing different values of their parameter results in fluctuations in the MSE, see two tables below. For instance, increasing the $K$ for the SAS method to 1.2 times the default value led to an approximately 1.6-fold increase in MSE. Similarly, increasing the bandwidth $h$ of the network histogram method by a factor of 0.5 leads to a twofold increase in the MSE. This indeed highlights the importance and difficulty of parameter selection in practice. In contrast, our method does not require parameter tuning, eliminating this complication.

The MSE of the SAS (Chan and Airoldi, 2014) under different parameter selections.

The MSE of the Nethist (Olhede and Wolfe, 2014) under different parameter selections.

• ``I am following their proof correctly, it seems that alternative equations involving $\hat\lambda_1$ and $\hat\lambda_2$ could be derived for the rank-2 case. Specifically, it seems that $\sum_{i,j}A_{ij}A_{ji}$ should approximate up to rescaling $\lambda_1^2+\lambda_2^2$. Is this correct, and would it be useful? How does this approach compare with the method they propose?''

Answer: Thank you for your interest in the specifics. We kindly recall that, in this paper, we consider simple undirected graphs for modeling the network (i.e., without multiple edges or self-loops, and $A_{ij}=A_{ji}$). Consequently, we have $\sum_{i,j}A_{ij}A_{ji}=\sum_{i,j}A_{ij}$ which approximate up to rescaling $\lambda_1(\int_0^1 G_1(u)du)^2+\lambda_2(\int_0^1 G_2(u)du)^2$. We propose new methods instead of the previous sampling approach, for effectively estimating ASG(r) model utilizing lines and cycles (see Lines 223–227 on page 5 for a formal definition). We choose these subgraphs because they not only ensure theoretical guarantees (Theorems 3 and 4) but also maintain computational efficiency. As demonstrated in Section 3.3, our method achieves an asymptotically equivalent complexity to matrix multiplication, effectively avoiding the challenges of general subgraph counting. In practice, we also suggest eliminating paths with repetitive nodes subject to the computational constraint to improve the finite sample performance. In fact, when $r$ is small, the extra computational cost of this correction can be small, as confirmed by our simulation (see the running time in Table 2 on page 9).

• ``Difficulty in Extending Results: Although the authors discuss potential extensions, constructing estimators for lower-rank graphons appears challenging as the equations grow complex quickly. Furthermore, the order (in powers of E) of the quantities needed increases with the rank, complicating the computation.''

Answer: Thank you for your suggestion. We have now extended the model to arbitrary $r$, provided $r$ is bounded. The theoretical complexity of our revised model remains $O(n^{2.373})$ (see Remark 3). Our simulation results also confirm that the proposed method is computationally efficient in practice (see Table 2). Additionally, we have conducted a thorough theoretical analysis for the case $r > 2$, which is detailed in Lines 203–338 on pages 4–7 of the revised version. In this more general setting, we obtain the same consistent results as before.

• ``Although the subsampling technique seems feasible, it lacks a thorough theoretical foundation. For example, in Remark 2, the justification only considers asymptotics, and I suspect finite sample considerations could introduce additional variance. ''

Answer: Thank you for your suggestions, which have significantly enhanced our work. After we updated the theory for general $r$, we find a new approach for computing to replace the original sampling methods in Lines 338-381 on pages 7-8. Specifically, we approximate simple paths with paths allowing node repetition and prove their asymptotic equivalence, ensuring that this approximation does not affect any theoretical results. This approach enables our method to achieve a time complexity matching that of matrix multiplication ($O(n^{2.373})$), eliminating the need for subsampling techniques. In practice, we also suggest eliminating paths with repetitive nodes subject to the computational constraint to improve the finite sample performance. In fact, when $r$ is small, the extra computational cost of this correction can be small, as confirmed by our simulation. Additionally, we have included a simulation result in section 4 for $r=3$. With this new approach, subsampling is no longer required. Once again, thank you for your insightful question, as it motivated us to develop this new method.

• ``Somewhat Unfair Experimental Comparisons: In their experiments, the authors compare their method to more generic methods that work under broader assumptions (not strictly rank 1 or 2). For example, USVT can adapt to various low-rank situations. I consider this comparison somewhat unfair, as the proposed method is tailored to specific ranks in the application, effectively using additional information. Additionally, the authors note that methods are used with default values. Spectral methods with low-rank information might perform as efficiently as the proposed method.''

Answer: Thank you for your suggestion. We appreciate the reviewer's insightful comment regarding the fairness of the comparison between methods. We conducted additional simulations where we used the true $r$ for USVT (i.e., retaining only the first $r$ eigenpairs). We report the results in Table 2 on page 9. Overall, the conclusion remains, namely, that our method performs on par with USVT for estimating the connection probability matrix. This is also expected, as both methods achieve the minimax optimal for estimating the connection probability matrix. Notably, our method offers a slight advantage in computational speed. Finally, we would like to emphasize that our method is capable of estimating the graphon function $f(u,v)$ over $[0,1]\times [0,1]$, and USVT cannot estimate the graphon function. We also compare with other estimating methods for the graphon function, such as the network histogram method proposed by Olhede and Wolfe, 2014. The network histogram algorithm relies on a greedy search to find the maximum likelihood estimator, leading to both high theoretical computational complexity and significant practical computational cost. Furthermore, our algorithm demonstrates superior accuracy with respect to both error metrics compared to the network histogram method, as shown in Table 2.

• ``Limited Applicability: Limiting the methodology to rank-1 and rank-2 graphons restricts its applicability. Additionally, for rank-2 graphons, Assumption 1 feels restrictive.''

Answer: Thank you for your nice comments! In the revision, we have extended our method to address the model with a general $r$. For Assumption 1, admittedly, our condition rules out $\int_0^1G_k(u)du=0$ for some $k$. We note that similar assumptions has made in other works, such as in Theorem 3 of [1]. Exploring how to relax this condition is an important and influential future work.

[1] Bickel, P. J., Chen, A., & Levina, E. (2011). The method of moments and degree distributions for network models. The Annals of Statistics, 39(5), 2280–2301.

• ``How they estimate the eigenvectors for ranks larger than one? Did they use the deflation technique?''

Answer: Thank you for your question. We confirm that we have employed the deflation technique, as outlined in Stoer et al. (1980). To clarify this approach, we have added the following remark in the appendix, where we describe how we estimate the eigenpairs for ranks greater than one, with a particular emphasis on the deflation technique.

``Suppose we want to estimate the leading $r$ eigenpairs of a matrix $A$. The power iteration algorithm can be used by iteratively applying the power method and deflating the matrix. Beginning with $A_0 = A$, each step $k$ (from $1$ to $r$) involves estimating the dominant eigenpair $(\hat \lambda_k, \hat v_k)$ of $A_{k-1}$ using the power iteration method (with normalization of $\hat v_k$) and then deflating the matrix via $A_k = A_{k-1} - \hat \lambda_k \hat v_k \hat v_k^\top$. After $r$ steps, the estimated eigenpairs are used to reconstruct an approximation $\tilde P = \sum_{k=1}^r \hat \lambda_k \hat v_k \hat v_k^\top$. Finally, the estimated probability connection matrix $\hat P$ is obtained by applying an element-wise thresholding operation to $\tilde P: \hat p_{ij} = 1 \wedge (0 \vee \tilde p_{ij})$, ensuring entries fall within the $[0, 1]$ probability range. The whole algorithm is

Algorithm: Power Iteration for $ASG(r)$.

Input: adjacency matrix $A$, rank $r$, maximum number of iterations $N$, tolerance $\epsilon$

Output: estimated probability connection matrix $\hat{P}$.

1. Initialize $A_0 \gets A$.

2. For $k = 1$ to $r$:

3. Initialize vector $\mathbf x_0$ of length $n$ with all ones. Normalize $\mathbf x_0 \gets \frac{\mathbf x_0}{\|\mathbf x_0\|_2}$.

4. For $i = 1$ to $N$:

5. $\mathbf x_{i} \gets A_{k-1} \mathbf x_{i-1}$.

6. Normalize $\mathbf x_{i} \gets \frac{\mathbf {x_{i}} }{\|\mathbf x_{i}\|_2}$.

7. If $\|\mathbf x_{i} - \mathbf x_{i-1}\|_2 < \epsilon$: break; End If

8. End For

9. With deflation technique, compute eigenvalue $\hat \lambda_{k} \gets \mathbf x_i^T A_{k-1} \mathbf x_i$. Set eigenvector $\mathbf{\hat v_k} \gets \mathbf{x_i}$. Update $A_k \gets A_{k-1} - \hat \lambda_k \mathbf{\hat{v}_k} \mathbf {\hat{v}_k}^\top$.

10. End For

11. Compute $\tilde P = \sum_{k=1}^r \hat \lambda_k \mathbf{\hat v_k} \mathbf{\hat v_k}^\top$, and take $\hat p_{ij} = 1 \wedge (0 \vee \tilde P_{ij})$ for $1\le i\neq j\le n, \hat p_{ii}=0$ for $1\le i\le n$.

12. Return estimated probability connection matrix $\hat P =( \{\hat p_{ij}\})_{1 \leq i, j \leq n}$. ''

Thank you again for your suggestion. We will include these statements in the appendix of the final version.

• ``Did they apply the same clipping to 1 in the other estimators they tested, apart from the one they proposed? In my initial comment, I pointed out: "Discounting the clipping on 1, which one can always do knowing that the matrix to be estimated has entries below 1". Given that their proposed estimator explicitly leverages this additional information (that the matrix has entries in $[0,1]$), it seems that any general-purpose algorithm could also be adapted to incorporate this knowledge.''

Answer: Thank you for your question. We would like to clarify that in our current version, all the methods (both ours and others') include this adjustment, i.e., using $1 \wedge (0 \vee \hat{p}_{ij})$ to calculate the final MSE and maximum error. We have added the statement at the beginning of Section 4.

We regret being rejected with a score of 666—reflecting all positive scores—after substantial improvements. We fully respect the decision, however, there are some points we need to clarify:

1. During the rebuttal period, we extended our work to cover cases with arbitrary fixed $r$. The algorithms and theory are, of course, firmly grounded in the methods and theoretical framework originally submitted for $r\le 2$. Therefore, the correctness is guaranteed. Furthermore, the reviewers acknowledged this extension without raising any concerns about its correctness. Finally, considering that the ICLR website clearly states that "reviewers are not required to read the appendix", I have reservations about the statement "the changes warrant a new review of the paper for correctness".

2. Our method and the spectral method differ significantly in motivation, assumptions, and performance in sparse scenarios, as clearly laid out in Remark 3. Moreover, recovering the graphon function from an estimate of the connection probability matrix is, indeed, far from trivial. First, the latent variables associated with the graphon are both unknown and unordered, which prevents the connection probability matrix from being treated as a lattice sampling of the underlying graphon. Second, further analyzing the estimated connection probability matrix becomes complex, particularly when the goal is to achieve sup-norm consistency for graphon estimation. It is truly a pity that there are still misunderstandings regarding this problem, and, somehow, regarding our methods and theoretical framework.