Eigenvalues and Eigenvectors Algorithms¶
We now discuss the numerical algorithms for computing eigenvalues and eigenvectors. We will mainly focus on the symmetric matrices with real eigenvalues and eigenvectors. For the general case, most of the tasks in data science need to conduct the singular value decomposition (SVD).
Recall the symmetric matrix
Sensitivity of Eigenvalues and Eigenvectors¶
There are two famous theorems about the sensitivity of eigenvalues and eigenvectors.
Theorem(Weyl's inequality): If
Theorem(Davis-Kahan Theorem): If
From the above two theorems, we can see that the eigenvalues are stable but the eigenvectors will be sensitive when the eigengap
Be aware that the above theorems are the simplified version. We can refer to the Matrix Perturbation Theory for more details.
Power Method¶
Here we consider the positive definite matrix
The power method applies the power of the matrix
Let us take a vector
where
Therefore, it follows that:
This indicates that, when
is a very good approximation of an eigenvector of
This motivates the power method to find the largest eigenvalue and its corresponding eigenvector.
Notice that we normalize the vector
import numpy as np
A = np.array([[4, 1],
[2, 3]])
tol = 1e-6
max_iter = 1000
# Initialize a random vector
n = A.shape[0]
x = np.random.rand(n)
for _ in range(max_iter):
# Multiply A with the vector
x = A @ x
x_norm = np.linalg.norm(x)
# Normalize the vector
x = x / x_norm
# Check for convergence
if np.linalg.norm(x_new - x) < tol:
break
Convergence and Sensitivity: We can see from the above analysis that the power method convergence speed is determined by
Other eigenvectors: The power method above only finds the largest eigenvalue and its corresponding eigenvector. As
We can see that the inverse power method needs to solve a linear system. You can improve the efficiency by apply the LU decomposition first. Similar to the above analysis, the convergence speed of the inverse power method is determined by
We can further apply the inverse power method to
QR Method¶
Though the power method is simple and easy to implement, it depends on the eigenvalue distribution. Besides, we can only find one eigenvalue and its corresponding eigenvector using the power method. We now discuss how to get the full eigen-decomposition using the QR method.
QR Decomposition: For any matrix
QR decomposition can be applied to any matrix
where
QR Iteration: We can use the QR iteration to find the eigenvalues and eigenvectors: Given
Theorem(Convergence of QR Iteration): If
In fact, we can see from the QR iteration that
Therefore, the QR algorithm can be seen as a more sophisticated variation of the basic "power" eigenvalue algorithm.
Here's a simple implementation of the QR iteration method:
import numpy as np
# Create a symmetric test matrix
n = 4
A = np.random.rand(n, n)
A = (A + A.T) / 2
# Initialize
A_k = A.copy()
V = np.eye(n) # To accumulate eigenvectors
tol = 1e-8
max_iter = 1000
# QR iteration
for _ in range(max_iter):
# QR decomposition
Q, R = np.linalg.qr(A_k)
# Update A_k
A_k = R @ Q
# Accumulate eigenvectors
V = V @ Q
# Check convergence of off-diagonal elements
if np.sum(np.abs(A_k - np.diag(np.diag(A_k)))) < tol:
break
# Eigenvalues are on the diagonal of final A_k
eigenvals = np.diag(A_k)
Practical Algorithms¶
We will briefly discuss the practical algorithms for computing eigenvalues and eigenvectors below. Please refer to the specialized textbooks, e.g., Matrix Computations, for more details.
The QR iteration method above is just the basic version. It can be further accelerated by applying the shift technique like the power method. In fact, the state-of-the-art Francis algorithm applies the famous implicit double-shift with no QR decomposition being explicitly performed. In fact, instead of QR decomposition, the Francis algorithm decomposes matrix into
When
For SVD decomposition
Python Implementation¶
Like the linear equation, both numpy
and scipy
provide the linalg
module for the eigenvalue and eigenvector problems.
You may need to consider to use different functions based on:
- Whether
is symmetric - Whether you only need the eigen/singular values
- Whether you only need the top
or the complete decomposition
import numpy as np
A = np.array([[1, 2],
[2, 3]])
# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)
# If you only need eigenvalues
eigenvalues = np.linalg.eigvals(A)
# Eigendecomposition for symmetric matrices
eigenvalues, eigenvectors = np.linalg.eigh(A)
# If you only need eigenvalues
eigenvalues = np.linalg.eigvalsh(A)
# Perform SVD
U, S, VT = np.linalg.svd(A, full_matrices=False) #If full_matrices = True (default), U and VT are of shape (M, M), (N, N). If False, the shapes are (M, K) and (K, N), where K = min(M, N).
# If you only need singular values
singular_values = np.linalg.svdvals(A)
# Compute the largest k eigenvalues and eigenvectors, when k << n
from scipy.sparse.linalg import eigs, eigsh, svds
k = 1
eigenvalues, eigenvectors = eigs(A, k=k, which='LM') # 'LM' selects largest magnitude eigenvalues
# For symmetric matrices
eigenvalues, eigenvectors = eigsh(A, k=k, which='LA') # 'LA' selects largest algebraic eigenvalues
# Top k singular values and vectors
U, S, VT = svds(A, k=k)
The scipy.sparse.linalg
package provides more efficient implementations for large sparse matrices and it implements the specialized algorithms if you only need top scipy.sparse.linalg
for the top