Eigenvalue and Eigenvector Statistics in Time Series Analysis

TL;DR

This paper develops a supersymmetric theoretical framework to analyze the eigenvalue and eigenvector statistics of cross-correlation matrices in correlated time series, providing a universal benchmark for complex system analysis.

Contribution

It introduces a novel supersymmetric approach to derive analytical results for eigenvector statistics in correlated time series, filling a gap in existing theoretical understanding.

Findings

Analytical expressions for eigenvector statistics derived

Benchmark results for correlated signal analysis established

Framework applicable to diverse complex systems

Abstract

The study of correlated time-series is ubiquitous in statistical analysis, and the matrix decomposition of the cross-correlations between time series is a universal tool to extract the principal patterns of behavior in a wide range of complex systems. Despite this fact, no general result is known for the statistics of eigenvectors of the cross-correlations of correlated time-series. Here we use supersymmetric theory to provide novel analytical results that will serve as a benchmark for the study of correlated signals for a vast community of researchers.

Figures

Equations

$$W=\sqrt{C_L}(W_0+W_1)\sqrt{C_R},$$

$$P(W_1)=(2\pi {\sigma }^2{)}^{-pn/2}\exp \left[-\frac{1}{2{\sigma }^2}\mathrm{tr}{W}_1{W}_1^T\right],\sigma >0.$$

$${\mathcal{I}}_b(\mu |\lambda )=\frac{1}{p\rho (\lambda )}{\sum }_{a=1}^p⟨\delta (\mu -|U_{ab}{|}^2)\delta (\lambda -{\lambda }_a)⟩,$$

$$\rho (\lambda )=\frac{1}{p}{\sum }_{a=1}^p⟨\delta (\lambda -{\lambda }_a)⟩$$

$${\mathcal{I}}_b^{(\mathrm{Haar})}(\mu |\lambda )=\sqrt{\frac{p}{2\pi \mu }}\exp \left[-\frac{p\mu }{2}\right],$$

$$I_{q,b}(\lambda )=⟨{\mu }^q⟩=\frac{1}{p\rho (\lambda )}\sum _{a=1}^p⟨|U_{ab}{|}^{2q}\delta ({\lambda }_a-\lambda )⟩,$$

$$I^{\prime}_{1,b}(\lambda)=-\frac{1}{\sqrt{\lambda_{+}}}\left\langle\left\{\left[\begin{array}[]{cc}\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}&W\\ W^{T}&\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}\end{array}\right]^{-1}\right\}_{bb}\right\rangle,$$

$$Z_{b}(\lambda)=\left\langle\exp\left[i\alpha^{2}e_{b}^{T}\left[\begin{array}[]{cc}\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}&W\\ W^{T}&\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}\end{array}\right]^{-1}e_{b}\right]\right\rangle,$$

$$I_{1,b}(\lambda )=\frac{\int d[T]\exp \left[F(T_{+},T_{-})\right]G_{1b,1b}\sqrt{\mathrm{Sdet}G}}{\sqrt{{\lambda }_{+}}\int d[T]\exp \left[F(T_{+},T_{-})\right]},$$

$$G=\left[\begin{array}[]{cc}\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}\otimes LJ-C_{L}\otimes(T_{+}-iT_{-}+\sigma^{2}L)LJ&\sqrt{C_{L}}W_{0}\sqrt{C_{R}}\otimes LJ\\ \sqrt{C_{R}}W_{0}^{T}\sqrt{C_{L}}\otimes LJ&\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}\otimes LJ-C_{R}\otimes(T_{+}+iT_{-}+\sigma^{2}L)LJ\end{array}\right]^{-1},$$

$$\begin{array}{rl}& \frac{r_0}{{\sigma }^2}=\mathrm{tr}{\left[A_L-W_0{A}_R^{-1}{W}_0^T\right]}^{-1},\\ & \frac{s_0}{{\sigma }^2}=\mathrm{tr}{\left[A_R-W_0^T{A}_L^{-1}{W}_0\right]}^{-1},\mathrm{with}\\ & A_L=\sqrt{\lambda }{C}_L^{-1}-s_0\text{leavevmode}11_p\mathrm{and}{A}_R=\sqrt{\lambda }{C}_R^{-1}-r_0\text{leavevmode}11_n.\end{array}$$

