The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
78,452 characters
Hypothesis testing on invariant subspaces of non-diagonalizable matrices with applications to network statistics
\global\long\def\uwrite#1#2{\underset{#2}{\underbrace{#1}} }
\global\long\def\blw#1{\ensuremath{\underline{#1}}}
\global\long\def\abv#1{\ensuremath{\overline{#1}}}
\global\long\def\vect#1{\mathbf{#1}}
\global\long\def\smlseq#1{\{#1\} }
\global\long\def\seq#1{\left\{ #1\right\} }
\global\long\def\smlsetof#1#2{\{#1\mid#2\} }
\global\long\def\setof#1#2{\left\{ #1\mid#2\right\} }
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long
\global\long\global\long\global\long\global\long\global\long\global\long\def\Ellp#1{\ensuremath{\mathcal{L}^{#1}}}
\global\long\global\long\global\long\global\long\global\long
\global\long\def\abs#1{\ensuremath{\left|#1\right|}}
\global\long\def\smlabs#1{\ensuremath{\lvert#1\rvert}}
\global\long\def\bigabs#1{\ensuremath{\bigl|#1\bigr|}}
\global\long\def\Bigabs#1{\ensuremath{\Bigl|#1\Bigr|}}
\global\long\def\biggabs#1{\ensuremath{\biggl|#1\biggr|}}
\global\long\def\norm#1{\ensuremath{\left\Vert #1\right\Vert }}
\global\long\def\smlnorm#1{\ensuremath{\lVert#1\rVert}}
\global\long\def\bignorm#1{\ensuremath{\bigl\|#1\bigr\|}}
\global\long\def\Bignorm#1{\ensuremath{\Bigl\|#1\Bigr\|}}
\global\long\def\biggnorm#1{\ensuremath{\biggl\|#1\biggr\|}}
\global\long\global\long\global\long\global\long\global\long\global\long\def\clsr#1{\ensuremath{\overline{#1}}}
\global\long\global\long\global\long\global\long
\global\long\def\smlinprd#1#2{\ensuremath{\langle#1,#2\rangle}}
\global\long\def\inprd#1#2{\ensuremath{\left\langle #1,#2\right\rangle }}
\global\long\global\long
\global\long\global\long\global\long\global\long
\global\long\global\long\global\long\def\sigf#1{\mathcal{#1}}
\global\long\global\long\global\long\def\flt#1{\mathcal{#1}}
\global\long\global\long\global\long\global\long\global\long\global\long\global\long
\global\long\global\long\global\long\global\long\global\long\global\long
\global\long\global\long\global\long\global\long\global\long
\global\long\def\independenT#1#2{\mathrel{\rlap{$#1#2$}\mkern2mu {#1#2}}}
\global\long\global\long\global\long\global\long\global\long\global\long\def\inprobu#1{\ensuremath{\overset{#1}{\ensuremath{\rightarrow}}}}
\global\long\global\long\global\long\def\inLp#1{\ensuremath{\overset{\Ellp{#1}}{\ensuremath{\rightarrow}}}}
\global\long\global\long\global\long\global\long\def\wkcu#1{\overset{#1}{\ensuremath{\rightsquigarrow}}}
\global\long
\global\long\global\long\global\long\global\long\global\long\global\long\global\long
\global\long\global\long\global\long\global\long\global\long\global\long\def\cv#1{\left\langle #1\right\rangle }
\global\long\def\smlcv#1{\langle#1\rangle}
\global\long\def\qv#1{\left[#1\right]}
\global\long\def\smlqv#1{[#1]}
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\def\smlfloor#1{\lfloor#1\rfloor}
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\def\spc#1{\mathcal{#1}}
\global\long\def\set#1{\mathscr{#1}}
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\def\largedec#1{\mathbf{#1}}
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long
\global\long
\href{https://arxiv.org/pdf/2303.18233.pdf}{Please click here for the arxiv version.}
\title{Hypothesis testing on invariant subspaces of non-diagonalizable matrices
with applications to network statistics}
\author{Jérôme R. Simons}
\date{8 October `25}
\begin{abstract}
We generalise the inference procedure for eigenvectors of symmetrizable
matrices of \citet{Tyler1981} to that of invariant and singular subspaces
of non-diagonalizable matrices. Wald tests for invariant vectors and
$t$-tests for their individual coefficients perform well in simulations,
despite the matrix being not symmetric. Using these results, it is
now possible to perform inference on network statistics that depend
on eigenvectors of non-symmetric adjacency matrices as they arise
in empirical applications from directed networks. Further, we find
that statisticians only need control over the first-order Davis-Kahan
bound to control convergence rates of invariant subspace estimators
to higher-orders. For general invariant subspaces, the minimal eigenvalue
separation dominates the first-order bound potentially slowing convergence
rates considerably. In an example, we find that accounting for uncertainty
in network estimates changes empirical conclusions about the ranking
of nodes' popularity.
\thanks{An earlier version of this paper was entitled \textquotedblleft Inference
on non-symmetric subspaces and network statistics\textquotedblright}\thanks{I would like to thank James A. Duffy, Steve Bond, Michael Leung, Richard
Samworth, David E. Tyler, Eric French, Alexei Onatskiy, Oliver Linton,
Richard Smith, Andrew Harvey, Patrick Allmis, Christian Ghiglino,
Carsten-Andreas Schulz, and Sam Gee for their comments and suggestions.
I am also grateful to the organizers and participants of Encounters
in Econometric Theory, seminars at Cambridge, Oxford, and the Summer
Meeting of the Econometric Society. I acknowledge funding from the
Keynes Fund at the Faculty of Economics, Cambridge.}
\end{abstract}
\maketitle
\section{Introduction}
This paper contributes hypothesis tests for both invariant subspace
and singular vectors of non-symmetric matrices that are not diagonalisable.
As an application of our theory, we specialise tests for a selection
of centrality and clustering statistics as they arise as functions
of network adjacency matrices. While network statistics are perhaps
most empirically relevant, our results are in the form of general
$t$- and Wald tests with the latter reducing to the procedure developed
in \citet{Tyler1981}, when the matrix has a real spectrum and is
diagonalizable.
The source of randomness in the network context is uncertainty about
the extent of weights or links where errors propagate to these statistics.
Allowing for non-symmetric adjacency matrices opens up many empirically
relevant applications. Directed, weighted networks for example arise
when weights depend on the flow direction between nodes. Specifically,
trade, input-output, and food chain networks trigger directed graphs
where direction matters. In this context, we assume a researcher has
a network adjacency matrix estimator at hand. For example, social
interaction models such as those described in \citet{depaula2023identifying,rothenhausler2015backshift,manresa}
treat adjacency matrix entries as estimands. Another variant is the
sampling of graphons that leads to noisy network matrices, developed
among others by \citet{10.1093/biomet/asac032,parise2023graphon}.
On the basis of such models, our results let us construct standard
errors for derived network statistics so that researchers can quantify
the propagated uncertainty.
In an application, we examine how confidence intervals for network
centralities arising in a simple network model provide a cautionary
tale about ranking nodes' popularity: reordering based on the upper
ends of the confidence intervals reorders the nodes' popularity in
one example but leaves the ordering undisturbed in another.
Besides, we also offer Monte Carlo evidence for the quality of the
distributional approximations. They perform well, but do depend on
the quality of the underlying matrix estimator. We also study the
performance of the $t$-test for a data-generating process that starts
with a random graph model, which experiences normally distributed
disturbances.
Beyond networks, there are many statistical applications that require
researchers to find eigenvectors of matrices estimated with error.
For example, companion matrices of vector auto-regressions are not
symmetric yet their spectrum carries information about the dynamics.
Eigenvectors associated with unit eigenvalues of these matrices identify
cointegrating relations, for which our inference methods are also
useful. Similarly, eigenvectors are used to estimate functional diversity
in ecology.
We calculate convergence rates of subspace-based estimators whenever
the convergence speed in the form of the Frobenius norm $\smlnorm{\hat{M}-M}_{\text{F}}$
is known. We also approximate higher-order bounds and learn that for
invariant subspaces, the eigenvalue separation dominates all higher-order
terms. To control convergence rates of eigenspaces, statisticians
only need control over the first-order bound as all higher-order bounds
are powers of the former. To first-order, we recover the version of
the \citet{davis1969some} bound found in \citet{demetrius}. These
results are helpful to strengthen consistency results to explicit
convergence rates.
Invariant subspaces of adjacency matrices also appear in latent space
graph models, where the latent space is either an invariant subspace
directly of the adjacency matrix or graph Laplacian or a hidden Euclidean
space that can be estimated via invariant subspaces as shown in \citet{ZhangXuZhu2022}.
Invariant subspaces are also used to approximate the latent spaces
in random dot product graph models \citep{10.1007/978-3-540-77004-6_11}
where \citet{XieXu2020} propose a method to estimate these spaces
and \citet{AthreyaEtAl2018} outline how to perform inference in random
dot product graphs.
An important distinction arises in the dimensionality of the $p\times p$
matrix $M$. In many scientific disciplines, network size $p$ is
modest and multiple observations of $M$ are available such as network
detection and modelling \citep{cattuto2010dynamics,krivitsky2014separable,prawesh2019small}
or health dynamics \citep{rothenberg1998social,cornwell2009network,christakis2010social},
or biology \citep{krause2009social,isella2011sociopatterns}. This
paper is predominantly about this case, which arises e.g. from the
measurement of many small networks, sampled over time. Importantly,
none of the Jacobians, the convergence results, and perturbation estimations
are sensitive to the size of $M$ whereas the hypothesis tests in
\prettyref{sec:hyp-tests} use fixed size CLTs and require estimation
of a $p^{2}\times p^{2}$ covariance matrix by standard methods.
Intuitively, symmetric matrices are convenient because small perturbations
to the entries correspond to small perturbations to the eigenvectors
and the eigenvector map is smooth. As one moves away from symmetric
matrices, eigenvalues acquire imaginary parts and come in complex
conjugate pairs. Therefore, eigenvalues may lie closely together though
differentiability still holds. Finally, if multiple eigenvalues correspond
to a single eigenvector, the eigenvector map is non-differentiable.
In this case, we can only distinguish groups of eigenvalues. Our methods
specialise in this latter case.
\section{Setup}
\label{sec:setup}
\subsection{General framework for subspace inference\label{subsec:General-framework}}
We define the framework for the invariant and singular subspace inference
problems. For a $p\times p$ matrix $M$, the columns of the $p\times q$
matrix $R$ span a right-invariant subspace of dimension $q$ iff
the relation
\begin{equation}
MR=R\Lambda\label{eq:inv-subs-rel}
\end{equation}
holds for some not necessarily full-rank $q\times q$ matrix $\Lambda$
with the analogous relationship $L^{\mathsf{T}}M=\Lambda L^{\mathsf{T}}$
and $R^{-1}=L^{\mathsf{T}}$. In this case, the columns of $R$ span
a right-invariant subspace of $M$. For inference on eigenvectors,
a requirement by \citet{Tyler1981} is that there exists a positive
definite symmetric matrix $\Gamma$ such that $\Gamma M$ is symmetric,
which we relax. We illustrate this condition in
\begin{example}
\label{exa:symmetrisable}The matrix $M_{1}$ is symmetric in the
metric of $\Gamma$ for
\begin{align*}
\Gamma & =\begin{bmatrix}2 & 1\\
1 & 6
\end{bmatrix} & & M_{1}=\begin{bmatrix}1 & 3\\
1 & 1
\end{bmatrix}
\end{align*}
whereas for parameters $\lambda,a$
\begin{align*}
M_{2}=\begin{bmatrix}\lambda & a\\
0 & \lambda
\end{bmatrix}
\end{align*}
is not symmetric in the metric of $\Gamma$. A $2\times2$ matrix
$\Gamma$ is not guaranteed to be positive-definite for $\Gamma M_{2}$
to be symmetric and $M_{2}$ is non-diagonalizable and only has a
single eigenvector $\left[\begin{smallmatrix}1 & 0\end{smallmatrix}\right]^{\mathsf{T}}$.
Its single eigenvector however spans an invariant subspace. Therefore,
Tyler's procedure applies to $M_{1}$ but not to $M_{2}$.
\end{example}
Essentially, the requirement we relax stipulates that a matrix $M$
be symmetric with respect to some positive definite inner product
which implies diagonalizability and real eigenvalues.
We assume that we observe a sample of $n$ matrices $M_{1},\dots,M_{n}$,
whence we estimate the mean $\hat{M}_{n}$ so that the columns of
$\sqrt{n}(\hat{M}_{n}-M)$ converge weakly to a multivariate normal
distribution centered at zero. Denote the span of the columns of
interest in $R$ by $\mathcal{S}_{I}\left(M\right)$ with associated set of eigenvalues $\mathcal{L}_{I}$.
Then, our null hypothesis is
\begin{equation}
H_{0}:\upsilon\in\mathcal{S}_{I}\left(M\right)\label{eq:basic-hypothesis-random-matrix}
\end{equation}
against the one-sided alternative $H_{1}:\upsilon\notin\mathcal{S}_{I}(M)$.
The elements of $\mathcal{L}_{J}$ index all other directions of interest
so that if $M$ is diagonalizable, eigenvalues $\mathcal{L}$ equal $\mathcal{L}_{I}\ensuremath{\cup}\mathcal{L}_{J}$.
To attain a statistic equalling zero under $H_{0}$, we focus on the
orthocomplement $\upsilon_{\perp}$, which is an element of the left-invariant
subspace $\mathcal{S}_{J}\left(M\right)$ formed by the columns of $(R)^{-\mathsf{T}}$ so that
\begin{equation}
H_{0}^{\ast}\,:\,\upsilon_{\perp}\in\mathcal{S}_{J}\left(M\right),\label{eq:perp-inference}
\end{equation}
which is equivalent to \prettyref{eq:basic-hypothesis-random-matrix}.
\prettyref{app:projections-and-gen-inverses} provides details. Therefore,
a researcher specifies up to $m\leq q$ columns of $\upsilon$ leading
to an orthocomplement with $p-m$ columns. It is, of course, possible
to only select a subset of the orthocomplement if $p$ is large as
long as the hypothesised vector of interest is orthogonal to $\upsilon_{\perp}$.
We also consider an extension to inference on the left-singular\footnote{An eigendecomposition of $M=R\Lambda^{\mathsf{T}}$ implies that $R$
contains right eigenvectors, while for a singular value decomposition,
$M=U\Sigma V^{\mathsf{T}}$, $U$ contains left singular vectors. Our
study focusses on those matrices appearing `on the left.'} subspace of $M$ spanned by the columns of $U$ in
\[
M=U\Sigma V^{\mathsf{T}}.
\]
Suppose a researcher observes a sequence of matrices $\{M_{t}\}_{t=1}^{n}$
that constitute noisy measurements from an underlying model
\[
M_{t}=M+\varepsilon_{t},
\]
for an idiosyncratic error $\varepsilon_{t}$. The estimator $\hat{M}_{T}\coloneqq T^{-1}\sum_{t=1}^{T}M_{t}$
is consistent.
The hypothesis tests covered in the next section focus mainly on the
case where the size of $M$ is fixed although we make no formal restrictions
on it. For some of the network applications in \prettyref{sec:centralities},
the observations are rows and columns of $M$ which consequently grows
in size. To make our setup formal, we have
\begin{assumption}
\label{assu:general-ass}~Let $M\in\mathbb{R}^{p\times p}$ be a general
matrix and let $\Omega\in\mathbb{R}^{p^{2}\times p^{2}}$ denote a positive-definite
covariance matrix. Then,
\end{assumption}
\begin{enumerate}
\item \label{enu:tightness}If $p=n$, is growing, the estimator $\hat{M}$
is $1/r_{n}$ consistent (tight), i.e. $\text{\ensuremath{\smlnorm{\hat{M}_{n}-M_{n}}}}_{\text{F}}=O_{p}(r_{n})$
for the Frobenius norm.
\item \label{enu:asymptotically-normal-estimator}The model $M_{t}=M+\varepsilon_{t}$
generates data $\left\{ M_{t}\right\} _{t=1}^{T}$ where $\varepsilon_{t}$
are identical with general covariance matrix $\Omega=\ensuremath{\mathbb{E}}\text{vec\ }\varepsilon_{i}\,\text{vec\ }\varepsilon_{j}^{\mathsf{T}}$
for all $i,j=1,\dots,T$.
\item \label{enu:cov-mat-estimation}We can consistently estimate the covariance
matrix by some covariance estimator $\hat{\Omega}$ so that $\hat{\Omega}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}}\Omega$.
\end{enumerate}
\prettyref{assu:general-ass}\prettyref{enu:tightness} implies
that $\hat{M}$ converges weakly and covers the setup where a network
researcher estimates network links based on a single sample of $M$.
Further, this assumption is the foundation for invariant subspaces
of large, random matrices. In those cases, special care has to be
taken. \citet[Example 2.3]{benaychgeorges2018lectureslocalsemicirclelaw}
show that eigenvectors of length $n$ tend to be scaled by $n^{-1/2}$.
Our focus in this study is simply how the convergence rate propagates
to estimators of invariant subspaces. \citet{cai2021network} provides
details.
\prettyref{assu:general-ass}\prettyref{enu:asymptotically-normal-estimator}
implies
\[
\sqrt{T}\,\text{vec\ }\left(\hat{M}_{T}-M\right)\ensuremath{\rightsquigarrow} N\left(0,\Omega\right),
\]
which admits a general covariance structure within $\Omega$ and imposes
a stationary error distribution over the index $t$ when $M$ has
fixed size $p\times p$. Finally, \prettyref{assu:general-ass}\prettyref{enu:cov-mat-estimation}
ensures that we can estimate the covariance matrix of the residuals.
The following assumption for invariant subspaces defines the set of
matrices we work with.
\begin{assumption}
\label{assu:invar-sub-ass}Let $\mathscr{M}\subset\mathbb{R}^{p\times p}$
to be the set of matrices such that for every $M\in\mathscr{M}$,
\end{assumption}
\begin{enumerate}
\item \label{enu:semisimple}There exists an invariant subspace of interest
spanned by the columns of $R_{I}\in\mathbb{R}^{p\times q}$ such that
$MR_{I}=R_{I}\Lambda_{I}$ for some matrix $\Lambda_{I}\in\mathbb{R}^{q\times q}$.
This requirement is equivalent to the existence of a non-singular
matrix $R\coloneqq\left[\begin{smallmatrix}R_{I} & R_{J}\end{smallmatrix}\right]\in\mathbb{R}^{p\times p}$
with equivalent matrices $L=R^{-\mathsf{T}}$ such that
\begin{equation}
R^{-1}MR=\begin{bmatrix}\Lambda_{I} & \Lambda_{IJ}\\
0 & \Lambda_{J}
\end{bmatrix}=:\Lambda\label{eq:basic-rel}
\end{equation}
where $\Lambda_{I}\in\mathbb{R}^{q\times q}$ and $\Lambda_{J}\in\mathbb{R}^{r\times r}$
for $r=p-q$.
\item \label{enu:no-evs}The eigenvalues of $\Lambda_{I}$ and $\Lambda_{J}$
denoted by $\mathcal{L}_{I}$ and $\mathcal{L}_{J}$ obey $\mathcal{L}_{I}\ensuremath{\cap}\mathcal{L}_{J}=\emptyset$
and $\lambda\in\mathcal{L}_{I}$ implies that $\lambda^{\ast}\in\mathcal{L}_{I}$,
i.e. $\mathcal{L}_{I}$ is closed under conjugation. Furthermore, $\abs{\lambda_{q}}>\abs{\lambda_{q+1}}$.
\item \label{enu:rankRlu}$\operatorname{rk} G^{\mathsf{T}}R_{I}=q$ for a full-rank normalizing
matrix $G\in\mathbb{R}^{q\times p}$.
\end{enumerate}
\prettyref{assu:invar-sub-ass}\prettyref{enu:semisimple} defines
general invariant subspaces. The most important special case of those
subspaces are eigenspaces which obtain if $\Lambda_{IJ}=0$ and $\Lambda_{I}$
and $\Lambda_{J}$ are diagonal. An intermediate case is that of a
partially diagonalizable matrix, which we could achieve by requiring
$\Lambda_{I}$ to be diagonal with distinct eigenvalues, $\Lambda_{IJ}=0$,
and $\Lambda_{J}$ to be in Jordan normal form with blocks of arbitrary
size. The appeal of generic invariant subspaces is that they always
exist even in the most adverse circumstances. We have defined $\Lambda$
matrices to be real because even if eigenvalues appear as possibly
defective complex conjugates of another, it is always possible to
define the Jordan real form. For a conjugate pair, we can then find
an invariant subspace of twice the dimension of the multiplicity of
the eigenvalue associated with it.
Finally, \prettyref{assu:invar-sub-ass}\prettyref{enu:no-evs} ensures
that the map from $M$ to its invariant subspaces spanned by the columns
of $M$ is differentiable and that we can discriminate between vectors
of interest in sets $I$ and those in $J$. In particular, we do
not assume that we can discriminate among eigenvalues in $\mathcal{L}_{I}$
or $\mathcal{L}_{J}$.
An example of $\hat{M}$ that is consistent and has asymptotically
normal columns is the social interactions model \citep{depaula2023identifying,manresa}
for estimands $m_{ij}$ and $\gamma$, where the outcome
\[
y_{it}=\gamma\sum_{j\neq i}m_{ij}x_{jt}+\varepsilon_{it},
\]
for unit $i=1,\dots,p$ at time $t=1,\dots,T$ depends on the values
of individual-specific covariates $x_{it}$. Collecting estimands
$\hat{m}_{ij}$ into a matrix results in an estimated network adjacency
matrix whence all derived statistics such as centrality scores inherit
the uncertainty.\footnote{We consider statistics that can be derived from the spectrum of the
adjacency matrix but our methods extend straightforwardly to those
originating with graph Laplacians, too.} Furthermore, the estimation techniques in \citet{vliet} and \citet[Section 4.2]{rothenhausler2015backshift}
correspond to OLS estimation of $M$. Generally, OLS methods work
well if $T>>p$, which we explore in \prettyref{subsec:application}.
\subsection{Singular subspaces}
For singular subspaces, our setup changes ever so slightly. If $M$
is instead an $m\times l$ general, wide matrix so that $m<l$ and
of rank $q\leq\min\left(m,l\right)$, with $\Sigma_{I}\in\mathbb{R}^{q\times q}$.
The columns of $V_{I}\in\mathbb{R}^{l\times q}$ and $U_{I}\in\mathbb{R}^{m\times q}$
span left- and right singular subspaces of dimension $q$ iff
\begin{equation}
M=U_{I}\Sigma_{I}V_{I}^{\mathsf{T}}+U_{J}\Sigma_{J}V_{J}^{\mathsf{T}}\label{eq:svd}
\end{equation}
holds for $U_{J}\in\mathbb{R}^{m\times\left(l-q\right)}$, $V_{J}\in\mathbb{R}^{\left(l-q\right)\times l}$,
and $\Sigma_{J}\in\mathbb{R}^{\left(l-q\right)\times\left(l-q\right)}$.
It is also instructive to consider that we could implement inference
on $U_{I}$ via the procedure in \citet{Tyler1981} applied to $MM^{\mathsf{T}}$.
While using eigenvectors of $MM^{\mathsf{T}}$ may make no difference
asymptotically, it may be less efficient because $\operatorname{var}\text{vec\ } MM^{\mathsf{T}}\geq\operatorname{var}\text{vec\ } M$.
In analogy to \prettyref{assu:invar-sub-ass}, we have for the singular
subspace decomposition
\begin{assumption}
\label{assu:svd-ass}Define $\mathscr{S}$ such that for every $M\in\mathscr{S}$,
$M$ is an $m\times l$ real matrix with $m\leq l$, and of rank $q\leq m$.
Then,
\begin{enumerate}
\item \label{enu:svd-split}The singular subspaces split according to $M=U_{I}\Sigma_{I}V_{I}^{\mathsf{T}}+U_{J}\Sigma_{J}V_{J}^{\mathsf{T}}$
where $\operatorname{rk} M=q$ implies that $\Sigma_{I}\in\mathbb{R}^{q\times q}$
is a diagonal, square matrix and $\Sigma_{J}=0$.
\item For a singular subspace of dimension $q$, we have $U_{I}\in\mathbb{R}^{m\times q}$
where $q\leq F$.
\end{enumerate}
\end{assumption}
\prettyref{assu:svd-ass} defines a general setting that also covers
graph models, where $M$ may be sparse so that $q<<\min\left(m,l\right)$.
\section{Hypothesis tests}
\label{sec:hyp-tests}
We present Wald and $t$-tests for inference on basis vectors of invariant
and singular subspaces. Proofs of results appear in \prettyref{app:results-hyp-tests}.
To test $H_{0}$ in \eqref{eq:basic-hypothesis-random-matrix}, we
construct the orthocomplement to the candidate vectors $\upsilon_{\perp}$
so that under the null hypothesis, $\upsilon_{\perp}^{\mathsf{T}}R_{I}L_{I}^{\mathsf{T}}=0$.
Let $m$ denote the column dimension of $\upsilon_{\perp}$ which
is chosen by the researcher based on hypothesised vectors. Normally,
for a hypothesised eigenspace of dimension $q$, we have $m=p-q$
columns in $\upsilon_{\perp}$ although we could add additional columns
so that $m$ could exceed $p-q$.
Letting $\hat{\Omega}_{W}^{+}$ denote a generalized inverse of $\hat{\Omega}_{W}$,
defined in \prettyref{eq:gen-inv-eig} and \prettyref{eq:gen-inv-svd},
for $\hat{\Omega}_{W}:=\hat{B}\hat{\Omega}\hat{B}$, we have
\begin{equation}
\hat{W}_{n}\left(\upsilon_{\perp}\right)\coloneqq T\text{vec\ }\left(\upsilon_{\perp}^{\mathsf{T}}\hat{R}_{I}\hat{L}_{I}^{\mathsf{T}}\right)^{\mathsf{T}}\hat{\Omega}_{W}^{+}\text{vec\ }\left(\upsilon_{\perp}^{\mathsf{T}}\hat{R}_{I}\hat{L}_{I}^{\mathsf{T}}\right)\label{eq:waldstat-1}
\end{equation}
and
\begin{equation}
B\coloneqq\left(L_{I}\otimes\upsilon_{\perp}^{\mathsf{T}}R_{J}\right)\left\{ \left(\Lambda_{I}^{\mathsf{T}}\otimes I_{J}\right)-\left(I_{I}\otimes\Lambda_{J}\right)\right\} ^{-1}\left(R_{I}^{\mathsf{T}}\otimes L_{J}^{\mathsf{T}}\right).\label{eq:jacobian}
\end{equation}
into which we can insert sample analogues. The Wald test is based
on a first-order approximation of how the estimated subspace $R_{I}$
changes in response to perturbations in the matrix estimate $\hat{M}$,
with the Jacobian $B$ capturing this response. In our next result,
we focus on inference but \prettyref{sec:jacobians-ho-perturbations}
lays out Jacobians and further analytic results in more detail. For
eigenspaces, $\Lambda_{IJ}=0$ which does not affect the test statistic
but the speed with which test statistics converge as we see in \prettyref{subsec:higher-order-dk}.
Therefore, the approximations to test statistics are robust to $\Lambda_{IJ}\neq0$
and converge fastest for eigenspaces.
\begin{thm}
\label{thm:waldstatdistrib}Suppose \prettyref{assu:general-ass}
holds. Then, $\hat{W}_{n}\left(\upsilon_{\perp}\right)\ensuremath{\rightsquigarrow}\chi_{qm}^{2}.$
\end{thm}
\noindent The proof appears in \prettyref{app:distributions}.
Analogously to \eqref{eq:basic-hypothesis-random-matrix}, we write
hypotheses for singular vectors as $H_{0}^{\ast}:\upsilon_{\perp}^{\mathsf{T}}U_{I}=0$
for a specified $m\times(m-h)$ matrix $\upsilon_{\perp}$ corresponding
to some $\upsilon_{0}$. For the covariance, define
\begin{align}
B_{\text{SVD}}= & \left(\Sigma_{I}^{-1\mathsf{T}}V_{I}^{\mathsf{T}}\otimes\upsilon_{\perp}^{\mathsf{T}} U_{J}U_{J}^{\mathsf{T}}\right)\label{eq:svd-jacobian}\\
& +\left(I_{F}\otimes\upsilon_{\perp}^{\mathsf{T}} U_{I}\right)\left[\left(\ensuremath{\mathbf{1}}^{\mathsf{T}}\otimes\text{vec\ } D_{\text{d}}\right)\cdot\left(\left(\Sigma_{I}^{\mathsf{T}}V_{I}^{\mathsf{T}}\otimes U_{I}^{\mathsf{T}}\right)+\left(U_{I}^{\mathsf{T}}\otimes\Sigma_{I}V_{I}^{\mathsf{T}}\right)\right)\right].\nonumber
\end{align}
In analogy to inference on invariant subspaces above, $\upsilon_{\perp,0}^{\mathsf{T}}\in\mathbb{R}^{m\times\left(m-h\right)}$
where $h\leq F$ is the size of the hypothesized singular subspace.
We define the covariance matrix $\hat{\Omega}_{S}:=\widehat{B_{\text{SVD}}}\hat{\Omega}\widehat{B_{\text{SVD}}}^{\mathsf{T}}$
and denote by $\hat{\Omega}_{S}^{+}$ a generalized inverse of $\hat{\Omega}_{S}$.
The statistic is then
\begin{equation}
\hat{W}_{\text{SVD},n}\left(\upsilon_{\perp}\right)\coloneqq T\text{vec\ }\left(\upsilon_{\perp}^{\mathsf{T}}\hat{U}_{I}\right)^{\mathsf{T}}\hat{\Omega}_{S}^{+}\text{vec\ }\left(\upsilon_{\perp}^{\mathsf{T}}\hat{U}_{I}\right)\label{eq:waldstat-svd}
\end{equation}
for which we obtain
\begin{thm}
\label{thm:svd-wald-test}Suppose \prettyref{assu:svd-ass} holds,
then $\hat{W}_{\text{SVD},n}\left(\upsilon_{\perp}\right)\ensuremath{\rightsquigarrow}\chi_{q\left(m-h\right)}^{2}$.
\end{thm}
\prettyref{thm:svd-wald-test} is directly analogous to \prettyref{thm:waldstatdistrib}
and the proof is identical and available in \prettyref{app:results-hyp-tests}.
\noindent Although straightforward to compute, a Wald test statistic
is only required when a researcher specifies a full vector or matrix
hypothesis. We imagine that this scenario is most likely of interest
when testing whether localized basis vectors belong to an invariant
or singular subspace of $M$. Empirical applications may require inference
on individual coefficients of invariant or singular vectors, in which
case a scalar $t$-test is useful.
\noindent To make invariant subspace coordinates unique, ensuring
any associated estimates are consistent across experiments, we normalise
$\upsilon_{\perp}=:\begin{smallmatrix}[\upsilon_{1} & \upsilon_{2}]\end{smallmatrix}$
where $\upsilon_{1}\in\mathbb{R}^{q\times q}$ and $\upsilon_{2}\in\mathbb{R}^{r\times q}$.
Then letting $D^{\mathsf{T}}\coloneqq\upsilon_{2}\upsilon_{1}^{-1}\in\mathbb{R}^{r\times q}$,
we have
\begin{equation}
\upsilon_{\perp}\coloneqq\begin{bmatrix}I_{q}\\
-D^{\mathsf{T}}
\end{bmatrix}\label{eq:eigenvector-normalization}
\end{equation}
\noindent If an eigenvector represents a node's centrality score,
this normalization corresponds to choosing a numeraire node with unit
centrality so that coordinates are comparable across samples. However,
this normalization may cause outliers making it potentially unstable.
In our simulations, we did not find any evidence of any problems.
Analogous to \prettyref{eq:eigenvector-normalization}, we can derive
a closed-form expression of $\hat{D}_{n}$ appearing in $R_{I}$.
Let $R_{I,1}\in\mathbb{R}^{r\times q}$ and $R_{I,2}\in\mathbb{R}^{q\times q}$
with $\operatorname{rk} R_{I,2}=q$ so that
\begin{equation}
\begin{bmatrix}R_{I,1}\\
R_{I,2}
\end{bmatrix}\coloneqq R_{I}.\label{eq:partitionforahat-1}
\end{equation}
Consequently,
\begin{equation}
\hat{D}_{I,n}^{\mathsf{T}}\coloneqq\hat{R}_{I,1,n}\hat{R}_{I,2,n}^{-1}.\label{eq:ahat-1}
\end{equation}
so the normalized vector is $\begin{smallmatrix}[-D & I_{q}]^{\mathsf{T}}\end{smallmatrix}$.
Observe that \prettyref{eq:ahat-1} defines a unique estimator.
\noindent For inference on individual entries of $D$, define $d_{ij}\coloneqq e_{i}^{\mathsf{T}}De_{j}$
for $e_{i}\in\mathbb{R}^{q\times1}$ and $e_{j}\in\mathbb{R}^{r\times1}$
where the vectors $e_{i}$ have unit entries at $i$ and zero elsewhere.
Let $\hat{R}_{I,2}^{-1}$ be the empirical analogue of $R_{I,2}^{-1}$
defined in \prettyref{eq:partitionforahat-1} and define
\begin{equation}
B_{ij}=\left(e_{i}^{\mathsf{T}}R_{2,I}^{\mathsf{T}-1}\otimes e_{j}^{\mathsf{T}}\upsilon_{\perp}^{\mathsf{T}}R_{J}\right)\left[\left(\Lambda_{I}^{\mathsf{T}}\otimes I_{J}\right)-\left(I_{I}\otimes\Lambda_{J}\right)\right]^{-1}\left(R_{I}^{\mathsf{T}}\otimes L_{J}^{\mathsf{T}}\right),\label{eq:jac-coeff}
\end{equation}
in analogy to \eqref{eq:jacobian}. We write the scalar variance
of $\sqrt{T}(\hat{d}_{ij}-d)$ as $\sigma_{ij}^{2}\coloneqq B_{ij}\Omega B_{ij}^{\mathsf{T}}$
and define the $t$-test statistic for inference on the $ij$th coefficient
of $D$
\begin{equation}
t_{ij,n}\left(d_{0}\right)\coloneqq\frac{\hat{d}_{n,ij}-d_{0}}{\sqrt{\hat{\sigma}_{ij}^{2}/T}}.\label{eq:tstat-1-1}
\end{equation}
We construct the estimator of the variance, $\hat{\sigma}_{ij}$ by
replacing $B_{ij}$ with $\hat{B}_{ij}$ which contains sample analogues
$\hat{R}_{i}$ and $\hat{L}_{i}$ for $i\in\left\{ I,J\right\} $.
The distribution is given in
\begin{thm}
\label{thm:tdist}Suppose \prettyref{assu:general-ass} holds, then
$t_{ij,n}\left(d_{0}\right)\ensuremath{\rightsquigarrow} N\left(0,1\right).$
\end{thm}
\noindent Eigenvector-based centralities are based on the absolute
values $|\hat{d}_{ij}|$. We define for this result the folded normal
distribution with cumulative distribution function
\begin{equation}
F_{G}\left(x;d_{ij},\sigma_{ij}\right)\coloneqq\Phi\left(\frac{x-d_{ij}}{\sigma_{ij}}\right)-\Phi\left(\frac{-x-d_{ij}}{\sigma_{ij}}\right)\label{eq:cdf-folded}
\end{equation}
where $\Phi\left(\frac{x-\mu}{\sigma}\right)$ denotes the normal
cumulative distribution function with mean $\mu$ and variance $\sigma^{2}$.
For the asymptotic distribution of the absolute values of the normalized
eigenvector entries, we obtain
\begin{cor}[Folded normal distribution]
\label{cor:folded-normal}Suppose \prettyref{assu:general-ass} holds,
then $\sqrt{n}|\hat{d}_{ij}|\ensuremath{\rightsquigarrow} G$ where $G$ is a random variable
that has a folded normal distribution with c.d.f. \textup{given in
\eqref{eq:cdf-folded}.}
\end{cor}
\prettyref{subsec:Network-centralities-as} gives an application of
this result.
\subsection{Estimators for basis vectors\label{sec:estimators}}
Our next result characterizes the distribution of $\hat{D}_{I,n}^{\mathsf{T}}$.
Define
\begin{equation}
B_{D^{\mathsf{T}}}\coloneqq\left(R_{I,2}^{\mathsf{T}-1}\otimes\upsilon_{\perp}^{\mathsf{T}}R_{J}\right)\left\{ \left(\Lambda_{I}^{\mathsf{T}}\otimes I_{J}\right)-\left(I_{I}\otimes\Lambda_{J}\right)\right\} ^{-1}\text{\ensuremath{\left(R_{I}^{\mathsf{T}}\otimes L_{J}^{\mathsf{T}}\right)}}.\label{eq:d-jac}
\end{equation}
Then, we obtain
\begin{thm}[Distribution of basis vectors]
\label{thm:estdist}Let \prettyref{eq:ahat-1} define $\hat{D}_{n}$.
Then, \prettyref{assu:general-ass} implies $n^{1/2}\text{vec\ }\{\hat{D}_{n}^{\mathsf{T}}-D^{\mathsf{T}}\}\ensuremath{\rightsquigarrow} N(0,B_{D^{\mathsf{T}}}\Omega B_{D^{\mathsf{T}}}^{\mathsf{T}}).$
\end{thm}
\noindent The above results rely on consistent covariance matrix estimators
which we discuss alongside other proofs in \prettyref{app:distributions}.
To make precise statements about smoothness of invariant and singular
subspace maps, we define $\psi\left(M;\upsilon\right)$ that represents
invariant subspace statistics and has informative null distributions
on spaces spanned by a subset of the columns of $R$ or $U$. Let
$\psi\,:\,\set M\rightarrow\mathbb{R}^{r\times q}$ for $M\in\mathscr{M}$
in \prettyref{assu:invar-sub-ass} by
\begin{equation}
\psi\left(M,\upsilon_{\perp}^{\mathsf{T}}\right)\coloneqq\upsilon_{\perp}^{\mathsf{T}}R_{I}L_{I}^{\mathsf{T}}\left(M\right)\label{eq:inv-subspace-map}
\end{equation}
and analogously $\Psi(M,\upsilon_{\perp}^{\mathsf{T}})\coloneqq\upsilon_{\perp}^{\mathsf{T}}U_{I}$.
We define the map $\Psi:\mathscr{S\rightarrow\mathbb{R}}^{r\times q}$
via
\begin{equation}
\Psi\left(M,\upsilon_{\perp}^{\mathsf{T}}\right)\coloneqq\upsilon_{\perp}^{\mathsf{T}} U_{I}\left(M\right),\label{eq:sin-subspace-map}
\end{equation}
so that $\Psi\left(M,\upsilon_{0,\perp}^{\mathsf{T}}\right)=0$ under
the null hypothesis.\footnote{Without loss of generality, we restrict ourselves to left-singular
vectors as inference on $\operatorname{col} V$ follows from $M^{\mathsf{T}}$ instead
of $M$.}
The zero-level set $\psi(M,\upsilon_{0,\perp}^{\mathsf{T}})=0$ defines
non-rejection regions for $\upsilon_{\perp}$ and $\upsilon$. We
have set estimators
\begin{equation}
\widehat{\upsilon}_{\perp,n}\coloneqq\left\{ \upsilon_{\perp}\in\mathbb{R}^{p\times m}\,:\,\hat{\Phi}\left(\upsilon_{\perp}^{\mathsf{T}}\right)=0\right\} \label{eq:betahat-1}
\end{equation}
where
\[
\hat{\Phi}\left(\upsilon_{\perp}^{\mathsf{T}}\right)\coloneqq\begin{cases}
\psi\left(\hat{M}_{n},\upsilon_{\perp}^{\mathsf{T}}\right) & \text{invariant subspace,}\\
\Psi\left(\hat{M}_{n},\upsilon_{\perp}^{\mathsf{T}}\right) & \text{singular subspace.}
\end{cases}
\]
Estimators of $\psi\left(M\right)$ obtain from evaluating $\psi$
at the estimate $\hat{M}$, so that
\[
\psi\left(\hat{M}_{n}\right)=\upsilon_{\perp}^{\mathsf{T}}\hat{P}_{I}\left(\hat{M}_{n}\right).
\]
Using \prettyref{lem:jacobian}\prettyref{enu:invariant}, we obtain
the distribution of the sample analogue of $\psi(M)$ in
\begin{lem}
\label{lem:proj-distrib}Suppose \prettyref{assu:general-ass} holds
for invariant and additionally \prettyref{assu:svd-ass} holds for
singular subspaces. Let $B$ be as in \prettyref{eq:jacobian} or
\prettyref{eq:svd-jacobian} and let $\upsilon_{\perp}$ be as in
\prettyref{eq:per-def}. Then,\textup{
\[
\sqrt{T}\,\text{vec\ }\hat{\Phi}\left(\upsilon_{\perp}^{\mathsf{T}}\right)\ensuremath{\rightsquigarrow} N\left(0,B\Omega B^{\mathsf{T}}\right).
\]
}
\end{lem}
For general eigenvectors, i.e. if the set of eigenvalues $\mathcal{L}_{I}$
is not closed und\ensuremath{\le}er conjugation, $\sqrt{n}\,\text{vec\ }\psi(\hat{M}_{n})$
converges to a multivariate complex normal distribution, with covariance
$B\Omega B^{\mathsf{T}}$ and relation $B\Omega B^{\prime}$. For now,
we shall continue to assume that $\mathcal{L}_{I}$ is closed under conjugation
or, in the case of $\abs{\mathcal{L}_{I}}=1$, that the Perron-Frobenius
theorem applies although our results also accommodate these cases.
\section{Jacobians and smoothness\label{sec:jacobians-ho-perturbations}}
\subsection{\label{subsec:expansion-psi}Invariant subspace map $\psi$}
In the background, the hypothesis tests in \prettyref{sec:hyp-tests}
relied on a first-order Taylor expansion argument to find standard
errors. To make our arguments about higher-order perturbations to
invariant subspaces precise, we use the map $\psi$ that assigns such
spaces to non-symmetric matrices. \prettyref{assu:invar-sub-ass}\prettyref{enu:no-evs}
ensures that eigenvalues are semi-simple within their respective groups,
which allows differentiability of the invariant subspaces spanned
by columns of $R$ belonging to $\mathcal{L}_{I}$. This setup allows eigenvalues
to coincide within $\mathcal{L}_{I}$ so that our setup covers matrices
whose Jordan blocks have non-unit size, i.e. are not simple. Generally,
such eigenvalues imply that the corresponding eigenvectors are not
differentiable because small perturbations in the underlying matrix
have unpredictable consequences for the corresponding eigenvectors.
Our baseline setup covers exactly this case while some results in
\prettyref{sec:centralities} reduce to the case of simple eigenvalues.
In this section, we establish that $\psi$ and its singular subspace
companion $\Psi$ are infinitely differentiable: as their arguments
change slightly, the maps respond in a controlled and predictable
way enabling expansions of arbitrary order. Smoothness is therefore
not just a technical regularity: infinite differentiability enables
refined inference, including higher-order likelihood and Edgeworth
expansions. Our main novelty lies in achieving smoothness and deriving
Jacobians for generic invariant subspaces without assuming symmetry
or diagonalisability, extending the classical results. Unsurprisingly,
smoothness applies to singular subspace maps, too. Furthermore, we
characterise higher-order bounds of invariant subspaces. For estimations
of bounds of singular subspace maps beyond second order, we recommend
adapting invariant subspace expressions\footnote{Invariant subspaces of $MM^{\mathsf{T}}$ are singular subspaces of $M$.}
because of the tractability of the perturbation theory for invariant
subspaces.
To summarize our results succinctly, we define the resolvent map $\mathscr{P}\rightarrow\mathbb{R}^{p\times p}$,
for
\begin{equation}
X\mapsto\left(X-\alpha I\right)^{-1}\label{eq:res}
\end{equation}
which forms the basis of all spectral-based statistics. The main result
is
\begin{thm}
\label{thm:resolvent}The resolvent from $\mathscr{P}\rightarrow\mathbb{R}^{p\times p}$
in \prettyref{eq:res} defines a smooth map.
\end{thm}
This result guarantees that network statistics based on \prettyref{eq:res},
or, equivalently, those based on spectral data of $X$, have existent
higher-order expansions and hence inherit a certain stability.\footnote{If $X$ is less than full rank, then a generalized inverse can be
used without altering the result.} Furthermore, all statistics based on \eqref{eq:res} are consistently
estimable. Because singular vectors of $M$ are eigenvectors of $MM^{\mathsf{T}}$,
\prettyref{thm:resolvent} applies to $\Psi$ as well. A proof of
\prettyref{thm:resolvent} appears in \prettyref{app:smooth-maps-perturb-expansions}.
Jacobians appear in
\begin{lem}[Jacobians]
\label{lem:jacobian}~
\end{lem}
\begin{enumerate}
\item \label{enu:invariant}Let \prettyref{eq:inv-subspace-map} define
$\psi$. Then, the Jacobian matrix with respect to $M$ in the direction
of $E\in\mathbb{R}^{p\times p}$ of $\psi$ is $B\in\mathbb{R}^{p^{2}\times rq}$
such that$\ensuremath{\,\ensuremath{\mathrm{d}}}\text{vec\ }\psi(M,E)=B\text{vec\ } E$ for $B$ in \eqref{eq:jacobian}.
\item \label{enu:singular}Let \prettyref{eq:sin-subspace-map} define $\Psi$.
Then, the Jacobian matrix with respect to $M$ in the direction of
$E\in\mathbb{R}^{m\times l}$ of $\Psi$ is $B_{\text{SVD}}\in\mathbb{R}^{F\left(m-h\right)\times ml}$
such that $\ensuremath{\,\ensuremath{\mathrm{d}}}\text{vec\ }\Psi\left(M;E\right)=B_{\text{SVD}}\text{vec\ } E$ for $B_{\text{SVD}}$
in \eqref{eq:svd-jacobian}.
\end{enumerate}
\noindent\prettyref{app:Additional-details-on-invar} contains a
proof. Importantly, \eqref{eq:jacobian} reduces to $C_{w}$ in \citet[4.3]{Tyler1981}
when $M$ is symmetrizable. In this case, the columns of $R_{I}$
and $L_{I}$ satisfy $r_{i}l_{i}^{\mathsf{T}}\equiv r_{i}r_{i}^{\mathsf{T}}\Gamma$
where $\Gamma$ is a real positive definite symmetric matrix such
that $M\Gamma$ is symmetric. If such $\Gamma$ is found, the distinction
between left and right eigenvectors disappears as $\Gamma r_{i}$
takes the place of $l_{i}$ for $i\in I$.
\subsection{Extension of Davis-Kahan bound to higher-order perturbations}
\label{subsec:higher-order-dk}
The Jacobian in \eqref{eq:jacobian} reduces to the bound in \citet{demetrius},
which asserts that small perturbations produce changes bounded by
the ratio of the perturbation and the size of the spectral gap. For
a $2\times2$ matrix, the first-order term
\begin{align}
\norm{\ensuremath{\,\ensuremath{\mathrm{d}}}\psi\left(M,E\right)}_{\text{F}} & \leq\frac{\norm E_{\text{F}}}{\lambda_{1}-\lambda_{2}}\label{eq:davis-kahan-disguise}
\end{align}
where $\text{d}\psi$ is the infinitesimal response to perturbation
$E$. This bound differs from the original found in \citet{davis1969some}
by the fact that only population eigenvalues appear in the denominator.
The relation \eqref{eq:davis-kahan-disguise} inspires the question
whether we can use perturbative arguments to derive higher-order bounds.
The answer is affirmative and to simplify the exposition, we define
the size of the perturbation
\begin{align*}
a & \coloneqq\max\left\{ \norm{\hat{M}-M}_{\text{F}},\norm{L_{j}^{\mathsf{T}}\left(M-\hat{M}\right)R_{i}}\right\} \\
& =\norm{\hat{M}-M}_{\text{F}}
\end{align*}
for $i,j\in\left\{ I,J\right\} $ where the equality follows from
the proof of \prettyref{thm:stoch-bound-prop} and the discussion
in \prettyref{app:inequality-block-diagonal}. We define the eigenvalue
separation as
\[
s\coloneqq\norm{\left\{ \left(\Lambda_{I}^{\mathsf{T}}\otimes I_{J}\right)-\left(I_{I}\otimes\Lambda_{J}\right)\right\} ^{-1}}_{\text{F}}.
\]
Similarly, we have the distance
\[
\alpha_{12}\coloneqq\norm{\left(L^{\mathsf{T}}MR\right)_{1:p-q,q+1:p}}_{\text{F}}
\]
which measures how far $L^{\mathsf{T}}MR$ is from being block-diagonal.
For eigenvectors, higher-order bounds are simply the first-order bound
raised to a higher power so that control over $as$ implies control
over all orders. For general invariant subspaces $\alpha_{12}\neq0$
which amplifies higher-order effects, revealing a qualitative difference
in robustness.
Denote by $h(s^{x},a^{y})$ a polynomial with leading powers $x$
and $y$ in $s$ and $a$, respectively. Then, the first-order bound
$\gamma_{1}$ is given by \eqref{eq:davis-kahan-disguise}. While
first-order error is known to scale with $a/s$, it it is unclear
how higher-order terms behave for eigenspaces or invariant subspaces.
For general invariant subspaces, $\alpha_{12}>0$ and the eigenvalue
separation $s$ has a progressively stronger effect on higher-order
error terms than the estimation precision. This relationship contrasts
sharply with the eigenspace case where $\alpha_{12}=0$, where higher-order
bounds are simply powers of the first-order term.
\begin{thm}[Higher-order Davis-Kahan bounds]
\label{thm:higher-order-dk}The first order bound in \eqref{eq:davis-kahan-disguise}
extends to higher orders according to
\begin{align*}
\gamma_{2} & \approx h\left(s^{2},a^{2}\right)+\ensuremath{\mathbf{1}}\left\{ \alpha_{12}>0\right\} \alpha_{12}s\gamma_{2}\approx h\left(s^{2},a^{2}\right)+\alpha_{12}\ensuremath{\mathbf{1}}\left\{ \alpha_{12}>0\right\} h\left(s^{3},a^{2}\right)\\
\gamma_{3} & \approx h\left(s^{3},a^{3}\right)+\ensuremath{\mathbf{1}}\left\{ \alpha_{12}>0\right\} \alpha_{12}^{2}\gamma_{3}\approx h\left(s^{3},a^{3}\right)+\alpha_{12}^{2}\ensuremath{\mathbf{1}}\left\{ \alpha_{12}>0\right\} s^{2}h\left(s^{3},a^{3}\right)\\
\gamma_{4} & \approx h\left(s^{4},a^{4}\right)+\ensuremath{\mathbf{1}}\left\{ \alpha_{12}>0\right\} \alpha_{12}^{3}h\left(s^{7},a^{4}\right)\\
& \vdots\\
\gamma_{k} & \approx h\left(s^{k},a^{k}\right)+\ensuremath{\mathbf{1}}\left\{ \alpha_{12}>0\right\} \alpha_{12}^{k-1}h\left(s^{2k-1},a^{k}\right).
\end{align*}
For eigenspaces $\alpha_{12}=0$, in which case the higher order bounds
are homogeneous polynomials in $as$.
\end{thm}
The approximations to higher-order terms only list the highest power
of $\alpha_{12}$ but a term linear in $\alpha_{12}$ nevertheless
appears and is a polynomial with leading powers $k$ in $s$ and $a$,
($h(s^{k},a^{k})$). The first two equalities follow from $\gamma_{2}=O(s^{2}a^{2})$
multiplying $s$and $\gamma_{3}=O(s^{3}a^{3})$ multiplying $s^{2}$.
The lesson of \prettyref{thm:higher-order-dk} is that invariant subspace
estimation converges slowly if $\alpha_{12}>1$. Moreover, a quickly
converging estimator with $a\ensuremath{\rightarrow}0$ cannot easily balance out a
larger reciprocal eigenvalue separation $s$ because its higher-order
contributions rise faster in $s$ than they do in $a$. Specifically,
its powers dominate with $2k-1$ relative to $k$ for $a$ so that
statisticians estimating invariant subspaces are advised to check
$\alpha_{12}$ whenever Davis-Kahan arguments are applied. A proof
is given in \prettyref{app:perturb-to-deriv}.
\section{Network statistics\label{sec:centralities}}
In this section, we study the scenario where a researcher has an estimator
$\hat{M}$ of $M$ that fulfills $\smlnorm{\hat{M}-M}_{\text{F}}=O_{p}\left(r_{n}\right)$
for some sequence $r_{n}$. Using $r_{n}$ as an input, we find the
stochastic order of $\smlnorm{\hat{\psi}-\psi}$ for subspace-based
network statistics $\psi$ as well as the node-wise and overall clustering
coefficients. Moreover, we apply the inference methods of \prettyref{sec:hyp-tests}
to the problem of finding standard errors for network centrality measures,
all of which are based on invariant subspace decompositions of the
network adjacency matrix. Proofs appear in \prettyref{app:results-network-centralities}.
\subsection{\label{subsec:Network-centralities-as}Convergence rates of network
statistics}
For our results to apply, the leading eigenvalue of the adjacency
matrix has to be distinct from all the others. In this case, \prettyref{assu:invar-sub-ass}\prettyref{enu:no-evs}
becomes the singleton $\mathcal{L}_{I}=\{\lambda_{\text{max}}\}$. Both
this assumption as well as the Perron-Frobenius theorem in the case
of adjacency matrices ensure that this condition is met. The application
of the Perron-Frobenius theorem here is crucial because it guarantees
that the top eigenvalue is simple so that the single top eigenvector
map is differentiable.
\subsection*{Convergence rate of $\psi$ as a network statistic for matrices with
potentially growing size}
The map $X\mapsto(X-\alpha I)^{-1}$ forms the basis of a large
class of network centrality measures, which we term \emph{resolvent-based
centrality scores }and in \prettyref{thm:resolvent}, we establish
that it is infinitely differentiable. The immediate consequence is
that all network statistics admit expansions to arbitrary order. For
$g\in\{n,T\}$, recall \prettyref{assu:general-ass}\prettyref{enu:tightness},
which states that
\begin{equation}
\text{\ensuremath{\norm{\hat{M}-M}}}_{\text{F}}=O_{p}\left(r_{g}\right)\label{eq:stoch-bd}
\end{equation}
and suppose that a researcher has shown that \prettyref{eq:stoch-bd}
holds for some estimator. The question that naturally arises is to
what extent
\[
\norm{\psi\left(\hat{M}\right)-\psi\left(M\right)}_{\text{F}}
\]
is then bounded in probability. Naturally, $\psi(\hat{M})\ensuremath{\overset{p}{\ensuremath{\rightarrow}}}\psi(M)$
as $r_{g}\ensuremath{\rightarrow}0$ follows from \prettyref{thm:resolvent} as $\alpha\ensuremath{\rightarrow}\lambda_{\text{max}}$.\footnote{As $\alpha\ensuremath{\rightarrow}\lambda_{\text{max}}$, the resolvent approaches
the top eigenvector. See the proof of \prettyref{lem:equivalence-centr}
in Appendix \ref{subsec:network-centrality-proofs} for an explanation.} But we wish to quantify the exact convergence rate. Denote by $\psi_{1}$
a statistic based on the top eigenvector of $M$. Then, we have
\begin{thm}
\label{thm:stoch-bound-prop}Suppose the stochastic bound \eqref{eq:stoch-bd}
(\prettyref{assu:general-ass}\prettyref{enu:tightness}) holds and
that the Perron-Frobenius theorem applies to $M$. Then, the following
apply:
\begin{enumerate}
\item Statistics based on the principal component of $M$, $\psi_{1}$,
obey
\begin{equation}
\norm{\psi_{1}\left(\hat{M}\right)-\psi_{1}\left(M\right)}_{\text{F}}=O_{p}\left(\frac{r_{g}}{\lambda_{1}-\lambda_{2}}\right).\label{eq:principal-comp-bound}
\end{equation}
\item If the statistic encompasses more than one eigenvector (invariant
vector), let $1\leq r\leq s\leq\operatorname{rk} A$ and $I=\left\{ r,r+1,\dots,s\right\} $
with ordered eigenvalues $\lambda_{r}<\lambda_{r+1}<\dots<\lambda_{s}$
where $\operatorname{sp}\psi\left(\hat{M}\right)=\operatorname{sp}\left\{ R_{I}\right\} +o_{p}\left(1\right)$.
In this case, the bound changes to
\[
O_{p}\left(\frac{r_{g}}{\min\left\{ \left(\lambda_{r-1}-\lambda_{r}\right),\left(\lambda_{s}-\lambda_{s+1}\right)\right\} }\right).
\]
\end{enumerate}
\end{thm}
\noindent The bound attained in \prettyref{thm:stoch-bound-prop}
resembles the Davis-Kahan theory \citep{davis1969some,demetrius}
except that we derived it from perturbative arguments.
\begin{rem}
The result in \prettyref{thm:stoch-bound-prop} is valid for $M$
with fixed size $p\times p$ or $n\times n$ where $n$ is understood
to be possibly divergent.
\end{rem}
\noindent For eigenvector-based statistics, we can set $\alpha_{12}=0$
in \prettyref{thm:higher-order-dk} and combine the result with \prettyref{thm:stoch-bound-prop}
to a bound to $k$th order which reads
\begin{equation}
\sum_{i=1}^{k}O_{p}\left(\frac{r_{g}}{\min\left\{ \left(\lambda_{r-1}-\lambda_{r}\right),\left(\lambda_{s}-\lambda_{s+1}\right)\right\} }\right)^{i}.\label{eq:higher-order-bound}
\end{equation}
An important case arises if the Perron-Frobenius theorem does not
apply or if eigenvalues are not simple, so that individual eigenvectors
are no longer differentiable. In this case, we have to focus the analysis
on generic invariant subspaces, where $\alpha_{12}\neq0$ in \prettyref{thm:higher-order-dk}.
It is possible to extend \prettyref{thm:stoch-bound-prop} to these
cases.
\subsection*{Clustering coefficient}
Follow \citet{jackson2008social} and define the node-wise clustering
coefficient $\text{cl}_{i}$ via
\begin{equation}
\text{cl}_{i}\coloneqq\frac{\sum_{j\neq i,k\neq j,k\neq i}m_{ij}m_{ik}m_{jk}}{\sum_{j\neq i,k\neq j,k\neq i}m_{ij}m_{ik}},\label{eq:clustering-coeff}
\end{equation}
where the sums run over all indices but $i$. If we are dealing with
the overall clustering coefficient we have instead
\[
\text{cl}\coloneqq\frac{\sum_{i,j\neq i,k\neq j,k\neq i}m_{ij}m_{ik}m_{jk}}{\sum_{i,j\neq i,k\neq j,k\neq i}m_{ij}m_{ik}}.
\]
Estimators obtain from replacing $m_{ij}$ with elements of $\hat{M}$.
\begin{thm}[Convergence rate of clustering coefficients]
\label{thm:convergence-rate-clustering-coeff}Suppose the stochastic
bound \eqref{eq:stoch-bd} (\prettyref{assu:general-ass}\prettyref{enu:tightness})
holds. For the node-wise and individual clustering coefficients, $\smlabs{\hat{\text{\text{cl}}}_{i}-\text{cl}_{i}}=O_{p}\left(r_{n}\right)$
and $\smlabs{\hat{\text{\text{cl}}}-\text{cl}}=O_{p}\left(r_{n}\right).$
\end{thm}
\prettyref{thm:convergence-rate-clustering-coeff} shows that the
node-wise and overall clustering coefficients converge at the same
rate as the adjacency matrix estimates.
\subsection{Network centrality measures}
In the following, we shall use standard errors for centrality scores
to construct $t$-tests. Note that these all centrality scores satisfy
the convergence rate given in \prettyref{thm:stoch-bound-prop}. Recall
that standardisation by $g\in\{n,T\}$ depends on whether $M$ has
growing or fixed size where \citet[Example 2.3]{benaychgeorges2018lectureslocalsemicirclelaw}
provide details on required normalization.
\subsection*{Eigenvector centrality}
This centrality score $c_{i}$ assigns popularity to node $i$ as
the sum of the popularity of its neighbors. Intuitively speaking,
if someone is connected to very popular nodes, one is themselves very
popular so that $c_{i}\coloneqq\frac{1}{\lambda}\sum_{j\in N_{i}}c_{j}.$
The sum runs over all $j$ that are connected to node $i$, denoted
by $N_{i}$ and $\lambda$ is a normalizing constant. For the full
vector of scores $c$, we make the substitution $\sum_{j\in N_{i}}c_{j}=\sum_{j=1}^{N}m_{ij}c_{j}$
and obtain the eigenvector problem $Mc=\lambda c$. Because $c$ is
a basis vector belonging to the eigenspace of $M$ associated with
the largest eigenvalue of $M$, we can apply the results of \prettyref{sec:estimators}
directly to this network statistic. Importantly, this measure is not
invariant to normalizations. We thus suggest applying the normalization
in \prettyref{eq:eigenvector-normalization} which defines centralities
uniquely by denoting one node as the reference and use the test in
\prettyref{thm:tdist} for inference on individual coefficients. Formally,
we have
\begin{thm}
\label{thm:eigenvec-smooth}The eigenvector centrality $c\left(M\right)$
is a smooth function of $M$. Therefore, $\hat{M}_{n}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}} M$
implies $\hat{c}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}} c$ and for an adjacency matrix with distribution
$N\left(\text{vec\ } M,\Omega\right)$, we have
\[
\frac{\hat{c}_{i}-c_{i}}{\sqrt{B_{1j}\Omega B_{1j}^{\mathsf{T}}/g}}\ensuremath{\rightsquigarrow} N\left(0,1\right)
\]
for the Jacobian of $B_{1j}$ defined in \eqref{eq:jac-coeff}.
\end{thm}
The above result is a direct consequence of \prettyref{thm:estdist}
and allows one to construct standard errors. In applied work, we wish
to test hypotheses of the form $H_{0}:c_{i}=c_{0}$ against the alternative
that $c_{i}\neq c_{0}$. We use the absolute value version of the
coefficient-wise $t$-test based on \prettyref{cor:folded-normal}
to construct one-sided confidence intervals and define the one-sided
confidence interval with level $\alpha$ for $\abs{c_{i}}$\footnote{This centrality measure is usually implemented as the abs. value of
the coefficients, see e.g. \href{https://github.com/JuliaGraphs/Graphs.jl/blob/master/src/centrality/eigenvector.jl}{https://github.com/JuliaGraphs/Graphs.jl/blob/master/src/centrality/eigenvector.jl.}.} via$\mathcal{C}_{i}\coloneqq\left[\abs{\hat{c}_{i}},x_{0}\right],$where
the upper limit $x_{0}>0$ solves $F_{G}\left(x_{0};\abs{c_{i}},\sigma_{i}\right)=1-\alpha$
for $x_{0}>0$. Expression \prettyref{eq:cdf-folded} defines $F_{G}$
with consistent point estimates $\abs{\hat{c}_{i}}$ and $\hat{\sigma}_{i}$
substituted for $\abs{s_{i}}$ and $\sigma_{i}$ for $i=1,\dots,p-1$.
Consequently, $\forall c\in\mathcal{C}$, we cannot reject $H_{0}\,:\,\abs{c_{i}}=c$.
\subsection*{PageRank, Katz, and Diffusion centralities}
A related family of centrality measures is conceptually very similar
to the eigenvector centrality but performs better for dense graphs.
Formally, the PageRank centrality \citep{brin1998anatomy} $c_{P}$
is given by
\begin{align}
c_{P} & =(I-\alpha M)^{-1}\beta\label{eq:pagerank-def}
\end{align}
where $\alpha\in\mathbb{R}$ and $\beta\in\mathbb{R}^{p}$ are constants.
The damping factor $\alpha$ should be chosen to lie between $0$
and $1/\lambda_{\text{max}}\left(M\right)$. We can easily derive
\prettyref{eq:pagerank-def} from the eigenvector centrality by adding
a small amount of ``centrality'' to each node in the eigenvector
relation in the form of $\beta$ to obtain $c=\alpha Mc+\beta$ for
some choice of $\alpha$ and solving for $c$. To motivate \prettyref{eq:pagerank-def}
further, we imagine a walk on a graph for $S$ periods with the entries
of $c$ providing a measure of which nodes would be visited most often
in the iterative scheme $c^{\left(t\right)}=\alpha Mc^{\left(t-1\right)}+\beta$
where $c^{\left(t\right)}$ converges to $c_{P}$ for any initial
choice $c^{\left(0\right)}$ as long as $\alpha$ is chosen appropriately.
\citet{banerjee2013diffusion} introduce this idea formally via the
diffusion centrality as
\begin{equation}
c_{D}\coloneqq\left(\sum_{s=0}^{S}\alpha^{s}M^{s}\right)\ensuremath{\mathbf{1}},\label{eq:diffusion-def}
\end{equation}
which we can interpret as a ``finite'' walk PageRank centrality
with exogenous centrality $\beta=\ensuremath{\mathbf{1}}$. Its limiting case is the
Katz centrality defined via
\begin{align}
c_{K} & \coloneqq(I-\alpha M)^{-1}\ensuremath{\mathbf{1}},\label{eq:katz-def}
\end{align}
which we obtain by setting $\beta=\ensuremath{\mathbf{1}}$ in the PageRank centrality
\eqref{eq:pagerank-def}. We summarize the relationship between eigenvector,
Katz, PageRank, and Diffusion centralities in
\begin{lem}
\label{lem:equivalence-centr}~
\begin{enumerate}
\item \label{enu:diff-katz-equi}As the chains in the diffusion centrality
become infinitely long, it converges to the Katz centrality, i.e.
as $S\ensuremath{\rightarrow}\infty$, $c_{D}\ensuremath{\rightarrow} c_{K}$
\item \label{enu:eigenvec-equiv}As $\alpha\ensuremath{\rightarrow}1/\lambda_{1}$, Katz
and PageRank centralities converge to the eigenvector centrality,
i.e. $c_{K}\ensuremath{\rightarrow} c$ and $c_{P}\ensuremath{\rightarrow} c$.
\end{enumerate}
\end{lem}
For inference on PageRank and Katz centralities, we derive the Jacobian
$\ensuremath{\,\ensuremath{\mathrm{d}}} c_{P,i}=B_{P,i}\ensuremath{\,\ensuremath{\mathrm{d}}}\text{vec\ } M$ as
\begin{equation}
B_{P,i}=\left(-\alpha\beta^{\mathsf{T}}\left(I-\alpha M^{\mathsf{T}}\right)^{-1}\right)\otimes\left(e_{i}^{\mathsf{T}}\left(I-\alpha M\right)^{-1}\right),\label{eq:jac-pr}
\end{equation}
while for the diffusion centrality we have
\begin{equation}
B_{D,i}=\sum_{s=0}^{S}\alpha^{-s}\sum_{j=1}^{s}\left(\ensuremath{\mathbf{1}}^{\mathsf{T}}M^{\mathsf{T}}\right)^{s-j}\otimes e_{i}^{\mathsf{T}}M^{j-1}.\label{eq:jac-diff}
\end{equation}
The following result allows conducting inference on eigenvector-based
centrality measures.
\begin{thm}
\label{thm:pagerank-katz-smooth}For an adjacency matrix estimator
that satisfies $\sqrt{n}\,\text{vec\ }(\hat{M}_{n}-M)\ensuremath{\rightsquigarrow} N\left(0,\Omega\right)$,
we have the following results.
\begin{enumerate}
\item \label{enu:pr-katz}The PageRank and Katz centralities in \prettyref{eq:pagerank-def}
and \prettyref{eq:katz-def} are smooth functions of $M$. Furthermore,
$\hat{M}_{n}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}} M$ implies $\hat{c}_{P,i}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}} c_{P,i}$ and
for an adjacency matrix with distribution $N\left(\text{vec\ } M,\Omega\right)$,
we have the $t$-statistic
\[
\sqrt{g}\frac{\left(\hat{c}_{P,i}-c_{P,i}\right)}{\sigma_{P,i}}\ensuremath{\rightsquigarrow} N\left(0,1\right)
\]
for \textup{$\sigma_{P,i}=B_{P,i}\Omega B_{P,i}^{\mathsf{T}}$} defined
in \prettyref{eq:jac-pr}. Similarly, setting $\beta=\ensuremath{\mathbf{1}}$ in $\sigma_{P,i}$
lets us obtain the same result for the Katz centrality.
\item \label{enu:diffusion}The diffusion centrality in \prettyref{eq:diffusion-def}
is a smooth function of $M$. Furthermore, $\hat{M}_{n}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}} M$
implies $\hat{c}_{D,i}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}} c_{D,i}$. Further, we have
\[
\sqrt{g}\frac{\left(\hat{c}_{D,i}-c_{D,i}\right)}{\sigma_{D,i}}\ensuremath{\rightsquigarrow} N\left(0,1\right)
\]
for $\sigma_{D,i}=B_{D,i}\Omega B_{D,i}^{\mathsf{T}}$, defined in \prettyref{eq:jac-diff}.
\end{enumerate}
\end{thm}
From a computational point of view, the evaluation of the Jacobian
$B_{P,i}$ can be memory-intensive and sometimes researchers may prefer
to avoid inversion altogether. It is straightforward to construct
a variance estimate based on the Moore-Penrose inverse instead. See
\citet[Ch. 8, Thm. 5]{magnus2019matrix}.
\subsection*{Degree centrality}
Lastly, we examine degree centrality, defined as $c_{N}\coloneqq M\ensuremath{\mathbf{1}}$
for a vector of ones $\ensuremath{\mathbf{1}}$ of length $p$. It is straightforward
to compute the distribution of the degree centrality via
\begin{thm}
\label{thm:deg-centrality}The degree centrality satisfies $c_{N}=\left(\ensuremath{\mathbf{1}}^{\mathsf{T}}\otimes I_{p}\right)\text{vec\ } M$
and $\hat{M}_{n}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}} M$ implies that $\hat{c}_{N}\ensuremath{\overset{p}{\ensuremath{\rightarrow}}} c_{N}$.
For an adjacency matrix estimator that satisfies $\sqrt{n}\,\text{vec\ }\left(\hat{M}_{n}-M\right)\ensuremath{\rightsquigarrow} N\left(0,\Omega\right)$,
we have for the $i$th coefficient
\[
\frac{\hat{c}_{N,i}-c_{N,i}}{\sqrt{\sigma_{N,i}^{2}/g}}\ensuremath{\rightsquigarrow} N\left(0,1\right),
\]
for $\sigma_{N,i}^{2}=\left(\ensuremath{\mathbf{1}}\otimes e_{i}\right)^{\mathsf{T}}\Omega\left(\ensuremath{\mathbf{1}}\otimes e_{i}\right)$.
\end{thm}
\section{Applications and Simulations\label{sec:Simulation-study}}
Our results cover directed networks with potential self-loops, which
imply non-symmetric adjacency matrices. For the simplest case, the
entries $M_{ij}$ are equal to one if a connection exists between
$i$ and $j$ and zero otherwise. For directed graphs, $M$ is not
symmetric and in the weighted case, its entries are proportional to
the intensity of the connection. The applications demonstrate the
utility of our methods for eigenvector centralities but equally apply
to related statistics used, e.g. by clustering algorithms.
\subsection{\label{subsec:application}Network models of international trade
and input-output}
Suppose that an empirical researcher wishes to quantify the uncertainty
in the associated centrality estimates using the results of \prettyref{sec:centralities},
in particular \prettyref{thm:eigenvec-smooth}. To demonstrate versatility,
we follow computational convention applying \prettyref{cor:folded-normal}
to $\smlabs{d_{i}}$ implying one-sided intervals for one application,
and rely on the PF theorem for our second, leading to two-sided intervals.
\prettyref{fig:one-sided-intervals} shows plots of intervals obtains
for both examples for centralities of (a) trade and (b) input-output
networks.
\begin{figure}
\caption{\label{fig:one-sided-intervals}Confidence intervals ($95\%$) for
eigenvector centralities of (a) trade network (one-sided) and (b)
input-output network (two-sided).}
\centering
\begin{tabular}{cc}
\includegraphics[width=7cm,height=6cm,keepaspectratio]{centrality_intervals} & \includegraphics[width=7cm,height=6cm,keepaspectratio]{centrality_intervals_io_raw}\tabularnewline
(a) Centralities for EU trade network (logarithms). & (b) Centralities for input-output network.\tabularnewline
\end{tabular}
\justifying
\noindent{\footnotesize\emph{Notes}}{\footnotesize : The plots show
estimated centrality scores based on weighted, directed trade networks
estimated with error. The left panel shows how accounting for uncertainty
in the estimated trade links provides upper bounds for the centrality
scores, though leaves the importance ordering based on the point estimates
intact. In the right panel (b), we present estimated input-output
links based on the sectors of the American economy, where we ordered
the nodes by the upper end of the confidence interval. We see that
the point estimates vary substantially and could have led to misreporting
of node importance.}{\footnotesize\par}
\end{figure}
\prettyref{tab:sector-codes} links the abbreviations in \prettyref{fig:one-sided-intervals}
(b) to the relevant sector. We see from panel (a) that accounting
for uncertainty in the centrality score estimates leaves the importance
ranking among the countries intact, whereas it reorders it for the
sectors of the US economy when considering the upper bounds.
\begin{table}
\caption{\label{tab:sector-codes}Sector codes for the input-output network
model.}
\centering
\begin{tabular}{l|l}
\toprule Code & Description \\
\hline 23 & Construction \\ 31G & Manufacturing \\ 42 & Wholesale trade \\ 44RT & Retail trade \\ 48TW & Transportation and warehousing \\ FIRE & Finance, insurance, real estate, rental, and leasing \\ PROF & Professional and business services \\ 7 & Arts, entertainment, recreation, accommodation, and food services \\ 81 & Other services, except government \\ G & Government \\ \bottomrule
\end{tabular}
\justifying
\noindent{\scriptsize\emph{Notes}}{\scriptsize : This table shows
the sectors used in estimating the input-output network adjacency
matrix. Data originate with the Bureau of Economic Analysis.}{\scriptsize\par}
\end{table}
\prettyref{fig:estimated-networks} visualizes the estimated networks,
where arrow thickness denotes the weight and node size corresponds
to centrality. Panel (a) displays the EU trade network and gives some
intuition why Portugal is such a central node, despite having lower
weights than the others. What it lacks in trade volume, it compensates
with diversity of connections and being connected to other very central
nodes, as is intended for eigenvector centrality. In panel (b), we
see the estimated network of sectors of the US economy. There, trade
tends to be fairly balanced shown by equally thick arrows. Pointwise,
the retail sector (44RT) is the largest, but it is estimated with
quite a bit of noise. Arts, food, and entertainment (7) could be much
more central than retail despite its lower point estimate. The lesson
here is clear: it is likely that arts, food, and entertainment is
more volatile than retail, so that perhaps a time-dependent graph
model may be better able to capture the dependencies and reduce the
uncertainty thus quantified.
\begin{figure}
\caption{\label{fig:estimated-networks}Estimated networks from (a) trade data
and (b) input-output of sectors of the US economy. Arrow thickness
indicates trade volume while node size indicates estimated centrality
score.}
\centering
\begin{tabular}{c}
\resizebox{11.5cm}{!}{\includegraphics[width=9cm,height=9cm,keepaspectratio]{graph_w_centralities_trade_raw_simple_gr}}\tabularnewline
(a) Digraph of estimated EU pharmaceutical trade network. Date range
is from '04 to '18, EUROSTAT.\tabularnewline
\begin{cellvarwidth}[t]
\centering
\resizebox{11.5cm}{!}{
\includegraphics[width=9cm,height=9cm,keepaspectratio]{graph_w_centralities_io_raw_simple_larger}}
\end{cellvarwidth}\tabularnewline
(b) Digraph of estimated sectoral trade of US economy. Date range
is from '17 to '22, BEA.\tabularnewline
\end{tabular}
\justifying
\noindent{\footnotesize\emph{Notes: }}{\footnotesize The plots show
estimated networks where the arrow indicates the direction of trade.}{\footnotesize\par}
\end{figure}
\subsection{Digraph model}
We studied the finite sample performance of the $t$-test based on
a random, binary directed graph with six nodes and ten connections
displayed in panel (a) of \prettyref{fig:eigenvector-centrality}.
The Q-Q plot for the $t$-test statistic in panel (b) shows that the
approximation performs very well. Our method is therefore well-suited
for binary digraphs, too.
\begin{figure}
\caption{\label{fig:eigenvector-centrality}Example graph measured with noise
and quality of the asymptotic approximation for inference on eigenvector
centralities. $1000$ MC repetitions were used for a sample size of
500. Q-Q plots are theoretical ($y$) vs. empirical ($x$).}
\centering
\begin{tabular}{cc}
\includegraphics[height=6cm]{graph_plot} & \includegraphics[height=6cm]{qq_graph}\tabularnewline
(a) Digraph with adjacency matrix $M$. & (b) Quality of approximation.\tabularnewline
\end{tabular}
\justifying
\noindent{\footnotesize Notes: The left panel shows the graph that
gives rise to adjacency matrix $M$. It has six nodes and ten connections.
On the right panel, we see the Q-Q plot of the $t$-statistic for
inference on a single centrality score on the $x$-axis and the theoretical
quantiles on the $y$-axis.}{\footnotesize\par}
\end{figure}
\subsection{\label{subsec:Monte-Carlo-Evidence}Generic invariant subspaces}
To check how well our methods worked for inference on generic invariant
subspaces, we used a dense matrix with $p=2$ whose stacked columns
have simplified covariance matrix $\Omega=I_{p}\otimes\Omega_{M}$.
Full details on the data-generating process appear in \prettyref{alg:wald-and-t}
in \prettyref{app:results-network-centralities}. The online supplement
to this paper contains further DGPs as well as examples on how to
include heteroskedasticity and autocorrelation-robust covariance matrix
estimators. We ran further simulations to study the performance of
the $t$- and Wald tests using Q-Q plots as well as the empirical
cumulative distribution functions compared with their theoretical
counterparts. These visual aids demonstrate the quality of the asymptotic
approximations found in \prettyref{sec:hyp-tests}.\footnote{Complementary to the results presented in this section, we refer the
reader to the extra material hosted at \href{https://github.com/jsimons8/networkmodelssubspaces}{https://github.com/jsimons8/networkmodelssubspaces.}} The panels in \prettyref{fig:quality-approx-invar-sub} show that
both the Wald and $t$-test statistics perform well in simulation
exercises. Left panels let us judge the approximation made in \prettyref{thm:waldstatdistrib}
while right ones display the quality of the $t$-test approximation
of \prettyref{thm:tdist}. The overall performance is very good with
only few outliers. Histograms overlain with densities displayed similar
results and appear in the online supplement. The bottom panels show
the same pattern where the cumulative distribution functions track
their empirical counterparts well.
\begin{figure}
\caption{\label{fig:quality-approx-invar-sub}Quality of asymptotic approximation.
Left column: Wald test, right column: t-test. $5000$ MC repetitions
were used for a sample size of $100$. : Q-Q plots are theoretical
($y$) vs. empirical ($x$).}
\centering
\begin{tabular}{ccc}
\includegraphics[width=8cm,height=6cm,keepaspectratio]{qqplot2_largefont} & & \includegraphics[width=8cm,height=6cm,keepaspectratio]{qqt_bigfont}\tabularnewline
(a) Wald test quantiles. & & (b) $t$-test quantiles.\tabularnewline
\includegraphics[width=8cm,height=6cm,keepaspectratio]{cdfwald_largefont} & & \includegraphics[width=8cm,height=6cm,keepaspectratio]{cdft_bigfont}\tabularnewline
(c) CDF comparisons of Wald test statistic. & & (d) CDF comparisons of $t$-test statistic.\tabularnewline
\end{tabular}
\justifying
\noindent{\footnotesize\emph{Notes}}{\footnotesize : The figures show
Q-Q plots and CDFs for the Wald and $t$-test statistics to verify
the accuracy of \prettyref{thm:waldstatdistrib} and \prettyref{thm:tdist}.
The Q-Q plots visualise the outliers in the tails of the Wald tests
which are less pronounced for the $t$-test. However, even for the
Wald-test, the number of outliers is only moderate in light of the
$5,000$ MC repetitions. The bottom panels, (c) and (d), show empirical
and reference CDFs which show close tracking across the entire support.}
\end{figure}
\subsection{Singular subspace inference: Monte Carlo evidence}
We considered a data-generating process using a covariance matrix
$\Omega=I_{m}\otimes\Omega_{W}$ for a positive-definite $\Omega_{W}\in\mathbb{R}^{m\times m}$
for $m=3$. \prettyref{alg:wald-and-t-svd} in \prettyref{app:results-network-centralities}
details the steps of how we generate samples for the $t$- and Wald
statistics.
\begin{figure}
\caption{\label{fig:quality-approx-svd}Quality of asymptotic approximation
for Wald test-based inference on singular vectors. $2000$ MC repetitions
were used for a sample size of $500$. : Q-Q plots are theoretical
($y$) vs. empirical ($x$).}
\centering
\begin{tabular}{cc}
\includegraphics[width=8cm,height=6cm,keepaspectratio]{qqplot2_svd} & \includegraphics[width=8cm,height=6cm,keepaspectratio]{cdfwald-svd}\tabularnewline
(a) Wald test quantiles. & (b) CDF comparisons of Wald test statistic.\tabularnewline
\end{tabular}
\justifying
\noindent{\footnotesize\emph{Notes}}{\footnotesize : The plots show
performance of the Wald test, both in terms of quantiles and the CDF.
In panel (a), we see only a few outliers and generally good agreement
between empirical and reference distributions throughout the support
in panel (b).}{\footnotesize\par}
\end{figure}
The Q-Q plot in panel (a) of \prettyref{fig:quality-approx-svd} shows
that the asymptotic approximation for the SVD performs well with few
outliers. Similarly, panel (b) shows the case for $q=m=1$ and that
the empirical CDF tracks the implied $\chi_{1}^{2}$ benchmark well.
These results are encouraging that our first-order expansion of the
SVD map delivers a good approximation for conducting inference on
singular subspace vectors. Details on the DGP appear in \prettyref{alg:wald-and-t-svd}.
\section{Conclusion\label{sec:Conclusion}}
This paper has extended the inferential theory of \citet{Tyler1981}
to cover non-diagonalizable matrices and applied the results to network
statistics. In addition to the Wald test for full vector hypotheses,
a $t$-test is practically useful because it allows inference on individual
coefficients. The method of smooth eigenvector estimation in \citet{10.1093/biomet/asad018}
presents a useful extension to possibly avoid introducing outliers
through the proposed normalization and to take advantage of infinite
differentiability.
Regarding the underlying perturbation theory, there are strong resemblances
between \citet{Sun1991} and \citet{Kato} although the former also
discusses the extensions to Gateaux derivatives. We are hopeful that
the present exposition can aid researchers in similar settings and
present a way to find Jacobians of maps that are related to invariant
or singular subspaces. We refer readers new to the literature on perturbation
theory to \citet{2greenbaum:2019} who offer a pedagogic and detailed
treatment of the subject. We are hopeful that our arguments are easily
adaptable for statistics that depend on matrix subspaces in a more
general way both in the graph domain and others.
The leading application of our results is to invariant subspaces of
estimated network adjacency matrices, which inherit the measurement
error of the estimated network. Similarly, our results on singular
vectors may be applied to network clustering algorithms to allow quantifying
the uncertainty in low-dimensional representations of networks. Generally,
invariant subspaces find applications in the analysis of not only
adjacency matrices but also graph Laplacians where they enable algorithms
for spectral clustering, community detection, or the finding of mixing
rates for random walks on graphs. In this vein, using covariate-level
information as suggested in \citet{BinkiewiczVogelsteinRohe2017}
could be useful to constrain invariant subspaces and potentially tighten
confidence intervals.
Network centralities are often used to identify interventions. However,
point estimates may be misleading if the uncertainty in the network
link identification is large. Therefore, we advocate reporting confidence
bands.
The higher-order Davis-Kahan bounds reveal important distinctions
between general invariant subspace and eigenspace perturbations. In
the case of the former, the eigengap dominates the perturbation estimations
while for the latter, higher-order bounds simply consist of the first-order
Davis-Kahan bound raised to a higher power. The consequence is that,
for eigenspaces, statisticians only need control over the Davis-Kahan
bound while for invariant subspaces, the eigenvalue gap is more important
than the estimation precision.
We have also shown in \prettyref{thm:higher-order-dk} and \prettyref{thm:stoch-bound-prop}
that the eigengap governs the convergence speed of subspace-based
network statistics and explicitly calculated convergence rates of
clustering coefficients whenever the network is estimated using a
single large matrix. A useful extension would be to apply methods
similar to those in \citet{10.1093/biomet/asac032} where a considered
network statistic is the maximum eigenvalue of the adjacency matrix
corresponding to a random eigengap.
The inference results assume fixed matrix size but the convergence
rate and perturbation bound calculations are agnostic about matrix
size. Therefore, an extension of this study is to consider our results
for large, random matrices in more detail, which we leave for future
work.
\bibliographystyle{chicago}
\bibliography{eigenvector-inference,biometrika-one-jasa-references,biometrika_spectral_graph,empirical-applications}