# Learning joint and individual structure in network data with covariates

Carson James , Dongbang Yuan , Irina Gaynanova , Jesús Arroyo Department of
Statistics, Texas A&M UniversityMetaDepartment of Biostatistics, University of
MichiganDepartment of Statistics, Texas A&M UniversityCorresponding author
(Email: jarroyo@tamu.edujarroyo@tamu.edu).

###### Abstract

Datasets consisting of a network and covariates associated with its vertices
have become ubiquitous. One problem pertaining to this type of data is to
identify information unique to the network, information unique to the vertex
covariates and information that is shared between the network and the vertex
covariates. Existing techniques for network data and vertex coviariates focus
on capturing structure that is shared but are usually not able to
differentiate structure that is unique to each dataset. This work formulates a
low-rank model that simultaneously captures joint and individual information
in network data with vertex covariates. A two-step estimation procedure is
proposed, composed of an efficient spectral method followed by a refinement
optimization step. Theoretically, we show that the spectral method is able to
consistently recover the joint and individual components under a general
signal-plus-noise model. Simulations and real data examples demonstrate the
ability of the methods to recover accurate and interpretable components. In
particular, the application of the methodology to a food trade network between
countries with economic, developmental and geographical country-level
indicators as covariates yields joint and individual factors that explain the
trading patterns.

##  1 Introduction

Network data is ubiquitous in many disciplines and application domains,
including computer science, statistics, biology, and physics. These data,
encoding relationships between units represented as nodes, are often
accompanied by additional information about the nodes, usually referred to as
node covariates, attributes, or metadata (Newman and Clauset, 2016; Liu, 2019;
Chunaev, 2020). In these situations, a common goal is to understand the
associations between the network connectivity and the node covariates. In our
example, we consider international food commodity trade data represented as a
network, where the nodes correspond to different countries and edge weights
encode food commodity trade volumes between corresponding countries. The
covariates at each node consist of economic and geographic information for
each country, such as gross domestic product (GDP) per capita, birth rate and
region. We wish to exploit that both datasets contain information about the
nodes in order to better understand the structure of the network, node
covariates and their relationship. Specifically, we seek to understand how
economic and geographic factors explain the observed trade between countries,
and identify additional information in the network that cannot be explained
solely by these variables.

There has been substantial work that incorporates network and node covariate
information. Some examples include methods that use node covariates to improve
community detection (Binkiewicz et al., 2017; Huang et al., 2023),
dimensionality reduction (Zhao et al., 2022), regression with network
information (Li et al., 2019) and mixed effect models for network edges (Hoff,
2005). These existing methods aim to exploit the common information in the
network and its node covariates but do not address whether the covariates
contain information that is not captured by the network or vice-versa. Such
methods may also discard information contained in only network or covariate
data.

When integrating multi-view datasets with matched units, a common task is to
identify joint and individual structures across these datasets, for example,
by partitioning the data into components that either capture shared
information or information that is unique to each individual dataset (Lock et
al., 2013; Feng et al., 2018). Existing work in this area has focused on
multi-view data matrices. In constrast, network data provides information
about interactions or relationships between units rather than unit-level data,
making integration of such data more challenging. Recent studies have explored
integration problems in the context of multiple network datasets with matched
units (MacDonald et al., 2022). To the best of our knowledge, the task of
integrating network data with matrix data (representing edge or node
covariates) with the goal of identifying joint and individual signals across
these two distinct sources has not been explored in prior works.

In this work, we study the problem of integrating network data with node
covariates making the following contributions:

  * •

We propose a new unified model that partitions the information among node
connectivity and node covariate data into joint and individual components.
These node components can be interpreted as node embeddings that capture
either joint or individual structure, and we discuss the connections between
these embeddings and latent space models (Hoff et al., 2002).

  * •

We develop an estimation procedure based on a computationally efficient
spectral method followed by an optimization-based refining step. We derive a
fast optimization algorithm for the latter based on closed-form block-
coordinate descent updates.

  * •

Under a general signal-plus-noise error model, we derive an explicit
characterization of the error rate for the joint and individual components
obtained from the spectral method and prove consistency of the spectral
estimates under reasonable restrictions on the model parameters.

  * •

We showcase that the proposed method effectively captures joint and individual
information via simulations and the analysis of an international food
commodity trade dataset paired with econometric data. We find that the joint
components for the network and node covariates capture joint trading structure
related to the growth-domestic product, whereas the individual network
components capture local and global trading patterns that cannot be explained
by the covariates alone.

To further motivate our methodology, consider a synthetic network data example
in Figure 1(a) (the full data generation process is discussed in Section 4.1)
with $n=40$ units (nodes) divided into four groups (indicated by color and
shape in Figure 1(a)). The nodes are divided into three communities (indicated
by shape), with groups 3 and 4 merged into one community. The node covariate
data consists of a mixture model with three clusters (indicated by color),
where groups 1 and 2 are merged into the same cluster. In this example,
neither the connectivity data nor the covariate data alone can differentiate
all four groups. Our methodology allows to discover joint components
differentiating the nodes into two larger clusters that are common across two
datasets (Figure 1(c)), and individual components that identify the remaining
groups that cannot be distinguished by a single view (Figure 1(d)).

(a) Network connectivity plot

(b) Top node covariate principal components

(c) Joint components

(d) Individual components

Figure 1: Synthetic dataset composed of a network with $40$ nodes (top left)
and three node-level covariates represented via their leading principal
components (top right). The nodes are partitioned into four groups, but each
dataset can only detect three groups. Shape indicates nodes distinguished by
the network data. Color indicates nodes distinguished by covariate data. The
joint (bottom left) and individual (bottom right) components estimated by our
methodology partition the information into the joint and individual group
structures across the datasets. See Section 4.1 for details.

###  1.1 Related Work

In the analysis of network data, it is common to represent the nodes via a
low-dimensional latent space (Hoff et al., 2002). This can be achieved by
estimating node embeddings through low-rank decompositions of the adjacency
matrix or the Laplacian (Rohe et al., 2011; Athreya et al., 2017). The
resulting low-dimensional representation is often used for downstream
inference tasks such as community detection. We follow a similar approach in
this work by obtaining low-dimensional representations of the nodes that
capture joint and individual variation across network data and covariates.

In the context of matrix datasets with matched units, there has been
substantial work on data integration techniques that appeal to the notion of
joint structure. Classical approaches include principal component analysis
(PCA) on combined multiple covariate datasets (Wold et al., 1996; Shen et al.,
2009) and canonical correlation analysis (Harold, 1936; Witten and Tibshirani,
2009). More recent approaches decompose the data into joint and individual
components (Lock et al., 2013; Feng et al., 2018; Gaynanova and Li, 2019;
Murden et al., 2022; Chen et al., 2022; Palzer et al., 2022). However, to our
knowledge, there persists a gap in the literature for similar integration
techniques for network data problems.

In the context of the analysis of network data with node covariates, much of
the existing literature is concerned with two types of approaches. On one
hand, the information in one of the datasets is modeled as a function of the
other. For example, the network may be treated as auxiliary information in
prediction problems with linked units (Li et al., 2019; Le and Li, 2022), or
alternatively, the edges of the network can be modeled as a function of node
and edge covariate information (Hoff, 2005). The second type of approach,
which is the direction we follow in this paper, studies both datasets jointly
to identify joint structures between the network and its nodal attributes.
Common tasks in this approach include community detection with node features
(Zhang et al., 2016; Binkiewicz et al., 2017; Deshpande et al., 2018; Yan and
Sarkar, 2021), or testing for dependence between a network and its nodal
attributes(Fosdick and Hoff, 2015; Lee et al., 2019). These methods are only
concerned with identifying common information across the datasets, but are
unable to distinguish between the structure that is shared or unique to each
dataset, which is our main goal in our analysis. Preparing this work, we
became aware of recent concurrent work that also introduced a notion of joint
and individual in network data with covariates (Li et al., 2024), but the
model and methodology are substantially different from ours.

Another related line of work is the study of multiple network datasets with
matched nodes, often referred to as multilayer or multiplex networks. A common
approach in this setup is to identify joint structure across the networks via
node embeddings (Eynard et al., 2015; Wang et al., 2019; Zhang et al., 2018;
Paul and Chen, 2020; Wang et al., 2019; Arroyo et al., 2021), but these
methods are typically unable to identify individual structures. Recent
progress in this direction has been made by MacDonald et al. (2022) who
partition the model structure into joint and individual parts, but extensions
to network data with covariates are not immediate. Although our model shares
some similarities with MacDonald et al. (2022), the approach we follow
considers modeling the column spaces of the matrices instead of directly
modeling the latent factors as in their work, resulting in more flexibility to
handle differences in scaling or rotation factors across the datasets.

From an algorithmic perspective, our methodology is based on a combination of
an efficient spectral algorithm to jointly embed the network and the
covariates, followed by an iterative optimization procedure as a refinement of
the initial estimates. Our spectral method shares some similarities with the
AJIVE method of Feng et al. (2018) for matrix data, but our model is
fundamentally different as it treats network data with a connection to a
latent space model. Furthermore, we provide theoretical guarantees on the
error rate of the spectral method, which are still not available in the matrix
setup. The spectral method is also related to algorithms for finding common
low-rank embeddings (Crainiceanu et al., 2011; Arroyo et al., 2021), but the
latter methods are unable to recover individual components.

###  1.2 Notation.

The column space of a matrix $M\in\mathbb{R}^{n\times m}$ is denoted by
$\mathcal{C}(M)$, and the orthogonal complement of the subspace
$\mathcal{C}(M)$ is denoted by $\mathcal{C}(M)^{\perp}$. The set of $n\times
p$ matrices with orthonormal columns is denoted by
$\mathbb{O}_{n,p}=\\{Q\in\mathbb{R}^{n\times p}:Q^{\top}Q=I\\}$, and we write
$\mathbb{O}_{n}$ in place of $\mathbb{O}_{n,n}$. Given matrices
$U,V\in\mathbb{R}^{n\times m}$, we denote by $U\cong_{O}V$ if there exists
$W\in\mathbb{O}_{p}$ such that $U=VW$. Given a symmetric matrix
$A\in\mathbb{R}^{n\times n}$ with eigendecomposition $U\Lambda U^{\top}=A$,
where $U=(u_{1},\ldots,u_{n})\in\mathbb{O}_{n}$ contains the eigenvectors and
$\Lambda=\text{diag}(\lambda_{1},\ldots,\lambda_{n})$ contains the eigenvalues
in its diagonal with $|\lambda_{1}|\geq\cdots\geq|\lambda_{n}|$, we denote the
matrix of $k$ leading eigenvectors of $A$ ordered by eigenvalue magnitude by
$\operatorname{EIG}(A,k)=(u_{1},\ldots,u_{k})\in\mathbb{R}^{n\times k}$. Given
a matrix $X\in\mathbb{R}^{n\times p}$, let $X=U\Sigma V^{\top}$ be its
singular value decomposition, where
$U=(u_{1},\ldots,u_{m})\in\mathbb{O}_{n\times m}$, $V\in\mathbb{O}_{p\times
m}$, and $\Sigma=\text{diag}(\sigma_{1},\ldots,\sigma_{m})$ with
$\sigma_{1}\geq\cdots\geq\sigma_{m}>0$. The $j$-th largest singular value of
$X$ is denoted by $\sigma_{j}(X)$, and the $k$ left leading singular vectors
of $X$ ordered by singular value magnitude are denoted by
$\operatorname{SV}(X,k)=(u_{1},\ldots,u_{k})\in\mathbb{R}^{n\times k}$. Let
$f,g:\mathbb{N}\rightarrow[0,\infty)$. We write $f=O(g)$ and $g=\Omega(f)$ if
there exists $M>0$ and $n_{0}\in\mathbb{N}$ such that for each
$n\in\mathbb{N}$, $n\geq n_{0}$ implies that $f(n)\leq Mg(n)$. We write
$f\asymp g$ if $f=O(g)$ and $f=\Omega(g)$. We write $f=o(g)$ and $g=\omega(f)$
if for each $\epsilon>0$, there exist $n_{0}\in\mathbb{N}$ such that for each
$n\in\mathbb{N}$, $n\geq n_{0}$ implies that $f(n)<\epsilon g(n)$. We denote
the Frobenius and spectral norms by $\|\cdot\|_{F}$ and $\|\cdot\|$
respectively.

##  2 Methodology

In this section, we propose a model for network data with covariates, and
develop estimation procedures for the model parameters.

###  2.1 Modeling joint and individual structure

We consider the problem of jointly modeling the structure of a network (a
graph with $n$ nodes) and node-level covariate information. We represent the
network by an observed adjacency matrix $A\in\mathbb{R}^{n\times n}$, where
the entry $A_{ij}$ encodes the value of the interaction between nodes $i$ and
$j$. We represent node-level covariates by a matrix $X\in\mathbb{R}^{n\times
p}$ so that the $i$-th row of $X$ represents $p$ variables measured at node
$i$. For simplicity, we focus on the undirected network setup, in which case
$A$ is symmetric.

We adopt a signal-plus-noise matrix model for the adjacency matrix
$A\in\mathbb{R}^{n\times n}$ and the covariate matrix $X\in\mathbb{R}^{n\times
p}$ so that

$\displaystyle A=P+E^{A},\quad\quad X=W+E^{X},$ (1)

where $P=\mathbb{E}[A]\in\mathbb{R}^{n\times n}$ and
$W=\mathbb{E}[X]\in\mathbb{R}^{n\times p}$ are the expected values, $E^{A}$
and $E^{X}$ are random matrices such that the collection of upper triangular
entries of $E^{A}$ and the entries of $E^{P}$ are mutually independent random
variables, and $P$ and $E^{A}$ are symmetric matrices. The matrix $P$
represents the expected node connection strengths (network signal), whereas
$W$ represents the signal of the node covariates. The additive decomposition
for $A$ incorporates many existing models. For instance, when elements of $A$
are binary, it is common to model its entries as independent Bernoulli random
variables with $A_{ij}{\sim}\text{Bern}(P_{ij})$, where $P_{ij}\in[0,1]$
indicates the probability of connection between nodes $i$ and $j$.

To model joint and individual structures across the network and covariate
signals, we propose to partition the column spaces of signal matrices $P$ and
$W$ such that the joint structure is the intersection of these spaces, and the
individual structures are the residuals after removing the intersection (Lock
et al., 2013; Feng et al., 2018). We formalize this with the definition below.

######  Definition 2.1.

Let $P\in\mathbb{R}^{n\times n}$ and $W\in\mathbb{R}^{n\times p}$ be the
expected values of the network adjacency matrix and node covariate matrix from
Equation (1). We define the _joint subspace_ , _individual subspace for the
network_ and the _individual subspace for the covariates_ as

Joint: $\displaystyle\mathcal{M}$
$\displaystyle=\mathcal{C}(P)\cap\mathcal{C}(W),$ (2) Network Individual:
$\displaystyle\mathcal{R}^{(1)}$
$\displaystyle=\mathcal{C}(\mathcal{P}_{\mathcal{M}^{\perp}}P),$ (3) Covariate
Individual: $\displaystyle\mathcal{R}^{(2)}$
$\displaystyle=\mathcal{C}(\mathcal{P}_{\mathcal{M}^{\perp}}W).$ (4)

We denote the respective dimensions of these subspaces by $r_{M}$, $r_{1}$ and
$r_{2}$. If $M\in\mathbb{O}_{n,r_{M}}$, $R^{(1)}\in\mathbb{O}_{n,r_{1}}$ and
$R^{(2)}\in\mathbb{O}_{n,r_{2}}$ satisfy $\mathcal{C}(M)=\mathcal{M}$,
$\mathcal{C}(R^{(1)})=\mathcal{R}^{(1)}$ and
$\mathcal{C}(R^{(2)})=\mathcal{R}^{(2)}$, we refer to $M$, $R^{(1)}$ and
$R^{(2)}$ as the _joint components_ , _network individual components_ and
_covariate individual components_ respectively.

The joint and individual subspaces in Definition 2.1 are uniquely defined by
construction and satisfy $\operatorname{rank}(P)=r_{M}+r_{1}$ and
$\operatorname{rank}(W)=r_{M}+r_{2}$. While the subspaces are unique, the
joint and individual components $M,R^{(1)},R^{(2)}$ are only unique up to
multiplication by orthogonal matrix. We formalize the uniqueness conditions as
well as decomposition of $P$ and $W$ through these components below under an
additional nontriviality assumption.

######  Assumption 1.