$$\begin{split}\rho(\lambda)\propto&{\rm Im}\left[{\rm tr\,}\left(Q^{-1}{\rm diag\,}(\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p},0)\right)\right],\\ Q=&\left[\begin{array}[]{cc}\sqrt{\lambda}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}-s_{0}C_{L}&\sqrt{C_{L}}W_{0}\sqrt{C_{R}}\\ \sqrt{C_{R}}W_{0}^{T}\sqrt{C_{L}}&\sqrt{\lambda}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}-r_{0}C_{R}\end{array}\right],\end{split}$$

$$I_{1,b}(\lambda )=\frac{\mathrm{Im}\left[\mathrm{tr}\left(Q^{-1}\mathrm{diag}({\hat{E}}_b,0)\right)\right]}{\mathrm{Im}\left[\mathrm{tr}\left(Q^{-1}\mathrm{diag}(\text{leavevmode}11_p,0)\right)\right]},$$

$$\begin{split}\frac{r_{0}}{\sigma^{2}}=&\sum_{j}\frac{1}{\sqrt{\lambda}l_{j}-s_{0}}\\ &+\frac{\sum_{t}x_{t}^{2}\sum_{j}w_{j}^{2}/(\sqrt{\lambda}l_{j}-s_{0})^{2}}{\sqrt{\lambda}-r_{0}-\sum_{t}x_{t}^{2}\sum_{j}w_{j}^{2}/(\sqrt{\lambda}l_{j}-s_{0})},\\ \frac{s_{0}}{\sigma^{2}}=&\frac{1}{\sqrt{\lambda}-r_{0}}\biggl{(}n\\ &+\frac{\sum_{t}x_{t}^{2}\sum_{j}w_{j}^{2}/(\sqrt{\lambda}l_{j}-s_{0})}{\sqrt{\lambda}-r_{0}-\sum_{t}x_{t}^{2}\sum_{j}w_{j}^{2}/(\sqrt{\lambda}l_{j}-s_{0})}\biggl{)}.\end{split}$$

$$\begin{array}{rl}(Q^{-1}{)}_{ii}=& \frac{l_i}{\sqrt{\lambda }{l}_i-s_0}+\frac{1}{(\sqrt{\lambda }{l}_i-s_0{)}^2}\\ & \times \frac{l_i{w}_i^2{\sum }_t{x}_t^2}{\sqrt{\lambda }-r_0-{\sum }_t{x}_t^2{\sum }_j{w}_j^2/(\sqrt{\lambda }{l}_j-s_0)}.\end{array}$$

$$I_{q,b}(\lambda )=\frac{(2q)!}{2^{q}q!}[I_{1,b}(\lambda ){]}^q,$$

$${\mathcal{I}}_b(\mu |\lambda )=\frac{1}{\sqrt{2\pi I_{1,b}(\lambda )\mu }}\exp \left[-\frac{\mu }{2I_{1,b}(\lambda )}\right].$$

$$Z_{b}(\lambda)=\left\langle\exp\left[i\alpha^{2}e_{b}^{T}\left[\begin{array}[]{cc}\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}&W\\ W^{T}&\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}\end{array}\right]^{-1}e_{b}\right]\right\rangle$$

$$\begin{split}&\exp\left[i\alpha^{2}e_{b}^{T}\left[\begin{array}[]{cc}\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}&W\\ W^{T}&\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}\end{array}\right]^{-1}e_{b}\right]\int d[\psi]d[\phi]\exp\left[-{\rm Str\,}(\phi,\psi)^{T}(\phi,\psi)J\right]\\ =&\int d[\psi]d[\phi]\exp\left[i{\rm Str\,}(\phi,\psi)^{T}\left[\begin{array}[]{cc}\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}&W\\ W^{T}&\sqrt{\lambda_{+}}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}\end{array}\right](\phi,\psi)J+2\alpha{\rm Str\,}(\phi,\psi)^{T}(e_{b},0,0,0)\right],\end{split}$$

$$(\phi,\psi)=\left(\begin{array}[]{c}V_{L}\\ V_{R}\end{array}\right).$$

$$\begin{array}{rl}{Z}_b(\lambda )& =\frac{⟨\int d[V_R,V_L]\exp \left[S(V_R,V_L,W)\right]⟩}{\int d[V_R,V_L]\exp \left[S_0(V_R,V_L)\right]},\\ S(V_R,V_L,W)& =-i\sqrt{{\lambda }_{+}}\mathrm{Str}LJ(V_L^T{V}_L+V_R^T{V}_R)-i\mathrm{Str}LJ(V_L^{T}W{V}_R+V_R^T{W}^T{V}_L)+2\mathrm{Str}(V_L,V_R{)}^T{E}_b,\\ S_0(V_R,V_L)& =-\mathrm{Str}J(V_L^T{V}_L+V_R^T{V}_R).\end{array}$$

