plain

# On varimax asymptotics in network models and spectral methods for
dimensionality reduction

J. Cape<EMAIL_ADDRESS>Department of Statistics, University of
Wisconsin–Madison,
1300 University Avenue, Madison, Wisconsin 53706 U.S.A.

###### Abstract

Varimax factor rotations, while popular among practitioners in psychology and
statistics since being introduced by H. Kaiser, have historically been viewed
with skepticism and suspicion by some theoreticians and mathematical
statisticians. Now, work by K. Rohe and M. Zeng provides new, fundamental
insight: varimax rotations provably perform statistical estimation in certain
classes of latent variable models when paired with spectral-based matrix
truncations for dimensionality reduction. We build on this newfound
understanding of varimax rotations by developing further connections to
network analysis and spectral methods rooted in entrywise matrix perturbation
analysis. Concretely, this paper establishes the asymptotic multivariate
normality of vectors in varimax-transformed Euclidean point clouds that
represent low-dimensional node embeddings in certain latent space random graph
models. We address related concepts including network sparsity, data
denoising, and the role of matrix rank in latent variable parameterizations.
Collectively, these findings, at the confluence of classical and contemporary
multivariate analysis, reinforce methodology and inference procedures grounded
in matrix factorization-based techniques. Numerical examples illustrate our
findings and supplement our discussion.

###### keywords:

Data embedding; factor analysis; latent variable; network; random graph;
varimax rotation.

††journal: Biometrika

##  section Introduction

###  section.subsection Background

Factor analysis has a long, rich history as a tenet of multivariate data
analysis in statistics and psychology. Stemming from foundational work by F.
Galton, K. Pearson, and C. Spearman more than a century ago, factor analysis
developed under the auspices of numerous researchers including L. Thurstone,
H. Kaiser, H. Rubin, M.G. Kendall, and T.W. Anderson (Vincent, 1953;
Bartholomew, 1995). In the twenty-first century, conducting factor analysis
remains a mainstay of applied scientific domains, finding widespread use
throughout psychology (Fabrigar et al., 1999; Brown, 2015), economics
(Chamberlain & Rothschild, 1983; Bai & Ng, 2008), and biology (Alter et al.,
2000; Argelaguet et al., 2018).

Much has been written about factor analysis. According to some:

> Factor analysis is a branch of statistical theory concerned with the
> resolution of a set of descriptive variables in terms of a small number of
> categories or factors. (Holzinger & Harman, 1941, p. 3)

According to others:

> Factor analysis is not a purely statistical technique; there is always a
> certain amount of guesswork in it which no fair-minded person would attempt
> to deny; but factor analysis is chiefly a research tool, and there is of
> necessity a lot of guesswork in all research […] (Vincent, 1953, p. 117)

Ultimately, matrices lie at the heart of this apparent dichotomy. In their
role as linear transformations, different matrices are used in practice to
obtain differently rotated estimated factors or loadings, with the goal of
obtaining scientifically meaningful, interpretable coordinate directions
(Thurstone, 1935). Practitioners frequently select from a collection of
diverse optimization criteria, leading to differently transformed point cloud
configurations in low-dimensional Euclidean space via orthogonal or oblique
matrices, with the most common criteria being promax, quartimax, equimax, and
finally varimax (Kaiser, 1958), the focus of this paper.

Historically, researchers have sought to develop theory that coherently
justifies the aforementioned process of rotating estimates and the widespread
subsequent empirical validation of exploratory analyses. Such efforts have
largely centered around classical Gaussian data generating models and maximum
likelihood-based procedures (Anderson & Rubin, 1956; Anderson, 1962). In the
case of varimax factor rotation matrices, which are orthogonal matrices, the
process of rotating estimates neither alters the behavior of standard Gaussian
likelihoods nor does it affect overall least-squares reconstruction error.
Taken together, these observations have led some to believe that factor
analysis is plagued by a fundamental conundrum.

At long last, Rohe & Zeng (2023) sheds new light on the statistical properties
of varimax factor rotations. The authors place their results in the context of
certain so-called semiparametric factor models which produce data matrices in
high-dimensional ambient space that possess low-dimensional population-level
latent structure. Under certain identifiability conditions and assumptions
that are markedly non-Gaussian, it is shown that pairing varimax rotations
with the truncated singular value decomposition of an observed data matrix
allows one to accurately estimate underlying generative variables. Put
differently, Rohe & Zeng (2023) establishes that varimax rotations provably
perform statistical estimation in certain classes of latent variable models
when applied to spectral-based matrix truncations for dimensionality
reduction.

The novel theoretical insights provided in Rohe & Zeng (2023) have been
enabled by the contemporary emergence and proliferation of entrywise matrix
perturbation analysis (Abbe et al., 2020; Cape et al., 2019b, a; Fan et al.,
2018, 2021), partially summarized in the monograph Chen et al. (2021). At a
high level, this line of work provides fine-grained high-resolution analysis
of structured data matrices and their perturbations arising from statistical
models. For such settings, matrix perturbations pertaining to estimation
procedures can be precisely analyzed and significantly differ from
deterministic worst-case behavior reflected by the classical theorems of C.
Davis, W. Kahan, P.-Å. Wedin, and H. Weyl found in matrix analysis textbooks.

Random graph and network models, typified by the study of large-dimensional
adjacency matrices that represent pairwise interactions between vertices, or
nodes, via the presence or absence of edges, or links, are among the
statistical settings where modern matrix perturbation tools are advancing
methodology and inference capabilities. As addressed in the present paper,
entrywise and indeed row-wise perturbation analysis can be used to obtain
high-resolution uniform error rates and asymptotic fluctuations that hold at
the level of individual node embeddings; such results improve upon coarser,
average error rates available from classical matrix theory. Further still,
certain random graph models that possess latent structure can be viewed as
cousins of general factor models, thereby motivating concurrent interest in
varimax factor rotations for obtaining interpretable coordinate
representations in network settings.

###  section.subsection Contributions of this paper

This paper studies the large-sample theoretical properties of varimax-
transformed network embeddings in the context of latent variable stochastic
graph models. We establish asymptotic multivariate normality results for low-
dimensional node-specific embedding vectors which estimate generative latent
variables. Our results hold for models of undirected networks, directed
networks, and degree-corrected networks. Taken together, this paper builds
upon the uniform error analysis of varimax conducted in Rohe & Zeng (2023).
Our results contribute to the contemporary proliferation of theory, methods,
and applications for spectral graph clustering, and for multivariate analysis
more generally. We provide further evidence that varimax factor rotations,
when paired with spectral methods for dimensionality reduction, provide a
coherent spectral estimation strategy. Spectral methods for dimensionality
reduction offer more than merely exploratory data analysis capabilities: they
do, in fact, possess significance criteria and quantifiable estimation
accuracy by virtue of our being able to discern the behavior of their
asymptotic fluctuations. Moreover, given the simplicity, flexibility, and
computational scalability afforded by truncated matrix factorizations,
spectral methods serve as a useful alternative to likelihood-based methods in
practice.

This paper provides numerical examples but deliberately avoids real data
analysis with varimax, since hundreds of existing papers have already
demonstrated its empirical success and show post hoc improved interpretability
of rotated estimated loadings and factors. Rather, this work is one of only a
few existing papers to date that studies the theoretical properties of varimax
factor rotations in statistical latent variable models. Interested readers can
find further discussion of related twenty-first century real data applications
and methodological advances in, for example, Rohe & Zeng (2023); Chen & Rohe
(2020).

##  section Preliminaries

The generic factor model posits that observed data
$X=[x_{1},\dots,x_{n}]\in\mathbb{R}^{p\times n}$ take the form

$X=\mu 1_{n}^{{\mathrm{\scriptscriptstyle T}}}+LF^{{\mathrm{\scriptscriptstyle
T}}}+E,$ (equation)

where $\mu\in\mathbb{R}^{p}$ is an intercept vector, $1_{n}\in\mathbb{R}^{n}$
is a vector of all ones, $L=[l_{1},\dots,l_{p}]^{{\mathrm{\scriptscriptstyle
T}}}\in\mathbb{R}^{p\times d}$ is a matrix of unobserved factor loadings,
$F=[f_{1},\dots,f_{n}]^{{\mathrm{\scriptscriptstyle T}}}\in\mathbb{R}^{n\times
d}$ is a matrix of unobserved individual factors, and $E\in\mathbb{R}^{p\times
n}$ is a matrix of error terms having mean zero, typically presumed to be
uncorrelated with or independent of $L$ and $F$. The problem of interest is to
estimate $\mu$, $L$, and $F$, given knowledge of the true dimension $d\ll p,n$
or with an estimate thereof. Additional identifiability conditions are
commonly assumed to enable inference in Eq. equation, e.g., see Anderson &
Rubin (1956); Fan et al. (2021).

Classical matrix perturbation theory, based on Weyl’s inequality for singular
values (Weyl, 1912) and Wedin’s sine theta theorem for singular subspaces
(Wedin, 1972), establishes that the data matrix $X$, after centering, is well-
approximated by $LF^{{\mathrm{\scriptscriptstyle T}}}$ when the error term $E$
has a comparatively small effect, for example, when its operator norm is small
relative to that of $LF^{{\mathrm{\scriptscriptstyle T}}}$. Given that the
matrix $LF^{{\mathrm{\scriptscriptstyle T}}}$ is necessarily low rank since
$d\ll p,n$, a natural candidate estimator is to take a singular value
truncation of $X$. Here, the challenge remains how to construct estimates
$\widehat{L}$ and $\widehat{F}$ from $X$ so as to then quantify their
statistical properties at a more refined level of granularity than mere
operator norm or Frobenius norm bounds permit.

In the above display, Eq. equation is an example of a signal-plus-noise matrix
model, in which an observed matrix $X$ additively decomposes into a signal
term $\mu 1_{n}^{{\mathrm{\scriptscriptstyle
T}}}+LF^{{\mathrm{\scriptscriptstyle T}}}$ and a noise term $E$. Spectral
methods are particularly well suited for such models which naturally arise in
numerous statistical applications including network analysis, image denoising,
and covariance matrix estimation. In these settings, computing an appropriate
truncated singular value decomposition often has the effect of denoising data,
in an effort to approximate the underlying, unobserved signal. Our study of
varimax factor rotations is rooted in this context.

Under the semiparametric factor model of Rohe & Zeng (2023), though here using
slightly different notation, one observes a data matrix
$A\in\mathbb{R}^{n\times m}$ consisting of independent entries such that, for
latent factor matrices $Z\in\mathbb{R}^{n\times d}$, $Y\in\mathbb{R}^{m\times
d}$ and a latent full rank matrix of coefficients $B\in\mathbb{R}^{d\times
d}$, at the population level, in expectation,

$\operatorname{E}(A\mid Z,Y)=ZBY^{{\mathrm{\scriptscriptstyle T}}}.$
(equation)

Upon observing $A$, the goal is to consistently estimate the matrices $Z$ and
$Y$. Conditional on $Z$ and $Y$, writing $A=ZBY^{{\mathrm{\scriptscriptstyle
T}}}+\left(A-ZBY^{{\mathrm{\scriptscriptstyle T}}}\right)$ illustrates signal-
plus-noise matrix model form. Focusing here on $Z$ for the sake of discussion,
Rohe & Zeng (2023) proposes and studies estimators of the form
$\widehat{Z}=n^{1/2}\widehat{U}R_{\widehat{U}}$, where the orthonormal matrix
$\widehat{U}$ consists of the top left singular vectors of $A$ and the
orthogonal matrix $R_{\widehat{U}}$ maximizes the varimax objective function

$\operatorname{Varimax}\left(R;\widehat{U}\right)=\sum_{j=1}^{d}\frac{1}{n}\sum_{i=1}^{n}\left[(\widehat{U}R)_{ij}^{4}-\left\\{\frac{1}{n}\sum_{k=1}^{n}(\widehat{U}R)_{kj}^{2}\right\\}^{2}\right]$
(equation)

over the set of $d$-dimensional orthogonal matrices, denoted by
$R\in\mathcal{O}(d)=\left\\{R\in\mathbb{R}^{d\times
d}:RR^{{\mathrm{\scriptscriptstyle T}}}=R^{{\mathrm{\scriptscriptstyle
T}}}R=I_{d}\right\\}$. Solutions $R_{\widehat{U}}$, deemed varimax factor
rotations, are by inspection only unique up to multiplication by signed
permutation matrices, the latter of which are denoted by
$\Pi\in\mathcal{P}(d)=\left\\{\Pi\in\mathcal{O}(d):\Pi_{ij}\in\\{-1,0,1\\}\right\\}$.
Moreover, since $\widehat{U}$ consists of orthonormal columns, maximizing Eq.
equation is equivalent to maximizing the expression
$\|\widehat{U}R\|_{\ell_{4}}^{4}$ as a function of $R$, again over the set of
$d\times d$ orthogonal matrices.