The joint and individual subspaces $\mathcal{M},\mathcal{R}^{(1)}$ and
$\mathcal{R}^{(2)}$ from Definition 2.1 are non-trivial, that is,
$\mathcal{M},\mathcal{R}^{(1)},\mathcal{R}^{(2)}\neq\\{0\\}$.

Assumption 1 implies that $r_{M}$, $r_{1}$, $r_{2}>0$, meaning that each
subspace in Definition 2.1 is non-empty.

######  Lemma 2.1.

Let $P\in\mathbb{R}^{n\times n}$ and $W\in\mathbb{R}^{n\times p}$ be as in
Equation (1) and let $r_{M},r_{1},r_{2}$ be the ranks from Definition 2.1.
Under Assumption 1, there exist $M\in\mathbb{O}_{n,r_{M}}$,
$R^{(k)}\in\mathbb{O}_{n,r_{k}}$ for $k=1,2$, as in Definition 2.1, and full
rank matrices $\Gamma^{(1)}\in\mathbb{R}^{(r_{M}+r_{1})\times n}$ and
$\Gamma^{(2)}\in\mathbb{R}^{(r_{M}+r_{2})\times p}$ such that

  1. (a)

$P=\begin{pmatrix}M&R^{(1)}\end{pmatrix}\Gamma^{(1)}$ and
$W=\begin{pmatrix}M&R^{(2)}\end{pmatrix}\Gamma^{(2)}$,

  2. (b)

$M\perp R^{(k)}$, $k=1,2$,

  3. (c)

$M,R^{(1)}$ and $R^{(2)}$ are unique up to right multiplication by an
orthogonal matrix.

Lemma 2.1 establishes the identifiability of the components up to an
orthogonal matrix multiplication, and shows that the joint and individual
components are orthogonal to each other. Related results for matrix-valued
data have been studied in Feng et al. (2018); Yuan and Gaynanova (2022). The
decomposition in (a) characterizes the column spaces of the matrices in terms
of such components. As a consequence, the network components
$\begin{pmatrix}M&R^{(1)}\end{pmatrix}$ are equivalent to the non-trivial
eigenvectors of $P$ up to an orthogonal transformation, whereas the covariate
components $\begin{pmatrix}M&R^{(2)}\end{pmatrix}$ are the left singular
vectors of $W$ up to an orthogonal transformation, which motivates us to
utilize a spectral method for estimation (see Section 2.2.1).

Definition 2.1 relates to several existing methods. When the network adjacency
matrix is replaced by a non-symmetric matrix of covariates (no network data),
Definition 2.1 coincides with existing models for joint and individual signals
for matrix-valued data (Lock et al., 2013; Feng et al., 2018). In comparison
to existing approaches for network data, the special case of $r_{1}=r_{2}=0$
(column spaces of $P$ and $W$ are the same) is similar to models that consider
common latent factors or community structure (Binkiewicz et al., 2017; Arroyo
et al., 2021); thus, our model is more general. A recent work by MacDonald et
al. (2022) considers a similar decomposition in terms of both joint and
individual components for multiplex networks. However, their model imposes
this decomposition in the factorization of the latent positions directly. In
contrast, we use the column spaces, which allow for more flexibility via
different scaling and rotation factors across the components of $P$ and $W$.

For estimation to be feasible, we assume that the signal matrices $P$ and $W$
are low rank, that is $n\gg r_{M}+r_{1}$, $\min(n,p)\gg r_{M}+r_{2}$, which is
a common assumption for network and matrix data (Athreya et al., 2017; Han et
al., 2023; Udell and Townsend, 2019). However, for network data, the low-rank
assumption may not hold exactly but only approximately. For a simple example,
if the network has no self-loops, the diagonal entries of $P$ are zero, which
can violate the low-rank assumption. In these cases, Definition 2.1 can be
modified to work with the leading eigenspace of $P$ and the leading singular
subspace of $W$. Formally, if $V^{(1)}\in\mathbb{O}_{n,r_{P}}$ is the matrix
of leading eigenvectors of $P$ ($V^{(1)}=\operatorname{EIG}(P,r_{P})$) and
$V^{(2)}\in\mathbb{O}_{n,r_{W}}$ is the matrix of left-leading singular
vectors of $W$ ($V^{(2)}=\operatorname{SV}(W,r_{W})$) for any given values of
$r_{P}$ and $r_{W}$, the joint and individual subspaces can be defined in an
analogous manner as
$\widetilde{\mathcal{M}}=\mathcal{C}(V^{(1)})\cap\mathcal{C}(V^{(2)})$,
$\widetilde{\mathcal{R}}^{(1)}=\mathcal{C}(\mathcal{P}_{{\widetilde{\mathcal{M}}}^{\perp}}V^{(1)})$
and
$\widetilde{\mathcal{R}}^{(2)}=\mathcal{C}(\mathcal{P}_{{\widetilde{\mathcal{M}}}^{\perp}}V^{(2)})$,
respectively. This definition coincides with Definition 2.1 when
$r_{P}=\operatorname{rank}(P)$ and $r_{W}=\operatorname{rank}(W)$. We opted to
work with a low-rank assumption for simplicity of notation, but all our
results continue to hold if the matrices are approximately low-rank.

####  2.1.1 Implications for latent space models with covariates

While our model in Lemma 2.1 is stated in terms of the column space of the
network, it has a connection to latent space models (Hoff et al., 2002) and
implies a decomposition of the latent positions of the network into joint and
individual parts. Concretely, consider the _generalized random dot product
graph_ (GRDPG) model (Rubin-Delanchy et al., 2022). This model represents the
nodes of the network with latent positions in some low-dimensional space and
characterizes the probability of two nodes being connected by an edge as an
indefinite inner product between their corresponding latent positions. This
model is particularly suitable to characterize networks with low-rank expected
adjacency matrices, and encompasses multiple models of interest in the
literature, including the stochastic blockmodel (Holland et al., 1983) and its
extensions. For a survey, see Athreya et al. (2017).

Specifically, for a binary and symmetric adjacency matrix $A$ whose edges
$(A_{ij})_{i<j}$ are independent, the GRDPG model assigns to each node
$i=1,\ldots,n,$ a corresponding latent position denoted by
$y_{i}\in\mathbb{R}^{d}$, and the probability of a connection between nodes
$i$ and $j$ is given by the indefinite inner product of their associated
latent positions $P_{ij}=y_{i}^{\top}I_{a,b}y_{j}$, where
$I_{a,b}=\text{diag}(1,\ldots,1,-1,\ldots,-1)$ is a $d\times d$ diagonal
matrix with $a$ ones and $b$ minus ones on its diagonal such that $a+b=d$. We
collect these latent positions into a matrix
$Y=\begin{pmatrix}y_{1},\cdots,y_{n}\end{pmatrix}^{\top}\in\mathbb{R}^{n\times
d}$ so that the expected adjacency matrix can be expressed as
$P=YI_{a,b}Y^{\top}$. We denote $A\sim\operatorname{GRDPG}(Y,a,b)$ for a
network adjacency matrix generated under this model. If $b=0$, this model
reduces to the _random dot product graph_ (RDPG) (Young and Scheinerman, 2007)
and we write $A\sim\operatorname{RDPG}(Y)$.

When a graph is modeled jointly with node covariates as in Definition 2.1, the
latent positions of the GRDPG model can be decomposed into joint and
individual components. In particular, the column space of $P$ is the same as
that of the latent position matrix, that is, $\mathcal{C}(P)=\mathcal{C}(Y)$.
Using the factorization for $P$ given in Lemma 2.1, we can decompose the
latent positions $Y$ into joint and individual parts as demonstrated next.

######  Corollary 2.2.

Let $P$ and $W$ be matrices as in Equation (1) such that $P$ is the
probability matrix of a generalized random dot product graph model
$\operatorname{GRDPG}(Y,a,b)$ with latent position matrix
$Y\in\mathbb{R}^{n\times d}$. Suppose that $M,R^{(1)}$ and $R^{(2)}$ are the
joint and individual components from Definition 2.1. Then, there exists a
matrix $\widetilde{\Gamma}^{(1)}\in\mathbb{R}^{(r_{M}+r_{1})\times d}$ such
that

$Y=\begin{pmatrix}M&R^{(1)}\end{pmatrix}\widetilde{\Gamma}^{(1)}=\mathcal{P}_{\mathcal{C}(M)}Y+\mathcal{P}_{\mathcal{C}(R^{(1)})}Y.$
(5)

This corollary is a direct consequence of Lemma 2.1. By virtue of the
factorization in Equation (5), we can interpret the rows of $M$ as the node
embeddings that explain the joint structure with the covariate information,
whereas the rows of $R^{(1)}$ correspond to latent positions individual to the
network. For example, in the stochastic blockmodel (Holland et al., 1983)
(SBM), which is a special case of the GRDPG model, the rows of $Y$ take only
$K$ different values, corresponding to the clusters or communities of the
model, and the decomposition in Equation (5) breaks these $K$ groups according
to the structure in the covariates. In particular, if the covariate data also
presents clustering structure (as illustrated in Figure 1(a)), the joint and
individual components can detect clustering structure that is common or
individual across the two datasets. We elaborate an example of this situation
in Section 4.1.

###  2.2 Estimation

In this section, we describe our methodology for estimating joint and
individual components in network data with covariates as defined in Definition
2.1. We first obtain estimates via a simple and efficient spectral method. We
then refine these initial estimates via joint optimization to obtain the final
estimates. Throughout, we assume that the ranks $r_{M}$, $r_{1}$ and $r_{2}$
are known. In our application, we estimate the ranks by finding elbows in
scree plots (Section 5), but alternative methods for rank estimation based on
singular value thresholding can be used (Chatterjee, 2015; Gavish and Donoho,
2017; Prothero et al., 2024).

####  2.2.1 Spectral Method

To motivate the spectral method, we start by focusing on the noiseless case,
in which the expectation matrices $P$ and $W$ are directly observed. We set
$V^{(1)}=\operatorname{SV}(P,r_{M}+r_{1})$ and
$V^{(2)}=\operatorname{SV}(W,r_{M}+r_{2})$ as the leading eigenvector and left
singular vector matrices of $P$ and $W$, respectively. Since $V^{(k)}$ and
$\begin{pmatrix}M&R^{(k)}\end{pmatrix}$ have the same column space,
$V^{(k)}\cong_{O}\begin{pmatrix}M&R^{(k)}\end{pmatrix}$. The key idea behind
the spectral estimation procedure is to combine the datasets in such a way
that the joint subspace has a stronger signal in the combination than in
either of the original datasets. To this end, we consider the matrix
$U=\begin{pmatrix}V^{(1)}&V^{(2)}\end{pmatrix}$ and observe that its $r_{M}$
leading left singular values correspond to $M$ (modulo some orthogonal
matrix). Once we have $M$, we can obtain $R^{(k)}$ using the fact that

$\operatorname{SV}(\mathcal{P}_{\mathcal{C}(M)^{\perp}}V^{(k)},r_{k})\cong_{O}\operatorname{SV}(\begin{pmatrix}0&R^{(k)}\end{pmatrix},r_{k})\cong_{O}R^{(k)}.$

In practice, as the matrices $P$ and $W$ are not observed, we apply this
process to observed $A$ and $X$ as detailed in Algorithm 1. When applied to
the noiseless $P$ and $W$, the algorithm indeed recovers the true model
parameters $M,R^{(1)},R^{(2)}$, as formalized below.

######  Proposition 2.3.

Let $P\in\mathbb{R}^{n\times n}$ and $W\in\mathbb{R}^{n\times p}$ be the
signal matrices as in Equation (1), and let $r_{M},r_{1},r_{2}$ be the ranks
as in Definition 2.1. Suppose that Algorithm 1 is started with
$(P,W,r_{M},r_{1},r_{2})$ as the input. Then, the algorithm returns the
matrices $M,R^{(1)},R^{(2)}$ that satisfy Lemma 2.1, up to right
multiplication by an orthogonal matrix.

Algorithm 1 Spectral method for joint and individual component estimation

$A\in\mathbb{R}^{n\times n}$, $X\in\mathbb{R}^{n\times p}$, $r_{M}$, $r_{1}$,
$r_{2}$

$\widehat{V}^{(1)}=\operatorname{EIG}(A,r_{M}+r_{1})$

$\widehat{V}^{(2)}=\operatorname{SV}(X,r_{M}+r_{2})$

$\widehat{U}=\begin{pmatrix}\widehat{V}^{(1)}&\widehat{V}^{(2)}\end{pmatrix}$

$\widehat{M}=\operatorname{SV}(\widehat{U},r_{M})$

$\widehat{R}^{(1)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(\widehat{M})^{\perp}}\widehat{V}^{(1)},r_{1})$

$\widehat{R}^{(2)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(\widehat{M})^{\perp}}\widehat{V}^{(2)},r_{2})$

$\widehat{M},\widehat{R}^{(1)},\widehat{R}^{(2)}.$

####  2.2.2 Optimization Algorithm Refinement

While the spectral method in Section 2.2.1 is simple to implement and
computationally fast, it can discard nontrivial signals in the initial low-
rank approximation step. To remedy this, we propose a refinement estimation
procedure based on optimizing an objective function associated with the model.
To measure the loss for the covariates, we consider the standard squared
Frobenius norm error, $\|X-W\|_{F}^{2}$. For the network, due to the symmetry
of the adjacency matrix, we consider non-symmetric matrices
$A^{\prime}=(|\widehat{\lambda}_{1}|^{1/2}\widehat{u}_{1},\ldots,|\widehat{\lambda}_{n}|^{1/2}\widehat{u}_{n})$
and $P^{\prime}=(|\lambda_{1}|^{1/2}u_{1},\ldots,|\lambda_{n}|^{1/2}u_{n})$
where $\widehat{\lambda}_{1},\ldots,\widehat{\lambda}_{n}$ and
$\widehat{u}_{1},\ldots,\widehat{u}_{j}$ are the eigenvalues and eigenvectors
of $A$ ordered by magnitude, and $\lambda_{1},\ldots,\lambda_{n}$ and
$u_{1},\ldots,u_{j}$ are the eigenvalues and eigenvectors of $P$ ordered by
magnitude. Since $A$ and $A^{\prime}$ have the same column space (similarly,
$P$ and $P^{\prime}$) the conditions as in Lemma 2.1 still apply. When $A$ is
positive-definite, $A^{\prime}$ and $P^{\prime}$ are the square roots of $A$
and $P$, respectively. We consider the Frobenius norm error
$\|A^{\prime}-P^{\prime}\|_{F}^{2}$ to measure the loss for the network.

To further refine the initial estimates obtained from Algorithm 1, we propose
to minimize the sum of the losses associated with the network and covariate
signals subject to identifiability conditions in Lemma 2.1. Specifically,
since $\mathcal{C}(P^{\prime})=\mathcal{C}(P)=\mathcal{C}(M\ R^{(1)})$ and
$\mathcal{C}(W)=\mathcal{C}(M\ R^{(2)})$, under the identifiability conditions
the signals minimizing Frobenius norm error loss must satisfy

$\widehat{P^{\prime}}=\mathcal{P}_{\mathcal{C}(\widehat{M})}A^{\prime}+\mathcal{P}_{\mathcal{C}(\widehat{R}^{(1)})}A^{\prime},\quad\quad\widehat{W}=\mathcal{P}_{\mathcal{C}(\widehat{M})}X+\mathcal{P}_{\mathcal{C}(\widehat{R}^{(2)})}X.$

This leads us to the following constrained optimization problem

$\displaystyle\operatorname*{minimize}_{M,R^{(1)},R^{(2)}}$
$\displaystyle\left\\{\|A^{\prime}-\mathcal{P}_{\mathcal{C}(M)}A^{\prime}-\mathcal{P}_{\mathcal{C}(R^{(1)})}A^{\prime}\|_{F}^{2}+\|X-\mathcal{P}_{\mathcal{C}(M)}X-\mathcal{P}_{\mathcal{C}(R^{(2)})}X\|_{F}^{2}\right\\}$
(6) s.t.
$\displaystyle\quad\begin{pmatrix}M&R^{(k)}\end{pmatrix}\in\mathbb{O}_{n,r_{M}+r_{k}}.$
(7)

The optimization problem above seeks the best low-rank approximations to the
observed data matrices simultaneously under the joint subspace constraint
imposed by the model. To ensure equal weights of both terms in the loss
function, we scale $A^{\prime}$ and $X$ by dividing all elements by the
Frobenius norm of $\operatorname{SV}(A^{\prime},r_{M}+r_{1})$ and
$\operatorname{SV}(X,r_{M}+r_{2})$, respectively. In the remainder of this
section, we assume that $A^{\prime}$ and $X$ have already been scaled.