$$\begin{array}{rl}{Z}_b(\lambda )=& \frac{\int d[V_R,V_L]\exp \left[S_1(V_R,V_L)+S_2(V_R,V_L)\right]}{\int d[V_R,V_L]\exp \left[S_0(V_R,V_L)\right]},\\ S_1(V_R,V_L)=& -i\sqrt{{\lambda }_{+}}\mathrm{Str}LJ(V_L^T{V}_L+V_R^T{V}_R)-i\mathrm{Str}LJ(V_L^T\sqrt{C_L}{W}_0\sqrt{C_R}{V}_R+V_R^T\sqrt{C_R}{W}_0^T\sqrt{C_L}{V}_L)\\ & +2\mathrm{Str}(V_L,V_R{)}^T{E}_b,\\ S_2(V_R,V_L)=& -2{\sigma }^2\mathrm{Str}{V}_L^T{C}_L{V}_{L}LJ{V}_R^T{C}_R{V}_{R}LJ.\end{array}$$

$$\begin{array}{rl}{Z}_b(\lambda )\propto & \int d[V_R,V_L]\int d[T]\exp \left[S_1(V_R,V_L)+S_3(V_R,V_L,T_{+},T_{-})\right],\\ S_3(V_R,V_L,T_{+},T_{-})=& i{\sigma }^2\mathrm{Str}J(V_L^T{C}_L{V}_L+V_R^T{C}_R{V}_R)-\frac{1}{2{\sigma }^2}\mathrm{Str}(T_{+}^2+T_{-}^2)\\ & +i\mathrm{Str}{T}_{+}JL(V_L^T{C}_L{V}_L+V_R^T{C}_R{V}_R+iJ)+\mathrm{Str}{T}_{-}JL(V_L^T{C}_L{V}_L-V_R^T{C}_R{V}_R).\end{array}$$

$$T_{+}=\left[\begin{array}[]{cc}B_{1}+iC(B_{2})&\eta_{1}\hat{\tau}_{2}\\ -\eta_{1}^{T}&iF_{1}\end{array}\right]\ {\rm and}\ T_{-}=\left[\begin{array}[]{cc}B_{2}&\eta_{2}\hat{\tau}_{2}\\ -\eta_{2}^{T}&iF_{2}\end{array}\right]$$

$$d[T]=d[B_1]d[B_2]d[F_1]d[F_2]d[{\eta }_1]d[{\eta }_2].$$

$$B_1=-O_1{b}_1{O}_1^{-1}\mathrm{and}{B}_2=-O_2{b}_2{O}_2^{-1},$$

$$C(B_2)=-O_2\sqrt{\text{leavevmode}11_2+b_2^2}{O}_2^{-1}.$$

$$Z_b(\lambda )\propto \int d[T]\exp \left[-\frac{1}{2{\sigma }^2}\mathrm{Str}(T_{+}^2+T_{-}^2)-\mathrm{Str}{T}_{+}L-i{\alpha }^2{G}_{1b,1b}\right]\sqrt{\mathrm{Sdet}G},$$

$$\begin{array}{rl}{Z}_b(\lambda )=& \frac{\int d[T]\exp \left[-\frac{1}{2{\sigma }^2}\mathrm{Str}(T_{+}^2+T_{-}^2)-\mathrm{Str}{T}_{+}L-i{\alpha }^2{G}_{1b,1b}\right]\sqrt{\mathrm{Sdet}G}}{\int d[T]\exp \left[-\frac{1}{2{\sigma }^2}\mathrm{Str}(T_{+}^2+T_{-}^2)-\mathrm{Str}{T}_{+}L\right]}.\end{array}$$

$$I_{1,b}(\lambda )=\frac{1}{\sqrt{{\lambda }_{+}}}\frac{\int d[T]\exp \left[-\frac{1}{2{\sigma }^2}\mathrm{Str}(T_{+}^2+T_{-}^2)-\mathrm{Str}{T}_{+}L\right]G_{1b,1b}\sqrt{\mathrm{Sdet}G}}{\int d[T]\exp \left[-\frac{1}{2{\sigma }^2}\mathrm{Str}(T_{+}^2+T_{-}^2)-\mathrm{Str}{T}_{+}L\right]},$$

