Joint Embedding of Graphs

TL;DR

The paper introduces a joint embedding method for multiple graphs that captures shared features and vertex information, improving graph classification and interpretability in applications like brain network analysis.

Contribution

It proposes a novel joint embedding technique for multiple graphs and a new random graph model, with theoretical guarantees and practical applications.

Findings

Produces low-error parameter estimates under the model

Achieves state-of-the-art graph classification performance

Extracts interpretable features from brain graphs

Abstract

Feature extraction and dimension reduction for networks is critical in a wide variety of domains. Efficiently and accurately learning features for multiple graphs has important applications in statistical inference on graphs. We propose a method to jointly embed multiple undirected graphs. Given a set of graphs, the joint embedding method identifies a linear subspace spanned by rank one symmetric matrices and projects adjacency matrices of graphs into this subspace. The projection coefficients can be treated as features of the graphs, while the embedding components can represent vertex features. We also propose a random graph model for multiple graphs that generalizes other classical models for graphs. We show through theory and numerical experiments that under the model, the joint embedding method produces estimates of parameters with small errors. Via simulation experiments, we demonstrate that the joint embedding method produces features which lead to state of the art performance in classifying graphs. Applying the joint embedding method to human brain graphs, we find it extracts interpretable features with good prediction accuracy in different tasks.

Figures

Tables

Table 1. TABLE I: Objective function and running time of four initialization approaches. The standard error is shown in parenthesis. SVD and truth initializations are significantly better than random initializations on this scenario.

Initialization	Objective	Running time (sec.)
SVD	$375.22\left(1.21\right)$	$8.3\left(1.0\right)$
$1$ Random	$383.29\left(1.60\right)$	$9.2\left(1.4\right)$
Best of $10$ Random	$379.63\left(1.39\right)$	$96.5\left(5.3\right)$
Truth	$374.69\left(1.22\right)$	$7.8\left(1.0\right)$

Table 2. TABLE II: Clustering Performance on Wikipedia Graphs. The adjusted rand index (ARI) and the purity of clustering of different spectral clustering approaches are shown. The first two scenarios consider each graph separately, and the third considers them jointly. The standard error (included in parenthesis) is estimated through repeatedly clustering bootstrapped latent positions. The joint embedding estimates latent positions which combine the information in both graphs with good clustering performance.

Method	ARI	Purity
ASE+EN	0.152 (0.02)	0.555 (0.018)
JE+EN	0.154 (0.019)	0.55 (0.019)
ASE+FR	0.114 (0.019)	0.521 (0.017)
JE+FR	0.154 (0.019)	0.549 (0.017)
ASE+(EN+FR)	0.155 (0.02)	0.547 (0.016)
JE+(EN, FR)	0.156 (0.02)	0.551 (0.021)

Equations

$${\mathbf{A}}_{st}\sim Bernoulli(x_s^T{x}_t).$$

$$P(\mathbf{A}|\mathbf{X})=\prod _{s<t}(x_s^T{x}_t{)}^{{\mathbf{A}}_{st}}(1-x_s^T{x}_t{)}^{1-{\mathbf{A}}_{st}}.$$

$$\{{\mathbf{A}}_i{\}}_{i=1}^m\sim \mathrm{MREG}(\{{\lambda }_i{\}}_{i=1}^m,h_1,...,h_d),$$

$${\mathbf{A}}_i[s,t]\sim Bernoulli(\sum _{k=1}^d{\lambda }_i[k]h_k[s]h_k[t]).$$

$$\{({\lambda }_i,{\mathbf{A}}_i){\}}_{i=1}^m\sim MREG(F,h_1,...,h_d).$$

$$P(\mathbf{A}|\tau ,\mathbf{B})=\prod _{i<j}{\mathbf{B}}_{{\tau }_s,{\tau }_t}^{{\mathbf{A}}_{s,t}}(1-{\mathbf{B}}_{{\tau }_s,{\tau }_t}{)}^{(1-{\mathbf{A}}_{s,t})}$$

$$({\hat{\lambda }}_1,...,{\hat{\lambda }}_m,{\hat{h}}_1,...,{\hat{h}}_d)=\underset{{\lambda }_i,\Vert h_k\Vert =1}{\mathrm{argmin}}\sum _{i=1}^m\Vert {\mathbf{A}}_i-\sum _{k=1}^d{\lambda }_i[k]h_k{h}_k^T{\Vert }^2.$$

$$(\hat{\mathbf{\Lambda }},\hat{\mathbf{H}})=\underset{\mathbf{\Lambda },\Vert h_k\Vert =1}{\mathrm{argmin}}\sum _{i=1}^m\Vert {\mathbf{A}}_i-\sum _{k=1}^d{\mathbf{\Lambda }}_{ik}{h}_k{h}_k^T{\Vert }^2.$$

$$\underset{{\mathbf{D}}_i,\Vert h_k\Vert =1}{\mathrm{argmin}}$$

$${\mathbf{D}}_i\text{ being diagonal.}$$

$$\underset{\mathbf{\Lambda },\Vert h_k\Vert =1}{\mathrm{argmin}}\Vert \mathbb{A}-\sum _{k=1}^d{\mathbf{\Lambda }}_{\ast k}\otimes h_k\otimes h_k{\Vert }^2,$$

$$({\hat{\lambda }}_1,...,{\hat{\lambda }}_m,{\hat{\mathbf{S}}}_1,...,{\hat{\mathbf{S}}}_d)=\underset{{\lambda }_i,\mathbf{S}}{\mathrm{argmin}}\sum _{i=1}^m\Vert {\mathbf{A}}_i-\sum _{k=1}^d{\lambda }_i[k]{\mathbf{S}}_k{\Vert }^2,$$

$$({\hat{\mathbf{\Lambda }}}_{\ast k_0},{\hat{h}}_{k_0})=$$

$$\sum _{i=1}^m\Vert {\mathbf{A}}_i-\sum _{k=1}^{k_0-1}{\hat{\mathbf{\Lambda }}}_{ik}{\hat{h}}_k{\hat{h}}_k^T-{\mathbf{\Lambda }}_{i{k}_0}{h}_{k_0}{h}_{k_0}^T{\Vert }^2.$$

$$({\hat{\mathbf{\Lambda }}}_{\ast 1},...,{\hat{\mathbf{\Lambda }}}_{\ast k_0})=\underset{\mathbf{\Lambda }\in {\mathbb{R}}^{m\times k_0}}{\mathrm{argmin}}\sum _{i=1}^m\Vert {\mathbf{A}}_i-\sum _{k=1}^{k_0}{\mathbf{\Lambda }}_{ik}{\hat{h}}_k{\hat{h}}_k^T{\Vert }^2.$$

$$\frac{\partial f}{\partial h_{k_0}}=-4\sum _{i=1}^m{\mathbf{\Lambda }}_{i{k}_0}({\mathbf{R}}_{ik}-{\mathbf{\Lambda }}_{i{k}_0}{h}_{k_0}{h}_{k_0}^T)h_{k_0}.$$

$$\frac{\partial f}{\partial {\mathbf{\Lambda }}_{i{k}_0}}=-2⟨{\mathbf{R}}_{ik}-{\mathbf{\Lambda }}_{i{k}_0}{h}_{k_0}{h}_{k_0}^T,h_{k_0}{h}_{k_0}^T⟩.$$

$${\hat{\mathbf{\Lambda }}}_{i{k}_0}=⟨{\mathbf{R}}_{ik},h_{k_0}{h}_{k_0}^T⟩,$$

$$({\hat{\mathbf{\Lambda }}}_{i1},...,{\hat{\mathbf{\Lambda }}}_{i{k}_0})$$

$$\mathbf{\Gamma }{\hat{\mathbf{\Lambda }}}_{\ast ,1:k_0}^T=\mathbf{\Psi }.$$

$$\frac{\partial f}{\partial h_{k_0}}$$

$${\hat{\mathbf{\Lambda }}}_{i{k}_0}$$

$$(\hat{\mathbf{\Lambda }},\hat{\mathbf{H}})=\underset{\mathbf{\Lambda },\Vert h_k\Vert =1}{\mathrm{argmin}}\sum _{i=1}^m\Vert {\mathbf{A}}_i-\sum _{k=1}^d{\mathbf{\Lambda }}_k{h}_k{h}_k^T{\Vert }^2,$$

$$(\hat{\mathbf{\Lambda }},\hat{\mathbf{H}})=\underset{\mathbf{\Lambda },\Vert h_k\Vert =1}{\mathrm{argmin}}\Vert \frac{1}{m}\sum _{i=1}^m{\mathbf{A}}_i-\sum _{k=1}^d{\mathbf{\Lambda }}_k{h}_k{h}_k^T{\Vert }^2.$$

$$(\hat{\mathbf{\Lambda }},\hat{\mathbf{H}})=\underset{\mathbf{\Lambda },\Vert h_k\Vert =1}{\mathrm{argmin}}\sum _{i=1}^m\Vert {\mathbf{A}}_i-\sum _{k=1}^d{\mathbf{\Lambda }}_{y_{i}k}{h}_k{h}_k^T{\Vert }^2.$$