To solve the optimization problem (6), we use the initial estimates obtained
from Algorithm 1, and perform iterative updates. That is, we fix $M$ and
optimize for $R^{(k)}$, $k=1,2$, and subsequently fix $R^{(k)}$ and optimize
for $M$. These updates are alternated until convergence. We describe each
update below.

When $M$ is fixed, by Lemma 2 in Yuan and Gaynanova (2022), the updates for
$R^{(1)}$ and $R^{(2)}$ are closed form and are given by

$\widehat{R}^{(1)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(M)^{\perp}}A^{\prime},r_{1}),\quad\quad\widehat{R}^{(2)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(M)^{\perp}}X,r_{2}).$

When $R^{(1)}$ and $R^{(2)}$ are fixed, we derive the closed-form update for
$M$ based on the following lemma.

######  Lemma 2.4.

Consider optimization problem (6) with respect to $M$ with fixed
$R^{(k)}\in\mathbb{O}_{n,r_{k}}$ for $k=1,2$ and fixed joint rank $r_{M}$. Let
$Y=\begin{pmatrix}\mathcal{P}_{\mathcal{C}(R^{(1)})^{\perp}}A^{\prime}&\mathcal{P}_{\mathcal{C}(R^{(2)})^{\perp}}X\end{pmatrix}$
and $R=\begin{pmatrix}R^{(1)}&R^{(2)}\end{pmatrix}$. Set
$\widehat{M}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(R)^{\perp}}Y,r_{M})$.
If $\operatorname{rank}(\mathcal{P}_{\mathcal{C}(R)^{\perp}}Y)\geq r_{M}$,
then $\widehat{M}$ is a global minimizer of (6).

In Lemma 2.4, we assume $\mathcal{P}_{\mathcal{C}(R)^{\perp}}Y$ is of rank at
least $r_{M}$. If $R=\begin{pmatrix}R^{(1)}&R^{(2)}\end{pmatrix}$ contains a
complete basis for $\mathbb{R}^{n}$, then
$\mathcal{P}_{\mathcal{C}(R)^{\perp}}Y=0$, which means there is no solution
for $M$ with given rank $r_{M}$. In practice, we never encounter such
situations as the signal ranks are smaller than the maximal ranks, and our
starting values for $M$, $R^{(1)}$ and $R^{(2)}$ always satisfy the
assumption.

Algorithm 2 summarizes overall iterative updates for problem (6). While each
block update has a closed-form solution, the final output can be affected by
the initial starting points as the overall problem is nonconvex. We
purposefully use our initial estimates obtained from Algorithm 1 as a starting
point, and thus Algorithm (6) can be viewed as an optimization-based
refinement of the spectral estimates. In Section 3, we further prove that the
initial estimates are consistent.

Algorithm 2 Iterative descent method for joint and individual component
estimation

Scaled data matrices, $A^{\prime}\in\mathbb{R}^{n\times n}$,
$X\in\mathbb{R}^{n\times p}$; initial estimate $M_{(0)}\in\mathbb{R}^{n\times
r_{M}}$; $r_{1},r_{2},t_{max},\epsilon>0$.

$t=0$; $L_{0}=0$

$R^{(1)}_{(0)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(M_{(0)})^{\perp}}A^{\prime},r_{1})$

$R^{(2)}_{(0)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(M_{(0)})^{\perp}}X,r_{2})$

repeat

Update $M$:

$Y_{(t)}=\begin{pmatrix}\mathcal{P}_{\mathcal{C}(R^{(1)}_{(t)})^{\perp}}A^{\prime}&\mathcal{P}_{\mathcal{C}(R^{(2)}_{(t)})^{\perp}}X\end{pmatrix}$

$R_{(t)}=\begin{pmatrix}R^{(1)}_{(t)}&R^{(2)}_{(t)}\end{pmatrix}$

$M_{(t+1)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(R_{(t)})^{\perp}}Y_{(t)},r_{M})$

Update $\widetilde{R}^{(1)}$ and $\widetilde{R}^{(2)}$:

$R^{(1)}_{(t+1)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(M_{(t+1)})^{\perp}}A^{\prime},r_{1})$

$R^{(2)}_{(t+1)}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(M_{(t+1)})^{\perp}}X,r_{2})$

Update iteration and loss:

$t=t+1$

$V^{(k)}_{(t)}=\begin{pmatrix}M_{(t)}&R^{(k)}_{(t)}\end{pmatrix}$ $k=1,2$

$L_{(t)}=\|A^{\prime}-\mathcal{P}_{V^{(1)}_{(t)}}A^{\prime}\|_{F}^{2}+\|X-\mathcal{P}_{V^{(2)}_{(t)}}X\|_{F}^{2}$

until $t>t_{max}$ or $|L_{(t)}-L_{(t-1)}|\leq\epsilon$

$\widehat{M}=M_{(t)}$

$\widehat{R}^{(k)}=R^{(k)}_{(t)}$ $k=1,2$

$\widehat{M}$, $\widehat{R}^{(1)}$, $\widehat{R}^{(2)}$

##  3 Estimation Error of Spectral Method

In this section, we characterize the estimation error of the spectral method
from Section 2.2.1. Throughout the analysis, we assume the true ranks
$r_{M},r_{1},r_{2}$ are known. To quantify the estimation error, we use the
Procrustes distance which is invariant to orthogonal transformations (Lemma
2.1). For matrices with orthonormal columns $U,V\in\mathbb{O}_{n,r}$, it is
defined as

$d(U,V)=\inf_{Q\in\mathbb{O}_{r}}\|U-VQ\|_{F}.$ (8)

Our goal is to characterize the expected distance between the true and
estimated components, that is $\mathbb{E}[d(\widehat{M},M)]$ and
$\mathbb{E}[d(\widehat{R^{(k)}},R^{(k)})]$, $k=1,2$. We will obtain such
characterization by leveraging perturbation theory for eigenspaces and
singular subspaces (Yu et al., 2015; Cai and Zhang, 2018) under the following
distributional assumptions.

######  Assumption 2.

Let $A\in\mathbb{R}^{n\times n}$ and $X\in\mathbb{R}^{n\times p}$ be the
observed adjacency and covariate matrices with expectation matrices $P$ and
$W$ as in Equation (1) The following distributional assumptions hold.

  1. (a)

$\|A-P\|$ is a subexponential random variable with parameter $\kappa>0$, that
is, for each $q\geq 1$, $\mathbb{E}(\|A-P\|^{q})^{1/q}\leq\kappa q$.

  2. (b)

The entries of $X$ are independent subgaussian random variables with
$\text{Var}(X_{ij})=\tau^{2}$, that is, there exists $c>0$ such that for each
$t\geq 0$, $\mathbb{P}(\|X_{ij}\|\geq t)\leq 2e^{-t^{2}/c}$.

Many existing network models satisfy this assumption, including the GRDPG with
binary entries discussed in Section 2.1.1.

One of the main difficulties in estimating the components from Definition 2.1
is distinguishing the joint components from individuals, as the latter are not
necessarily orthogonal. When the individual subspaces are closely aligned,
they can be mistaken for joint subspaces in the presence of noise, making
estimation difficult. To provide a theoretical characterization of this
difficulty, we use the first principal angle between the individual subspaces
$\mathcal{R}^{(1)}$ and $\mathcal{R}^{(2)}$ to quantify their separation as

$\delta=\delta(\mathcal{R}^{(1)},\mathcal{R}^{(2)})=1-\cos\\{\theta_{1}(\mathcal{R}^{(1)},\mathcal{R}^{(2)})\\}=1-\sigma_{1}({R^{(1)}}^{\top}R^{(2)}),$
(9)

where $R^{(1)}$, $R^{(2)}$ are as in Definition 2.1. Since
$\mathcal{R}^{(1)}\cap\mathcal{R}^{(2)}=\\{0\\}$ by definition, and
$\mathcal{R}^{(1)},\mathcal{R}^{(2)}\neq\\{0\\}$ by Assumption 1, we have a
strictly positive separation, that is $\delta>0$.

We now state our main result.

######  Theorem 3.1.

Let $A\in\mathbb{R}^{n\times n}$ and $X\in\mathbb{R}^{n\times p}$ be the
observed adjacency and covariate matrices as in Equation (1) and let
$r_{M},r_{1},r_{2}\in\mathbb{N}$ be the joint and individual subspace
dimensions as in Definition 2.1. Define
$(\widehat{M},\widehat{R}^{(1)},\widehat{R}^{(2)})$ to be the output of
$Algorithm\ref{alg:alg:model:spectral:1}$ when supplied
$(A,X,r_{M},r_{1},r_{2})$ as arguments. Suppose that Assumption 1 and
Assumption 2 hold, and set

$\epsilon^{(1)}=\frac{\kappa\sqrt{r_{M}+r_{1}}}{\lambda_{r_{M}+r_{1}}(P)},\quad\quad\epsilon^{(2)}=\frac{\tau\sqrt{n(r_{M}+r_{2})(\sigma_{r_{M}+r_{2}}^{2}(W)+p)}}{\sigma_{r_{M}+r_{2}}^{2}(W)}\wedge\sqrt{r_{M}+r_{2}}$
(10)

where $\kappa$, $\tau$ are as in Assumption 2. Then, the estimators satisfy

Joint components: $\displaystyle\mathbb{E}[d(\widehat{M},M)]$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}}}{\delta}[\epsilon^{(1)}+\epsilon^{(2)}]\bigg{)},$
Individual components: $\displaystyle\mathbb{E}[d(\widehat{R}^{(k)},R^{(k)})]$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}r_{k}}}{\delta}[\epsilon^{(1)}+\epsilon^{(2)}]\bigg{)},\hskip
8.5359ptk=1,2.$

From Theorem 3.1, we notice that the error rates are influenced by several
factors. First, the terms $\epsilon^{(1)}$ and $\epsilon^{(2)}$ quantify the
estimation error for the leading eigenvectors and left singular vectors of the
network and the covariates, respectively (Yu et al., 2015; Cai and Zhang,
2018).

Second, the larger is the separation $\delta$ between the individual subspaces
$\mathcal{R}^{(1)}$ and $\mathcal{R}^{(2)}$, the smaller is the estimation
error, matching our intuition. When $\delta$ is close to zero in the presence
of noise, the individual components are easily mistaken for joint.

Theorem 3.1 provides a general error rate that depends on the spectral norm of
the noise in the network adjacency matrix. We now turn our attention to the
more concrete setting of independent Bernoulli noise, which has the GRDPG
model introduced in Section 2.1.1 as a special case. To quantify the noise in
the Bernoulli model, we denote by $\mu(P)$ the largest expected node degree,
that is, $\mu(P)=\max_{i\in[n]}\sum_{j=1}^{n}P_{ij},$ and consider the
following assumption.

######  Assumption 3.

Let $A\in\mathbb{R}^{n\times n}$ be the observed adjacency and
$P\in\mathbb{R}^{n\times n}$ be its expected value as in Equation (1). Then,
the entries $(A_{ij})_{i\leq j}$ of $A$, are independently distributed with
$A_{ij}\sim\text{Ber}(P_{ij})$, and there exists $\rho_{n}>0$ that only
depends on $n$ such that the following conditions hold:

  1. (a)

$\lambda_{r_{M}+r_{1}}(P)=\Omega(\rho_{n}n)$,

  2. (b)

$\rho_{n}=\omega(n^{-1}\log n)$,

  3. (c)

$\mu(P)=\Omega(\rho_{n}n)$.

The constant $\rho_{n}>0$ is a _sparsity factor_ which affects the asymptotic
behavior of the largest expected degree $\mu(P)$ and
$\lambda_{r_{M}+r_{1}}(P)$. Condition (b) ensures a lower bound for the rate
of the sparsity factor so that under (a), the network signal is sufficiently
strong for estimation and under (b) and (c), the network is connected with
high probability (Chung and Lu, 2006). We now formulate a special case of the
above result for this specific Bernoulli model.

######  Corollary 3.2.

Let $A\in\mathbb{R}^{n\times n}$ and $X\in\mathbb{R}^{n\times p}$ be the
observed adjacency and covariate matrices as in Equation (1) and let
$r_{M},r_{1},r_{2}\in\mathbb{N}$ be the joint and individual subspace
dimensions as in Definition 2.1. Define
$(\widehat{M},\widehat{R}^{(1)},\widehat{R}^{(2)})$ to be the output of
Algorithm 1 when supplied $(A,X,r_{M},r_{1},r_{2})$ as arguments. Suppose that
$\sigma_{r_{M}+r_{2}}(W)=\Omega(\sqrt{np})$, $r_{M},r_{1},r_{2}=O(1)$ and
$\delta=\Omega(1)$. Then under Assumptions 1, 2 and 3 we have that

$\displaystyle\mathbb{E}[d(\widehat{M},M)]$
$\displaystyle=O\bigg{(}\frac{1}{\sqrt{\rho_{n}n}}+\frac{\tau}{\sqrt{p}}\bigg{)},$
$\displaystyle\mathbb{E}[d(\widehat{R}^{(k)},R^{(k)})]$
$\displaystyle=O\bigg{(}\frac{1}{\sqrt{\rho_{n}n}}+\frac{\tau}{\sqrt{p}}\bigg{)},\quad\quad
k=1,2.$

Corollary 3.2 provides an explicit characterization of the error rate in terms
of the dimensions and magnitude of the noise. Qualitatively, the network
sparsity ($\rho_{n}$) increases the error in the estimation, as well as the
variance in the covariate data ($\tau$). The assumptions on the signals
$\lambda_{r_{M}+r_{1}}(P)$ and $\sigma_{r_{M}+r_{2}}(W)$ are typical for
signal-plus-noise and network models (Athreya et al., 2017).

##  4 Simulations

In this section, we use simulated data to evaluate the performance of the
proposed method for estimating the joint and individual components, as well as
compare with alternative approaches.

###  4.1 Example: data with group structure

We start with an illustration of the decomposition into joint and individual
components in latent space models with group structure in the nodes,
previously shown in Figure 1. We generate a synthetic dataset with $n=40$
nodes divided into $K=4$ groups. For simplicity, we enumerate the nodes so
that for $k=1,\ldots,4$, nodes $(k-1)n/4+1\ldots,kn/4$ belong to group $k$. We
demonstrate how the methodology can recover joint and individual group
structures in the network and the covariates.

The network data is simulated from a stochastic blockmodel with $K=3$
communities. Formally, the SBM is defined using a connection probability
matrix $B\in[0,1]^{K\times K}$ so that if nodes $i$ and $j$ belong to
communities $k$ and $l$ respectively, then the probability of a connection
between nodes $i$ and $j$ is $P_{ij}=B_{kl}$. This probability matrix is given
by

$B=\begin{pmatrix}0.6&0.05&0.05\\\ 0.05&0.6&0.05\\\
0.05&0.05&0.6\end{pmatrix}.$

The community labels are encoded in a vector $z^{(1)}\in\\{1,\ldots,4\\}^{n}$,
so that $P_{ij}=B_{z^{(1)}_{i}z^{(1)}_{j}}$, and we define these memberships
as

