SANVI: A Fast Spectral-Assisted Network Variational Inference Method with an Extended Surrogate Likelihood Function

TL;DR

SANVI introduces a fast, spectral-assisted variational inference method for network models, significantly reducing computational costs while maintaining accuracy, by using an extended surrogate likelihood function within the generalized random dot product graph framework.

Contribution

It proposes a novel extended surrogate likelihood function and a spectral-assisted variational inference algorithm that improves computational efficiency in Bayesian network analysis.

Findings

Comparable accuracy to MCMC in estimating latent positions

Significantly reduced computational costs

Asymptotic efficiency of the estimator and posterior distribution

Abstract

Bayesian inference has been broadly applied to statistical network analysis, but suffers from the expensive computational costs due to the nature of Markov chain Monte Carlo sampling algorithms. This paper proposes a novel and computationally efficient Spectral-Assisted Network Variational Inference (SANVI) method within the framework of the generalized random dot product graph. The key idea is a cleverly designed extended surrogate likelihood function that enjoys two convenient features. Firstly, it decouples the generalized inner product of latent positions in the random graph model. Secondly, it relaxes the complicated domain of the original likelihood function to the entire Euclidean space. Leveraging these features, we design a computationally efficient Gaussian variational inference algorithm via stochastic gradient descent. Furthermore, we show the asymptotic efficiency of the maximum extended surrogate likelihood estimator and the Bernstein-von Mises limit of the variational posterior distribution. Through extensive numerical studies, we demonstrate the usefulness of the proposed SANVI algorithm compared to the classical Markov chain Monte Carlo algorithm, including comparable estimation accuracy for the latent positions and less computational costs.

Tables

Table 1. Table 1: The sums of squared errors (and their standard errors over 100 repetitions, in parenthesis) of ASE, OSE, BE, and SANVI, respectively, in the example of stochastic block model.

Estimate	ASE	OSE	BE	SANVI
$n=1000$	8.352	19.611	7.658	7.739
	(0.367)	(11.116)	(0.358)	(0.359)
$n=3000$	7.880	41.459	7.535	7.564
	(0.210)	(28.124)	(0.199)	(0.203)
$n=5000$	7.834	39.936	7.588	7.593
	(0.145)	(29.470)	(0.142)	(0.141)
$n=7000$	7.913	39.742	7.727	7.721
	(0.129)	(38.149)	(0.130)	(0.130)
$n=9000$	7.833	23.114	7.690	7.685
	(0.115)	(18.278)	(0.116)	(0.116)
$n=10000$	7.774	19.486	7.650	7.643
	(0.094)	(17.361)	(0.092)	(0.093)

Table 2. Table 2: The p p -values of paired two-sample t t -tests of the sums of squared errors of ASE with BE and SANVI, respectively, in the example of the stochastic block model.

$n$	$n=1000$	$n=3000$	$n=5000$	$n=7000$
ASE vs BE	$2.7\times {10}^{-82}$	$7.1\times {10}^{-87}$	$5.5\times {10}^{-81}$	$3.7\times {10}^{-84}$
ASE vs SANVI	$2.2\times {10}^{-76}$	$4.9\times {10}^{-81}$	$6.9\times {10}^{-75}$	$8.1\times {10}^{-85}$
	$n=9000$	$n=10000$
ASE vs BE	$7.3\times {10}^{-84}$	$3.4\times {10}^{-78}$
ASE vs SANVI	$4.4\times {10}^{-76}$	$1.1\times {10}^{-76}$

Table 3. Table 3: The sums of squared errors (and their standard errors over 100 repetitions, in parenthesis) of ASE, OSE, BE, and SANVI, respectively, in the example of degree-corrected stochastic block model.

Estimate	ASE	OSE	BE	SANVI
$n=1000$	3.323	3.307	3.097	3.082
	(0.133)	(0.443)	(0.137)	(0.138)
$n=3000$	3.340	3.157	3.161	3.157
	(0.074)	(0.068)	(0.069)	(0.069)
$n=5000$	3.317	3.140	3.148	3.146
	(0.059)	(0.058)	(0.058)	(0.058)
$n=7000$	3.302	3.124	3.134	3.133
	(0.049)	(0.046)	(0.046)	(0.046)
$n=9000$	3.330	3.149	3.160	3.160
	(0.040)	(0.038)	(0.038)	(0.038)
$n=10000$	3.291	3.117	3.128	3.128
	(0.042)	(0.040)	(0.041)	(0.041)

Table 4. Table 4: The p p -values of paired two-sample t t -tests of the sums of squared errors of ASE with OSE, BE, and SANVI, respectively, in the example of the degree-corrected stochastic block model.

	$n=1000$	$n=3000$	$n=5000$	$n=7000$
ASE vs OSE	$0.7$	$3.9\times {10}^{-87}$	$1.4\times {10}^{-101}$	$1.1\times {10}^{-105}$
ASE vs BE	$1.1\times {10}^{-65}$	$1.2\times {10}^{-85}$	$3.8\times {10}^{-98}$	$1.8\times {10}^{-102}$
ASE vs SANVI	$9.8\times {10}^{-66}$	$2.1\times {10}^{-86}$	$9.6\times {10}^{-99}$	$2.6\times {10}^{-102}$
	$n=9000$	$n=10000$
ASE vs OSE	$4.2\times {10}^{-115}$	$8.1\times {10}^{-111}$
ASE vs BE	$1.3\times {10}^{-108}$	$1.7\times {10}^{-106}$
ASE vs SANVI	$3.2\times {10}^{-110}$	$3.5\times {10}^{-107}$

Table 5. Table 5: The sums of squared errors (and their standard errors over 100 repetitions, in parenthesis) of ASE, OSE, BE, and SANVI, respectively, in the example of rank-two latent curve.

Estimate	ASE	OSE	BE	SANVI
$n=1000$	31.800	41.128	28.713	28.655
	(0.381)	(11.628)	(0.490)	(0.546)
$n=3000$	64.216	86.080	60.203	60.119
	(0.371)	(19.760)	(0.521)	(0.556)
$n=5000$	85.716	114.552	80.924	80.711
	(16.277)	(31.286)	(16.267)	(16.077)
$n=7000$	28.950	76.529	25.113	25.420
	(0.995)	(51.900)	(0.848)	(0.872)
$n=9000$	25.198	73.932	21.978	22.260
	(0.644)	(52.166)	(0.559)	(0.572)
$n=10000$	24.153	71.077	21.112	21.383
	(0.517)	(51.425)	(0.441)	(0.452)

Table 6. Table 6: The p p -values of paired two-sample t t -tests of the sums of squared errors of ASE with BE and SANVI, respectively, in the example of the rank-two latent curve.

	$n=1000$	$n=3000$	$n=5000$	$n=7000$
ASE vs BE	$4.4\times {10}^{-91}$	$1.2\times {10}^{-99}$	$1.9\times {10}^{-105}$	$1.1\times {10}^{-127}$
ASE vs SANVI	$9.2\times {10}^{-83}$	$1.8\times {10}^{-95}$	$6.0\times {10}^{-101}$	$1.5\times {10}^{-125}$
	$n=9000$	$n=10000$
ASE vs BE	$1.7\times {10}^{-136}$	$2.2\times {10}^{-139}$
ASE vs SANVI	$8.0\times {10}^{-133}$	$9.6\times {10}^{-135}$

Table 7. Table 7: The sums of squared errors (and their standard errors over 100 repetitions, in parenthesis) of ASE, OSE, BE, and SANVI, respectively, in the example of rank-three latent curve.

Estimate	ASE	OSE	BE	SANVI
$n=1000$	74.147	87.363	69.920	70.058
	(2.581)	(12.521)	(2.559)	(2.589)
$n=3000$	58.578	107.439	55.517	55.835
	(8.160)	(37.719)	(7.940)	(7.954)
$n=5000$	46.232	107.456	44.070	44.260
	(1.039)	(38.946)	(0.993)	(0.996)
$n=7000$	42.373	112.582	40.593	40.747
	(0.590)	(51.598)	(0.558)	(0.565)
$n=9000$	40.291	123.629	38.719	38.860
	(0.448)	(63.700)	(0.413)	(0.424)
$n=10000$	39.572	132.974	38.079	38.205
	(0.369)	(75.372)	(0.354)	(0.358)

Table 8. Table 8: The p p -values of paired two-sample t t -tests of the sums of squared errors of ASE with BE and SANVI, respectively, in the example of rank-three latent curve.

	$n=1000$	$n=3000$	$n=5000$	$n=7000$
ASE vs BE	$3.0\times {10}^{-95}$	$1.6\times {10}^{-100}$	$2.1\times {10}^{-112}$	$1.7\times {10}^{-118}$
ASE vs SANVI	$5.0\times {10}^{-87}$	$1.5\times {10}^{-95}$	$1.1\times {10}^{-105}$	$2.3\times {10}^{-106}$
	$n=9000$	$n=10000$
ASE vs BE	$1.3\times {10}^{-115}$	$9.7\times {10}^{-124}$
ASE vs SANVI	$5.3\times {10}^{-107}$	$2.4\times {10}^{-114}$

Table 9. Table 9: The adjusted Rand indices computed from the four estimates and the corresponding computation time for Section 4.5 .

Estimate	ASE	OSE	BE	SANVI
Adjusted Rand Index	0.1321	0.0416	0.4374	0.3117
Computation Time (seconds)	0.05	0.05	45.20	16.64

Equations

$${\ell }_{\mathbf{A}}(\mathbf{X})=\sum _{1\le i\le j\le n}\{A_{ij}\log ({\mathbf{x}}_i^{\mathrm{T}}{\mathbf{I}}_{p,q}{\mathbf{x}}_j)+(1-A_{ij})\log ({\mathbf{x}}_i^{\mathrm{T}}{\mathbf{I}}_{p,q}{\mathbf{x}}_j)\}.$$

$${\ell }_{0in}({\mathbf{x}}_i)$$

$$+\{A_{ii}\log ({\mathbf{x}}_i^{\mathrm{T}}{\mathbf{I}}_{p,q}{\mathbf{x}}_i)+(1-A_{ii})\log (1-{\mathbf{x}}_i^{\mathrm{T}}{\mathbf{I}}_{p,q}{\mathbf{x}}_i)\}.$$

$${\ell }_{0in}({\mathbf{x}}_i)\approx \sum _{j=1}^n\{A_{ij}\log ({\mathbf{x}}_i^{\mathrm{T}}{\tilde{\mathbf{x}}}_j)+(1-A_{ij})\log (1-{\mathbf{x}}_i^{\mathrm{T}}{\tilde{\mathbf{x}}}_j)\}.$$

$${\ell }_{in}({\mathbf{x}}_i;\{{\mathbf{x}}_j{\}}_{j\ne i})=\sum _{j=1}^n\{A_{ij}\log ({\mathbf{x}}_i^{\mathrm{T}}{\mathbf{I}}_{p,q}{\mathbf{x}}_j)+(1-A_{ij})\log (1-{\mathbf{x}}_i^{\mathrm{T}}{\mathbf{I}}_{p,q}{\mathbf{x}}_j)\}.$$

$${\psi }_n(t)$$

$${\hat{\ell }}_{in}({\mathbf{x}}_i)=\sum _{j=1}^n\{A_{ij}{\psi }_n({\mathbf{x}}_i^{\mathrm{T}}{\tilde{\mathbf{x}}}_j)+(1-A_{ij}){\psi }_n(1-{\mathbf{x}}_i^{\mathrm{T}}{\tilde{\mathbf{x}}}_j)\}.$$

$${\pi }_{in}({\mathbf{x}}_i\mid \mathbf{A})\propto \exp \{{\hat{\ell }}_{in}({\mathbf{x}}_i)\}{\pi }_i({\mathbf{x}}_i).$$