Full text

Introduction

The theory of complex systems ultimately deals with the identification of patterns of simple behaviours accounting for the emergence of universal dynamics in the time series measured in a vast range of disciplines, including condensed matter physics, medicine, finance, signal transmission, biology, and more recently computational social sciences Strogatz (2018). A time series is a series of values scanned over time of a given observable of a system Chatfield (2018) such as the sea level Rahmstorf (2007), the temperature of a lake Sharma et al. (2015), the neuron activity in electroencephalography (EEG) Šeba (2003); Müller et al. (2006), the response in a unit of volume of a magnetic resonance imaging experiment Kwong et al. (1992), the gross domestic product of a country Lee (2005), the price or return of a stock Plerou et al. (1999); Barucca (2014), the volume of an order in the market Bouchaud et al. (2009); Chiarella et al. (2009), the infected individuals in a region affected by an epidemics Grenfell et al. (2001), and the online activity of a user O’Connor et al. (2010).

The basic analysis that is ubiquitously performed when dealing with multiple time series are covariance and correlation analysis, especially with the aim of identifying the main factors accounting for time variability and parsimoniously representing the state space of the system, through denoising and dimensionality reduction. The generality of this statistical approach constitutes the basis for Principal Component Analysis (PCA) Pearson (1901); Jolliffe (2011), clustering analysis, and many other data mining algorithms Lloyd (1982). In these techniques one distinguishes between eigenvalue and eigenvector statistics and both of them carry important information as we know, for example, from the theory of quantum disordered systems. Therefore it is even more surprising that only few results are available for the cross-statistics between eigenvalues and eigenvectors, when dealing with the covariance and correlation matrices of noisy time series.

The spectral density of the eigenvalues is up to now the major quantity where the theory provides robust and general results Laloux et al. (1999); Lillo and Mantegna (2005); Allez et al. (2012); Majumdar and Vivo (2012). For instance, the Marc̆enko–Pastur distribution (MPD) Marc̆enko and Pastur (1967) usually serves as a blueprint for describing the influence of white noise in the time series on the spectral density. Any deviation from the MPD, for instance outliers, can be considered as system specific information so that the MPD serves as a filter. However, some eigenvalues encoding relevant information might be obscured by the bulk of the spectrum described by the MPD. Then PCA may remove relevant data that should be taken into account. To distinguish those system specific eigenvalues from the eigenvalues of the MPD one needs to take into consideration the eigenvector statistics. An important step in this direction is made in the present Letter. We derive an analytical formula for the first moment of a fixed eigenvector component conditioned to a chosen eigenvalue. Moreover, we state a conjecture on their general moments and distributions for a correlation matrix of noisy time series. Our results provide insights and pave the way for a much more informative spectral decomposition in time series analysis, allowing not only to focus on the spectral density but also on the individual contribution of each component to the spectrum, leading to a much deeper understanding of a system’s dynamics.

Random Matrix Model

Specifically, we study the statistics of the eigenvectors and the eigenvalues of the matrix $C=WW^{T}$, with $W\in\mathbb{R}^{p\times n}$ representing $p$ time series of length $n$ or, in the case of PCA, $p$ descriptors with $n$ variants, and $W^{T}$ being the transpose of $W$. Thus $C$ can be interpreted as the covariance matrix between the time series aggregated in $W$ or the covariance between the descriptors respectively. The real rectangular matrix $W$ in our model is composed of four matrices

[TABLE]

where $W_{0}\in\mathbb{R}^{p\times n}$ is a deterministic real matrix and $W_{1}\in\mathbb{R}^{p\times n}$ is a Gaussian random matrix distributed by

[TABLE]

The two real symmetric matrices $C_{L}=C_{L}^{T}\in\mathbb{R}^{p\times p}$ and $C_{R}=C_{R}^{T}\in\mathbb{R}^{n\times n}$ are positive definite and represent a spatio-temporal correlation between the various time series. Here, the matrix $C_{L}$ can be identified with a time correlation, the matrix $C_{R}$ with the spatial correlations, and additionally, at difference with many common models, we include an offset $W_{0}$. Hence $W$ is a non-centred and doubly correlated Gaussian random matrix. This form allows the model to capture in detail the case of factor models ubiquitous in statistics and econometrics.