$${\hat{\mathbf{\Lambda }}}_{yk}=\left(\sum _{i=1}^m\mathbb{I}\{y_i=y\}⟨{\mathbf{R}}_{ik},h_{k_0}{h}_{k_0}^T⟩\right)/\left(\sum _{i=1}^m\mathbb{I}\{y_i=y\}\right).$$

$$(\hat{\mathbf{\Lambda }},\hat{\mathbf{H}})=\underset{\mathbf{\Lambda }\ge 0,\Vert h_k\Vert =1}{\mathrm{argmin}}\sum _{i=1}^m\Vert {\mathbf{A}}_i-\sum _{k=1}^d{\mathbf{\Lambda }}_{ik}{h}_k{h}_k^T{\Vert }^2.$$

$$\rho ({\mathbf{A}}_i,h)=\Vert {\mathbf{A}}_i-⟨{\mathbf{A}}_i,h{h}^T⟩h{h}^T{\Vert }^2,$$

$$D_m(h,h_1)=\frac{1}{m}\sum _{i=1}^m\rho ({\mathbf{A}}_i,h),$$

$$D(h,h_1)=\mathbb{E}(\rho ({\mathbf{A}}_i,h)).$$

$${\hat{h}}_1^m=\underset{\Vert h\Vert =1}{\mathrm{argmin}}\underset{{\lambda }_i}{\mathrm{argmin}}\sum _{i=1}^m\Vert {\mathbf{A}}_i-{\lambda }_{i}h{h}^T\Vert .$$

Full text

1 Introduction

In many problems arising in science and engineering, graphs arise naturally as data structure to capture complex relationships between a set of objects. Graphs have been used in various application domains as diverse as social networks [1], internet mapping [2], brain connectomics [3], political voting networks [4], and many others. The graphs are naturally high dimensional objects with complicated topological structure, which makes graph clustering and classification a challenge to traditional machine learning algorithms. Therefore, feature extraction and dimension reduction techniques are helpful in the applications of learning graph data. In this paper, we propose an algorithm to jointly embed multiple graphs into low dimensional space. We demonstrate through theory and experiments that the joint embedding algorithm produces features which lead to state of the art performance for subsequent inference tasks on graphs.

There exist a few unsupervised approaches to extract features from graphs. First, classical Principal Component Analysis can be applied by treating each edge of a graph as a raw feature[5]. This approach produces features which are linear combinations of edges, but it ignores the topological structure of graphs and the features extracted are not easily interpretable. Second, features can be extracted by computing summary topological and label statistics from graphs[6, 7]. These statistics commonly include number of edges, number of triangles, average clustering coefficient, maximum effective eccentricity, etc. In general, it is hard to know what intrinsic statistics to compute a priori and computing some statistics can be computationally expensive. Third, many frequent subgraph mining algorithms are developed [8]. For example, the fast frequent subgraph mining algorithm can identify all connected subgraphs that occur in a large fraction of graphs in a graph data set [9]. Finally, spectral feature selection can also be applied to graphs. It treats each graph as a node and constructs an object graph based on a similarity measure. Features are computed through the spectral decomposition of this object graph [10].

Adjacency Spectral Embedding (ASE) and Laplacian Eigenmap (LE) are proposed to embed a single graph observation [11, 12]. The inference task considered in these papers is learning of the block structure of the graph or clustering vertices. Given a set of graphs $\{G_{i}=(V_{i},E_{i})\}_{i=1}^{m}$, ASE and LE need to embed an adjacency matrix or Laplacian matrix of $G_{i}$ individually, and there is no easy way to combine multiple embeddings. The joint embedding method considers the set of graphs together. It takes a matrix factorization approach to extract features for multiple graphs. The algorithm manages to simultaneously identify a set of rank one matrices and project adjacency matrices into the linear subspace spanned by this set of matrices. The joint embedding can be understood as a generalization of ASE for multiple graphs. We demonstrate through simulation experiments that the joint embedding algorithm extracts features which lead to good performance for a variety of inference tasks. In the next section, we review some random graph models and present a model for generating multiple random graphs. In Section 3, we define the joint embedding of graphs and present an algorithm to compute it. In Section 4, we perform some theoretical analyses of our joint embedding. The theoretical results and real data experiments are explored in Section 5. We conclude the paper with a brief discussion of implications and possible future work.

2 Setting

We focus on embedding unweighted and undirected graphs for simplicity, although the joint embedding algorithm works on weighted graphs, and directed graphs with some modifications. Let $\{G_{i}=(V_{i},E_{i})\}_{i=1}^{m}$ be $m$ graphs, each with $n$ vertices, and $\mathbf{A}_{i}$ be the adjacency matrix of graph $G_{i}$. The vertices in these graphs should be matched, which means that all the graphs have a common vertex set $V$. The joint embedding algorithm embeds all $G_{i}$s simultaneously into $\mathbb{R}^{d}$ and represents $G_{i}$ by a vector $\lambda_{i}\in\mathbb{R}^{d}$. Before discussing the joint embedding algorithm, we need a random graph model on multiple graphs, on which the theoretical analysis is based. Let us first recall a model on a single graph: Random Dot Product Graph [13].

• Definition

Random Dot Product Graph (RDPG). Let $\mathcal{X}$ be a subset of $\mathbb{R}^{d}$ such that $x^{T}y\in[0,1]$ for all $x,y\in\mathcal{X}$. Let $\mathbf{X}=[x_{1}^{T},x_{2}^{T},...,x_{n}^{T}]\in\mathcal{X}^{n}$ be a $n\times d$ matrix, and given $\mathbf{X}$, suppose that $\mathbf{A}$ is a random $n\times n$ adjacency matrix such that

[TABLE]

Alternatively,

[TABLE]

Also, define $\mathbf{P}:=\mathbf{X}\mathbf{X}^{T}$ to be edge probability matrix. We write $\mathbf{A}\sim\operatorname{RDPG}(\mathbf{X})$ to denote the distribution of a random dot product graph with latent positions $\mathbf{X}$. When the rows of $\mathbf{X}$ are not fixed, but instead are random variables with a distribution $F$ on $\mathcal{X}$, $(\mathbf{X},\mathbf{A})\sim\operatorname{RDPG}(F)$ denotes the distribution of a random dot product graph with latent positions distributed according to $F$.

The RDPG is a convenient model which is designed to capture the relationship between the vertices of a graph using latent positions, and some extensions have been proposed to capture more general connectivity structures [14]. Moreover, other popular models, including the stochastic block model (SBM) [15] and mixed membership SBM [16] are special cases of the RDPG and its generalizations.. The RDPG can be further generalized to a Latent Position Graph by replacing the inner product by a kernel [17]. The Adjacency Spectral Embedding of a RDPG adjacency matrix is well studied [18]. Next, we propose a new random graph model which generalizes the RDPG to multiple graphs.

• Definition

Multiple Random Eigen Graphs (MREG). Let $h_{1},\ldots,h_{d}$ be vectors in $\mathbb{R}^{n}$ with $\|h_{i}\|_{2}=1,i=1,\ldots,d$, and denote by $\mathcal{X}\subset\mathbb{R}^{d}$ the set of vectors satisfying $\sum\limits_{k=1}^{d}\lambda[k]h_{k}h_{k}^{T}\ \in[0,1]^{n\times n}$ for all $\lambda\in\mathcal{X}$, where $\lambda[k]$ is the $k$th entry of vector $\lambda$. The random adjacency matrices $\mathbf{A}_{1},\ldots,\mathbf{A}_{m}$ follow a $d$-dimensional multiple random eigen graphs model, denoted by

[TABLE]

if the entries of $A_{i}$ are independent Bernoulli random variables,

[TABLE]

$\mathbf{P}_{i}:=\sum_{k=1}^{d}\lambda_{i}[k]h_{k}h_{k}^{T}$ is defined to be the edge probability matrix for graph $i$. In cases that $\{\lambda_{i}\}_{i=1}^{m}$ are of primary interest, they are treated as parameters. When the $\{\lambda_{i}\}_{i=1}^{m}$ are random variables, $F$ denotes their distribution on $\mathcal{X}$, and the model is expressed as

[TABLE]

Compared to the RDPG model, MREG is designed to model multiple graphs. The vectors $\{h_{k}\}_{k=1}^{d}$ are shared across graphs and represent joint latent positions of the vertices, where $(h_{j}[1],\ldots,h_{j}[d])\in\mathbb{R}^{d}$ represents the position of vertex $j$; a $\lambda_{i}$ represents the parameter of an individual graph relative to the latent positions $\{h_{k}\}_{k=1}^{d}$. Note that for a single graph, the edge probability matrix can be written as $\mathbf{P}_{i}=\mathbf{H}\mathbf{\Lambda}_{i}\mathbf{H}^{T}$, where $\mathbf{H}=[h_{1}\cdots h_{d}]\in\mathbb{R}^{n\times d}$ and $\mathbf{\Lambda}_{i}=\text{diag}(\lambda_{i})$. When the entries of $\lambda_{i}$ are non-negative, then $\mathbf{X}_{i}=\mathbf{H}\mathbf{\Lambda}_{i}^{1/2}$ are the latent positions of graph $i$, and hence, on a single graph, RDPG and MREG are equivalent if the edge probability matrix is positive semidefinite. In MREG, we allow self loops to happen. This is mainly for theoretical convenience. Next, we introduce another random graph model: Stochastic Block Model [15], which generalizes the Erdos-Renyi (ER) model [19] that corresponds to a SBM with only one block. SBM is a widely used model to study the community structure of a graph [20, 21].