For certain latent variable statistical matrix models, including but not
limited to semi-sparse stochastic blockmodels to be specified in Section
section, Rohe & Zeng (2023) establishes that the $n\times d$ matrix of latent
factors $Z$ is well-estimated by $\widehat{Z}$ in the sense that the uniform
row-wise estimation error asymptotically vanishes with high probability as
$n\rightarrow\infty$ when $m\asymp n$, i.e.,

$\underset{\Pi\in\mathcal{P}(d)}{\inf}\hskip 1.99997pt\underset{1\leq i\leq
n}{\max}\|\widehat{Z}_{i}-\Pi^{{\mathrm{\scriptscriptstyle T}}}\cdot
Z_{i}\|_{\ell_{2}}=\operatorname{o}_{\operatorname{P}}(1).$ (equation)

In fact, more is true, since Rohe & Zeng (2023) provides explicit, context-
specific big-O bounds that hold with high probability. To be clear, these
bounds control the maximum Euclidean row norm, denoted alternatively at times
by $\ell_{2,\infty}$ or $\|\cdot\|_{2,\infty}$ or
$\|\cdot\|_{2\rightarrow\infty}$. Additional comments are provided in Section
section.subsection.

Our point of departure is Eq. equation. We seek to quantify the asymptotic
distribution and covariance structure of individual vectors in varimax-
transformed Euclidean point clouds corresponding to spectral embeddings of
networks. Our subsequent usage of the breve symbol to write $\breve{Z}_{i}$ in
Section section plays the role of $Z_{i}$ in Section section provided the
latter is assumed to be appropriately scaled for identifiability.

##  section Main results

###  section.subsection Setup

In an effort to streamline the readability of the main results, the following
three subsections are devoted to stand-alone statements of three main
theorems. Our asymptotic results hold for undirected, directed, and degree-
corrected stochastic blockmodel graphs (Holland et al., 1983; Karrer & Newman,
2011) as specified below. Throughout, latent cluster memberships are
represented using the canonical set of $k$-dimensional standard basis vectors,
denoted by $\\{e_{1},\dots,e_{k}\\}$.

Concretely, this paper considers the following algorithm for obtaining low-
dimensional varimax-rotated data embeddings from large observed data matrices.

  1. Algorithm 1: Dimensionality reduction for blockmodel-type data using varimax.

  2. –

Input: An adjacency matrix $A\in\\{0,1\\}^{n\times m}$ and embedding dimension
value $r\ll\min(n,m)$.

  3. enumi.

Compute the best rank-$r$ approximation of $A$. Let
$\widehat{U}\in\mathbb{R}^{n\times r}$ denote the top-$r$ left singular
vectors of $A$, and let $\widehat{V}\in\mathbb{R}^{m\times r}$ denote the
top-$r$ right singular vectors of $A$.

  4. enumi.

Separately for $\widehat{U}$ and $\widehat{V}$, compute an orthogonal matrix
that maximize the varimax objective function, i.e.,

$R_{\widehat{U}}\in\underset{R\in\mathcal{O}(r)}{\arg\max}\;\operatorname{Varimax}(R;\widehat{U}),\qquad
R_{\widehat{V}}\in\underset{R\in\mathcal{O}(r)}{\arg\max}\;\operatorname{Varimax}(R;\widehat{V}).$

  5. –

Output: The estimated latent factor matrices
$\widehat{Z}=n^{1/2}\widehat{U}R_{\widehat{U}}$ and
$\widehat{Y}=m^{1/2}\widehat{V}R_{\widehat{V}}$.

  6. (Optional: The value $r$ could instead be selected in a data-driven manner, for example, based on inspecting the singular values of $A$.)

In this paper, $\widehat{Z}$ and $\widehat{Y}$ are viewed as statistical
estimators, and their row-wise distributional properties are quantified in the
large-data limit. In practice, after obtaining $\widehat{Z}$ and
$\widehat{Y}$, analysts might choose to visualize them using pair-pair plots,
to apply a clustering algorithm to their rows, or to use them as plug-in
estimates, depending on the problem at hand.

The above algorithm is sufficient for the purposes of this paper. We emphasize
that Algorithm 1 represents a special case of the more general vsp algorithm
developed in Rohe & Zeng (2023). Therein, the authors provide supporting
discussion and guidance regarding the utility of optional centering,
recentering, and scaling steps that are appropriate for more general data
settings but are beyond the scope of the present paper.

###  section.subsection Undirected stochastic blockmodels

Fix $k\geq 2$. Let $\pi$ be a $k$-dimensional probability vector with positive
entries summing to unity. Let $B\in(0,1)^{k\times k}$ be a full-rank symmetric
matrix. Given any positive integer $n$, for all $i\in\\{1,\dots,n\\}$ generate
independent, identically distributed $k$-dimensional latent indicator vectors
$Z_{i}\sim\mathcal{G}=\operatorname{Multinomial}(1;\pi)$, arranged in the
matrix $Z=[Z_{1},\dots,Z_{n}]^{{\mathrm{\scriptscriptstyle
T}}}\in\\{0,1\\}^{n\times k}$.

Conditional on $Z$, sample a symmetric adjacency matrix
$A\in\\{0,1\\}^{n\times n}$ whose upper triangular entries are independent
Bernoulli random variables with edge probabilities given by
$\operatorname{E}(A\mid Z)=\rho ZBZ^{{\mathrm{\scriptscriptstyle T}}}$, where
$\rho=\rho_{n}\in(0,1]$ denotes a sparsity term. Define $\widehat{Z}_{i}$ as
the $i$-th row of $\widehat{Z}=n^{1/2}\widehat{U}R_{\widehat{U}}$ obtained
from Algorithm 1. Let $\breve{Z}_{i}=\operatorname{diag}(\pi)^{-1/2}Z_{i}$.

######  Theorem section.theorem.

Assume the undirected stochastic blockmodel setting of Section
section.subsection with $n\rho=\omega(\log^{c}n)$ for some sufficiently large
constant $c>1$ where either $\rho=1$ or $\rho\rightarrow 0$ with limit
$\rho_{\infty}$. It suffices to set $c=20$. There exists a sequence of signed
permutation matrices $(\Pi_{Z})\subset\mathcal{P}(k)$ such that for any choice
of fixed index $i$, as $n\rightarrow\infty$, the sequence of random vectors of
the form

$(n\rho)^{1/2}\left\\{\Pi_{Z}\cdot\widehat{Z}_{i}-\breve{Z}_{i}\mid
Z_{i}=e_{\ell}\right\\}$ (equation)

converges in distribution to a multivariate Gaussian random vector with mean
zero and covariance matrix given by

$\left\\{\operatorname{diag}(\pi)^{-1/2}\cdot T_{B}^{-1}\cdot
J\cdot\Delta_{X}^{-1}\right\\}\cdot\operatorname{E}_{\xi}\\{g_{\rho_{\infty}}(\xi,T_{B}e_{\ell})\\}\cdot\left\\{\Delta_{X}^{-1}\cdot
J\cdot(T_{B}^{-1})^{{\mathrm{\scriptscriptstyle
T}}}\cdot\operatorname{diag}(\pi)^{-1/2}\right\\}.$ (equation)

In Eq. equation, $T_{B}=|\Lambda_{B}|^{1/2}U_{B}^{{\mathrm{\scriptscriptstyle
T}}}$ is defined in terms of the spectral decomposition of $B$ such that
$B=T_{B}^{{\mathrm{\scriptscriptstyle T}}}JT_{B}$ where
$J=\operatorname{diag}(1_{p},-1_{q})\in\mathbb{R}^{k\times k}$ when $B$ has
$p$ positive eigenvalues and $q$ negative eigenvalues. Furthermore,
$\Delta_{X}=\operatorname{E}(\xi\xi^{{\mathrm{\scriptscriptstyle T}}})$ where
$\xi\sim T_{B}\mathcal{G}$, and
$g_{\rho_{\infty}}(\xi,x)=(x^{{\mathrm{\scriptscriptstyle
T}}}J\xi)\cdot(1-\rho_{\infty}x^{{\mathrm{\scriptscriptstyle
T}}}J\xi)\cdot\xi\xi^{{\mathrm{\scriptscriptstyle T}}}$ where
$\rho_{\infty}\in\\{0,1\\}$.

###  section.subsection Directed stochastic blockmodels

Fix $k\geq 2$. Let $\pi_{Z}$ and $\pi_{Y}$ be two $k$-dimensional probability
vectors each having positive entries summing to unity. Let $B\in(0,1)^{k\times
k}$ be a full-rank matrix though not necessarily symmetric. Given any positive
integer $n$, for all $i\in\\{1,\dots,n\\}$ generate independent, identically
distributed $k$-dimensional latent indicator vectors
$Z_{i}\sim\mathcal{G}=\operatorname{Multinomial}(1;\pi_{Z})$, arranged in the
matrix $Z=[Z_{1},\dots,Z_{n}]^{{\mathrm{\scriptscriptstyle
T}}}\in\\{0,1\\}^{n\times k}$. Similarly, for all $j\in\\{1,\dots,n\\}$,
generate independent, identically distributed $k$-dimensional latent indicator
vectors $Y_{j}\sim\mathcal{H}=\operatorname{Multinomial}(1;\pi_{Y})$, arranged
in the matrix $Y=[Y_{1},\dots,Y_{n}]^{{\mathrm{\scriptscriptstyle
T}}}\in\\{0,1\\}^{n\times k}$.

Conditional on $Z$ and $Y$, sample a directed adjacency matrix
$A\in\\{0,1\\}^{n\times n}$ consisting of independent Bernoulli entries with
edge probabilities given by $\operatorname{E}(A\mid Z,Y)=\rho
ZBY^{{\mathrm{\scriptscriptstyle T}}}$, where $\rho=\rho_{n}\in(0,1]$ denotes
a sparsity term. Define $\widehat{Z}_{i}$ as the $i$-th row of
$\widehat{Z}=n^{1/2}\widehat{U}R_{\widehat{U}}$ obtained from Algorithm 1.
Similarly, define $\widehat{Y}_{j}$ as the $j$-th row of
$\widehat{Y}=n^{1/2}\widehat{V}R_{\widehat{V}}$ obtained from Algorithm 1. Let
$\breve{Z}_{i}=\operatorname{diag}(\pi_{Z})^{-1/2}Z_{i}$ and
$\breve{Y}_{j}=\operatorname{diag}(\pi_{Y})^{-1/2}Y_{j}$.

######  Theorem section.theorem.

Assume the directed stochastic blockmodel setting of Section
section.subsection with $n\rho=\omega(\log^{c}n)$ for some sufficiently large
constant $c>1$ where either $\rho=1$ or $\rho\rightarrow 0$ with limit
$\rho_{\infty}$. It suffices to set $c=20$. There exists a sequence of signed
permutation matrices $(\Pi_{Z})\subset\mathcal{P}(k)$ such that for any choice
of fixed index $i$, as $n\rightarrow\infty$, the sequence of random vectors of
the form

$(n\rho)^{1/2}\left\\{\Pi_{Z}\cdot\widehat{Z}_{i}-\breve{Z}_{i}\mid
Z_{i}=e_{\ell}\right\\}$ (equation)

converges in distribution to a multivariate Gaussian random vector with mean
zero and covariance matrix given by

$\left\\{\operatorname{diag}(\pi_{Z})^{-1/2}\cdot
T_{\Xi}^{-1}\cdot\Delta_{\Upsilon}^{-1}\right\\}\cdot\operatorname{E}_{\upsilon}\left\\{g_{\rho_{\infty}}(\upsilon,T_{\Xi}e_{\ell})\right\\}\cdot\left\\{\Delta_{\Upsilon}^{-1}\cdot(T_{\Xi}^{-1})^{{\mathrm{\scriptscriptstyle
T}}}\cdot\operatorname{diag}(\pi_{Z})^{-1/2}\right\\}.$ (equation)

Simultaneously, there exists a sequence of signed permutation matrices
$(\Pi_{Y})\subset\mathcal{P}(k)$ such that for any choice of fixed index $j$,
as $n\rightarrow\infty$, the sequence of random vectors of the form

$(n\rho)^{1/2}\left\\{\Pi_{Y}\cdot\widehat{Y}_{j}-\breve{Y}_{j}\mid
Y_{j}=e_{\ell}\right\\}$ (equation)

converges in distribution to a multivariate Gaussian random vector with mean
zero and covariance matrix given by

