Network Representation Using Graph Root Distributions
Jing Lei



§ INTRODUCTION

In recent years network analysis has been the focus of many theoretical and applied research efforts in the scientific community, due to the increasing popularity of relational data.  Generally speaking, a network records the presence and absence of pairwise interactions among a group of individuals, and statistical network analysis aims at recovering properties of the underlying population of individuals from their pairwise interactions.  There is a vast literature on network analysis, and we refer to <cit.> for more detailed review of this field from a statistical perspective.

Exchangeable random graphs <cit.> are an important class of probabilistic models for network data.  The exchangeability requirement is quite natural: The individuals recorded in the network are somewhat like random sample points, and the data distribution remains unchanged under permutation of the nodes.  Many popularly studied network models are special cases of exchangeable random graphs, including the stochastic block model <cit.>, the degree-corrected block model <cit.>, the mixed-membership block model <cit.>, the random dot-product graph model <cit.>, and random geometric graphs <cit.>.

A central piece of the theoretical foundation of exchangeable random graphs is the celebrated Aldous-Hoover Theorem <cit.>, which says that any exchangeable random graph of infinite size can be generated  by first sampling independent node variables (s_i:i≥ 1) uniformly on [0,1], and then connect each pair of nodes (i,j) independently with probability W(s_i,s_j), for some symmetric W:[0,1]^2↦ [0,1].  The representation theory also says that two functions W_1, W_2 lead to the same distribution of exchangeable random graphs if and only if there exist measure-preserving mappings h_1, h_2, both from [0,1] to [0,1], such that W_1(h_1(s),h_1(t))=W_2(h_2(s),h_2(t)) almost everywhere over (s,t)∈ [0,1]^2.  
Other random graph models have been proposed, such as exchangeable random measures <cit.>, edge exchangeability <cit.>, and ideas using a Bayesian framework <cit.>. A notable difference is that these models can cover sparse networks while the Aldous-Hoover type exchangeable arrays must be dense.  Nevertheless, when a finite sample is concerned, one can add a sparsity parameter to generate a sparse network using the Aldous-Hoover type exchangeable array.  See <Ref> for further discussion on sparse networks.  

The purpose of this work is to develop an alternative parametrization for a subclass of exchangeable random graphs and corresponding estimators. This new parameterization is based on a representation of nodes as independent random vectors in a separable space 𝒦, and the edge probability is the inner product of the two node vectors. A space is isomorphic to a Hilbert space, but has an indefinite inner product.  With such a space node embedding, we shift the information hidden in the function W to the probability measure on 𝒦, which can exhibit salient structures in a transparent and geometric manner.
  

We highlight a few key contributions.

  
  * We provide a constructive proof for the correspondence between a subclass of exchangeable random graphs and probability distributions on the space 𝒦.  Our construction starts from viewing W as an integral operator and considering its spectral decomposition. The variable s∼ Unif(0,1) is treated as the input variable in an infinite dimensional inverse transform sampling. The induced measure is thus invariant under measure-preserving transforms of s. The existence and identifiability of such a representation is established for a wide class of exchangeable random graphs. In our construction, it becomes apparent that the induced measure is closely related to the square root of the integral operator W, with appropriate treatment of negative eigenvalues.  Thus we call this induced measure the graph root distribution (GRD).
  
  * We show that the Wasserstein distance between two graph root distributions on 𝒦 provides an upper bound of the cut-distance between sampling distributions of the corresponding exchangeable random graphs.  This result is further extended to a modified version of Wasserstein distance that is suitable for measuring the distance between two equivalence classes of graph root distributions.
  
  * We show that a truncated adjacency spectral embedding weighted by the square roots of the absolute eigenvalues can approximate the empirical distribution of the latent node vectors in 𝒦 with vanishing Wasserstein distance error when the network size n goes to infinity, under suitable regularity conditions.  This in turn implies that such a weighted truncated spectral embedding can also consistently estimate the underlying graph root distribution when the truncation dimension is chosen appropriately.
  
  * As demonstrated in our numerical examples, the new parameterization allows for a simple, theoretically justifiable estimator, which can reveal salient geometric features in the network data. 




 §.§ Related work.

The GRD parameterization is closely related to latent space network models.  The idea of modeling the edge probability between a pair of nodes by the inner product of the corresponding latent vectors is studied by <cit.> under the name of “latent eigenmodel”, and further developed by <cit.> under the name of “random dot-product graph”.   The GRD framework extends this idea to a population perspective and connects it with the graphon literature.

The GRD parameterization and estimation are also related to spectral methods in random graphs. The spectral approach has been used in graphon estimation by <cit.>.  The GRD estimation method considered in <Ref> uses a similar singular value thresholding of the adjacency matrix, but takes a further step of weighted spectral embedding to recover the latent vectors corresponding to each node, with the target parameter being a distribution in the latent space.  Besides spectral methods, graphon estimation has also been studied using histogram and smoothing methods, including theoretical analysis <cit.>, and practical algorithms <cit.>. GRD estimation and graphon (or probability matrix) estimation are different, as the target parameters are in different spaces with different error metrics (see <Ref> for further discussion and comparison).  Some potential benefits of using GRD are discussed in <Ref>, and visualized in some simulated and real data examples in <Ref>.   



§ BACKGROUND



 §.§ Exchangeable random graphs and the graphon parameterization


Consider a random symmetric two-way binary array

    𝐀 = (A_ij: 1≤ i<j)

such that A_ii=0 and satisfies the row-column joint exchangeability

    (A_ij: i≥ 1, j≥ 1) d= (A_σ(i)σ(j): i≥ 1,  j≥ 1)