• Definition

Stochastic Block Model (SBM). Let $\pi$ be a prior probability vector for block membership which lies in the unit $K-1$-simplex. Denote by $\tau=(\tau_{1},\tau_{2},...,\tau_{n})\in[K]^{n}$ the block membership vector, where $\tau$ is a multinomial sequence with probability vector $\pi$. Denote by $\mathbf{B}\in[0,1]^{K\times K}$ the block connectivity probability matrix. Suppose $\mathbf{A}$ is a random adjacency matrix given by,

[TABLE]

Then, $\mathbf{A}$ is an adjacency matrix of a $K$-block stochastic block model graph, and the notation is $\mathbf{A}\sim SBM(\pi,\mathbf{B})$. Sometimes, $\tau$ may also be treated as the parameter of interest, in this case the notation is $\mathbf{A}\sim SBM(\tau,\mathbf{B})$.

The top panel of Figure 1 shows the relationships between three random graph models defined above and the ER model on $1$ graph. The models considered are those conditioned on latent positions, that is $\tau$, $\mathbf{X}$ and $\lambda$ in SBM, RDPG and MREG respectively are treated as parameters; furthermore, loops are ignored in MREG. If an adjacency matrix $\mathbf{A}\sim SBM(\tau,\mathbf{B})$ and the block connectivity matrix $\mathbf{B}$ is positive semidefinite, $\mathbf{A}$ can also be written as an $RDPG(\mathbf{X})$ with $\mathbf{X}$ having at most $K$ distinct rows. If an adjacency matrix $\mathbf{A}\sim RDPG(\mathbf{X})$, then it is also a $1$-graph $MREG(\lambda_{1},h_{1},...,h_{d})$ with $h_{k}$ being the normalized $k$th column of $\mathbf{X}$ and $\lambda_{1}$ being the vector containing the squared norms of columns of $\mathbf{X}$. However, a $1$-graph $MREG(\lambda_{1},h_{1},...,h_{d})$ is not necessarily an RDPG graph since $\lambda_{1}$ could contain negative entries which may result in an indefinite edge probability matrix.

The bottom panel of Figure 1 shows the relationships between the models on multiple graphs. For RDPG, the graphs are sampled i.i.d. with the same parameters. MREG has the flexibility to have $\lambda$ differ across graphs, which leads to a more generalized model for multiple graphs. A $d$ dimensional MREG can represent any SBM with $K\leq d$ blocks, in which the block memberships are the same across all the graphs in the population, but possible different connectivity matrices $\mathbf{B}_{i}$ for each graph, a common assumption in modeling multilayer and time-varying networks [22, 23]. Other models in this setting also allow some vertices of the SBM to change their block memberships over time [22, 24]; in this case, the MREG can still represent those models, but the dimension $d$ may need to increase. Actually, it turns out that if $d$ is allowed to be as large as $\frac{n(n-1)}{2}$, MREG can represent any distribution on binary graphs, which includes distributions in which edges are not independent.

Theorem 2.1

Given any distribution $\mathcal{F}$ on graphs and a random adjacency matrix $\mathbf{A}\sim\mathcal{F}$, there exists a dimension $d$, a distribution $F$ on $\mathbb{R}^{d}$, and a set of vectors $\{h_{k}\}_{k=1}^{d}$, such that $\mathbf{A}\sim MREG(F,h_{1},...,h_{d})$.

Theorem 2.1 implies that MREG is really a semi-parametric model, which can capture any distribution on graphs. One can model any set of graphs by MREG with the guarantee that the true distribution is in the model with $d$ being large enough. However, in practice, a smaller $d$ may lead to better inference performance due to reduction in the dimensionality.

In the next section, we consider the joint embedding algorithm which can be understood as a parameter estimation procedure for MREG.

3 Joint Embedding

3.1 Joint Embedding of Graphs

The joint embedding method considers a collection of vertex-aligned graphs, and estimates a common embedding space across all graphs and a loading for each graph. Specifically, it simultaneously identifies a subspace spanned by a set of rank one symmetric matrices and projects each adjacency matrix $\mathbf{A}_{i}$ into the subspace. The coefficients obtained by projecting $\mathbf{A}_{i}$ are denoted by $\hat{\lambda}_{i}\in\mathbb{R}^{d}$, which is called the loading for graph $i$. To estimate rank one symmetric matrices and loadings for graphs, the algorithm minimizes the sum of squared Frobenius distances between adjacency matrices and their projections as described below.

• Definition

Joint Embedding of Graphs (JE). Given $m$ graphs $\{G_{i}\}_{i=1}^{m}$ with $\mathbf{A}_{i}$ being the corresponding adjacency matrix, the $d$-dimensional joint embedding of graphs $\{G_{i}\}_{i=1}^{m}$ is given by

[TABLE]

Here, $\|\cdot\|$ denotes the Frobenius norm and $\lambda_{i}[k]$ is the $k$th entry of vector $\lambda_{i}$.

To make sure that the model is identifiable and avoid the scaling factor, $h_{k}$ is required to have norm $1$. In addition, $\{h_{k}h_{k}^{T}\}_{k=1}^{d}$ must be linearly independent to avoid identifiability issues in estimating $\lambda_{i}$; however, $\{h_{k}\}_{k=1}^{d}$ needs not to be linearly independent or orthogonal. To ease the notations, let us introduce two matrices $\mathbf{\Lambda}\in\mathbb{R}^{m\times d}$ and $\mathbf{H}\in\mathbb{R}^{n\times d}$, where $\lambda_{i}$ is the $i$th row of $\Lambda$ and $h_{k}$ is the $k$th row of $\mathbf{H}$; that is, $\mathbf{\Lambda}=[\lambda_{1}^{T},...,\lambda_{m}^{T}]$ and $\mathbf{H}=[h_{1},...,h_{d}]$. The equation (1) can be rewritten using $\mathbf{\Lambda}$ and $\mathbf{H}$ as

[TABLE]

Denote the function on the left hand side of the equation by $f(\mathbf{\Lambda},\mathbf{H})$ which is explicitly a function of $\lambda_{i}$s and $h_{k}$s. There are several alternative ways to formulate the problem. If vector $\lambda_{i}$ is converted into a diagonal matrix $\mathbf{D}_{i}\in\mathbb{R}^{d\times d}$ by putting entries of $\lambda_{i}$ on the diagonal of $\mathbf{D}_{i}$, then solving equation (1) is equivalent to solving

[TABLE]

Equation (1) can be also viewed as a tensor factorization problem. If $\{\mathbf{A}_{i}\}_{i=1}^{m}$ are stacked in a 3-D array ${\mathbb{A}}\in\mathbb{R}^{m\times n\times n}$, then solving equation (1) is also equivalent to

[TABLE]

where $\otimes$ denotes the tensor product and $\mathbf{\Lambda}_{*k}$ is the $k$th column of $\mathbf{\Lambda}$. It is well known in the tensor factorization community that the solution to Equation 1 may not necessarily exist for $d\geq 2$. This phenomenon was first found by Bini et al. [25], and Silva and Lim gives a characterization all such tensors in the order-$3$ rank-$2$ case [26]. Although there may not exist a global minimum, finding the local solution in a compact region still provide significant insights to the data. We design an algorithm which is guaranteed to converge, and provide analysis under the $d=1$ case.

The joint embedding algorithm assumes the graphs are vertex-aligned and undirected. The vertex-aligned graphs are common in many applications of interest such as neuroimaging [27], multilayer networks [28] or time-varying graphs [7]. In case that the graphs are not aligned, graph matching should be performed before the joint embedding [29, 30]. The mis-alignments of some vertices will have adverse effects in estimating corresponding latent positions in $\mathbf{H}$; however, a small number of mis-aligned vertices should not have a big impact in estimating $\mathbf{\Lambda}$. If the graphs have weighted edges, the joint embedding can still be applied. Also, the MREG model can be easily extended to weighted graphs by replacing the Bernoulli distribution with other proper distributions. In fact, in the experiment of section 5.3, the graphs are weighted, where the edge weights are the log of fiber counts across regions of brains. In case of directed graphs, to apply the joint embedding, one can symmetrize the graph by removing the direction of edges. Alternatively, $h_{k}h_{k}^{T}$ in equation (1) can be replaced by $h_{k}g_{k}^{T}$, with $h_{k}$ and $g_{k}$ representing the in and out latent positions respectively. With this modification, equation (1) becomes the tensor factorization problem [31].