$\left\\{\operatorname{diag}(\pi_{Y})^{-1/2}\cdot
T_{\Upsilon}^{-1}\cdot\Delta_{\Xi}^{-1}\right\\}\cdot\operatorname{E}_{\xi}\left\\{g_{\rho_{\infty}}(\xi,T_{\Upsilon}e_{\ell})\right\\}\cdot\left\\{\Delta_{\Xi}^{-1}\cdot(T_{\Upsilon}^{-1})^{{\mathrm{\scriptscriptstyle
T}}}\cdot\operatorname{diag}(\pi_{Y})^{-1/2}\right\\}.$ (equation)

In Eqs. equation and equation, the matrices of coefficients
$T_{\Xi}=S_{B}^{1/2}U_{B}^{{\mathrm{\scriptscriptstyle T}}}$ and
$T_{\Upsilon}=S_{B}^{1/2}V_{B}^{{\mathrm{\scriptscriptstyle T}}}$ are defined
in terms of the singular value decomposition of $B$, written
$B=U_{B}S_{B}V_{B}^{{\mathrm{\scriptscriptstyle T}}}$. Here,
$\Delta_{\Xi}=\operatorname{E}(\xi\xi^{{\mathrm{\scriptscriptstyle T}}})$
where $\xi\sim T_{\Xi}\mathcal{G}$,
$\Delta_{\Upsilon}=\operatorname{E}(\upsilon\upsilon^{{\mathrm{\scriptscriptstyle
T}}})$ where $\upsilon\sim T_{\Upsilon}\mathcal{H}$, and
$g_{\rho_{\infty}}(\alpha,\beta)=(\beta^{{\mathrm{\scriptscriptstyle
T}}}\alpha)\cdot(1-\rho_{\infty}\beta^{{\mathrm{\scriptscriptstyle
T}}}\alpha)\cdot\alpha\alpha^{{\mathrm{\scriptscriptstyle T}}}$ with
$\rho_{\infty}\in\\{0,1\\}$ takes $k$-dimensional input vectors $\alpha$ and
$\beta$.

In general, the asymptotic covariance matrices differ for Eq. equation and Eq.
equation. The asymptotic distribution of Eq. equation depends on the
distribution $\mathcal{H}$, and similarly, Eq. equation depends on the
distribution $\mathcal{G}$, manifest in the covariance expressions Eqs.
equation and equation.

For simplicity, Theorem section.theorem is stated and proved for large, square
adjacency matrices. It can be shown that asymptotic normality continues to
hold for rectangular adjacency matrices $A\in\\{0,1\\}^{n\times m}$ in the
large-system asymptotic regime $n,m\rightarrow\infty$ when $n$ and $m$ grow at
the same rate, requiring only minor modifications to the current notation and
proof.

###  section.subsection Degree-corrected stochastic blockmodels

Fix $k\geq 2$. Let $\pi$ be a $k$-dimensional probability vector with positive
entries summing to unity. Let $B\in(0,1)^{k\times k}$ be a full-rank symmetric
matrix. Given any positive integer $n$, for all $i\in\\{1,\dots,n\\}$ generate
independent, identically distributed $k$-dimensional latent indicator vectors
$Z_{i}\sim\mathcal{G}=\operatorname{Multinomial}(1;\pi)$, arranged in the
matrix $Z=[Z_{1},\dots,Z_{n}]^{{\mathrm{\scriptscriptstyle
T}}}\in\\{0,1\\}^{n\times k}$. Independently of $Z$, generate independent,
identically distributed scalar degree parameters $\theta_{i}\sim\mathcal{H}$
from a distribution $\mathcal{H}$ with compact support contained in the open
unit interval, to form the vector of degree-heterogeneity parameters
$\theta=(\theta_{1},\dots,\theta_{n})^{{\mathrm{\scriptscriptstyle T}}}$.

Conditional on $Z$ and $\theta$, sample a symmetric adjacency matrix
$A\in\\{0,1\\}^{n\times n}$ whose upper triangular entries are independent
Bernoulli random variables with edge probabilities given by
$\operatorname{E}(A\mid
Z,\theta)=\rho\operatorname{diag}(\theta)ZBZ^{{\mathrm{\scriptscriptstyle
T}}}\operatorname{diag}(\theta)$, where $\rho=\rho_{n}\in(0,1]$ denotes a
sparsity term. Define $\widehat{Z}_{i}$ as the $i$-th row of
$\widehat{Z}=n^{1/2}\widehat{U}R_{\widehat{U}}$ obtained from Algorithm 1. Let
$\breve{Z}_{i}^{\prime}=\operatorname{diag}(\eta)^{-1/2}\cdot(\theta_{i}Z_{i})$
where $\eta=\operatorname{E}(\theta_{1}^{2}Z_{1}^{2})\in\mathbb{R}^{k}$.

######  Theorem section.theorem.

Assume the degree-corrected stochastic blockmodel setting of Section
section.subsection with $n\rho=\omega(\log^{c}n)$ for some sufficiently large
constant $c>1$ where either $\rho=1$ or $\rho\rightarrow 0$ with limit
$\rho_{\infty}$. It suffices to set $c=20$. There exists a sequence of signed
permutation matrices $(\Pi_{Z})\subset\mathcal{P}(k)$ such that, for any
choice of fixed index $i$ and any vector $a\in\mathbb{R}^{k}$, it holds that

$\operatorname{pr}\left\\{(n\rho)^{1/2}\left(\Pi_{Z}\cdot\widehat{Z}_{i}-\breve{Z}_{i}^{\prime}\right)\leq
a\right\\}\longrightarrow\int_{\operatorname{supp}\mathcal{F}}\Phi\left\\{a,\Gamma_{\rho_{\infty}}(b)\right\\}\operatorname{d}\mathcal{F}(b)$
(equation)

in the large-$n$ limit, where $\mathcal{F}$ denotes the product distribution
$\mathcal{G}\mathcal{H}$ and $\Phi\left\\{\cdot,\cdot\right\\}$ denotes the
mean zero multivariate Gaussian cumulative distribution function evaluated at
the vector $a\in\mathbb{R}^{k}$ with covariance matrix
$\Gamma_{\rho_{\infty}}(b)$, for $b\in\operatorname{supp}\mathcal{F}$, given
by

$\left\\{\operatorname{diag}(\eta)^{-1/2}\cdot T_{B}^{-1}\cdot
J\cdot\Delta_{X^{\prime}}^{-1}\right\\}\cdot\operatorname{E}_{\xi^{\prime}}\\{g_{\rho_{\infty}}(\xi^{\prime},b)\\}\cdot\left\\{\Delta_{X^{\prime}}^{-1}\cdot
J\cdot(T_{B}^{-1})^{{\mathrm{\scriptscriptstyle
T}}}\cdot\operatorname{diag}(\eta)^{-1/2}\right\\}.$ (equation)

Above, $T_{B}=|\Lambda_{B}|^{1/2}U_{B}^{{\mathrm{\scriptscriptstyle T}}}$
satisfies $B=T_{B}^{{\mathrm{\scriptscriptstyle T}}}JT_{B}$ where
$J=\operatorname{diag}(1_{p},-1_{q})\in\mathbb{R}^{k\times k}$ when $B$ has
$p$ positive eigenvalues and $q$ negative eigenvalues. Further,
$\Delta_{X^{\prime}}=\operatorname{E}(\xi^{\prime}\xi^{\prime{\mathrm{\scriptscriptstyle
T}}})$ where $\xi^{\prime}\sim T_{B}\mathcal{G}\mathcal{H}$, and
$g_{\rho_{\infty}}(\xi^{\prime},b)=(b^{{\mathrm{\scriptscriptstyle
T}}}J\xi^{\prime})\cdot(1-\rho_{\infty}b^{{\mathrm{\scriptscriptstyle
T}}}J\xi^{\prime})\cdot\xi^{\prime}\xi^{\prime{\mathrm{\scriptscriptstyle
T}}}$ where $\rho_{\infty}\in\\{0,1\\}$.

###  section.subsection Coda to main results

Several remarks are in order to further contextualize and clarify the main
results. First, regarding notation, the literature on stochastic blockmodels
commonly writes $Z_{i}$ to denote the unscaled community membership vector for
node $i$ taking values in the set of $k$-dimensional standard basis vectors.
Above, our notation $\breve{Z}_{i}$ indicates that the community membership
vector is scaled for the sake of identifiability, though for simplicity we
drop the breve symbol in some of the surrounding discussion here. Second,
writing the entrywise matrix statement $\widehat{Z}\approx Z\Pi$ corresponds
to writing the column vector statement
$\widehat{Z}_{i}\approx\Pi^{{\mathrm{\scriptscriptstyle T}}}Z_{i}$ and hence
$\Pi\widehat{Z}_{i}\approx Z_{i}$ uniformly for all vector indices
$i\in\\{1,\dots,n\\}$.

In the main theorems, each sequence of signed permutation matrices
$(\Pi)\subset\mathcal{P}(k)$ is implicitly indexed by $n$ and does not depend
on the choice of node index $i$. The need to account for these matrices is
unavoidable in general and discussed further in Section section.

The big-O rate established in Rohe & Zeng (2023) for stochastic blockmodels,
which implies the uniform error bound in Eq. equation, is approximately
$\operatorname{O}_{\operatorname{P}}\left\\{(n\rho)^{-1/4}\right\\}$, ignoring
logarithmic terms. In contrast, the results in this paper suggest the non-
uniform row-wise $\ell_{2}$ norm bound
$\operatorname{O}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$, ignoring
logarithmic terms, noting that asymptotic normality is driven by the
stabilization of a linear combination of independent latent bounded random
variables. To be clear, this paper obtains asymptotic normality for individual
low-dimensional embedding vectors, thereby building on the results in Rohe &
Zeng (2023), but without obtaining an improved uniform error bound for Eq.
equation.

The independence assumption between $Z_{i}$ and $Y_{i}$ in each tuple
$(Z_{i},Y_{i})$ in Section section.subsection is conventional, though it is
perhaps possible that versions of Eqs. equation and equation continue to hold
under certain forms of dependence. The choice of constant $c=20$ in the growth
condition $n\rho=\omega(\log^{c}n)$, i.e.,
$n\rho/(\log^{c}n)\rightarrow\infty$ where $\log^{c}n=(\log n)^{c}$, is made
for convenience and is based on the perturbation analysis collectively found
in Cape et al. (2019a); Tang et al. (2022); Rohe & Zeng (2023). The asymptotic
covariance matrices depend on the graph sparsity regime
$\rho_{\infty}\in\\{0,1\\}$ and convey larger asymptotic variability when
$\rho_{\infty}=0$ compared to the dense graph regime $\rho_{\infty}=1$.

A key challenge in going from the $\ell_{2,\infty}$ bound in Rohe & Zeng
(2023) to a row-wise distributional guarantee is to decompose the matrix
difference between the sample quantity and the population-level quantity,
namely the difference between the estimated latent factors and ground-truth
latent factors, into tractable dominant and residual terms that reflect the
underlying signal-plus-noise structure. Beyond simply identifying the dominant
leading-order behavior in the perturbation analysis, it must then be
appropriately oriented and decomposed so as to obtain explicit formulas for
the asymptotic mean vector and covariance matrix in terms of the data-
generating model parameters. Non-trivial bookkeeping and analysis involving
various orthogonal transformations is required. Further details and
discussions can be found in the proofs which are provided in the Supplementary
Material.

##  section Numerical examples

###  section.subsection Four-block undirected stochastic blockmodel graphs

Here, we simulate $100$ undirected stochastic blockmodel graphs as defined in
Section section.subsection, each having $n=5000$ nodes, sparsity term
$\rho=1$, and with block connectivity matrix and membership probability vector
given by

$B=\begin{bmatrix}0.6\phantom{0}&0.2\phantom{0}&0.1\phantom{0}&0.1\phantom{0}\\\
\cdot&0.7\phantom{0}&0.05\phantom{0}&0.05\phantom{0}\\\
\cdot&\cdot&0.6\phantom{0}&0.25\phantom{0}\\\
\cdot&\cdot&\cdot&0.6\phantom{0}\end{bmatrix},\qquad\pi=(0.25,0.25,0.25,0.25)^{{\mathrm{\scriptscriptstyle
T}}},$