$z_{i}^{(1)}=\left\\{\begin{array}[]{cl}1&i=1,\ldots,10,\\\
2&i=11,\ldots,20\\\ 3&i=21,\ldots,40.\end{array}\right.$

In this way, the network is only able to identify three of the groups in the
nodes, and the nodes from groups 3 and 4 are combined into the same community.

We consider a similar model for the node covariates by assuming the existence
of $L=3$ clusters in a Gaussian mixture model. For each cluster
$l=1,\ldots,L$, there is a corresponding cluster mean
$\mu_{l}\in\mathbb{R}^{p}$. We assume that for each node $i=1,\ldots,n$ the
expected value of its node covariates, $w_{i}=\mathbb{E}(x_{i})$, satisfies
$w_{i}=\mu_{l}$ iff node $i$ belongs to cluster $l$. The rows of the covariate
matrix are generated from a 3-dimensional Gaussian mixture model, where the
cluster means are given by
$\mu_{1}=\begin{pmatrix}-30&-60&30\end{pmatrix}^{\top}$,
$\mu_{2}=\begin{pmatrix}16&8&16\end{pmatrix}^{\top}$,
$\mu_{3}=\begin{pmatrix}-20&40&20\end{pmatrix}^{\top}$. For all clusters, we
use the same covariance matrix equal to the identity. The cluster memberships
are encoded in a vector $z^{(2)}\in\\{1,\ldots,3\\}^{n}$ in a similar way as
in the network, but here the nodes from groups 1 and 2 are combined in the
same cluster, resulting in

$z_{i}^{(2)}=\left\\{\begin{array}[]{cl}1&i=1,\ldots,20,\\\
2&i=21,\ldots,30\\\ 3&i=31,\ldots,40.\end{array}\right.$

Under this SBM model with covariates, each node $i$ has two different
memberships, $z_{i}^{(1)}$ and $z_{i}^{(2)}$. Using our model to identify
joint and individual components, we distinguish common and individual group
structure. In particular, two nodes appear in two different joint clusters if
their memberships are different in both datasets, while two nodes appear in
different individual clusters if they can be separated by one dataset but not
the other. Using Definition 2.1, the dimensions of the model described above
are $r_{M}=2$, $r_{1}=r_{2}=1$. Figure 1(c) shows the initial estimates for
the joint components following our spectral estimation procedure described
previously. The joint components separate the nodes into the two clusters that
are common to the network and the covariates. This aligns with the design of
‘triangle’ cluster in network and ‘green’ cluster in the covariates, namely,
groups 1-2 and 3-4. Figure 1(d) shows the estimates for the individual
components, each one on a different axis. Looking at the x-axis, the values of
$\widehat{R}^{(1)}$ are able to distinguish between clusters 1, 2 and 3-4, as
this information is individual to the network, whereas the values of
$\widehat{R}^{(2)}$ in the y-axis can separate groups 3, 4 and 1-2.

###  4.2 Numerical comparisons

Here we evaluate the performance of our methods and compare with other
embedding techniques for network data and covariates.

####  4.2.1 Data generation

We consider $n=200$ nodes with $p=10$ covariates for each node and use the
model in Definition 2.1 for joint and individual signal components with
$r_{M}=r_{1}=r_{2}=1$. We generate the adjacency matrix
$A\sim\operatorname{RDPG}(\sqrt{\alpha}Y)$ of an undirected network with no
self-loops (zero diagonal) as in Section 2.1.1, where
$Y=\begin{pmatrix}M&R^{(1)}\end{pmatrix}\widetilde{\Gamma}^{(1)}$ is the
latent position matrix and $\alpha$ controls the density of the network to
make the average expected degree of the nodes equal to $20$. For covariates,
we generate $X_{ij}\overset{iid}{\sim}N(W_{ij},\tau^{2})$, with $\tau=0.1$ and
$W=\begin{pmatrix}M&R^{(2)}\end{pmatrix}\Gamma^{(2)}$. We set

$\widetilde{\Gamma}^{(1)}=\begin{pmatrix}\sqrt{nq_{1}}&0\\\
0&\sqrt{nq_{2}}\end{pmatrix},\quad\quad\quad\Gamma^{(2)}=\begin{pmatrix}\sqrt{ns_{1}}&0\\\
0&\sqrt{ns_{2}}\end{pmatrix}Q^{\top},$

where $q_{1},q_{2},s_{1},s_{2}>0$ are parameters that control the strength of
each component, and $Q\in\mathbb{O}_{p,2}$ is a random matrix obtained by
taking the top two right singular vectors from a $p\times p$ matrix of i.i.d.
$N(0,1)$ entries.

To construct $M$, $R^{(1)}$ and $R^{(2)}$, we consider three orthonormal
vectors:

$T^{(0)}=\frac{1}{\sqrt{n}}1_{n},\quad\quad
T^{(1)}=\frac{1}{\sqrt{n}}1_{n/2}\otimes\begin{pmatrix}1\\\
-1\end{pmatrix},\quad\quad
T^{(2)}=\frac{1}{\sqrt{n}}1_{n/4}\otimes\begin{pmatrix}1_{2}\\\
-1_{2}\end{pmatrix},$

where $1_{n}$ is the vector of dimension $n$ with all of its entries equal to
one. We consider two settings in terms of the relative strength of joint and
individual signals by defining the components as follows, where
$\delta=1-\langle R^{(1)},R^{(2)}\rangle$ controls the separation between the
individual $R^{(1)}$ and $R^{(2)}$:

  1. a)

Strong joint signal: $M=T^{(0)}$, $R^{(1)}=T^{(1)}$, $R^{(2)}=\delta
T^{(1)}+\sqrt{1-\delta^{2}}T^{(2)}$, $q_{1}=0.5$, $q_{2}=0.3$, $s_{1}=0.6$.

  2. b)

Weak joint signal: $M=T^{(1)}$, $R^{(1)}=T^{(0)}$, $R^{(2)}=\delta
T^{(0)}+\sqrt{1-\delta^{2}}T^{(2)}$ and $q_{1}=0.2$, $q_{2}=0.6$, $s_{1}=0.2$.

For each setting, we consider how the performance varies as a function of
$\delta$ and as a function of the signal strength of covariates individual
component:

• Varying Individual Subspace Separation:

When the joint signal is strong, we set $s_{2}=0.2$, otherwise $s_{2}=.7$. We
generate the data by varying $\delta$ over a grid of $10$ values in $[0,1]$.

• Varying Individual Component Signal Strength:

We set $\delta=1$ and generate the data by varying $s_{2}$ over a grid of $10$
values in $[0,1]$.

We consider 50 replications for each combination of parameters and report the
average error.

####  4.2.2 Methods for comparison

We consider: (i) the spectral estimator (Spectral) described in Algorithm 1,
(ii) the optimization method from Algorithm 2 initialized at the spectral
estimate obtained from Algorithm 1 (Spectral + Optimization), (iii) the
covariate assisted spectral embedding (CASE) from (Binkiewicz et al., 2017),
(iv) the top $r_{M}$ leading eigenvectors of $A$ (Top_SV_Net) and (v) the
first $r_{M}$ left singular vectors of $X$ (Top_SV_Cov). Methods (i) and (ii)
return estimates of three different components (joint, network individual, and
covariate individual), whereas methods (iii)-(v) only return a single set of
components. While methods (iii)-(v) do not intend to specifically estimate the
joint and individual components of our model, we use them as a baseline since
they still capture the true components in the low-rank structure of the data
for various parameter settings. We compare the output of each such method to
the true joint and individual components separately to investigate under which
conditions the components recovered by these methods contain signal specific
to the joint or individual components. Explicitly, the components obtained
from method (iii) are given by
$\operatorname{SV}(L_{\tau}^{2}+\alpha_{0}XX^{\top},1)$ where $L_{\tau}$ is
the regularized graph Laplacian and $\alpha_{0}$ is the initial choice of
tuning parameter given in Binkiewicz et al. (2017). The components obtained
from methods (iv) and (v) are given by $\operatorname{SV}(A,1)$ and
$\operatorname{SV}(X,1)$ respectively.

To evaluate the performance, we calculate the Procrustes distance from
Equation (8) between the estimated components and the truth separately for
joint components $M$ and individual components $R^{(1)}$, $R^{(2)}$. Note that
for matrices with orthonormal columns $U,V\in\mathbb{O}_{n,d}$, this distance
is bounded by $\sqrt{2d}$, so in our simulations this bound is $\sqrt{2}$.

####  4.2.3 Results

Figures 2-3 show the average Procrustes distance over 50 replications for each
method, setting and component type. In general, our Spectral and Spectral +
Optimization methods show good performance in estimating the correct
components when the signal is sufficiently strong. For instance, we observe in
Figure 2 that the Procrustes error is low whenever the separation between the
individual components is sufficiently large (top row) or the signal in the
joint components is strong (bottom row). In Figure 3, our methods perform
generally well except for the estimation of the covariate individual
components when the signal in this component is low, which is reasonable.
These results confirm our observations from Theorem 3.1 regarding the effect
of the $\delta$ and the magnitude of the eigenvalues and singular values of
the signal matrices. While both the Spectral and Spectral + Optimization
methods show similar performance, we see that in cases where the joint
component signal is weak, the optimization method gives an improvement to the
results, particularly noticeable in the estimation of $M$.

The alternative approaches (CASE, Top_SV_Net and Top_SV_Cov) aim to estimate
low-rank components that capture the most variability in a given matrix, and
as such, they do not differentiate between joint and individual components.
Nevertheless, they coincide in some cases with the joint or individual
components, which allows us to compare with our methods. When the joint
components have strong signal (bottom row of each figure), all of these
methods recover $M$ quite well (left column), and our methods show similar
performance. On the other hand, when the joint signal is weak (top row of each
figure), Top_SV_Net and Top_SV_Cov estimate the individual components
$R^{(1)}$ and $R^{(2)}$, respectively (middle and right columns), and our
methods once again show similar performance as long as the subspace separation
is large enough. CASE aims to capture the joint signal, but the performance
depends on how strong the individual signals are. For instance, CASE generally
estimates the individual components when the joint signal is weak (top row of
each figure), and the method would either capture $R^{(1)}$ or $R^{(2)}$
depending on which component has the strongest signal. In contrast, our
methods can always differentiate between the joint or individual components
regardless of the signal strength.

Figure 2: Joint and individual component estimate error (measured by
Procrustes distance) averaged over 50 replications as a function of the
separation between the individual subspaces ($\delta$) while maintaining the
individual component signal strength constant. Figure 3: Joint and individual
component estimate error (measured by Procrustes distance) averaged over 50
replications as a function of the covariate individual component signal
strength $(s_{2})$ while maintaining the individual components to be
orthogonal.

##  5 Application to trade network data with covariates

We illustrate the methodology in the analysis of the food trade network
between countries reported by the Food and Agriculture Organization of the
United Nations (FAO). In addition to the network, we consider country-level
economic and development indicators from the World Bank open data database
that could explain part of the variation in the trade data. We use our
methodology to explore the joint and individual structure in the trade data
and the covariates, with the goal of understanding trading patterns that can
be explained by country-level indicators, as well as additional information
that is unique to each dataset.

Data description and processing. The version of the FAO trade data we used was
collected by De Domenico et al. (2015) and represents the total food trade
from one country to another of a given food product during the year 2010. For
the node covariate data, we collected economic and development country data
from the same year reported in the World Bank open data database through the
wbstats package (Piburn, 2020).

The original network data reports multiplex-directed food trade between
countries. Each of the 364 layers of the multiplex network corresponds to a
food product. Nodes correspond to countries and edges correspond to the import
and export relations between them for the associated food product, weighted by
the traded quantity of the particular product. We aggregate the adjacency
matrices over all layers to obtain a single weighted and directed adjacency
matrix representing the total food import/export trade between countries. We
then symmetrize this adjacency matrix by adding its transpose to itself.
Finally, the weights are log-transformed by $\log(x+1)$.

The covariate data consists of the following economic and developmental
indicators: (i) gross domestic product per capita in US dollars (GDP), (ii)
agricultural, forestry and fishing value added as a percentage of the GDP
(AFF), (iii) total natural resource rents as a percentage of the GDP (NAT),
(iv) government expenditure on education as a percentage of the GDP (EDU) and
(v) crude birth rate per 1000 people (BIR). The values of these covariates
were log-transformed similarly to the trade network weights. In addition,
categorical covariates indicating the continents to which each country belongs
(Africa, Europe, the Americas, Asia and Oceania) are represented via
additional columns in the covariate matrix with dummy variables. The columns
of the covariate matrix were then centered and standardized.

As the countries with available data from the FAO data are not the same as
those with available data from the World Bank dataset, we only consider the
146 countries with data present in both. In summary, we work with a symmetric,
weighted adjacency matrix $A\in\mathbb{R}^{146\times 146}$ and a covariate
matrix $X\in\mathbb{R}^{146\times 10}$.

Data analysis. We begin with an exploratory analysis of each dataset
independently by constructing separate node embeddings for the network and the
covariates. For the network, we compute the adjacency spectral embedding (ASE)
(Sussman et al., 2012) into four dimensions, which is defined as
$\widehat{Y}=\widehat{V}^{(1)}|\widehat{\Lambda}|^{1/2}$, where
$\widehat{V}^{(1)}=\operatorname{EIG}(A,4)$ is the matrix of leading
eigenvectors of $A$ and
$\Lambda=\text{diag}(\lambda_{1}(A),\ldots,\lambda_{4}(A))$ are its leading
eigenvalues. The rows of the ASE are a consistent estimator for the latent
positions of the GRDPG model defined in Section 2.1.1 (Rubin-Delanchy et al.,
2022). For the covariate data, we obtain the top four left leading singular
vectors, which are equivalent to the scores of the principal components.
Ideally, any nontrivial structure should be identified by these components as
they capture most of the variation in the data. We display these plots in
Figure 4. Upon inspection, the network embeddings in Figure 4(a) only display
part of the geographic structure related to the region (continents). In
contrast, the covariate principal components display a clustering structure
according to continent in Figure 4(b), as these indicator variables are
included in the covariate data. The PCA biplot shows two main directions of
variability corresponding to the size of the economy (GDP per capita and other
variables correlated to it), and the expenditure in education (EDU).

We now estimate joint and individual components from the data. To estimate the
ranks in the model $\operatorname{rank}(P)=r_{M}+r_{1}$ and
$\operatorname{rank}(W)=r_{M}+r_{2}$, we first estimate the embedding
dimensions via the elbow method on the scree plots of the singular values of
the adjacency and covariate matrices, resulting in values
$\operatorname{rank}(P)=4$ and $\operatorname{rank}(W)=4$. For the dimensions
of the joint and individual components, we first estimate the dimension of the
joint structure based on the scree plot of the singular values of the matrix
$\widehat{U}$ obtained in Algorithm 1 and look for an elbow. This results in
an estimated value of $r_{M}=2$, and therefore, we set $r_{1}=2$ and
$r_{2}=2$.

(a) ASE of the trade network

(b) Biplot of the covariate principal components

Figure 4: Network and covariate embeddings, constructed by separately
computing the network adjacency spectral embedding (ASE) and the covariate
principal components. Each point represents a different country (node),
colored by continent and size proportional to GDP per capita.

The joint and individual components are first estimated with the spectral
method (Algorithm 1), and these are used to initialize the optimization method
(Algorithm 2). The joint components obtained by each of these methods are
shown in Figures 5(a) and 5(b). While both methods show a very similar
pattern, the joint components from the optimization method show tighter
groupings than those obtained from the spectral method. We observed a similar
result in the individual components, and hence, we only report the results for
the network and covariate individual components obtained from the optimization
method in Figure 6.

(a) Spectral method

(b) Optimization method

Figure 5: Plots of estimated joint components ($\widehat{M}$) from Algorithm 1
(left) and Algorithm 2 (right). Each point represents a different country
(node), colored by continent and size corresponding to GDP value.

(a) Network individual components

(b) Covariate individual components

Figure 6: Plots of estimated individual components for the network
($\widehat{R}^{(1)}$) and covariates ($\widehat{R}^{(2)}$) using Algorithm 2.

(a) Europe

(b) Americas

Figure 7: Estimated joint components for countries in Europe and the Americas.
Each point size is proportional to the country GDP per capita.

(a) Europe

(b) Americas

Figure 8: Estimated network individual components for countries in Europe and
the Americas. Each point size is proportional to the country GDP per capita.

(a) total trade inside vs total trade outside Europe for European countries.

(b) total trade inside vs total trade outside the Americas for American
countries.

Figure 9: Total trade volume per country divided into trade with other
countries inside (x-axis) and outside (y-axis) the same region. The blue line
corresponds to the fitted regression value of a linear model without an
intercept; countries below this line have higher local trade volume than
average and vice-versa. The location of the countries in this figure closely
resembles the relative location in the individual network components (Figure
8).

To gain more insight into the information displayed in the estimated
components, we focus on the countries from specific regions, in particular,
the Americas and Europe. Figures 7 and 8 show the estimated joint and
individual components only for the countries in those regions, together with
the country names and the point size proportional to the GDP per capita. In
the joint components (Figure 7), the scatterplots show an association between
GPD per capita and the joint latent positions, as countries with higher GDP
are located closer to the bottom left of the corresponding plots. This
observation suggests that this covariate plays an important role in explaining
trade between countries. The individual components, on the other hand, do not
capture a clear association with GDP, as this was already captured in the
joint.

The individual components in Figure 8 reveal trading structure not captured by
the covariates. A further analysis of the data reveals that the components
capture global and local trade variability, as countries located closer to the
right tend to trade more with countries from the same region, whereas
countries on the left side have higher trade volume with countries from a
different region. This trend is confirmed in Figure 9, which shows the total
trade volume of each country with other countries from the same region against
the total trade volume with countries from different regions. For each country
$i$, these quantities are defined as

$w^{\text{inside}}_{i}=\sum_{j:\text{region}(i)=\text{region}(j)}A_{ij},\quad\quad\quad\quad
w^{\text{outside}}_{i}=\sum_{j:\text{region}(i)\neq\text{region}(j)}A_{ij}.$

