Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
78,452 characters · 23 sections · 36 citation commands
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{#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.}
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 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 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 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 davis1969some bound found in 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 ZhangXuZhu2022. Invariant subspaces are also used to approximate the latent spaces in random dot product graph models 10.1007/978-3-540-77004-6_11 where XieXu2020 propose a method to estimate these spaces and 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 cattuto2010dynamics,krivitsky2014separable,prawesh2019small or health dynamics rothenberg1998social,cornwell2009network,christakis2010social, or biology 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.
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
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 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
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
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
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
\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. 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. 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.
\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 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 vliet and rothenhausler2015backshift correspond to OLS estimation of $M$. Generally, OLS methods work well if $T>>p$, which we explore in \prettyref{subsec:application}.
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
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 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
\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)$.
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 (ref), 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
and
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.
The proof appears in \prettyref{app:distributions}.
Analogously to (ref), 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
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
for which we obtain
\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}.
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.
To make invariant subspace coordinates unique, ensuring any associated estimates are consistent across experiments, we normalise $\upsilon_{\perp}=:
$ 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
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
Consequently,
so the normalized vector is $
$. Observe that \prettyref{eq:ahat-1} defines a unique estimator.
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
in analogy to (ref). 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$
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
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
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
\prettyref{subsec:Network-centralities-as} gives an application of this result.
Our next result characterizes the distribution of $\hat{D}_{I,n}^{\mathsf{T}}$. Define
Then, we obtain
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
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
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
where \[ \hat{\Phi}\left(\upsilon_{\perp}^{\mathsf{T}}\right)\coloneqq
\] 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
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.
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
which forms the basis of all spectral-based statistics. The main result is
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 (ref) 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
\prettyref{app:Additional-details-on-invar} contains a proof. Importantly, (ref) reduces to $C_{w}$ in 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$.
The Jacobian in (ref) reduces to the bound in 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
where $\text{d}\psi$ is the infinitesimal response to perturbation $E$. This bound differs from the original found in davis1969some by the fact that only population eigenvalues appear in the denominator. The relation (ref) 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
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 (ref). 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.
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}.
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}.
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.
The map $X\mapsto(X-\alpha I)^{-1}$ forms the basis of a large class of network centrality measures, which we term 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
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) 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
The bound attained in \prettyref{thm:stoch-bound-prop} resembles the Davis-Kahan theory davis1969some,demetrius except that we derived it from perturbative arguments.
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
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.
Follow jackson2008social and define the node-wise clustering coefficient $\text{cl}_{i}$ via
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}$.
\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.
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 benaychgeorges2018lectureslocalsemicirclelaw provide details on required normalization.
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
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$.
A related family of centrality measures is conceptually very similar to the eigenvector centrality but performs better for dense graphs. Formally, the PageRank centrality brin1998anatomy $c_{P}$ is given by
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. banerjee2013diffusion introduce this idea formally via the diffusion centrality as
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
which we obtain by setting $\beta=\ensuremath{\mathbf{1}}$ in the PageRank centrality (ref). We summarize the relationship between eigenvector, Katz, PageRank, and Diffusion centralities in
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
while for the diffusion centrality we have
The following result allows conducting inference on eigenvector-based centrality measures.
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 magnus2019matrix.
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
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.
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.
\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.
\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.
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.
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.
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.
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}.
This paper has extended the inferential theory of 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 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 Sun1991 and 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 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 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 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.