inspired by the simulation study in Chen & Rohe (2020). We separately embed
each graph adjacency matrix $A^{(j)}$ into dimension $k=4$ to obtain a
varimax-transformed matrix estimate $\widehat{Z}^{(j)}\in\mathbb{R}^{n\times
k}$ using the function varimax in the R package `stats`. Given each graph
embedding, we compute estimated class-conditional covariance matrices
$\widehat{\Sigma}^{(j)}_{\ell}$ for each cluster label $\ell\in\\{1,2,3,4\\}$
by fitting a one-component Gaussian density to the collection of vectors
$(n\rho)^{1/2}\left\\{\Pi_{Z}\cdot\widehat{Z}^{(j)}_{i}-\breve{Z}_{i}\right\\}$
satisfying $Z_{i}=e_{\ell}$. Our implementation uses the
expectation–maximization algorithm found in the R package `mclust`. In this
simulation setting, the signed permutation matrix $\Pi_{Z}$ can be computed
using knowledge of the true latent factor matrix $Z$, though it is only needed
to ensure proper ordering of coordinates and signs for comparison with the
theoretical asymptotic covariance matrices. Here, Section section.subsection
reports the block-wise averages of the relative errors
$\|\widehat{\Sigma}^{(j)}_{\ell}-\Sigma_{\ell}\|_{\operatorname{F}}/\|\Sigma_{\ell}\|_{\operatorname{F}}$
computed in Frobenius norm, $1\leq\ell\leq k$, together with sample standard
errors, each scaled by a factor of one hundred for readability. As a further
illustration, compare for example $\widehat{\Sigma}^{(1)}_{1}$ with
$\Sigma_{1}$ and $\widehat{\Sigma}^{(1)}_{2}$ with $\Sigma_{2}$ below.

$\displaystyle\Sigma_{1}\overset{\cdot}{=}\begin{bmatrix}\phantom{-}14.7&-5.84&-1.71&-1.71\\\
\cdot&\phantom{-}7.42&\phantom{-}0.196&\phantom{-}0.196\\\
\cdot&\cdot&\phantom{-}7.11&-4.65\\\
\cdot&\cdot&\cdot&\phantom{-}7.11\end{bmatrix},$ $\displaystyle\hskip
10.00002pt\widehat{\Sigma}^{(1)}_{1}\overset{\cdot}{=}\begin{bmatrix}\phantom{-}14.4&-6.00&-1.75&-1.69\\\
\cdot&\phantom{-}7.75&\phantom{-}0.372&\phantom{-}0.176\\\
\cdot&\cdot&\phantom{-}6.67&-4.29\\\
\cdot&\cdot&\cdot&\phantom{-}6.72\end{bmatrix},$
$\displaystyle\Sigma_{2}\overset{\cdot}{=}\begin{bmatrix}\phantom{-}10.2&-5.25&-1.07&-1.07\\\
\cdot&\phantom{-}9.08&\phantom{-}0.054&\phantom{-}0.054\\\
\cdot&\cdot&\phantom{-}3.77&-2.43\\\
\cdot&\cdot&\cdot&\phantom{-}3.77\end{bmatrix},$ $\displaystyle\hskip
10.00002pt\widehat{\Sigma}^{(1)}_{2}\overset{\cdot}{=}\begin{bmatrix}\phantom{-}11.4&-5.92&-1.04&-1.43\\\
\cdot&\phantom{-}9.41&\phantom{-}0.215&\phantom{-}0.079\\\
\cdot&\cdot&\phantom{-}3.32&-2.10\\\
\cdot&\cdot&\cdot&\phantom{-}3.64\end{bmatrix}.$

Strictly speaking, for each stochastic blockmodel graph, Theorem
section.theorem does not simultaneously hold jointly for all order
$\operatorname{O}(n)$ embedding vectors in each block $\ell$. In particular,
there is dependence in the order $\operatorname{O}(n)$ vectors used to
construct $\widehat{\Sigma}^{(j)}_{\ell}$, and we do not claim that the
estimation approach used in the current numerical example is optimal for
estimating the asymptotic variances and covariances. Nevertheless, embedding
the graphs and then aggregating the node embeddings follows conventional
spectral graph clustering methodology, and doing so produces reasonably close
agreement between empirical and theoretical covariances as shown in Section
section.subsection.

Estimation of asymptotic covariance matrices in Section section.subsection
Block 1 Block 2 Block 3 Block 4 Rel.0Frob.0error (sample stand.0err.) $\times
100$ 6.198 (0.2491) 6.383 (0.2505) 5.506 (0.2352) 5.424 (0.2460)

###  section.subsection Two-block directed stochastic blockmodel graphs

Here, we simulate directed stochastic blockmodel graphs according to Section
section.subsection, with left block membership vector
$\pi_{Z}=(0.25,0.75)^{{\mathrm{\scriptscriptstyle T}}}$, right block
membership vector $\pi_{Y}=(0.667,0.333)^{{\mathrm{\scriptscriptstyle T}}}$,
sparsity term $\rho=1$, and block connectivity matrix defined entrywise as
$B_{11}=0.4,B_{12}=0.6,B_{21}=0.3,B_{22}=0.7$. For two independent graphs of
size $n=10^{3}$ and $n=10^{4}$, respectively, Fig. figure depicts the varimax-
transformed graph embeddings for the left and right estimated latent factor
matrices. Crucially, the point cloud centroids lie along the standard
coordinate axes, however each block exhibits elliptic radial covariance
structure whose axes do not simultaneously align with the standard axes. In
addition, each theoretical covariance matrix is reasonably well estimated in
relative Frobenius norm by its sample counterpart.

(a) Vectors $\widehat{Z}_{i}$ when $n=10^{3}$. (b) Vectors $\widehat{Y}_{i}$
when $n=10^{3}$.

(c) Vectors $\widehat{Z}_{i}$ when $n=10^{4}$. (d) Vectors $\widehat{Y}_{i}$
when $n=10^{4}$.

Figure figure: Varimax-rotated left and right embeddings for two directed
graphs. Embedding vectors are uncentered. Latent vectors are scaled in
agreement with Section section.subsection. Dotted lines denote circles with
radii given by the reciprocals of the square-roots of the block membership
probabilities. Solid black circles denote the sample mean vectors computed for
each block. Additional details are provided in Section section.subsection.

###  section.subsection Degree-corrected stochastic blockmodels with affinity
network structure

Here, we simulate a single dense three-block degree-corrected stochastic
blockmodel graph per Section section.subsection, with two thousand nodes,
equal membership probabilities, and block connectivity matrix defined
entrywise as $B_{ii}=0.2$ for all $1\leq i\leq 3$ and $B_{ij}=0.1$ for all
$i\neq j$. The within-block edge probabilities are larger than the between-
block edge probabilities which can be interpreted as producing affinity
network structure. The node-specific degree parameters are drawn independently
in the manner $\theta_{i}\sim\operatorname{Uniform}[0.25,0.75]$ for all $1\leq
i\leq n$. For each panel in Fig. figure, two of the three latent factor
dimensions are isolated along the standard coordinate axes, with the one
remaining dimension yielding coordinate values near the origin. This behavior
generalizes to the same model having $k\geq 3$ blocks, illustrating the
concept of “simple structure” (Thurstone, 1947) as well as the concept of
approximate sparsity in embedding space after performing dimensionality
reduction.

(a) View of dimensions one and two (b) View of dimensions two and three

Figure figure: Varimax-rotated degree-corrected stochastic blockmodel
embedding for one graph. Dotted circles have radius equal to the magnitude of
each block’s non-zero theoretical coordinate centroid, while solid black
circles denote empirical coordinate centroids. Additional details are provided
in Section section.subsection.

###  section.subsection Asymptotic variance and covariance expressions

The asymptotic covariance expressions in Theorems section.theorem,
section.theorem and section.theorem are rather complicated in general. For
simplicity and concreteness, consider undirected stochastic blockmodel graphs
generated from the three-parameter model

$B=\begin{bmatrix}a&b\\\ \cdot&a\end{bmatrix},\quad
0<b<a<1;\qquad\pi=(\pi_{1},1-\pi_{1})^{{\mathrm{\scriptscriptstyle T}}},\quad
0<\pi_{1}<1.$

In this setting, it happens to be possible to analytically compute the block-
specific asymptotic covariance matrices. In particular, for the first block,
namely for the first community of nodes, a calculation reveals that

$\Sigma_{1}=\begin{bmatrix}\frac{a^{3}(1-a)(1-\pi_{1})+b^{3}(1-b)\pi_{1}}{(a^{2}-b^{2})^{2}\pi_{1}^{2}(1-\pi_{1})}&\frac{-ab\\{a(1-a)+(a-b)(a+b-1)\pi_{1}\\}}{(a^{2}-b^{2})^{2}\\{\pi_{1}(1-\pi_{1})\\}^{3/2}}\\\
\cdot&\frac{a^{2}b\pi_{1}+ab^{2}(1-a-\pi_{1})}{(a^{2}-b^{2})^{2}\pi_{1}(1-\pi_{1})^{2}}\end{bmatrix}.$

For the second block, a similar calculation reveals that

$\Sigma_{2}=\begin{bmatrix}\frac{ab^{2}\pi_{1}+a^{2}b(1-b-\pi_{1})}{(a^{2}-b^{2})^{2}\pi_{1}^{2}(1-\pi_{1})}&\frac{-ab\left\\{b(1-b)+(a-b)(1-a-b)\pi_{1}\right\\}}{(a^{2}-b^{2})^{2}\left\\{\pi_{1}(1-\pi_{1})\right\\}^{3/2}}\\\
\cdot&\frac{a^{3}(1-a)\pi_{1}+b^{3}(1-b)(1-\pi_{1})}{(a^{2}-b^{2})^{2}\pi_{1}(1-\pi_{1})^{2}}\end{bmatrix}.$

The determinants of the above covariance matrices, namely the generalized
variances of the respective distributions, are given by

$\det(\Sigma_{1})=\det(\Sigma_{2})=\frac{a(1-a)\cdot
b(1-b)}{(a^{2}-b^{2})^{2}\cdot\pi_{1}^{2}(1-\pi_{1})^{2}}.$

The above expression succinctly reflects the near-degeneracy of the asymptotic
covariance matrices in situations where the underlying stochastic blockmodel
is itself nearly degenerate, namely when either $b\approx 0$, $b\approx a$,
$a\approx 1$, $\pi_{1}\approx 0$, or $\pi_{1}\approx 1$.

Evaluating the above expressions for several choices of model parameters
yields the following examples of diverse covariance structure.

$\displaystyle a=\tfrac{3}{4},b=\tfrac{1}{4},\pi_{1}=\tfrac{1}{4},$
$\displaystyle\qquad\Sigma_{1}=\begin{bmatrix}7&-3^{1/2}\\\
\cdot&\phantom{-}1\phantom{{}^{1/2}}\end{bmatrix},$
$\displaystyle\Sigma_{2}=\Sigma_{1};$ $\displaystyle
a=\tfrac{4}{5},b=\tfrac{3}{5},\pi_{1}=\tfrac{1}{2},$
$\displaystyle\qquad\Sigma_{1}\overset{\cdot}{=}\begin{bmatrix}9.63&-9.79\\\
\cdot&\phantom{-}10.77\end{bmatrix},$
$\displaystyle\Sigma_{2}\overset{\cdot}{=}\begin{bmatrix}10.77&-9.79\\\
\cdot&\phantom{-}9.63\end{bmatrix};$ $\displaystyle
a=\tfrac{8}{10},b=\tfrac{3}{10},\pi_{1}=\tfrac{4}{10},$
$\displaystyle\qquad\Sigma_{1}\overset{\cdot}{=}\begin{bmatrix}2.37&-1.21\\\
\cdot&\phantom{-}1.43\end{bmatrix},$
$\displaystyle\Sigma_{2}\overset{\cdot}{=}\begin{bmatrix}2.97&-1.28\\\
\cdot&\phantom{-}1.20\end{bmatrix}.$

##  section Extensions

This paper considers stochastic blockmodel graphs having full-rank $k\times k$
connectivity matrix $B$ which is a widespread assumption in the literature. If
instead $\operatorname{rank}(B)=r<k$, then modified versions of the results in
Section section possibly still hold, though in such cases the $k$-dimensional
singular subspaces of the latent factor matrices $Z$ and $Y$ are mapped into
lower $r$-dimensional space for the population-level expectation term
$ZBY^{{\mathrm{\scriptscriptstyle T}}}$, resulting in quantitatively different
varimax factor rotations and estimation properties. Broadly speaking,
understanding such rank-degenerate blockmodeling is desirable because it
permits flexibility in the true number of embedding dimensions $r$ and the
true number of blocks or clusters $k$.

This paper establishes the asymptotic normality of individual vectors
representing node embeddings. Recent related work establishes that
standardized averages of spectral estimates for the entries of $B$ also
exhibit asymptotic normality but with bias in the presence of rank degeneracy
that is exacerbated by network sparsity (Tang et al., 2022). Similar behavior
is anticipated in the present paper when estimating the non-zero components of
membership vectors via the standardized block-conditional averages of the
vectors $\widehat{Z}_{i}^{(j)}$; indeed, for the data in Section
section.subsection, conducting a per-block Shapiro–Wilk test for normality
yields $\operatorname{p}$-values larger than $0.2$.

