{'fraction_non_alphanumeric': 0.05124445447026092, 'fraction_numerical': 0.01915557560718851, 'mean_word_length': 4.0959862055752465, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 32, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical inference of random graphs with a surrogate likelihood function

Spectral estimators have been broadly applied to statistical network analysis but they do not incorporate the likelihood information of the network sampling model. This paper proposes a novel surrogate likelihood function for statistical inference of a class of popular network models referred to as random dot product graphs. In contrast to the structurally complicated exact likelihood function, the surrogate likelihood function has a separable structure and is log-concave yet approximates the exact likelihood function well. From the frequentist perspective, we study the maximum surrogate likelihood estimator and establish the accompanying theory. We show its existence, uniqueness, large sample properties, and that it improves upon the baseline spectral estimator with a smaller sum of squared errors. A computationally convenient stochastic gradient descent algorithm is designed for finding the maximum surrogate likelihood estimator in practice. From the Bayesian perspective, we establish the Bernstein--von Mises theorem of the posterior distribution with the surrogate likelihood function and show that the resulting credible sets have the correct frequentist coverage. The empirical performance of the proposed surrogate-likelihood-based methods is validated through the analyses of simulation examples and a real-world Wikipedia graph dataset. An R package implementing the proposed computation algorithms is publicly available at https://fangzheng-xie.github.io./publication/ .


Introduction
In the contemporary world of data science, network data are pervasive in a broad range of applications such as sociology (Lacetera et al., 2016;Young and Scheinerman, 2007), econometrics (Mele, 2017;Mele et al., 2022), and neuroscience (Tang et al., 2019). Statistical network analysis is also an interdisciplinary area of research connected with many other fields, including computer science, machine learning, combinatorics, applied mathematics, and physics. To model and analyze network data, various random graph models have been proposed in the literature, including the Erdös-Rényi graph (Erdős et al., 1960), the stochastic block model (Holland et al., 1983), and the latent space model (Hoff et al., 2002).
In this paper, we focus on random dot product graphs (Young and Scheinerman, 2007), a class of random graph models that are popular due its simple architecture and flexibility. On one hand, the edge probability matrix of a random dot product graph has a low-rank structure, which motivates, among others, the use of spectral methods in statistical network analysis. On the other hand, the random dot product graph model is quite flexible because it not only encompasses the popularly used stochastic block models (Holland et al., 1983;Abbe et al., 2016) and their offspring (Airoldi et al., 2008;Binkiewicz et al., 2017;Lyzinski et al., 2017;Sengupta and Chen, 2018), but can also approximate general latent position graphs when the rank of the edge probability matrix grows with the number of vertices at a certain rate (Gao et al., 2015;Tang et al., 2013).
Because the adjacency matrix has a low expected rank, spectral quantities such as the leading eigenvectors of the adjacency matrix and those of the normalized Laplacian matrix, have been extensively used for lowrank random graph inference. In particular, it is well known that the rows of these eigenvectors encode the cluster membership information when the underlying graph is generated from a stochastic block model (Abbe et al., 2020;Lyzinski et al., 2014;Lei and Rinaldo, 2015;Rohe et al., 2011;Sussman et al., 2012).
There has been substantial recent development on the theory for spectral methods and the corresponding subsequent inference tasks in random dot product graphs. For an incomplete list of reference, see Athreya et al. (2016); Sussman et al. (2014); Sarkar and Bickel (2015); Tang and Priebe (2018); Tang et al. (2013Tang et al. ( , 2017a. The readers are also referred to the survey paper Athreya et al. (2017) for a review of the recent advances in this topic.
It has been pointed out (Xie andXu, 2020, 2021;Xie, 2022) that, although the spectral methods for random dot product graphs have gained marvelous success and been broadly applied, the Bernoulli likelihood information contained in the graph distribution has been neglected. This observation has motivated the development of likelihood-based inference for random dot product graphs. Xie and Xu (2020) proposed a fully Bayesian approach for estimating the latent positions in random dot product graphs, referred to as posterior spectral embedding, and established its global minimax optimality. Xie and Xu (2021) proposed a novel one-step procedure, which lead to a one-step estimator that took advantage of the Bernoulli likelihood information of the sampling model through the score function and the Fisher information matrix, to estimate random dot product graphs from the frequentist perspective. There, the authors further established the asymptotic efficiency of the one-step estimator and its smaller asymptotic sum of squared errors compared to that of the spectral estimators. Later, Tang et al. (2022) applied the idea of the one-step refinement of spectral methods to stochastic block models when the block probability matrix is rank deficient.
Despite the success of the one-step estimator, a central question regarding likelihood-based inference for random dot product graphs remains open: What is the behavior of the frequentist maximum likelihood estimator? Also, a related question is: What is the behavior of the Bayes estimator? Efforts attempting to address these two questions aim to gain deeper insight into the likelihood-based inference for random dot product graphs from the frequentist and the Bayesian perspective, respectively. These two questions are also closely related through the Bernstein-von Mises phenomenon (see, for example, Section 10.2 in Van der Vaart, 2000). Here, the major technical barrier is the complicated structure of the parameter space for the latent positions. In this paper, we partially answer the aforementioned two questions by resorting to a cleverly-designed surrogate likelihood function that simplifies the parameter space enormously. Our work features the following novel contributions: Firstly, the surrogate likelihood function has a separable structure, is log-concave, and the associated parameter space for the latent positions is a convex relaxation of the original latent space. These features greatly facilitate both the theoretical analyses and the related practical computations. Secondly, we establish the existence, uniqueness, and the asymptotic efficiency of the frequentist maximum surrogate likelihood estimator under the minimal sparsity condition. In particular, similar to the one-step estimator, the maximum surrogate likelihood estimator improves upon the baseline spectral estimators with a smaller sum of squared errors. Thirdly, we design a computationally efficient stochastic gradient descent algorithm for the maximum surrogate likelihood estimator with adaptive step sizes. Fourthly, regarding the Bayes procedure, we establish the Bernstein-von Mises theorem for the posterior distribution with the surrogate likelihood function and show that the resulting credible sets have the correct frequentist coverage probabilities.
The remaining part of the article is structured as follows. In Section 2, we review the background of random dot product graphs and introduce the surrogate likelihood function. Section 3 elaborates on the theoretical properties and the computational algorithm of the maximum surrogate likelihood estimation. In Section 4, we establish the large sample properties of the Bayes procedure with the surrogate likelihood function. Section 5 demonstrates the empirical performance of the proposed methods through simulation examples and the analysis of a real-world Wikipedia network dataset. We conclude the paper with a discussion in Section 6.
Notations: Let [n] denote the set of consecutive integers from 1 to n: [n] = {1, . . . , n}. The symbol δ means an inequality up to a constant depending on δ, that is, a δ b if a ≤ C δ b for some constant C δ > 0 depending on δ; a similar definition also applies to the symbol δ . The notation x denotes the Euclidean norm of a vector x = [x 1 , . . . , x d ] T ∈ R d , that is, x = ( d k=1 x 2 k ) 1/2 . The d × d identity matrix is denoted as I d . The notation O(n, d) = {U ∈ R n×d : U T U = I d } denotes the set of all orthonormal d-frames in R n , where d ≤ n, and we write O(d) = O (d, d). For a matrix X = [x ik ] n×d , σ k (X) denotes its kth largest singular value, and when X is square and symmetric, λ k (X) denotes its kth largest eigenvalue in magnitude.
Matrix norms with following definitions are used: the spectral norm X 2 = σ 1 (X), the Frobenius norm ik ) 1/2 , the matrix infinity norm X ∞ = max i∈[n] d k=1 |x ik |, and the two-to-infinity norm X 2→∞ = max i∈[n] ( d k=1 x 2 ik ) 1/2 . In particular, these norm notations apply to any Euclidean vector x ∈ R d viewed as a d × 1 matrix. Given two symmetric positive semidefinite matrices A, B of the same dimension, we write A B (A B, respectively) if B − A (A − B, respectively) is positive semidefinite.
2 Background and the Surrogate Likelihood 2.1 Background on random dot product graphs We begin by briefly reviewing the background on random dot product graphs and adjacency spectral embedding. Consider a graph with n vertices labeled as [n] = {1, . . . , n}. Let X be a subset of R d such that x T 1 x 2 ∈ (0, 1) for all x 1 , x 2 ∈ X , where d is fixed and d ≤ n, and let ρ n ∈ (0, 1] be a sparsity factor. Each ver-tex i ∈ [n] is associated with a vector x i ∈ X , referred to as the latent position for vertex i. We say that a symmetric random matrix A = [A ij ] n×n ∈ {0, 1} n×n is an adjacency matrix generated by a random dot product graph with latent position matrix X = [x 1 , . . . , x n ] T and sparsity factor ρ n , denoted by A ∼ RDPG(ρ The sparsity factor ρ n fundamentally controls the overall average graph expected degree through nρ n when the entries of XX T are bounded away from 0 and ∞. Remark 1 (Deterministic versus stochastic latent positions). In this work, we consider the setup where the latent positions x 1 , . . . , x n are deterministic parameters to be estimated. Another slightly different modeling approach is to consider x 1 , . . . , x n as independent and identically distributed latent random variables (see, for example, Athreya et al., 2016;Tang et al., 2017b;Tang and Priebe, 2018). This random formulation of the latent positions introduces implicit homogeneity and is connected to the infinite exchangeable random graphs (Janson and Diaconis, 2008). The same homogeneity condition was retained in Xie and Xu (2021) using a Glivenko-Cantelli type condition when x 1 , . . . , x n are deterministic. The latter Glivenko-Cantelli type condition is also relaxed in the current work as we only require that σ d (X) > 0 (see Remark 2 below).
Remark 2 (Nonidentifiability). The latent position matrix X is not uniquely identified in the following two senses. Firstly, any low-rank positive semidefinite connection probability matrix P = XX T can have different factorizations because for any orthogonal matrix W ∈ O(d), XX T = (XW)(XW) T . Secondly, for any d > d and any latent position matrix X ∈ R d , there exists another matrix X ∈ R d such that XX T = X (X ) T . The latter source of non-identifiability can be removed by requiring that σ d (X) > 0, while the former source is inevitable without further constraints. Thus, any estimator of the latent position matrix X can only recover it up to an orthogonal transformation.
Example 1 (Stochastic block model). Random dot product graphs have connections with the popular stochastic block model (Holland et al., 1983). Consider a graph with n vertices that are partitioned into K communities, where K is assumed to be much smaller than n. Let τ : [n] → [K] be a cluster assignment function that assigns each vertex to a unique community. Let B = [B kl ] K×K ∈ (0, 1) K×K be a symmetric probability matrix and A ij be the binary indicator of the existence of an edge between vertices i and j. Then the stochastic block model specifies that A ij ∼ Bernoulli(B τ (i)τ (j) ) independently for all i, j ∈ [n], i ≤ j, and A ij = A ji for all i > j. By converting the community assignment to a matrix Z = [1{τ (i) = k}] n×K , we see that the expected adjacency matrix ZBZ T is symmetric and of low rank. Furthermore, if B is positive semidefinite with rank d ≤ K and can be factorized as B = VV T for a K × d matrix V, then A can be seen as an adjacency matrix generated by the random dot product graph with latent position matrix X = ZV, Motivated by the low-rank structure of random dot product graphs, Sussman et al. (2012) proposed to estimate the latent position matrix X by solving the least squares problem X = arg min X∈R n×d A−XX T 2 F . The interpretation is that X X T can be viewed as the projection of the data matrix A onto the space of all n × n rank-d positive semidefinite matrices with regard to the Frobenius norm distance. The solution X is referred to as the adjacency spectral embedding of A into R d , and can be computed as the matrix of eigenvectors associated with the top d eigenvalues of A, scaled by the square roots of the corresponding eigenvalues (Eckart and Young, 1936). The asymptotic properties of X have been established in the literature Athreya et al., 2016;Tang and Priebe, 2018). Notably, Athreya et al. (2016), Tang andPriebe (2018), and Xie and Xu (2021) have shown that each row of the adjacency spectral embedding converges to a mean-zero multivariate normal distribution after appropriate standardization.

The surrogate likelihood function
In this subsection, we derive the surrogate likelihood function for the random dot product graph model.
The motivation is that the exact likelihood function has a complicated structure, bringing challenges for developing the theory of the maximum likelihood estimation. The difficulty partially comes from the fact that the random dot product graph model belongs to a curved exponential family, and the theory of the maximum likelihood estimation is much more difficult in curved exponential families than in the canonical ones (see, for example, Section 2.3 in Bickel and Doksum, 2007).
To distinguish a generic latent position x i ∈ R d and its true value associated with the data generating distribution, let x 0i denote the ground truth of x i , i ∈ [n], and X 0 = [x 01 , . . . , x 0n ] T . Let us begin by considering the log-likelihood function of a single x i when the remaining latent positions (x 0j ) j =i are accessible: (2.1) We refer to 0in (x i ) in (2.1) as the oracle log-likelihood function because it requires the true values of the remaining x j 's with j = i. Theorem 2 in Xie and Xu (2021) established the consistency and asymptotic normality of the maximizer of the oracle log-likelihood function 0in (x i ) in (2.1). Nevertheless, the oracle log-likelihood is not computable because (x 0j ) are not accessible in practice. Following the idea in Xie and Xu (2021), we consider replacing the unknown latent positions by the corresponding rows of the adjacency spectral embedding. Formally, let x j be the jth row of the adjacency spectral embedding X, j ∈ [n]. Then we obtain the following approximation to the oracle log-likelihood: Note that the last term in 0in is replaced by A ii log(ρ for convenience, which is immaterial. This approximation step is motivated by the uniform consistency of the adjacency spectral embedding: There exists a d × d orthogonal W such that XW − ρ 1/2 n X 0 2 = O{ (log n)/n} with high probability (Lyzinski et al., 2014;Xie, 2022).
With the approximation in (2.2), the constraints for the latent position x i become a system of linear . Geometrically, these constraints correspond to a convex polyhedron. Namely, given any vector x i ∈ R d , checking whether x i is in such a convex polyhedron requires O(n) operations, so that the relevant computation could be cumbersome. We now resolve this issue by applying a quadratic Taylor approximation to the terms log(ρ 1/2 n x T i x j ) and relax the parameter space for x i .
Here we can drop the sparsity factor ρ n without loss of generality. Formally, write g ij (x i ) = A ij log(x T i x j ). Then a quadratic Taylor approximation to g ij at x i = x i leads to Hence, ignoring the constant terms that are free of x i , combining the approximations in (2.2), (2.3), and (2.4) leads to the following surrogate log-likelihood function (2.5) Therefore, by Cauchy-Schwarz inequality, the the parameter space for x i associated with the surrogate likelihood of vertex i can be taken as the unit ball {x i ∈ R d : x i ≤ 1} for all i ∈ [n] when max j∈[n] x j 2 < 1 (which holds with high probability). Consequently, we relax the original complicated parameter space } to a simple unit ball, which is much more tractable to work with. Moreover, the surrogate likelihood function has a separable structure because in (x i ) does not involve x j for j = i. This convenience enables parallelization when related computation is requested. In addition, a simple algebra shows that the surrogate likelihood function is log-concave, a highly desired feature when optimization and Monte Carlo sampling are needed.

Comparison with the one-step estimator
Recently, Xie and Xu (2021) proposed a one-step estimator X OS = [ x OS 1 , . . . , x OS n ] T for random dot product graphs that improves upon the adjacency spectral embedding: (2.6) The one-step estimator originates from a one-step updating scheme of the Newton-Raphson method for maximizing the log-likelihood function with the initial guess being the adjacency spectral embedding (see, e.g., Section 5.7 in Van der Vaart, 2000). It is clear from the construction that the one-step estimator takes advantage of the likelihood information of the sampling distribution through the Fisher information matrix and the score function.
In Section 2.2, we have shown the derivation of the surrogate log-likelihood function by applying a quadratic Taylor approximation to the logarithm function log(x T i x j ). We now show that the same approximation treatment applied to the entire function in (2.2) results in an approximate log-likelihood function whose maximizer is exactly the one-step estimator. Formally, applying a second-order Taylor expansion to Following the idea in (2.4), we can also conceive the following approximation: . (2.8) We thus obtain the following quadratic approximation to (2.2) modulus a constant term from (2.5), (2.7), and (2.8): . (2.9) Then a simple algebra shows that the one-step estimator x Clearly, the surrogate log-likelihood function in (2.5) is constructed by applying a Taylor expansion to the term log(x T i x j ), whereas the one-step estimator is obtained by applying the Taylor expansion to the entire function. Thus, intuitively, the surrogate likelihood retains more likelihood information than the one-step procedure does. Below, we visualize this heuristic using a toy numerical example.
Example 2. Consider the following random dot product graph model. Let n = 300, (t i ) n i=1 be equidistant points over [0, 1], x 0i = 0.2+0.6 sin(πt i ), i ∈ [n], and X 0 = [x 01 , . . . , x 0n ] T ∈ R n×1 . Suppose A ∼ RDPG(X 0 ) and we focus on the likelihood function for x i with i = 100. Figure 1 visualizes the comparison among the oracle log-likelihood 0in (x i ), the surrogate log-likelihood in (x i ), and the approximate log-likelihood (OS) in (x i ) associated with the one-step estimator. The constant terms of these functions have been added to make them comparable. The vertical lines mark the maximizers of the three functions, respectively. It is visually clear that the maximizer of the surrogate log-likelihood estimate is closer to that of the oracle log-likelihood than the one-step estimate is, suggesting that the maximum surrogate likelihood estimator may outperform the one-step estimator in some practical finite sample problems.

Theoretical properties
This subsection elaborates on the theoretical properties of the frequentist inference with the surrogate likelihood. Below, Theorem 3.1 establishes the existence and uniqueness of the maximum surrogate likelihood estimator.
Theorem 3.1. Suppose A ∼ RDPG(ρ 1/2 n X 0 ) and (log n)/(nρ n ) → 0 as n → ∞. Assume λ d (X T 0 X 0 /n) ≥ λ for some constant λ > 0 for all n > d, and min i,j∈ [n]   in (x i ) associated with the one-step estimator. The three vertical lines mark the one-step estimate, the maximum surrogate likelihood estimate, and the oracle maximum likelihood estimate, respectively. δ > 0. Let i ∈ [n] be a fixed vertex and consider the maximum surrogate likelihood estimator x i = arg max xi: xi 2≤1 in (x i ). Then for any c > 0, there exists some constant N c,δ,λ ∈ N + depending on c, δ, λ such that P 0 ( x i exists and is unique) ≥ 1 − n −c for all n ≥ N c,δ,λ .
Remark 3 (Sparsity condition). The sparsity condition that (log n)/(nρ n ) → 0 required in Theorem 3.1 and in the asymptotic normality of Theorem 3.2 is minimal in the following sense. It is well known that the random adjacency matrix A no longer concentrates around its expected value E 0 (A) when (log n)/(nρ n ) → ∞ (Tang and Priebe, 2018). Furthermore, Abbe et al. (2020) and Xie (2022) showed that in order to have ρ −1/2 n XW − X 0 2→∞ = o(1) with high probability, which is an indispensable ingredient in our employed proof, it is necessary that (log n)/(nρ n ) → 0.
Remark 4 (Comparison with the adjacency spectral embedding and the one-step estimator). Athreya et al. (2016), Tang and Priebe (2018), and Xie and Xu (2021) have establish the large sample properties of the adjacency spectral embedding and the one-step estimator as the following. Let X = [ x 1 , . . . , x n ] T and in distribution. Also, . Theorem 3.2 thus suggests that the maximum surrogate likelihood estimator improves upon the adjacency spectral embedding and is (first-order) asymptotically equivalent to the one-step estimator. This phenomenon is also known as the local efficiency (Xie and Xu, 2021) because the asymptotic covariance matrix for a single latent position x i is the same as that of the oracle maximum likelihood estimator.

Computation details
This subsection discusses the detailed algorithm for computing the maximum surrogate likelihood estimator.
For a given vertex i ∈ [n], the estimator x i = arg max xi in (x i ) can be computed separately for each vertex i ∈ [n]. Thus, it is sufficient to design an algorithm for solving the optimization problem max xi 2≤1 (1/n) in (x i ).
Observe that the objective function (1/n) in (x i ) is concave and can be written in a sample average fashion, which motivates us to adopt the stochastic gradient descent algorithm (Robbins and Monro, 1951). Let j 1 , . . . , j s be independent Unif(1, . . . , n) random variables, where s ∈ {1, . . . , n} is the so-called batch size, and for any j ∈ [n], let can be viewed as a noisy measurement of the objective function Then given a sequence of step sizes {α t } t≥1 and a initial guess s )} t≥1 are independent copies of (j 1 , . . . , j s ). The advantage of the stochastic gradient descent method over the classical gradient descent algorithm is that, with a comparatively small batch size s, one only needs to compute s gradient measurements of m i (x i , j) rather than all the gradient measurements . This computational convenience is especially desired when the network contains large number of vertices. To implement the algorithm with adaptive step sizes, we follow the suggestion given by Duchi et al. (2011) and Li and Orabona (2019) and take where a 0 , b 0 > 0 and 0 < ≤ 1/2 are constants.
The key difference between our algorithm and the standard stochastic gradient descent algorithm is that shows the convergence of such a step-halving stochastic gradient descent with adaptive step sizes.
For sufficiently large n, the constraint is satisfied by requiring that x i 2 ≤ 1 since the adjacency spectral embedding X = [ x 1 , . . . , x n ] T satisfies max j∈[n] x j 2 < 1 with high probability. However, this requirement may not hold in certain finite sample problems, in which case the surrogate log-likelihood function in (x i ) is no longer well-defined. This numerical issue can be practically addressed by the following smooth concatenation technique. Roughly speaking, for a fixed j ∈ [n], when 1 − x T i x j drops below a small threshold, we replace the objective function (1/n) in (x i ) by a quadratic function such that the two pieces of functions are concatenated smoothly. Formally, let τ > 0 be a small threshold and define Then the objective function (1/n) in (x i ) can be replaced by (1/n) n j=1 h i (x i , j) and the aforementioned stochastic gradient descent algorithm applies with ∂m i (x i , j)/∂x i replaced by ∂h i (x i , j)/∂x i .