Comparing Figure 8(b) and Figure 8(a), we see that the network individual
components are capturing much of the same information as the total trade
volume inside and outside of the respective region. Moreover, this information
is not associated with GDP per capita, as countries with high GDP appear on
both sides of the figure.

To gain more interpretation of the role of the covariates in explaining the
components, we decompose the covariate matrix by projecting the data onto the
joint and individual subspaces as

$X=\mathcal{P}_{\widehat{M}}X+\mathcal{P}_{\widehat{R}^{(2)}}X+(I-\mathcal{P}_{\widehat{M}}-\mathcal{P}_{\widehat{R}^{(2)}})X.$
(11)

For the joint and individual projections, we do PCA to obtain the loadings of
each of the covariates onto the joint and individual components. Figures 10(a)
and 10(b) show the biplots representing the country scores obtained from each
of the projections, together with the variable loadings (arrows in the
figure). We see in Figure 10(a) that there is a large correlation between the
economic covariates (GDP, BIR, AFF) and the first principal component in the
joint PCA, suggesting that the effect captured by both the network and
covariate data is related to these factors. This observation coincides with
the results in Figure 7. On the other hand, the loadings for education are low
in the joint components, which is reasonable given that we do not expect that
government expenditure on education would impact the trading network. The
projection of the covariate data onto the individual components captures the
remaining variability with respect to the economic factors, and the loadings
for education and natural resources show higher values, indicating that these
variables do not have a strong impact on the network trading patterns.

(a) PCA biplot for $\mathcal{P}_{\widehat{M}}X$

(b) PCA biplot for $\mathcal{P}_{\widehat{R}^{(2)}}X$

Figure 10: PCA biplots for the joint and individual parts of the covariates

Finally, we measure the proportion of variation explained by the joint and
individual components. For the covariates, we use the decomposition in
Equation (11) to calculate the proportion of variation in the covariates
explained by the joint, individual and residual components, denoted
$\operatorname*{Var}_{J}(X)$, $\operatorname*{Var}_{I}(X)$ and
$\text{Var}_{\text{res}}(X)$ respectively, as

$\text{Var}_{J}(X)=\frac{\|\mathcal{P}_{\mathcal{C}(\widehat{M})}X\|_{F}^{2}}{\|X\|_{F}^{2}},\quad\quad\text{Var}_{I}(X)=\frac{\|\mathcal{P}_{\mathcal{C}(\widehat{R}^{(2)})}X\|_{F}^{2}}{\|X\|_{F}^{2}},\quad\quad\text{Var}_{\text{res}}(X)=1-(\text{Var}_{J}(X)+\text{Var}_{I}(X)).$

For the network adjacency matrix, we first decompose the variation into parts
explained by the signal and the residual as

$A=\widehat{P}+(A-\widehat{P}),$

where
$\widehat{P}=\mathcal{P}_{(\widehat{M},\widehat{R}^{(1)})}A\mathcal{P}_{(\widehat{M},\widehat{R}^{(1)})}^{\top}$
is the projection of $A$ onto the joint and individual components. We denote
the proportion of signal and residual variation as
$\text{Var}_{\widehat{P}}(A)$ and $\text{Var}_{\text{res}}(A)$ respectively,
defined by

$\text{Var}_{\widehat{P}}(A)=\frac{\|\widehat{P}\|_{F}^{2}}{\|A\|_{F}^{2}},\quad\quad\text{Var}_{\text{res}}(A)=1-\text{Var}_{\Sigma}(A).$

We then define the proportion of variation in the network signal explained by
the joint and individual components in terms of the latent positions of the
model as in Equation (5), denoted $\text{Var}_{J|\widehat{P}}(A)$ and
$\text{Var}_{I|\widehat{P}}(A)$, respectively. Specifically, let
$\widehat{Y}=\operatorname{ASE}(\widehat{P},4)$ be an estimate of the latent
positions of the model. We set

$\quad\text{Var}_{J|\widehat{P}}(A)=\frac{\|\mathcal{P}_{\widehat{M}}\widehat{Y}\|^{2}_{F}}{\|\widehat{Y}\|^{2}_{F}},\quad\quad\quad\quad\text{Var}_{I|\widehat{P}}(A)=\frac{\|\mathcal{P}_{\widehat{R}^{(1)}}\widehat{Y}\|^{2}_{F}}{\|\widehat{Y}\|^{2}_{F}}$

respectively. We then define the proportion of variation in the network
explained by the joint and individual components, denoted $\text{Var}_{J}(A)$
and $\text{Var}_{J}(A)$ respectively, by

$\text{Var}_{J}(A)=\text{Var}_{J|\widehat{P}}(A)\text{Var}_{\widehat{P}}(A),\quad\quad\text{Var}_{I}(A)=\text{Var}_{I|\widehat{P}}(A)\text{Var}_{\widehat{P}}(A).$

The proportions of variation explained are reported in Table 1. Overall, the
low-rank models capture a large portion of the variation in the data (88% for
the network and 79.75% for the covariates). However, when looking at the joint
and individual components, the proportions vary significantly in the network
and the covariates. For instance, while most of the variation in the
covariates (51.59%) is captured by the joint components, these only capture
15.6% of the variation in the network. This observation suggests that most of
the variation in the covariates is associated with the network trading
structure, but this information is not enough to capture the variation in the
trading patterns.

Dataset | Joint | Individual | Residual
---|---|---|---
Network | 15.60% | 72.40% | 12.00%
Covariates | 51.59% | 28.16% | 20.25%
Table 1: Proportion of variation explained by component for network and
covariate datasets

##  6 Discussion

In this work, we have introduced a model for joint and individual information
present in network data with node covariates, developed an estimation
procedure to obtain joint and individual components of the model, and
established theoretical guarantees under a general signal-plus-noise model for
the covariates. Illustrations in simulated and real networks show the ability
of the methodology to recover meaningful components from the data.

Even though the theory assumes an undirected network with independent
Bernoulli entries, similar theoretical guarantees can be obtained when
imposing other distributional assumptions on $A$.

This work is limited to studying only one network and its node covariates. In
many situations, more complex data types are also available, such as multiplex
or multilayer networks with their corresponding node covariates, or even edge-
covariate information. Having access to multiple datasets of this form can
yield improvements in the estimation of the joint components. Extensions of
the methodology and theory to these scenarios are left as future work. In this
paper, we limit our theoretical analysis to studying the estimation error in
Frobenius norm of the joint and individual embeddings. Many applications
require further subsequent inferences, including node clustering, prediction
or hypothesis testing, for which a more refined analysis of the statistical
error and its distributional properties are required, which we leave as future
work.

## Acknowledgement

IG was supported by NSF DMS 2422478.

##  Appendix A Proofs of main results

###  A.1 Identifiability

###### Proof.

[Proof of Lemma 2.1] We recall that, by definition,
$\mathcal{M}=\mathcal{C}(P)\cap\mathcal{C}(W)$,
$\mathcal{R}^{(1)}=\mathcal{C}(\mathcal{P}_{\mathcal{M}^{\perp}}P)$,
$\mathcal{R}^{(2)}=\mathcal{C}(\mathcal{P}_{\mathcal{M}^{\perp}}W)$,
$r_{M}=\dim(\mathcal{M})$ and $r_{k}=\dim(\mathcal{R}^{(k)})$. Therefore, by
construction, $\operatorname{rank}(P)=r_{M}+r_{1}$ and
$\operatorname{rank}(W)=r_{M}+r_{2}$. Now, suppose that for $k=1,2$,
$\mathcal{M},\mathcal{R}^{(k)}\neq\\{0\\}$.

  1. (a)

Existence:
Define $J^{(1)}=P_{\mathcal{M}}P$, $J^{(2)}=P_{\mathcal{M}}W$,
$L^{(1)}=\mathcal{P}_{\mathcal{M}^{\perp}}P$ and
$L^{(2)}=\mathcal{P}_{\mathcal{M}^{\perp}}W$. By definition,
$\mathcal{R}^{(k)}=\mathcal{C}(L^{(k)})$. Then for $k=1,2$,

    1. (1)

$P=J^{(1)}+L^{(1)}$ and $W=J^{(2)}+L^{(2)}$ ,

    2. (2)

$\mathcal{C}(J^{(k)})=\mathcal{M}$ and
$\mathcal{M}\subset\mathcal{C}(P)\cap\mathcal{C}(W)$,

    3. (3)

$\mathcal{M}\perp\mathcal{C}(L^{(k)})$,

    4. (4)

$\mathcal{C}(L^{(1)})\cap\mathcal{C}(L^{(2)})=\\{0\\}$.

The proof of Lemma 1 in Yuan and Gaynanova (2022) implies that the properties
$(1)-(4)$ uniquely identify $J^{(1)}$, $J^{(2)}$, $L^{(1)}$ and $L^{(2)}$.
Since $\mathcal{M}\neq\\{0\\}$, there exists $M\in\mathbb{O}_{n,r_{M}}$ such
that $\mathcal{C}(M)=\mathcal{M}$. Since $\mathcal{R}^{(k)}\neq\\{0\\}$, there
exists $R^{(k)}\in\mathbb{O}_{n,r_{k}}$ such that
$\mathcal{C}(R^{(k)})=\mathcal{R}_{k}$. We set
$\Gamma^{(k)}_{1}=M^{\top}J^{(k)}$, $\Gamma^{(k)}_{2}=(R^{(k)})^{\top}L^{(k)}$
and $\Gamma^{(k)}=\begin{pmatrix}\Gamma^{(k)}_{1}\\\
\Gamma^{(k)}_{2}\end{pmatrix}$. Then

$\displaystyle\begin{pmatrix}M&R^{(k)}\end{pmatrix}\Gamma^{(k)}$
$\displaystyle=\begin{pmatrix}M&R^{(k)}\end{pmatrix}\begin{pmatrix}\Gamma^{(k)}_{1}\\\
\Gamma^{(k)}_{2}\end{pmatrix}$
$\displaystyle=M\Gamma^{(k)}_{1}+R^{(k)}\Gamma^{(k)}_{2}$
$\displaystyle=MM^{\top}J^{(k)}+R^{(k)}(R^{(k)})^{\top}L^{(k)}$
$\displaystyle=\mathcal{P}_{\mathcal{M}}J^{(k)}+\mathcal{P}_{\mathcal{R}^{(k)}}L^{(k)}$
$\displaystyle=J^{(k)}+L^{(k)}$ $\displaystyle=\begin{cases}P,&k=1,\\\
W,&k=2.\end{cases}$

Since $\begin{pmatrix}M&R^{(1)}\end{pmatrix}\Gamma^{(1)}=P$ and
$\mathcal{C}(P)=\mathcal{C}\begin{pmatrix}M&R^{(1)}\end{pmatrix}$, we have
that

$\displaystyle\operatorname{rank}(\Gamma^{(1)})$
$\displaystyle=\operatorname{rank}[\begin{pmatrix}M&R^{(1)}\end{pmatrix}^{\top}P]$
$\displaystyle\geq\operatorname{rank}[\begin{pmatrix}M&R^{(1)}\end{pmatrix}\begin{pmatrix}M&R^{(1)}\end{pmatrix}^{\top}P]$
$\displaystyle=\operatorname{rank}(\mathcal{P}_{\mathcal{C}(P)}P)$
$\displaystyle=r_{M}+r_{1}.$

On the other hand, $\Gamma^{(1)}\in\mathbb{R}^{(r_{M}+r_{1})\times n}$. Hence
$\operatorname{rank}(\Gamma^{(1)})\leq r_{M}+r_{1}$. Thus
$\operatorname{rank}(\Gamma^{(1)})=r_{M}+r_{1}$. Similarly,
$\operatorname{rank}(\Gamma^{(2)})=r_{M}+r_{2}$.

  2. (b)

Orthogonality:
Since $\mathcal{C}(R^{(k)})\subset\mathcal{C}(M)^{\perp}$, we have that
$M\perp R^{(k)}$.

  3. (c)

Uniqueness:
Suppose that there are $N\in\mathbb{O}_{n\times r_{M}}$,
$S^{(k)}\in\mathbb{O}_{n\times r_{k}}$, $k=1,2$, such that

$\displaystyle\mathcal{C}(N)=\mathcal{M},$
$\displaystyle\mathcal{C}(S^{(k)})=\mathcal{R}^{(k)}.$

Since $\mathcal{C}(M)=\mathcal{M}$, there exists full rank
$W\in\mathbb{R}^{r_{M}\times r_{M}}$ such that $M=NW$. This implies that
$N^{\top}M=W$, and $I_{r_{M}}=M^{\top}NW$. Therefore,
$W^{-1}=M^{\top}N=(N^{\top}M)^{\top}=W^{\top}$, that is, the inverse of $W$ is
its transpose, which implies that $W$ is an orthogonal matrix. Hence, $M$ is
unique up to an orthogonal matrix right multiplication. By similar arguments,
$S^{(k)}=R^{(k)}W^{\prime}$ for some $W^{\prime}\in\mathbb{O}_{r_{k}}$, so
$R^{(k)}$ is unique up to right multiplication by an orthogonal matrix.

∎

###  A.2 Spectral method solution

###### Proof.

[Proof of Proposition 2.3] Suppose that $P$ and $W$ satisfy the assumptions in
Lemma 2.1 with joint and individual components $(M,R^{(1)},R^{(2)})$ and ranks
$(r_{M},r_{1},r_{2})$. Define
$(\widehat{M},\widehat{R}^{(1)},\widehat{R}^{(2)})$ to be the output of
Algorithm 1 when supplied $(P,W,r_{M},r_{1},r_{2})$ as arguments. Lemma 2.1
implies that there exist $\Gamma^{(1)}\in\mathbb{R}^{(r_{M}+r_{1})\times n}$
and $\Gamma^{(2)}\in\mathbb{R}^{(r_{M}+r_{2})\times p}$ such that

$P=\begin{pmatrix}M&R^{(1)}\end{pmatrix}\Gamma^{(1)},\quad
W=\begin{pmatrix}M&R^{(2)}\end{pmatrix}\Gamma^{(2)}.$

Let $\widehat{V}^{(1)},\widehat{U}$ be the matrices in the intermediate steps
of Algorithm 1. Then, by construction,
$\widehat{V}^{(1)}\cong_{O}\begin{pmatrix}M&R^{(1)}\end{pmatrix}$ and
$\widehat{V}^{(2)}\cong_{O}\begin{pmatrix}M&R^{(2)}\end{pmatrix}$. Set
$R=\begin{pmatrix}R^{(1)}&R^{(2)}\end{pmatrix}$. By construction, this
$n\times(r_{1}+r_{2})$ matrix has full rank, and write the eigenvalue
decomposition of $RR^{\top}$ as $RR^{\top}=LDL^{\top}$, where
$L\in\mathbb{O}_{n,r_{1}+r_{2}}$ and
$D\in\mathbb{R}_{(r_{1}+r_{2})\times(r_{1}+r_{2})}$ is the diagonal matrix
containing the non-zero eigenvalues. Since $M\perp R^{(i)}$, we have that
$M\perp L$ and therefore $\widehat{U}\widehat{U}^{\top}$ has eigenvalue
decomposition

$\displaystyle\widehat{U}\widehat{U}^{\top}$
$\displaystyle=\widehat{V}^{(1)}(\widehat{V}^{(1)})^{\top}+\widehat{V}^{(2)}(\widehat{V}^{(2)})^{\top}$
$\displaystyle=2MM^{\top}+RR^{\top}$
$\displaystyle=\begin{pmatrix}M&L\end{pmatrix}\begin{pmatrix}2I_{r_{M}}&0\\\
0&D\end{pmatrix}\begin{pmatrix}M&L\end{pmatrix}^{\top}.$

Let $\delta=1-\sigma_{1}((R^{(1)})^{\top}R^{(2)})$ as in Equation (9). Denote
the leading eigenvalue of $D$ in magnitude by $\lambda_{1}(D)$. We note that
we have

$\displaystyle\lambda_{1}(D)$ $\displaystyle=\lambda_{1}(RR^{\top})$
$\displaystyle=\lambda_{1}(R^{\top}R)$
$\displaystyle=\lambda_{1}\bigg{[}\begin{pmatrix}(R^{(1)})^{\top}R^{(1)}&(R^{(1)})^{\top}R^{(2)}\\\
(R^{(2)})^{\top}R^{(1)}&(R^{(2)})^{\top}R^{(2)}\\\ \end{pmatrix}\bigg{]}$
$\displaystyle=\lambda_{1}\bigg{[}\begin{pmatrix}I_{r_{1}}&(R^{(1)})^{\top}R^{(2)}\\\
(R^{(2)})^{\top}R^{(1)}&I_{r_{2}}\\\ \end{pmatrix}\bigg{]}$
$\displaystyle=\lambda_{1}\bigg{[}I_{r_{1}+r_{2}}+\begin{pmatrix}0&(R^{(1)})^{\top}R^{(2)}\\\
(R^{(2)})^{\top}R^{(1)}&0\\\ \end{pmatrix}\bigg{]}$
$\displaystyle=1+\lambda_{1}\bigg{[}\begin{pmatrix}0&(R^{(1)})^{\top}R^{(2)}\\\
(R^{(2)})^{\top}R^{(1)}&0\\\ \end{pmatrix}\bigg{]}.$