The findings in Section section possibly extend to non-blockmodel special
cases of the more expressive generalized random dot product graph setting in
Rubin-Delanchy et al. (2022) which permits more flexible latent variable
modeling of networks beyond stochastic blockmodels. A challenge with such
general, flexible models regarding varimax factor rotations, however, is the
typical absence of even approximately sparse latent factor structure and the
difficulty of establishing identifiability conditions for estimation.

Rohe & Zeng (2023) develops theory and methods for more general data-
generating models than stochastic blockmodels. Additional steps such as data
centering, recentering, and row-normalized varimax are discussed therein. In
future work, it would be interesting to quantify the asymptotic distributional
properties of varimax-based estimators in such settings, which foreseeably
would introduce additional complexities and technical challenges.

##  section Discussion

Stochastic blockmodels have attracted significant research interest in recent
years at the time of writing this paper. Despite their seemingly
straightforward structure, the statistics community is only now developing a
more complete understanding of spectral-based inference for these baseline
statistical network models (Rohe et al., 2011; Lei & Rinaldo, 2015; Lei, 2016;
Cape et al., 2019a; Tang et al., 2022; Athreya et al., 2022). While stochastic
blockmodels are of moderate interest in and of themselves, studying them
ultimately provides points of departure for investigating more complex
relational data and for moving beyond node clustering or the so-called
community detection problem.

Readers familiar with the statistical network analysis literature on spectral
embedding will notice that the results in this paper somewhat resemble those
found in several existing works, e.g., Athreya et al. (2016); Tang & Priebe
(2018); Rubin-Delanchy et al. (2022). These related works and others primarily
establish the uniform asymptotic recovery of node embeddings and asymptotic
normality for sequences of vectors acted on by unknown linear transformations,
typically indefinite or orthogonal transformations, which admit no general
convergence properties or obvious interpretability. Consequently, such results
are more useful for informing clustering methodology than for direct
estimation in practice. In contrast, our results incorporate varimax factor
rotations with graph embeddings to enable more straightforward estimation
modulo only the set of signed permutation matrices $\mathcal{P}(\cdot)$.
Advantageously, the collection $\mathcal{P}(\cdot)$ reflects the minimal
possible necessary statistical non-uniqueness in general, since singular
vectors in left-right pairs are uniquely computed only up to a choice of sign
and since blockmodels are unique only up to globally permuting the latent
cluster labels or the block coordinates.

In the undirected and degree-corrected stochastic blockmodel numerical
examples, nodes within the same graph share a common expected degree, hence
preventing trivial degree-based node clustering. In the directed stochastic
blockmodel numerical example, separately applying varimax to the estimated
left and right singular vectors is shown to successfully recover left-right
directional membership heterogeneity. Our simulations make use of the true
latent block memberships which are asymptotically perfectly recoverable via
$k$-means clustering in the present semi-sparse regime
$n\rho=\omega(\log^{c}n)$. Moreover, the embedding dimension can be
asymptotically perfectly recovered via an elbow in the scree plot, namely via
a demonstrable spectral gap.

This paper considers computing a truncated singular value decomposition of a
given matrix $A$ to obtain estimated latent factor matrices $\widehat{Z}$ and
$\widehat{Y}$ for $Z$ and $Y$ appearing in $\operatorname{E}(A\mid Z,Y)$, per
Eq. equation. As such, the methods and results in this paper correspond to the
so-called problem of matrix denoising. Matrix denoising is distinct from the
problem of covariance estimation in factor models via principal components
analysis. In particular, this paper allows for heteroskedasticity in the
underlying noise matrix, whereas it is well-known that heteroskedastic noise
can pose significant challenges for classical principal components analysis
and factor analysis. Further discussion and results to this effect can be
found in, for example, Zhang et al. (2022); Yan et al. (2021).

Unanswered questions abound. For example, the statistical properties of other
common factor rotation methods, such as promax, quartimax, and equimax, are
not fully developed, even for semiparametric factor models. The present paper
does not consider the optional row-normalization step in varimax which remains
popular with practitioners and whose theoretical properties are incompletely
understood. For networks specifically, this paper considers low-dimensional
embeddings derived from adjacency matrices, while the widespread success and
popularity of graph Laplacians makes them natural candidates to consider
together with varimax rotations. In this direction, it has been shown that
Laplacian-based analyses perform well empirically, even for high-dimensional,
relatively sparse data matrices, e.g., see Rohe & Zeng (2023); Chen & Rohe
(2020). That being said, their theoretical properties can be quite challenging
to analyze, e.g., see Tang & Priebe (2018). These are but a few of the many
open problems in contemporary multivariate analysis.

## Acknowledgement

The author thanks the Editor, the Associate Editor, and two referees for their
valuable feedback and suggestions. The author also thanks Karl Rohe and Minh
Tang for stimulating discussions. This research is supported in part by the
National Science Foundation under grants DMS-1902755 and SES-1951005. The
author gratefully acknowledges support from the University of
Wisconsin–Madison, Office of the Vice Chancellor for Research and Graduate
Education with funding from the Wisconsin Alumni Research Foundation.

## References

  * Abbe et al. (2020) Abbe, E., Fan, J., Wang, K. & Zhong, Y. (2020).  Entrywise eigenvector analysis of random matrices with low expected rank.  Annals of Statistics 48, 1452–1474.
  * Alter et al. (2000) Alter, O., Brown, P. O. & Botstein, D. (2000).  Singular value decomposition for genome-wide expression data processing and modeling.  Proceedings of the National Academy of Sciences 97, 10101–10106.
  * Anderson (1962) Anderson, T. W. (1962).  An introduction to multivariate statistical analysis.  Wiley, 3rd ed.
  * Anderson & Rubin (1956) Anderson, T. W. & Rubin, H. (1956).  Statistical inference in factor analysis.  In Proceedings of the third Berkeley symposium on mathematical statistics and probability, vol. 5.
  * Argelaguet et al. (2018) Argelaguet, R., Velten, B., Arnol, D., Dietrich, S., Zenz, T., Marioni, J. C., Buettner, F., Huber, W. & Stegle, O. (2018).  Multi-omics factor analysis: a framework for unsupervised integration of multi-omics data sets.  Molecular Systems Biology 14, e8124.
  * Athreya et al. (2022) Athreya, A., Cape, J. & Tang, M. (2022).  Eigenvalues of stochastic blockmodel graphs and random graphs with low-rank edge probability matrices.  Sankhya A 84, 36–63.
  * Athreya et al. (2016) Athreya, A., Priebe, C. E., Tang, M., Lyzinski, V., Marchette, D. J. & Sussman, D. L. (2016).  A limit theorem for scaled eigenvectors of random dot product graphs.  Sankhya A 78, 1–18.
  * Bai & Ng (2008) Bai, J. & Ng, S. (2008).  Large dimensional factor analysis.  Now Publishers Inc.
  * Bartholomew (1995) Bartholomew, D. J. (1995).  Spearman and the origin and development of factor analysis.  British Journal of Mathematical and Statistical Psychology 48, 211–220.
  * Brown (2015) Brown, T. A. (2015).  Confirmatory factor analysis for applied research.  Guilford Publications.
  * Cape et al. (2019a) Cape, J., Tang, M. & Priebe, C. E. (2019a).  Signal-plus-noise matrix models: eigenvector deviations and fluctuations.  Biometrika 106, 243–250.
  * Cape et al. (2019b) Cape, J., Tang, M. & Priebe, C. E. (2019b).  The two-to-infinity norm and singular subspace geometry with applications to high-dimensional statistics.  Annals of Statistics 47, 2405–2439.
  * Chamberlain & Rothschild (1983) Chamberlain, G. & Rothschild, M. (1983).  Arbitrage, factor structure, and mean-variance analysis on large asset markets.  Econometrica 51, 1281–1304.
  * Chen & Rohe (2020) Chen, F. & Rohe, K. (2020).  A new basis for sparse PCA.  arXiv preprint arXiv:2007.00596v1 .
  * Chen et al. (2021) Chen, Y., Chi, Y., Fan, J. & Ma, C. (2021).  Spectral methods for data science: a statistical perspective.  Foundations and Trends in Machine Learning 14, 566–806.
  * Fabrigar et al. (1999) Fabrigar, L. R., Wegener, D. T., MacCallum, R. C. & Strahan, E. J. (1999).  Evaluating the use of exploratory factor analysis in psychological research.  Psychological Methods 4, 272–299.
  * Fan et al. (2021) Fan, J., Wang, K., Zhong, Y. & Zhu, Z. (2021).  Robust high-dimensional factor models with applications to statistical machine learning.  Statistical Science 36, 303–327.
  * Fan et al. (2018) Fan, J., Wang, W. & Zhong, Y. (2018).  An $\ell_{\infty}$ eigenvector perturbation bound and its application.  Journal of Machine Learning Research 18, 1–42.
  * Holland et al. (1983) Holland, P. W., Laskey, K. B. & Leinhardt, S. (1983).  Stochastic blockmodels: first steps.  Social Networks 5, 109–137.
  * Holzinger & Harman (1941) Holzinger, K. J. & Harman, H. H. (1941).  Factor analysis: a synthesis of factorial methods.  University of Chicago Press.
  * Kaiser (1958) Kaiser, H. F. (1958).  The varimax criterion for analytic rotation in factor analysis.  Psychometrika 23, 187–200.
  * Karrer & Newman (2011) Karrer, B. & Newman, M. E. (2011).  Stochastic blockmodels and community structure in networks.  Physical Review E 83, 016107.
  * Lei (2016) Lei, J. (2016).  A goodness-of-fit test for stochastic block models.  Annals of Statistics 44, 401–424.
  * Lei & Rinaldo (2015) Lei, J. & Rinaldo, A. (2015).  Consistency of spectral clustering in stochastic block models.  Annals of Statistics 43, 215–237.
  * Rohe et al. (2011) Rohe, K., Chatterjee, S. & Yu, B. (2011).  Spectral clustering and the high-dimensional stochastic blockmodel.  Annals of Statistics 39, 1878–1915.
  * Rohe & Zeng (2023) Rohe, K. & Zeng, M. (2023).  Vintage factor analysis with varimax performs statistical inference.  Journal of the Royal Statistical Society Series B: Statistical Methodology (with Discussion) 85, 1037–1060.
  * Rubin-Delanchy et al. (2022) Rubin-Delanchy, P., Cape, J., Tang, M. & Priebe, C. E. (2022).  A statistical interpretation of spectral embedding: the generalised random dot product graph.  Journal of the Royal Statistical Society, Series B 84, 1446–1473.
  * Tang et al. (2022) Tang, M., Cape, J. & Priebe, C. E. (2022).  Asymptotically efficient estimators for stochastic blockmodels: the naive MLE, the rank-constrained MLE, and the spectral.  Bernoulli 28, 1049–1073.
  * Tang & Priebe (2018) Tang, M. & Priebe, C. E. (2018).  Limit theorems for eigenvectors of the normalized Laplacian for random graphs.  Annals of Statistics 46, 2360–2415.
  * Thurstone (1935) Thurstone, L. L. (1935).  The vectors of mind: multiple-factor analysis for the isolation of primary traits.  University of Chicago Press.
  * Thurstone (1947) Thurstone, L. L. (1947).  Multiple-factor analysis: a development and expansion of “The Vectors of Mind”.  University of Chicago Press.
  * Vincent (1953) Vincent, D. F. (1953).  The origin and development of factor analysis.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 2, 107–117.
  * Wedin (1972) Wedin, P.-Å. (1972).  Perturbation bounds in connection with singular value decomposition.  BIT Numerical Mathematics 12, 99–111.
  * Weyl (1912) Weyl, H. (1912).  Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Differentialgleichungen (mit einer Anwendung auf die Theorie der Hohlraumstrahlung).  Mathematische Annalen 71, 441–479.
  * Yan et al. (2021) Yan, Y., Chen, Y. & Fan, J. (2021).  Inference for heteroskedastic PCA with missing data.  arXiv preprint arXiv:2107.12365 .
  * Zhang et al. (2022) Zhang, A. R., Cai, T. T. & Wu, Y. (2022).  Heteroskedastic PCA: algorithm, optimality, and applications.  Annals of Statistics 50, 53–80.

##  section Supplementary Material

###  section.subsection Asymptotic normality in varimax-transformed undirected
graph embeddings

######  Proof section.theorem (of Theorem 1).