Though our model is quite general, it is still not the most general Gaussian random matrix model. We assume that the spatio-temporal correlations of the multivariate time series factorize in the two matrices, $C_{L}$ and $C_{R}$. Therefore time-dependent spatial correlations, like the two epoch model Akemann et al. (2016), are not considered here.

The random matrix model defined above can be also considered as a simple deformation of the standard real Wishart ensemble of random matrices, in which the orthogonal invariance is broken in several ways. Such non-invariant deformations of the standard random matrix ensembles were introduced and studied in different contexts including wireless communication Couillet and Debbah (2011), vibration analysis Soize (2003), signal processing Nadakuditi and Edelman (2008) and neural networks Ahmadian et al. (2015). There is a growing interest to the statistical properties of the eigenvectors in these ensembles. While there are some recent results about the statistics of the eigenvectors in the deformed Gaussian Orthogonal and Unitary ensembles Allez et al. (2014); Truong and Ossipov (2016a, b); Bun et al. (2017); Bourgade and Yau (2017); Benigni (2017); Truong and Ossipov (2018), we are not aware of similar results for the Wishart ensemble except for Ref.Bourgade and Yau (2017), in which the ergodicity of the eigenvectors was proven for the special case $C_{L}=\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}$, $C_{R}=\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}$.

In the following, we will not simply focus on the computation of the spectral density of the eigenvalues, analysed in Recher et al. (2010, 2012) with the same supersymmetric (SUSY) approach as in the present work, but also calculate a detailed eigenvector statistics of the matrix $WW^{T}=U\Lambda U^{T}$, whose eigenvalues represented by the diagonal matrix $\Lambda={\rm diag\,}(\lambda_{1},\ldots,\lambda_{p})$ and the eigenvectors represented by the columns of the matrix $U=\{U_{ab}\}\in\mathrm{O}(p)$. The full information about the statistics of the eigenvector components is contained in the conditional density

[TABLE]

where $b=1,\dots,p$ refers to a particular eigenvector component and

[TABLE]

is the mean density of the eigenvalues and $\langle.\rangle$ stands for the ensemble average over the distribution of $W_{1}$. In the case of a factorisation of the eigenvector and eigenvalue statistics, as in the Wishart ensemble, one finds the Porter–Thomas distribution Porter and Thomas (1956)

[TABLE]

which is independent of the component $b$ and the eigenvalue $\lambda$ due to the Haar distributed eigenvectors. This simplification cannot be expected to hold in our non-trivial model as well as in a realistic situation. The computation of (3) or its arbitrary moments

[TABLE]

where $q$ is a positive integer, is technically a very challenging problem. In this Letter we focus on the analytical derivation of the first moment $I_{1,b}(\lambda)$ and make a conjecture about an arbitrary moment $I_{q>1,b}(\lambda)$ and $\mathcal{I}_{b}(\mu|\lambda)$ in the conclusions.

The moments of the eigenvectors are also a standard tool to characterise properties of complex quantum systems and are used to distinguish different phases in condensed matter physics Evers and Mirlin (2008). Hence, we expect that it may give valuable insights for time series as well.

Before we start with the analytical calculation of $I_{1,b}$, we want to point out that the eigenvector components $U_{ab}$ are basis dependent. Thus the conditional distribution $\mathcal{I}_{b}(\mu|\lambda)$ strongly depends on the reference frame. In this work such a frame is chosen as the eigenbasis of $C_{L}$, allowing us to investigate the broadening of the eigenvectors due to the white noise $W_{1}$ and its strength $\sigma$. Another natural and valuable reference frame could be the eigenbasis of $\sqrt{C_{L}}W_{0}C_{R}W_{0}^{T}\sqrt{C_{L}}$ which we do not consider here for simplicity.

Eigenvector Statistics with SUSY

The first moment of the eigenvectors, see (6) for $q=1$, can be computed by taking the imaginary part and the limit of a regularization $\epsilon\to 0$ of the quantity

[TABLE]

where $\lambda_{+}=\lambda+i\varepsilon$. Defining the $(p+n)$-dimensional unit vector $e_{b}$ with unity at the position $b$ and zero otherwise, this quantity can be generated by differentiating

[TABLE]