$$\mathcal{Q}=\left\{\prod _{i=1}^n\mathcal{N}({\mathbf{x}}_i\mid {\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i):{\mathbf{\mu }}_i\in {\mathbb{R}}^d,{\mathbf{\Sigma }}_i\in {\mathbb{M}}_{+}(d),i\in [n]\right\},$$

$$\underset{{\mathbf{\mu }}_i\in {\mathbb{R}}^d,{\mathbf{\Sigma }}_i\in {\mathbb{M}}_{+}(d)}{\min }{D}_{\mathrm{KL}}\left({\varphi }_d({\mathbf{x}}_i\mid {\mathbf{\mu }}_i,n^{-1}{\mathbf{L}}_i{\mathbf{L}}_i^{\mathrm{T}})\Vert {\pi }_{in}({\mathbf{x}}_i\mid \mathbf{A})\right),i\in [n].$$

$$D({\varphi }_d({\mathbf{x}}_i\mid {\mathbf{\mu }}_i,n^{-1}{\mathbf{L}}_i{\mathbf{L}}_i^{\mathrm{T}})\Vert {\pi }_{in}({\mathbf{x}}_i\mid \mathbf{A}))$$

$$\displaystyle=-\log\det(\mathbf{L}_{i})-\frac{d}{2}\log(2\pi)-\mathbb{E}_{\mathbf{z}_{i}}\bigg{(}\frac{1}{2}\|\mathbf{z}_{i}\|_{2}^{2}\bigg{)}$$

$$-{\mathbb{E}}_{{\mathbf{z}}_i}\left\{{\hat{\ell }}_{in}\left({\mathbf{\mu }}_i+\frac{1}{\sqrt{n}}{\mathbf{L}}_i{\mathbf{z}}_i\right)+\log {\pi }_i\left({\mathbf{\mu }}_i+\frac{1}{\sqrt{n}}{\mathbf{L}}_i{\mathbf{z}}_i\right)\right\}+\log (d_{in}),$$

$$F_{in}({\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i)=-\log \det ({\mathbf{L}}_i)-{\mathbb{E}}_{{\mathbf{z}}_i}\left\{{\hat{\ell }}_{in}\left({\mathbf{\mu }}_i+\frac{1}{\sqrt{n}}{\mathbf{L}}_i{\mathbf{z}}_i\right)+\log {\pi }_i\left({\mathbf{\mu }}_i+\frac{1}{\sqrt{n}}{\mathbf{L}}_i{\mathbf{z}}_i\right)\right\},$$

$$\underset{{\mathbf{\mu }}_i\in {\mathbb{R}}^d,{\mathbf{L}}_i\in \mathcal{L}}{\min }{F}_{in}({\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i).$$

$$f_{in}({\mathbf{\mu }}_i,{\mathbf{L}}_i,{\mathbf{z}}_i)=-\log \det ({\mathbf{L}}_i)-{\hat{\ell }}_{in}\left({\mathbf{\mu }}_i+\frac{1}{\sqrt{n}}{\mathbf{L}}_i{\mathbf{z}}_i\right)-\log {\pi }_i\left({\mathbf{\mu }}_i+\frac{1}{\sqrt{n}}{\mathbf{L}}_i{\mathbf{z}}_i\right).$$

$$\frac{\partial f_{in}}{\partial {\mathbf{\mu }}_i}({\mathbf{\mu }}_i,{\mathbf{L}}_i,{\mathbf{z}}_i)$$

$$\frac{\partial f_{in}}{\partial {\mathbf{L}}_i}({\mathbf{\mu }}_i,{\mathbf{L}}_i,{\mathbf{z}}_i)$$

$${-\frac{1}{\sqrt{n}}\mathrm{tril}\left\{\frac{\partial }{\partial {\mathbf{x}}_i}\log {\pi }_i\left({\mathbf{x}}_i\right){\mathbf{z}}_i^{\mathrm{T}}\right\}|}_{{\mathbf{x}}_i={\mathbf{\mu }}_i+\frac{1}{\sqrt{n}}{\mathbf{L}}_i{\mathbf{z}}_i},$$

$${\tilde{h}}_n(x)=\{\begin{array}{rlr}& \frac{c_n}{c_{n}x+1}& \text{mbox}ifx>0,\\ & -c_n^{2}x+c_n& \text{mbox}ifx\le 0.\end{array}$$

$${\tilde{\nabla }}_{{\mathbf{L}}_i}{f}_{in}({\mathbf{\mu }}_i,{\mathbf{L}}_i,{\mathbf{z}}_i)$$

$$\tilde{\mathbf{X}}=[{\tilde{\mathbf{x}}}_1,\dots ,{\tilde{\mathbf{x}}}_n{]}^{\mathrm{T}}=\left[\mathrm{sign}({\hat{\lambda }}_1)|{\hat{\lambda }}_1{|}^{1/2}{\hat{\mathbf{u}}}_1,\dots ,\mathrm{sign}({\hat{\lambda }}_d)|{\hat{\lambda }}_d{|}^{1/2}{\hat{\mathbf{u}}}_d\right].$$

$$m_{{\mathbf{\mu }}_i,1}^{(t)}$$

$$m_{{\mathbf{L}}_i,1}^{(t)}$$

$$m_{{\mathbf{\mu }}_i,2}^{(t)}$$

$$m_{{\mathbf{L}}_i,2}^{(t)}$$

$${\mathbf{\mu }}_i^{(t)}$$

$$\mathbb{P}\left({\hat{\mathbf{x}}}_{i}existsandisuniqueforalli\in [n]\right)>1-n^{-c},$$

$$\mathbb{P}\left\{\underset{i\in [n]}{\max }{\Vert {\mathbf{W}}^{\mathrm{T}}{\hat{\mathbf{x}}}_i-{\rho }_n^{1/2}{\mathbf{x}}_{0i}\Vert }_2<C_{c,\delta ,\lambda }\sqrt{\frac{\log n}{n}}\right\}>1-n^{-c},$$

$${\int }_{{\mathbb{R}}^d}\left|{\pi }_{in}({\mathbf{x}}_i\mid \mathbf{A})-{\varphi }_d({\mathbf{x}}_i\mid {\hat{\mathbf{x}}}_i,(n\mathbf{W}{\mathbf{G}}_{0in}{\mathbf{W}}^{\mathrm{T}}{)}^{-1})\right|\mathrm{d}{\mathbf{x}}_i{\lesssim }_{c,\delta ,\lambda }\frac{1}{\log n}.$$

$$q_{in}^{\ast }({\mathbf{x}}_i)=\underset{q\in {\mathcal{Q}}_d}{\mathrm{argmin}}{D}_{\mathrm{KL}}(q({\mathbf{x}}_i)\Vert {\pi }_{in}({\mathbf{x}}_i|\mathbf{A}))$$

$${\int }_{{\mathbb{R}}^d}\left|q_{in}^{\ast }({\mathbf{x}}_i)-{\varphi }_d({\mathbf{x}}_i\mid {\hat{\mathbf{x}}}_i,(n\mathbf{W}{\mathbf{G}}_{0in}{\mathbf{W}}^{\mathrm{T}}{)}^{-1})\right|\mathrm{d}{\mathbf{x}}_i{\lesssim }_{c,\delta ,\lambda }\sqrt{\frac{1}{\log n}}$$

$$\tilde{\mathbf{X}}\mathbf{W}-{\rho }_n^{1/2}{\mathbf{X}}_0{\mathbf{I}}_{p,q}={\rho }_n^{-1/2}(\mathbf{A}-\mathbf{P}){\mathbf{X}}_0({\mathbf{X}}_0^{\mathrm{T}}{\mathbf{X}}_0{)}^{-1}+{\mathbf{R}}_{\tilde{\mathbf{X}}},$$

$$\Vert \tilde{\mathbf{X}}\mathbf{W}-{\rho }_n^{1/2}{\mathbf{X}}_0{\mathbf{I}}_{p,q}{\Vert }_{2\to \infty }$$

$${\mathbf{\Delta }}_n=\left[\begin{array}{cc}\frac{1}{n}{\mathbf{X}}_{0+}^{\mathrm{T}}{\mathbf{X}}_{0+}& \\ & \frac{1}{n}{\mathbf{X}}_{0-}^{\mathrm{T}}{\mathbf{X}}_{0-}\end{array}\right].$$

Full text

**Keywords: **generalized random dot product graphs, extended surrogate likelihood, variational Bayes, stochastic gradient descent

1 Introduction

Using graphs, a mathematical abstraction of real-world networks, to represent relational data, with the vertices denoting entities and the edges encoding relationships between connected entities, has been attracting attention in a broad range of applications, such as social networks (Girvan and Newman,, 2002; Wasserman and Faust,, 1994; Young and Scheinerman,, 2007), biological networks (Girvan and Newman,, 2002; Tang et al.,, 2019), and computer networks (Neil et al.,, 2013; Rubin-Delanchy et al.,, 2016), among others. Network analysis also connects to other fields beyond statistics, including computer science, machine learning, probability, and physics. A variety of network models that are conformable to statistical analyses have been developed, including the renowned stochastic block model (Holland et al.,, 1983) as well as its offspring (Airoldi et al.,, 2008; Karrer and Newman,, 2011; Lyzinski et al.,, 2017), the (generalized) random dot product graph model (Rubin-Delanchy et al.,, 2022; Young and Scheinerman,, 2007), the latent space model (Hoff et al.,, 2002), exchangeable random graphs (Caron and Fox,, 2017; Lei,, 2021), and graphons (Lovász,, 2012). Meanwhile, there has also been substantial progress on the subsequent inference tasks for the latent structures of network models, such as community detection (Abbe,, 2018; Abbe et al.,, 2016; Lei and Rinaldo,, 2015; Sussman et al.,, 2012), vertex classification (Sussman et al.,, 2014; Tang et al.,, 2013), and network hypothesis testing (Lei,, 2016; Tang et al., 2017a, ; Tang et al., 2017b, ).

In this paper, we focus on the generalized random dot product graph (GRDPG). Informally, GRDPG assigns each vertex a low-dimensional vector called the latent position, and the connection probability between any pair of vertices is given by the generalized inner product of the associated latent positions. We defer the formal definition to Section 2.1. GRDPG has been attracting attention because it not only has a simple low-rank structure but also is versatile as it encompasses several popular network models, such as stochastic block models (Holland et al.,, 1983), degree-corrected stochastic block models (Karrer and Newman,, 2011), mixed membership stochastic block models (Airoldi et al.,, 2008), and degree-corrected mixed membership models (Jin et al.,, 2023). GRDPG also provides building blocks for approximating general latent position random graphs (Lei,, 2021; Tang et al.,, 2013).

Graph data is usually represented in the form of an adjacency matrix. Due to the low expected rank of the adjacency matrix generated from a GRDPG, spectral methods have been widely applied in statistical analysis of graph data, among which the adjacency spectral embedding (ASE) is a popular one. The random dot product graph (RDPG) community has been developing theory and methods based on ASE. The readers are referred to Athreya et al., (2016, 2021); Koo et al., (2023); Levin and Levina, (2025); Levin et al., (2021); Li et al., (2020); Lyzinski et al., (2014); Rubin-Delanchy et al., (2022); Sengupta and Chen, (2017); Sussman et al., (2014, 2012); Tang et al., 2017a ; Tang et al., 2017b ; Tang and Priebe, (2018); Xie, (2023, 2024); Xie and Wu, (2023); Xie and Xu, (2020, 2023); Young and Scheinerman, (2007) for an incomplete list of references. However, it has also been observed in Xie and Xu, (2020) that spectral estimators do not take advantage of the likelihood information of the network adjacency matrix, and likelihood-based methods for (generalized) RPDG are comparatively unexplored. This research theme aims to develop a novel likelihood-based method for learning GRDPG that is computationally efficient, numerically stable for finite-sample problems, and theoretically solid and optimal.

Recently, Xie and Xu, (2023) discovered a striking fact: Spectral estimators are sub-optimal for estimating the latent positions due to the negligence of the graph likelihood structure. Specifically, Xie and Xu, (2023) proposed a one-step estimator (OSE) that absorbs the network likelihood information and established that OSE improves upon ASE. Better estimation of the latent positions is not only interesting by itself but also useful for more effective subsequent inference methods, such as more powerful hypothesis testing of the equality of latent positions (Xie,, 2024) or membership profiles in mixed membership models (Fan et al.,, 2022).

Despite the large-sample optimality, OSE typically requires comparatively large network sizes to outperform ASE (Xie and Xu,, 2023). For small-network problems, OSE can be numerically unstable because the estimated Hessian matrix may contain negative eigenvalues. Subsequently, Wu and Xie, (2025) developed a Bayesian method for RDPG based on a cleverly-designed surrogate likelihood that retains more likelihood information than OSE does. The Bayes estimate based on the surrogate likelihood is not only asymptotically efficient but also exhibits superior numerical stability compared to OSE and ASE, even for moderately small network sizes.

Nevertheless, the Bayesian methods, although theoretically solid and numerically competitive, are practically inconvenient. This is largely due to the expensive computational cost associated with Markov chain Monte Carlo (MCMC) sampling methods. Compared to classical MCMC methods, variational inference (VI) methods (Blei et al.,, 2017) have emerged as a popular alternative. Unlike MCMC, VI is optimization-based, which tends to be faster while still having comparable numerical performance. VI methods have been gaining rapid development recently. Readers are referred to Bhattacharya et al., (2025); Katsevich and Rigollet, (2024); Zhang and Yang, (2024); Wang and Blei, (2019); Han and Yang, (2019); Hinton and van Camp, (1993); Jordan et al., (1998); Peterson and Anderson, (1987) and references therein for the recent advances of VI methods in general. For VI methods in the context of network models, Loyal, (2024) and Zhao et al., (2024) developed structured mean-field VI methods for dynamic networks that were built upon latent space models, and they do not apply to the GRDPG framework.

In this paper, we propose a computationally efficient spectral-assisted network variational inference (SANVI) method through a pivotal extended surrogate likelihood (ESL) function in the context of GRDPG. Given a fixed vertex, SANVI minimizes the Kullback–Leibler (KL) divergence between a candidate Gaussian distribution and the posterior distribution of the latent position of interest, where the posterior distribution is computed based on the ESL function. Note that the algorithm can be parallelized thanks to the separable structure of the ESL function. We also establish the corresponding large sample theory, including the asymptotic efficiency of the proposed estimator and the Bernstein-von Mises theorem of the variational posterior distribution, thereby generalizing and popularizing the existing framework in Wu and Xie, (2025).

The remaining part of this paper is structured as follows. In Section 2, we review the background of GRDPG and ASE and introduce the ESL function. In Section 3, we introduce SANVI based on the Gaussian VI by leveraging the ESL function, establish the asymptotic properties of the variational posterior distribution, and discuss the stochastic gradient descent algorithm for the computation of SANVI. In Section 4, we demonstrate the empirical finite-sample performance of SANVI through some simulated examples and the analysis of a real-world network dataset. We conclude the paper with a discussion in Section 5.

Notations: Most of the notations that we mainly use in this paper are explained in the following. The notation $[n]$ stands for the set of consecutive integers from $1$ to $n$, that is, $[n]=\{1,\ldots,n\}$. The symbol $\lesssim$ means an inequality up to a constant, that is, $a\lesssim b$ if $a\leq Cb$ for some constant $C>0$. The constant $C$ can depend on some other constants, of which we use subscripts to denote the dependency, e.g., $C_{\delta,\lambda}$ showing the dependency of $C$ on $\delta$ and $\lambda$. A similar definition also applies to the symbol $\gtrsim$. The notation $\|\mathbf{x}\|_{2}$ denotes the Euclidean norm of a vector $\mathbf{x}=[x_{1},\ldots,x_{d}]^{\mathrm{T}}\in\mathbb{R}^{d}$, that is, $\|\mathbf{x}\|_{2}=(\sum_{k=1}^{d}x_{k}^{2})^{1/2}$. The $d\times d$ identity matrix is denoted by ${\mathbf{I}}_{d}$. The notation $\mathbb{O}(n,d)=\{\mathbf{U}\in\mathbb{R}^{n\times d}:\mathbf{U}^{\mathrm{T}}\mathbf{U}=\mathbf{I}_{d}\}$ denotes the set of all orthonormal $d$-frames in $\mathbb{R}^{n}$, where $d\leq n$, and we write $\mathbb{O}(d)=\mathbb{O}(d,d)$. For a matrix $\mathbf{X}=[x_{ik}]_{n\times d}$, $\sigma_{k}(\mathbf{X})$ denotes its $k$th largest singular value. Matrix norms with following definitions are used: the spectral norm $\|\mathbf{X}\|_{2}=\sigma_{1}(\mathbf{X})$, the Frobenius norm $\|\mathbf{X}\|_{\mathrm{F}}=(\sum_{i=1}^{n}\sum_{k=1}^{d}x_{ik}^{2})^{1/2}$, the matrix infinity norm $\|\mathbf{X}\|_{\infty}=\max_{i\in[n]}\sum_{k=1}^{d}|x_{ik}|$, and the two-to-infinity norm $\|\mathbf{X}\|_{2\to\infty}=\max_{i\in[n]}(\sum_{k=1}^{d}x_{ik}^{2})^{1/2}$. In particular, these norm notations apply to any Euclidean vector $\mathbf{x}\in\mathbb{R}^{d}$ viewed as a $d\times 1$ matrix. Given two symmetric positive semidefinite matrices $\mathbf{A},\mathbf{B}$ of the same dimension, we write $\mathbf{A}\preceq\mathbf{B}$ ($\mathbf{A}\succeq\mathbf{B}$, respectively) if $\mathbf{B}-\mathbf{A}$ ($\mathbf{A}-\mathbf{B}$, respectively) is positive semidefinite. For a vector $\mathbf{x}\in\mathbb{R}^{d}$, the notation $[\mathbf{x}]_{k}=x_{k}$ denotes its $k$th coordinate. For a matrix $\mathbf{X}\in\mathbb{R}^{n\times d}$, the notation $\mathbf{X}_{i*}$ denotes its $i$th row, $\mathbf{X}_{*j}$ its $j$th column, and $x_{ij}$ its $(i,j)$th entry. We use $\{\mathbf{W}_{n}\}_{n=1}^{\infty}$ to denote the sequence of orthogonal matrices aligning a sequence of estimators $\{\widehat{\mathbf{X}}_{n}\}_{n=1}^{\infty}$ and the true value, and we may drop the subscript $n$ for simplicity of notation.

2 Generalized random dot product graphs and the extended surrogate likelihood

2.1 Background on generalized random dot product graphs

We begin by briefly reviewing GRDPG and ASE.

Definition 2.1* (Generalized random dot product graph).*

Let $n,d\in\mathbb{N}_{+}$, $n\geq d$, $p,q\in\mathbb{N}$ with $p+q=d$, and ${\mathbf{I}}_{p,q}=\mathrm{diag}(1,\ldots,1,-1,\ldots,-1)$ with $p$ positive ones followed by $q$ negative ones on its diagonal. Given an $n\times d$ matrix $\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]^{\mathrm{T}}$, where $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\in\mathbb{R}^{d}$, with first $p$ columns orthogonal to last $q$ columns, such that $\mathbf{x}_{i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{j}\in[0,1]$ for all $i,j\in[n]=\{1,\ldots,n\}$, we say that $\mathbf{A}=[A_{ij}]_{n\times n}$ is the adjacency matrix of a generalized random dot product graph, denoted as $\mathbf{A}\sim\mathrm{GRDPG}(\mathbf{X})$ with signature $(p,q)$ if $A_{ij}\sim\mathrm{Bernoulli}(\mathbf{x}_{i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{j})$ independently for all $i\leq j$, and $A_{ij}=A_{ji}$ if $i>j$. The matrix $\mathbf{X}$ is referred to as the latent position matrix, and the $d$-dimensional vector $\mathbf{x}_{i}$ is referred to as the latent position of vertex $i$. When $q=0$, a GRDPG is also called a random dot product graph (RDPG).

In this paper, we consider the latent positions $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ to be deterministic parameters to be estimated. Another slightly different modeling approach is to consider $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ as independent and identically distributed latent random variables following some distribution $F$ supported on the latent space $\mathcal{X}$ (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, (2023) using a Glivenko–Cantelli type condition when $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ are deterministic. The latter Glivenko–Cantelli type condition is also relaxed in the current work, as we only require that $\sigma_{d}(\mathbf{X})>0$ (see Remark 1 below).

Remark 1* (Nonidentifiability).*

For convenience, in this work, we follow the setup in Xie, (2024) and require the first $p$ columns of the latent position matrix $\mathbf{X}$ to be orthogonal to the last $q$ columns. For GRDPG with more general latent position matrices, please see Rubin-Delanchy et al., (2022). The latent position matrix $\mathbf{X}$ is not uniquely identified in the following two senses. First, any low-rank connection probability matrix $\mathbf{P}=\mathbf{X}{\mathbf{I}}_{p,q}\mathbf{X}^{\mathrm{T}}$ can have different factorizations because for any orthogonal matrix with a $(p,q)$ block structure, $\mathbf{W}=\mathrm{diag}(\mathbf{W}_{p},\mathbf{W}_{q})$, where $\mathbf{W}\in\mathbb{O}(d)$, $\mathbf{W}_{p}\in\mathbb{O}(p)$, $\mathbf{W}_{q}\in\mathbb{O}(q)$, we have $\mathbf{X}{\mathbf{I}}_{p,q}\mathbf{X}^{\mathrm{T}}=(\mathbf{X}\mathbf{W}){\mathbf{I}}_{p,q}(\mathbf{X}\mathbf{W})^{\mathrm{T}}$. Second, for any $d^{\prime}>d$ and any latent position matrix $\mathbf{X}\in\mathbb{R}^{n\times d}$, there exists another matrix $\mathbf{X}^{\prime}\in\mathbb{R}^{n\times d^{\prime}}$ such that $\mathbf{X}{\mathbf{I}}_{p,q}\mathbf{X}^{\mathrm{T}}=\mathbf{X}^{\prime}{\mathbf{I}}_{p+(d^{\prime}-d),q}(\mathbf{X}^{\prime})^{\mathrm{T}}$. The latter source of non-identifiability can be removed by requiring that $\sigma_{d}(\mathbf{X})>0$, while the former source is inevitable without further constraints. Thus, any estimator of the latent position matrix $\mathbf{X}$ can only recover it up to an orthogonal transformation.

We consider undirected and unweighted graphs, so the adjacency matrices are binary and symmetric. We allow self-loops, so the adjacency matrices may have non-zero diagonal elements. One convenient feature of GRDPG is that the edge probability matrix $\mathbf{P}=\mathbb{E}\mathbf{A}=\mathbf{P}=\mathbf{X}{\mathbf{I}}_{p,q}\mathbf{X}^{\mathrm{T}}$ is low-rank. This motivates spectral decomposition methods for learning the latent position matrix $\mathbf{X}$ (Rubin-Delanchy et al.,, 2022; Sussman et al.,, 2012).

Definition 2.2* (Adjacency spectral embedding).*

Given $\mathbf{A}\sim\mathrm{GRDPG}(\mathbf{X})$, let $\mathbf{A}$ yield spectral decomposition $\mathbf{A}=\sum_{i=1}^{n}\lambda_{i}(\mathbf{A})\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\mathrm{T}}$, where $|\lambda_{1}(\mathbf{A})|\geq\ldots\geq|\lambda_{n}(\mathbf{A})|$, arranged in decreasing order of absolute value, are the eigenvalues of $\mathbf{A}$, $\widehat{\mathbf{u}}_{i}$ is the eigenvector associated with $\lambda_{i}(\mathbf{A})$, $\widehat{\mathbf{u}}_{i}^{\mathrm{T}}\widehat{\mathbf{u}}_{j}=0$ for all $i\neq j$, and $\|\mathbf{u}_{i}\|_{2}=1$ for all $i\in[n]$. Pick the first $d$ eigenvalues and the correspongding eigenvectors, and rearrange them in the decreasing order of the eigenvalues as real numbers, $\lambda_{k_{1}}(\mathbf{A})\geq\ldots\geq\lambda_{k_{d}}(\mathbf{A})$. Then the adjacency spectral embedding of $\mathbf{A}$ into $\mathbb{R}^{n\times d}$ is defined as $\breve{\mathbf{X}}=[\breve{\mathbf{x}}_{1},\ldots,\breve{\mathbf{x}}_{n}]^{\mathrm{T}}=\mathbf{U}_{\mathbf{A}}|\mathbf{S}_{\mathbf{A}}|^{1/2}$, where $\mathbf{U}_{\mathbf{A}}=[\widehat{\mathbf{u}}_{k_{1}},\ldots,\widehat{\mathbf{u}}_{k_{d}}]$, $\mathbf{S}_{\mathbf{A}}=\mathrm{diag}\{\lambda_{k_{1}}(\mathbf{A}),\ldots,\lambda_{k_{d}}(\mathbf{A})\}$, and $|\mathbf{S}_{\mathbf{A}}|=\mathrm{diag}\{|\lambda_{k_{1}}(\mathbf{A})|,\ldots,|\lambda_{k_{d}}(\mathbf{A})|\}$. Also, the signature-adjusted adjacency spectral embedding of $\mathbf{A}$ into $\mathbb{R}^{n\times d}$ is defined as $\widetilde{\mathbf{X}}=[\widetilde{\mathbf{x}}_{1},\ldots,\widetilde{\mathbf{x}}_{n}]^{\mathrm{T}}=\mathbf{U}_{\mathbf{A}}|\mathbf{S}_{\mathbf{A}}|^{1/2}\mathrm{sgn}(\mathbf{S}_{\mathbf{A}})$, where $\mathrm{sgn}(\cdot)$ is the sign function and $\mathrm{sgn}(\mathbf{S}_{\mathbf{A}})$ applies entrywise on the diagonals of $\mathbf{S}_{\mathbf{A}}$.

2.2 The extended surrogate likelihood function

We now introduce the extended surrogate likelihood function for GRDPG. The motivation is that the exact likelihood function has a complicated structure, bringing challenges for developing the theory and computation of maximum likelihood estimation. The difficulty partially comes from the fact that GRDPG belongs to a curved exponential family, and the theory of the maximum likelihood estimation is much more difficult in curved exponential families than in canonical ones (see, for example, Section 2.3 in Bickel and Doksum,, 2015).

Consider the log-likelihood function of $\mathbf{A}\sim\mathrm{GRDPG}(\mathbf{X})$:

[TABLE]

The parameter space is defined by $\{\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]^{\mathrm{T}}\in\mathbb{R}^{n\times d}:0<\mathbf{x}_{i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{j}<1\text{ for all }i,j\}$, which is a complicated set whose boundary renders the maximum likelihood estimation intractable, both computationally and analytically. In addition, the log-likelihood function has an unbounded gradient over the boundary.

For the sake of generality, we introduce an $n$-dependent sparsity factor $\rho_{n}\in(0,1]$ that governs the average expected degree of GRDPG through the quantity $n\rho_{n}$. Note that by taking $\rho_{n}\to 0$ as $n\to\infty$, we allow the modeling of sparse random graphs that are more practical for real-world network data. To distinguish a generic latent position $\mathbf{x}_{i}\in\mathbb{R}^{d}$ and its true value associated with the data generating distribution, let $\rho_{n}^{1/2}\mathbf{x}_{0i}$ denote the ground truth of $\mathbf{x}_{i}$, $i\in[n]$, and $\mathbf{X}_{0}=[\mathbf{x}_{01},\ldots,\mathbf{x}_{0n}]^{\mathrm{T}}$. We first consider the log-likelihood function of a single $\mathbf{x}_{i}$ when the remaining latent positions $\{\mathbf{x}_{0j}\}_{j\neq i}$ are accessible:

[TABLE]

We refer to $\ell_{0in}(\mathbf{x}_{i})$ in (2.1) as the oracle log-likelihood function. Theorem 2 in Xie and Xu, (2023) established the consistency and asymptotic normality of the maximizer of the oracle log-likelihood function $\ell_{0in}(\mathbf{x}_{i})$ in (2.1). Nevertheless, the oracle log-likelihood is not computable because $\{\mathbf{x}_{0j}\}_{j\neq i}$ are not accessible in practice. Following the idea in Wu and Xie, (2025), we replace the unknown latent positions together with the signature by the corresponding rows of the signature-adjusted adjacency spectral embedding. Formally, let $\widetilde{\mathbf{x}}_{j}$ be the $j$th row of the signature-adjusted adjacency spectral embedding $\widetilde{\mathbf{X}}$, $j\in[n]$. Then, we obtain the following approximation to the oracle log-likelihood:

[TABLE]

Note that the last term in $\ell_{0in}$ is replaced by $A_{ii}\log(\mathbf{x}_{i}^{\mathrm{T}}\widetilde{\mathbf{x}}_{i})+(1-A_{ii})\log(1-\mathbf{x}_{i}^{\mathrm{T}}\widetilde{\mathbf{x}}_{i})$ for convenience, which is asymptotically unimportant. The above approximation can be made precise by the uniform consistency of the adjacency spectral embedding: There exists a $d\times d$ orthogonal $\mathbf{W}$ such that $\|\widetilde{\mathbf{X}}\mathbf{W}-\rho_{n}^{1/2}\mathbf{X}_{0}\|_{2\to\infty}=O\{\sqrt{(\log n)/n}\}$ with high probability (Lyzinski et al.,, 2014; Xie,, 2024).

The complication of the oracle log-likelihood function primarily comes from the constraint that $\mathbf{x}_{i}^{\mathrm{T}}\mathbf{x}_{j}\in(0,1)$ for all $i,j\in[n]$. Nonetheless, the approximation step (2.2) does not fully resolve it. To address this technical challenge, Wu and Xie, (2025) proposed a surrogate likelihood method for learning the latent position matrix $\mathbf{X}$ for RDPG by applying a Taylor’s expansion of the term $\log(\mathbf{x}_{i}^{\mathrm{T}}\widetilde{\mathbf{x}}_{j})$ (note that in RDPG, ${\mathbf{I}}_{p,q}={\mathbf{I}}_{p}$), such that the constraint $\mathbf{x}_{i}^{\mathrm{T}}\mathbf{x}_{j}\in(0,1)$ is relaxed to $\|\mathbf{x}_{i}\|_{2}\leq 1$ due to the uniform consistency of ASE. The surrogate log-likelihood function for the entire graph is formed by taking the sum of the individual surrogate log-likelihood functions for each vertex. In particular, by doing so, the surrogate log-likelihood has a separable structure for each individual latent position $\mathbf{x}_{i}$ and provides immediate convenience for both theoretical analysis and practical computation.

Wu and Xie, (2025) observed that Bayesian methods are typically comparable and sometimes outperform the frequentist point estimator for RPDG, such as ASE and OSE. Additionally, Bayesian methods offer a natural and principled approach for uncertainty quantification. One disadvantage of the Bayesian method proposed in Wu and Xie, (2025) is the computational expense due to the nature of MCMC. This practical inconvenience motivates us to develop a computationally efficient VI method for GRDPG and to provide it with the necessary theoretical guarantee.

The surrogate likelihood proposed Wu and Xie, (2025) is only defined on a compact subset $\prod_{i=1}^{n}\{\mathbf{x}_{i}\in\mathbb{R}^{d}:\|\mathbf{x}_{i}\|_{2}\leq 1\}$ of the Euclidean space. In this work, we adopt the Gaussian VI and take the variational distribution family to be the space of all Gaussian distributions. The detailed formulation is deferred to Section 3, but one requirement of Gaussian VI is that the target posterior distribution needs to be supported over $\mathbb{R}^{d}$. Hence, extending the feasible set of the surrogate likelihood to the entire Euclidean space is necessary. In addition, the surrogate likelihood derivation relies on Taylor’s expansion argument, which is necessary to drop the constraint that $\mathbf{x}_{i}^{\mathrm{T}}\widetilde{\mathbf{x}}_{j}>0$ for all $j\in[n]$. However, this approximation step could still be rough in moderate and small network problems. The extended surrogate likelihood to be developed will address the above issues while preserving the attractive features of the surrogate likelihood, including separability and log-concavity.

For $\mathbf{A}\sim\mathrm{GRDPG}(\mathbf{X})$, the local log-likelihood for a single latent position $\mathbf{x}_{i}$ is

[TABLE]

Rather than using the $j$th row of ASE $\breve{\mathbf{x}}_{j}$ directly, we use the signature-adjusted adjacency spectral embedding, $\widetilde{\mathbf{x}}_{j}=\mathrm{sgn}(\mathbf{S}_{\mathbf{A}})\breve{\mathbf{x}}_{j}$, to replace ${\mathbf{I}}_{p,q}\mathbf{x}_{j}$, since $(p,q)$ may be unknown in practice. Next, observe that the functions $\log(t)$ and $\log(1-t)$ are both well defined over $[\tau,1-\tau]$ for a small threshold $\tau>0$ , but not the entire $\mathbb{R}$. It is desirable that these functions can be extended beyond this interval while certain regularities, such as differentiability and smoothness, are preserved. For this purpose, let $\tau_{n}$ be a small positive number that may depend on $n$, and we define $\psi_{n}(t)$ with the following properties: $\psi_{n}(t)=\log(t)$ on $[\tau_{n},1]$; $\psi_{n}$ is a quadratic function for $t<\tau_{n}$ and $t>1$; $\psi_{n}$ is twice continuously differentiable over $\mathbb{R}$. Formally,

[TABLE]

With these modifications of $\log(t)$ and $\log(1-t)$, we then define the local extended surrogate log-likelihood (ESL) function for the latent position of a single vertex $\mathbf{x}_{i}$ for GRDPG as

[TABLE]

The global ESL function for the entire graph is defined as $\widehat{\ell}_{n}(\mathbf{X})=\sum_{i=1}^{n}\widehat{\ell}_{in}(\mathbf{x}_{i})$. Here, the term “extended” means that the domain of the target function (2.3) is extended to the entire $\mathbb{R}^{d}$ without constraint, and the term “surrogate” means that the unknown latent positions are replaced by their signature-adjusted ASE.

3 Spectral-assisted network variational inference

We now leverage the ESL function to develop SANVI. To begin with, we first consider the posterior distribution of the latent positions $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ associated with the ESL function (2.3). Note that the spectral assistance occurs directly in the formulation of (2.3) since ($\widetilde{\mathbf{x}}_{1},\ldots,\widetilde{\mathbf{x}}_{n}$) are spectral estimators.

It should be noted that we do not use the exact likelihood for the entire graph, as it is difficult to analyze and not separable in $i\in[n]$. Additionally, we do not use the oracle likelihood, since the true latent positions are unknown. Instead, we substitute the signature-adjusted ASE for the unknown latent positions together with the signature in the oracle likelihood (hence the term “surrogate”). This idea of using a general statistical criterion function to replace the likelihood in the Bayes formula when the exact likelihood function is not available or intractable for analysis or computation is not new, and among the literature, an influential work is Chernozhukov and Hong, (2003).

Formally, given independent prior distributions with densities $\pi_{i}(\mathbf{x}_{i})$ over $\mathbb{R}^{d}$, $i\in[n]$ for the latent positions $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$, the posterior distribution of the latent position $\mathbf{x}_{i}$ of a single vertex $i$ associated with the ESL function (2.3) has the following density function up to a normalizing constant:

[TABLE]

The joint posterior density of the latent position matrix $\mathbf{X}$ of the entire graph takes the product form $\pi_{n}(\mathbf{X}\mid\mathbf{A})=\prod_{i=1}^{n}\pi_{in}(\mathbf{x}_{i}\mid\mathbf{A})$ thanks to the separable structure of $\widehat{\ell}_{n}(\mathbf{X})$. The exact computation of the posterior distribution in (3.1) typically relies on MCMC and is generally inconvenient, even though the separable structure permits parallelization. Instead, we resort to VI methods and focus on the Gaussian VI.

The goal of VI is to find a distribution $q(\mathbf{X})\in\mathcal{Q}$ for the latent position matrix $\mathbf{X}$, where $\mathcal{Q}$ is a collection of candidate distributions over $\mathbf{X}$ that are tractable to compute, such that the Kullback-Leibler (KL) divergence between $q(\mathbf{X})$ and the posterior distribution $\pi_{n}(\mathbf{X}\mid\mathbf{A})$ is minimized. Formally, VI solves $\min_{q\in\mathcal{Q}}D_{\mathrm{KL}}(q(\cdot)\|\pi_{n}(\cdot\mid\mathbf{A}))$, where $D_{\mathrm{KL}}$ denotes the KL divergence. Since $\pi_{n}(\mathbf{X}\mid\mathbf{A})$ factorizes as $\prod_{i}\pi_{in}(\mathbf{x}_{i}\mid\mathbf{A})$, it is also reasonable to require that $\mathcal{Q}$ reduces to the class of all product distributions of $\mathbf{x}_{i}$’s: $\mathcal{Q}=\{\prod_{i=1}^{n}q_{i}(\mathbf{x}_{i}):q_{i}\in\mathcal{P}\}$, where $\mathcal{P}$ is some distribution class for $\mathbf{x}_{i}$.

Specialized to the Gaussian VI, we take $\mathcal{P}$ as the class of all $d$-dimensional (non-degenerate) multivariate Gaussian distributions, namely,

[TABLE]

where $\mathcal{N}(\mathbf{x}_{i}\mid{\bm{\mu}}_{i},{\bm{\Sigma}}_{i})$ denotes the multivariate Gaussian distribution of $\mathbf{x}_{i}$ with mean ${\bm{\mu}}_{i}$ and covariance matrix ${\bm{\Sigma}}_{i}$, and $\mathbb{M}_{+}(d)$ denotes the class of all $d\times d$ symmetric positive definite matrices. Notationally, we use $\phi_{d}(\mathbf{x}_{i}\mid{\bm{\mu}}_{i},{\bm{\Sigma}}_{i})$ to denote the density function of the Gaussian distribution $\mathcal{N}(\mathbf{x}_{i}\mid{\bm{\mu}}_{i},{\bm{\Sigma}}_{i})$. It then follows that the Gaussian VI for the posterior distribution of vertex $i$ associated with the ESL function solves the following optimization problem:

[TABLE]

We call the Gaussian distribution with parameters being the solution to the above optimization problem the variational posterior distribution, and we call its mean parameter the variational inference estimator. We next introduce the computation and the theory of SANVI.

3.1 Computation algorithm

Following the idea in Xu and Campbell, (2023) and Kucukelbir et al., (2017), we reparameterize the covariance matrix ${\bm{\Sigma}}_{i}$ of $\mathcal{N}(\mathbf{x}_{i}\mid{\bm{\mu}}_{i},{\bm{\Sigma}}_{i})$ using the Cholesky factorization ${\bm{\Sigma}}_{i}=(1/n)\mathbf{L}_{i}\mathbf{L}_{i}^{\mathrm{T}}$, where $\mathbf{L}_{i}$ is a lower triangular matrix with positive diagonal entries. With the change of variable $\mathbf{x}_{i}={\bm{\mu}}_{i}+n^{-1/2}\mathbf{L}_{i}\mathbf{z}_{i}$ where $\mathbf{z}_{i}\sim\mathrm{N}_{d}({\bm{0}}_{d},{\mathbf{I}}_{d})$, a simple algebra shows that the objective function of VI is

[TABLE]

where $d_{in}=\int_{\mathbb{R}^{d}}\exp\{\widehat{\ell}_{in}(\mathbf{x}_{i})\}\pi_{i}(\mathbf{x}_{i})\mathrm{d}\mathbf{x}_{i}$ is the marginal density of the data matrix $\mathbf{A}$. Dropping the terms that do not depend on ${\bm{\mu}}_{i}$ and $\mathbf{L}_{i}$, we define the Gaussian VI objective function

[TABLE]

where $\widehat{\ell}_{in}(\mathbf{x}_{i})$ is the ESL function defined in (2.3), and $\mathcal{L}$ is the class of all $d\times d$ lower-triangular matrices with positive diagonals. Then the optimization problem (3.2) is equivalent to

[TABLE]

Denote the noisy version of the objective function in (3.3) by

[TABLE]

Then, a simple algebra shows that

[TABLE]

where $\mathrm{tril}(\mathbf{B})$ replaces the upper triangular entries (excluding diagonals) of a $d\times d$ matrix $\mathbf{B}$ with zeros.

Below, Theorem 3.1 establishes the strong convexity of $F_{in}({\bm{\mu}}_{i},\mathbf{L}_{i})$.

Theorem 3.1**.**

Suppose Assumption 1 and Assumption 2 hold. Then $F_{in}({\bm{\mu}}_{i},\mathbf{L}_{i})$ viewed as a function from $\mathbb{R}^{d}\times\mathcal{L}_{d\times d}$ to $\mathbb{R}$ is strongly convex with probability at least $1-n^{-c}$ for all $n\geq N_{c,\delta,\lambda}$ depending on $c,\delta,\lambda$, and $\nabla_{{\bm{\mu}}_{i},\mathbf{L}_{i}}F_{in}({\bm{\mu}}_{i},\mathbf{L}_{i})=\mathbb{E}_{\mathbf{z}}\left[\nabla_{{\bm{\mu}}_{i},\mathbf{L}_{i}}f_{in}({\bm{\mu}}_{i},\mathbf{L}_{i},\mathbf{z})\right]$.

The derivative of $\log\det\mathbf{L}_{i}$ with respect to $\mathbf{L}_{i}$ is $\mathrm{diag}(L_{11}^{-1},\ldots,L_{dd}^{-1})$, and $L_{kk}^{-1}$ is unbounded as $L_{kk}$ approaches 0 from the right. In practical implementation, to avoid unbounded gradients and improve the numerical stability of our gradient-based algorithm, we borrow the idea in Xu and Campbell, (2023) to modify the gradient of $\mathbf{L}_{i}$ as follows. Let $c_{n}$ be a positive number that depends on $n$, and consider the function

[TABLE]

The function $\widetilde{h}_{n}(x)$ has a continuous derivative at $x=0$ and asymptotically equals $\frac{1}{x}$ as $x$ goes to positive infinity. With this modification, the scaled gradient of $\log\det\mathbf{L}_{i}$ with respect to $\mathbf{L}_{i}$ is defined as the $d\times d$ diagonal matrix whose $k$th diagonal element is $\widetilde{h}_{n}(L_{kk})$. Then the scaled gradient of $f_{in}$ with respect to $\mathbf{L}_{i}$ is defined as the $d\times d$ matrix

[TABLE]

We adopt the Adam scheme in Kingma and Ba, (2015) to define the update step for stochastic gradient descent. See Algorithm 1 for the detailed SANVI computation algorithm.

3.2 Theoretical properties

We now introduce and establish the asymptotic properties of the variational posterior distribution whose parameters solve the optimization problem (3.3). Several assumptions are necessary before we state the main results.

Assumption 1*.*

The following conditions hold:

1. (a)

$d$, $p$, and $q$ are constant integers with $d\geq{}1$, $p\geq{}1$, $q\geq{}0$, and $d=p+q$. 2. (b)

$\|\mathbf{x}_{0i}\|_{2}\in[\sqrt{\delta},\sqrt{1-\delta}]$ for all $i\in[n]$, and $\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}\in[\delta,1-\delta]$ for all $i,j\in[n]$, for a constant $\delta\in(0,1/2)$. 3. (c)

The eigenvalues of $(1/n)\sum_{i=1}^{n}\mathbf{x}_{0i}\mathbf{x}_{0i}^{\mathrm{T}}$, $\sigma_{0,1}\geq\sigma_{0,2}\geq\ldots\geq\sigma_{0,d}$, satisfy either $\sigma_{0,k}=\sigma_{0,k+1}$ or $\sigma_{0,k}-\sigma_{0,k+1}>\lambda$ where $1\leq k\leq d-1$, and $\sigma_{0,d}>\lambda$, for a positive constant $\lambda$ for all $n\geq{}d$. 4. (d)

$\rho_{n}\in(0,1]$ for all $n$, $\lim_{n\to\infty}\rho_{n}$ exists with $(\log n)/(n\rho_{n})\to 0$ as $n\to\infty$. 5. (e)

$\delta^{2}<\tau_{n}/\rho_{n}<\delta/2$ for all $n$. 6. (f)

The first $p$ columns of $\mathbf{X}_{0}$ are orthogonal to the last $q$ columns of $\mathbf{X}_{0}$, where $\mathbf{X}_{0}=[\mathbf{x}_{01},\ldots,\mathbf{x}_{0n}]^{\mathrm{T}}$. 7. (g)

$\mathbf{A}\sim\mathrm{GRDPG}(\rho_{n}^{1/2}\mathbf{X}_{0},\,{\mathbf{I}}_{p,q})$.

In Assumption 1 above, items (b) and (c) are standard, and item (d) is a weak requirement on the network sparsity (also see Xie,, 2024). Item (f) can be made without loss of generality by Sylvester’s law of inertia. Item (e) guarantees that the true values of $\mathbb{E}_{0}A_{ij}$’s stay inside the truncated interval $[\tau_{n},1-\tau_{n}]$ and requires that the truncation level $\tau_{n}$ is not too small.

Assumption 2*.*

The prior densities $\pi_{i}(\mathbf{x}_{i})$, $i\in[n]$, which are independent, satisfy the following conditions, where $C,c>0$ are absolute constants:

1. (a)

$0<\pi_{i}(\mathbf{x}_{i})\leq C$, for all $\mathbf{x}_{i}\in\mathbb{R}^{d}$, and $\pi_{i}(\rho_{n}^{1/2}\mathbf{W}\mathbf{x}_{0i})\geq c$, for all $i\in[n]$; 2. (b)

$\|\partial\log\pi_{i}/\partial\mathbf{x}_{i}(\rho_{n}^{1/2}\mathbf{W}\mathbf{x}_{0i})\|_{2}\leq C$, for all $i\in[n]$; 3. (c)

$\log\pi_{i}(\mathbf{x}_{i})$ is concave in $\mathbf{x}_{i}$, and $\|\partial^{2}\log\pi_{i}/\partial\mathbf{x}_{i}\partial\mathbf{x}_{i}^{\mathrm{T}}(\mathbf{x}_{i})\|_{2}\leq C$, for all $\mathbf{x}_{i}\in\mathbb{R}^{d}$, for all $i\in[n]$.

Assumption 2 above lists several standard requirements for the prior distribution, such as the prior thickness in a neighborhood of the truth and log-concavity. In particular, $\pi_{i}(\mathbf{x}_{i})$ can be taken as a multivariate Gaussian distribution with a bounded mean vector and a covariance matrix (in spectra).

We now establish the asymptotic properties of the maximum extended surrogate likelihood estimator (MESLE). The MESLE provides the theoretical foundations for the Gaussian VI, and it is also theoretically appealing by itself. For convenience, denote by $\mathbf{G}_{in}(\mathbf{x}_{i})=(1/n)\sum_{j=1}^{n}{\rho_{n}\mathbf{x}_{0j}\mathbf{x}_{0j}^{\mathrm{T}}}/\{\rho_{n}^{1/2}\mathbf{x}_{i}^{\mathrm{T}}\mathbf{x}_{0j}(1-\rho_{n}^{1/2}\mathbf{x}_{i}^{\mathrm{T}}\mathbf{x}_{0j})\}$ and let $\mathbf{G}_{0in}=\mathbf{G}_{in}(\rho_{n}^{1/2}\mathbf{x}_{0i})$. Note that $\mathbf{G}_{0in}$ is precisely the Fisher information matrix for $\mathbf{x}_{i}$ at $\rho_{n}^{1/2}\mathbf{x}_{0i}$.

Theorem 3.2**.**

Suppose Assumption 1 holds. For each $i\in[n]$, let $\widehat{\mathbf{x}}_{i}=\operatorname*{arg\,max}_{\mathbf{x}_{i}\in\mathbb{R}^{d}}\widehat{\ell}_{in}(\mathbf{x}_{i})$ be the maximizer of the ESL function. Then, there exists an orthogonal matrix $\mathbf{W}\in\mathbb{O}(d)$ depending on $n$, and for any $c>0$, there exist a constant integer $N_{c,\delta,\lambda}$ and a constant $C_{c,\delta,\lambda}$ depending on $c,\,\delta,\,\lambda$, such that for all $n\geq N_{c,\delta,\lambda}$,

[TABLE]

and $\sqrt{n}\mathbf{G}_{0in}^{1/2}(\mathbf{W}^{\mathrm{T}}\widehat{\mathbf{x}}_{i}-\rho_{n}^{1/2}\mathbf{x}_{0i})\overset{\mathcal{L}}{\to}\mathrm{N}_{d}({\bm{0}}_{d},\,{\mathbf{I}}_{d})$ as $n\to\infty$. If furthermore $(\log n)^{4}/(n\rho_{n})\to 0$ as $n\to\infty$, then $\sum_{i=1}^{n}\|\mathbf{W}^{\mathrm{T}}\widehat{\mathbf{x}}_{i}-\rho_{n}^{1/2}\mathbf{x}_{0i}\|_{2}^{2}-\sum_{i=1}^{n}\mathrm{tr}(\mathbf{G}_{0in}^{-1})/n\overset{{\mathbb{P}}}{\to}0$ as $n\to\infty$.

Theorem 3.3 below is a Bernstein-von-Mises theorem for the posterior distribution (3.1) associated with the ESL function. It says that the posterior distribution converges in total variation distance to a normal distribution centered at the MESLE with covariance being the inverse Fisher information matrix scaled by $1/{n}$. For technical considerations, we impose the condition $\rho_{n}=1$, although it is possible to relax it and let $\rho_{n}\to 0$ as $n\to\infty$.

Theorem 3.3**.**

Suppose Assumption 1 and Assumption 2 hold, and assume $\rho_{n}=1$ for all $n$. For each $i\in[n]$, let $\widehat{\mathbf{x}}_{i}=\operatorname*{arg\,max}_{\mathbf{x}_{i}\in\mathbb{R}^{d}}\widehat{\ell}_{in}(\mathbf{x}_{i})$ be the maximizer of the ESL function (MESLE), and let $\mathbf{W}$ be the orthogonal alignment matrix in Theorem 3.2. Then, with probability at least $1-n^{-c}$,

[TABLE]

Since the variational posterior distribution is a minimizer of the KL divergence from the family of Gaussian distributions to the true posterior distribution, intuitively, the distance between the variational posterior distribution (i.e., solution to the problem (3.2)) and the posterior distribution defined in (3.1) should be small. With Theorem 3.3, we can make this intuition precise and establish the following Bernstein-von Mises theorem for VI.

Theorem 3.4**.**

Suppose the conditions in Theorem 3.3 hold. Let

[TABLE]

be the variational posterior distribution, where $\mathcal{Q}_{d}$ denotes the family of all $d$-dimensional Gaussian distributions. Then

[TABLE]

with probability at least $1-n^{-c}$.

We also provide the asymptotic normality of the variational posterior mean as a point estimator in Theorem 3.5 below.

Theorem 3.5**.**

Suppose the conditions in Theorem 3.3 hold. Let $\mathbf{x}_{i}^{*}$ be the variational posterior mean of $q_{in}^{*}(\mathbf{x}_{i})$. Then, $\sqrt{n}\mathbf{G}_{0in}^{1/2}(\mathbf{W}^{\mathrm{T}}\mathbf{x}_{i}^{*}-\rho_{n}^{1/2}\mathbf{x}_{0i})\overset{\mathcal{L}}{\to}\mathrm{N}_{d}({\bm{0}}_{d},\,{\mathbf{I}}_{d})$ as $n\to\infty$.

4 Numerical examples

In this section, we study the finite-sample numerical performance of SANVI in several simulated examples of GRDPG. For comparison, we implement the following competing estimates: ASE, OSE developed by Xie and Xu, (2023), Bayes estimate (BE) as the posterior mean from MCMC with the ESL function, and SANVI. For both BE and SANVI, the improper uniform prior distribution over the Euclidean space is used. We evaluate the performance of an estimator $\widehat{\mathbf{X}}$ by computing the sum of squared errors (the global error) between the aligned estimated latent positions and the true value counterparts, defined as $\mathrm{SSE}(\widehat{\mathbf{X}},\,\mathbf{X}_{0})=\inf_{\mathbf{W}\in\mathbb{O}(d)}\|\widehat{\mathbf{X}}\mathbf{W}-\mathbf{X}_{0}\|_{\mathrm{F}}^{2}$. Besides the simulated examples, we also apply the proposed method to a real-world graph dataset. We then discuss briefly the computation time of the proposed algorithm.

We consider four examples of GRDPG to investigate the numerical performance of the proposed estimate in various scenarios: a rank-two stochastic block model, a rank-two degree-corrected stochastic block model, a generic rank-two RDPG, and a generic rank-three GRDPG. For each example, several sample sizes are considered: $n=$ 1000, 3000, 5000, 7000, 9000, and 10000. We take the truncation parameter in the ESL function to be $\tau=\mathrm{min}(0.001,e^{1.5}/n)$, where $n$ denotes the number of vertices in the graph. For the parameters in the stochastic gradient descent for SANVI, we take the batch size $s=2$, the step size $\alpha_{0}=0.01$, decay rates for the moments of gradients $\beta_{1}=0.01$, $\beta_{2}=0.95$, and the maximum number of iterations to be $1000$. For the MCMC sampling, we use the Metropolis-Hastings algorithm with a Gaussian random walk. The total length of the chain is set to be $3000$, with a thinning of 2 and then a burn-in of 500, giving a set of $1000$ draws to compute the posterior mean. The covariance of the random walk is tuned so that the acceptance rates of the chains for most vertices among the $n$ vertices of a graph lie between $20\%$ and $30\%$. The experiments are repeated for $100$ times for each simulated scenario.

4.1 A stochastic block model example

As the first simulated example, consider a rank-two stochastic block model in the context of random dot product graphs with five blocks with latent positions $\mathbf{v}_{1}=[0.3,0.3]$, $\mathbf{v}_{2}=[0.5,0.5]$, $\mathbf{v}_{3}=[0.7,0.7]$, $\mathbf{v}_{4}=[0.3,0.7]$, $\mathbf{v}_{5}=[0.7,0.3]$. Each vertex is randomly assigned to one of the five blocks with equal probability. Arrange the five latent positions as the rows of a $5\times 2$ matrix $\mathbf{B}$, and let $\mathbf{Z}$ be an $n\times{}5$ matrix whose $i$th row $\mathbf{z}_{i}^{\mathrm{T}}$ encodes the block membership of vertex $i$, i.e., the $k$th entry of $\mathbf{z}_{i}$ is $1$ if vertex $i$ belongs to block $k$ and [math] otherwise. Conditional on the block assignments of all the vertices, we have $A_{ij}{\overset{\mathrm{i.i.d.}}{\sim}}\mathrm{Bernoulli}(\mathbf{v}_{k_{i}}^{\mathrm{T}}\mathbf{v}_{k_{j}})$ for all $i,j\in[n]$, and $\mathbf{A}\sim\mathrm{RDPG}(\mathbf{Z}\mathbf{B})$.

The sums of squared errors of the four estimates are summarized in Table 1. While OSE is numerically unstable and has sums of squared errors larger than all other estimates due to the latent position $\mathbf{v}_{3}=[0.7,\,0.7]$ being close to the unit circle (two vertices in this block have an edge probability of $0.98$), BE and SANVI are nevertheless numerically stable and perform better than ASE. The paired two-sample $t$-tests between the sums of squared errors of ASE and those of BE, and those of SANVI, respectively, are performed, and the $p$-values are listed in Table 2, from which we can see that the two estimates that are based on the ESL function indeed have smaller sums of squared errors.

4.2 A degree-corrected stochastic block model example

In this example, we consider a rank-two degree-corrected stochastic block model in the context of RDPG with two blocks. Specifically, let $\mathbf{v}_{1}=[{3\sqrt{10}}/{10},{\sqrt{10}}/{10}]$ and $\mathbf{v}_{2}=[{\sqrt{10}}/{10},{3\sqrt{10}}/{10}]$. Each vertex is randomly assigned to a block with equal probability and then assigned a degree-corrected parameter (weight) $\theta_{i}$ that follows $\mathrm{Uniform}(0.05,0.95)$. Arrange the latent positions of the two blocks as the rows of a $2\times 2$ matrix $\mathbf{B}$, let $\mathbf{Z}$ be an $n\times{}2$ matrix whose $i$th row $\mathbf{z}_{i}^{\mathrm{T}}$ encodes the block membership of vertex $i$, i.e., the $k$th entry of $\mathbf{z}_{i}$ is $1$ if vertex $i$ belongs to block $k$ and [math] otherwise, and let ${\bm{\Theta}}$ be an $n\times{}n$ diagonal matrix whose $(i,i)$th entry $\theta_{i}$ is the degree-corrected parameter of vertex $i$. Then, conditional on the block assignments and the degree corrections of all the vertices, we have the edge indicator $A_{ij}{\overset{\mathrm{i.i.d.}}{\sim}}\mathrm{Bernoulli}(\theta_{i}\theta_{j}\mathbf{v}_{k_{i}}^{\mathrm{T}}\mathbf{v}_{k_{j}})$ for all $i,j\in[n]$, and $\mathbf{A}\sim\mathrm{RDPG}({\bm{\Theta}}\mathbf{Z}\mathbf{B})$.

The sums of squared errors of the four estimates are summarized in Table 3. We can see that for large samples, the three likelihood-based estimates (OSE, BE, and SANVI) all have smaller sums of squared errors than ASE does, and in the case of $n=1000$, the two estimates based on the ESL function still perform well. This phenomenon empirically validates the statement that the likelihood-based estimates improve upon spectral-based ASE, since the latter does not incorporate likelihood information in the graph. The $p$-values given by the paired $t$-tests on the sums of squared errors of ASE and the other three estimates are listed in Table 4, which quantitatively verifies the smaller errors given by the likelihood-based estimates.

4.3 A two-dimensional latent curve example

Now we consider a rank-two generic RDPG whose latent positions are drawn from a latent curve in $\mathbb{R}^{2}$, parameterized as $[0.15\sin(\pi{}t)+0.6,\,0.15\cos(\pi{}t)+0.6]^{\mathrm{T}}$, for $0<t\leq{}1$, where the $n$ latent positions $\mathbf{x}_{i}$ for $i\in[n]$ are then obtained by taking $t={i}/{n}$ for $i\in[n]$. Then, we take $A_{ij}{\overset{\mathrm{i.i.d.}}{\sim}}{}\mathrm{Bernoulli}(\mathbf{x}_{0i}^{\mathrm{T}}\mathbf{x}_{0j})$ for all $i,j\in[n]$. Writing $\mathbf{X}_{0}=[\mathbf{x}_{1},\,\ldots,\,\mathbf{x}_{n}]^{\mathrm{T}}$, we then have $\mathbf{A}\sim\mathrm{RDPG}(\mathbf{X}_{0})$.

The sums of squared errors of the four estimates are summarized in Table 5. We can see that the errors are obviously larger than those in the previous simulated examples, due to the complex nature of the latent positions obtained from a curve, in comparison to the relatively simple nature of the latent positions in a stochastic block model. Similar to the case in the stochastic block model example, while OSE is numerically unstable in the presence of latent positions close to the unit circle, the two estimates based on the ESL function are nevertheless numerically stable and perform relatively well compared to ASE. From Table 4.3, we observe that the sums of squared errors of BE and SANVI are approximately $5\%$ to $15\%$ less than those of ASE. This indicates that the likelihood information in the ESL function facilitates the estimation of latent positions with complex structures. The $p$-values given by the paired two-sample $t$-tests between the sums of squared errors of ASE and those of the two ESL-based estimates are listed in Table 6, which quantitatively supports this observation.

4.4 A three-dimensional latent curve example

In this example, we consider a rank-three generic GRDPG, with signature $(2,1)$, whose latent positions are drawn from a latent curve in $\mathbb{R}^{3}$, parameterized as $[0.15\sin(2\pi{}t)+0.6,\,0.15\cos(2\pi{}t)+0.6,\,0.15\cos(4\pi{}t)]^{\mathrm{T}}$, for $0<t\leq{}1$, where the $n$ latent positions are then obtained by taking $t={i}/{n}$ for $i\in[n]$. In particular, the resulting edge probability matrix $\mathbf{P}=\mathbf{X}_{0}{\mathbf{I}}_{2,1}\mathbf{X}_{0}^{\mathrm{T}}$, where $\mathbf{X}_{0}=[\mathbf{x}_{1},\,\ldots,\,\mathbf{x}_{n}]^{\mathrm{T}}$ and ${\mathbf{I}}_{2,1}=\mathrm{diag}(1,\,1,\,-1)$, is an indefinite matrix. We then take $A_{ij}{\overset{\mathrm{i.i.d.}}{\sim}}{}\mathrm{Bernoulli}(\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{2,1}\mathbf{x}_{0j})$ for all $i,j\in[n]$, and $\mathbf{A}\sim\mathrm{GRDPG}(\mathbf{X}_{0})$ with signature $(2,1)$.

The sums of squared errors of the four estimates are summarized in Table 7. As in the previous example of the rank-two latent curve, the errors are relatively large due to the complex nature of the latent positions obtained from a curve. BE and SANVI still give relatively smaller sums of squared errors compared to ASE, and the $p$-values from the paired two-sample $t$-tests of their sums of squared errors are listed in Table 8.

4.5 Analysis of a real-world graph dataset

We finally apply the proposed algorithm on a real-world network of political blogs (Adamic and Glance,, 2005). The network corresponds to the hyperlinks of blogs regarding U.S. politics after the 2004 presidential election. These blogs are manually classified as either liberal or conservative, which we use as the true value of community labels. After following the rule of thumb by extracting the largest connected component and converting the resulting network with undirected edges, we obtain an $1224\times{}1224$ adjacency matrix with 33430 entries being 1 and others being 0. We apply MCMC and stochastic gradient Algorithm in 1 to compute the BE and SANVI estimates, together with ASE and OSE as the competitors. We choose the embedding dimension to be $d=2$ since there are two true communities in the network. These latent position estimates are then applied to the Gaussian-mixture-model-based clustering to obtain the set of estimated community labels, which we compare against the true community labels via the adjusted Rand index (ARI). The results are listed in Table 9, along with the corresponding computation time. Clearly, the two estimates that are based on the ESL function, i.e., BE and SANVI, are more accurate in terms of recovering the liberal-versus-conservative community structure in the network of these political blogs.

4.6 Discussion on computation time

A graph that has $n$ vertices gives an $n\times{}n$ adjacency matrix that has $n^{2}$ entries. Estimating GRDPG involves finding a $d$-dimensional representation for each vertex, resulting in an $n\times{}d$ latent position matrix. We decompose this into $n$ subproblems, each of which finds a $d$-dimensional latent position. Each subproblem is $O(n)$ in time, so the entire problem is $O(n^{2})$ in time.

Among the four estimates considered above, ASE and OSE require little time in computation, since the former is just the spectral decomposition truncated at the first $d$ dimensions of the adjacency matrix, and the latter is just a one-step update of the former. Methods based on the likelihood function typically require optimization and/or sampling, which are often computationally intensive. MCMC sampling is useful and often yields good results in various statistical problems; however, it is known to suffer from a long mixing time in some high-dimensional cases. Identifying a stopping criterion for an MCMC sampler is also nontrivial. VI tries to deal with this issue by turning the sampling problem into an optimization problem that requires relatively less computation time.

Specialized to the simulated examples above, while both BE and SANVI perform relatively as well as each other, VI requires less computation time than MCMC sampling. The relationship of computation time and sample size in the example of the stochastic block is given in Figure 1, with two quadratic curves fitted for the points corresponding to the MCMC sampler for BE and the stochastic gradient descent algorithm for SANVI, respectively. We can see that although both algorithms are $O(n^{2})$ in time, the optimization-based stochastic gradient descent requires around only $20\%$ of that of the sampling-based MCMC algorithm.

5 Discussion

In this paper, we propose an ESL function for GRDPG and leverage it to develop a computationally efficient spectral-assisted network variational inference method (SANVI). We establish the asymptotic properties of the point estimator, namely, MESLE, including its existence and uniqueness, its consistency at the rate of $\sqrt{n}$ with high probability, its asymptotic normality and efficiency, and, globally for all vertices, the consistency of the global error. The MESLE is theoretically interesting and computationally appealing in its own right. For SANVI, we establish the Bernstein-von Mises theorem of the variational posterior distribution and the asymptotic normality and efficiency of the variation inference estimator.

We also provide a stochastic gradient descent algorithm for implementing the computation of SANVI. Numerical study shows that, measured in terms of the global error, the point estimate of SANVI (variational posterior mean) is numerically comparable to BE, the latter of which is computed via a classical MCMC sampler. The computation time of the stochastic gradient descent algorithm for SANVI, although still in $O(n^{2})$ and is the same as the MCMC sampler due to the intrinsic properties of our setting, is indeed much less than the computation time of the classical MCMC algorithm.

Appendix A Preliminary Results

This section contains some preliminary results that will be used in the proofs of the main results.

Theorem A.1**.**

Suppose Assumption 1 holds. Let $\widetilde{\mathbf{X}}$ denote the signature-adjusted adjacency spectral embedding. Then

[TABLE]

where, for any $c>0$, there exists a constant $N_{c,\delta,\lambda}\in\mathbb{N}_{+}$ depending on $c,\delta,\lambda$, such that for all $n\geq N_{c,\delta,\lambda}$,

[TABLE]

with probability at least $1-n^{-c}$.

Remark 2*.*

Theorem A.1 is a generalization of Corollary 4.1 in Xie, (2024) to generalized random dot product graphs in our settings. The proof is mostly identical to its original version, with slight modifications such as the presence of the signature matrix ${\mathbf{I}}_{p,q}$ and the different definition of the orthogonal alignment matrix $\mathbf{W}$. We provide a proof here. For more theory on the entrywise limit theorems for eigenvectors of signal-plus-noise matrix models with weak signals, please refer to Xie, (2024).

Proof.

We clarify some notations first. Let $\mathbf{X}_{0+}$ denote the first $p$ columns of $\mathbf{X}_{0}$ (those corresponding to the positive part of the signature), and $\mathbf{X}_{0-}$ the last $q$ columns of $\mathbf{X}_{0}$ (those corresponding to the negative part of the signature), that is, $\mathbf{X}_{0}=[\mathbf{X}_{0+},\;\mathbf{X}_{0-}]$. Define $\Delta_{n}=(1/n)\mathbf{X}_{0}^{\mathrm{T}}\mathbf{X}_{0}$, then by the assumption that $\mathbf{X}_{0+}$ is orthogonal to $\mathbf{X}_{0-}$, we have

[TABLE]

Perform the eigendecomposition of $\mathbf{P}=\rho_{n}\mathbf{X}_{0}{\mathbf{I}}_{p,q}\mathbf{X}_{0}^{\mathrm{T}}$, we have $\mathbf{P}=\mathbf{U}_{\mathbf{P}}\mathbf{S}_{\mathbf{P}}\mathbf{U}_{\mathbf{P}}^{\mathrm{T}}$. Group by positive eigenvalues and negative eigenvalues, we have

[TABLE]

Let $\mathbf{W}_{\mathbf{X}_{0+}}\in\mathbb{O}(p)$ and $\mathbf{W}_{\mathbf{X}_{0-}}\in\mathbb{O}(q)$ be the orthogonal matrices such that $\rho_{n}^{1/2}\mathbf{X}_{0+}=\mathbf{U}_{\mathbf{P}+}\mathbf{S}_{\mathbf{P}+}^{1/2}\mathbf{W}_{\mathbf{X}_{0+}}$ and $\rho_{n}^{1/2}\mathbf{X}_{0-}=\mathbf{U}_{\mathbf{P}-}(-\mathbf{S}_{\mathbf{P}-})^{1/2}\mathbf{W}_{\mathbf{X}_{0-}}$. It is easy to see that $\rho_{n}\mathbf{X}_{0}{\mathbf{I}}_{p,q}\mathbf{X}_{0}^{\mathrm{T}}=\mathbf{U}_{\mathbf{P}}\mathbf{S}_{\mathbf{P}}\mathbf{U}_{\mathbf{P}}^{\mathrm{T}}$. Let $\hat{\lambda}_{1}\geq\hat{\lambda}_{2}\geq\ldots\geq\hat{\lambda}_{d}$ be the first $d$ eigenvalues of $\mathbf{A}$ that are largest in absolute value, and let $\hat{\sigma}_{1}\geq\hat{\sigma}_{2}\geq\ldots\geq\hat{\sigma}_{d}$ be the first $d$ singular values of $\mathbf{A}$. Note that the numbering for $\hat{\lambda}_{k}$ and $\hat{\sigma}_{k}$ are different in order. Recall the definition of the adjacency spectral embedding $\breve{\mathbf{X}}=\mathbf{U}_{\mathbf{A}}|\mathbf{S}_{\mathbf{A}}|^{1/2}$, and the definition of the signature-adjusted adjacency spectral embedding $\widetilde{\mathbf{X}}=\mathbf{U}_{\mathbf{A}}|\mathbf{S}_{\mathbf{A}}|^{1/2}\mathrm{sgn}(\mathbf{S}_{\mathbf{A}})$, where $\mathbf{S}_{\mathbf{A}}$ is the diagonal matrix with $\hat{\lambda}_{k}$, $k\in[d]$, arranged in the order of real numbers, $\mathrm{sgn}(\mathbf{S}_{\mathbf{A}})$ is the diagonal matrix with the signs (+1 and -1) of the corresponding eigenvalues, and $\mathbf{U}_{\mathbf{A}}$ is the matrix with the corresponding eigenvectors as columns. By grouping the positive eigenvalues and negative eigenvalues respectively, we can write

[TABLE]

Define $\breve{\mathbf{X}}_{+}=\mathbf{U}_{\mathbf{A}+}\mathbf{S}_{\mathbf{A}+}^{1/2}$ and $\breve{\mathbf{X}}_{-}=\mathbf{U}_{\mathbf{A}-}(-\mathbf{S}_{\mathbf{A}-})^{1/2}$, and $\widetilde{\mathbf{X}}_{+}=\mathbf{U}_{\mathbf{A}+}\mathbf{S}_{\mathbf{A}+}^{1/2}$ and $\widetilde{\mathbf{X}}_{-}=-\mathbf{U}_{\mathbf{A}-}(-\mathbf{S}_{\mathbf{A}-})^{1/2}$, we have

[TABLE]

Note that $\widetilde{\mathbf{X}}=\breve{\mathbf{X}}{\mathbf{I}}_{p,q}$. Let $\mathbf{U}_{\mathbf{P}+}^{\mathrm{T}}\mathbf{U}_{\mathbf{A}+}$ and $\mathbf{U}_{\mathbf{P}-}^{\mathrm{T}}\mathbf{U}_{\mathbf{A}-}$ yield the singular value decompositions $\mathbf{W}_{1+}{\bm{\Sigma}}_{+}\mathbf{W}_{2+}^{\mathrm{T}}$ and $\mathbf{W}_{1-}{\bm{\Sigma}}_{-}\mathbf{W}_{2-}^{\mathrm{T}}$, respectively, and define $\mathbf{W}_{+}^{*}=\mathbf{W}_{1+}\mathbf{W}_{2+}^{\mathrm{T}}$ and $\mathbf{W}_{-}^{*}=\mathbf{W}_{1-}\mathbf{W}_{2-}^{\mathrm{T}}$, and let $\mathbf{W}^{*}=\mathrm{diag}(\mathbf{W}_{+}^{*},\;\mathbf{W}_{-}^{*})$. Then the orthogonal alignment matrix between $\breve{\mathbf{X}}_{+}$ and $\rho_{n}^{1/2}\mathbf{X}_{0+}$ is selected as $(\mathbf{W}_{+}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0+}}$, and the same for that between $\widetilde{\mathbf{X}}_{+}$ and $\rho_{n}^{1/2}\mathbf{X}_{0+}$; the orthogonal alignment matrix between $\breve{\mathbf{X}}_{-}$ and $\rho_{n}^{1/2}\mathbf{X}_{0-}$ is selected as $(\mathbf{W}_{-}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0-}}$, and that between $\widetilde{\mathbf{X}}_{-}$ and $\rho_{n}^{1/2}\mathbf{X}_{0-}$ is selected as $-(\mathbf{W}_{-}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0-}}$. So the orthogonal alignment matrix between the adjacency spectral embedding $\breve{\mathbf{X}}$ and $\rho_{n}^{1/2}\mathbf{X}_{0}$ is the block diagonal matrix $\mathrm{diag}((\mathbf{W}_{+}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0+}},\;(\mathbf{W}_{-}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0-}})$, and that between $\widetilde{\mathbf{X}}$ and $\rho_{n}^{1/2}\mathbf{X}_{0}$ is $\mathrm{diag}((\mathbf{W}_{+}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0+}},\;-(\mathbf{W}_{-}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0-}})$. The $\mathbf{W}$ in the statement of this lemma is $\mathrm{diag}((\mathbf{W}_{+}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0+}},\;(\mathbf{W}_{-}^{*})^{\mathrm{T}}\mathbf{W}_{\mathbf{X}_{0-}})$ because we are aligning $\widetilde{\mathbf{X}}$ and $\rho_{n}^{1/2}\mathbf{X}_{0}{\mathbf{I}}_{p,q}$.