The optimization problem in equation (1) is similar to Principal Component Analysis (PCA) in the sense of minimizing squared reconstruction error to recover loadings and components [5]. However, there are extra symmetries and rank constraints on the components. Specifically, if $h_{k}h_{k}^{T}$ is replaced by a matrix $\mathbf{S}_{k}$ in equation (1)

[TABLE]

the problem can be solved by applying PCA on vectorized adjacency matrices. In this case, there is a $\mathbf{S}_{k}$ to estimate for each latent dimension which has $\frac{n(n+1)}{2}$ parameters. Compared to PCA, the joint embedding estimates a rank one matrix $h_{k}h_{k}^{T}$ for each latent dimension which has $n$ parameters, and $h_{k}$ can be treated as latent positions for vertices, but the joint embedding yields a larger approximation error due to the extra constraints. Similar optimization problems have also been considered in the simultaneous diagonalization literature [32, 33]. The difference is that the joint embedding is estimating an $n$-by-$d$ matrix $\mathbf{H}$ by minimizing reconstruction error instead of finding a $n$-by-$n$ non-singular matrix by trying to simultaneously diagonalize all matrices. The problem in equation (1) has considerably fewer parameters to optimize, which makes it more stable and applicable with $n$ being moderately large. In case of embedding only one graph, the joint embedding is equivalent to the Adjacency Spectral Embedding solved by singular value decomposition [11]. Next, we describe an algorithm to optimize the objective function $f(\mathbf{\Lambda},\mathbf{H})$.

3.2 Alternating Descent Algorithm

The joint embedding of $\{G_{i}\}_{i=1}^{m}$ is estimated by solving the optimization problem in equation (1). There are a few methods proposed to solve similar problems. Alternating least squares (ALS) is a popular method to solve similar problems [34, 31], but often ignore symmetric constraints. Gradient approaches have also been considered for similar problems [35, 36]. We develop an alternating descent algorithm to minimize $f(\mathbf{\Lambda},\mathbf{H})$ that combines ideas from both approaches [37]. The algorithm can also be understood as a block coordinate descent method with $\mathbf{\Lambda}$ and $\mathbf{H}$ being the two blocks [38, 39]. The algorithm iteratively updates one of $\mathbf{\Lambda}$ and $\mathbf{H}$ while treating the other parameter as fixed. Optimizing $\mathbf{\Lambda}$ when fixing $H$ is straightforward, since it is essentially a least squares problem. However, optimizing $\mathbf{H}$ when fixing $\mathbf{\Lambda}$ is hard due to the fact that the problem is non-convex and there is no closed form solution available. In this case, the joint embedding algorithm utilizes gradient information and takes an Armijo backtracking line search strategy to update $\mathbf{H}$ [40].

Instead of optimizing all columns $\mathbf{\Lambda}$ and $\mathbf{H}$ simultaneously, we consider a greedy algorithm which solves the optimization problem by only considering one column of $\mathbf{H}$ at a time. Specifically, the algorithm fixes all estimates for the first $k_{0}-1$ columns of $\mathbf{\Lambda}$ and $\mathbf{H}$ at iteration $k_{0}$, and then the objective function is minimized by searching through only the $k_{0}$th column of $\mathbf{\Lambda}$ and $\mathbf{H}$. That is,

[TABLE]

When $m=1$, by the Eckart-Young theorem, the solution for $(\hat{\mathbf{\Lambda}}_{1,k_{0}},\hat{h}_{k_{0}})$ is given by the leading eigenvalue and eigenvector of $\mathbf{A}_{1}-\sum\limits_{k=1}^{k_{0}-1}\hat{\mathbf{\Lambda}}_{1k}\hat{h}_{k}\hat{h}_{k}^{T}$ , and hence the solution obtained by our greedy approach coincides with the eigendecomposition of $\mathbf{A}_{1}$. If $m>1$ there is no closed-form solution in general, so we develop an optimization method next.

Let $f(\mathbf{\Lambda}_{*k_{0}},h_{k_{0}})$ denote the sum on the left hand side of the equation. To compute a $d$-dimensional joint embedding $(\hat{\mathbf{\Lambda}},\hat{\mathbf{H}})$, the algorithm iteratively finds one component at a time by solving the optimization for the one dimensional embedding problem (3). Once a new component $k_{0}$ is found, the matrix of loadings $\hat{\mathbf{\Lambda}}$ is updated by minimizing the optimization problem (1) with the first $k_{0}$ components fixed at $\hat{h}_{1},\ldots,\hat{h}_{k_{0}}$, resulting in the least squares problem

[TABLE]

This step finds the projection of the graphs into the space spanned by the first $k_{0}$ components $\{\hat{h}_{k}\hat{h}_{k}^{T}\}_{k=1}^{k_{0}}$. The algorithm solves the $d$-dimensional embedding problem by iteratively updating $\hat{\mathbf{H}}$ and $\hat{\mathbf{\Lambda}}$ using equations (3) and (4).

There are a few advantages in iteratively solving one dimensional embedding problems. First, there are fewer parameters to fit at each iteration, since the algorithm is only allowed to vary $h_{k_{0}}$ at iteration $k_{0}$. This makes initialization and optimization steps much easier compared to optimizing all columns of $\mathbf{H}$ simultaneously. Second, it implicitly enforces an ordering on the columns of $\mathbf{H}$. This ordering allows us to select the top few columns of $\mathbf{\Lambda}$ and $\mathbf{H}$ in cases where model selection is needed after the joint embedding. Third, it allows incremental computation. If $d$ and $d^{\prime}$ dimensional joint embeddings are both computed, the first $\min(d,d^{\prime})$ columns of $\hat{\mathbf{H}}$ will be the same. Fourth, although the global rank-$d$ tensor approximation problem does not always have a solution for $d\geq 2$, a one dimensional embedding that corresponds to a rank-1 tensor approximation always have a solution [26]. Finally, based on numerical experiments, the difference between optimizing iteratively and optimizing all the parameters when $d$ is small is negligible; however, the iterative algorithm yields a slightly smaller objective function when $d$ is large. The disadvantage of optimizing each column separately is that the algorithm is more likely to end up at a local minimum when the objective function is structured not in favor of embedding iteratively. In practice, this problem can be mitigated by running the joint embedding algorithm several times with random initializations.

To find $\mathbf{\Lambda}_{*k_{0}}$ and $h_{k_{0}}$ in equation (3), the algorithm needs to evaluate two derivatives: $\frac{\partial f}{\partial h_{k_{0}}}$ and $\frac{\partial f}{\partial\mathbf{\Lambda}_{ik_{0}}}$. Denote by $\mathbf{R}_{ik_{0}}$ the residual matrix after iteration $k_{0}-1$ which is $\mathbf{A}_{i}-\sum\limits_{k=1}^{k_{0}-1}\hat{\mathbf{\Lambda}}_{ik}\hat{h}_{k}\hat{h}_{k}^{T}$. The gradient of the objective function with respect to $h_{k_{0}}$ is given by

[TABLE]

The derivative of the objective function with respect to $\mathbf{\Lambda}_{ik_{0}}$ is given by

[TABLE]

Setting the derivative to [math] yields

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the inner product.

Once a new component $\hat{h}_{k_{0}}\hat{h}_{k_{0}}^{T}$ is obtained, the algorithm proceeds to update $\mathbf{\Lambda}$ by solving the optimization problem (4). Note that the problem can be split into $m$ least square subproblems, one for each graph $\mathbf{A}_{i}$, of the form

[TABLE]

for each $i=1,\ldots,m$. By obtaining the derivative with respect to $\mathbf{\Lambda}$ for all $i=1,\ldots,m$ and setting them equal to zero, $\mathbf{\Lambda}$ is updated by solving $m$ systems of linear equations with $k_{0}$ variables each. Let $\mathbf{\Gamma}$ be a $k_{0}\times k_{0}$ matrix and $\mathbf{\Psi}$ be a ${k_{0}\times m}$ matrix such that $\mathbf{\Gamma}_{kj}=(\hat{h}_{j}^{T}\hat{h}_{k})^{2}$ and $\mathbf{\Psi}_{ki}=(\hat{h}_{k}^{T}\mathbf{A}_{i}\hat{h}_{k})$. The algorithm updates the first $k_{0}$ columns of $\mathbf{\Lambda}$ by solving

[TABLE]

The joint embedding algorithm alternates between updating $\hat{\mathbf{\Lambda}}_{*k_{0}}$ and $\hat{h}_{k_{0}}$ according to equation (5) and (6), after which the values of all the loadings $(\hat{\mathbf{\Lambda}}_{*1},\ldots,\hat{\mathbf{\Lambda}}_{*k_{0}})$ are updated using equation (7). Algorithm 1 describes the general procedure to compute the $d$-dimensional joint embedding of graphs $\{G_{i}\}_{i=1}^{m}$. The algorithm outputs two matrices: $\hat{\mathbf{\Lambda}}$ and $\hat{\mathbf{H}}$. The rows of $\hat{\mathbf{\Lambda}}$ denoted by $\{\hat{\lambda}_{i}\}_{i=1}^{m}$ can be treated as estimates of $\{\lambda_{i}\}_{i=1}^{m}$ in MREG and features for graphs. Columns of $\hat{\mathbf{H}}$ denoted by $\{\hat{h}_{k}\}_{k=1}^{d}$ are estimates of $\{h_{k}\}_{k=1}^{d}$. If a new graph $G$ is observed with adjacency matrix $\mathbf{A}$, $\mathbf{A}$ can be projected into the linear space spanned by $\{\hat{h}_{k}\hat{h}_{k}^{T}\}_{k=1}^{d}$ to obtain features for the graph.

