Let $M$ be an \(n \times n\) matrix. A scalar \(\lambda \) is an eigenvalue of $M$ if there is a non-zero vector $x$ (called eigenvector) s.t.:
$$M x = \lambda x$$
Eigenvalues and eigenvectors are crucial in many fields of science. For example, consider a discrete-time and discrete states Markov chain, whose transition matrix $M$ is defined as follows:
Let the initial state vector \(x_1\) be:
We know that from $M$ and \(x\_1\) we could compute all the successive states and it's true that:
$$x\_2 = M x \_1$$
$$x\_3 = M x\_2$$
and in general
$$x\k = M x\{k-1}$$
We may want to find a vector $x$ s.t.
$$Mx = x$$
Vectors with this property as known as steady-state vectors. It can be demonstrated that finding steady-state vectors equals finding any eigenvector $x$ with eigenvalue 1.
For example, the steady-state vector for the matrix $M$ is:
and one can easily show that
Finding eigenvalues and eigenvectors is not always easy to do by hand, and there are some algorithms to compute them. Unfortunately, this calculation may be expensive, especially with large matrices, and the result may be inaccurate due to approximations.
However, some algorithms perform better than others, and I want to discuss some of them in this article.
Solving characteristic equation
We can rewrite \(M x = \lambda x\) as
$$M x-\lambda x = 0$$
$$( M-\lambda I)x = 0$$
This system has a non-trivial solution (i.e. \(x \neq 0\)) only if \(det(M-\lambda I) =0\). \(det(M-\lambda I) =0\) is known as the characteristic equation.
Expanding \(det(M-\lambda I) =0\) we obtain a polynomial of degree $n$, whose roots are the eigenvalues of $M$. Computing eigenvectors from eigenvalues is trivial: for each eigenvalue \(\lambda\), we just need to find the null space of the matrix \(M-\lambda I\).
This is how we compute eigenvalues and eigenvectors by hand, but following this approach on a computer leads to some problems:
it depends on the computation of the determinant, which is a time-consuming process (due to the symbolic nature of the computation);
there is no formula for solving polynomial equations of degree higher than 4. Even though some techniques exist, like Newton's method, it's tough to find all the roots.
Therefore we need a different approach.
Iterative methods
Unfortunately, there is no simple algorithm to directly compute eigenvalues and eigenvectors for general matrices (there are special cases of matrices where it's possible, but I won't cover them in this article).
However, there are iterative algorithms that produce sequences that converge to eigenvectors or eigenvalues. There are several variations of these methods, I'll just cover two of them: the power method and the QR algorithm.
The power method
This method applies to matrices that have a dominant eigenvalue \(\lambda\_d\) (i.e. an eigenvalue that is larger in absolute value than the other eigenvalues).
Let $M$ be an \(n \times n\) matrix, the power method approximates a dominant eigenvector in the following steps:
$$x\_1 = Mx\_0$$
$$x\_2 = Mx\_1$$
$$x\k = Mx\{k-1}$$
And the more steps we take (i.e. the bigger $k$ is) the more accurate will be our approximation. This is expressed in the following formula
Once we have an approximation of the dominant eigenvector \(x\_d\) we find the corresponding dominant eigenvalue \(\lambda\_d\) with the Rayleigh quotient
$$\frac{(Mx)x}{xx} = \frac{(\lambda\_d x)x}{xx} = \frac{\lambda\_d (xx)}{xx}\lambda\_d$$
Once we have \(\lambda\_d\), we use the observation that if \(\lambda\) is an eigenvalue of $M$, \(\lambda - \beta\) is an eigenvalue of \(M-\beta I\) for any scalar \(\beta\). We can then apply the power method to compute a second eigenvalue. Repeating this process will allow us to compute all of the eigenvalues.
In Python this is:
import numpy as np
def power_method(M, n_iter = 100):
n = M.shape[0]
x_d = np.repeat(.5, n)
lambda_d = n
for i in range(n_iter):
x_0 = x_d
x_d = np.matmul(M, x_0)
lambda_d = np.matmul(np.matmul(M, x_d), x_d) / np.matmul(x_d, x_d)
h = np.zeros((n, n), int)
np.fill_diagonal(h, lambda_d)
N = M - h
x_1 = np.array([1, 0])
lambda_1 = n
for j in range(n_iter):
x_0 = x_1
x_1 = np.matmul(N, x_0)
lambda_1 = np.matmul(np.matmul(M, x_1), x_1) / np.matmul(x_1, x_1)
return [[x_d, lambda_d], [x_1, lambda_1]]
The function above works only for \(2 \times2\) matrices, but can easily be modified to \(n\times n\) matrices. We now test the function:
Matr = np.array([[1, 3], [2, 1]])
power_method(Matr)
We can even prove that those values represent a good approximation by checking the equation
$$Mx=\lambda x$$
Since this is an approximation, the == operator is not suited, we define instead the is_close function.
def is_close(x, y):
if all(abs(x-y) < 1e-5):
return True
else:
return False
Matr = np.array([[.7, .2], [.3, .8]])
sol = power_method(Matr)
lambda_a = sol[0][1]
lambda_b = sol[1][1]
x_a = sol[0][0]
x_b = sol[1][0]
print(is_close(np.matmul(Matr, x_a), lambda_a * x_a))
print(is_close(np.matmul(Matr, x_b), lambda_b * x_b))
Above we defined the algorithm as follows
$$x\k = Mx\{k-1}$$
We can notice that if
$$x\{k-1} = Mx\{k-2}$$
then we can substitute
$$x\k = MMx\{k-2}$$
By induction, we can prove that
$$x\k = M^kx{0}$$
We now use this formula to update the Python function above. The new function is the following:
def power_method_2(M, n_iter = 100):
n = M.shape[0]
x_d = np.array([1, 0])
M_k = np.linalg.matrix_power(M, n_iter)
M_k = M_k / np.max(M_k)
x_d = np.matmul(M_k, x_d)
x_d = x_d / np.max(x_d)
lambda_d = np.matmul(np.matmul(M, x_d), x_d) / np.matmul(x_d, x_d)
D = np.zeros((n, n), float)
np.fill_diagonal(D, lambda_d)
N = M - D
x_nd = np.array([1,0])
N_k = np.linalg.matrix_power(N, n_iter)
N_k= N_k / np.max(N_k)
x_nd = np.matmul(N_k, x_nd)
x_nd = x_nd/np.max(x_nd)
lambda_nd = np.matmul(np.matmul(N, x_nd), x_nd) / np.matmul(x_nd, x_nd)
lambda_nd = lambda_nd + lambda_d
return [[x_d, lambda_d], [x_nd, lambda_nd]]
Again we test the function:
Matr = np.array([[.7, .2], [.3, .8]])
sol_2 = power_method(Matr)
lambda_a = sol_2[0][1]
lambda_b = sol_2[1][1]
x_a = sol_2[0][0]
x_b = sol_2[1][0]
print(is_close(np.matmul(Matr, x_a), lambda_a * x_a))
print(is_close(np.matmul(Matr, x_b), lambda_b * x_b))
Once we are sure both the functions work correctly, we can now test which has a better performance.
import timeit
%timeit power_method(Matr)
%timeit power_method_2(Matr)
And we have a winner: the second function is 3.875 times faster than the first one.
The QR algorithm
One of the best methods for approximating the eigenvalues and the eigenvectors of a matrix applies the QR factorization and for this reason is known as the QR algorithm.
Let $M$ be an \(n\times n\) matrix, first of all, we need to factor it as
$$M = Q\_0R\_0$$
then we set
$$M\_1 = R\_0Q\_0$$
We then factor \(M\_1 = Q\_1R\_1\) and define \(M\_2 = R\_1Q\_1\) and so on.
It can be proven that $M$ is similar to \(M\_1,M\_1, \dots, M\_k\), which means $M $ and \(M\_1,M\_1, \dots, M\_k\) have the same eigenvalues.
It can also be shown that the matrices \(M\_k\) converge to a triangular matrix $T$ and that the elements on the diagonal are the eigenvalues of \(M\_k\).
In Python this is:
import numpy as np
def QR_argo(M, n_iter = 100):
n = M.shape[1]
Q_k = np.linalg.qr(M)[0]
R_k = np.linalg.qr(M)[1]
e_values = []
for i in range(n_iter):
M_k = np.matmul(R_k, Q_k)
Q_k = np.linalg.qr(M_k)[0]
R_k = np.linalg.qr(M_k)[1]
for j in range(M_k.shape[1]):
e_values.append(M_k[j, j])
return e_values
We can now test the function and compare it to the power method.
def is_close(x, y):
if abs(x-y) < 1e-5:
return True
else:
return False
Matr = np.array([[1, 3], [2, 1]])
pow_lambda_a = power_method(Matr)[0][1]
pow_lambda_b = power_method(Matr)[1][1]
QR_lambda_a = QR_argo(Matr)[0]
QR_lambda_b = QR_argo(Matr)[1]
is_close(QR_lambda_a, pow_lambda_a)
is_close(QR_lambda_b, pow_lambda_b)
Once we have eigenvalues \(\lambda\_i\), computing eigenvectors is easy: they are the non-trivial solution of
$$(M-\lambda\_i I) x=0$$
And that's it for this article.
Thanks for reading.
For any questions or suggestions related to what I covered in this article, please add them as a comment. In case of more specific inquiries, you can contact me here.