To prove the theorem, we follow the proofs in Xie, (2024). We first present a useful result for random graphs.

Result A.1*.*

Suppose Assumption 1 holds. Let $\mathbf{P}=\rho_{n}\mathbf{X}_{0}{\mathbf{I}}_{p,q}\mathbf{X}_{0}^{\mathrm{T}}$. Then for any $c>0$, there exists some constant $K_{c}>$ depending on $c$, such that $\|\mathbf{A}-\mathbf{P}\|_{2}\leq K_{c}(n\rho_{n})^{1/2}$ with probability at least $1-n^{-c}$. This follows exactly from Theorem 5.2 in Lei and Rinaldo, (2015).

We need to verify the Assumptions 1-5 in Xie, (2024) for our setup. By our definition of generalized random dot product graphs, Assumptions 1-3 in Xie, (2024) automatically holds. We now verify Assumptions 4 in Xie, (2024). Fix an arbitrary constant $c\geq 1$. Write $\mathbf{A}=\mathbf{P}+\mathbf{E}$. Let $\mathbf{e}_{i}$ denote the unit basis vector whose $i$th coordinate is one and the rest of coordinates are zeros. Let $\mathbf{E}^{(m)}$ denote the matrix constructed by replacing the $m$th row and $m$th column of $\mathbf{E}$ by their expected values which are zeros. Define the function $\phi(x)=(2+\beta_{c})(\max(\log(1/x),\,1))^{-1}\lambda_{d}({\bm{\Delta}}_{n})^{-1}$ for a constant $\beta_{c}>0$ that satisfies $\beta_{c}n\rho_{n}\geq(c+2)\log{n}$. By Lemma S6.1 in Xie, (2024), for any deterministic $\mathbf{V}\in\mathbb{R}^{n\times d}$, we have