In case of $\mathbf{A}_{i}$ being large, the updating equations (5) and (6) are not practical due to $h_{k}h_{k}^{T}$ and $\mathbf{R}_{ik}$ being large and dense. However, they can be rearranged to avoid explicit computation of $h_{k}h_{k}^{T}$ and $\mathbf{R}_{ik}$. Equation (5) becomes

[TABLE]

Similarly, equation (6) can be rewritten as

[TABLE]

Based on the rearranged equations, efficiently evaluating matrix vector product $\mathbf{A}_{i}h_{k_{0}}$ is needed to calculate the derivatives. This can be completed for a variety of matrices, in particular, sparse matrices [41].

The Algorithm 1 is guaranteed to converge to a stationary point. Specifically, at the termination of iteration $k_{0}$, $\frac{\partial f}{\partial h_{k_{0}}}\approx 0$ and $\frac{\partial f}{\partial\mathbf{\Lambda}_{ik_{0}}}\approx 0$. First, $\frac{\partial f}{\partial\mathbf{\Lambda}_{ik_{0}}}\approx 0$ is ensured due to exact updating by equation (6). Second notice that updating according to equation (5) and (6) always decreases the objective function. Due to the fact that $h_{k_{0}}$ lies on the unit sphere and the objective is twice continuous differentiable, $\frac{\partial f}{\partial h_{k_{0}}}$ is Lipschitz continuous. This along with Armijo backtracking line search guarantees a ”sufficient” decrease $c\|\frac{\partial f}{\partial h_{k_{0}}}\|^{2}$ in objective function each time when the algorithm updates $h_{k_{0}}$ [40], where $c$ is a constant independent of $h_{k_{0}}$. Since the objective function is bounded below by [math], this implies convergence of gradient, that is $\frac{\partial f}{\partial h_{k_{0}}}\rightarrow 0$.

The joint embedding of graphs model requires for identifiability that the basis of the embedding space $\{h_{k}h_{k}^{T}\}_{k=1}^{m}$ is linearly independent. This condition also ensures that equation (4) has a unique solution. Although the optimization problem (3) does not directly enforce this constraint, the next theorem guarantees that the components obtained by our method are linearly independent. The proof is given in the appendix.

Theorem 3.1

Suppose that $\hat{\mathbf{\Lambda}}=(\hat{\mathbf{\Lambda}}_{*1},...,\hat{\mathbf{\Lambda}}_{*k_{0}})$ and $\hat{\mathbf{H}}=[\hat{h}_{1},...,\hat{h}_{k_{0}}]$ are obtained by iteratively solving (3) and (4) for $k_{0}$ components, and the first $k_{0}-1$ components $\{\hat{h}_{k}\hat{h}_{k}^{T}\}_{k=1}^{k_{0}-1}$ are linearly independent. Then $\{\hat{h}_{k}\hat{h}_{k}^{T}\}_{k=1}^{k_{0}}$ are also linearly independent, or $\hat{\mathbf{\Lambda}}_{ik_{0}}=0$ for all $i=1,\ldots,m$.

In general, the objective function may have multiple stationary points due to non-convexity. Therefore, the joint embedding algorithm is sensitive to initializations. Actually, like many of the problems in tensor factorization, finding the global minimum in joint embedding is NP-Hard [42]. When time permits, we recommend running the joint embedding several times including many random initializations. In Section 5.1, we study the effects of different initialization approaches through a numerical simulation experiment. For other simulation and real experiments, we initialize $\hat{h}_{k_{0}}$ using the top left singular vector of the average residual matrix $\sum\mathbf{R}_{ik_{0}}/m$. Computing this vector is computationally cheap compared to the computational cost of the joint embedding when the number of graphs $m$ is larger than 1. When $m=1$, the SVD initialization and the joint embedding solution coincide, and when $m>1$ the SVD initialization gives a good approximation to the solution if all the graphs are identically distributed. The optimization algorithm described above may not be the fastest approach to solving the problem; however, numerical optimization is not the focus of this paper. Based on results from numerical applications, our approach works well in estimating parameters and extracting features for subsequent statistical inference. Next, we discuss some variations of the joint embedding algorithm.

3.3 Variations

The joint embedding algorithm described in the previous section can be modified to accommodate several different settings.

Variation 1. When all graphs come from the same distribution, we can force estimated loadings $\hat{\lambda}_{i}$ to be equal across all graphs. This is useful when the primary inference task is to extract features for vertices. Since all graphs share the same loadings, with slightly abusing notations, let $\mathbf{\Lambda}$ be a vector in $\mathbb{R}^{d}$ and the optimization problem becomes

[TABLE]

which is equivalent to

[TABLE]

Therefore, the optimization problem can be solved exactly by finding the singular value decomposition of the average adjacency matrix $\frac{1}{m}\sum\limits_{i=1}^{m}\mathbf{A}_{i}$.

Variation 2. When there is a discrete label $y_{i}\in\mathbb{Y}$ associated with the graph $G_{i}$ available, we may require all loadings $\hat{\lambda}_{i}$ to be equal within class. Let $\mathbf{\Lambda}\in\mathbb{R}^{|\mathbb{Y}|\times d}$, the optimization problem becomes

[TABLE]

In this case, when updating $\mathbf{\Lambda}$ as in equation (6), the algorithm should average $\mathbf{\Lambda}_{yk}$ within the same class, that is

[TABLE]

Variation 3. In some applications, we may require all $\mathbf{\Lambda}_{ik}$ to be greater than [math], as in non-negative matrix factorization. One advantage of this constraint is that graph $G_{i}$ may be automatically clustered based on the largest entry of $\hat{\lambda}_{i}$. In this case, the optimization problem is

[TABLE]

To guarantee nonnegativity, the algorithm should use nonnegative least squares in updating $\mathbf{\Lambda}$ [43]. Furthermore, a constraint on the number of non-zero elements in $i$th row of $\mathbf{\Lambda}$ can be added as in K-SVD [44], and a basis pursuit algorithm could be used to update $\mathbf{\Lambda}$ [45, 46]. Next, we discuss some theoretical properties of the MREG model and joint embedding when treated as a parameter estimation procedure for the model.

4 Theoretical Results

In this section, we consider a simple setting where graphs follow a $1$-dimensional MREG model, that is $\{(\lambda_{i},\mathbf{A}_{i})\}_{i=1}^{m}\sim MREG(F,h_{1})$. The $1$-dimensional joint embedding is well defined in this case, that is $\hat{\lambda}_{i}$ and $\hat{h}_{1}$ defined in Equation 1 is guaranteed to exist. Under this MREG model, the joint embedding of graphs can be understood as estimators for parameters of the model. Specifically, $\hat{\lambda}_{i}$ and $\hat{h}_{1}$ are estimates of $\lambda_{i}$ and $h$. We prove two theorems concerning the asymptotic behavior of estimator $\hat{h}_{1}$ produced by joint embedding.

Let $\hat{h}_{1}^{m}$ denote the estimates based on $m$ graphs and define functions $\rho$, $D_{m}$ and $D$ as below:

[TABLE]

[TABLE]

[TABLE]

One can understand $D_{m}$ and $D$ as sample and population approximation errors respectively. By equation (1),

[TABLE]

By equation (6),

[TABLE]

Therefore,

[TABLE]

The first theorem states that $\hat{h}_{1}^{m}$ converges almost surely to a global minimum of $D(h,h_{1})$. Alternatively, the theorem implies the sample minimizer converges to the population minimizer.

Theorem 4.1

The estimator $\hat{h}_{1}^{m}$ converges almost surely to the set of global minimizers of $D(h,h_{1})$ as $m$ goes to infinity. That is,

[TABLE]

Theorem 4.1 ensures that, in the limit, $\hat{h}^{m}$ is arbitrarily close to parameters in the set of global minimizer of $D(h,h_{1})$. However, the global minimizer is definitely not unique due to the symmetry up to sign flip of $h$, that is $D(h,h_{1})=D(-h,h_{1})$ for any $h$. This problem can be addressed by forcing an orientation of $\hat{h}_{1}^{m}$ or stating that the convergence is up to a sign flip. In this case, Theorem 4.1 does not apply. We are currently only certain that when all graphs are from the Erdos-Renyi random graph model, the global minimizer is unique up to a sign flip. The next theorem concerns the asymptotic bias of $h^{\prime}$. It gives a bound on the difference between the population minimizer $h^{\prime}$ and the truth $h_{1}$.

Theorem 4.2

If $h^{\prime}$ is a minimizer of $D(h,h_{1})$, then

[TABLE]