This proof proceeds via matrix eigendecompositions, or equivalently, spectral
decompositions. This choice is primarily for convenience; we emphasize that
the analysis also holds for matrix singular value decompositions, with only
minor notational modifications, since the matrices $A$ and $B$ of primary
interest are real-valued and symmetric.

Given the full-rank symmetric matrix of probabilities $B\in(0,1)^{k\times k}$,
write its spectral decomposition as

$B=U_{B}\Lambda_{B}U_{B}^{{\mathrm{\scriptscriptstyle
T}}}=U_{B}\cdot\left(|\Lambda_{B}|^{1/2}\cdot
J\cdot|\Lambda_{B}|^{1/2}\right)\cdot U_{B}^{{\mathrm{\scriptscriptstyle
T}}},$

where the $k\times k$ orthogonal matrix $U_{B}$ consists of columns that are
unit-norm eigenvectors of $B$, the $k\times k$ diagonal matrix $\Lambda_{B}$
contains the eigenvalues of $B$ arranged in decreasing order, $|\cdot|$
denotes the entrywise absolute value operation, and $J$ is the $k\times k$
diagonal matrix $J=\operatorname{diag}(1_{p},-1_{q})$ when $B$ has $p$
positive eigenvalues and $q$ negative eigenvalues, hence $p+q=k$. The matrix
$U_{B}$ is not uniquely defined, modulo column-wise signflips and up to
subspace basis non-uniqueness if $B$ has repeated eigenvalues, but this non-
uniqueness poses no appreciable problem throughout the proof due to
appropriate invariances, the consideration of certain matrix products, and the
bounding of subspace distances.

By hypothesis, the sparse latent factor matrix $Z\in\\{0,1\\}^{n\times k}$
consists of independent, identically distributed $k$-dimensional random latent
factor vectors $Z_{i}\sim\mathcal{G}=\operatorname{Multinomial}(1;\pi)$ whence
$Z_{i}\in\\{e_{1},\dots,e_{k}\\}$ and where $\pi\in\mathbb{R}^{k}$ is a fixed
probability vector in the interior of the unit simplex. Now, define the matrix
of coefficients $T_{B}=|\Lambda_{B}|^{1/2}\cdot
U_{B}^{{\mathrm{\scriptscriptstyle T}}}$, and for each index $i$ define the
$k$-dimensional random vector $\xi_{i}=T_{B}\cdot Z_{i}$. Write
$X=[\xi_{1},\dots,\xi_{n}]^{{\mathrm{\scriptscriptstyle
T}}}\in\mathbb{R}^{n\times k}$ as the matrix whose rows are the vectors
$\xi_{i}$, hence $ZBZ^{{\mathrm{\scriptscriptstyle
T}}}=XJX^{{\mathrm{\scriptscriptstyle T}}}$ almost surely. At the population
level, $\operatorname{E}(A\mid Z)=\operatorname{E}(A\mid X)$.

For semi-sparse undirected stochastic blockmodel graphs satisfying
$n\rho=\omega(\log^{c}n)$ with positive-definite connectivity matrix $B$, it
holds as in (Cape et al., 2019a, Theorem 3) that for any fixed index $i$, the
sequence of random vectors of the form

$n\rho^{1/2}\left\\{W_{Z}^{{\mathrm{\scriptscriptstyle
T}}}\cdot\left(W_{U}\widehat{U}_{i}-U_{i}\right)\mid Z_{i}=e_{\ell}\right\\}$

converges in distribution to a mean zero multivariate Gaussian random vector
as $n\rightarrow\infty$, where (i) the explicit asymptotic covariance matrix
depends on the block membership of the $i$-th node, indicated by $e_{\ell}$,
and (ii) the orthogonal matrices $W_{Z}$ and $W_{U}$, implicitly indexed by
$n$, form sequences denoted by $(W_{Z}),(W_{U})\subset\mathcal{O}(k)$. Here,
the column vector $\widehat{U}_{i}$ is the $i$-th row of the matrix
$\widehat{U}$ consisting of the top-$k$ unit norm eigenvectors of $A$.
Similarly, the column vector $U_{i}$ is the $i$-th row of the matrix $U$
consisting of the top-$k$ unit norm eigenvectors of $\operatorname{E}(A\mid
Z)$. This existing result does not involve varimax factor rotations.

In this proof, we shall first generalize the above asymptotic normality result
to hold for all full-rank stochastic blockmodels, not only those having
positive-definite connectivity matrices $B$. In other words, we shall first
extend the existing result (Cape et al., 2019a, Theorem 3); the minor
distinction between fixed versus generated block memberships is immaterial
here. Our extension will be achieved by revisiting and rewriting matrix
factorizations involving $U$ and $\Lambda^{-1}$ as equivalent, properly
oriented matrix factorizations involving $X$ and $J$. To this end, it will be
convenient to work with the reparametrization and minor change of notation
established at the start of this proof.

Let $P$ denote the matrix $\operatorname{E}(A\mid Z)$, equivalently
$\operatorname{E}(A\mid X)$, which in terms of $X$ is written as $P=\rho
XJX^{{\mathrm{\scriptscriptstyle T}}}$ and which has rank $k$ almost surely.
Denote the skinny spectral decomposition of $P$ by $U\Lambda
U^{{\mathrm{\scriptscriptstyle T}}}$, whence it follows that
$P^{\dagger}=U\Lambda^{-1}U^{{\mathrm{\scriptscriptstyle T}}}$ is the
necessarily unique Moore–Penrose inverse of $P$. Consequently, by standard
properties of the Moore–Penrose inverse,

$U\Lambda^{-1}U^{{\mathrm{\scriptscriptstyle T}}}=P^{\dagger}=(\rho
XJX^{{\mathrm{\scriptscriptstyle
T}}})^{\dagger}=\rho^{-1}(X^{{\mathrm{\scriptscriptstyle
T}}})^{\dagger}J^{\dagger}X^{\dagger}=\rho^{-1}X(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}J(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}X^{{\mathrm{\scriptscriptstyle T}}}.$

Squaring the matrix $U\Lambda^{-1}U^{{\mathrm{\scriptscriptstyle T}}}$ yields

$U\Lambda^{-2}U^{{\mathrm{\scriptscriptstyle
T}}}=\rho^{-2}X(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}J(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}J(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}X^{{\mathrm{\scriptscriptstyle T}}},$

which is a positive-semidefinite matrix, hence there exists a $k\times k$
orthogonal matrix $W_{X}$ such that

$U\Lambda^{-1}W_{X}=\rho^{-1}X(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}J(X^{{\mathrm{\scriptscriptstyle T}}}X)^{-1/2}.$

Let $E=A-P$, hence it follows that

$EU\Lambda^{-1}W_{U}=E\left(U\Lambda^{-1}W_{X}\right)W_{X}^{{\mathrm{\scriptscriptstyle
T}}}W_{U}=\rho^{-1}EX(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}J(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1/2}W_{X}^{{\mathrm{\scriptscriptstyle T}}}W_{U}.$

By taking the above display equation and proceeding via the same perturbation
analysis approach as in (Cape et al., 2019a, Supplement), we obtain the key
matrix identity

$\left(\widehat{U}W_{U}^{{\mathrm{\scriptscriptstyle
T}}}-U\right)W_{X}=\rho^{-1}EX(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}J(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1/2}+E^{\prime}W_{U}^{{\mathrm{\scriptscriptstyle T}}}W_{X},$

where $E^{\prime}\in\mathbb{R}^{n\times k}$ denotes a residual matrix of
smaller order in maximum Euclidean row norm. Namely, the $i$-th row vector in
the above display equation, when written as a column vector and scaled by
$n\rho^{1/2}$, satisfies

$n\rho^{1/2}W_{X}^{{\mathrm{\scriptscriptstyle
T}}}\left(W_{U}\widehat{U}_{i}-U_{i}\right)=(n^{-1}X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1/2}J(n^{-1}X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1}\left\\{(n\rho)^{-1/2}(EX)_{i}\right\\}+\operatorname{o}_{\operatorname{P}}(1),$
(equation)

where the notation $\chi_{1}=\chi_{2}+\operatorname{o}_{\operatorname{P}}(1)$
conveys
$\|\chi_{1}-\chi_{2}\|_{\ell_{2}}=\operatorname{o}_{\operatorname{P}}(1)$ as
$n\rightarrow\infty$. Hence, by an application of the multivariate central
limit theorem and Slutsky’s theorem, the sequence of random vectors of the
form

$n\rho^{1/2}\left\\{W_{X}^{{\mathrm{\scriptscriptstyle
T}}}\cdot\left(W_{U}\widehat{U}_{i}-U_{i}\right)\mid Z_{i}=e_{\ell}\right\\}$
(equation)

converges in distribution to a multivariate normal random vector with mean
zero and covariance matrix
$\Delta_{X}^{-1/2}J\Delta_{X}^{-1}\cdot\operatorname{E}_{\xi}\\{g_{\rho_{\infty}}(\xi,T_{B}e_{\ell})\\}\cdot\Delta_{X}^{-1}J\Delta_{X}^{-1/2}$
as $n\rightarrow\infty$. Here,
$\Delta_{X}=\operatorname{E}(\xi\xi^{{\mathrm{\scriptscriptstyle
T}}})\in\mathbb{R}^{k\times k}$ where $\xi\sim T_{B}\mathcal{G}$, while
$g_{\rho_{\infty}}(\xi,x)=(x^{{\mathrm{\scriptscriptstyle
T}}}J\xi)\cdot(1-\rho_{\infty}x^{{\mathrm{\scriptscriptstyle
T}}}J\xi)\cdot\xi\xi^{{\mathrm{\scriptscriptstyle T}}}$. Recall that by
assumption, $n\rho=\omega(\log^{c}n)$ with $\rho=\rho_{n}\in(0,1]$ and
$\rho_{n}\rightarrow\rho_{\infty}\in\\{0,1\\}$. In words, the asymptotic
covariance matrix differs between the dense graph regime $\rho=1$ and the
semi-sparse graph regime $\rho\rightarrow 0$, exhibiting larger asymptotic
variability in the latter case.

The next step of the proof is to synthesize the above perturbation analysis
together with the effect of varimax factor rotations. Our approach will make
use of existing strong consistency results, namely $\ell_{2,\infty}$ uniform
row-wise error bounds. We shall also need to demonstrate suitable orthogonal
alignment of certain matrix terms appearing in perturbation expressions.

Towards this end, first observe that $UW_{X}=(U\Lambda
U^{{\mathrm{\scriptscriptstyle
T}}})\cdot(U\Lambda^{-1}W_{X})=X(X^{{\mathrm{\scriptscriptstyle
T}}}X)^{-1/2}$, hence for each fixed index $i$ it holds almost surely that

$W_{X}^{{\mathrm{\scriptscriptstyle T}}}\cdot
n^{1/2}U_{i}=(n^{-1}X^{{\mathrm{\scriptscriptstyle T}}}X)^{-1/2}X_{i}\\\
=(n^{-1}X^{{\mathrm{\scriptscriptstyle T}}}X)^{-1/2}\cdot T_{B}\cdot Z_{i}.$

At the same time, $W_{U}$ denotes the closest $k$-dimensional orthogonal
matrix approximation to $U^{{\mathrm{\scriptscriptstyle T}}}\widehat{U}$ which
satisfies the operator norm bound $\|U^{{\mathrm{\scriptscriptstyle
T}}}\widehat{U}-W_{U}\|_{\operatorname{op}}=\operatorname{O}_{\operatorname{P}}\left\\{(n\rho)^{-1}\right\\}$
per (Cape et al., 2019a, Eq. (S2) in Supplement). Furthermore, applying (Cape
et al., 2019a, Theorem 1) yields the maximum Euclidean row norm bound
$\|\widehat{U}-UW_{U}\|_{2,\infty}=\operatorname{o}_{\operatorname{P}}\left(n^{-1/2}\right)$.
Consequently,
$n^{1/2}\|\widehat{U}\|_{2,\infty}=\operatorname{O}_{\operatorname{P}}(1)$
since $n^{1/2}\|U\|_{2,\infty}$ is bounded with overwhelming probability.
Thus, since
$\widehat{Z}_{i}=n^{1/2}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\widehat{U}_{i}$ and $X=ZT_{B}^{{\mathrm{\scriptscriptstyle T}}}$, it
follows that