[TABLE]

where $\xi=c\geq 1$ and $c_{0}=2$. To show that the same concentration bound also holds for $\|\mathbf{e}_{i}^{\mathrm{T}}\mathbf{E}^{(m)}\mathbf{V}\|_{2}$, we simply observe that $[\mathbf{E}^{(m)}]_{im}$ can be viewed as a centered Bernoulli random variable whose success probability is zero. Then applying Lemma S6.1 in Xie, (2024) leads to

[TABLE]

where $\xi=c\geq 1$ and $c_{0}=2$. We now verify Assumptions 5 in Xie, (2024). By Result A.1, there exists a constant $K_{c}\geq 1$ that depends on $c$ such that ${\mathbb{P}}(\|\mathbf{E}\|_{2}\leq K_{c}(n\rho_{n})^{1/2})\geq 1-n^{-c}$. Let $\kappa({\bm{\Delta}}_{n})=\lambda_{1}({\bm{\Delta}}_{n})/\lambda_{d}({\bm{\Delta}}_{n})$ be the condition number of ${\bm{\Delta}}_{n}$. Then with

[TABLE]

we immediately see that

[TABLE]

which shows that the Assumption 5 in Xie, (2024) holds with $\zeta=c\geq 1$ and $c_{0}=1$. The five Assumptions in Xie, (2024) are thus verified for our setting.

