Eigenvectors of a small norm adjustment

I have a dataset that is slowly changing, and I need to keep track of eigenvectors/eigenvalues of its covariance matrix.

I've been using scipy.linalg.eigh, but it's too expensive, and it doesn't use the fact that I already have a decomposition that's only slightly incorrect.

Can anyone suggest a better approach to deal with this problem?

• How large is your data? Do you need the complete eigensystem, or only some of the largest eigenvalues? Do you need them exactly, or would an approximation do?
– cfh
Apr 24 '17 at 10:17
• I need complete eigensystem. I found an algorithm for updating inverse of a matrix after small norm update using regression interpretation of inverse of covariance matrix, so I assumed something similar should exist for eigenvectors. Apr 24 '17 at 16:11
• What do you do with that full eigendecomposition? There might be a better shortcut that does not go through it... And I reiterate cfh's question: "how large"? Apr 29 '17 at 10:56
• I have 8k features and millions of datapoints, so covariance is approximate. This is to implement this algorithm. Gradient update depends on eigenvalues of certain covariance matrix, and this covariance matrix changes at each step Apr 29 '17 at 17:19

A naive approach is to use the eigenvalue solution of your matrix $A(t)$ as the initial guess of an iterative eigensolver for matrix $A(t + \delta t)$. You might use QR if you need the full spectrum, or the power method otherwise. This isn't an entirely robust approach, however, as the eigenvalues of a matrix aren't necessarily close to a nearly neighboring matrix (1), especially if its poorly conditioned (2).

A subspace tracking method is apparently more useful (3). An excerpt from (4):

The iterative computation of an extreme (maximal or minimum) eigen pair (eigenvalue and eigenvector) can date back to 1966 . In 1980, Thompson proposed a LMS-type adaptive algorithm for estimating eigenvector, which correspond to the smallest eigenvalue of sample covariance matrix, and provided the adaptive tracking algorithm of the angle/frequency combing with Pisarenko’s harmonic estimator . Sarkar et al.  used the conjugate gradient algorithm to track the variation of the extreme eigenvector which corresponds to the smallest eigenvalue of the covariance matrix of the slowly changing signal and proved its much faster convergence than Thompson’s LMS-type algorithm. These methods were only used to track single extreme value and eigenvector with limited application, but later they were extended for the eigen-subspace tracking and updating methods. In 1990, Comon and Golub  proposed the Lanczos method for tracking the extreme singular value and singular vector, which is a common method designed originally for determining some big and sparse symmetrical eigen problem $Ax = kx$ .

: Comon, P., & Golub, G. H. (1990). Tracking a few extreme singular values and vectors in signal processing. In Processing of the IEEE (pp. 1327–1343).

: Thompson, P. A. (1980). An adaptive spectral analysis technique for unbiased frequency

: Bradbury, W. W., & Fletcher, R. (1966). New iterative methods for solutions of the eigenproblem. Numerical Mathematics, 9(9), 259–266.

: Sarkar, T. K., Dianat, S. A., Chen, H., & Brule, J. D. (1986). Adaptive spectral estimation by the conjugate gradient method. IEEE Transactions on Acoustic, Speech, and Signal Processing, 34(2), 272–284.

: Golub, G. H., & Van Load, C. F. (1989). Matrix computation (2nd ed.). Baltimore: The John Hopkins University Press.

I should also mention that solutions to symmetric matrices, such as what you must be solving given your use of scipy.linalg.eigh, are somewhat cheap. If you are only interested in a few eigenvalues, you might find speed improvements in your method as well. The Arnoldi method is often used in such situations.

• thanks for the pointer, the QR algorithm seems like a good starting point Apr 26 '17 at 20:23
• I don't think that the magnitude of the perturbation in the eigenvalues is related to the condition number. This because $A$ has the same eigenvalues as $A+\lambda I$, but a different condition number. Apr 29 '17 at 11:01
• ps: linalg.eigh on a 4k-by-4k matrix is taking about 20 seconds (it only uses single core for some reason). I need about 0.25 second per update May 3 '17 at 16:47

There exist special techniques for updating the eigen-decomposition of time-dependent covariance matrices. Given a "prior" eigenvalue decomposition (say at some initial time $t^0$), these recursive algorithms lower the complexity of the spectrum update from $\mathcal{O}(N^3)$ (essentially the cost of a new eigendecomposition) to $\mathcal{O}(k N^2)$ where $N$ is the size of your matrix and $k$ is the rank of your update.

Here's a couple of relevant references:

Adaptive Eigendecomposition of Data Covariance Matrices Based on First-Order Perturbations (Champagne, IEEE TSP 42(10) 1994)

Recursive updating the eigenvalue decomposition of a covariance matrix (Yu, IEEE TSP, 39(5) 1991)

A Stable And Fast Algorithm For Updating The Singular Value Decomposition (Gu and Eisenstadt, 1994)

• unfortunately I don't have small rank updates, I have small norm updates of full rank Apr 24 '17 at 16:06
• @YaroslavBulatov I'm not aware of an efficient algorithm that can handle small-norm full rank updates - best I could find was this reference, but it does not look very promising. There is of course a large body of literature on eigenvalue perturbation that you may want to look at (see the other answer). Apr 24 '17 at 19:12