To see an application of Theorem 4.2, let us consider the case in which all graphs are Erdos-Renyi graphs with $100$ vertices and edge probability of $0.5$. Under this setting, Theorem 4.2 implies $\|h^{\prime}-h_{1}\|\in[0,0.04]\cup[1.28,1.52]$. The second interval is disturbing. It is due to the fact that when $h_{1}^{T}h^{\prime}$ is small, the bound is useless. We provide some insights as to why the second interval is there and how we can get rid of it with additional assumptions. In the proof of Theorem 4.2, we show that the global optimizer $h^{\prime}$ satisfies

[TABLE]

Taking a closer look at $E(\langle\mathbf{A}_{i},hh^{T}\rangle^{2})$,

[TABLE]

Therefore,

[TABLE]

We can see that $\mathbb{E}(\lambda_{i}^{2})(h_{1}^{T}h)^{4}$ is maximized when $h=h_{1}$; however, the noise term $\mathbb{E}((h^{T}(\mathbf{A}_{i}-\mathbf{P}_{i})h)^{2})$ is generally not maximized at $h=h_{1}$. If $n$ is large, we can apply a concentration inequality to $(h^{T}(\mathbf{A}_{i}-\mathbf{P}_{i})h)^{2}$ and have an upper bound on $\mathbb{E}((h^{T}(\mathbf{A}_{i}-\mathbf{P}_{i})h)^{2})$. If we further assume $A_{i}$ is not too sparse, that is $\mathbb{E}(\lambda_{i}^{2})$ grows with $n$ fast enough, then the sum of these two terms is dominated by the first term. This provides a way to have a lower bound on $h_{1}^{T}h^{\prime}$. We may then replace the denominator of the bound in Theorem 4.2 by the lower bound. In general, if $n$ is small, the noise term may cause $h^{\prime}$ to differ from $h_{1}$ by a significant amount. In this chapter, we focus on the case that $n$ is fixed. The case that $n$ goes to infinity for Random Dot Product Graph is considered in [47].

The two theorems above concern only the estimation of $h_{1}$, but not $\lambda_{i}$. Based on equation (6), the joint embedding estimates $\lambda_{i}$ by

[TABLE]

When $m$ goes to infinity, we can apply Theorem 4.1,

[TABLE]

Then, applying the bound on $\|h^{\prime}-h_{1}\|$ derived in Theorem 4.2 and utilizing the fact that $h^{T}\mathbf{A}_{i}h$ is continuous in $h$, we can obtain an upper bound on $|\hat{\lambda}_{i}^{m}-h_{1}^{T}\mathbf{A}_{i}h_{1}|$. When $\mathbf{A}_{i}$ is large, $h_{1}^{T}\mathbf{A}_{i}h_{1}$ is concentrated around $\lambda_{i}$ with high probability. As a consequence, with high probability $|\hat{\lambda}_{i}^{m}-\lambda_{1}|$ is small. In the next section, we demonstrate properties and utilities of the joint embedding algorithm through experiments.

5 Experiments

Before going into details of our experiments, we want to discuss how to select the dimensionality $d$ of the joint embedding. Estimating $d$ is an important model selection question that has been studied for years under various settings [48]. Model selection is not the focus of this paper, but we still face this problem in numerical experiments. In the simulation experiments of this section, we assume $d$ is known to us and simply set the dimensionality estimate $\hat{d}$ equal to $d$. In the real data experiment, we recommend two approaches to determine $\hat{d}$. Both approaches require first running the $d^{\prime}$-dimensional joint embedding algorithm, where $d^{\prime}$ is sufficiently large. We then plot the objective function versus dimension, and determine $\hat{d}$ to be where the objective starts to flatten out. Alternatively, we can plot $\{\hat{\mathbf{\Lambda}}_{ik}\}_{i=1}^{m}$ for $k=1,...,d^{\prime}$, and select $\hat{d}$ when the loadings start to look like noise with [math] mean. These two approaches should yield a similar dimensionality estimate of $\hat{d}$.

5.1 Simulation Experiment 1: Joint Embedding Under a Simple Model

In the first experiment, we present a simple numerical example to demonstrate some properties of the joint embedding procedure as the number of graphs grows. We repeatedly generate graphs with $20$ vertices from $3$-dimensional MREG, where $\lambda_{i}[1]\sim\text{Uniform}(8,16)$, $\lambda_{i}[2]\sim\text{Uniform}(0,2)$ and $\lambda_{i}[3]\sim\text{Uniform}(0,1)$, with

[TABLE]

We keep doubling the number of graphs $m$ from $2^{4}$ to $2^{12}$. At each value of $m$, we compute the $3$-dimensional joint embedding of graphs. Let the estimated parameters based on $m$ graphs be denoted by $\hat{\lambda}_{i}^{m}$ and $\hat{h}_{k}^{m}$. Two quantities based on $\hat{h}_{k}^{m}$ are calculated. The first is the norm difference between the current $h_{k}$ estimates and the previous estimates, namely $\|\hat{h}_{k}^{m}-\hat{h}_{k}^{m/2}\|$. This provides numerical evidence for the convergence of our principled estimation procedure. The second quantity is $\|\hat{h}^{m}_{k}-h_{k}\|$. This investigates whether $\hat{h}_{k}$ is an unbiased estimator for $h_{k}$. The procedure described above is repeated $20$ times. Figure 2 presents the result.

Based on the plot, the norm of differences $\|\hat{h}^{m}_{k}-\hat{h}_{k}^{m/2}\|$ seem to converge to [math] as $m$ increases. This suggests the convergence of $\hat{h}_{1}^{m}$. Second, we notice that the bias $\|\hat{h}^{m}_{2}-h_{2}\|$ and $\|\hat{h}^{m}_{3}-h_{3}\|$ do not converge to [math]; instead, it stops decreasing at around $0.1$ and $0.2$ respectively. This suggests that $\hat{h}_{k}$ is an asymptotically biased estimator for $h_{k}$. Actually, this is as to be expected: when there are infinitely many nuisance parameters present, Neyman and Scott demonstrate that maximum likelihood estimator is inconsistent [49]. In our case, there are infinitely many $\lambda_{i}$ as $m$ grows; therefore, we do not expect the joint embedding to provide an asymptotic consistent estimate of $h_{k}$.

In applications such as clustering or classifying multiple graphs, we may be not interested in $\hat{h}_{k}$. $\hat{\lambda}_{i}$ is of primary interest, which provides information specifically about the graphs $G_{i}$. Here, we consider two approaches to estimate $\lambda_{i}[1]$. The first approach is estimating $\lambda_{i}[1]$ through joint embedding, that is

[TABLE]

The second approach estimates $\lambda_{i}$ by assuming $h_{1}$ is known. In this case, equation (6) gives

[TABLE]

$\hat{\lambda}_{i}[1]$ calculated this way can be thought as the ’oracle’ estimate. Figure 3 shows the differences in estimates provided by two approaches. Not surprisingly, the differences are small due to the fact that $\hat{h}_{1}^{m}$ and $h_{1}$ are close.

Next, we investigate the effects of four different initialization approaches. The first approach utilizes SVD on average residual matrix to initialize $h_{k}$ at each iteration. The second approach randomly samples independent Gaussian variable for each entry of $h_{k}$. The third approach takes the best initialization among $10$ random initializations. The fourth approach initializes $h_{k}$ using the truth. To compare these approaches, we generate $16$ graphs from the MREG model and jointly embed them with four different initializations. Then, another $16$ graphs are generated and the objective function on these graphs are evaluated using $\hat{\mathbf{H}}$ estimated by joint embedding. This procedure is repeated $100$ times. Mean objective function and total running time with standard error of these four approaches are shown in Table I.

Based on Wilcoxon signed-rank test, SVD and truth initializations are significantly better than random initializations on this scenario. For the rest experiments, the initialization is completed by SVD.

5.2 Simulation Experiment 2: Joint Embedding to Classify Graphs

In this experiment, we consider the inference task of classifying graphs. We have $m$ pairs $\{(\mathbf{A}_{i},y_{i})\}_{i=1}^{m}$ of observations. Each pair consists of an adjacency matrix $\mathbf{A}_{i}\in\{0,1\}^{n\times n}$ and a label $y_{i}\in[K]$. Furthermore, all pairs are assumed to be independent and identically distributed according to an unknown distribution $\mathbb{F}_{\mathbf{A},y}$, that is

[TABLE]

The goal is to find a classifier $g$ which is a function $g:\{0,1\}^{n\times n}\rightarrow[K]$ that has a small classification error $L_{g}=P(g(\mathbf{A})\neq y)$.

We consider a binary classification problem where $y$ takes value $1$ or $2$. $200$ graphs with $100$ vertices are independently generated. The graphs are sampled from a $2$-dimensional MREG model. Let $h_{1}$ and $h_{2}$ be two vectors in $\mathbb{R}^{100}$, and

[TABLE]

Here, $h_{2}$ has $-0.1$ as its first $50$ entries and $0.1$ as its last $50$ entries. Graphs are generated according to the MREG model,

[TABLE]

where $F$ is a mixture of two point masses with equal probability,

[TABLE]

