# How to estimate the largest eigenvalue of a correlation matrix from one observation of underlying data matrix?

Suppose that I have $$N$$ time series $$x_{1t},x_{2t},\dots,x_{Nt},$$, that are correlated with each other. A $$N\times N$$ correlation matrix is $$R=\rho_{ij}$$. It can be represented with eigen value dcomposition: $$R=\sum_{j=1}^N\lambda_j\xi_j\xi_j'$$ Assuming the eigenvalues are sorted in descending order $$\lambda_1$$ is the largest eigenvalue. It represents the overall correlation of series. My question: How to estimate $$\lambda_1$$ when you have only one observation at time $$t$$ of each series $$x_{it}$$? We have only $$N$$ values $$x_{it}$$ where $$i\in[1,N]$$.

If it helps, I suspect that the series might have $$K$$ common factors, so that $$x_{it}=\sum_{j=1}^KB_{ij}f_{jt}+\varepsilon_i,$$ where the noise is independent, i.e. $$corr[\varepsilon_i,\varepsilon_j]=\delta_{ij}$$ and $$corr[\varepsilon_i,f_j]=0$$.

Me thinks that maybe when $$K=1$$, we can somehow estimate $$\hat\lambda_1$$. However, I'd like a more general solution if possible.

# Idea

Eigen value $$\lambda_1$$ represents overall correlation. Let's write a covariance of a sum of series: $$cov(\sum_i x_i)=\sum_{ij}\sigma_i\sigma_j\rho_{ij}=\sum_i\sigma^2+\sum_{i\ne j}\sigma_i\sigma_j\rho_{ij}$$

Assuming $$\rho_{ij}=\rho$$ we get the following metric: $$r=\frac{(\sum_ix_i)^2-\sum_ix_i^2}{\sum_{i\ne j}x_ix_j}$$

Suppose that one day I have two series observations (1,3), then another day they are (5,-1). Here's how we calculate the metrics: $$r_1=\frac{4^2-(1^2+3^2)}{2\times1\times 3}=\frac{6}{6}=1$$ $$r_2=\frac{4^2-(5^2+1^2)}{2\times5\times (-1)}=\frac{-10}{-10}=-1$$

In the first day both series moved in the same direction, and although the moves were smaller than in the second day, they still produced the same average move 2 because they were more correlated. Hence, this metric seems to be related to $$\lambda_1$$.

The idea's inspired by Eqs.6 and 7 in Bouchaud's paper.

• Would you be wanting to deconvolve the different temporal frequencies across the N series (wavelet or Fourier transform?), then rank them by variance? – ReneBt Feb 24 at 20:36

I suggest you a very simple online update solution:

At each time $$t$$ you provide me a point of a time series $$x$$ of size $$N$$, then

(1) $$\hat\mu=\hat\mu+\frac{1}{1+t}(x-\hat\mu)$$

(2) $$\hat\Sigma=\frac{t}{1+t}\hat\Sigma+\frac{t}{(1+t)^2}(x-\hat\mu)^\prime (x-\hat\mu)$$

(3) $$v^{new}=\hat\Sigma v$$

(4) $$\hat\lambda=\frac{1}{N}\sum_{i=1}^{N} v_{i}^{new}/v_{i}$$

(5) $$v=v^{new}$$

(6) $$v=v/||v||$$

Explanation:

1. Eq. (1) and (2) are basic online update formula for the first and second moment. See for instance this paper.

2. For understanding Eq. (3), note que $$Av=\lambda v$$. Thus, the right side is original eigenvector multiplied by $$\lambda$$. If $$A$$ is a constant matrix, one can prove that it converges to the eigenvector associated with the largest eigenvalue.

3. It is just an estimation of lambda. I average the division of each element of $$\lambda v$$ by each element of $$v$$, which is an estimation of $$\lambda$$.
4. Eq. (6) is just to avoid numerical instability.
5. Steps (3), (4), (5) and (6) are based on the basic algorithm to evaluate the numerical value of the largest eigenvalue of a matrix. See for instance Scientific Computing: An Introductory Survey - Michael T. Heath.

See the code below to note that it converges for the desired value. In this code I create a time series that come from a multinormal distribution with known mean and known covariance matrix with known eigenvalues and eigenvectors.

import numpy as np
from scipy.linalg import norm
def rotationMatrix(theta):
return np.array([[np.cos(theta), -np.sin(theta)],[np.sin(theta), np.cos(theta)]])

def generateCovarianceMatrix(theta,lambda1,lambda2):
P=rotationMatrix(theta)
eigMatrix=np.array([[lambda1, 0],[0, lambda2]])
Sigma=P.T.dot(eigMatrix.dot(P))
return Sigma

def generateData(mu,sigma):
return np.random.multivariate_normal(mu, sigma, 1)

if __name__=="__main__":
sigma=generateCovarianceMatrix(np.pi/6,2.0,7.0)
mu = [1, 2]
T=1000
hatMu=[0,0]
hatSigma=np.array([[1, 0],[0, 1]])
hatLargeVector=[1,1] # eigenvector associated with the largest eigenvalue
for i in range(T):
x=generateData(mu,sigma)
hatMu=hatMu+(1/(i+1))*(x-hatMu)
if(i==0):
hatSigma=np.array([[1, 0],[0, 1]])
else:
hatSigma=(i/(i+1))*hatSigma+i*(1/(i+1))*(1/(i+1))*(x-hatMu).T.dot((x-hatMu))
oldHatLargeVector=hatLargeVector
hatLargeVector=hatSigma.dot(hatLargeVector)

print('x',x)
print("mu",hatMu)
print("sigma",hatSigma)
print("v",hatLargeVector)
epsilon=0.001
hatLargestLambda=np.sum([(hatLargeVector[i]+epsilon)/(oldHatLargeVector[i]+epsilon) for i in range(len(hatLargeVector))])/len(hatLargeVector)
print(hatLargestLambda)
hatLargeVector=hatLargeVector/norm(hatLargeVector)

• @Aksakal does this solution not work for you? It really works. – DanielTheRocketMan Mar 3 at 0:28