Bayesian Estimation With Surrogate Likelihood
This section explores Bayesian estimation of random dot product graphs with the proposed surrogate likelihood. Although Xie and Xu (2020) has established the minimax optimality of the Bayesian random dot product graph model with the exact likelihood, the asymptotic shape of the posterior distribution is yet to be characterized because of the complicated structure of the exact likelihood function. In contrast, thanks to the separable and log-concave properties of the surrogate likelihood, we are able to completely characterize the asymptotic posterior distribution of the latent positions with the exact likelihood replaced by the surrogate. Formally, for any fixed vertex i ∈ [n] and a prior distribution π(·) supported on {x ∈ R d : x 2 ≤ 1}, the posterior distribution of x i given A with the surrogate log-likelihood function in (x i ) can be written as Then the joint posterior density of the entire latent position matrix X = [x 1 , . . . , x n ] T is taken as the When the exact likelihood function is not available or intractable for analysis or computation, the idea of using a general statistical criterion function to replace the likelihood in the Bayes formula is not entirely new, among which an influential work is Chernozhukov and Hong (2003). There have also been several recent works addressing the large sample properties of the so-called quasi-posterior or Gibbs posterior distributions (Kleijn and van der Vaart, 2012;Miller, 2021;Martin, 2018, 2022). One key difference is that unlike the well-specified exact posterior distributions, the frequentist coverage of the credible sets of the quasi-posterior distributions may not agree with their credibility level (Kleijn and van der Vaart, 2012).
Below, we show that, with the surrogate likelihood, the posterior distribution produces credible sets that have the correct frequentist coverage. This is achieved through the following Bernstein-von Mises theorem.
Theorem 4.1. Suppose the conditions of Theorem 3.1 hold and the embedding dimension d is fixed. Let π(·) be a prior density satisfying c ≤ π(x i ) ≤ C and |π(x) − π(y)| ≤ C x − y 2 for any x, y with x 2 , y 2 ≤ 1 for some constants 0 < c, C, C < ∞. Let W be the d × d orthogonal matrix in Theorem 3.2. For any fixed be the density of t induced from (4.1). Then for any α > 0, in probability.
Below, Corollary 4.1 discusses the effect of Theorem 4.1 on subsequent inference. In particular, it shows that for each vertex i ∈ [n], the posterior mean has the same asymptotic distribution as the maximum surrogate likelihood estimator, and the asymptotic level-α credible set has the correct frequentist coverage probability.
dx i as the posterior mean and the posterior covariance matrix of x i , respectively. Let q 1−α be the (1 − α) quantile of the χ 2 d distribution and C in (α) = {x i : in distribution, and In practice, the posterior distribution based on the surrogate likelihood can be computed using a standard Metropolis-Hastings algorithm with parallelization over the vertices i ∈ [n]. The detailed algorithm is provided in the Supplementary Material. Note that in practice, we can also apply the smooth concatenation technique discussed in Remark 5 to the posterior computation by simply replacing the surrogate log-likelihood 5 Numerical Examples

A latent curve example
In this subsection, we study the empirical performance of the proposed estimation procedures through a simulated random dot product graph example, where the latent positions are generated from a onedimensional curve. Consider a random dot product graph with n vertices and latent dimension d = 1.
For each vertex i ∈ [n], the latent position x 0i for the ith vertex is set to x 0i = 0.8 sin{π(i − 1)/(n − 1)} + 0.1. Let X 0 = [x 01 , . . . , x 0n ] T , n = 1000. Given A ∼ RDPG(X 0 ), we consider the following four estimation procedures for X 0 : the adjacency spectral embedding, the one-step estimate, the maximum surrogate likelihood estimate obtained using the step-halving stochastic gradient descent algorithm, and the Specifically, for the first objective, given a generic estimateX for X 0 , we use the sum of squared errors inf W∈{±1} X W − X 0 2 F as the evaluation metric. For the second objective, we compute the vertex-wise asymptotic 95% frequentist confidence intervals and Bayesian credible intervals. The vertex-wise 95% confidence intervals based on the maximum surrogate likelihood estimate are computed as follows: Denote the 1 − α/2 quantile of the standard normal distribution by z 1−α/2 . Then by Theorem 3.2, for each i ∈ [n],  For the first objective, the sum of squared errors of the estimates are shown in Table 1. We can see that the sum of squared errors of the adjacency spectral embedding is comparatively larger than those of the remaining competitors, while the likelihood-based estimates have smaller sum of squared errors. This phenomenon empirically validate the conclusion that the likelihood-based estimates, namely, the one-step estimate, the maximum surrogate likelihood estimate, and the Bayes estimate, improve upon the the spectralbased adjacency spectral embedding.

A stochastic block model example
We now consider a stochastic block model in the context of a random dot product graph. The latent dimension is d = 2, the number of communities is K = 5, and the unique latent positions are The cluster assignments of the vertices (z i ) n i=1 are drawn from a categorical distribution with probability vector [1/K, . . . , 1/K] T and we set x 0i = v zi , i ∈ [n]. Note that v 3 is very close to the boundary of the parameter space. Let X 0 = [x 01 , . . . , x 0n ] T and suppose an adjacency matrix A is generated from RDPG(X 0 ). The number of vertices n is 2000.
We consider the performance of the same estimates as in Section 5.1 given a realization A ∼ RDPG(X 0 ): the adjacency spectral embedding, the one-step estimate, the maximum surrogate likelihood estimate and the Bayesian estimation with the surrogate likelihood. For the maximum surrogate likelihood estimate, we implement the step-halving stochastic gradient descent algorithm with the batch size set to s = 500 and s = n (giving rise to the classical gradient descent algorithm) to compare the computational costs. The setup of the Metropolis-Hastings sampler for the Bayesian estimation is the same as in Section 5.1, and the convergence diagnostics are provided in the Supplementary Material. We take the posterior mean as the point estimate as before. The same experiment is repeated for 1000 independent Monte Carlo replicates.
Similar to Section 5.1, given a generic estimateX, we compute the sum of squared errors of the estimates inf W∈O(2) X W−X 0 2 F to measure the estimation accuracy. These results are summarized in Table 2. We see that the one-step estimate is numerically unstable because v 3 is close to the boundary of the parameter space.
Overall, the Bayes estimate outperforms the other competitors with the least errors, while the adjacency spectral embedding and the maximum surrogate likelihood estimate have similar performance in terms of the estimation error. This phenomenon suggests that, when some latent positions are close to the boundary of the parameter space, the Bayesian estimation method based on the Markov chain Monte Carlo sampler is numerically more stable than the optimization-based frequentist adjacency spectral embedding and the maximum surrogate likelihood estimate.
The computation times of a single experiment for different estimation procedures are reported in Table   2. We see that the adjacency spectral embedding and the one-step estimate are faster to compute, whereas the maximum surrogate likelihood estimate obtained through the classical gradient descent algorithm and the Bayes estimate are more computationally expensive. We also observe that the stochastic gradient descent algorithm is significantly faster than the classical gradient descent algorithm for finding the maximum surrogate likelihood estimate and it well balances between the computational efficiency and the estimation accuracy.

Analysis of Wikipedia Graph Dataset
In this section, we apply the proposed surrogate likelihood estimation methods to a real-world Wikipedia graph dataset. The network data is structured as follows: The vertices represent 1382 Wikipedia articles that are connected to the article named Algebraic Geometry within two hyperlinks, and an edge is assigned to link two articles if they are connected by a hyperlink. Besides the network itself, each Wikipedia article is also assigned with one of the following six class labels: people, places, dates, things, math and category.   The dataset is publicly available at at http://www.cis.jhu.edu/~parky/Data/data.html.
The goal is to study the clustering accuracy using different estimates when the embedding dimension varies. Given a selected embedding dimension d ≥ 1, we consider the following four estimates: the adjacency spectral embedding, the one-step estimate, the maximum surrogate likelihood estimate computed using the step-halving stochastic gradient descent algorithm, and the Bayes estimate based on the surrogate likelihood Rather, only the class labels of the vertices are available to us. To this end, we follow the suggestion in Tang and Priebe (2018) and apply the Gaussian-mixture-model-based clustering to the aforementioned four estimates. Namely, these estimates are regarded as the input for learning the clustering structure of the Wikipedia article network. We report the clustering accuracy using the Rand index (Rand, 1971) as the evaluation metric.
The Rand indices of the clustering results using different estimates across different embedding dimensions d ∈ {1, 2, . . . , 10} are shown in Table 3. On one hand, we can see that when d ≤ 2, the adjacency spectral embedding yields better clustering accuracy with a higher Rand index value than the remaining competitors.
On the other hand, as the embedding dimension d increases from 2 to 10, the maximum surrogate likelihood estimate and the Bayes estimate with the surrogate likelihood outperform the other two competitors. A plausible explanation of this phenomenon could be that the eigenvectors of the adjacency matrix with smaller eigenvalues are noisier than the top two eigenvectors, but this source of noise is reduced through the additional information introduced by the surrogate likelihood function.

Discussion
In Section 2.3, we have seen that the one-step estimator also corresponds to the maximizer of an approximate likelihood function, but it has a worse approximation quality than the proposed surrogate likelihood near the oracle maximum likelihood estimator. Surprisingly, under a framework of generalized estimating equations proposed by Xie and Wu (2022), the gradients of both the surrogate log-likelihood function and the approximate log-likelihood function associated with the one-step estimator can be viewed as some generalized estimating equations that take advantage of the likelihood function information. This intuition conforms to the fact that the estimators based on approximation of likelihood are asymptotically equivalent up to the first order. However, we have also found in some finite sample problems that the maximum surrogate likelihood estimator outperforms the one-step estimator. This difference may be caused by the difference in higher order terms, which is an interesting topic that we defer to future research.

Supplementary material
This supplementary material contains the proofs of the main results of the manuscript, the proof of the convergence of the step-halving stochastic gradient descent algorithm, additional implementation details including the algorithm details, an additional simulation example, and the convergence diagnostics of the Metropolis-Hastings sampler.

Supplementary Material for "Statistical inference of random graphs with a surrogate likelihood"
Abstract This supplementary material contains the proofs of the main results of the manuscript, the proof of the convergence of the stochastic gradient descent, and additional implementation details, including the algorithm details, an additional simulation example, and the convergence diagnostics of the Metropolis-Hastings sampler.

A Preliminary Results
Then for all c > 0, there exists some constant N c,λ ∈ N + depending on c, λ, such that for all n ≥ N c,λ , with probability at least 1 − n −c .
This completes the proof.
Lemma A.2 (Some frequently used results). Suppose A ∼ RDPG(ρ 1/2 n X 0 ) and assume the conditions of Theorem 3.1 hold. Denote p ij = x T i x j , i, j ∈ [n]. Then for any c > 0, there exists a constant N c,δ,λ ∈ N + depending on c, δ, λ such that for all n ≥ N c,δ,λ , the following hold with probability at least 1 − n −c : c,λ ρ 1/2 n log n n .
Proof. For the first result, by Lemma A.1 and the condition that log n nρn → 0, we can pick a constant N c,δ,λ ∈ N + depending on c, δ, λ such that for all n ≥ N c,δ,λ , with probability at leat 1 − n −c , For the second result, over the same event as above, we have x j 2 + ρ 1/2 n ) XW − ρ 1/2 n X 0 2→∞ c,λ ρ 1/2 n log n n .
For the third result, over the same event as above, we have ρ n x T 0i x 0j ≤ C c,λ ρ 1/2 n log n n + ρ n (1 − δ).
Lemma A.3 (Concentration of Hessian matrices). Suppose A ∼ RDPG(ρ 1/2 n X 0 ) and assume the conditions of Theorem 3.1 hold. Denote p ij = x T i x j , i, j ∈ [n] and let > 0 be sufficiently small. Then for any c > 0, there exists a constant N c,δ,λ ∈ N + depending on c, δ, λ such that for all n ≥ N c,δ,λ , the following hold with probability at least 1 − n −c : Proof. For simplicity of notation, denote p 0ij = ρ n x T 0i x 0j . The large probability bounds below are with regard to n ≥ N c,δ,λ for some large constant N c,δ,λ depending on c, δ, λ.
We show the first conclusion first. Write For the first term, with probability at least 1 − n −c , ( where in the second inequality we use Lemma A.2, in the third inequality triangle inequality and Cauchy-Schwarz inequality, and in the fourth inequality Lemma A.1 and Lemma A.2.
For the second term, with probability at least 1 − n −c , by Lemma A.2 and the result that A ∞ c nρ n with probability at least 1−n −c which follows from triangle inequality and Bernstein's inequality.
For the third term, for a typical (k, l)th entry, by Bernstein's inequality and a union bound over i ∈ [n], for any t > 0, where K δ > 0 is a constant depending on δ. Taking t = C (ρ 3 n log n)/n for an appropriate constant C > 0, we see that c,δ ρ 3 n log n n with probability at least 1 − n −c . Since d is fixed (it implicitly depends on λ), we have (1 − p 0ij ) 2 2 c,δ,λ ρ 3 2 n log n n with probability at least 1 − n −c .
For the fourth term, with probability at least 1 − n −c , For the fifth term, with probability at least 1 − n −c , c,δ,λ log n nρ n by Lemma A.2. So the first conclusion is shown by combining the above five bounds.
Next, we show the second conclusion. Write The first term is O(ρ 3 2 n log n n ) with probability at least 1 − n −c as previously shown. For the second term, with probability at least 1 − n −c , c,δ,λ 1 n ρ n ρ 1/2 n log n n + ρ 1/2 n log n n ρ n c,δ,λ ρ 3 2 n log n n by Cauchy-Schwarz inequality, Lemma A.2, and For the third term, with probability at least 1 − n −c , x j 2 2 c,δ,λ ρ 3 2 n log n n by Lemma A.2. So the second conclusion is shown by combining the above three bounds.
Lemma A.4. Suppose A ∼ RDPG(ρ 1/2 n X 0 ) and assume the conditions of Theorem 3.2 hold. Denote Proof. Denote The proof of Lemma A.4 here is almost identical to those of Lemma D.5 in Xie and Xu (2021) and Lemma A.6 in Tang et al. (2017a) and we include the proof here only for the sake of completeness. The key idea relies on the following logarithmic Sobolev concentration inequality: Lemma A.5 (Theorem 6.7 in Boucheron et al. (2013)). Let A, A ∈ {0, 1} n×n be two symmetric random adjacency matrices, Z = Z(A) be a measurable function of of A. Denote by A (kl) the adjacency matrix obtained by replacing the (k, l) and (l, k) entries of A by those of A , and Z kl = Z(A (kl) ). If there exists a constant v > 0 such that then for all t > 0, P(|Z − E(Z)| > t) ≤ 2 exp{−t 2 /(2v)} + η.
Let A be another symmetric random adjacency matrix. Denote by A (kl) the adjacency matrix obtained by replacing the (k, l) and (l, k) entries of A by those of A and Z kl = Z(A (kl) ) Observe that A and A (kl) are the same except for the (k, l) and (l, k) entries. Also note that A and A are binary random matrices.
Therefore, when Z − Z kl = 0, Invoking Lemma A.5, we obtain that for all > 0. The proof is thus completed.
Theorem A.6 (Theorem 4.7 in Xie, 2022). Suppose A ∼ RDPG(ρ 1/2 n X 0 ) and assume the conditions of Theorem 3.1 hold. Define the one-step estimator x , and for all c > 0, there exists a constant N c,δ,λ ∈ N + depending on c, δ, λ, such that for all n ≥ N c,δ,λ , with

B Proofs of the Main Results
B.1 Proof of Theorem 3.1

Proof.
We first prove existence. For any c > 0, there exists N c,δ,λ ∈ N + such that x j 2 ≤ ρ n (1 − δ 2 ) < 1 with probability at least 1 − n −c , where the first inequality follows from Cauchy-Schwarz inequality, the second from Lemma A.2. By definition of M in (x i ), it is continuous over the closed unit ball {x i ∈ R d : x i 2 ≤ 1} over this event. Hence the maximizer x i of M in (x i ) exists with probability at least 1 − n −c .
Next we prove uniqueness. By definition, with probability at By Theorem 5.2 in Lei and Rinaldo (2015) and Weyl's inequality, there exists a constant depending on c, such that with probability at least 1 − n −c , Therefore, for any c > 0, there exists N c,δ,λ ∈ N + such that for all n ≥ N c,δ,λ , with probability at least 1 − n −c , M in (x i ) is strictly concave. Hence it has a unique maximizer x i .
We first establish the following consistency result: For any c > 0, there exists some constant N c,δ,λ ∈ N + depending on c, δ, λ such that for all n ≥ N c,δ,λ ∈ N + , there exists an orthogonal matrix Simple calculation shows that For simplicity of notation, in what follows the large probability bounds are with regard to n ≥ N c,δ,λ for some large constant N c,δ,λ depending on c, δ, λ.
Claim I (identifiability): For all > 0, Claim II (uniform convergence): With probability at least 1 − n −c , Now we show Claim I. It is obvious that ∂Min ∂xi (ρ which implies that M in (x i ) is strictly concave with ρ 1/2 n x 0i as a unique maximizer. By Taylor's theorem, Thus Claim I is shown. Now we show Claim II. By triangle inequality, For the second term, log n nρ n with probability at least 1 − n −c . For the third term, For the fourth term, log n nρ n with probability at least 1 − n −c .
In order to bound the first term, a maximal inequality is required. We use the results in Chapter 8 of Kosorok (2008). Define a stochastic process on {y ∈ R d : y 2 ≤ 1} for each k ∈ [d], Then for any y, y with y 2 ≤ 1, y 2 ≤ 1, where y = θy + (1 − θ)y for some θ ∈ [0, 1]. By Hoeffding's inequality, where C δ > 0 is a constant depending on δ, which indicates that J ink (y) is a sub-Gaussian process on {y ∈ R d : y 2 ≤ 1} with respect to the metric d n (y, y ) = y − y 2 √ C δ ρ 2 n /n. The metric entropy of the metric space ({y ∈ R d : y 2 ≤ 1}, d n ) can be bounded by where K δ is a constant depending on δ. Recall that the ψ 2 -Orlicz norm (sub-Gaussian norm) of a random variable X is defined as where ψ 2 (x) = e x 2 − 1 (see Chapter 8 of Kosorok, 2008).
By Theorem 8.4 in Kosorok (2008), where we note that d depends on λ implicitly. Then by Lemma 8.1 in Kosorok (2008) and a union bound over i ∈ [n], max i∈[n] sup xi 2 ≤1 |J ink (x i )| c,δ,λ (ρ 2 n log n)/n with probability at least 1 − n −c . So |J ink (x i )| c,δ,λ ρ n log n n with probability at least 1 − n −c . Thus Claim II is shown.
where the first inequality follows from x i being the unique zero of ∂ M in /∂x i (x i ) 2 with probability at least 1 − n −c , the second inequality from triangle inequality, and the third inequality from Claim II.
We next establish the asymptotic normality. We utilize the asymptotic normality of the one-step estimator x (OS) i (Theorem A.6) to establish the asymptotic normality of the maximum surrogate likelihood estimator x i . By the previous part of the theorem, we know that with probability at least 1 − n −c , x i is in the interior of the closed unit ball , we apply Taylor's theorem to where in the last inequality we applied the fact that A − P 2 c √ nρ n with probability at least 1 − n −c (Theorem 5.2 in Lei and Rinaldo, 2015). By Lemma A.1 and the previous part of the theorem, with probability at least 1 − n −c , So the Taylor expansion of (∂ M in )/(∂x i ) mentioned above can be written as where R in1 ∈ R d×d is a random matrix with R in1 2 c,δ,λ ρ 1/2 n (log n)/n with probability at least 1−n −c .
By definition of M in (x i ) and Lemma A.3, where R in2 ∈ R d×d is a random matrix with R in2 2 c,δ,λ ρ 1/2 n (log n)/n with probability at least 1 − n −c . Similarly as in the proof of Theorem 3.1, i.e., G in is finite and positive definite with probability at least 1 − n −c .
We finally show the convergence of the sum of squares errors, that is By the previous result, we have By Lemma A.4, the first term equals For the second term, by Theorem 4.7 in Xie (2022), we have r in 2 2 c,δ,λ ρ n (log n) 2 n + (log n) 4 nρ n with probability at least 1 − n −c for all n ≥ N c,δ,λ , so n i=1 r in 2 2 = o P0 (1) by the condition that (log n) 4 = o(nρ n ). For the third term, by Cauchy-Schwarz inequality, we have Hence, we conlcude that

B.3 Proof of Theorem 4.1
Proof. Similar to the earlier proofs, the large probability bounds below are with regard to n ≥ N c,δ,λ for some large constant N c,δ,λ depending on c, δ, λ.
Denote the parameter space of t by Θ in = {t ∈ R d : It is sufficient to show that To see this, note that (4.2) in the manuscript can be rewritten as Since (B.2) implies that max i∈[n] |d in − det(2πG −1 0in ) 1/2 π(Wρ 1/2 n x 0i )| = o P0 (1) (by taking α = 0), it can be seen that (B.2) implies that the two terms on the right hand side of the previous display are o P0 (1). Hence, we are left with establishing (B.2).
Let {η n } ∞ n=1 be a sequence to be determined later with 0 < η n → ∞ and consider the following partition of R d : We first consider the integral of (B.2) over A 3 . By definition of 1(t ∈ Θ in ), the integral over A 3 can be bounded by since A 3 is shrinking to empty set and min i∈[n] (G 0in ) ≥ λ has been shown in the proof of Theorem 3.2 (see diaplay (B.1)). We next consider the integral of (B.2) over A 2 . Define the event Note that by Lemma A.2, Theorem 5.2 in Lei and Rinaldo (2015), and Weyl's inequality, with probability This shows that P 0 (E 2n ) ≥ 1 − n −c for all n ≥ N c,δ,λ . By Taylor's expansion, for any t ∈ Θ in , we have (1 + t α 2 ) exp (1 + t α 2 ) e −t T G0int/2 π(ρ 1/2 n Wx 0i )dt (1 + t α 2 ) e −λ t 2 2 /4 dt.

B.4 Proof of Corollary 4.1
Proof. We first show the convergence of the mean and covariance of π * in (t | A), which is a direct consequence of Theorem 4.1: then by Theorem 3.2 and Slutsky's Theorem, √ nG Note that G 0in is finite and positive definite. By continuous mapping theorem, We now focus on the last assertion. By the previous proof, we know that max i∈ [n] x * − x i 2 2 = o P0 (1/n). It follows directly that Therefore, by Theorem 3.2 and Cauchy-Schwarz inequality, we have where ·, · F denotes the Frobenius inner product between matrices. The proof is thus completed.

C Proof of the Convergence of the Stochastic Gradient Descent
Lemma C.1 (Lemma A.5 in Mairal, 2013). Let (a t ) t≥1 , (b t ) t≥1 be two non-negative real sequences. Assume that ∞ t=1 a t b t converges and ∞ t=1 a t diverges, and |b t+1 − b t | ≤ Ka t for some constant K ≥ 0. Then b t converges to 0.
Lemma C.2 (Lemma 2 in Li and Orabona, 2019). Let a 0 > 0, a i ≥ 0, i = 1, . . . , T and β > 1. Then Li and Orabona, 2019). Let f : X ⊂ R d → R be twice continuously differentiable whose minimum is attained at x = x * and suppose there exists a constant L > 0, such that for all x, y ∈ X , Suppose g(x, z) is a function of a random vector z, such that E z g(x, z) = ∂f (x)/∂x. Let (z t ) t≥1 be a sequence of independent and identically distributed (i.i.d.) copies of z. Consider a sequence of iterates x (t) generated by where H t ∈ R d×d is a step-size matrix for the tth iteration. Then the sequence (x (t) ) t≥1 satisfies the following inequality: of Theorem 3.3. The proof is similar to Theorem 1 in Li and Orabona (2019), with some slight modifications.
In the setting here, the expectation is taken with respect to the randomness of the stochastic gradient descent conditioned on the adjacency matrix, that is, the data and the ASE are viewed as deterministic. Here, we suppress the subscript i ∈ [n] and use x (t) to denote the tth iterate in the optimization, and x the maximizer of the average surrogate log-likelihood function M in (x) := (1/n) in (x).
For the surrogate log-likelihood function, by the computation of the gradient and Hessian of M in in the proof of Theorem 3.2, they are bounded over {x i : x i 2 ≤ 1} when max j x j 2 < 1. So both M in (x) and its gradient are Lipschitz in {x ∈ R d : x 2 ≤ 1} by the mean value theorem. Let C 1 and C 2 be the Lipschitz constants for M in (x) and its gradient, respectively. In the context of Section 3.2, the random vector z t corresponds to the randomly generated indices (j s ) in a single iteration of the mini-batch SGD algorithm, and g(x (t) , z t ) takes the form It is clear that E zt g(x (t) , z t ) coincides with the gradient of M in (x (t) ). Also, for the stochastic gradient

Additional details of the algorithms
This subsection provides the detailed Metropolis-Hastings sampler for computing the joint posterior distribution π n (X | A) using the surrogate likelihood function. For each i ∈ [n], we use the normal random walk truncated in the unit ball as the proposal distribution, with the covariance matrix being the inverse of .

D.2 An additional simulation example
In this subsection, we provide an additional simulation example to explore the performance of the proposed estimation methods in a comparatively small sample regime with n = 30 to supplement the example in Section 5.1. The simulation setup is the same as that in Section 5.1 but we set the number of vertices to n = 30.
Here, we focus on the performance of different estimates using the sum of square errors as the evaluation metric. Besides the four estimates considered in Section 5.1 of the manuscript (the adjacency spectral embedding, the one-step estimate, the maximum surrogate likelihood estimate, and the Bayes estimate with the surrogate likelihood), we also consider the maximum likelihood estimate. Note that although the theory of the maximum likelihood estimation is still open, it is always possible to find a local maximizer of the likelihood function using any optimization toolkit. Here we use the R built-in optim function in practice.
We repeat the same numerical experiment for 1000 independent Monte Carlo replicates and visualize the boxplots of the sum of squares errors in Figure 3. We can see that in this small sample regime, the maximum surrogate likelihood estimate and the onestep estimate do not outperform the baseline adjacency spectral embedding and the maximum likelihood estimate, while the Bayes estimate has the least sum of squares errors. This observation shows the potential advantage of the Bayesian estimation procedure based on the Markov chain Monte Carlo sampling algorithm over the classical optimization-based estimation methods for finite-sample problems in practice.

D.3 Convergence diagnostics of the Metropolis-Hastings sampler
In this subsection, we provide some convergence diagnostics of Metropolis-Hastings sampler. Specifically, Trace plot of a parameter whose Metropolis-Hastings sample gives a p-value less than 0.05. Bottom right panel: ACF plot of a parameter whose Metropolis-Hastings sample gives a p-value less than 0.05.