We let the class label $y_{i}$ indicate which point mass $\lambda_{i}$ is sampled from. In terms of SBM, this graph generation scheme is equivalent to

[TABLE]

[TABLE]

To classify graphs, we first jointly embed $200$ graphs. The first two dimensional loadings are shown in Figure 4. We can see two classes are separated after being jointly embedded.

Then, a $1$-nearest neighbor classifier (1-NN) is constructed based on loadings $\{\hat{\lambda}_{i}\}_{i=1}^{m}$.

We compare classification performances of using the joint embedding to extract features to five other feature extraction approaches: Adjacency Spectral Embedding, Laplacian Eigenmap, Graph Statistics, Graph Spectral Statistics, and PCA. For Adjacency Spectral Embedding (ASE) and Laplacian Eigenmap (LE), we first embed each adjacency matrix or normalized Laplacian matrix and then compute the Procrustes distance between embeddings. For Graph Statistics (GS), we compute topological statistics of graphs considered by Park et al. in [7]. For Graph Spectral Statistics (GSS), we compute the eigenvalues of adjacency matrices and treat them as features [50]. For PCA, we vectorize the adjacency matrices and compute the factors through SVD. After the feature extraction step, we also apply a $1$-NN rule to classify graphs. We let the number of graphs $m$ increase from $4$ to $200$. For each value of $m$, we repeat the simulation $100$ times. Figure 5 shows the result. The joint embedding takes advantage of increasing sample size and outperforms other approaches when given more than $10$ graphs.

5.3 Real Data Experiment 1: Subject classification on HNU1 data

In this section, we use the HNU1 data [51] to classify graphs from different subjects. These data consist of diffusion tensor imaging (DTI) records of the brain of 30 healthy different subjects, each of which was scanned 10 times over a period of one month. Each scan was processed with the NeuroData’s MRI to Graphs (NDMG) pipeline [52] using the Talairach atlas [53] to register the vertices, and we obtained a sample of 300 graphs (10 graphs per subject), each of which is composed of 1105 vertices.

It has been suggested that the information encoded in patterns of brain connectivity can uniquely identify different subjects [54, 55], and there is some evidence of low-rank structure in those differences [56, 57]. We study how low-dimensional representations can capture inter-individual variability by using the HNU1 data to classify subject scans. The joint embedding of graphs is an ideal method for this task, since it provides a low-rank representation of each of the graphs which is jointly learned from the sample.

The information encoded in the adjacency matrices can be used to accurately discriminate between different individuals. This can be confirmed by a 1-NN classifier using the vectorized adjacency matrices, which gives an almost perfect classification accuracy. However, our goal here is to study the accuracy of low-dimensional representations of the graphs to discriminate between subjects. Thus, we measure how does our joint dimensional embedding method perform as the number of embedding dimensions grows. We compare the performance of our method with PCA, as another dimensionality reduction method. Note that the joint embedding of graphs directly imposes a rank-one restriction on the components, enforcing a low-rank representation of the graphs. PCA instead represents the adjacency matrices using components that are usually full rank, and requires many more parameters to represent those components.

Before computing the embedding, all graphs are centered by the mean, that is, we compute $\tilde{\mathbf{A}}_{i}=\mathbf{A}_{i}-\frac{1}{300}\sum_{j=1}^{300}\mathbf{A}_{j}$. Using a 5-fold cross-validation, the data of each subject is divided into training and test sets. We measure the effect of the sample size using two different scenarios, the first uses $30$ graphs on the training set (one scan per subject) and the second uses $240$ graphs (eight scans per subject). The rest of the graphs are used in the test set to evaluate the accuracy, and the average over the 5 folds is computed. We follow the same procedure as in experiment 1. Using the training set, we compute an embedding of the graphs into $d$ dimensions, either by doing JEG or PCA on the vectorized adjacency matrices. After that, all the data is projected into the $d$-dimensional embedding, and we use 1-NN to estimate the labels of the test data.

Figure 6 shows a comparison of the average classification accuracy and model complexity for both methods. We measure model complexity as the number of embedding dimensions $d$, and as the total number of parameters used in each instance ($d(n+m)$ for JEG and $d(n^{2}+m)$ for PCA). In general, both methods are able to perform very accurate classification when $d$ is large enough, but JEG uses far fewer parameters in the components for the same number of dimensions (see middle plot). The advantage of JEG is especially remarkable when the sample size is small (first row). In this case, JEG shows a better performance than PCA when $d$ is not large. For values close to $30$, PCA shows almost perfect performance, which is not surprising, since the maximum number of principal components that can be obtained by PCA is $m=30$, and thus PCA is not really performing any dimensionality reduction when $d$ is close to this value, but rather projecting onto the training data itself. JEG on the other hand is able to provide very accurate classification error in all scenarios using low-rank representations of the graphs with a fewer number of parameters that are interpretable. These can be observed in Figure 7, which shows the latent positions of the vertices obtained by JEG for the first seven dimensions. To construct these plots, only the vertices that are labeled as left (507 of them) or right (525) were considered. As it can be observed, several dimensions of the latent positions show a structure that is clearly related to the hemisphere side. This structure is specially highlighted on the 7th dimension of the embedding.

5.4 Real Data Experiment 2: Joint Embedding to Cluster Vertices

In the previous experiments, we focus on feature extraction for graphs. Here, we consider a different task, that is spectral clustering through the joint embedding. In general, spectral clustering first computes (generalized) eigenvalues and eigenvectors of adjacency matrix or Laplacian matrix, then clusters the latent positions into groups [11, 12]. The cluster identities of latent positions become the cluster identities of vertices of the original graph. Adjacency Spectral Embedding (ASE) is one of the spectral clustering algorithms used to find the latent positions of the vertices to fit the SBM and RDPG [58], which are special cases of our MREG model. When applied to one graph, the joint embedding is equivalent to Adjacency Spectral Embedding (ASE). When given multiple graphs, the joint embedding can estimate latent positions for graph $i$ as $[\hat{\lambda}_{i}[1]^{\frac{1}{2}}\hat{h}_{1},\hat{\lambda}_{i}[2]^{\frac{1}{2}}\hat{h}_{2},...,\hat{\lambda}_{i}[d]^{\frac{1}{2}}\hat{h}_{d}]$ or equivalently $\hat{\mathbf{H}}\hat{\mathbf{D}}_{i}^{\frac{1}{2}}$. Alternatively, the matrix of components $\hat{H}=[\hat{h}_{1},\ldots,\hat{h}_{d}]$ gives joint latent positions for all the graphs. Then, a clustering algorithm can be applied to the latent positions to obtain communities.

We apply this spectral clustering approaches to two Wikipedia graphs [59]. The vertices of these graphs represent Wikipedia article pages. Two vertices are connected by an edge if either of the associated pages hyperlinks to the other. Two graphs are constructed based on English webpages and French webpages. The full graph has $1382$ vertices which represents articles within $2$-neighborhood of ”Algebraic Geometry”. Based on the content of the associated articles, they are grouped by hand into $6$ categories: people, places, dates, things, math terms, and categories.

We consider a subset of vertices from $3$ categories: people, things, and math terms. After taking the induced subgraph of these vertices and removing isolated vertices, there are $n=704$ vertices left. Specifically, $326$, $181$, and $197$ vertices are from people, things and math terms respectively. We consider approaches to embed the graphs to obtain latent positions. First, we consider clustering of each graph separately by doing ASE on the English graph $\mathbf{A}_{en}$ (ASE+EN), or equivalently, JE on $\mathbf{A}_{en}$, and ASE on the French Graph $\mathbf{A}_{fr}$ (ASE+FR), and compare with the individual latent positions obtained by JE $\hat{\mathbf{H}}\hat{\mathbf{D}}_{en}^{\frac{1}{2}}$ (JE+EN) and $\hat{\mathbf{H}}\hat{\mathbf{D}}_{fr}^{\frac{1}{2}}$ (JE+FR). We also consider methods to estimate joint latent positions, by doing ASE on the mean of both graphs $\bar{\mathbf{A}}=(\mathbf{A}_{en}+\mathbf{A}_{fr})/2$ (ASE+(EN+FR)) [56], and the matrix $\hat{\mathbf{H}}$ obtained by JE on both graphs (JE+(EN,FR)). The dimension $d$ is set to $3$ for all approaches, and the latent positions are scaled to have norm $1$ for degree correction. Then, we apply $3$-means algorithm to the latent positions.

The joint latent positions of the graphs $\hat{\mathbf{H}}$ estimated by the joint embedding are provided on Figure 8. The latent positions of math terms are separated from the other two clusters. However, the latent positions of people and things are mixed. Table II shows the clustering results measured by adjusted rand index and the purity of clustering [60, 61]. The standard error is estimated through repeatedly clustering bootstrapped latent positions. All methods yield clustering results which are significantly better than random. The English graph demonstrates clearer community structure than the French graph. The joint embedding produces latent positions that combine the information in both graphs, and leads to better clustering performance. The JE and ASE give similar results on the English graph, but JE is able to improve the clustering performance on the French graph significantly. The JE also shows better performance than ASE on $\bar{\mathbf{A}}$, which also uses the information of both graphs. These results show that our method provide interpretable representations for the vertices, with an accuracy comparable to state of the art methods for spectral clustering.