Write

[TABLE]

which can be view as a sum of two terms. Then by Lemma S2.1 in Xie, (2024) we have

[TABLE]

for all $n\geq N_{c}$ that depends on $c$ with probability at least $1-n^{-c}$, and by Theorem 3.2 in Xie, (2024) we have

[TABLE]

for all $n\geq N_{c}$ that depends on $c$ with probability at least $1-n^{-c}$, and we also have

[TABLE]

So, with the assumption that $(\log n)/(n\rho_{n})\to 0$, we have

[TABLE]

for all $n\geq N_{c}$ that depends on $c,\lambda$ with probability at least $1-n^{-c}$. ∎

Lemma A.2** (Some frequently used results).**

Suppose Assumption 1 holds. Let $\breve{\mathbf{X}}$ denote the adjacency spectral embedding, and $\widetilde{\mathbf{X}}$ the signature-adjusted adjacency spectral embedding. Let $p_{0ij}=\rho_{n}\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}$, and $\widetilde{p}_{ij}=\breve{\mathbf{x}}_{i}^{\mathrm{T}}\widetilde{\mathbf{x}}_{j}$, $i,j\in[n]$. Then for any $c>0$, there exists a constant $N_{c,\delta,\lambda}\in\mathbb{N}_{+}$ depending on $c,\delta,\lambda$, such that for all $n\geq N_{c,\delta,\lambda}$, the following hold with probability at least $1-n^{-c}$:

[TABLE]

Proof.

We prove the results one by one. For simplicity of notation, in the proof of this lemma, the results are stated to hold with probability at least $1-n^{-c}$ for all $n\geq{}{N}_{c,\delta,\lambda}$ for some large constant integer ${N}_{c,\delta,\lambda}$ that depends on $c,\delta,\lambda$, where $c>0$ is an arbitrary positive constant. Also, the results that hold for a single $i\in[n]$ with probability at least $1-n^{-c}$ can be strengthened to hold for all $i\in[n]$ by taking a union bound over $i\in[n]$.

For (a), by assumption we have $(\log n)/(n\rho_{n})\to 0$, so we can pick an ${N}_{c,\delta,\lambda}$ large enough such that for all $n\geq{}{N}_{c,\delta,\lambda}$, we have $C_{c,\delta,\lambda}(\log n)/(n\rho_{n})<1-\delta/2-\sqrt{1-\delta}$ (this is because $(1-\delta/2)^{2}=1-\delta+\delta^{2}/4>1-\delta$ and recall that $\delta\in(0,1/2)$), and we also have $C_{c,\delta,\lambda}(\log n)/(n\rho_{n})<\sqrt{\delta}-\delta/2$ (this is because $1-\delta/2-\sqrt{1-\delta}=1-\sqrt{1-\delta}-\sqrt{\delta}+\sqrt{\delta}-\delta/2=1-\sqrt{(\sqrt{1-\delta}+\sqrt{\delta})^{2}}+\sqrt{\delta}-\delta/2=1-\sqrt{1+2\sqrt{\delta(1-\delta)}}+\sqrt{\delta}-\delta/2\leq\sqrt{\delta}-\delta/2$ and recall that $\delta\in(0,1/2)$). Then by triangle inequality and Theorem A.1,