with respect to the auxiliary parameter $\alpha$, at $\alpha=0$. $\alpha$ is chosen to be real to guarantee convergence later on. Following the standard steps of the SUSY method Recher et al. (2010, 2012), we represent first the generating function $Z_{b}(\lambda)$ by the supersymmetric Gaussian integral, average over the random matrix $W_{1}$ and finally apply the Hubbard-Stratonovich transformation Fyodorov et al. (2008). In this way, we derive the following supersymmetric representation for $I_{1,b}(\lambda)$ (see the Supplemental Material sup for details),

[TABLE]

where $F(T_{+},T_{-})=-\frac{1}{2\sigma^{2}}{\rm Str\,}(T_{+}^{2}+T_{-}^{2})-{\rm Str\,}T_{+}L$, $L={\rm diag\,}(-\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2};\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2})$. The $(2|2)\times(2|2)$ supermatrices $T_{\pm}$ are symmetric in the boson-boson block and self-dual in the fermion-fermion block and their eigenvalues run along complex contours that are detailed in the Supplemental Material sup . The supersymmetric Green function $G$ has the form

[TABLE]

with $J={\rm diag\,}(\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2};\tau_{2})$. The representation (9) is exact, but rather involved and technical. An expression for the mean level density can be obtained by summing over $b=1,\ldots,p$ and should be compared with the corresponding result in Recher et al. (2010, 2012); Waltner et al. (2015). The above expression simplifies a lot in the limit $n,p\to\infty$, which is considered next.

Macroscopic level density and limiting statistics

In most applications, one is interested in the limit $n,p\to\infty$. In this limit the integral in Eq.(9) can be evaluated using the saddle-point approximation. To derive the saddle-point equation, it is convenient to introduce the supermatrices $S=T_{+}-iT_{-}+\sigma^{2}L$ and $R=T_{+}+iT_{-}+\sigma^{2}L$, which can be considered as independent. The saddle-point solution contributing most to the integral is given by the diagonal matrices $S_{0}=s_{0}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2|2}$ and $R_{0}=r_{0}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2|2}$ with the complex parameters $s_{0}$ and $r_{0}$ that satisfy the coupled equations sup

[TABLE]

The mean level density is up to a normalisation constant given by

[TABLE]

where we assume $p\leq n$ without loss of generality. The case $p>n$ only yields an additional Dirac delta function at the origin. The formula (12) reduces to the MPD Marc̆enko and Pastur (1967) in the case of the Wishart ensemble, i.e., $C_{L}=\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{p}$, $C_{R}=\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}$ and $W_{0}=0$. We illustrate the result for $\rho(\lambda)$ in Fig. 1 for the one-factor model, which is described in the next subsection.

The result for $I_{1,b}(\lambda)$ can be expressed in terms of the same matrix $Q$ and reads

[TABLE]

which constitutes the main result of the present Letter. The normalisation is fixed by the condition $\sum_{b=1}^{p}I_{1,b}(\lambda)=1$. We note that for a Haar distributed vector one has $I_{1,b}^{\rm(Haar)}(\lambda)=1/p$.

One-factor model

To illustrate our findings we apply our general results to the one-factor model supplemented with Gaussian noise. Specifically, we set $W_{0}=wx^{T}$, where $w$ and $x$ are column vectors of length $p$ and $n$, respectively. The correlation matrices are chosen to be diagonal $C_{L}={\rm diag\,}(l_{1}^{-1},\ldots,l_{p}^{-1})$ and $C_{R}=\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n}$. The vector $x$ represents a common factor, e.g. the market mode in financial time series analysis, and the component $w_{j}$ quantifies the relative weight of the common factor on the $j$th time series, before normalization.

We plug the matrices of the one-factor model into the saddle-point equation (11) and simplify the resulting expression via the Sherman-Morrison identity for the inverse matrices Sherman and Morrison (1950), i.e. $(A+uv^{T})^{-1}=A^{-1}-\frac{A^{-1}uv^{T}A^{-1}}{1+v^{T}A^{-1}u}$. This leads to the coupled equations

[TABLE]

Solving these saddle-point equations we can derive the spectral density and the moments of the eigenvectors, simply by plugging the following matrix elements in Eqs. (12)-(13)

[TABLE]

We illustrate these results in Figs. 1 and 2, where we also compare them with Monte-Carlo simulations. The deviations from the Porter-Thomas distribution (5), which yields for the first moment the constant $I_{1,b}^{\rm(Haar)}=1/p$, can be readily seen for some components of the eigenvectors. They indicate that the corresponding eigenvalues still carry a lot of information on the matrix $C_{L}$, although these eigenvalues are evidently inside the bulk of the spectrum, cf., Fig 1. This simple example demonstrates the strength of the combined statistics of eigenvalues and eigenvectors.