6 Conclusion

In summary, we developed a joint embedding method that can simultaneously embed multiple graphs into low dimensional space. The joint embedding can be utilized to estimate features for inference problems on multiple vertex matched graphs. Learning on multiple graphs has significant applications in diverse fields and our results have both theoretical and practical implications for the problem. As the real data experiment illustrates, the joint embedding is a practically viable inference procedure. We also proposed a Multiple Random Eigen Graphs model. It can be understood as a generalization of the Random Dot Product Graph model or the Stochastic Block Model for multiple random graphs. We analyzed the performance of joint embedding on this model under simple settings. We demonstrated that the joint embedding method provides estimates with bounded error. Our approach is intimately related to other matrix and tensor factorization approaches such as singular value decomposition and CP decomposition. Indeed, the joint embedding and these algorithms all try to estimate a low dimensional representation of high dimensional objects through minimizing a reconstruction error. We are currently investigating the utility of joint embedding with more or less regularizations on parameters and under different set ups. We are optimistic that our method provides a viable tool for analyzing multiple graphs and can contribute to a deeper understanding of the joint structure of networks.

Appendix A

Proof:

Denote the probability of observing a particular adjacency matrix $\mathbf{A}_{i}$ under distribution $\mathcal{F}$ by $p_{i}$. It suffices to show that there is a set of parameters for MREG such that observing $\mathbf{A}_{i}$ under MREG is also $p_{i}$.

For undirected graphs with loops on $n$ vertices, there are ${n+1}\choose{2}$ possible edges. Let $\mathbf{A}_{1},\mathbf{A}_{2},...,\mathbf{A}_{2^{{n+1}\choose{2}}}$ be all the possible adjacency matrices. Since real symmetric matrix of size $n$ has ${n+1}\choose{2}$ free entries which lies in a linear space, if there exists ${n+1}\choose{2}$ linearly independent rank one symmetric matrices, they form a basis for this space. It turns out that the rank one symmetric matrices generated by vectors $\{e_{i}\}_{i=1}^{n}\cup\{e_{i}+e_{j}\}_{i<j}$ are linearly independent, where $\{e_{i}\}_{i=1}^{n}$ is the standard basis for $n$-dimensional Euclidean space.

Next, we construct parameters for the MREG. Let $d$ be ${n+1}\choose{2}$ and

[TABLE]

Since $\{h_{k}h_{k}^{T}\}_{k=1}^{d}$ forms a basis for real symmetric matrices, for each adjacency matrix $\mathbf{A}_{i}$, there exists a vector $\lambda_{i}$, such that

[TABLE]

Let $F$ be a finite mixture distribution on points $\{\lambda_{i}\}_{i=1}^{2^{{n+1}\choose{2}}}$, that is

[TABLE]

Under this MREG model, for any adjacency matrix $\mathbf{A}_{i}$

[TABLE]

This concludes that the distribution $\mathcal{F}$ and $MREG(F,h_{1},...,h_{d})$ are equal. ∎

Proof:

Suppose that $\hat{h}_{k_{0}}\hat{h}_{k_{0}}^{T}=\sum_{k=1}^{k_{0}-1}\alpha_{k}\hat{h}_{k}\hat{h}_{k}^{T}$. Therefore, from equation (3) it can be noted that the solution for $\hat{\mathbf{\Lambda}}_{*k_{0}}$ is given by

[TABLE]

Note that equation (4) implies that for any $i=1,\ldots,m$,

[TABLE]

and hence, $\hat{\mathbf{\Lambda}}_{ik_{0}}\sum\limits_{k=1}^{k_{0}-1}\alpha_{k}\hat{h}_{k}\hat{h}_{k}^{T}=0$, so $\hat{\mathbf{\Lambda}}_{ik_{0}}=0$ for all $i=1,\ldots,m$. ∎

Proof:

First, we show that $|D_{n}(h,h_{1})-D(h,h_{1})|$ converges uniformly to [math]. To begin with, notice three facts:

• (1)

the set $\{h:\|h\|=1\}$ is compact;

• (2)

for all $h$, the function $\rho(\cdot,h)$ is continuous

• (3)

for all $h$, the function $\rho(\cdot,h)$ is bounded by $n^{2}$.

Therefore, by the uniform law of large numbers [62], we have

[TABLE]

To prove the claim of the theorem, we use a technique similar to that employed by Bickel and Doksum [63]. Let $h^{\prime}$ be a global minimizer of $D(h,h_{1})$. By definition, we must have $D_{m}(\hat{h}_{1}^{m},h_{1})\leq D_{m}(h^{\prime},h_{1})$ and $D(h^{\prime},h_{1})\leq D(\hat{h}_{1}^{m},h_{1})$. From these two inequalities,

[TABLE]

Therefore,

[TABLE]

This implies

[TABLE]

Hence, $|D_{m}(\hat{h}_{1}^{m},h_{1})-D(h^{\prime},h_{1})|$ must converge almost surely to [math], that is

[TABLE]

If $\hat{h}_{1}^{m}$ does not converge almost surely to a global minimizer of $D(h,h_{1})$, then for all $h^{\prime}\in\arg\min D(h,h_{1})$, $\|\hat{h}_{1}^{m}-h^{\prime}\|\geq\epsilon$ for some $\epsilon$ and infinitely many values of $m$. Since $h^{\prime}$ is a global minimum, $|D(\hat{h}_{1}^{m},h_{1})-D(h^{\prime},h_{1})|>\epsilon^{\prime}$ for infinitely many values of $m$ and some $\epsilon^{\prime}$. This contradicts with the previous equation. Therefore, $\hat{h}_{1}^{m}$ must converge almost surely to $\arg\min D(h^{\prime},h_{1})$. ∎

Proof:

The proof of theorem relies on two lemmas. The first lemma shows that $h^{\prime}$ is the eigenvector corresponding to the largest eigenvalue of $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle\mathbf{A}_{i})$. The second lemma shows that $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle\mathbf{A}_{i})$ is close to $E(\lambda_{i}^{2})(h_{1}^{T}h^{\prime})^{2}h_{1}h_{1}^{T}$ under Frobenius norm. Then, we apply Davis-Kahan theorem [64] to establish the result of theorem.

Lemma 6.1

The vector $h^{\prime}$ is the eigenvector corresponding to the largest eigenvalue of $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle\mathbf{A}_{i})$.

We notice that

[TABLE]

Therefore,

[TABLE]

Taking the derivative of $E(\langle\mathbf{A}_{i},hh^{T}\rangle^{2})+c(h^{T}h-1)$ with respect to $h$,

[TABLE]

Setting this expression to [math] yields,

[TABLE]

Using the fact that $\|h^{\prime}\|=1$, we can solve for $c$:

[TABLE]

Then, substituting for $c$,

[TABLE]

Therefore, we see that $h^{\prime}$ is an eigenvector of $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle\mathbf{A}_{i})$ and the corresponding eigenvalue is $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle^{2})$. Furthermore, $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle^{2})$ must be the eigenvalue with the largest magnitude. For if not, then there exists an $h^{\prime\prime}$ with norm $1$ such that

[TABLE]

however, by Cauchy-Schwarz inequality we must have

[TABLE]

implying $E(\langle\mathbf{A}_{i},h^{\prime\prime}h^{\prime\prime T}\rangle^{2})>E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle^{2})$, which contradicts equation (9) of this appendix. Thus, we conclude that $h^{\prime}$ is the eigenvector corresponding to the largest eigenvalue of $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle\mathbf{A}_{i})$.

Lemma 6.2

[TABLE]

We compute $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle\mathbf{A}_{i})$ by conditioning on $P_{i}$.

[TABLE]

Here, $\operatorname{diag}()$ means only keep the diagonal of the matrix; $*$ means the Hadamard product, and $\mathbf{J}$ is a matrix of all ones. Using the fact that $\mathbf{P}_{i}=\lambda_{i}h_{1}h_{1}^{T}$, we have

[TABLE]

If we consider the norm difference between $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle\mathbf{A}_{i})$ and $E(\lambda_{i}^{2})(h_{1}^{T}h^{\prime})^{2}h_{1}h_{1}^{T}$, we have

[TABLE]

This finishes the proof for the lemma.

Notice that the only non-zero eigenvector of $E(\lambda_{i}^{2})(h_{1}^{T}h^{\prime})^{2}h_{1}h_{1}^{T}$ is $h_{1}$ and the corresponding eigenvalue is $E(\lambda_{i}^{2})(h_{1}^{T}h^{\prime})^{2}$. We apply the Davis-Kahan theorem [64] to the eigenvector corresponding to the largest eigenvalue of matrices $E(\langle\mathbf{A}_{i},h^{\prime}h^{\prime T}\rangle\mathbf{A}_{i})$ and $E(\lambda_{i}^{2})(h_{1}^{T}h^{\prime})^{2}h_{1}h_{1}^{T}$, yielding

[TABLE]

∎