$\displaystyle W_{X}^{{\mathrm{\scriptscriptstyle T}}}W_{U}\cdot
n^{1/2}\widehat{U}_{i}$ $\displaystyle\hskip
10.00002pt=W_{X}^{{\mathrm{\scriptscriptstyle
T}}}\left(U^{{\mathrm{\scriptscriptstyle T}}}\widehat{U}\right)\cdot
n^{1/2}\widehat{U}_{i}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip 10.00002pt=\left(X^{{\mathrm{\scriptscriptstyle
T}}}X\right)^{-1/2}X^{{\mathrm{\scriptscriptstyle
T}}}\cdot\widehat{U}\left(R_{\widehat{U}}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot
n^{1/2}\widehat{U}_{i}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip 10.00002pt=\left(n^{-1}X^{{\mathrm{\scriptscriptstyle
T}}}X\right)^{-1/2}\cdot T_{B}\cdot\left(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\widehat{Z}\right)\cdot\widehat{Z}_{i}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}.$

Given the asymptotic normality of the eigenvector components in Eq. equation,
an application of Slutsky’s theorem together with several steps of algebra
yields that

$(n\rho)^{1/2}\left\\{\left(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\widehat{Z}\right)\cdot\widehat{Z}_{i}-Z_{i}\mid Z_{i}=e_{\ell}\right\\}$

is asymptotically normal with mean zero. In particular, necessarily
$(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\widehat{Z})\cdot\widehat{Z}_{i}\rightarrow Z_{i}$ in probability given
$Z_{i}=e_{\ell}$, since otherwise we obtain a contradiction to the established
convergence in distribution to a multivariate normal random vector.

For any choice of $\Pi_{Z}\in\mathcal{P}(k)$, adding and subtracting
$Z\operatorname{diag}(\pi)^{-1/2}\Pi_{Z}$ below gives

$\left(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\widehat{Z}\right)=\left(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}Z\right)\operatorname{diag}(\pi)^{-1/2}\Pi_{Z}+n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\left(\widehat{Z}-Z\operatorname{diag}(\pi)^{-1/2}\Pi_{Z}\right).$
(equation)

The law of large numbers guarantees that
$\left(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}Z\right)\rightarrow\operatorname{diag}(\pi)$ almost surely as
$n\rightarrow\infty$. Furthermore, (Rohe & Zeng, 2023, Corollary C.1)
establishes that there exists a sequence of signed permutation matrices
$\left(\Pi_{Z}\right)\subset\mathcal{P}(k)$ such that

$\max_{1\leq i\leq
n}\left\|\left(\widehat{Z}-Z\operatorname{diag}(\pi)^{-1/2}\Pi_{Z}\right)_{i}\right\|_{\ell_{2}}=\operatorname{o}_{\operatorname{P}}(1).$

For this same sequence of signed permutation matrices, properties of the
$\ell_{2,\infty}$ norm yield

$\left\|n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\left(\widehat{Z}-Z\operatorname{diag}(\pi)^{-1/2}\Pi_{Z}\right)\right\|_{2,\infty}\leq\left\|n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\right\|_{\infty,\infty}\cdot\left\|\widehat{Z}-Z\operatorname{diag}(\pi)^{-1/2}\Pi_{Z}\right\|_{2,\infty},$

hence,

$\left\|n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\left(\widehat{Z}-Z\operatorname{diag}(\pi)^{-1/2}\Pi_{Z}\right)\right\|_{2,\infty}=\operatorname{o}_{\operatorname{P}}(1),$

where $\|\cdot\|_{2,\infty}$ again denotes the maximum Euclidean row norm of a
matrix, where $\|\cdot\|_{\infty,\infty}$ denotes the maximum absolute row sum
of a matrix, and where we have used the fact that
$\left\|n^{-1}Z^{{\mathrm{\scriptscriptstyle T}}}\right\|_{\infty,\infty}\leq
1$ almost surely. Moreover,
$\Pi_{Z}\widehat{Z}_{i}\rightarrow\operatorname{diag}(\pi)^{-1/2}Z_{i}$ in
probability, recalling that
$\operatorname{diag}(\pi)^{-1/2}Z_{i}=\breve{Z}_{i}$ with $Z_{i}=e_{\ell}$.
Taken together and by writing
$Z_{i}=\operatorname{diag}(\pi)^{1/2}\operatorname{diag}(\pi)^{-1/2}Z_{i}$,
these deductions imply that

$(n\rho)^{1/2}\left\\{\Pi_{Z}\cdot\widehat{Z}_{i}-\breve{Z}_{i}\mid
Z_{i}=e_{\ell}\right\\}$

is asymptotically multivariate normal as $n\rightarrow\infty$. By tracing
through the proof and again appealing to Slutsky’s theorem, the asymptotic
covariance matrix is given by

$\left\\{\operatorname{diag}(\pi)^{-1/2}\cdot T_{B}^{-1}\cdot
J\cdot\Delta_{X}^{-1}\right\\}\cdot\operatorname{E}_{\xi}\\{g_{\rho_{\infty}}(\xi,T_{B}e_{\ell})\\}\cdot\left\\{\Delta_{X}^{-1}\cdot
J\cdot(T_{B}^{-1})^{{\mathrm{\scriptscriptstyle
T}}}\cdot\operatorname{diag}(\pi)^{-1/2}\right\\}.$

To reiterate, here $\pi$ is the block membership probability vector for the
latent factor vectors, here
$T_{B}=|\Lambda_{B}|^{1/2}U_{B}^{{\mathrm{\scriptscriptstyle T}}}$ is defined
such that $B=T_{B}^{{\mathrm{\scriptscriptstyle T}}}JT_{B}$ where
$J=\operatorname{diag}(1_{p},-1_{q})\in\mathbb{R}^{k\times k}$, here
$\Delta_{X}=\operatorname{E}(\xi\xi^{{\mathrm{\scriptscriptstyle T}}})$ where
$\xi$ has distribution given by $T_{B}\mathcal{G}$, and
$g_{\rho_{\infty}}(\xi,x)=(x^{{\mathrm{\scriptscriptstyle
T}}}J\xi)\cdot(1-\rho_{\infty}x^{{\mathrm{\scriptscriptstyle
T}}}J\xi)\cdot\xi\xi^{{\mathrm{\scriptscriptstyle T}}}$ where
$\rho_{\infty}\in\\{0,1\\}$. This concludes the proof.

###  section.subsection Asymptotic normality in varimax-transformed directed
graph embeddings

######  Proof section.theorem (of Theorem 2).

By hypothesis, consider directed stochastic blockmodels having edge
probabilities given by the matrix $\operatorname{E}(A\mid Z,Y)=\rho
ZBY^{{\mathrm{\scriptscriptstyle T}}}$, consisting of independent, identically
distributed left latent factor vectors
$Z_{i}\sim\mathcal{G}=\operatorname{Multinomial}(1;\pi_{Z})$ and independent,
identically distributed right latent factor vectors
$Y_{i}\sim\mathcal{H}=\operatorname{Multinomial}(1;\pi_{Y})$ and where $Z_{i}$
is independent of $Y_{i}$ for all $1\leq i\leq n$. Write the singular value
decomposition of the full-rank $k\times k$ matrix $B$ as
$B=U_{B}S_{B}V_{B}^{{\mathrm{\scriptscriptstyle T}}}$. Define the random
vectors $\xi_{i}=T_{\Xi}\cdot
Z_{i}=S_{B}^{1/2}U_{B}^{{\mathrm{\scriptscriptstyle T}}}\cdot Z_{i}$ and
$\upsilon_{i}=T_{\Upsilon}\cdot
Y_{i}=S_{B}^{1/2}V_{B}^{{\mathrm{\scriptscriptstyle T}}}\cdot Y_{i}$ to form
the rows of $\Xi$ and $\Upsilon$, whence $ZBY^{{\mathrm{\scriptscriptstyle
T}}}=\Xi\Upsilon^{{\mathrm{\scriptscriptstyle T}}}$ almost surely.
Consequently, $\operatorname{E}(A\mid
Z,Y)=\operatorname{E}(A\mid\Xi,\Upsilon)$ which for convenience shall be
denoted simply by $P$.

By extending the proof analysis found in Section section.subsection and in the
supplementary material of Cape et al. (2019a) to the case of directed graphs
and their truncated singular value decompositions, we obtain the key
expression

$n\rho^{1/2}\cdot W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}\left(W_{V}\widehat{U}_{i}-W_{V}W_{U}^{{\mathrm{\scriptscriptstyle
T}}}U_{i}\right)=(n^{-1}\Xi^{{\mathrm{\scriptscriptstyle
T}}}\Xi)^{-1/2}(n^{-1}\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1}\left\\{(n\rho)^{-1/2}(E\Upsilon)_{i}\right\\}+\operatorname{o}_{\operatorname{P}}(1)$
(equation)

which generalizes Eq. equation. Eq. equation can alternatively be established
with the help of Theorem 9 in Yan et al. (2021). In the above display
equation, $W_{U}$ minimizes $\|\widehat{U}-UW\|_{\operatorname{F}}$ over all
$k\times k$ orthogonal matrices $W$. Separately, $W_{V}$ minimizes
$\|\widehat{V}-VW\|_{\operatorname{F}}$ over all $k\times k$ orthogonal
matrices $W$. Here, $W_{\Upsilon}$ denotes the orthogonal matrix satisfying

$VS^{-1}W_{\Upsilon}=\rho^{-1}(\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1}(\Xi^{{\mathrm{\scriptscriptstyle T}}}\Xi)^{-1/2},$

hence $UW_{\Upsilon}=(USV^{{\mathrm{\scriptscriptstyle
T}}})\cdot(VS^{-1}W_{\Upsilon})=\Xi(\Xi^{{\mathrm{\scriptscriptstyle
T}}}\Xi)^{-1/2}$. By analogy, for the orthogonal matrix $W_{\Xi}$, it holds
that $US^{-1}W_{\Xi}=\rho^{-1}(\Xi^{{\mathrm{\scriptscriptstyle
T}}}\Xi)^{-1}(\Upsilon^{{\mathrm{\scriptscriptstyle T}}}\Upsilon)^{-1/2}$ and
so $VW_{\Xi}=\Upsilon(\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1/2}$.

A routine computation shows $\|U^{{\mathrm{\scriptscriptstyle
T}}}\widehat{U}-W_{U}\|_{\operatorname{op}}\leq\|\widehat{U}\widehat{U}^{{\mathrm{\scriptscriptstyle
T}}}-UU^{{\mathrm{\scriptscriptstyle
T}}}\|_{\operatorname{op}}^{2}=\operatorname{O}_{\operatorname{P}}\left\\{(n\rho)^{-1}\right\\}$
and similarly $\|V^{{\mathrm{\scriptscriptstyle
T}}}\widehat{V}-W_{V}\|_{\operatorname{op}}=\operatorname{O}_{\operatorname{P}}\left\\{(n\rho)^{-1}\right\\}$,
where the stated bounds hold by an application of Wedin’s theorem, see for
example (Chen et al., 2021, Theorem 2.9).

Next, consider the expansion

$\displaystyle W_{\Upsilon}^{{\mathrm{\scriptscriptstyle T}}}W_{V}$
$\displaystyle\hskip 10.00002pt=W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}V^{{\mathrm{\scriptscriptstyle
T}}}\widehat{V}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip 10.00002pt=W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}W_{\Xi}W_{\Xi}^{{\mathrm{\scriptscriptstyle
T}}}V^{{\mathrm{\scriptscriptstyle
T}}}\widehat{V}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip 10.00002pt=W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}W_{\Xi}(\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1/2}\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\widehat{V}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip 10.00002pt=W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}W_{\Xi}\cdot(n^{-1}\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1/2}\cdot T_{\Upsilon}\cdot
n^{-1}\left(Y^{{\mathrm{\scriptscriptstyle T}}}\cdot
n^{1/2}\widehat{V}R_{\widehat{V}}\right)R_{\widehat{V}}^{{\mathrm{\scriptscriptstyle
T}}}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip 10.00002pt=W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}W_{\Xi}\cdot(n^{-1}\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1/2}\cdot
T_{\Upsilon}\cdot\left(n^{-1}Y^{{\mathrm{\scriptscriptstyle
T}}}\widehat{Y}\right)R_{\widehat{V}}^{{\mathrm{\scriptscriptstyle
T}}}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}.$

Crucially, on the right-hand side, the leftmost matrix product can be expanded
as