Exercise $1.3.17$ in Tao (2023) implies that

$\displaystyle
1+\lambda_{1}\bigg{[}\begin{pmatrix}0&(R^{(1)})^{\top}R^{(2)}\\\
(R^{(2)})^{\top}R^{(1)}&0\\\ \end{pmatrix}\bigg{]}$ $\displaystyle\leq
1+\sigma_{1}[(R^{(2)})^{\top}R^{(1)}]$ $\displaystyle=1+\cos(\theta)$
$\displaystyle=2-\delta.$

Since $\mathcal{C}(\mathcal{R}^{(1)})\cap\mathcal{C}(R^{(1)})=0$ and
$r_{M},r_{1},r_{2}>0$, we have that $\delta>0$. Therefore $\lambda_{1}(D)<2$
and thus

$\displaystyle\widehat{M}$
$\displaystyle=\operatorname{SV}(\widehat{U},r_{M})$
$\displaystyle\cong_{O}\operatorname{EIG}(\widehat{U}\widehat{U}^{\top},r_{M})$
$\displaystyle\cong_{O}M.$

Then, we have that

$\displaystyle\widehat{R}^{(1)}$
$\displaystyle=\operatorname{SV}((I-\widehat{M}\widehat{M}^{\top})\widehat{V}^{(1)},r_{1})$
$\displaystyle=\operatorname{SV}(R^{(1)},r_{1})$
$\displaystyle\cong_{O}R^{(1)}.$

Similarly, $\widehat{R}^{(2)}\cong_{O}R^{(2)}$. ∎

###  A.3 Optimization method solution

###### Proof.

[Proof of Lemma 2.4] Writing $V^{(k)}=\begin{pmatrix}M&R^{(k)}\end{pmatrix}$,
the minimization in (6) can be rewritten as

$\displaystyle\operatorname*{minimize}_{M\in\mathbb{O}_{n,r_{M}}}$
$\displaystyle\\{\|A^{\prime}-\mathcal{P}_{\mathcal{C}(V^{(1)})}A^{\prime}\|_{F}^{2}+\|X-\mathcal{P}_{\mathcal{C}(V^{(2)})}X\|_{F}^{2}\\}$
(12) such that $\displaystyle\quad\quad
V^{(k)}=\begin{pmatrix}M&R^{(k)}\end{pmatrix}\quad k=1,2.$

Equivalently, we only need to minimize

$\|\mathcal{P}_{\mathcal{C}(V^{(1)})^{\perp}}A^{\prime}\|_{F}^{2}+\|\mathcal{P}_{\mathcal{C}(V^{(2)})^{\perp}}X\|_{F}^{2}.$

Due to orthogonality of $R^{(k)}$ and $M$, this is equivalent to minimizing

$\|\mathcal{P}_{\mathcal{C}(R^{(1)})^{\perp}}A^{\prime}-\mathcal{P}_{\mathcal{C}(M)}\mathcal{P}_{\mathcal{C}(R^{(1)})^{\perp}}A^{\prime}\|_{F}^{2}+\|\mathcal{P}_{\mathcal{C}(R^{(2)})^{\perp}}X-\mathcal{P}_{\mathcal{C}(M)}\mathcal{P}_{\mathcal{C}(R^{(2)})^{\perp}}X\|_{F}^{2}.$

Set

$Y=\begin{pmatrix}\mathcal{P}_{\mathcal{C}(R^{(1)})^{\perp}}A^{\prime}&\mathcal{P}_{\mathcal{C}(R^{(2)})^{\perp}}X\end{pmatrix},\quad\quad
R=\begin{pmatrix}R^{(1)}&R^{(2)}\end{pmatrix}.$

Then the minimization problem can be equivalently written as

$\displaystyle\operatorname*{minimize}_{M\in\mathbb{O}_{n,r_{M}}}$
$\displaystyle\left\\{{\|Y-\mathcal{P}_{\mathcal{C}(M)}Y\|_{F}^{2}}\right\\}$
s.t. $\displaystyle\quad M^{\top}R=0.$

Because $M$ is orthogonal to $R$, it is equivalent to minimize

$\|\mathcal{P}_{\mathcal{C}(R)^{\perp}}Y-\mathcal{P}_{\mathcal{C}(M)}(\mathcal{P}_{\mathcal{C}(R)^{\perp}}Y)\|_{F}^{2}.$

The Eckart-Young-Mirsky theorem implies that the optimal $\widehat{M}$ is
given by
$\widehat{M}=\operatorname{SV}(\mathcal{P}_{\mathcal{C}(R)^{\perp}}Y,r_{M})$.
∎

###  A.4 Error rate of spectral method

###### Proof.

[Proof of Theorem 3.1]
We split the proof intro three parts, covering the error rate for each of the
joint, network individual and covariate individual components.

  * •

Rate for Joint Components:
We recall from Definition 5.13 in Vershynin (2010), the $\psi_{1}$-norm of a
random variable $T$ is given by $\|T\|_{\psi_{1}}=\sup\limits_{q\geq
1}[q^{-1}(\mathbb{E}|T|^{q})^{1/q}]$. We recall that

$\displaystyle\epsilon^{(1)}=\frac{\sqrt{r_{M}+r_{1}}\kappa}{\lambda_{r_{M}+r_{1}}(P)},$
$\displaystyle\epsilon^{(2)}=\frac{\tau\sqrt{n(\sigma_{r_{M}+r_{2}}^{2}(W)+p)(r_{M}+r_{2})}}{\sigma_{r_{M}+r_{2}}^{2}(W)}\wedge\sqrt{r_{M}+r_{2}},$
$\displaystyle\delta=1-\cos(\theta),$

where $r_{M}=\dim(\mathcal{M})$, $r_{k}=\dim(\mathcal{R}^{(k)})$ and $\theta$
is the first principal angle between $\mathcal{C}(R^{(1)})$ and
$\mathcal{C}(R^{(2)})$ so that
$\cos(\theta)=\sigma_{1}((R^{(2)})^{\top}R^{(1)})$. According to Lemma 2.1,
there exist $M\in\mathbb{R}^{n\times r_{M}}$, $R^{(1)}\in\mathbb{R}^{n\times
r_{1}}$, $R^{(2)}\in\mathbb{R}^{n\times r_{2}}$,
$\Gamma^{(1)}\in\mathbb{R}^{(r_{M}+r_{1})\times n}$ and
$\Gamma^{(2)}\in\mathbb{R}^{(r_{M}+r_{2})\times p}$ such that

$\displaystyle P=\begin{pmatrix}M&R^{(1)}\end{pmatrix}\Gamma^{(1)},$
$\displaystyle W=\begin{pmatrix}M&R^{(2)}\end{pmatrix}\Gamma^{(2)},$ (13)

$M\perp R^{(k)}$ and $M$, $\mathcal{R}^{(1)}$, $R^{(2)}$ are unique modulo
$\cong_{O}$. For notational convenience, we set $A^{(1)}=A$, $A^{(2)}=X$ and
$P^{(k)}=\mathbb{E}(A^{(k)})$. Following Algorithm 1 we define

$\displaystyle\widehat{V}^{(k)}=\operatorname{SV}(A^{(k)},r_{M}+r_{i}),$
$\displaystyle V^{(k)}=\operatorname{SV}(P^{(k)},r_{M}+r_{i}),$
$\displaystyle\widehat{U}=\begin{pmatrix}\widehat{V}^{(1)}&\widehat{V}^{(2)}\end{pmatrix},$
$\displaystyle U=\begin{pmatrix}V^{(1)}&V^{(2)}\end{pmatrix},$
$\displaystyle\widehat{M}=\operatorname{SV}(\widehat{U},r_{M}),$
$\displaystyle\widehat{R}^{(k)}=\operatorname{SV}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\widehat{V}^{(k)},r_{k}).$

From the proof of Proposition 2.3, we know that
$V^{(k)}\cong_{O}\begin{pmatrix}M&R^{(k)}\end{pmatrix}$,
$M\cong_{O}\operatorname{SV}(U,r_{M})$ and
$R^{(k)}\cong_{O}\operatorname{SV}(\mathcal{P}_{\mathcal{M}^{\perp}}V^{(k)},r_{k})$.
Since we need only consider the projection matrix onto a components respective
column space and such a projection matrix is invariant under right
multiplication by an orthogonal matrix, we may assume that
$V^{(k)}=\begin{pmatrix}M&R^{(k)}\end{pmatrix}$,
$M=\operatorname{SV}(U,r_{M})$ and
$R^{(k)}=\operatorname{SV}(\mathcal{P}_{\mathcal{M}^{\perp}}V^{(k)},r_{k})$ so
that $P=V^{(1)}\Gamma^{(1)}$ and $W=V^{(2)}\Gamma^{(2)}$. We will work with
corresponding projection matrices

$\displaystyle\mathcal{P}_{\widehat{\mathcal{V}}^{(k)}}=\widehat{V}^{(k)}(\widehat{V}^{(k)})^{\top},$
$\displaystyle\mathcal{P}_{\mathcal{V}^{(k)}}=V^{(k)}(V^{(k)})^{\top},$
$\displaystyle\mathcal{P}_{\widehat{\mathcal{M}}}=\widehat{M}\widehat{M}^{\top},$
$\displaystyle\mathcal{P}_{\mathcal{M}}=MM^{\top},$
$\displaystyle\mathcal{P}_{\widehat{\mathcal{R}}^{(k)}}=\widehat{R}^{(k)}(\widehat{R}^{(k)})^{\top},$
$\displaystyle\mathcal{P}_{\mathcal{R}^{(k)}}=R^{(k)}(R^{(k)})^{\top},$

and we define

$\mathcal{Q}_{\widehat{\mathcal{U}}}=\widehat{U}\widehat{U}^{\top},\quad\quad\mathcal{Q}_{\mathcal{U}}=UU^{\top}.$

With those definitions in mind, we now begin deriving our error bound. From
Equation (13), we know that
$V^{(k)}\cong_{O}\begin{pmatrix}M&R^{(k)}\end{pmatrix}$. By defintion,

$\displaystyle\mathcal{Q}_{\widehat{\mathcal{U}}}$
$\displaystyle=\widehat{U}\widehat{U}^{\top}$
$\displaystyle=\begin{pmatrix}\widehat{V}^{(1)}&\widehat{V}^{(2)}\end{pmatrix}\begin{pmatrix}(\widehat{V}^{(1)})^{\top}\\\
(\widehat{V}^{(2)})^{\top}\end{pmatrix}$
$\displaystyle=\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}+\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}.$

Similarly,
$\mathcal{Q}_{\mathcal{U}}=\mathcal{P}_{\mathcal{V}^{(1)}}+\mathcal{P}_{\mathcal{V}^{(2)}}$.
Finally, we set

$R=\begin{pmatrix}R^{(1)}&R^{(2)}\end{pmatrix}$

and we write the eigendecomposition of $RR^{\top}$ as

$RR^{\top}=LDL^{\top}.$

As in the proof of Proposition 2.3, we have that the eigenvalue decomposition
of $\mathcal{Q}_{\mathcal{U}}$ is given by

$\displaystyle\mathcal{Q}_{\mathcal{U}}$
$\displaystyle=\mathcal{P}_{\mathcal{V}^{(1)}}+\mathcal{P}_{\mathcal{V}^{(2)}}$
$\displaystyle=2MM^{\top}+LDL^{\top}$
$\displaystyle=\begin{pmatrix}M&L\end{pmatrix}\begin{pmatrix}2I&0\\\ 0&D\\\
\end{pmatrix}\begin{pmatrix}M&L\end{pmatrix}^{\top}$

so that

$\displaystyle\lambda_{1}(D)$
$\displaystyle=1+\lambda_{1}\bigg{[}\begin{pmatrix}0&(R^{(1)})^{\top}R^{(2)}\\\
(R^{(2)})^{\top}R^{(1)}&0\\\ \end{pmatrix}\bigg{]}$

and

$\displaystyle
1+\lambda_{1}\bigg{[}\begin{pmatrix}0&(R^{(1)})^{\top}R^{(2)}\\\
(R^{(2)})^{\top}R^{(1)}&0\\\ \end{pmatrix}\bigg{]}$ $\displaystyle=2-\delta.$

Since $\delta>0$, we see that $\lambda_{1}(D)<2$ and thus
$\operatorname{EIG}(\mathcal{Q}_{\mathcal{U}},r_{M})=M$. In addition,

$\displaystyle\lambda_{r_{M}}(\mathcal{Q}_{\mathcal{U}})-\lambda_{r_{M}+1}(\mathcal{Q}_{\mathcal{U}})$
$\displaystyle=2-\lambda_{1}(D)$ $\displaystyle=2-(2-\delta)$
$\displaystyle=\delta.$

Therefore

$\frac{1}{\lambda_{r_{M}}(\mathcal{Q}_{\mathcal{U}})-\lambda_{r_{M}+1}(\mathcal{Q}_{\mathcal{U}})}\leq\frac{1}{\delta}.$

From here, Lemma $2.6$ in Chen et al. (2021) and Theorem $2$ in Yu et al.
(2015) imply that that for some $b_{0}>0$,

$\displaystyle\mathbb{E}[\min\limits_{Q\in\mathbb{O}_{r_{M}}}\|\widehat{M}-MQ\|_{F}]$
$\displaystyle\leq\frac{b_{0}\sqrt{r_{M}}\mathbb{E}\|\mathcal{Q}_{\widehat{\mathcal{U}}}-\mathcal{Q}_{\mathcal{U}}\|_{F}}{\lambda_{r_{M}}(\mathcal{Q}_{\mathcal{U}})-\lambda_{r_{M}+1}(\mathcal{Q}_{\mathcal{U}})}$
$\displaystyle\leq\frac{b_{0}\sqrt{r_{M}}}{\delta}\mathbb{E}\|\mathcal{Q}_{\widehat{\mathcal{U}}}-\mathcal{Q}_{\mathcal{U}}\|_{F}.$

To bound
$\|\mathcal{Q}_{\widehat{\mathcal{U}}}-\mathcal{Q}_{\mathcal{U}}\|_{F}$, we
note that

$\displaystyle\|\mathcal{Q}_{\widehat{\mathcal{U}}}-\mathcal{Q}_{\mathcal{U}}\|_{F}$
$\displaystyle=\bigg{\|}\sum_{i=1}^{2}\mathcal{P}_{\widehat{\mathcal{V}}^{(i)}}-\mathcal{P}_{\mathcal{V}^{(i)}}\bigg{\|}_{F}$
$\displaystyle\leq\sum_{i=1}^{2}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(i)}}-\mathcal{P}_{\mathcal{V}^{(i)}}\|_{F}.$

We make some useful observations about the singular vectors of
$\mathbb{E}(XX^{\top})$ and $W$. First, by definition,
$\operatorname{EIG}(P,r_{M}+r_{1})\cong_{O}V^{(1)}$. Second, from Equation
(13), we know that

$\displaystyle\mathbb{E}(XX^{\top})$
$\displaystyle=\mathbb{E}((X-W)(X-W)^{\top})+WW^{\top}$
$\displaystyle=p\tau^{2}I_{n}+WW^{\top}.$

Let $\Gamma^{(2)}(\Gamma^{(2)})^{\top}$ have eigendecomposition
$\Gamma^{(2)}(\Gamma^{(2)})^{\top}=Q\Lambda Q^{\top}$. Then we write

$\displaystyle WW^{\top}$
$\displaystyle=V^{(2)}\Gamma^{(2)}(V^{(2)}\Gamma^{(2)})^{\top}$
$\displaystyle=V^{(2)}Q\Lambda(V^{(2)}Q)^{\top}$
$\displaystyle=\begin{pmatrix}V^{(2)}&V^{(2)}_{\perp}\end{pmatrix}\begin{pmatrix}Q&0\\\
0&I_{n-r_{M}-r_{2}}\end{pmatrix}\begin{pmatrix}\Lambda&0\\\
0&0\end{pmatrix}\bigg{[}\begin{pmatrix}V^{(2)}&V^{(2)}_{\perp}\end{pmatrix}\begin{pmatrix}Q&0\\\
0&I_{n-r_{M}-r_{2}}\end{pmatrix}\bigg{]}^{\top}.$

