Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
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?
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 [72]. 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 [14]. Sarkar et al. [73] 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 [6] 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$ [74].
[6]: 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).
[14]: Thompson, P. A. (1980). An adaptive spectral analysis technique for unbiased frequency
[72]: Bradbury, W. W., & Fletcher, R. (1966). New iterative methods for solutions of the eigenproblem. Numerical Mathematics, 9(9), 259–266.
[73]: 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.
[74]: 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.
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)
Online Principal Component Analysis in High Dimension: Which Algorithm to Choose? (Cardot and Degras)
A Stable And Fast Algorithm For Updating The Singular Value Decomposition (Gu and Eisenstadt, 1994)