[TABLE]

For (b), by triangle inequality, Cauchy–Schwarz inequality, and Theorem A.1,

[TABLE]

For (c), similar to (a), we can pick an ${N}_{c,\delta,\lambda}$ such that for all $n\geq{}{N}_{c,\delta,\lambda}$, we have

$C_{c,\delta,\lambda}(\log n)/(n\rho_{n})<\delta/2$. Then by triangle inequality and the previous result,

[TABLE]

For (d), by triangle inequality, Cauchy–Schwarz inequality, and Theorem A.1,

[TABLE]

For (e), by triangle inequality, Cauchy–Schwarz inequality, and Theorem A.1,

[TABLE]

so

[TABLE]

With a union bound, the result holds with maximum over $i\in[n]$.

For (f), similar to (a), we can pick an ${N}_{c,\delta,\lambda}$ such that for all $n\geq{}{N}_{c,\delta,\lambda}$, we have

$C_{c,\delta,\lambda}(\log n)/(n\rho_{n})<\delta/2$. Then by triangle inequality and the previous result,

[TABLE]

With a union bound, the result holds with minimum over $i\in[n]$.

For (g), similar to (a), we can pick an ${N}_{c,\delta,\lambda}$ such that for all $n\geq{}{N}_{c,\delta,\lambda}$, we have

$C_{c,\delta,\lambda}(\log n)/(n\rho_{n})<\delta/2$. Then by triangle inequality and the previous result,

[TABLE]

With a union bound, the result holds with maximum over $i\in[n]$. ∎

Lemma A.3** (Some results with $A_{ij}$).**

Suppose Assumption 1 holds. Denote by $p_{0ij}=\rho_{n}\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}$. Let $\alpha_{ijn}$ be a two-dimensional array of real numbers such that $|\alpha_{ijn}|\leq C_{\delta}\alpha_{n}$ for all $n$ where $C_{\delta}$ is a constant that depends on $\delta$ and $\alpha_{n}$ is a function of $n$. Then for any $c>0$, there exists a constant $N_{c,\delta,\lambda}\in\mathbb{N}_{+}$ depending on $c,\delta,\lambda$, such that for all $n\geq N_{c,\delta,\lambda}$, the following hold with probability at least $1-n^{-c}$:

[TABLE]

With a union bound over $i\in[n]$, the results above hold with maximum over $i\in[n]$.

Proof.

For simplicity of notation, in the proof of this lemma, the results are stated to hold with probability at least $1-n^{-c}$ for all $n\geq{}{N}_{c,\delta,\lambda}$ for some large constant integer ${N}_{c,\delta,\lambda}$ that depends on $c,\delta,\lambda$, where $c>0$ is an arbitrary positive constant. For (a), by Bernstein’s inequality,

[TABLE]

Let $c>0$ be given, and by the assumption that $(\log n)/(n\rho_{n})\to 0$, we have $\sqrt{\log n}\leq C_{2}\sqrt{n\rho_{n}}$ for a constant $C_{2}$ for all sufficiently large $n$. Take $t=C_{1}C_{\delta}\alpha_{n}\rho_{n}^{1/2}\sqrt{(\log n)/n}$, where $C_{1}$ is a constant that depends on $c$ and $C_{2}$ that satisfies $-3C_{1}^{2}/(6+2C_{1}C_{2})<-(\log 2)/(\log n)-c$ for all sufficiently large $n$, we have

[TABLE]

so $|(1/n)\sum_{j=1}^{n}(A_{ij}-p_{0ij})\alpha_{ijn}|\lesssim_{c,\delta}\alpha_{n}\rho_{n}\sqrt{{\log n}/(n\rho_{n})}$ for all $n\geq N_{c}$ with probability at least $1-n^{-c}$. With a union bound, we can take maximum over $i\in[n]$ and the bound still holds.

For (b), by triangle inequality, the previous result, and the assumption that $(\log n)/(n\rho_{n})\to 0$ as $n\to\infty$,

[TABLE]

For (c),

[TABLE]

For (d), by the assumption that $\|\mathbf{x}_{0j}\|_{2}\in[\sqrt{\delta},\sqrt{1-\delta}]$ for all $j\in[n]$, and the previous result (a), we have

[TABLE]

in which we note that the chosen embedding dimension $d$ implicitly depends on $\lambda$.

For (e), similar to (c), we have

[TABLE]

For (f), we have

[TABLE]

for all $n$ large enough such that $C_{c,\delta,\lambda}\sqrt{(\log n)/(n\rho_{n})}<(1/2)\delta\lambda$, by assumption and by the previous result (e).

For (g), with the previous result (f), we have

[TABLE]

where

[TABLE]

by the previous result (b) and Lemma A.2, so $(1/n)\sum_{j=1}^{n}A_{ij}\widetilde{\mathbf{x}}_{j}\widetilde{\mathbf{x}}_{j}^{\mathrm{T}}\succeq(1/4)\delta\lambda\rho_{n}^{2}{\mathbf{I}}_{d}$ for all $n$ large enough such that $C_{c,\delta,\lambda}\sqrt{(\log n)/(n\rho_{n})}<(1/4)\delta\lambda$. ∎

Lemma A.4** (Concentration of the gradient).**

Suppose Assumption 1 holds. Then for any $c>0$, there exists a constant $N_{c,\delta,\lambda}\in\mathbb{N}_{+}$ depending on $c,\delta,\lambda$, such that for all $n\geq N_{c,\delta,\lambda}$,

[TABLE]

with probability at least $1-n^{-c}$, where

[TABLE]

and $\widetilde{\mathbf{X}}$ denotes the signature-adjusted adjacency spectral embedding.

Proof.

In the proof of this lemma, the large probability bounds with probability at least $1-n^{-c}$ are stated with respect to all $n\geq{}{N}_{c,\delta,\lambda}$ for some large constant integer ${N}_{c,\delta,\lambda}$ that depends on $c,\delta,\lambda$, where $c>0$ is an arbitrary positive constant. Let $p_{0ij}=\rho_{n}\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}$, and let $g_{n}(t)=\psi_{n}^{\prime}(t)$ for simplicity of notation. Then it is easy to see that $g_{n}(t)>0$, that $g^{\prime}_{n}(t)=-\tau_{n}^{-2}\mathbbm{1}(t<\tau_{n})-t^{-2}\mathbbm{1}(t\in[\tau_{n},1])-\mathbbm{1}(t>1)$, which implies that $g_{n}(t)$ is decreasing in $t$, that $g_{n}^{\prime}(t)$ is constant on $t<\tau_{n}$ or $t>1$ and increasing on $t\in[\tau_{n},1]$, and that $-1\leq\tau_{n}^{2}g^{\prime}_{n}(t)\leq-\tau_{n}^{2}$.

Note that $\sqrt{1-\delta}<1-\delta/2$, which is because $1-\delta<1-\delta+\delta^{2}/4=(1-\delta/2)^{2}$ and recall that $\delta\in(0,1/2)$; and also note that $\delta/2\leq 1-(1-\delta/2)\rho_{n}$, which is because $\delta/2=(\delta/2)(1-\rho_{n}+\rho_{n})=(\delta/2)(1-\rho_{n})+(\delta/2)\rho_{n}\leq 1-\rho_{n}+(\delta/2)\rho_{n}=1-(1-\delta/2)\rho_{n}$; so we have $\tau_{n}<1-(1-\delta/2)\rho_{n}$, which is because $\tau_{n}<(\delta/2)\rho_{n}\leq\delta/2$ by assumption.

We have $\max_{j\in[n]}\sup_{\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}}|\rho_{n}^{1/2}\mathbf{x}_{i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}|\leq\sqrt{1-\delta}\rho_{n}<(1-\delta/2)\rho_{n}$ by assumption,

$\max_{j\in[n]}\sup_{\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}}|\mathbf{x}_{i}^{\mathrm{T}}\mathbf{W}^{\mathrm{T}}\widetilde{\mathbf{x}}_{j}|\leq\rho_{n}^{1/2}\max_{j\in[n]}\|\widetilde{\mathbf{x}}_{j}\|_{2}\leq(1-\delta/2)\rho_{n}$ for all $n\geq N_{c,\delta,\lambda}$ that depends on $c,\delta,\lambda$ with probability at least $1-n^{-c}$ by Lemma A.2,

$\max_{j\in[n]}\sup_{\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}}|\mathbf{x}_{i}^{\mathrm{T}}\mathbf{W}^{\mathrm{T}}\widetilde{\mathbf{x}}_{j}-\rho_{n}^{1/2}\mathbf{x}_{i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}|\leq C_{c,\delta,\lambda}\rho_{n}^{1/2}\sqrt{(\log n)/n}$ for all $n\geq N_{c,\delta,\lambda}$ and a constant $C_{c,\delta,\lambda}$ that depend on $c,\delta,\lambda$ with probability at least $1-n^{-c}$ by Theorem A.1, and $\max_{j\in[n]}\sup_{\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}}|\rho_{n}^{1/2}\mathbf{x}_{i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}+\theta(\mathbf{x}_{i}^{\mathrm{T}}\mathbf{W}^{\mathrm{T}}\widetilde{\mathbf{x}}_{j}-\rho_{n}^{1/2}\mathbf{x}_{i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j})|<(1-\delta/2)\rho_{n}$ for all $n\geq N_{c,\delta,\lambda}$ such that $C_{c,\delta,\lambda}\sqrt{(\log n)/(n\rho_{n})}<1-\delta/2-\sqrt{1-\delta}$ with probability at least $1-n^{-c}$ by the results above and by the assumption that $(\log n)/(n\rho_{n})\to 0$, where $\theta\in(0,1)$.

Write

[TABLE]

which can be viewed as a sum of six terms. For simplicity of notation, in the remaining of the proof of this lemma, the large probability bounds are stated with respect to all $n\geq{}{N}_{c,\delta,\lambda}$ for some large constant integer ${N}_{c,\delta,\lambda}$ that depends on $c,\delta,\lambda$.

For the first term, with probability at least $1-n^{-c}$,

[TABLE]

by Cauchy-Schwarz inequality, mean value theorem, the properties of the function $g_{n}(t)$, the assumption that $\delta^{2}\rho_{n}<\tau_{n}$, Theorem A.1, Lemma A.2, and Lemma A.3.

For the second term, with probability at least $1-n^{-c}$,

[TABLE]

by Cauchy-Schwarz inequality, mean value theorem, the properties of the function $g_{n}(t)$, the result that $\tau_{n}<(\delta/2)\rho_{n}\leq\delta/2\leq 1-(1-\delta/2)\rho_{n}$ shown above, Theorem A.1, and Lemma A.2.

For the third term, with probability at least $1-n^{-c}$,

[TABLE]

where the inequalities and equalities follow from Cauchy-Schwarz inequality, triangle inequality, the properties of the function $g_{n}(t)$, the assumption that $\delta^{2}<\tau_{n}/\rho_{n}<\delta/2$, Theorem A.1, Lemma A.2, and Lemma A.3.

For the fourth term, with probability at least $1-n^{-c}$,

[TABLE]

where the inequalities and equalities follow from Cauchy-Schwarz inequality, triangle inequality, the properties of the function $g_{n}(t)$, the assumption that $\delta^{2}<\tau_{n}/\rho_{n}<\delta/2$, and Theorem A.1.

For the fifth term and the sixth term, some methods in empirical processes are needed. We first define a stochastic process indexed by $\mathbf{x}_{i}$ in the closed ball that is centered at origin and of radius $\rho_{n}^{1/2}$, then use the results in Chapter 2.2 of van der Vaart and Wellner, (2023) to compute the bounds on the Orlicz $\psi_{1}$ norm for the supremum of the process, and then use a Bernstein-type inequality (Theorem 12.2 in Boucheron et al., (2013)) to obtain a tail probability bound for the supremum of the process. We now show the fifth term. Let $J_{ijk}(\mathbf{x}_{i})=(A_{ij}-p_{0ij})g_{n}(\rho_{n}^{1/2}\mathbf{x}_{i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j})\rho_{n}^{1/2}x_{0jk}$, and for each $k\in[d]$, define a stochastic processes $J_{ink}(\mathbf{x}_{i})=(1/n)\sum_{j=1}^{n}J_{ijk}(\mathbf{x}_{i})$, where $\mathbf{x}_{i}\in B({\bm{0}}_{d},\rho_{n}^{1/2})=\{\mathbf{x}_{i}\in\mathbb{R}^{d}:\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}\}$. Then we have $\mathbb{E}\left[J_{ijk}(\mathbf{x}_{i})-J_{ijk}(\mathbf{x}_{i}^{\prime})\right]=0$,

[TABLE]

and similarly, $\mathbb{E}|J_{ijk}(\mathbf{x}_{i})-J_{ijk}(\mathbf{x}_{i}^{\prime})|^{2}\leq\delta^{-8}\rho_{n}^{-1}\|\mathbf{x}_{i}-\mathbf{x}_{i}^{\prime}\|_{2}^{2}$ for all $\mathbf{x}_{i}\,\mathbf{x}_{i}^{\prime}\in B({\bm{0}}_{d},\rho_{n}^{1/2})$ and all $j\in[n]$. Then for any $\mathbf{x}_{i},\,\mathbf{x}_{i}^{\prime}\in B({\bm{0}}_{d},\rho_{n}^{1/2})$, by Bernstein’s inequality,

[TABLE]

where we take $C_{\delta}=\delta^{-4}$. We then consider the case where $|J_{ink}(\mathbf{x}_{i})-J_{ink}(\mathbf{x}_{i}^{\prime})|\leq 3C_{\delta}\|\mathbf{x}_{i}-\mathbf{x}_{i}^{\prime}\|_{2}$ (the sub-Gaussian part) and the case where $|J_{ink}(\mathbf{x}_{i})-J_{ink}(\mathbf{x}_{i}^{\prime})|>3C_{\delta}\|\mathbf{x}_{i}-\mathbf{x}_{i}^{\prime}\|_{2}$ (the sub-exponential part) separately, noting that both the two bounds hold for all $t>0$, i.e.,

[TABLE]

Recall that the Orlicz $\psi_{p}$ norm of a random variable $X$ is $\|X\|_{\psi_{p}}=\inf\{c>0:\mathbb{E}[\psi_{p}(|x|/c)]\leq 1\}$ with $\psi_{p}(x)=e^{x^{p}}-1$. Then, by sub-additivity, the fact that $\|X\|_{\psi_{1}}\leq(1/\sqrt{\log 2})\|X\|_{\psi_{2}}$ (Problem 2.2.5 in van der Vaart and Wellner, (2023)), and Lemma 2.2.1 in van der Vaart and Wellner, (2023), we can bound the Orlicz $\psi_{1}$ norm of $J_{ink}(\mathbf{x}_{i})-J_{ink}(\mathbf{x}_{i}^{\prime})$, i.e.,

[TABLE]

where $C_{\delta}=1/\delta^{4}$ (note that $C_{\delta}>1$) and $K(n)=\sqrt{(12/\log 2)}(n\rho_{n})^{-1/2}+4(n\rho_{n})^{-1}$.

Define a metric $d_{J}(\mathbf{x}_{i},\mathbf{x}_{i}^{\prime})=C_{\delta}K(n)\|\mathbf{x}_{i}-\mathbf{x}_{i}^{\prime}\|_{2}$ on $B({\bm{0}}_{d},\rho_{n}^{1/2})$, then the diameter of $B({\bm{0}}_{d},\rho_{n}^{1/2})$ under $d_{J}$ is $2C_{\delta}\rho_{n}^{1/2}{}K(n)$, and the packing number of the metric space $(B({\bm{0}}_{d},\rho_{n}^{1/2}),d_{J})$ satisfies $D(\epsilon,d_{J})\leq(2C_{\delta}\rho_{n}^{1/2}{}K(n)/\epsilon)^{d}$. Then by Corollary 2.2.5 of van der Vaart and Wellner, (2023),

[TABLE]

where $C$ is a constant (related to the function $\psi_{1}(x)=e^{x}-1$).

Consider $J_{ijk}({\bm{0}}_{d})=(A_{ij}-p_{0ij})2\tau_{n}^{-1}\rho_{n}^{1/2}x_{0jk}$. We then have $|J_{ijk}({\bm{0}}_{d})|\leq 2C_{\delta}^{1/2}\rho_{n}^{-1/2}$ and $\mathbb{E}|J_{ijk}({\bm{0}}_{d})|^{2}\leq 4C_{\delta}$ for all $j\in[n]$. Then by Bernstein’s inequality,

[TABLE]

from which, similar to the computation for $\|J_{ink}(\mathbf{x}_{i})-J_{ink}(\mathbf{x}_{i}^{\prime})\|_{\psi_{1}}$ shown above, we can obtain $\|J_{ink}({\bm{0}}_{d})\|_{\psi_{1}}\leq\sqrt{(48/\log 2)}C_{\delta}^{1/2}n^{-1/2}+8C_{\delta}^{1/2}n^{-1}\rho_{n}^{-1/2}=2C_{\delta}^{1/2}\rho_{n}^{1/2}K(n)$. By triangle inequality, monotonicity of integral, and sub-additivity of norm,

[TABLE]

and we also have

[TABLE]

since $\|X\|_{1}\leq\|X\|_{\psi_{1}}$ (recall that $x\leq e^{x}-1$ for $x\geq 0$).

Similar to the previous steps, it is straightforward to compute that $|(2C_{\delta})^{-1}\rho_{n}^{1/2}J_{ijk}(\mathbf{x}_{i})|<1$ for all $\mathbf{x}_{i}\in B({\bm{0}}_{d},\rho_{n}^{1/2})$ and for all $j\in[n]$, and also that

[TABLE]