Conclusions

The general result in Eq. (13) provides a powerful analytical methodology to quantify the expected value of the square of specific components in a given eigenvalue interval for a wide range of random matrices. We tested numerically these analytical results in detail for the one-factor model (see Figs. 1-2). Our general formulation allows an arbitrary number of factors to be added in the matrix $W_{0}$. Although our analytical results were derived in the limit $n,\>p\to\infty$, they show a very good agreement with the results of numerical simulations at finite $n$ and $p$. The rate of convergence to the limiting statistics will generally depend on the input $W_{0}$, $C_{L}$, and $C_{R}$.

In the present work we derived analytically a closed result only for the first moment $I_{1,b}(\lambda)$ of an eigenvector under the condition of a fixed eigenvalue. However we conjecture that all higher moments are related to the first moment as follows:

[TABLE]

which corresponds to a locally rescaled Porter-Thomas distribution

[TABLE]

A similar result has been also found for the conditioned eigenvector statistics of the deformed Gaussian Unitary Ensemble (GUE) in Truong and Ossipov (2016a, b). The only difference is the prefactor in (16), which is equal to $(2q)!/(2^{q}\,q!)$ in our case and given by $q!$ for the complex eigenvectors in the deformed GUE Truong and Ossipov (2016a, b). These numerical values result from the averaged moments of real and complex normalized vectors, respectively. We have tested this conjecture numerically for $q=2$ and found a nice agreement, see Fig. 3.

We are confident that our analytical results are of general relevance for the spectral decomposition of time series and could lead to unprecedented understanding of the full statistics of the eigen-components in signal analysis. A strong deviation of the moment $I_{q,b}(\lambda)$ from the constant $(2q)!/(2^{q}\,q!)$ hints at an eigenvector-eigenvalue pair that contains system specific information. This knowledge can improve PCA and other techniques to reduce highly dimensional data without loosing relevant information.

Acknowledgements.

MK acknowledges financial support by the German research council (DFG) through CRC 1283: “Taming uncertainty and profiting from randomness and low regularity in analysis, stochastics and their applications”. PB acknowledges support from the London Institute for Mathematical Sciences (LIMS).

I Derivation of the supersymmetric integral representation for the moments of the eigenvectors

The quantity $I^{\prime}_{1,b}$ defined in Eq.(7) can be computed by differentiating the generating function (8)

[TABLE]

with respect to $i\alpha^{2}$ and setting $\alpha=0$. The normalization is given as $\lim_{\lambda\to\infty}Z_{b}(\lambda)=1$. In order to construct a representation of $Z_{b}(\lambda)$ in terms of the supersymmetric integral we use the identity

[TABLE]

where we employed the matrix $\psi$ which is an $(n+p)\times 2$ dimensional matrix of real Grassmann variables and $\phi$ is an $(n+p)\times 2$ dimensional ordinary real matrix. The two matrices are introduced in order to cancel the resulting determinants from the Gaussian integral. To ensure integrability we have introduced the constant matrix $J={\rm diag\,}(\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2};\tau_{2})$, where $\tau_{2}$ is the second Pauli matrix.

To simplify the notation, we define the diagonal $(2|2)\times(2|2)$ supermatrix $L={\rm diag\,}(-\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2};\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2})$ and the $(n+p)\times(2|2)$ rectangular supermatrix $E_{b}=(\alpha e_{b},0;0,0)$. Moreover, we rearrange the matrices $\psi$ and $\phi$ in the $p\times(2|2)$ supermatrix $V_{L}$ and the $n\times(2|2)$ supermatrix $V_{R}$ as follows

[TABLE]

Both matrices are two real rectangular supermatrices $V_{L}=V_{L}^{*}$ and $V_{R}=V_{R}^{*}$ with dimensions $p\times(2|2)$ and $n\times(2|2)$ respectively. The first two columns of $V_{L}$ and $V_{R}$ are real variables while the last two columns are Grassmann variables. In this way, we find

[TABLE]

The average over $W_{1}$ yields

[TABLE]

Since the action contains a quartic term in the matrices $V_{L}$ and $V_{R}$, the next step is to perform the Hubbard-Stratonovich transformation, which allows one to decouple such terms. Up to the normalization the result reads

[TABLE]