for all finite index permutation mapping σ: for some 1≤ i_0 < j_0,

    σ(i) = {[                i if  i∉{i_0,j_0},;              j_0     if  i = i_0,;              i_0    if  i = j_0 . ].

Here “d=” means that two random objects have the same distribution.

Analogous to the de Finetti theorem, the Aldous-Hoover theorem <cit.> says that any symmetric exchangeable binary array 𝐀 of infinite size can be generated by sampling independent (s_i:i≥ 1) from Unif(0,1) (the uniform distribution on [0,1]), and sampling A_ij independently from a Bernoulli distribution with probability W(s_i,s_j) for a symmetric measurable function W(·,·):[0,1]^2↦ [0,1].  Here W is a random object measurable in the doubly-exchangeable σ-field.  For any given realization of 𝐀, we can simply treat W as a non-random parameter.




Once W is given, the distribution of 𝐀 is completely determined.  However, the converse is not true.  Let h:[0,1]↦ [0,1] be a measure-preserving mapping in the sense that

    μ(h^-1(B)) = μ(B) ,  ∀ B∈ℬ_[0,1] ,

where μ(·) denotes the Lebesgue measure and ℬ_[0,1] is the Borel σ-field.  




Two functions W_1 and W_2 generate the same distribution of exchangeable arrays if and only if there exist two measure-preserving mappings h_1,h_2 such that

    W_1(h_1(s),h_1(s'))=W_2(h_2(s),h_2(s')) ,  a.e. .

When (<ref>) holds, we say W_1 and W_2 are weakly isomorphic, denoted as

    W_1w.i.=W_2 .

The notion “w.i.=” defines an equivalence relation on the space of all symmetric functions that map [0,1]^2 to [0,1]. We use W̃ to denote the equivalence class containing W.

When W_1 and W_2 are not weakly isomorphic, then they lead to different distributions of exchangeable random graphs. In this case, the sub-graph counts have different distributions under W_1 and W_2. Such a sampling distribution difference can be linked to the cut-distance, defined as

    δ_□(W_1,W_2)
    
    =   inf_h_1,h_2sup_S× S':⊆[0,1]^2|∫_S× S'[W_1(h_1(s),h_1(s'))-W_2(h_2(s),h_2(s'))]dsds'| ,

where h_1, h_2 range over all measure-preserving mappings.
It can be shown that W_1w.i.=W_2 if and only if δ_□(W_1, W_2)=0.

Therefore, the cut-distance δ_□(·,·) can also be used to measure the distance between two equivalence classes W̃_1 and W̃_2.  

We will also adopt the terminology used in <cit.> to call the function W a graphon (abbreviation for “graph function”).  In the rest of this paper we will often use a graphon W to represent its equivalence class. 








 §.§ Graph root distributions

The focus of this paper is to develop an alternative characterization of a subclass of exchangeable random graphs.  The construction involves probability distributions on a separable space, which we introduce first.


  A space 𝒦=ℋ_+⊖ℋ_- is the direct sum of two Hilbert spaces ℋ_+ and ℋ_-. For each (x,y), (x',y')∈𝒦 with x,x'∈ℋ_+ and y,y'∈ℋ_-, the inner product is
  
    ⟨ (x,y), (x',y') ⟩_𝒦=
        ⟨ x,x'⟩_ℋ_+ - ⟨ y,y' ⟩_ℋ_- .

  The space 𝒦 is also a linear normed space isomorphic to the Hilbert space ℋ_+⊕ℋ_- equipped with norm
  
    (x,y)_𝒦=(x_ℋ_+^2+y_ℋ_-^2)^1/2 .


The notation “⊖” is used to emphasize the non-positive definite inner product associated with the space 𝒦, and
ℋ_+, ℋ_- represent the subspaces containing the positive and negative components, respectively.  This notation has been used in existing machine learning
 literature involving spaces <cit.>.  The traditional notation “⊕” is saved for the direct sum of Hilbert spaces in the usual sense, where the direct sum is still a Hilbert space with a positive definite inner product.


Now we define graph root distributions.

  We call a probability measure F on 𝒦 a graph root distribution if for two independent samples  Z_1 and Z_2 from F
  
    ℙ(⟨ Z_1,Z_2 ⟩_𝒦∈ [0,1])=1 .


Let F be a GRD on 𝒦. We can generate an exchangeable random graph by first generating independent random vectors (Z_i:i≥ 1) from F and then generating A_ij independently from a Bernoulli distribution with parameter ⟨ Z_i, Z_j ⟩_𝒦.
We call this sampling procedure the graph root sampling with F.  In contrast, the graphon based sampling scheme is called graphon sampling with W.


The embedding of network nodes in a space 𝒦 has a clear interpretation.  Suppose each node i (1≤ i≤ n) corresponds to a Z_i=(X_i,Y_i)∈𝒦, with X_i, Y_i being the positive and negative components, respectively.  Then two nodes i and j are more likely to connect if ⟨ X_i,X_j ⟩ is large, or equivalently, X_iX_j⟨X̃_i,  X̃_j⟩ is large (where X̃_i=X_i/X_i).  The quantities X_i, X_j measure how “active” the individuals i, j are, respectively, while the normalized inner product ⟨X̃_i,  X̃_j⟩ measures how well the two individuals match each other.  Analogous interpretations can be given to the negative components Y_i, Y_j.

 Now we list how some commonly considered network models fit in the framework of GRD.  Further explanations of the correspondence are given in the Supplementary Material <cit.>.  Simulations based on these models are reported in <Ref>.


  
Stochastic block models: point mass mixture.  A stochastic block model <cit.> with k blocks is parameterized by (π,B), where π is in the (k-1)-dimensional simplex and B is a k× k symmetric matrix with entries in [0,1]. The exchangeable random graph is generated by sampling (e_i:i≥ 1) independently from a multinomial distribution with parameter π, and connecting nodes i,j independently with probability B_e_i,e_j. 

The corresponding GRD is a mixture of no more than k point masses in a k-dimensional space 𝒦, with the point mass locations determined by B and the point mass weights determined by π.

For example, consider an SBM with k=3, π=(1/3,1/3,1/3), and

    B=([ 1/4 1/2 1/4; 1/2 1/4 1/4; 1/4 1/4 1/6 ]) ,

which corresponds to three blocks, but the rank is only 2, with one positive eigenvalue and one negative eigenvalue.
Then a corresponding GRD is (after rounding)

    F=   (1/3) δ_(0.61;0.35) + (1/3) δ_(0.61;-0.35) + (1/3) δ_(0.41;0) ,

where δ_z=δ_(x;y) denotes the point mass at z=(x;y) and the semicolon is used to delineate positive and negative components.



  
Degree corrected block models: 1-D subspace mixture. The degree-corrected block model <cit.> is parameterized by (π,B,Θ) where π,B are the same as in SBM and Θ is a distribution on (0,∞).  The random graph is generated similarly as in SBM, except that it also generates (θ_i:i≥ 1) independently from Θ and the connection probability of nodes i,j is θ_iθ_j B_e_i,e_j.  In this case, the corresponding GRD is a mixture of distributions, each supported on a line connecting one of the SBM point masses and the origin.

For example, if we use the same k, π and B as in the SBM example above, and set Θ to be the uniform distribution on [0,1], then a GRD for this DCBM is

    F=   (1/3) U(0,(0.61;0.35))+ (1/3)U(0,(0.61;-0.35)) + (1/3) U(0,(0.41;0)) ,

where U(z,z') denotes the uniform distribution on the line segment between z and z'.  Other distributions Θ are allowed, and can be chosen differently for each mixture component.  This will lead to different ending points of line segments and the distributions on them.



  
Mixed membership block models: convex polytope. The mixed membership block model <cit.> allows each node to have a mixture of memberships.  

Given a matrix B as in the SBM, and a Dirichlet distribution Dir(a) with parameter a∈ (0,∞)^k, the random graph is generated by sampling (g_i:i≥ 1) from Dir(a) independently, and connecting nodes (i,j) with probability g_i^T B g_j.  In this case the corresponding GRD is a distribution supported on the convex polytope with extreme points given by the SBM point masses determined by B.

For example, using the same B matrix as in the previous examples for SBM and DCBM, and choosing a=(1,1,1) so that the Dirichlet distribution is uniform on the simplex, a GRD for this MMBM is

    F= U((0.61;0.35),(0.61;-0.35),(0.41;0)) ,

where U(z_1,z_2,z_3) denotes the uniform distribution on the convex hull of {z_1,z_2,z_3}.



  
Random dot-product graphs: finite dimensional subspace. The random dot-product graph <cit.> and generalized random dot-product graph <cit.> generate the random graph by connecting nodes (i,j) independently with probability ⟨ X_i, X_j ⟩, where (X_i:1≤ i≤ n) are node covariate vectors in an Euclidean space.  The original random dot-product graph only considers positive semidefinite inner products, while the generalized model allows for indefinite inner products in a similar fashion as we have defined for spaces.  In principle, a generalized random dot-product graph can be viewed as a finite-sample realization of a GRD supported on a finite dimensional space.







 §.§ Potential features and benefits of the graph root parameterization

As we will see in the following section, GRDs can be used to parameterize exchangeable random graphs under some mild regularity conditions. 
Here we list some features and potential benefits of the GRD parameterization.



  
GRDs are identifiable up to orthogonal transforms. Roughly speaking, two GRDs that differ by a pair of orthogonal transforms can lead to the same distribution of exchangeable random graphs.  The choice of orthogonal transforms used in GRD estimation is straightforward: It diagonalizes the covariance operator of the embedded node vectors. Moreover, many important geometric features of the distribution, such as clusters, pairwise distances, and the shape of support, are invariant under orthogonal transforms. In contrast, graphon estimators that do not attempt to recover the ordering of the nodes cannot reveal the same geometric structures in the data. This difference is illustrated in our numerical examples in <Ref>.



  
GRD provides new tools for some inference problems.
The embedding of network nodes as independent realizations of a common latent distribution makes it possible to apply the methods and tools developed for iid data to network related problems.

For example, suppose we observe an exchangeable random graph with n nodes, where each node i is associated with a covariate vector U_i∈ℝ^d. Here the nodes can be students in a school, edges represent friendship, and covariates are demographic informations.  A question of interest is to test whether the covariate and the network are independent.  In the GRD parameterization, the problem can be formulated as testing independence of two random vectors (U_i, Z_i) in a paired sample, where Z_i is the space embedding of node i.  In a second example, suppose we have two exchangeable random graphs on two disjoint set of individuals. Here the two networks can be student friendship networks from different middle schools. Then one may want to test whether these two networks have the same distribution.  In GRD parameterization, this reduces to testing equality (up to orthogonal transform) of two distributions using independent samples.  Moreover, if the underlying GRDs are the same, it is even possible to aggregate two independently estimated GRDs, as well as to predict edge probabilities between nodes in different samples.



  
GRD provides connection between the graphon, spectral embedding, and latent space perspectives. Spectral embedding of network vertices has been an active research topic related to exchangeable random graphs, especially stochastic block models <cit.>. Examples of spectral embeddings beyond stochastic block models include the random dot product graph <cit.> and the latent eigenmodel <cit.>. The GRD parameterization shows that such embeddings exist in an infinite dimensional space for all trace-class graphons, partially reconciling the graphon, spectral embedding, and latent space model literature.



§ EXISTENCE, IDENTIFIABILITY, AND TOPOLOGY OF GRD




From now on we only consider separable spaces, where 𝒦=ℋ_+⊖ℋ_-, ℋ_+=ℋ_-={x∈ℝ^∞: ∑_j x_j^2< ∞} with inner product ⟨ x, x' ⟩_ℋ_±=∑_j≥ 1x_j x_j'. These spaces are associated with the Borel σ-field.




 §.§ Existence of GRD for exchangeable random graphs

To find an underlying GRD for an exchangeable random graph, we consider the spectral decomposition of the corresponding graphon.  We will show that such GRDs exist for graphons whose spectral series converge in a strong sense.


Recall that a graphon W is a symmetric function from [0,1]^2 to [0,1].  We can view W as an integral operator on L^2([0,1]):

    (Wf)(·) = ∫_[0,1] W(·,s)f(s)ds ,  ∀ f∈ L^2([0,1]) .

Since ∫_[0,1]^2 W(s,s')^2 dsds'≤ 1, W is a Hilbert-Schmidt operator and hence  admits a spectral decomposition

    W(s,s') = ∑_j=1^∞λ_j ϕ_j(s)ϕ_j(s')-∑_j=1^∞γ_j ψ_j(s)ψ_j(s')

where λ_1≥λ_2≥...> 0, γ_1≥γ_2≥...> 0, and (ϕ_j:j≥ 1)∪(ψ_j:j≥ 1) are orthonormal functions in L^2([0,1]).  The convergence in (<ref>) shall be interpreted as L^2-convergence.  In general, almost everywhere convergence does not hold without further assumptions.


  We say a graphon W admits strong spectral decomposition if the eigen-components (λ_j,ϕ_j)_j≥ 1, (γ_j,ψ_j)_j≥ 1 in (<ref>) satisfy
  
    ∑_j≥ 1[λ_jϕ_j^2(s)+γ_jψ_j^2(s)]<∞    a.e. .


Strong spectral decomposition implies, among other things, that the sum in (<ref>) converges almost everywhere.

The spectral decomposition of a graphon has been considered in the mathematical side of the literature, such as in <cit.>, and recently in graphon estimation in <cit.>.  Here we use the spectral decomposition to define a mapping from [0,1] to a pair of infinite sequences: [√(λ_j)ϕ_j(s), j≥ 1]
and [√(γ_j)ψ_j(s), j≥ 1], and our key object, the graph root distribution (GRD), is the corresponding induced probability measure on the infinite dimensional space.  Such an induced probability measure carries all the information about the corresponding exchangeable random graph, and removes the ambiguity caused by measure preserving transforms.


  Any exchangeable random graph generated by a graphon that admits 
 strong spectral decomposition can be generated by a GRD on 𝒦.


For a graphon W, consider its spectral decomposition (<ref>), and define Z(s)=(X(s),Y(s)) as

    X_j(s) =    λ_j^1/2ϕ_j(s) ,  ∀ j≥ 1 ,
    
      Y_j(s) =    γ_j^1/2ψ_j(s) ,  ∀ j≥ 1 .

If s∼ Unif(0,1), the resulting Z(s)=(X(s),Y(s)) is a random object.
By the strong spectral decomposition assumption, X_ℋ_+ and Y_ℋ_- are finite with probability one, so Z is a well-defined random vector in 𝒦.  Moreover,

    ⟨ Z(s), Z(s') ⟩_𝒦 =∑_j [λ_jϕ_j(s)ϕ_j(s') - γ_jψ_j(s)ψ_j(s')]

converges almost everywhere since for s, s' we have

    |λ_jϕ_j(s)ϕ_j(s')|≤   1/2(λ_jϕ_j^2(s)+λ_jϕ_j^2(s')) ,  
    
    |γ_jψ_j(s)ψ_j(s')|≤   1/2(γ_jψ_j^2(s)+γ_jψ_j^2(s')) ,

and the summability is ensured for all s and s' satisfying (<ref>).

Let F be the probability measure induced by Z(s):[0,1]↦𝒦 with s∼ Unif(0,1). By construction, W(s,s')=⟨ Z(s), Z(s') ⟩_𝒦 almost everywhere, so that the graphon W and GRD F lead to the same sampling distribution of exchangeable random graph. 

The examples given in <Ref> are special cases of <Ref>. We provide
more detailed explanation of the correspondence in the Supplementary Material <cit.>.



 §.§ How stringent is strong spectral decomposition?
  The requirement of strong spectral decomposition is, indeed, quite mild.  The following proposition states that trace-class integral operators admit strong spectral decomposition.


  If W is trace-class, in the sense that ∑_j≥ 1(λ_j+γ_j)<∞ with λ_j,γ_j defined in (<ref>), then W admits strong spectral decomposition, and the GRD constructed from the spectral decomposition of W is square-integrable.


Trace-class integral operators are well-studied in functional analysis <cit.>. Some important subclasses are the following.

  
  * Finite rank graphons.  If W has finite rank, then it only has finitely many non-zero eigenvalues, and hence strong spectral decomposition holds trivially.  Important examples covered by this case include the stochastic block models, the degree corrected block models, and the mixed membership block models.
  
  * Smooth graphons. If a graphon W, or an element in its equivalence class, is in α-Hölder class:
  
    |W(x,y)-W(x,y')|≤ C |y-y'|^α

  for constants C>0, α>1/2 and all x,y,y'∈ [0,1], then it is trace-class and hence admits strong spectral decomposition. 
  
  * Continuous positive graphons. If W is positive semidefinite and continuous, then W is trace-class and hence admits strong spectral decomposition.  This is the famous Mercer's theorem.  One can relax the requirement of positivity by instead requiring W=W_+-W_- with W_+, W_- both being positive semidefinite and continuous.  An example covered in this case is
  
    W(x,y) = 1/log(e/x)log(e/y) ,   if (x,y)∈(0,1]^2

  and W(x,y)=0 if x=0 or y=0.  This graphon is not in any Hölder class but is trace-class.


The second case in the list above has a useful consequence.
Even though not all continuous graphons are trace-class, one can approximate any continuous graphon arbitrarily well using trace-class graphons.

  Let W be a continuous graphon. For any ϵ>0 there exists a trace-class graphon W' such that
  
    δ_□(W,W')≤ϵ .

  Moreover, the set of trace-class continuous graphons is a dense subset of continuous graphons.

If a continuous graphon W does not satisfy strong spectral decomposition, the approximation W' given in <Ref> with a small approximation error ϵ may have a very large trace, and its eigenvalues may decay very slowly. This can pose challenges in estimation.  We make further discussion in <Ref>.
  















 §.§ Identifiability of GRD



When do two graph root distributions F_1 and F_2 on 𝒦 lead to the same exchangeable random graph distribution?  We first exclude some trivial sources of ambiguity.



  
Ambiguity by concatenation 

Let (X,Y) be a random vector on 𝒦, and let R be any random variable in an Euclidean space or separable Hilbert space, the random vector (X',Y') in the augmented space with X'=(X,R), Y'=(Y,R) leads to the same random graph sampling distribution as (X,Y).  





  
Ambiguity by rotation Let Q be an inner product preserving mapping from 𝒦 to 𝒦: 
    ⟨ z,z' ⟩_𝒦=⟨ Qz,Qz' ⟩_𝒦 ,  ∀ z,z'∈𝒦 .

Then, in terms of the generated exchangeable random graph, a GRD F is indistinguishable from F_Q, the measure induced by transforming Z∼ F↦ QZ.
An obvious example of Q is the direct sum of two orthogonal transforms
Q_+, Q_- on ℋ_+, ℋ_- respectively, such that Q(x,y)=(Q_+x,Q_-y). Such a Q preserves the inner product ⟨·,·⟩_𝒦 because it preserves the inner products in both the positive and negative components. However, due to the indefinite inner product, this is not the only type of inner product preserving transforms on 𝒦.  Other transforms, such as hyperbolic rotations, can also preserve the indefinite inner product. See <cit.> for some examples of hyperbolic rotations under the context of random dot-product graphs.



To resolve the identifiability issue, a key observation is that both concatenation and hyperbolic rotation necessarily mix up the positive and negative components. So these ambiguities can be precluded by the requiring uncorrelated positive and negative components.  In this subsection, we show that the direct sum of a pair of orthogonal transforms is the only possible ambiguity in identifying a square-integrable GRD with uncorrelated positive and negative components.

 We say two distributions F_1, F_2 on 𝒦 are equivalent up to orthogonal transform, written as F_1o.t.=F_2, if there exist orthogonal transforms Q_+ on ℋ_+ and Q_- on ℋ_-, such that (X,Y)∼ F_1⇔ (Q_+X,Q_-Y)∼ F_2.



 Two square-integrable GRDs F_1, F_2 with uncorrelated positive and negative components give the same exchangeable random graph sampling distribution if and only if
  F_1o.t.=F_2.


The main idea of the proof is to establish a direct connection between a GRD F and its corresponding graphon W.  Now F is a probability measure on 𝒦, while W is a function from [0,1]^2→ [0,1].  Our idea is to use an inverse transform sampling mapping to relate the distribution F to a measurable function on [0,1].


  Let F be a distribution on 𝒦.  A measurable function Z:[0,1]↦𝒦 is called an inverse transform sampling mapping of F if
  
    s∼ Unif(0,1)⇒ Z(s)∼ F .


In other words, an ITS induces the Lebesgue measure on [0,1] to F on 𝒦.  The mapping Z(·) given by (<ref>) in the proof of <Ref> is an example of an ITS of the GRD F. 
If 𝒦 is one-dimensional, then a well-known example of ITS is the inverse cumulative distribution function.  It is also straightforward to see that ITS's are not unique since if Z(·) is an ITS of F and h(·) is measure-preserving then Z(h(·)) is also an ITS of F.  The following result ensures that ITS's always exist for distributions on a separable Hilbert space.

  Let F be a distribution on a separable Hilbert space, then there exists an ITS of F.


Here we give a sketch of proof of <Ref>.
For i=1,2, let Z_i(s)=(X_i(s),Y_i(s)) be an ITS of F_i and define graphon 

    W_i(s,s')=⟨ Z_i(s),Z_i(s') ⟩_𝒦 .

  By assumption that F_1 and F_2 lead to the same exchangeable random graph sampling distribution, we have W_1w.i.=W_2.  Also, by choosing appropriate orthogonal rotations in the positive and negative components we can make the covariance of Z_i diagonal so that W_i(s,s')=⟨ Z_i(s),Z_i(s') ⟩_𝒦 indeed corresponds to the spectral decomposition of W_i.  Then the desired result follows by invoking an exchangeable array representation theorem in the form of spectral decompositions due to Kallenberg <cit.>.

We summarize our representation results in the following corollary.

  There exists a one-to-one correspondence between trace-class graphons (under the equivalence relation “w.i.=”) and square-integrable GRD's with uncorrelated positive and negative components (under the equivalence relation “o.t.=”).





  
Canonical GRD Since any square-integrable GRD with uncorrelated positive and negative components is identifiable up to a pair of orthogonal transforms, we can choose appropriate orthogonal transforms so that the covariance of the GRD is diagonalized.  Such a choice can be used as a canonical representation.  If all eigenvalues of the covariance operator have multiplicity one, then the canonical GRD F is determined up to the sign of each coordinate.   As we will see in <Ref> below, our estimator recovers one of the canonical GRD's.  




 §.§ Topology of the GRD space: orthogonal Wasserstein distance

Having established the GRD representation of exchangeable random graphs, we can study the closeness of graph sampling distributions by looking at the closeness of GRD's. To this end, we consider a metric on the quotient space of square-integrable distributions on 𝒦 with respect to the equivalence relation “o.t.=”, which we call the orthogonal Wasserstein metric. We will show that convergence of a sequence of GRD's in this metric implies convergence of corresponding graphons in cut-distance.

We start by recalling the Wasserstein distance. We will only use a special case of the Wasserstein distance suitable for our purpose. Given two probability distributions F_1,F_2 on 𝒦, the Wasserstein distance between F_1,F_2 is

    d_ w(F_1,F_2):=inf_ν∈𝒱(F_1,F_2)𝔼_(Z_1,Z_2)∼νZ_1-Z_2 ,

where 𝒱(F_1,F_2) is the collection of all distributions on 𝒦×𝒦 with F_1 and F_2 being its two marginal distributions.

The following lemma says that if two square-integrable GRD's are close in Wasserstein distance, then the corresponding graphons are close in cut-distance.

  Let F_1 and F_2 be two square-integrable GRD's on 𝒦, with corresponding graphons W_1, W_2 defined using ITS as in (<ref>).  Then
  
    δ_□(W_1,W_2)≤ (𝔼_Z∼ F_1Z+𝔼_Z∼ F_2Z) d_ w(F_1,F_2) .




Since we do not distinguish two GRD's differing only by orthogonal transforms on positive and negative components, we consider the orthogonal Wasserstein distance

    d_ ow(F_1,F_2):=inf_ν∈𝒱(F_1,F_2)inf_Q_+,Q_-𝔼_(Z_1,Z_2)∼νZ_1-(Q_+⊕ Q_-)Z_2 ,

where Q_+, Q_- range over all orthogonal transforms on ℋ_+, ℋ_-, respectively, and (Q_+⊕ Q_-) denotes the orthogonal transform as the direct sum of Q_+ and Q_-: (Q_+⊕ Q_-)(x,y)=(Q_+x,Q_-y).

We can improve <Ref> to the orthogonal Wasserstein distance, which says that orthogonal Wasserstein distance induces a stronger topology than cut-distance.

  Let F_1, F be two square-integrable GRD's on 𝒦, with corresponding graphons W_1, W as obtained by ITS in (<ref>).
  Then 
    δ_□(W_1,W)≤ (𝔼_Z∼ F_1Z+𝔼_Z∼ FZ) d_ ow(F_1,F) .

As a consequence, if (F_N:N≥ 1) are square-integrable GRD's on 𝒦 with corresponding graphons (W_N:N≥ 1), then
  
    d_ ow(F_N,F)→ 0  ⇒ δ_□(W_N,W)→ 0 .






§ ESTIMATION OF GRAPH ROOT DISTRIBUTIONS

Given n≥ 1, suppose we have observed an n× n block of 𝐀: 𝐀_n=(A_ij:1≤ i,j≤ n), where 𝐀 is generated from a GRD F.  In such finite sample scenarios, GRD and graphon are used as modeling tools and are no longer linked to the infinite exchangeability, since the Aldous-Hoover theorem is only applicable to the infinite case.  We consider the following two inference questions.

  
  * Node embedding: Can we recover the realized sample of node vectors Z_1,..,Z_n in 𝒦?
  
  * Distribution estimation: Can we recover the GRD F with small orthogonal Wasserstein distance?




  
Notation For an infinite vector x, x^(p) denotes the first p elements of x.  For a matrix 𝐌 with countably infinite number of columns, 𝐌^(p) denotes the submatrix consisting of the first p columns.  For a matrix 𝐙=(𝐗,𝐘) with n rows and each row taking value in 𝒦, 𝐙^(p_1,p_2)=(𝐗^(p_1),𝐘^(p_2)). 




 §.§ Truncated weighted spectral embedding

Write 𝐀_n in its eigen decomposition

    𝐀_n = ∑_j=1^n_1λ̂_j,Aâ_jâ_j^T - 
    ∑_j=1^n-n_1γ̂_j,Ab̂_jb̂_j^T ,

where λ̂_1,A≥λ̂_2,A≥ ...≥λ̂_n_1,A≥ 0 are the non-negative eigenvalues of A, and γ̂_1,A≥ ...≥γ̂_n-n_1,A>0 are the absolute negative eigenvalues of A.


Let p_1,p_2<n be nonnegative integers to be specified later. We consider the weighted (p_1+p_2)-dimensional spectral embedding of the nodes

    𝐙̂_A=    [𝐗̂^(p_1)_A, 𝐘̂^(p_2)_A] ,
    𝐗̂^(p_1)_A =    [λ̂_1,A^1/2â_1,...,λ̂_p_1,A^1/2â_p_1] = [â_1,...,â_p_1]Λ̂_p_1,A^1/2 ,
    𝐘̂^(p_2)_A =    [γ̂_1,A^1/2b̂_1,...,γ̂_p_2,A^1/2b̂_p_2] = [b̂_1,...,b̂_p_2]Γ̂_p_2,A^1/2 ,

where Λ̂_p_1,A is the p_1× p_1 diagonal matrix with diagonal entries being (λ̂_1,A,...,λ̂_p_1,A), and Γ̂_p_2,A is defined similarly.

We use the rows of 𝐙̂_A to estimate the realized sample points Z_1,...,Z_n as follows.

    X̂_i,A =     (λ̂_1,A^1/2â_1i,...,λ̂_p_1,A^1/2â_p_1,i,0,...)∈ℋ_+ ,
    Ŷ_i,A =    (γ̂_1,A^1/2b̂_1i,...,γ̂_p_2,A^1/2b̂_p_2,i,0,...)∈ℋ_- .

In other words, X̂_i,A, Ŷ_i,A are the ith row of 𝐗̂_A^(p_1), 𝐘̂_A^(p_2), padded with zeros in the tails.  The same weighted spectral embedding has been considered in random dot-product graphs in finite dimensional spaces and at the sample level <cit.>. Here we focus more on the infinite  dimensional case, where p_1, p_2 need to grow with n, and study the statistical properties of the embeddings at a population level with a goal of estimating the GRD F.

 



 §.§ Reconstruction error of sample points

In order to show that the estimated node vectors (X̂_i,A,Ŷ_i,A) are close to the true but hidden realized node vectors (X_i, Y_i), it is necessary to identify a particular orthogonal transform Q=Q_+⊕ Q_- to work with.  To this end, we make the following assumption to clear the identifiability issue.

  
  (A1) For all j,j'≥ 1, 𝔼_(X,Y)∼ F(X_j X_j')=λ_j 1(j=j'), 𝔼_(X,Y)∼ F(Y_j Y_j')=γ_j 1(j=j'), 𝔼_(X,Y)∼ F(X_j Y_j')=0.

This assumption is non-technical, it merely says that we pick a canonical element among all possible orthogonal transforms on the positive and negative spaces.

Our next assumption is a polynomial eigen decay and eigen-gap condition.

  
  (A2) There exist positive numbers c_1≤ c_2, 1<α≤β such that for all j≥ 1
  
    c_1 j^-α≤ (λ_j∧γ_j)≤ (λ_j∨γ_j) ≤ c_2 j^-α ,
       (λ_j-λ_j+1)∧ (γ_j-γ_j+1)≥ c_1 j^-β ,

  with λ_j,γ_j defined in assumption (A1).

Assumption (A2) is often used in the literature of functional data analysis, where one needs to control the estimation error of individual eigenvectors for random variables in Hilbert spaces, using a truncated empirical eigen decomposition <cit.>. When λ_j∝ j^-α, the eigengap condition usually holds with β=α+1. The random vector Z=(X,Y)∼ F is square-integrable if α>1. 

Assumption (A2) may seem a bit too stringent. Indeed we only need to consider the first p_1+1 (p_2+1) positive (negative) eigenvalues.  The estimation error bound can be given as a function of all individual gaps between these eigenvalues, where equal or nearly equal eigenvalues can be treated by considering the corresponding principal subspace.  This will make the presentation too cumbersome and will not change much of the nature of our argument.  Another simplification made in Assumption (A2) is that the positive eigenvalues (λ_j:j≥ 1) and negative eigenvalues (γ_j:j≥ 1) decay at the same speed. Our analysis does allow for different decay speeds for the positive and negative eigenvalues. That will require choosing p_1 and p_2 separately, which involves a heavier notation. We choose to work with the version of Assumption (A2) stated above for presentation simplicity.



As mentioned in the discussion after <Ref>, if we use a trace-class graphon W' to approximate a continuous graphon W that does not satisfy strong spectral decomposition, the eigenvalues of W' may decay slowly, which corresponds to a small value of α in Assumption (A2), and leads to a slower rate of convergence of the estimation error bound.

Finally, our procedure requires accurate estimation of the eigenvectors, which in turn requires accurate estimation of the covariance operator.  We assume that the GRD has finite fourth moment.

  
  (A3) 𝔼_Z∼ FZ^4<∞.



  Let 𝐀_n be generated from a GRD F satisfying (A1-A3). Let 𝐙=(𝐗,𝐘) be the hidden node data matrix with the ith row being (X_i,Y_i)∈𝒦 (1≤ i≤ n).
 If 
  
    p_1=p_2=p=o(n^1/2β+α) ,
 then the estimator 𝐙̂_A given in (<ref>) satisfies
  
    n^-1𝐙̂_A-𝐙^(p,p)_F^2 = O_P(n^-α-1/2β+p^2β+1n^-1)

  where ·_F denotes the Frobenius norm.  As a consequence, let
  F̂_A^(p) be the empirical distribution putting 1/n probability mass at each row of 𝐙̂_A, and F̂^(p) putting 1/n mass at each row of
  𝐙^(p,p), then
  
    d_ w(F̂_A^(p), F̂^(p))= O_P(n^-α-1/4β+p^β+1/2n^-1/2) .





The main technical task in the proof is to control the difference between the true realized random vectors (X_i,Y_i) and their projections on principal subspaces obtained from several approximations of the gram matrix.  Thus the tools used are similar to those in functional data analysis <cit.>.  However, the additional challenge here is that we do not observe any of the empirical covariance matrices, and the adjacency matrix we observe is actually a noisy version of the indefinite gram matrix, which is the difference of two positive semidefinite gram matrices, one in the positive space and one in the negative space.  This issue does not exist in the ordinary functional data analysis literature and requires more delicate spectral perturbation analysis.



 §.§ Estimating the GRD

The second part of <Ref> gives the possibility of estimating the GRD F using 
F̂_A^(p).  According to <Ref>, we only need to show that F̂^(p) 
is close to F.  To this end, we consider an intermediate object F^(p), the distribution of 
truncated vector (X^(p),Y^(p)) with (X,Y)∼ F.  The argument proceeds in two steps.

The first step is to compare F and F^(p), which is straightforward.

Under assumptions (A1-A2),

    d_ w(F,F^(p))≤ c p^-(α-1)/2 .











The second part is comparing the population truncated distribution F^(p) and its empirical version F̂^(p).  
We apply the result of <cit.> which provides Wasserstein error bounds for empirical distributions.  
Here we state a special case, which is suitable for our purpose.

Under assumptions (A1-A3), there exists a constant c independent of n,p, such that

    𝔼 d_ w(F̂^(p), F^(p)) ≤ c n^-1/(p ∨ 2)[1+(log n) 1(p=2)] ,

where 1(·) is the indicator function .



Combining the above two lemmas with <Ref>, we have the following result on estimating the GRD.

Under assumptions (A1-A3), we have, when p≥ 3 satisfies the conditions in <Ref>,

    d_ w(F̂_A^(p), F)=O_P[n^-α-1/4β+p^β+1/2n^-1/2+p^-(α-1)/2+n^-1/p] .

The right hand side is o_P(1) if p→∞ and p=o(log n).


This result seems to suggest that one must have p→∞ to have a vanishing error.
The term p^-(α-1)/2 comes from <Ref> which corresponds to the truncation error.  It is necessary only because we assumed a particular eigenvalue sequence in Assumption (A2).
What really matters here is the sum of absolute eigenvalues beyond p: (λ_j,γ_j:j> p).  When the graphon
is nearly low rank or the GRD is supported close to a finite dimensional space, there is no need to use a large value of p.



 §.§ Estimation for sparse graphs

One limitation of the theoretical framework of exchangeable random graphs is that they can only model dense graphs.  The total number of edges in 𝐀_n will concentrate around n^2∫_[0,1]^2W.  In reality the number of edges in a network rarely grows as the squared number of nodes.  Therefore, sparse networks are of greater practical interest.  To this end, for a given graphon W and a node sample size n one can consider adding a “sparsity parameter” to the network sampling scheme <cit.>:

    𝐀_n,i,j∼ Bernoulli(ρ_n W(s_i,s_j)) ,  ∀ 1≤ i<j≤ n .


This sparsity parameter can be carried over to the graph root sampling scheme.  Let F be a GRD. For a node sample size n and sparsity parameter ρ_n, the corresponding sparse graph root sampling scheme is equivalent to generating node sample points from a scaled distribution:

    𝐀_n,i,j∼ Bernoulli(⟨ρ_n^1/2Z_i,ρ_n^1/2Z_j⟩_𝒦) ,

  where Z_iiid∼F.
For notational simplicity, for scalar a and distribution F we use a F to denote the distribution obtained by scaling the distribution F by a factor of a: Z∼ F⇔ aZ∼ aF.

In the SBM and DCBM literature, it is well known that consistent estimation of network communities is possible only if nρ_n→∞.
  Our estimation theory developed in the previous subsections can be extended to cover sparse sampling schemes.  One technical challenge is that when nρ_n=o(log n), spectral methods tend to be sensitive to overly large node degrees.  In the spectral clustering literature <cit.> a common approach is to zero out rows and columns of 𝐀_n for which the degrees are too high. Some data-driven degree thresholding rules are developed, for example, in <cit.>. In the following we consider an adaptive trimmed spectral embedding method.

  Let d_i=∑_1≤ j≤ n𝐀_n,i,j be the degree of node i. 
Let I_n be the set of nodes whose degrees are among the ⌊n(n-1)/∑_i=1^n d_i⌋ largest, and 𝐀̃_n be the adjacency matrix obtained by zeroing out the columns and rows in I_n.  Let 𝐙̃_A^(p,p) and F̃_A be the corresponding embeddings and GRD estimate defined in <Ref> and <Ref>, respectively, with 𝐀_n replaced by 𝐀̃_n.

Under assumptions (A1-A3), assuming sparse sampling scheme (<ref>), the following hold.


  * If nρ_n→∞ and
p_1=p_2=p=o[n^1/(2β+α)∧ (nρ_n)^1/(2β)] then 
    n^-1ρ_n^-1/2𝐙̃_A^(p,p)-𝐙^(p,p)_F^2
    
      =   O_P(p^2β-α+1(nρ_n)^-1+ n^-α-1/2β+p^2β+1n^-1) .


  * If in addition we assume p≥ 3 then

    d_ w(ρ_n^-1/2F̃_A, F)
    
    =   O_P(p^β-(α-1)/2(nρ_n)^-1/2+ n^-α-1/4β+p^β+1/2n^-1/2+p^-α-1/2 +n^-1/p)   ,

where the error bound is o_P(1) if p→∞ and p=o(log n ∧ (nρ_n)^1/(2β)).





  
Comparison to graphon estimation using spectral methods. Graphon estimation using singular value thresholding has been considered in <cit.>. For specificity of discussion we focus on <cit.>.  The method first performs a singular value decomposition of the adjacency matrix 𝐀_n, and keeps only the components whose singular values exceed a threshold τ.  Then the remaining low-rank approximation is multiplied by ρ_n^-1 and entry-wise trimmed to [0,1] to obtain 𝐆̂_n as an approximation of the probability matrix 𝐆_n=ρ_n^-1𝔼𝐀_n, which can be further used as a piecewise constant approximation to the graphon. The error rates reported in <cit.> and discussed here refer to the probability matrix estimation error.

We first explain the difference  between the GRD estimation problem and probability matrix estimation problem.  In GRD estimation, an intermediate error quantity is the empirical distribution approximation error n^-1∑_i=1^nẐ_i-Z_i^2, where Z_i's are the true latent vectors sampled from the underlying GRD and Ẑ_i's are the estimated versions. The total GRD estimation error bound needs to add another term due to the empirical distribution approximation.
On the other hand, the probability matrix estimation is concerned with the error metric n^-2𝐆̂_n-𝐆_n_F^2.  Assuming that the GRD estimation and spectral probability matrix estimation use the same rules to select significant eigen components, and ignoring the trimming step in obtaining 𝐆̂_n, we have 𝐆̂_n,i,j=⟨Ẑ_i, Ẑ_j ⟩_𝒦. Using Cauchy-Schwartz we obtain

    n^-2𝐆̂_n-𝐆_n_F^2
    ≤    2[n^-1∑_i=1^nẐ_i-Z_i^2][n^-1∑_i=1^nZ_i^2+n^-1∑_i=1^n Ẑ_i^2] .

If we accept the assumption that n^-1∑_i=1^n Ẑ_i^2 is close to n^-1∑_i=1^n Z_i^2≈𝔼Z_1^2, then the GRD empirical distribution approximation error provides an upper bound of the probability matrix estimation error, up to a multiplicative factor.

Our requirement of polynomially decaying eigenvalues (first part of Assumption A2) implies the tail sum condition of eigenvalue sequence in <cit.>. In addition, we require a lower bound of the eigenvalue gap (second part of Assumption A2), which is not required in <cit.>.  This is because in GRD estimation we need to recover the eigenvectors, and the use of subspace perturbation theory (Davis–Kahan sinΘ theorem) involves a multiplicative factor of the inverse of eigenvalue gap near the threshold.  This multiplicative factor of inverse eigenvalue gap leads to a slower convergence rate.  When λ_j=c j^-α for some α>1, Assumption A2 holds with β=α+1.  In the moderately sparse case nρ_n is only a polynomial of log n, then <Ref> implies a GRD empirical distribution estimation error rate of (nρ_n)^-a-1/2(a+1), while the probability matrix estimation error rate in <cit.> is (nρ_n)^-2a-1/2a.
In <Ref> we will empirically observe that when the underlying graphon is low rank (e.g., a stochastic block model with a small number of blocks), then the GRD empirical distribution approximation error and probability matrix estimation error roughly differ by a multiplicative factor; and when the underlying graphon is high rank, such as a graphon whose eigenvalues decay polynomially, then the GRD empirical distribution approximation error exhibits a slower rate of convergence than the corresponding probability matrix estimation error.






 §.§ Choice of embedding dimensions


In practice, the value of p can affect the quality of the estimated GRD. If p is too small, the estimate may not have sufficient dimensionality to carry all useful structures in the GRD.  If p is too large, the estimation becomes less stable and there would be a waste on computing and storage resources.  Moreover, in many applications it may make sense to use different values of p_1 and p_2, since the effective dimensionality can be different for the positive and negative components as seen in the numerical examples in <Ref>.

One potential way of choosing (p_1,p_2) is to follow a common practice in functional data analysis, where one chooses the leading principal subspace that explain a certain fraction (such as 90%) of the total variance. In network data, this approach has limited success due to the low-rank and high-noise nature of the adjacency matrix. Real-world network data often have low rank structures, but are observed with an entry-wise Bernoulli noise.  

Another way to choose p_1,p_2 is singular value thresholding <cit.>, where p_1 and p_2 are the number of positive and negative eigenvalues whose absolute value exceeds a threshold, respectively.  In particular, <cit.> suggests the threshold 2.01√(nσ̅^2), where σ̅^2≥max_i,j Var(A_ij) is an upper bound of the maximum entry-wise variance of the adjacency matrix.

When σ̅^2 is unavailable, one may use the conservative bound Var(A_ij)≤ 1/4 ≡σ̅^2, which results in the conservative singular value threshold 1.005 √(n)≈√(n).

We find this simple rule of singular value thresholding working quite well for some real data sets.  









§ NUMERICAL EXAMPLES



 §.§ Simulation 1: dense SBM, DCBM, and MMBM

We apply the truncated weighted spectral embedding for SBM, DCBM, and MMBM. For comparison, we also apply the graphon estimation method based on stochastic block model approximation (SBA) <cit.>.  

The GRD estimator F̂ and the graphon estimator Ŵ are in different spaces, where F̂ is a probability measure on an Euclidean space and Ŵ is a symmetric function defined on the unit square.  Thus they may reveal different underlying structures about the network data.  We shall see that the GRD estimator is able to reveal the clustering and subspace clustering of the network nodes, while the graphon estimator is less visually informative without a correct ordering of the nodes.  Estimation errors and convergence rates in more general settings, including sparse networks and infinite dimensional GRDs, are considered in <Ref>.









Following the notations in <Ref>, we set k=3 and
B given in (<ref>).
The remaining parameters are set as follows.

  
  * π=(1/3,1/3,1/3) for the SBM and DCBM. 
  
  * Θ= Unif(0.7,1.4) for the DCBM, so the effects of node activeness parameter θ_i's range from halving to doubling the corresponding SBM edge probabilities.
  
  * a=(0.5,0.5,0.5) for the MMBM, so that the mixed memberships are not too close to the extreme points. 

For each model, we generate a random graph with n=1000 nodes, and apply the truncated weighted spectral embedding.  The number of eigen-components is determined by the singular value thresholding rule as described in <Ref> with threshold √(n), which chooses top two absolute eigenvalues in all three cases, with one positive component and one negative component.  The smooth graphon estimation method requires two independent realizations of the adjacency matrix on the same set of nodes.  To make it a fair comparison, we generate two adjacency matrices of size 708× 708 to use in the smooth graphon estimation algorithm, so that the number of independent observations is the same as a 1000× 1000 adjacency matrix.






















A typical output of GRD estimation and the SBA algorithm are visualized in <Ref>, <Ref>, and <Ref> for the SBM, DCBM, and MMBM, respectively. The SBA algorithm requires a tuning parameter δ for grouping similar nodes. Here we use δ=0.1 since it achieves the smallest probability matrix estimation error.
In each figure, the left plot shows the truncated and weighted spectral embedding of network nodes.  The red dots, line segments, and triangles are corresponding supports of the true GRD as theoretically predicted in <Ref>. 
The embedded empirical distributions exhibit reasonable approximations to the underlying graph root distributions.  On the other hand, the graphon estimation method SBA outputs an estimated probability matrix.










 §.§ Simulation 2: sparse and infinite dimensional graphons.

In this simulation study we demonstrate GRD estimation in sparse and infinite dimensional settings, and compare with the simulation results in corresponding graphon estimation using singular value thresholding (USVT, <cit.>).  We adopt two simulation settings in <cit.>: a stochastic block model with four communities and a smooth graphon.



  
Stochastic block model.
In the stochastic block model setting, we consider stochastic block models with k=4 equal sized communities, and the B matrices have randomly generated entries from the uniform distribution on [0,1] subject to symmetry. We consider four different values of ρ: 0.4, 0.2, 0.1, 0.05, and six values of n such that log (nρ/k) takes equally spaced values between 2.2 and 3.2. For each combination of (n,ρ) the simulation is repeated 30 times with independently generated B, community membership, and 𝐀_n.  The singular value thresholding algorithm for probability matrix estimation uses threshold 2.01√(nρ).  This is to make sure we can reproduce the results in <cit.>.  For GRD estimation, we choose p_1 and p_2 by thresholding the absolute eigenvalues at 2.01√(nρ(1-ρ)), following the suggestion in <cit.>.

The results are summarized in <Ref>.  The error metrics reported here are empirical GRD approximation error and the probability matrix estimation error as introduced in <Ref>.  The similar slopes between empirical GRD errors and probability matrix estimation errors seem to suggest that in this low-dimensional case, the two estimation errors roughly differ by a constant factor.



  
Smooth graphon. In the smooth graphon setting, we use W(x,y)=min(x,y), whose jth eigenvalue is 4/π^2(2j-1)^2. We consider the same values of ρ: 0.4, 0.2, 0.1, 0.05.  Given the small eigenvalues, we consider larger values of log(nρ_n), which are equally spaced between 4 and 8.5.  Due to the computer memory limit, we carry out the experiment when n<1.5× 10^4. That is, for ρ=0.4 the experiment covers n such that log (nρ)∈ [4,8.5]; For ρ=0.2 it covers log(nρ)∈[4,8]; For ρ=0.1 it covers log(nρ)∈[4,7]; For ρ=0.05 it covers log(nρ)∈[4,6.5].  The results are summarized in the left plot of <Ref>.  The plot seems to confirm a slower rate of convergence for GRD estimation error.

The estimation errors of both the empirical GRD and probability matrix exhibits a sharp drop when log(nρ)≈ 7.  To better understand this we decompose the total estimation error into two parts.

  
  * The finite dimensional estimation error: This is the error in approximating the low-rank component of the probability matrix 𝐆_n=ρ_n^-1𝔼𝐀_n.  
  For empirical GRD estimation, this corresponds to n^-1∑_i=1^nρ^-1/2Ẑ_i,A^(p_1,p_2)- Z_i^(p_1,p_2)^2, where Ẑ_i,A^(p_1,p_2) is the estimated latent vector Z_i truncated to retain p_1 and p_2 coordinates in the positive and negative parts, respectively, as given in (<ref>). For probability matrix estimation, this corresponds to 
  n^-2𝐆̂_n-𝐆_n^(p)_F^2, where 𝐆_n^(p) is the best rank-p approximation to 𝐆_n in Frobenius norm, and p is the number of singular values used in the USVT method.
  
  * The truncation error: This is the error incurred by ignoring the eigen-components of 𝐆_n with smaller absolute eigenvalues.  In empirical GRD estimation, this error is n^-1∑_i=1^n𝐙_i-𝐙_i^(p_1,p_2)^2.  In probability matrix estimation, this error is n^-2𝐆_n-𝐆_n^(p)_F^2.




The finite dimensional error and truncation error are plotted in the middle and right plots in <Ref>, respectively.  Near the point log (nρ)=7, the signal becomes strong enough to pick up the second eigenvalue of the underlying probability matrix, therefore the finite dimensional approximation error increases for both methods, because there are more eigen components to estimate.  After this increase, the finite dimensional approximation errors start dropping again with a similar linear slope.  The behavior of the truncation error matches the intuition, as it stays constant until a new eigen component is picked up when the signal strength increases.  In this example, the truncation error is larger and decays more slowly for the empirical GRD estimation than for the probability matrix estimation. 



 §.§ The political blogs data

The political blogs data <cit.> is one of the most widely studied network data sets with a well-believed degree-corrected community structure <cit.>.  The data set records undirected hyperlinks among 1222 political blogs during the 2004 presidential election, and the nodes have been manually classified as “liberal” and “conservative”.

Among many statistical methods applied to this data set, spectral methods are quite popular and have used the top two singular vectors of the adjacency matrix.  Here we apply the truncated and weighted spectral embedding to this data set.  The singular value thresholding rule suggests two significant eigen-components, both of which correspond to the positive component.  The embedded nodes in the two-dimensional space reflects a mixture of two components each on a one dimensional subspace, with each mixture component corresponding to a labeled class.  For this data set we only have one realization of the adjacency matrix so the SBA algorithm is no longer applicable.  For comparison, we apply the sorting-and-smoothing (SAS) estimator developed by <cit.>, which adapts the SBA method by sorting the nodes according to the degrees.  We also apply the USVT method to estimate the probability matrix, with the singular value threshold 1.005√(n).

The results are visualized in <Ref>, in a similar fashion as in the simulated examples in <Ref>.  The GRD node embedding scatter plot is colored according to the ground truth of manual labeling of the blogs. It clearly shows that each group is represented by a one-dimensional subspace on which the GRD is supported.  The SAS  estimator sorts the nodes according to the degrees, and misses the subspace clustering hidden in the data.  The USVT probability matrix estimation output is similar to that of SAS, with a different but random sorting of the nodes.





 §.§ The political books data

The political books data records undirected links among 105 political books with links defined by the co-purchase records on .  This data set, available on Mark Newman's website[<http://www-personal.umich.edu/ mejn/netdata/>], was collected by Krebs <cit.> during the 2004 presidential election.  The nodes have been manually labeled as one of the three categories: “neutral”, “liberal”, and “conservative”.

Given the three labeled classes, it seems natural to assume three significant eigen-components.  However, the singular value thresholding rule indicates only two significant components, both with positive eigenvalues.  Again, we also apply the SAS algorithm and the USVT probability matrix estimator with the threshold 1.005√(n) to this data set.

As shown in <Ref>,
the truncated weighted spectral embedding of the first two components (left plot of <Ref>) shows a two-component mixture with each component supported on a one-dimensional subspace, which strongly indicates a two-block DCBM.  The “neutral” class, plotted as green square points, appears near the intersection of the other two classes.  The SAS estimator and the USVT probability matrix estimator do not explicitly indicate such subspace clustering structure.





§ DISCUSSION




  
Kernel based learning  A side result of our theory is the relationship between the generating distribution of a random sample and the distribution of the corresponding kernel/gram matrix. Let F be a probability measure on a separable Hilbert space 𝒳.  Let (X_i:i≥ 1) be a sequence of independent samples from F, and 𝐆=(⟨ X_i, X_j ⟩, i,j≥ 1) be the (infinite size) gram matrix.  The perspective of viewing 𝐆 as an exchangeable random array allows us to establish the correspondence between F and the distribution of 𝐆.  The result essentially says that the gram matrix carries all information about F up to an orthogonal transform.  We believe that this result is elementary and highly intuitive, but are not able to find it in the literature.

Let F be a probability measure on a separable Hilbert space 𝒳.  Denote 𝒢_F the distribution of the corresponding infinite gram matrix 𝐆.  Then for two probability measures F_1, F_2 on 𝒳, 𝒢_F_1=𝒢_F_2 if and only if F_1o.t.=F_2, provided that one of the following holds:

  
  * 𝒳 is finite dimensional;
  
  * 𝔼_X∼ F_1X^2<∞;
  
  * 𝔼_X∼ F_2X^2<∞.
 

Here the equivalence relation “o.t.=” is defined as in <Ref> by treating ℋ_+=𝒳 and ℋ_-=∅.



  
Modeling and inference for relational data The framework of graph root representation can be extended in several interesting directions.  First, one can model the connection probability with a logistic link function so that the two nodes i,j connect with probability (1+e^-⟨ Z_i, Z_j ⟩_𝒦)^-1.  With such a logistic transform, each distribution on 𝒦 can be used to generate an exchangeable random graph, and therefore can model a wider collection of structures.  Moreover, one can also use the same framework to model relational data beyond binary observations.  For example, one may observe event counting between a pair of nodes, such as number of email correspondences and frequency of research article citations.  In applications such as multivariate time series and multimodal imaging, one may even observe a vector for each pair of nodes.

The graph root embedding also facilitates many subsequent inferences. We have discussed two examples in <Ref>.  There are other potential uses of GRD representations of networks. For example, in addition to clustering the embedded nodes as in SBM and DCBM, one can also test for specific structures of the graph root distribution, or compare the graph root distributions for multiple networks.  See <cit.> for an example of two-sample comparison for random dot-product graphs.  Another way to make use of the graph root embedding is to model the node movement in temporal networks.  A challenge is to find the orthogonal transforms to match the embeddings at different time points.  See <cit.> for an example using a similar idea with a different latent space model.