Let $B^{\prime}({\bm{0}}_{d},\rho_{n}^{1/2})$ be the set of all rational points in the closed ball that is centered at origin and of radius $\rho_{n}^{1/2}$, then $B^{\prime}({\bm{0}}_{d},\rho_{n}^{1/2})$ is countable and is a dense subset of $B({\bm{0}}_{d},\rho_{n}^{1/2})$. Since $J_{ink}(\mathbf{x}_{i})$ is continuous in $\mathbf{x}_{i}$ surely, for any $\mathbf{x}_{i}^{*}\in B({\bm{0}}_{d},\rho_{n}^{1/2})$, there exists a sequence $\{\mathbf{x}_{i}^{(m)}\}\subset B^{\prime}({\bm{0}}_{d},\rho_{n}^{1/2})$ with $\lim_{m\to\infty}\mathbf{x}_{i}^{(m)}=\mathbf{x}_{i}^{*}$ such that $\lim_{m\to\infty}J_{ink}(\mathbf{x}_{i}^{(m)})=J_{ink}(\mathbf{x}_{i}^{*})$. So $J_{ink}(\mathbf{x}_{i})$ indexed by $\mathbf{x}_{i}\in B({\bm{0}}_{d},\rho_{n}^{1/2})$ is a separable process, and we have

[TABLE]

Then, by Theorem 12.2 in Boucheron et al., (2013), we have

[TABLE]

and also

[TABLE]

from both of which we have

[TABLE]

Take $t=8C_{t}C_{\delta}\sqrt{(\log n)/n}$, and recall the bound on $\mathbb{E}[\sup_{\mathbf{x}_{i}\in B({\bm{0}}_{d},\rho_{n}^{1/2})}|J_{ink}(\mathbf{x}_{i})|]$, we have

[TABLE]

For any $c>0$, let $N_{c,\delta,\lambda}$ be an integer large enough such that $\rho_{n}^{1/2}{}K(n)<\sqrt{(\log n)/n}$ for all $n\geq N_{c,\delta,\lambda}$ (recall that $K(n)=\sqrt{(12/\log 2)}(n\rho_{n})^{-1/2}+4(n\rho_{n})^{-1}$ and that we assume $(\log n)/(n\rho_{n})\to 0$), and then choose $C_{t}$ large enough such that $(\log 2)/(\log n)-2C_{t}^{2}/(1+C_{t}\sqrt{(\log n)/(n\rho_{n})})<-c$ for all $n\geq N_{c,\delta,\lambda}$ (it is decreasing in $C_{t}$ and in $n$), we have

[TABLE]

with probability at least $1-n^{-c}$ (the embedding dimension implicitly depends on $\lambda$).

By the fact that $\inf_{x}f(x)\leq\sup_{x}f(x)$ and $-\inf_{x}f(x)=\sup_{x}-f(x)$, the previous computation also gives (omitting some details)

[TABLE]

which gives

[TABLE]

with probability at least $1-n^{-c}$. Then by the fact that $\sup_{x}|f(x)|\leq|\sup_{x}f(x)|+|\inf_{x}f(x)|$, we have

[TABLE]

with probability at least $1-n^{-c}$. So we have

[TABLE]

with probability at least $1-n^{-c}$. This shows the bound for the fifth term in A.1.

By a similar approach, for the sixth term, we have (omitting the details to save space)

[TABLE]

with probability at least $1-n^{-c}$. The conclusion of the lemma then follows from applying triangle inequality and combining the six large probability bounds above, then with a union bound over $i\in[n]$. ∎

Lemma A.5** (Concentration of the Hessian matrix).**

Suppose Assumption 1 holds. Then for any $c>0$, there exists a constant integer $N_{c,\delta,\lambda}\in\mathbb{N}_{+}$ and a positive constant $C_{c,\delta,\lambda}$ depending on $c,\delta,\lambda$, such that for all $n\geq N_{c,\delta,\lambda}$,

[TABLE]

with probability at least $1-n^{-c}$, where $\varepsilon_{n}=C_{c,\delta,\lambda}\sqrt{(\log n)/n}$.

Proof.

For simplicity of notation, Let $p_{0ij}=\rho_{n}\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}$. A simple algebra shows that

[TABLE]

which can be viewed as a sum of six terms. For simplicity of notation, in the remaining of the proof of this lemma, the large probability bounds are stated with respect to all $n\geq{}{N}_{c,\delta,\lambda}$ for some large constant integer ${N}_{c,\delta,\lambda}$ that depends on $c,\delta,\lambda$.

For the first term, with probability at least $1-n^{-c}$,

[TABLE]

by from Cauchy–Schwarz inequality, the properties of the function $\psi_{n}(t)$, mean value theorem, $(\log n)/(n\rho_{n})\to 0$, and Lemma A.2.

For the second term, with probability at least $1-n^{-c}$,

[TABLE]

where the inequalities follow from Cauchy–Schwarz inequality, the properties of the function $\psi_{n}(t)$, mean value theorem, $(\log n)/(n\rho_{n})\to 0$, and Lemma A.2.

For the third term, with probability at least $1-n^{-c}$,

[TABLE]

where the inequalities follow from Cauchy–Schwarz inequality, triangle inequality, Lemma A.2, and Lemma A.3.

For the fourth term, note that $\left\{-\psi_{n}^{\prime\prime}(p_{0ij})+\psi_{n}^{\prime\prime}(1-p_{0ij})\right\}\rho_{n}=\{p_{0ij}^{-2}-(1-p_{0ij})^{-2}\}\rho_{n}\leq C_{\delta}\rho_{n}^{-1},$ so

[TABLE]

with probability at least $1-n^{-c}$ by Lemma A.3.

For the fifth term, with probability at least $1-n^{-c}$,

[TABLE]

by triangle inequality and Lemma A.2.

For the sixth term, since $\tau_{n}<\delta\rho_{n}$ for all $n$ and $p_{0ij}\in[\delta\rho_{n},(1-\delta)\rho_{n}]$ for all $i,j\in[n]$ by assumption, we have

[TABLE]

The conclusion follows from applying triangle inequality, combining the six bounds above, and applying a union bound over $i\in[n]$. ∎

Lemma A.6** (Lipschitz property of the Hessian matrix).**

Suppose Assumption 1 holds. Let $\mathbf{x}_{i},\mathbf{x}_{i}^{\prime}\in\mathbb{R}^{d}$, then there exists a constant $N_{c,\delta,\lambda}\in\mathbb{N}_{+}$ depending on $c,\delta,\lambda$, such that for all $n\geq N_{c,\delta,\lambda}$,

[TABLE]

with probability at least $1-n^{-c}$.

Proof.

Write

[TABLE]

which can be viewed as a sum of two terms. For the first term,

[TABLE]

with probability at least $1-n^{-c}$. For the second term,

[TABLE]

with probability at least $1-n^{-c}$. So we have

[TABLE]

with probability at least $1-n^{-c}$. ∎

Theorem A.7** (One-step estimator).**

Suppose Assumption 1 holds. Let $\breve{\mathbf{X}}$ denote the adjacency spectral embedding, and $\widetilde{\mathbf{X}}$ the signature-adjusted adjacency spectral embedding. Let $p_{0ij}=\rho_{n}\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}$, and $\widetilde{p}_{ij}=\breve{\mathbf{x}}_{i}^{\mathrm{T}}\widetilde{\mathbf{x}}_{j}$, $i,j\in[n]$. For each $i\in[n]$, define the one-step estimator $\widehat{\mathbf{x}}_{i}^{\mathrm{(OS)}}$ by

[TABLE]

Then

[TABLE]

and for any $c>0$, there exist a constant integer $N_{c,\delta,\lambda}$ and a constant $C_{c,\delta,\lambda}$ that depend on $c,\delta,\lambda$ such that for all $1\leq t\leq C_{c,\delta,\lambda}\log n$ and for all $n\geq N_{c,\delta,\lambda}$, $\|\mathbf{r}_{in}^{\mathrm{(OS)}}\|_{2}\lesssim_{c,\delta,\lambda}{t^{2}}/{(n\rho_{n}^{1/2})}$ with probability at least $1-c_{0}n^{-c}-c_{0}e^{-t}$ for some absolute constant $c_{0}>0$. Furthermore, $\sqrt{n}\mathbf{G}_{0in}^{1/2}(\mathbf{W}^{\mathrm{T}}\widehat{\mathbf{x}}_{i}^{\mathrm{(OS)}}-\rho_{n}^{1/2}\mathbf{x}_{0i})\overset{\mathcal{L}}{\to}\mathrm{N}_{d}({\bm{0}}_{d},{\mathbf{I}}_{d})$ as $n\to\infty$.

Remark 3*.*

Theorem A.7 is a generalization of Theorem 4.7 in Xie, (2024) to generalized random dot product graphs in our settings. The proof is mostly identical to its original version, with slight modifications such as the presence of the signature matrix ${\mathbf{I}}_{p,q}$ and the different definition of the orthogonal alignment matrix $\mathbf{W}$. We omit the proof and use the theorem directly.

Lemma A.8** (Theorem 1 in de la Pena and Montgomery-Smith, (1995)).**

Let $\{X_{i}\}$ be a sequence of independent random variables on a measurable $(S,\mathscr{S})$ space and let $\{X_{i}^{(1)}\}$, $\{X_{i}^{(2)}\}$ be two independent copies of $\{X_{i}\}$. Let $f_{i_{1}i_{2}}$ be families of functions of two variables taking $(S\times S)$ into a Banach space $(B,\|\cdot\|)$. Then, for all $n\geq 2$, $t>0$, there exist a numerical constant $C$ such that

[TABLE]

Lemma A.9** (A weak law of large numbers).**

Suppose Assumption 1 holds. Let

[TABLE]

then $Z=\mathbb{E}[Z]+o_{{\mathbb{P}}}(1)$, where $\mathbb{E}[Z]=(1/n)\sum_{i=1}^{n}\mathrm{tr}(\mathbf{G}_{0in}^{-1})$.

Proof.

Recall $p_{0ij}=\rho_{n}\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}$ and $\mathbf{G}_{0in}=(1/n)\sum_{j=1}^{n}{\rho_{n}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}\mathbf{x}_{0j}^{\mathrm{T}}{\mathbf{I}}_{p,q}}/\{p_{0ij}(1-p_{0ij})\}$. For $i,j\in[n]$, let $E_{ij}=A_{ij}-p_{0ij}$, and ${\bm{\gamma}}_{ij}={\mathbf{G}_{0in}^{-1}\rho_{n}^{1/2}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}}/\{np_{0ij}(1-p_{0ij})\}$, then we can write

[TABLE]

We have $\mathbb{E}[E_{ij}]=0$, $\mathbb{E}[E_{ij}^{2}]=p_{0ij}(1-p_{0ij})$, and $\|{\bm{\gamma}}_{ij}\|_{2}\asymp_{\delta,\lambda}1/(n\rho_{n}^{0.5})$ by assumption and the result that $\|\mathbf{G}_{0in}^{-1}\|_{2}\asymp_{\delta,\lambda}1$ which is shown in (B.2) on page B.2 (in the proof of asymptotic normality part of Theorem 3.2). By Bernstein’s inequality,

[TABLE]

from which by taking $t=C_{c,\delta,\lambda}\sqrt{(\log n)/n^{2}\rho_{n}}$ we have

[TABLE]

with probability at least $1-n^{-c}$. We then use Lemma A.8 to deal with the sum of cross terms. Consider the sequence of random variables $\{E_{(i,j)}:(i,j)\in[n]^{2}\}$, with two independent copies $\{E_{ij}\}$ and $\{\bar{E}_{ij}\}$, and equip the index set $[n]^{2}$ with the lexicographic order, i.e., for $x=(x_{1},x_{2})$ and $y=(y_{1},y_{2})$ in $[n]^{2}$, we have $x<y$ if either $x_{1}<y_{1}$ or $x_{1}=y_{1}$ and $x_{2}<y_{2}$, and we have $x=y$ if $x_{1}=y_{1}$ and $x_{2}=y_{2}$. And consider the family of functions $f_{(i_{1},a),(i_{2},b)}(E_{(i_{1},a)},E_{(i_{2},b)})=E_{(i_{1},a)}E_{(i_{2},b)}{\bm{\gamma}}_{(i_{1},a)}^{\mathrm{T}}{\bm{\gamma}}_{(i_{2},b)}{\mathbbm{1}}(i_{1}=i_{2})$. Then

[TABLE]

By Lemma A.8, conditional probability, and Bernstein’s inequality,

[TABLE]

from which by taking $t=C_{c,\delta,\lambda}\sqrt{(\log n)/n\rho_{n}}$ we have

[TABLE]

with probability at least $1-n^{-c}$. Note that $\mathbb{E}[Z]=\mathbb{E}[\sum_{i=1}^{n}\sum_{j=1}^{n}E_{ij}^{2}\|{\bm{\gamma}}_{ij}\|_{2}^{2}]$ by independence of $E_{ij}$. So by combining several previous results, we have

[TABLE]

Finally, a simple algebra shows that

[TABLE]

∎

Appendix B Proofs of the Main Results

B.1 Proof of Theorem 3.2

Proof.

For simplicity of notation, in the proof of this theorem, the large probability bounds with probability at least $1-n^{-c}$ are stated with respect to all $n\geq{}{N}_{c,\delta,\lambda}$ for some large constant integer ${N}_{c,\delta,\lambda}$ that depends on $c,\delta,\lambda$, where $c>0$ is an arbitrary positive constant. Also, the results that hold for a single $i\in[n]$ with probability at least $1-n^{-c}$ can be strengthened to hold for all $i\in[n]$ by taking a union bound over $i\in[n]$.

Proof of existence and uniqueness. Note that we have the average ESL function for vertex $i$

[TABLE]

and define the population counterpart of the average of the ESL function

[TABLE]

Note that $\psi_{n}^{\prime\prime}(t)\in[-\tau_{n}^{-2},-1]$ for all $t\in\mathbb{R}$. Therefore,

[TABLE]

in which by Theorem 5.2 in Lei and Rinaldo, (2015) and Weyl’s inequality,

[TABLE]

with probability at least $1-n^{-c}$, so we have

[TABLE]

for all $\mathbf{x}_{i}\in\mathbb{R}^{d}$ with probability at least $1-n^{-c}$, i.e., $(1/n)\widehat{\ell}_{in}(\mathbf{x}_{i})$ is strongly concave over $\mathbb{R}^{d}$ with probability at least $1-n^{-c}$. This implies that $\operatorname*{arg\,max}_{\mathbf{x}_{i}\in\mathbb{R}^{d}}\widehat{\ell}_{in}(\mathbf{x}_{i})$ exists and is unique because $\widehat{\ell}_{in}(\mathbf{x}_{i})$ is clearly bounded from above.

Next, we let $\widehat{\mathbf{x}}_{i}=\operatorname*{arg\,max}_{\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}}\widehat{\ell}_{in}(\mathbf{x}_{i})$ be the local maximizer of $\widehat{\ell}_{in}(\mathbf{x}_{i})$ in the closed ball $\{\mathbf{x}_{i}\in\mathbb{R}^{d}:\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}\}$. In $\{\mathbf{x}_{i}:\,\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}\}$,

[TABLE]

It is easy to see that $\widehat{\ell}_{in}(\mathbf{x}_{i})$ is continuous over $\mathbb{R}^{d}$, so there exists a maximizer of $\widehat{\ell}_{in}(\mathbf{x}_{i})$ in $\{\mathbf{x}_{i}:\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}\}$ which is a compact set. Let $\widehat{\mathbf{x}}_{i}$ denote this maximizer of $\widehat{\ell}_{in}(\mathbf{x}_{i})$ in $\{\mathbf{x}_{i}:\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}\}$, and over this set, by Taylor’s theorem,

[TABLE]

and also by Taylor’s theorem and by Lemma A.4

[TABLE]

for all $i\in[n]$ with probability at least $1-n^{-c}$. Combining the two inequalities above, we have

[TABLE]

with probability at least $1-n^{-c}$. Then by taking $N_{c,\delta,\lambda}$ large enough such that

$C_{c,\delta,\lambda}((\log n)/(n\rho_{n}))^{1/4}<1-\delta/2-\sqrt{1-\delta}$ for all $n\geq N_{c,\delta,\lambda}$, we have

[TABLE]

which implies that $\widehat{\mathbf{x}}_{i}$ lies in the interior of the ball $\{\mathbf{x}_{i}:\|\mathbf{x}_{i}\|_{2}\leq\rho_{n}^{1/2}\}$, i.e., $\widehat{\mathbf{x}}_{i}$ is a local maximizer, with probability at least $1-n^{-c}$. Since $\widehat{\ell}_{in}(\mathbf{x}_{i})$ is strongly concave over $\mathbb{R}^{d}$ with probability at least $1-n^{-c}$, the local maximizer $\widehat{\mathbf{x}}_{i}$ is the unique global maximizer with probability at least $1-n^{-c}$. This implies that, for all $i\in[n]$ with probability at least $1-n^{-c}$, $\|\widehat{\mathbf{x}}_{i}\|_{2}<(1-\delta/2)\rho_{n}^{1/2}$.

Proof of Consistency. Over this event that happens with probability at least $1-n^{-c}$, we have $\|\widehat{\mathbf{x}}_{i}\|_{2}<(1-\delta/2)\rho_{n}^{1/2}$, and $\widehat{\mathbf{x}}_{i}$ is the unique global maximizer and thus $(\partial\widehat{\ell}_{in})/(\partial\mathbf{x}_{i})(\widehat{\mathbf{x}}_{i})={\bm{0}}_{d}$. Now, by mean value theorem,

[TABLE]

where $\bar{\mathbf{x}}_{i}=\theta\mathbf{W}^{\mathrm{T}}\widehat{\mathbf{x}}_{i}+(1-\theta)\rho_{n}^{1/2}\mathbf{x}_{0i}$ for some $\theta\in(0,1)$, then it follows that

[TABLE]

and by Lemma A.4,

[TABLE]

with probability at least $1-n^{-c}$. Combining the two inequalities above, we have

[TABLE]

with probability at least $1-n^{-c}$. The consistency is thus proved.