$\displaystyle W_{\Upsilon}^{{\mathrm{\scriptscriptstyle T}}}W_{\Xi}$
$\displaystyle\hskip 10.00002pt=W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}U^{{\mathrm{\scriptscriptstyle T}}}UV^{{\mathrm{\scriptscriptstyle
T}}}VW_{\Xi}$ $\displaystyle\hskip
10.00002pt=(\Xi^{{\mathrm{\scriptscriptstyle
T}}}\Xi)^{-1/2}\Xi^{{\mathrm{\scriptscriptstyle
T}}}\cdot\left(UV^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot\Upsilon(\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1/2}$ $\displaystyle\hskip
10.00002pt=(\Xi^{{\mathrm{\scriptscriptstyle
T}}}\Xi)^{-1/2}\cdot(T_{\Xi}Z^{{\mathrm{\scriptscriptstyle
T}}})\cdot\left(UV^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot(YT_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}})\cdot(\Upsilon^{{\mathrm{\scriptscriptstyle T}}}\Upsilon)^{-1/2}$
$\displaystyle\hskip 10.00002pt=(\Xi^{{\mathrm{\scriptscriptstyle
T}}}\Xi)^{-1/2}\cdot(T_{\Xi}Z^{{\mathrm{\scriptscriptstyle
T}}})\cdot\left(n^{-1}Z(Z^{{\mathrm{\scriptscriptstyle
T}}}Z)^{-1/2}\widetilde{R}_{U}^{{\mathrm{\scriptscriptstyle
T}}}\widetilde{R}_{V}(Y^{{\mathrm{\scriptscriptstyle
T}}}Y)^{-1/2}Y^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot(YT_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}})\cdot(\Upsilon^{{\mathrm{\scriptscriptstyle T}}}\Upsilon)^{-1/2},$

where the orthogonal matrix $\widetilde{R}_{U}^{{\mathrm{\scriptscriptstyle
T}}}$ consists of the left singular vectors of the matrix
$(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}Z)^{1/2}B(n^{-1}Y^{{\mathrm{\scriptscriptstyle T}}}Y)^{1/2}$ and
$\widetilde{R}_{V}$ consists of the corresponding right singular vectors.

The matrix product $(n^{-1}Y^{{\mathrm{\scriptscriptstyle
T}}}Y)^{-1/2}(n^{-1}Y^{{\mathrm{\scriptscriptstyle
T}}}Y)T_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}(n^{-1}\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1/2}(n^{-1}\Upsilon^{{\mathrm{\scriptscriptstyle
T}}}\Upsilon)^{-1/2}T_{\Upsilon}=(n^{-1}Y^{{\mathrm{\scriptscriptstyle
T}}}Y)^{-1/2}$ converges to $\operatorname{diag}(\pi_{Y})^{-1/2}$ almost
surely as $n\rightarrow\infty$ by the law of large numbers. Moreover, by the
same proof strategy used to analyze Eq. equation, there exists a sequence of
signed permutation matrices $(\Pi_{Y})\subset\mathcal{P}(k)$ indexed by $n$
such that $(n^{-1}Y^{{\mathrm{\scriptscriptstyle
T}}}\widehat{Y})\Pi_{Y}^{{\mathrm{\scriptscriptstyle T}}}$ converges to
$\operatorname{diag}(\pi_{Y})^{1/2}$. By putting the pieces together and with
a slight abuse of notation,

$\displaystyle W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}W_{V}\cdot(n^{1/2}\widehat{U}_{i})$ $\displaystyle\hskip
20.00003pt=\left\\{\Delta_{\Xi}^{-1/2}T_{\Xi}\operatorname{diag}(\pi_{Z})^{1/2}\right\\}\cdot\left(\widetilde{R}_{U}^{{\mathrm{\scriptscriptstyle
T}}}\widetilde{R}_{V}\Pi_{Y}R_{\widehat{V}}^{{\mathrm{\scriptscriptstyle
T}}}R_{\widehat{U}}\right)\cdot\left(n^{1/2}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\widehat{U}_{i}\right)+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}.$

Above, we are using the fact that
$\|\widehat{U}-UW_{U}\|_{2,\infty}=\operatorname{o}_{\operatorname{P}}\left(n^{-1/2}\right)$
and consequently
$n^{1/2}\|\widehat{U}\|_{2,\infty}=\operatorname{O}_{\operatorname{P}}(1)$
since $n^{1/2}\|U\|_{2,\infty}$ is bounded with overwhelming probability, as
in the setting of undirected stochastic blockmodels. Similarly,
$\|\widehat{V}-VW_{V}\|_{2,\infty}=\operatorname{o}_{\operatorname{P}}\left(n^{-1/2}\right)$
and
$n^{1/2}\|\widehat{V}\|_{2,\infty}=\operatorname{O}_{\operatorname{P}}(1)$.

Next, observe that for any choice of sequence
$(\Pi_{Z})\subset\mathcal{P}(k)$, for each value of $n$,

$\left(\widetilde{R}_{U}^{{\mathrm{\scriptscriptstyle
T}}}\widetilde{R}_{V}\Pi_{Y}R_{\widehat{V}}^{{\mathrm{\scriptscriptstyle
T}}}R_{\widehat{U}}\right)\Pi_{Z}^{{\mathrm{\scriptscriptstyle
T}}}=\Pi_{Z}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\cdot\left(R_{\widehat{U}}\Pi_{Z}^{{\mathrm{\scriptscriptstyle
T}}}\widetilde{R}_{U}^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot\left(R_{\widehat{V}}\Pi_{Y}^{{\mathrm{\scriptscriptstyle
T}}}\widetilde{R}_{V}^{{\mathrm{\scriptscriptstyle
T}}}\right)^{{\mathrm{\scriptscriptstyle T}}}\cdot
R_{\widehat{U}}\Pi_{Z}^{{\mathrm{\scriptscriptstyle T}}},$

hence choosing the sequences $(\Pi_{Z})$ and $(\Pi_{Y})$ per Rohe & Zeng
(2023) yields that the matrix product in the preceding display equation
converges in probability to the identity matrix in the large-$n$ limit.

In a similar fashion, and with slight abuse of notation, it also holds that

$\displaystyle W_{U}^{{\mathrm{\scriptscriptstyle
T}}}\cdot\left(n^{1/2}U_{i}\right)$ $\displaystyle\hskip
10.00002pt=\left(R_{\widehat{U}}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot\widehat{U}^{{\mathrm{\scriptscriptstyle
T}}}U\cdot\left(n^{1/2}U_{i}\right)+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip
10.00002pt=R_{\widehat{U}}\left(n^{1/2}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\widehat{U}^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot\left(UW_{\Upsilon}W_{\Upsilon}^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot
U_{i}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip
10.00002pt=R_{\widehat{U}}\left(n^{1/2}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\widehat{U}^{{\mathrm{\scriptscriptstyle
T}}}\right)\cdot\Xi\left(\Xi^{{\mathrm{\scriptscriptstyle
T}}}\Xi\right)^{-1}\cdot\Xi_{i}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip
10.00002pt=R_{\widehat{U}}\left(n^{1/2}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\widehat{U}^{{\mathrm{\scriptscriptstyle T}}}\right)\cdot
ZT_{\Xi}^{{\mathrm{\scriptscriptstyle
T}}}\left(T_{\Xi}Z^{{\mathrm{\scriptscriptstyle
T}}}ZT_{\Xi}^{{\mathrm{\scriptscriptstyle T}}}\right)^{-1}\cdot
T_{\Xi}Z_{i}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip
10.00002pt=R_{\widehat{U}}\left(n^{1/2}R_{\widehat{U}}^{{\mathrm{\scriptscriptstyle
T}}}\widehat{U}^{{\mathrm{\scriptscriptstyle T}}}\right)\cdot
Z\left(Z^{{\mathrm{\scriptscriptstyle T}}}Z\right)^{-1}\cdot
Z_{i}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip
10.00002pt=R_{\widehat{U}}\left(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\widehat{Z}\right)^{{\mathrm{\scriptscriptstyle
T}}}\cdot\left(n^{-1}Z^{{\mathrm{\scriptscriptstyle T}}}Z\right)^{-1}\cdot
Z_{i}+\operatorname{o}_{\operatorname{P}}\left\\{(n\rho)^{-1/2}\right\\}$
$\displaystyle\hskip
10.00002pt=R_{\widehat{U}}\Pi_{Z}^{{\mathrm{\scriptscriptstyle
T}}}\cdot\breve{Z}_{i}+\operatorname{o}_{\operatorname{P}}(1),$

where $\Pi_{Z}\cdot\left(n^{-1}Z^{{\mathrm{\scriptscriptstyle
T}}}\widehat{Z}\right)^{{\mathrm{\scriptscriptstyle T}}}\cdot Z_{i}$ converges
to $\operatorname{diag}(\pi)^{1/2}Z_{i}$ given $Z_{i}=e_{\ell}$ by an
analogous argument as in the preceding proof, leveraging the asymptotic
normality arising from Eq. equation. Consequently, there exists a sequence of
signed permutation matrices $(\Pi_{Z})\subset\mathcal{P}(k)$ such that for
each fixed index $i$, the sequence of random vectors

$(n\rho)^{1/2}\left\\{\Pi_{Z}\cdot\widehat{Z}_{i}-\breve{Z}_{i}\mid
Z_{i}=e_{\ell}\right\\}$

converges in distribution to a multivariate Gaussian random vector with mean
zero and covariance matrix given by

$\left\\{\operatorname{diag}(\pi_{Z})^{-1/2}\cdot
T_{\Xi}^{-1}\cdot\Delta_{\Upsilon}^{-1}\right\\}\cdot\operatorname{E}_{\upsilon}\left\\{g_{\rho_{\infty}}(\upsilon,T_{\Xi}e_{\ell})\right\\}\cdot\left\\{\Delta_{\Upsilon}^{-1}\cdot(T_{\Xi}^{-1})^{{\mathrm{\scriptscriptstyle
T}}}\cdot\operatorname{diag}(\pi_{Z})^{-1/2}\right\\}.$

By transposing the proof arguments to instead analyze $\widehat{Y}_{j}$, there
exists a sequence of signed permutation matrices
$(\Pi_{Y})\subset\mathcal{P}(k)$ such that for each fixed index $j$, the
sequence of random vectors of the form

$(n\rho)^{1/2}\left\\{\Pi_{Y}\cdot\widehat{Y}_{j}-\breve{Y}_{j}\mid
Y_{j}=e_{\ell}\right\\}$

converges in distribution to a multivariate Gaussian random vector with mean
zero and covariance matrix given by

$\left\\{\operatorname{diag}(\pi_{Y})^{-1/2}\cdot
T_{\Upsilon}^{-1}\cdot\Delta_{\Xi}^{-1}\right\\}\cdot\operatorname{E}_{\xi}\left\\{g_{\rho_{\infty}}(\xi,T_{\Upsilon}e_{\ell})\right\\}\cdot\left\\{\Delta_{\Xi}^{-1}\cdot(T_{\Upsilon}^{-1})^{{\mathrm{\scriptscriptstyle
T}}}\cdot\operatorname{diag}(\pi_{Y})^{-1/2}\right\\}.$ (equation)

To reiterate, here $\pi_{Z}$ and $\pi_{Y}$ denote the block membership
probability vectors for the left and right latent factor vectors $Z_{i}$ and
$Y_{j}$. The matrices of coefficients $T_{\Xi}$ and $T_{\Upsilon}$ defined in
terms of $B$ satisfy $B=T_{\Xi}^{{\mathrm{\scriptscriptstyle
T}}}T_{\Upsilon}$. Here,
$\Delta_{\Xi}=\operatorname{E}(\xi\xi^{{\mathrm{\scriptscriptstyle T}}})$
where $\xi\sim T_{\Xi}\mathcal{G}$ and
$\Delta_{\Upsilon}=\operatorname{E}(\upsilon\upsilon^{{\mathrm{\scriptscriptstyle
T}}})$ with $\upsilon\sim T_{\Upsilon}\mathcal{H}$, while
$g_{\rho_{\infty}}(\alpha,\beta)=(\beta^{{\mathrm{\scriptscriptstyle
T}}}\alpha)\cdot(1-\rho_{\infty}\beta^{{\mathrm{\scriptscriptstyle
T}}}\alpha)\cdot\alpha\alpha^{{\mathrm{\scriptscriptstyle T}}}$ depends on
$\rho_{\infty}\in\\{0,1\\}$ and takes $k$-dimensional input vectors
$\alpha,\beta$. This concludes the proof.

###  section.subsection Asymptotic normality in varimax-transformed degree-
corrected graph embeddings

######  Proof section.theorem (outline of Theorem 3).

The proof proceeds in the same manner as in Section section.subsection but
instead by replacing $Z_{i}$ with $\theta_{i}\cdot Z_{i}$, properly scaling,
and then finally integrating over all values
$b\in\operatorname{supp}\mathcal{F}$. The compactness assumption precludes
potential degenerate behavior near the boundary. Here, the integration step
arises due to the general inability to condition on events of the form
$\\{\theta_{i}\cdot Z_{i}=x\\}$ for fixed vectors $x$, since the distribution
$\mathcal{F}$ need not be discrete.