Since
$\operatorname{EIG}(\mathbb{E}(XX^{\top}),r_{m}+r_{2})=\operatorname{EIG}(WW^{\top},r_{m}+r_{2})$
with corresponding eigenvalues
$\lambda_{i}(\mathbb{E}(XX^{\top}))=p\tau^{2}+\lambda_{i}(WW^{\top})$, we have
that

$\displaystyle\operatorname{EIG}(\mathbb{E}(XX^{\top}),r_{M}+r_{2})$
$\displaystyle=\operatorname{EIG}(WW^{\top},r_{m}+r_{2})$
$\displaystyle=\operatorname{EIG}\bigg{(}\begin{pmatrix}V^{(2)}&V^{(2)}_{\perp}\end{pmatrix}\begin{pmatrix}Q&0\\\
0&I_{n-r_{M}-r_{i}}\end{pmatrix}\bigg{[}\begin{pmatrix}V^{(2)}&V^{(2)}_{\perp}\end{pmatrix}\begin{pmatrix}Q&0\\\
0&I_{n-r_{M}-r_{2}}\end{pmatrix}\bigg{]}^{\top},r_{M}+r_{2}\bigg{)}$
$\displaystyle=\operatorname{EIG}(V^{(2)}QQ^{\top}(V^{(2)})^{\top}+V^{(2)}_{\perp}(V^{(2)}_{\perp})^{\top},r_{M}+r_{2})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\mathcal{V}^{(2)}}+V^{(2)}_{\perp}(V^{(2)}_{\perp})^{\top},r_{M}+r_{2})$
$\displaystyle=\operatorname{EIG}\bigg{(}\begin{pmatrix}V^{(2)}&V^{(2)}_{\perp}\end{pmatrix}\begin{pmatrix}V^{(2)}&V^{(2)}_{\perp}\end{pmatrix}^{\top},r_{M}+r_{2}\bigg{)}$
$\displaystyle\cong_{O}V^{(2)}.$

Since $\operatorname{EIG}(P,r_{M}+r_{1})=V^{(1)}$, we apply the Davis-Kahan
theorem to obtain that for some $b_{1}>0$

$\displaystyle\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|_{F}$
$\displaystyle\leq\frac{b_{1}\sqrt{r_{M}+r_{1}}\|A-P\|}{\lambda_{r_{M}+r_{1}}(P)}.$

Since $\|A-P\|$ is subexponential and for each $q\geq 1$,

$\mathbb{E}(\|A-P\|^{q})^{1/q}\leq\bigg{\|}\|A-P\|\bigg{\|}_{\psi_{1}}q.$

By Assumption 2, for each $q\geq 1$,
$\bigg{\|}\|A-P\|\bigg{\|}_{\psi_{1}}\leq\kappa$. Therefore

$\displaystyle\mathbb{E}\bigg{(}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|_{F}^{q}\bigg{)}^{1/q}$
$\displaystyle\leq\mathbb{E}\bigg{[}\frac{b_{1}^{q}(r_{M}+r_{1})^{q/2}\|A-P\|^{q}}{(\lambda_{r_{M}+r_{1}}(P))^{q}}\bigg{]}^{1/q}$
$\displaystyle=\frac{b_{1}\sqrt{r_{M}+r_{1}}\mathbb{E}(\|A-P\|^{q})^{1/q}}{\lambda_{r_{M}+r_{1}}(P)}$
$\displaystyle\leq\frac{b_{1}\sqrt{r_{M}+r_{1}}\bigg{\|}\|A-P\|\bigg{\|}_{\psi_{1}}q}{\lambda_{r_{M}+r_{1}}(P)}$
$\displaystyle\leq\frac{b_{1}\sqrt{r_{M}+r_{1}}\kappa
q}{\lambda_{r_{M}+r_{1}}(P)}.$

This implies that

$\displaystyle\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|_{F}$
$\displaystyle\leq\bigg{\|}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|_{F}\bigg{\|}_{\psi_{1}}$
$\displaystyle=\sup_{q\geq
1}\bigg{[}\frac{\mathbb{E}(\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|_{F}^{q})^{1/q}}{q}\bigg{]}$
$\displaystyle\leq\frac{b_{1}\sqrt{r_{M}+r_{1}}\kappa}{\lambda_{r_{M}+r_{1}}(P)}$
$\displaystyle=O\bigg{(}\epsilon^{(1)}\bigg{)}.$

Jensen’s inequality implies that
$\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|_{F}\leq\bigg{(}\mathbb{E}[\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|_{F}^{2}]\bigg{)}^{1/2}$.
Under Assumption 2, Theorem 3 in Cai and Zhang (2018) implies that there
exists $c_{0}>0$ such that

$\displaystyle\mathbb{E}[\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|_{F}^{2}]$
$\displaystyle\leq\frac{c_{0}n(\sigma_{r_{M}+r_{2}}^{2}(\tau^{-1}W)+p)(r_{M}+r_{2})}{\sigma_{r_{M}+r_{2}}^{4}(\tau^{-1}W)}\wedge(r_{M}+r_{2})$
$\displaystyle=\frac{c_{0}\tau^{2}n(\sigma_{r_{M}+r_{2}}^{2}(W)+\tau^{2}p)(r_{M}+r_{2})}{\sigma_{r_{M}+r_{2}}^{4}(W)}\wedge(r_{M}+r_{2}).$

Therefore,

$\displaystyle\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|_{F}$
$\displaystyle\leq\frac{\tau\sqrt{c_{0}n(\sigma_{r_{M}+r_{2}}^{2}(W)+\tau^{2}p)(r_{M}+r_{2})}}{\sigma_{r_{M}+r_{2}}^{2}(W)}\wedge\sqrt{r_{M}+r_{2}}$
$\displaystyle=O\bigg{(}\epsilon^{(2)}\bigg{)}.$

We therefore obtain

$\displaystyle\mathbb{E}\|\mathcal{Q}_{\widehat{\mathcal{U}}}-\mathcal{Q}_{\mathcal{U}}\|_{F}$
$\displaystyle=\mathbb{E}\bigg{\|}\sum_{i=1}^{2}\mathcal{P}_{\widehat{\mathcal{V}}^{(i)}}-\mathcal{P}_{\mathcal{V}^{(i)}}\bigg{\|}_{F}$
$\displaystyle\leq\sum_{i=1}^{2}\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(i)}}-\mathcal{P}_{\mathcal{V}^{(i)}}\|_{F}$
$\displaystyle=O\bigg{(}\epsilon^{(1)}+\epsilon^{(2)}\bigg{)}.$

Summarizing, we have that

$\displaystyle\mathbb{E}[d(\widehat{M},M)]$
$\displaystyle=\mathbb{E}[\min\limits_{Q\in\mathbb{O}_{r_{M}}}\|\widehat{M}-MQ\|_{F}]$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}}}{\delta}\mathbb{E}\|\mathcal{Q}_{\widehat{\mathcal{U}}}-\mathcal{Q}_{\mathcal{U}}\|_{F}\bigg{)}$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}}}{\delta}[\epsilon^{(1)}+\epsilon^{(2)}]\bigg{)}.$

  * •

Rate for Network Individual Components:
We note that

$\displaystyle\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}$
$\displaystyle=\bigg{(}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\bigg{)}$
$\displaystyle\hskip
28.45274pt+\bigg{(}\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}\bigg{)}$
$\displaystyle=(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}})\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}$
$\displaystyle\hskip
28.45274pt+\mathcal{P}_{\mathcal{M}^{\perp}}(\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}).$

The expression in parentheses in the second term can be further expanded as

$\displaystyle\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}$
$\displaystyle=\bigg{(}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}\bigg{)}+\bigg{(}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}-\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}\bigg{)}$
$\displaystyle=\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}})+(\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}})\mathcal{P}_{\mathcal{M}^{\perp}}.$

Since the spectral norm of a projection matrix is equal to one, we have that

$\displaystyle\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|$
$\displaystyle\leq\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\|\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\|\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\|$
$\displaystyle\hskip
28.45274pt+\|\mathcal{P}_{\mathcal{M}^{\perp}}\|\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\|\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\|$
$\displaystyle\hskip
28.45274pt+\|\mathcal{P}_{\mathcal{M}^{\perp}}\|\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|\|\mathcal{P}_{\mathcal{M}^{\perp}}\|$
$\displaystyle\leq
2\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\|+\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|$
$\displaystyle=2\|\mathcal{P}_{\widehat{\mathcal{M}}}-\mathcal{P}_{\mathcal{M}}\|+\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|.$

Since

$\displaystyle\widehat{R}^{(1)}$
$\displaystyle=\operatorname{SV}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\widehat{V}^{(1)},r_{1})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\widehat{V}^{(1)}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\widehat{V}^{(1)})^{\top},r_{1})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}},r_{1}),$
$\displaystyle R^{(1)}$
$\displaystyle=\operatorname{SV}(\mathcal{P}_{\mathcal{M}^{\perp}}V^{(1)},r_{1})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\mathcal{M}^{\perp}}V^{(1)}(\mathcal{P}_{\mathcal{M}^{\perp}}V^{(1)})^{\top},r_{1})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}},r_{1})$

and

$\displaystyle\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}$
$\displaystyle=\mathcal{P}_{\mathcal{M}^{\perp}}(MM^{\top}+R^{(1)}(R^{(1)})^{\top})\mathcal{P}_{\mathcal{M}^{\perp}}$
$\displaystyle=R^{(1)}(R^{(1)})^{\top}$
$\displaystyle=\mathcal{P}_{\mathcal{R}^{(1)}},$

the Davis-Kahan theorem and the fact that orthogonal matrices have unit
nontrivial singular values imply that

$\displaystyle\|\widehat{R^{(1)}}(\widehat{R^{(1)}})^{\top}-R^{(1)}(R^{(1)})^{\top}\|_{F}$
$\displaystyle\leq
2\sqrt{r_{1}}\frac{\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|}{\lambda_{r_{1}}(\mathcal{P}_{R^{(1)}})}$
$\displaystyle=2\sqrt{r_{1}}\frac{\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|}{\sigma_{r_{1}}(R^{(1)})^{2}}$
$\displaystyle=2\sqrt{r_{1}}\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|$

and that

$\displaystyle\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|$
$\displaystyle\leq\frac{\sqrt{2}\mathbb{E}\|A-P\|}{\lambda_{r_{M}+r_{1}}(P)}$
$\displaystyle\leq\frac{\sqrt{2}\bigg{\|}\|A-P\|\bigg{\|}_{\psi_{1}}}{\lambda_{r_{M}+r_{1}}(P)}$
$\displaystyle\leq\frac{\sqrt{2}\sqrt{r_{M}+r_{1}}\kappa}{\lambda_{r_{M}+r_{1}}(P)}$
$\displaystyle=\sqrt{2}\epsilon^{(1)}.$

Hence

$\displaystyle\mathbb{E}[d(\widehat{R}^{(1)},R^{(1)})]$
$\displaystyle=\mathbb{E}\bigg{[}\min_{Q\in\mathbb{O}_{r_{1}}}\|R^{(1)}-\widehat{R}^{(1)}Q\|_{F}\bigg{]}$
$\displaystyle=O\bigg{(}\sqrt{r_{1}}\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(1)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|\bigg{)}$
$\displaystyle=O\bigg{(}\sqrt{r_{1}}\bigg{[}\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{M}}}-\mathcal{P}_{\mathcal{M}}\|+\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(1)}}-\mathcal{P}_{\mathcal{V}^{(1)}}\|\bigg{]}\bigg{)}$
$\displaystyle=O\bigg{(}\sqrt{r_{M}r_{1}}\bigg{[}\frac{\epsilon^{(1)}+\epsilon^{(2)}}{\delta}+\epsilon^{(1)}\bigg{]}\bigg{)}$
$\displaystyle=O\bigg{(}\sqrt{r_{M}r_{1}}\bigg{[}\frac{\epsilon^{(1)}+\epsilon^{(2)}}{\delta}\bigg{]}\bigg{)}$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}r_{1}}}{\delta}\bigg{[}\epsilon^{(1)}+\epsilon^{(2)}\bigg{]}\bigg{)}.$

  * •

Rate for Covariate Individual Components:
Similarly to the network individual components, we note that

$\displaystyle\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}$
$\displaystyle=\bigg{(}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\bigg{)}$
$\displaystyle\hskip
28.45274pt+\bigg{(}\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}\bigg{)}$
$\displaystyle=(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}})\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}$
$\displaystyle\hskip
28.45274pt+\mathcal{P}_{\mathcal{M}^{\perp}}(\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}).$

The expression in parentheses in the second term can be further expanded as

$\displaystyle\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}$
$\displaystyle=\bigg{(}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}\bigg{)}+\bigg{(}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}-\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}\bigg{)}$
$\displaystyle=\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}})+(\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}})\mathcal{P}_{\mathcal{M}^{\perp}}.$

Since the spectral norm of a projection matrix is $1$, we have that

$\displaystyle\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|$
$\displaystyle\leq\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\|\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\|\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\|$
$\displaystyle\hskip
28.45274pt+\|\mathcal{P}_{\mathcal{M}^{\perp}}\|\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\|\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\|$
$\displaystyle\hskip
28.45274pt+\|\mathcal{P}_{\mathcal{M}^{\perp}}\|\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|\|\mathcal{P}_{\mathcal{M}^{\perp}}\|$
$\displaystyle\leq
2\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\|+\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|$
$\displaystyle=2\|\mathcal{P}_{\widehat{\mathcal{M}}}-\mathcal{P}_{\mathcal{M}}\|+\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|.$

Since

$\displaystyle\widehat{R}^{(2)}$
$\displaystyle=\operatorname{SV}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\widehat{V}^{(2)},r_{2})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\widehat{V}^{(2)}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\widehat{V}^{(2)})^{\top},r_{2})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}},r_{2}),$
$\displaystyle R^{(2)}$
$\displaystyle=\operatorname{SV}(\mathcal{P}_{\mathcal{M}^{\perp}}V^{(2)},r_{2})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\mathcal{M}^{\perp}}V^{(2)}(\mathcal{P}_{\mathcal{M}^{\perp}}V^{(2)})^{\top},r_{2})$
$\displaystyle=\operatorname{EIG}(\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}},r_{2})$

and

$\displaystyle\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}$
$\displaystyle=\mathcal{P}_{\mathcal{M}^{\perp}}(MM^{\top}+R^{(2)}(R^{(2)})^{\top})\mathcal{P}_{\mathcal{M}^{\perp}}$
$\displaystyle=R^{(2)}(R^{(2)})^{\top}$
$\displaystyle=\mathcal{P}_{\mathcal{R}^{(2)}},$

the Davis-Kahan theorem and the fact that orthogonal matrices have unit
nontrivial singular values imply that

$\displaystyle\|\widehat{R^{(2)}}(\widehat{R^{(2)}})^{\top}-R^{(2)}(R^{(2)})^{\top}\|_{F}$
$\displaystyle\leq
2\sqrt{r_{2}}\frac{\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|}{\lambda_{r_{2}}(\mathcal{P}_{\mathcal{R}^{(2)}})}$
$\displaystyle=2\sqrt{r_{2}}\frac{\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|}{\sigma_{r_{2}}(R^{(2)})^{2}}$
$\displaystyle=2\sqrt{r_{2}}\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|.$

Theorem 3 in Cai and Zhang (2018) then implies that there exists $c_{0}>0$
such that

$\displaystyle\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|$
$\displaystyle\leq(\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|^{2})^{1/2}$
$\displaystyle\leq\bigg{[}\frac{c_{0}\tau^{2}n(\sigma_{r_{M}+r_{2}}(W)^{2}+p)}{\sigma_{r_{M}+r_{2}}(W)^{4}}\wedge
1\bigg{]}^{1/2}$
$\displaystyle\leq\sqrt{r_{M}+r_{2}}\bigg{[}\frac{c_{0}\tau^{2}n(\sigma_{r_{M}+r_{2}}(W)^{2}+p)}{\sigma_{r_{M}+r_{2}}(W)^{4}}\wedge
1\bigg{]}^{1/2}$ $\displaystyle=\sqrt{c_{0}}\epsilon^{(2)}.$

Hence