Proof of Asymptotic normality. Let $\breve{\mathbf{x}}_{i}$ denote the adjacency spectral embedding, $i\in[n]$, and let $\widetilde{p}_{ij}=\breve{\mathbf{x}}_{i}^{\mathrm{T}}\widetilde{\mathbf{x}}_{j}$, $p_{0ij}=\rho_{n}\mathbf{x}_{0i}^{\mathrm{T}}{\mathbf{I}}_{p,q}\mathbf{x}_{0j}$, $i,j\in[n]$. Define the and its plug-in estimate of $\mathbf{G}_{0in}$ with the adjacency spectral embedding by $\widetilde{\mathbf{G}}_{in}=(1/n)\sum_{j=1}^{n}{\widetilde{\mathbf{x}}_{j}\widetilde{\mathbf{x}}_{j}^{\mathrm{T}}}/\{\widetilde{p}_{ij}(1-\widetilde{p}_{ij})\}$. It is not hard to see that

[TABLE]

i.e., $\mathbf{G}_{0in}$ is a positive definite matrix with eigenvalues bounded away from [math] and $\infty$, and by Lemma A.2, Theorem 5.2 in Lei and Rinaldo, (2015), and Weyl’s inequality, we also have $\lambda_{1}(\widetilde{\mathbf{G}}_{in})\leq 2/\delta$ and $\lambda_{d}(\widetilde{\mathbf{G}}_{in})\geq\lambda/\{2(1-\delta/2)\}$, i.e., $\widetilde{\mathbf{G}}_{in}$ is a positive definite matrix with eigenvalues bounded away from [math] and $\infty$ with probability at least $1-n^{-c}$. Now we have

[TABLE]

with probability at least $1-n^{-c}$, where the inequalities follows from triangle inequality, mean value theorem, Cauchy-Schwarz inequality, and Lemma A.2. And we have

[TABLE]

with probability at least $1-n^{-c}$, by Theorem A.1 (bound for $\breve{\mathbf{x}}_{i}$) and Lemma A.5. So, by triangle inequality and the previous two bounds, we have

[TABLE]

with probability at least $1-n^{-c}$.

Now, we show the asymptotic normality of the maximizer of the ESL function by showing that $\widehat{\mathbf{x}}_{i}$ and $\widehat{\mathbf{x}}_{i}^{\mathrm{(OS)}}$ are close enough and then applying Slutsky’s theorem to utilize the asymptotic normality of the one-step estimator $\widehat{\mathbf{x}}_{i}^{\mathrm{(OS)}}$. For each $k\in[d]$, apply Taylor’s theorem to $(1/n)(\partial\ell_{in})/(\partial\mathbf{x}_{i})(\widehat{\mathbf{x}}_{i})=0$ at $\mathbf{x}_{i}=\breve{\mathbf{x}}_{i}$ to obtain

[TABLE]

where $\bar{\mathbf{x}}_{i}=\theta\widehat{\mathbf{x}}_{i}+(1-\theta)\breve{\mathbf{x}}_{i}$ for some $\theta\in(0,1)$. By triangle inequality, Theorem A.1, and the consistency result that has been shown above, $\|\mathbf{W}\bar{\mathbf{x}}_{i}-\rho_{n}^{1/2}\mathbf{x}_{0i}\|_{2}\lesssim_{c,\delta,\lambda}\sqrt{(\log n)/n}$ with probability at lease $1-n^{c}$. It is easy to see that, over such an event,

[TABLE]

Then, by Lemma A.2, Lemma A.3, and Cauchy-Schwarz inequality,

[TABLE]

with probability at lease $1-n^{c}$, and also by triangle inequality, Theorem A.1, and the consistency result that has been shown above,

[TABLE]

with probability at lease $1-n^{c}$. So the Taylor expansion mentioned above can be written as

[TABLE]

where $\mathbf{R}_{in1}$ is a random matrix with $\|\mathbf{R}_{in1}\|_{2}\lesssim_{c,\delta,\lambda}\sqrt{(\log n)/(n\rho_{n})}$ with probability at lease $1-n^{c}$. By the approximation of $(1/n)(\partial^{2}\widehat{\ell}_{in})/(\partial\mathbf{x}_{i}\partial\mathbf{x}_{i}^{\mathrm{T}})(\breve{\mathbf{x}}_{i})$ and $\widetilde{\mathbf{G}}_{in}$ that has been shown above, Lemma A.2 (for the large probability bounds on $\widetilde{p}_{ij}$), and the definition of $(1/n)(\partial\widehat{\ell}_{in})/(\partial\mathbf{x}_{i})(\breve{\mathbf{x}}_{i})$, we have

[TABLE]

where $\mathbf{R}_{in2}$ is a random matrix with $\|\mathbf{R}_{in2}\|_{2}\lesssim_{c,\delta,\lambda}\sqrt{(\log n)/(n\rho_{n})}$ with probability at lease $1-n^{c}$. Now we have $\|\widetilde{\mathbf{G}}_{in}^{-1}\mathbf{R}_{in2}\|_{2}\lesssim_{c,\delta,\lambda}\sqrt{(\log n)/(n\rho_{n})}$ with probability at lease $1-n^{c}$. Write

[TABLE]

then write

[TABLE]

and under the assumption that $(\log n)/(n\rho_{n})\to 0$, by Theorem A.1 and Theorem A.7 (take $t=C_{c,\delta,\lambda}\log(n\rho_{n})$), we have

[TABLE]

with probability at least $1-(n\rho_{n})^{-c}$, so $\|\widehat{\mathbf{x}}_{i}-\widehat{\mathbf{x}}_{i}^{\mathrm{(OS)}}\|_{2}\lesssim_{c,\delta,\lambda}\log n/(n\rho_{n}^{1/2})$ with probability at least $1-(n\rho_{n})^{-c}$. By Theorem A.7 and Slutsky’s theorem, we have

[TABLE]

where

[TABLE]

with probability at least $1-(n\rho_{n})^{-c}$, and $\sqrt{n}\mathbf{G}_{0in}^{1/2}(\mathbf{W}^{\mathrm{T}}\widehat{\mathbf{x}}_{i}-\rho_{n}^{1/2}\mathbf{x}_{0i})\overset{\mathcal{L}}{\to}\mathrm{N}_{d}({\bm{0}}_{d},\,{\mathbf{I}}_{d})$.

Proof of consistency of global error. Under the assumption that $(\log n)^{4}/(n\rho_{n})\to 0$, by Theorem A.1 and Theorem A.7 (take $t=C_{c,\delta,\lambda}\log n$ in Theorem A.7), we have

[TABLE]

with probability at least $1-n^{-c}$, so $\|\widehat{\mathbf{x}}_{i}-\widehat{\mathbf{x}}_{i}^{\mathrm{(OS)}}\|_{2}\lesssim_{c,\delta,\lambda}{\log n}/{(n\rho_{n}^{1/2})}$ with probability at least $1-n^{-c}$, then by Theorem A.7 and Slutsky’s theorem, we have

[TABLE]

where $\|\mathbf{r}_{in}\|_{2}\lesssim_{c,\delta,\lambda}{(\log n)^{2}}/{(n\rho_{n}^{1/2})}$ with probability at least $1-n^{-c}$. Now write

[TABLE]

which can be viewed as a sum of three terms. For the first term, by Lemma A.9,

[TABLE]

For the second term, recalling that $\mathbf{G}_{0in}$ is a positive definite matrix with eigenvalues bounded away from [math] and $\infty$, we have

[TABLE]

with probability at least $1-n^{-c}$. For the third term, by triangle inequality and Cauchy-Schwarz inequality (twice),

[TABLE]

Hence, we conclude that $\|\widehat{\mathbf{X}}\mathbf{W}-\rho_{n}^{1/2}\mathbf{X}_{0}\|_{\mathrm{F}}^{2}=(1/n)\sum_{i=1}^{n}\mathrm{tr}(\mathbf{G}_{0in}^{-1})+o_{{\mathbb{P}}}(1)$. ∎

B.2 Proof of Theorem 3.3

Proof.

It is sufficient to prove that

[TABLE]

with probability at least $1-n^{-c}$. Recall that the posterior density associated with the ESL function for vertex $i$ is $\pi_{in}(\mathbf{x}_{i}\mid\mathbf{A})=e^{\widehat{\ell}_{in}(\mathbf{x}_{i})-\widehat{\ell}_{in}(\widehat{\mathbf{x}}_{i})}\pi_{i}(\mathbf{x}_{i})/d_{in}$, where $d_{in}=\int_{\mathbb{R}^{d}}\exp\{\widehat{\ell}_{in}(\mathbf{x}_{i})-\widehat{\ell}_{in}(\widehat{\mathbf{x}}_{i})\}\pi_{i}(\mathbf{x}_{i})\mathrm{d}\mathbf{x}_{i}$ is the normalizing constant. Let $\eta_{n}=C_{c,\delta,\lambda}\sqrt{\log n}$, and partition $\mathbb{R}^{d}$ as $\mathcal{A}_{1in}\cup\mathcal{A}_{2in}$ where

[TABLE]

For the integral over $\mathcal{A}_{2in}$,

[TABLE]

of which for the first term, by (B.1) on page B.1 and the assumption that $\rho_{n}=1$, we have

[TABLE]

with probability at least $1-n^{-c}$, and for the second term, by (B.2) on page B.2, we have

[TABLE]

with probability at least $1-n^{-c}$. Over $\mathcal{A}_{1in}$,

[TABLE]

with probability at least $1-n^{-c}$ by Theorem 3.2 and the definition of $\eta_{n}$. For the integral over $\mathcal{A}_{1in}$, with the change of variable $\mathbf{t}_{i}=\sqrt{n}\mathbf{W}^{\mathrm{T}}(\mathbf{x}_{i}-\widehat{\mathbf{x}}_{i})$,

[TABLE]

which is a sum of two terms, in which $\bar{\mathbf{x}}_{i}$ is a convex combination of $\mathbf{x}_{i}$ and $\widehat{\mathbf{x}}_{i}$. For the first term of this integral over $\mathcal{A}_{1in}$ we have

[TABLE]

with probability at least $1-n^{-c}$ by Lemma A.5, and for the second term of this integral over $\mathcal{A}_{1in}$ we have

[TABLE]

with probability at least $1-n^{-c}$ by Assumption 2.

Hence, the integral in the statement of this theorem is bounded by

[TABLE]

with probability at least $1-n^{-c}$. ∎

B.3 Proof of Theorem 3.4

The proof of Theorem 3.4 can be done by the triangle inequality, Pinsker’s inequality, and Lemma B.1 below.

Lemma B.1**.**

*Suppose the conditions in Theorem 3.3 hold. Then *

[TABLE]

Proof.

Note that Theorem 3.3 implies that

[TABLE]

with probability at least $1-n^{-c}$.

[TABLE]

where $\mathbf{x}_{i}^{\prime}$ and $\mathbf{x}_{i}^{\prime\prime}$ lie between $\mathbf{x}_{i}$ and $\rho_{n}^{1/2}\mathbf{W}\mathbf{x}_{0i}$, and $\theta_{in}=O(1)$ with probability at least $1-n^{-c}$ by Theorem 3.3. Let $\eta_{n}=C_{c,\delta,\lambda}\sqrt{\log n}$, and partition $\mathbb{R}^{d}$ as $\mathcal{A}_{1in}\cup\mathcal{A}_{2in}$ where

[TABLE]

On $\mathcal{A}_{1in}$,

[TABLE]

with probability at least $1-n^{-c}$ by Theorem 3.2, and the same bound also holds for $\|\mathbf{W}^{\mathrm{T}}\mathbf{x}_{i}^{\prime}-\rho_{n}^{1/2}\mathbf{x}_{0i}\|_{2}$ and $\|\mathbf{W}^{\mathrm{T}}\mathbf{x}_{i}^{\prime\prime}-\rho_{n}^{1/2}\mathbf{x}_{0i}\|_{2}$. Then by Lemma A.5, Taylor’s theorem, and Assumption 2,

[TABLE]

with probability at least $1-n^{-c}$. So the integral over $\mathcal{A}_{1in}$

[TABLE]

with probability at least $1-n^{-c}$. We then consider the integral over $\mathcal{A}_{2in}$. Note that over $\mathcal{A}_{2in}$, by Theorem 3.2,

[TABLE]

with probability at least $1-n^{-c}$. Then by triangle inequality, Lemma A.5, and Lemma A.6, over $\mathcal{A}_{2in}$ we have

[TABLE]

with probability at least $1-n^{-c}$. And by Assumption 2 we have

[TABLE]

with probability at least $1-n^{-c}$. So the integral over $\mathcal{A}_{2in}$

[TABLE]

with probability at least $1-n^{-c}$. This completes the proof. ∎

B.4 Proof of Theorem 3.5

Proof.

Let $Q_{in}^{*}$ denote the variational posterior distribution $N(\mathbf{x}_{i}^{*},{\bm{\Sigma}}_{in}^{*})$, with density $q_{in}^{*}(\mathbf{x}_{i})$, and let $N_{in}^{*}$ denote the normal distribution $N(\widehat{\mathbf{x}}_{i},\{n\mathbf{W}\mathbf{G}_{0in}\mathbf{W}^{\mathrm{T}}\}^{-1})$, with density $\phi_{d}(\mathbf{x}_{i}|\widehat{\mathbf{x}}_{i},(n\mathbf{W}\mathbf{G}_{0in}\mathbf{W}^{\mathrm{T}})^{-1})$. Note that Theorem 3.3 implies that

[TABLE]

with probability at least $1-n^{-c}$, and Theorem 3.4 implies that

[TABLE]

with probability at least $1-n^{-c}$ by triangle inequality.

Let $\varphi_{Q}(\mathbf{t})$ denote the characteristic function of a distribution $Q$. We have

[TABLE]

with probability at least $1-n^{-c}$. Recall that both $Q_{in}^{*}$ and $N_{in}^{*}$ are $d$-dimensional normal distributions, so

[TABLE]

By triangle inequality and the fact that $|\exp\{ix\}|=1$ for all $x\in\mathbb{R}$, we have

[TABLE]

with probability at least $1-n^{-c}$. Also note that, for $\mathbf{u}\in\mathbb{R}^{d}$ with $\|\mathbf{u}\|_{2}=1$,

[TABLE]

since $\mathbf{G}_{0in}$ has all eigenvalues positive and bounded away from [math] and $+\infty$ by (B.2) on page B.2. Then,

[TABLE]

which, by taking $\mathbf{t}=\sqrt{n}\mathbf{u}$, implies that

[TABLE]

for all $\|\mathbf{u}\|_{2}=1$, $\mathbf{u}\in\mathbb{R}^{d}$, with probability at least $1-n^{-c}$. Namely, $\sqrt{n}(\widehat{\mathbf{x}}_{i}-\mathbf{x}_{i}^{*})=o_{{\mathbb{P}}}(1)$. By Theorem 3.2 and Slutsky’s theorem, $\sqrt{n}\mathbf{G}_{0in}^{1/2}(\mathbf{W}^{\mathrm{T}}\mathbf{x}_{i}^{*}-\rho_{n}^{1/2}\mathbf{x}_{0i})\overset{\mathcal{L}}{\to}\mathrm{N}_{d}({\bm{0}}_{d},\,{\mathbf{I}}_{d})$. ∎

B.5 Proof of Theorem 3.1

Proof.

We borrow the idea in the proof of Theorem 11 in Xu and Campbell, (2023) to proof the strong convexity. By Assumption 2, $-\log\pi_{i}(\mathbf{x}_{i})$ is convex in $\mathbf{x}_{i}$. By Exercise 12.21 in Abadir and Magnus, (2005), $-\log\det(\mathbf{L}_{i})$ is convex in $\mathbf{L}_{i}$. By (B.1) in the proof of Theorem 3.2 on page B.1, $-\widehat{\ell}(\mathbf{x}_{i})$ is strongly convex in $\mathbf{x}_{i}$ with strong convexity parameter $n\lambda\rho_{n}$ with probability at least $1-n^{-c}$. Let $\mathbf{D}_{n}=n\lambda\rho_{n}{\mathbf{I}}_{d}$, and note that

[TABLE]

Then by the linearity of expectation and the strong convexity of $-\widehat{\ell}_{in}(\mathbf{x}_{i})$, the function

[TABLE]

is convex in $({\bm{\mu}}_{i},\mathbf{L}_{i})\in\mathbb{R}^{d}\times\mathcal{L}_{d\times d}$ with probability at least $1-n^{-c}$. So the function

[TABLE]

is convex in $({\bm{\mu}}_{i},\mathbf{L}_{i})\in\mathbb{R}^{d}\times\mathcal{L}_{d\times d}$, with probability at least $1-n^{-c}$, which equivalently means that the function $F_{in}({\bm{\mu}}_{i},\mathbf{L}_{i})$ is strongly convex in $({\bm{\mu}}_{i},\mathbf{L}_{i})\in\mathbb{R}^{d}\times\mathcal{L}_{d\times d}$, with probability at least $1-n^{-c}$.

We next prove the interchange of derivative and expectation. By the definition of $\widehat{\ell}(\mathbf{x}_{i})$,

[TABLE]

which is polynomial in $\mathbf{x}_{i}$, so $\|\nabla_{{\bm{\mu}}_{i},\mathbf{L}_{i}}\widehat{\ell}_{in}\left({\bm{\mu}}_{i}+\frac{1}{\sqrt{n}}\mathbf{L}_{i}\mathbf{Z}\right)\|_{2}$ is dominated by a function that is integrable with respect to the standard normal measure. By Assumption 2,

[TABLE]

which is polynomial in $\mathbf{x}_{i}$, so $\|\nabla_{{\bm{\mu}}_{i},\mathbf{L}_{i}}\log\pi_{i}\left({\bm{\mu}}_{i}+\frac{1}{\sqrt{n}}\mathbf{L}_{i}\mathbf{Z}\right)\|_{2}$ is dominated by a function that is integrable with respect to the standard normal measure. Hence,

[TABLE]

So $\nabla_{{\bm{\mu}}_{i},\mathbf{L}_{i}}F_{in}({\bm{\mu}}_{i},\mathbf{L}_{i})=\mathbb{E}_{\mathbf{Z}}\left[\nabla_{{\bm{\mu}}_{i},\mathbf{L}_{i}}f_{in}({\bm{\mu}}_{i},\mathbf{L}_{i})\right]$. ∎