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.
 A: 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:


*

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

*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.

*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$.  

*Eq. (6) is just to avoid numerical instability.

*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)  