The parametrization of the two $(2|2)\times(2|2)$ supermatrices $T_{\pm}$ needs to be chosen carefully to guarantee the convergence of the integral. They are given by

[TABLE]

equipped with the flat Berezinian measure

[TABLE]

The ordinary matrices $B_{1}$ and $B_{2}$ are negative definite and symmetric and can be diagonalized with orthogonal matrices $O_{1},O_{2}\in{\rm O}(2)$ as follows

[TABLE]

with $b_{1},b_{2}$ two positive definite diagonal matrices. The matrix $C(B_{2})$ has the form

[TABLE]

The matrices $F_{1}$ and $F_{2}$ are Hermitian self-dual matrices and $\eta_{1}$ and $\eta_{2}$ are two $2\times 2$ rectangular matrices whose entries are independent real Grassmann variables.

The shift of $B_{1}$ in $T_{+}$ by the imaginary part $\sqrt{\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2}+B_{2}^{2}}$ solves a convergence problem in the Gaussian terms in (S6). In particular the Gaussian integrals over the supermatrices $V_{L}$ and $V_{R}$ are absolutely convergent and yield

[TABLE]

where $G$ is defined as in (10). Hence, $G_{\mu a,\nu b}$ has four indices with $\mu,\nu=1,\ldots,4$ and $a,b=1,\ldots,n+p$. To fix the normalization we take $\lambda\to\infty$ and notice that $G$ becomes approximately $\lambda_{+}^{-1/2}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{n+p}\otimes LJ$. Therefore we end up with the intermediate result

[TABLE]

Coming back to our original problem we notice that we are interested in the first derivative with respect to $i\alpha^{2}$ at $\alpha=0$. In particular, the quantity $I_{1,b}(\lambda)$ is given by

[TABLE]

which coincides with Eq.(9).

II Saddle-point equation

For deriving the saddle-point equation we only need to consider the exponential function and the superdeterminant in the integral (S13). The term $G_{1b,1b}$ is only a polynomial prefactor which does not influence the saddle-point solution. It is easier to study the saddle-point by introducing the supermatrices $S=T_{+}-iT_{-}+\sigma^{2}L$ and $R=T_{+}+iT_{-}+\sigma^{2}L$, which can be considered to be independent. Then the action, i.e. the function that need to be minimised, is

[TABLE]

Differentiating it with respect to $S$ and $R$ yields two coupled equations

[TABLE]

The operator ${\rm tr\,}_{1}$ is the partial trace over the first tensor space which is here the space of ordinary $n\times n$ and $p\times p$ matrices, respectively.

The saddle-point equation is rotation invariant, i.e., when $(S_{0},R_{0})$ is a solution then this is also true for $(R_{0}S_{0}R_{0}^{-1},R_{0})$ as well as $(S_{0},S_{0}R_{0}S_{0}^{-1})$ and any kind of combination. This can be seen by multiplying both equations from the left and the right with $R$ and $R^{-1}$, which is equivalent to replacing $S$ by $RSR^{-1}$. Assuming that the saddle-point solution $(S_{0},R_{0})$ is unique, we conclude then that $S_{0}$ and $R_{0}$ must commute. The uniqueness of the solution should follow from the fact the contour of integration, which was shifted by the term $i\epsilon$, can’t cross the poles and the fact that the Berezinian (the Jacobian in superspace), that is $|b_{1j}-b_{2j}|/[(b_{1j}-if_{j})^{2}(b_{1j}-if_{j})^{2}]$ for $j=1,2$, is not suppressed only when the multiplicity of the eigenvalues in the Fermion-Fermion blocks is equal to those in the Boson-Boson block. The Fermion-Fermion blocks are doubly degenerate due to their Hermitian self-duality. Thus also the Boson-Boson blocks are doubly degenerate, which implies for $(2|2)\times(2|2)$ supermatrices that we can diagonalize $S$ and $R$ simultaneously and the solution has to be diagonal and degenerate, i.e., $S_{0}=s_{0}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2|2}$ and $R_{0}=r_{0}\leavevmode\hbox{\small 1\kern-3.8pt\normalsize 1}_{2|2}$. Substituting this ansatz into Eq. (S15) we derive Eq.(11), which is

[TABLE]

The $\epsilon$ regularization only determines which saddle-point has to be chosen, especially which sign the imaginary part carries. Assuming the correct sign of the imaginary part we neglected this regularization in Eq. (11).