$\displaystyle\mathbb{E}[d(\widehat{R}^{(1)},R^{(2)})]$
$\displaystyle=\mathbb{E}\bigg{[}\min_{Q\in\mathbb{O}_{r_{2}}}\|R^{(2)}-\widehat{R}^{(2)}Q\|_{F}\bigg{]}$
$\displaystyle=O\bigg{(}\sqrt{r_{2}}\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}\mathcal{P}_{\widehat{\mathcal{M}}^{\perp}}-\mathcal{P}_{\mathcal{M}^{\perp}}\mathcal{P}_{\mathcal{V}^{(2)}}\mathcal{P}_{\mathcal{M}^{\perp}}\|\bigg{)}$
$\displaystyle=O\bigg{(}\sqrt{r_{2}}\bigg{[}\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{M}}}-\mathcal{P}_{\mathcal{M}}\|+\mathbb{E}\|\mathcal{P}_{\widehat{\mathcal{V}}^{(2)}}-\mathcal{P}_{\mathcal{V}^{(2)}}\|\bigg{]}\bigg{)}$
$\displaystyle=O\bigg{(}\sqrt{r_{M}r_{2}}\bigg{[}\frac{\epsilon^{(1)}+\epsilon^{(2)}}{\delta}+\epsilon^{(2)}\bigg{]}\bigg{)}$
$\displaystyle=O\bigg{(}\sqrt{r_{M}r_{2}}\bigg{[}\frac{\epsilon^{(1)}+\epsilon^{(2)}}{\delta}\bigg{]}\bigg{)}$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}r_{2}}}{\delta}\bigg{[}\epsilon^{(1)}+\epsilon^{(2)}\bigg{]}\bigg{)}.$

∎

###  A.5 Proof of Corollary 3.2

###### Proof.

[Proof of Corollary 3.2] From Vershynin (2010), we know that
$\bigg{\|}\|A-P\|\bigg{\|}_{\psi_{1}}$ differs from $\kappa$ by a
multiplicative constant. Since $\mu(P)=\Omega(\rho_{n}n)$ and
$\rho_{n}=\omega(n^{-1}\log n)$, we have that $\mu(P)=\omega(\log(n))$. The
proof of Lemma 13 in Arroyo et al. (2021) implies that

$\displaystyle\bigg{\|}\|P-A\|\bigg{\|}_{\psi_{1}}$
$\displaystyle=O\bigg{(}\sqrt{\mu(P)}\bigg{)}$
$\displaystyle=O\bigg{(}\sqrt{\rho_{n}n}\bigg{)}.$

Therefore

$\displaystyle\epsilon^{(1)}$
$\displaystyle=\frac{\sqrt{r_{M}+r_{1}}\kappa}{\lambda_{r_{M}+r_{1}}(P)}$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}+r_{1}}\bigg{\|}\|A-P\|\bigg{\|}_{\psi_{1}}}{\lambda_{r_{M}+r_{1}}(P)}\bigg{)}$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}+r_{1}}\sqrt{\rho_{n}n}}{\lambda_{r_{M}+r_{1}}(P)}\bigg{)}$
$\displaystyle=O\bigg{(}\frac{1}{\sqrt{\rho_{n}n}}\bigg{)}$

and

$\displaystyle\epsilon^{(2)}$
$\displaystyle=\frac{\tau\sqrt{n(r_{M}+r_{2})(\sigma_{r_{M}+r_{2}}^{2}(W)+p)}}{\sigma_{r_{M}+r_{2}}^{2}(W)}\wedge\sqrt{r_{M}+r_{2}}$
$\displaystyle=\bigg{[}\tau\bigg{(}\frac{n(r_{M}+r_{2})}{\sigma_{r_{M}+r_{2}}^{2}(W)}+\frac{n(r_{M}+r_{2})p}{\sigma_{r_{M}+r_{2}}^{4}(W)}\bigg{)}^{1/2}\bigg{]}\wedge\sqrt{r_{M}+r_{2}}$
$\displaystyle=O\bigg{(}\bigg{[}\tau\bigg{(}\frac{1}{p}+\frac{1}{np}\bigg{)}^{1/2}\bigg{]}\wedge
1\bigg{)}$ $\displaystyle=O\bigg{(}\frac{\tau}{\sqrt{p}}\bigg{)}.$

Hence
$\epsilon^{(1)}+\epsilon^{(2)}=O\bigg{(}\frac{1}{\sqrt{\rho_{n}n}}+\frac{\tau}{\sqrt{p}}\bigg{)}$,
so that

$\displaystyle\mathbb{E}[d(\widehat{M},M)]$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}}}{\delta}[\epsilon^{(1)}+\epsilon^{(2)}]\bigg{)}$
$\displaystyle=O\bigg{(}\frac{1}{\sqrt{\rho_{n}n}}+\frac{\tau}{\sqrt{p}}\bigg{)}$

and

$\displaystyle\mathbb{E}[d(\widehat{R}^{(k)},R^{(k)})]$
$\displaystyle=O\bigg{(}\frac{\sqrt{r_{M}r_{k}}}{\delta}[\epsilon^{(1)}+\epsilon^{(2)}]\bigg{)}$
$\displaystyle=O\bigg{(}\frac{1}{\sqrt{\rho_{n}n}}+\frac{\tau}{\sqrt{p}}\bigg{)}.$

Since $\rho_{n}n=\omega(\log n)$, we have that $(\rho_{n}n)^{-1/2}=o((\log
n)^{-1/2})$ and consistency follows. ∎

## References

  * Arroyo et al. (2021) Arroyo, J., A. Athreya, J. Cape, G. Chen, C. E. Priebe, and J. T. Vogelstein (2021).  Inference for Multiple Heterogeneous Networks with a Common Invariant Subspace.  Journal of Machine Learning Research 22(142), 1–49.
  * Athreya et al. (2017) Athreya, A., D. E. Fishkind, M. Tang, C. E. Priebe, Y. Park, J. T. Vogelstein, K. Levin, V. Lyzinski, and Y. Qin (2017).  Statistical Inference on Random Dot Product Graphs: a Survey.  The Journal of Machine Learning Research 18(1), 8393–8484.
  * Binkiewicz et al. (2017) Binkiewicz, N., J. T. Vogelstein, and K. Rohe (2017).  Covariate-assisted spectral clustering.  Biometrika 104(2), 361–377.
  * Cai and Zhang (2018) Cai, T. T. and A. Zhang (2018).  Rate-optimal perturbation bounds for singular subspaces with applications to high-dimensional statistics.  The Annals of Statistics 46(1), 60 – 89.
  * Chatterjee (2015) Chatterjee, S. (2015).  Matrix estimation by universal singular value thresholding.  The Annals of Statistics 43(1), 177–214.
  * Chen et al. (2022) Chen, H., B. Caffo, G. Stein-O’Brien, J. Liu, B. Langmead, C. Colantuoni, and L. Xiao (2022).  Two-stage linked component analysis for joint decomposition of multiple biologically related data sets.  Biostatistics 23(4), 1200–1217.
  * Chen et al. (2021) Chen, Y., Y. Chi, J. Fan, and C. Ma (2021).  Spectral methods for data science: A statistical perspective.  Foundations and Trends® in Machine Learning 14(5), 566–806.
  * Chunaev (2020) Chunaev, P. (2020).  Community detection in node-attributed social networks: a survey.  Computer Science Review 37, 100286.
  * Chung and Lu (2006) Chung, F. R. and L. Lu (2006).  Complex graphs and networks.  Number 107. American Mathematical Soc.
  * Crainiceanu et al. (2011) Crainiceanu, C. M., B. S. Caffo, S. Luo, V. M. Zipunnikov, and N. M. Punjabi (2011).  Population value decomposition, a framework for the analysis of image populations.  Journal of the American Statistical Association 106(495), 775–790.
  * De Domenico et al. (2015) De Domenico, M., V. Nicosia, A. Arenas, and V. Latora (2015).  Structural reducibility of multilayer networks.  Nature Communications 6(1), 6864.
  * Deshpande et al. (2018) Deshpande, Y., S. Sen, A. Montanari, and E. Mossel (2018).  Contextual Stochastic Block Models.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Eds.), Advances in Neural Information Processing Systems, Volume 31. Curran Associates, Inc.
  * Eynard et al. (2015) Eynard, D., A. Kovnatsky, M. M. Bronstein, K. Glashoff, and A. M. Bronstein (2015).  Multimodal Manifold Analysis by Simultaneous Diagonalization of Laplacians.  IEEE Transactions on Pattern Analysis and Machine Intelligence 37(12), 2505–2517.
  * Feng et al. (2018) Feng, Q., M. Jiang, J. Hannig, and J. Marron (2018).  Angle-based joint and individual variation explained.  Journal of Multivariate Analysis 166, 241–265.
  * Fosdick and Hoff (2015) Fosdick, B. K. and P. D. Hoff (2015).  Testing and modeling dependencies between a network and nodal attributes.  Journal of the American Statistical Association 110(511), 1047–1056.
  * Gavish and Donoho (2017) Gavish, M. and D. L. Donoho (2017).  Optimal shrinkage of singular values.  IEEE Transactions on Information Theory 63(4), 2137–2152.
  * Gaynanova and Li (2019) Gaynanova, I. and G. Li (2019).  Structural learning and integrative decomposition of multi-view data.  Biometrics 75(4), 1121–1132.
  * Han et al. (2023) Han, X., Q. Yang, and Y. Fan (2023).  Universal rank inference via residual subsampling with application to large networks.  The Annals of Statistics 51(3), 1109–1133.
  * Harold (1936) Harold, H. (1936).  Relations between two sets of variables.  Biometrika 28(3), 321–377.
  * Hoff (2005) Hoff, P. D. (2005).  Bilinear mixed-effects models for dyadic data.  Journal of the American Statistical Association 100(469), 286–295.
  * Hoff et al. (2002) Hoff, P. D., A. E. Raftery, and M. S. Handcock (2002).  Latent space approaches to social network analysis.  Journal of the American Statistical Association 97(460), 1090–1098.
  * Holland et al. (1983) Holland, P. W., K. B. Laskey, and S. Leinhardt (1983).  Stochastic blockmodels: First steps.  Social Networks 5(2), 109–137.
  * Huang et al. (2023) Huang, S., J. Sun, and Y. Feng (2023).  PCABM: Pairwise Covariates-Adjusted Block Model for Community Detection.  Journal of the American Statistical Association, 1–13.
  * Le and Li (2022) Le, C. M. and T. Li (2022).  Linear regression and its inference on noisy network-linked data.  Journal of the Royal Statistical Society Series B: Statistical Methodology 84(5), 1851–1885.
  * Lee et al. (2019) Lee, Y., C. Shen, C. E. Priebe, and J. T. Vogelstein (2019).  Network dependence testing via diffusion maps and distance-based correlations.  Biometrika 106(4), 857–873.
  * Li et al. (2024) Li, J., G. Xu, and J. Zhu (2024).  High-dimensional factor analysis for network-linked data.  arXiv preprint arXiv:2403.17121.
  * Li et al. (2019) Li, T., E. Levina, and J. Zhu (2019).  Prediction models for network-linked data.  The Annals of Applied Statistics 13(1), 132 – 164.
  * Liu (2019) Liu, B. (2019).  Statistical learning for networks with node features.  Ph. D. thesis.
  * Lock et al. (2013) Lock, E. F., K. A. Hoadley, J. S. Marron, and A. B. Nobel (2013).  Joint and individual variation explained (JIVE) for integrated analysis of multiple data types.  The Annals of Applied Statistics 7(1), 523.
  * MacDonald et al. (2022) MacDonald, P. W., E. Levina, and J. Zhu (2022).  Latent space models for multiplex networks with shared structure.  Biometrika 109(3), 683–706.
  * Murden et al. (2022) Murden, R. J., Z. Zhang, Y. Guo, and B. B. Risk (2022).  Interpretive JIVE: Connections with CCA and an application to brain connectivity.  Frontiers in Neuroscience 16, 969510.
  * Newman and Clauset (2016) Newman, M. E. and A. Clauset (2016).  Structure and inference in annotated networks.  Nature Communications 7(1), 11863.
  * Palzer et al. (2022) Palzer, E. F., C. H. Wendt, R. P. Bowler, C. P. Hersh, S. E. Safo, and E. F. Lock (2022).  sjive: Supervised joint and individual variation explained.  Computational Statistics & Data Analysis 175, 107547.
  * Paul and Chen (2020) Paul, S. and Y. Chen (2020).  Spectral and matrix factorization methods for consistent community detection in multi-layer networks.
  * Piburn (2020) Piburn, J. (2020).  wbstats: Programmatic Access to the World Bank API.  Oak Ridge, Tennessee: Oak Ridge National Laboratory.
  * Prothero et al. (2024) Prothero, J., M. Jiang, J. Hannig, Q. Tran-Dinh, A. Ackerman, and J. Marron (2024).  Data integration via analysis of subspaces (DIVAS).  TEST, 1–42.
  * Rohe et al. (2011) Rohe, K., S. Chatterjee, and B. Yu (2011).  Spectral clustering and the high-dimensional stochastic blockmodel.  The Annals of Statistics 39(4), 1878––1915.
  * Rubin-Delanchy et al. (2022) Rubin-Delanchy, P., J. Cape, M. Tang, and C. E. Priebe (2022, 06).  A Statistical Interpretation of Spectral Embedding: The Generalised Random Dot Product Graph.  Journal of the Royal Statistical Society Series B: Statistical Methodology 84(4), 1446–1473.
  * Shen et al. (2009) Shen, R., A. B. Olshen, and M. Ladanyi (2009).  Integrative clustering of multiple genomic data types using a joint latent variable model with application to breast and lung cancer subtype analysis.  Bioinformatics 25(22), 2906–2912.
  * Sussman et al. (2012) Sussman, D. L., M. Tang, D. E. Fishkind, and C. E. Priebe (2012).  A consistent adjacency spectral embedding for stochastic blockmodel graphs.  Journal of the American Statistical Association 107(499), 1119–1128.
  * Tao (2023) Tao, T. (2023).  Topics in random matrix theory, Volume 132.  American Mathematical Society.
  * Udell and Townsend (2019) Udell, M. and A. Townsend (2019).  Why are big data matrices approximately low rank?  SIAM Journal on Mathematics of Data Science 1(1), 144–160.
  * Vershynin (2010) Vershynin, R. (2010).  Introduction to the non-asymptotic analysis of random matrices.  arXiv preprint arXiv:1011.3027.
  * Wang et al. (2019) Wang, L., Z. Zhang, and D. Dunson (2019).  Common and individual structure of brain networks.  The Annals of Applied Statistics 13(1), 85–112.
  * Wang et al. (2019) Wang, S., J. Arroyo, J. T. Vogelstein, and C. E. Priebe (2019).  Joint embedding of graphs.  IEEE Transactions on Pattern Analysis and Machine Intelligence 43(4), 1324–1336.
  * Witten and Tibshirani (2009) Witten, D. M. and R. J. Tibshirani (2009).  Extensions of sparse canonical correlation analysis with applications to genomic data.  Statistical Applications in Genetics and Molecular Biology 8(1).
  * Wold et al. (1996) Wold, S., N. Kettaneh, and K. Tjessem (1996).  Hierarchical multiblock PLS and PC models for easier model interpretation and as an alternative to variable selection.  Journal of Chemometrics 10(5-6), 463–482.
  * Yan and Sarkar (2021) Yan, B. and P. Sarkar (2021).  Covariate Regularized Community Detection in Sparse Graphs.  Journal of the American Statistical Association 116(534), 734–745.
  * Young and Scheinerman (2007) Young, S. J. and E. R. Scheinerman (2007).  Random dot product graph models for social networks.  In Algorithms and Models for the Web-Graph: 5th International Workshop, WAW 2007, San Diego, CA, USA, December 11-12, 2007. Proceedings 5, pp. 138–149. Springer.
  * Yu et al. (2015) Yu, Y., T. Wang, and R. J. Samworth (2015).  A useful Variant of the Davis–Kahan theorem for Statisticians.  Biometrika 102(2), 315–323.
  * Yuan and Gaynanova (2022) Yuan, D. and I. Gaynanova (2022).  Double-matched matrix decomposition for multi-view data.  Journal of Computational and Graphical Statistics 31(4), 1114–1126.
  * Zhang et al. (2016) Zhang, Y., E. Levina, and J. Zhu (2016).  Community detection in networks with node features.  Electronic Journal of Statistics 10(2), 3153 – 3178.
  * Zhang et al. (2018) Zhang, Z., G. I. Allen, H. Zhu, and D. Dunson (2018).  Relationships between human brain structural connectomes and traits.  BioRxiv, 256933.
  * Zhao et al. (2022) Zhao, J., X. Liu, H. Wang, and C. Leng (2022).  Dimension reduction for covariates in network data.  Biometrika 109(1), 85–102.