4. Finding Eigenvalues and Eigenvectors#
In this chapter we implement basic methods for findinging eigenvalues and eigenvectors. The corresponding theory is covered in Chapter 4 of the lecture notes.
4.1. Norms#
We begin by implementing a few norms which will be useful in the methods we develop below.
import numpy as np
def norm_loo(x:np.array):
"""
l_\infty norm
"""
return np.max(np.abs(x))
def norm_l1(x:np.array):
"""
l_1 norm
"""
return np.sum(np.abs(x))
def norm_l2(x:np.array):
"""
l_2 norm
"""
return np.sqrt(sum(x**2))
4.2. Power method#
This is a basic implementation of the power method. In addition to the matrix and the starting vector, we pass the function \(\phi:\mathbb{C}^n\to \R\) we on the approximations of the eigenvectos, a norm, the maximum number of iterations, and the desider error bound.
def power_method(A:np.array,x0:np.array,phi,norm,max_iterations:int, epsilon:float):
"""
Implementation of the Power method.
Arguments
- A : n x n matrix (numpy array)
- x0 : the initial approximation fo the eigenvector
- phi: A function from C^n to R
- norm: A function from C^n to R
- max_iterations: the maximum number of iterations
- epsilon: the error tolerance
"""
# Inizializes the current vector (normalized)
x = x0/norm(x0)
# Initializes the error
err = float('inf')
# Initializes the number of iterations
k = 1
while err > epsilon and k <= max_iterations:
# Finds the new vector
x1 = np.matmul(A,x)
# Finds the current estimate of the eigenvalue
r = phi(x1)/phi(x)
x1 = x1/norm(x1)
# Updates the error
err = norm(x1-x)
# Prints a log
print(f"{k:<10} {str(x):<25} {str(x1):<25} {r:<10.6f}")
# Updates the current vector
x = x1
# Updates the iteration counter
k = k + 1
return r
4.2.1. Example 1#
A = np.array([[3,1],[1,3]])
x0 = np.array([0.5,1])
def phi(x:np.array):
"""
This function returns the sum of the component of the input vector.
"""
return np.sum(x)
e = power_method(A,x0,phi,norm_l1,50,10**-5)
print(f"The dominant eigenvalue is {e}")
1 [0.33333333 0.66666667] [0.41666667 0.58333333] 4.000000
2 [0.41666667 0.58333333] [0.45833333 0.54166667] 4.000000
3 [0.45833333 0.54166667] [0.47916667 0.52083333] 4.000000
4 [0.47916667 0.52083333] [0.48958333 0.51041667] 4.000000
5 [0.48958333 0.51041667] [0.49479167 0.50520833] 4.000000
6 [0.49479167 0.50520833] [0.49739583 0.50260417] 4.000000
7 [0.49739583 0.50260417] [0.49869792 0.50130208] 4.000000
8 [0.49869792 0.50130208] [0.49934896 0.50065104] 4.000000
9 [0.49934896 0.50065104] [0.49967448 0.50032552] 4.000000
10 [0.49967448 0.50032552] [0.49983724 0.50016276] 4.000000
11 [0.49983724 0.50016276] [0.49991862 0.50008138] 4.000000
12 [0.49991862 0.50008138] [0.49995931 0.50004069] 4.000000
13 [0.49995931 0.50004069] [0.49997965 0.50002035] 4.000000
14 [0.49997965 0.50002035] [0.49998983 0.50001017] 4.000000
15 [0.49998983 0.50001017] [0.49999491 0.50000509] 4.000000
16 [0.49999491 0.50000509] [0.49999746 0.50000254] 4.000000
The dominant eigenvalue is 4.0
4.3. Inverse power method#
To implement the inverse power method we need implementations of the forward and backward substitution methods, as well as of an LU factorization method. We provided implementations in Chapter 3. We need, however, to adjust them slightly, in order to be able to work with complex numbers. The adjusted implementations are provided below.
4.3.1. Adjusted code for linear systems#
import numpy as np
def forward_substitution(L,b):
'''
Implements the forward substitution algorithm to solve a lower triangular system.
- L: an n x n real lower triangular matrix (numpy array)
- b: an n dimensional real vector (numpy array)
'''
# Identifies the dimension of the matrix
n = L.shape[0]
# Initializes the solution as a vector of zeros
x = np.zeros(n,dtype="complex") # We ensure that we perform complex arithmetic operations.
for i in range(n):
# Initializes the numerator of the expression of x_i
numerator = b[i]
for j in range(i):
# Updates the numerator
numerator = numerator - L[i,j] * x[j]
# Divides by \ell_ii
x[i] = numerator / L[i,i]
# Returns the solution
return x
import numpy as np
def backward_substitution(U,b):
'''
Implements the backward substitution algorithm to solve an upper triangular system.
- U: an n x n real upper triangular matrix (numpy array)
- b: an n dimensional real vector (numpy array)
'''
# Identifies the dimension of the matrix U
n = U.shape[0]
# Initializes x with all zeros.
x = np.zeros(n,dtype="complex") # We ensure that we perform complex arithmetic operations.
# For i = n-1,...,0 with a step of -1
for i in range(n-1,-1,-1):
# Initializes the fraction that gives x_i
numerator = b[i]
# For j = i+1,...,n-1
for j in range(i+1,n):
numerator = numerator - U[i,j] * x[j]
x[i] = numerator / U[i,i]
# Returns x
return x
import numpy as np
def LUfactorization(A):
'''
Implements the LU factorization of a matrix A.
It uses a recursive algorithm.
- A: n x n numpy matrix
'''
# Get the length of the matrix
n=len(A)
# Initialize L and U with zeros. We ensure we perform complex arithmetic operations
L=np.zeros(shape=(n,n),dtype="complex")
U=np.zeros(shape=(n,n),dtype="complex")
# The operations on L_11 and U_11
# are always performed, independently of the size of the matrix.
L[0,0] = 1.0
U[0,0] = A[0,0]
# Base case
if n==1:
# If n = 1 we are done
return L,U
else:
L[0,0] = 1.0
U[0,0] = A[0,0]
# Populates row 1 (index 0) of U, from column 2 (index 1) to n (index n-1)
U[[0],1:n] = A[[0],1:n]/L[0,0]
# Populates column 1 (index 0) of L, from row 2 (index 1) to n (index n-1)
L[1:n,[0]] = A[1:n,[0]]/U[0,0]
# Recursion step on the matrix M = A22-\ell_1u_1
M = A[1:n,1:n] - np.matmul(L[1:n,[0]],U[[0],1:n])
(L22,U22) = LUfactorization(M)
# Puts the decomposed submatrices in their respective positions
L[1:n,1:n] = L22
U[1:n,1:n] = U22
return L,U
4.3.2. Inverse power method code#
def inverse_power_method(A:np.array,x0:np.array,phi,norm,max_iterations:int, epsilon:float):
"""
Implementation of the inverse power method.
The termination criteria is based on convergence of the eigenvalue.
- A: n x n matrix (np.array)
- x0: starting vector (np.array)
- phi: A linear function from C^n to R
- norm: A function from C^n to R
- max_iterations: the maximum number of iterations
- epsilon: the tolerance
"""
# Sets the initial approximation of the eigenvector
x = x0/norm(x0)
# Factorizes A
L,U = LUfactorization(A)
# Initializes the error
err = float('inf')
# Initializes r
r0 = float('inf')
# Initializes the number of iterations
k = 1
while err > epsilon and k <= max_iterations:
# Finds x_k+1 by solving Ax_k+1=x_k
y = forward_substitution(L,x)
x1 = backward_substitution(U,y)
# Finds the current estimate of the eigenvalue (should be the reciprocal for inverse power method)
r = phi(x1)/phi(x)
# Normalizes the vector using norm(x1)
x1 = x1/norm(x1)
# Updates the error
err = norm(r-r0)
# Prints a log
print(f"{k:<10} {r:<10.6f} {err:<10.6f}")
# Updates the current vector
x = x1
# Updates r
r0 = r
# Updates the iteration counter
k = k + 1
return r
4.3.3. Example#
A = np.array([[-1+1j,0,1/4],[1/4,1,1/4],[1,1,3]],dtype = 'complex')
x0 = np.array([1+1j,1+1j,1+1j])
def phi(x:np.array):
"""
e^\top x
"""
return np.sum(x)
e1 = inverse_power_method(A,x0,phi,norm_loo,50,10**-5)
print(1/e1)
1 0.242622-0.084922j inf
2 0.865504+0.755846j 1.046361
3 0.808433-0.537594j 1.294699
4 0.697444+0.205557j 0.751394
5 1.539751-0.010838j 0.869660
6 0.993897-0.121592j 0.556978
7 1.136335+0.135067j 0.293534
8 1.189864-0.058735j 0.201058
9 1.095477+0.008738j 0.116024
10 1.164797+0.021265j 0.070443
11 1.138086-0.013102j 0.043526
12 1.135787+0.012730j 0.025933
13 1.147652+0.002228j 0.015845
14 1.138055+0.002044j 0.009599
15 1.142213+0.006094j 0.005805
16 1.142041+0.002572j 0.003527
17 1.140665+0.004204j 0.002135
18 1.141950+0.004050j 0.001295
19 1.141316+0.003588j 0.000785
20 1.141406+0.004055j 0.000476
21 1.141559+0.003810j 0.000288
22 1.141390+0.003855j 0.000175
23 1.141484+0.003905j 0.000106
24 1.141463+0.003844j 0.000064
25 1.141447+0.003880j 0.000039
26 1.141469+0.003871j 0.000024
27 1.141455+0.003866j 0.000014
28 1.141459+0.003874j 0.000009
(0.8760613574842928-0.002972919686608703j)
4.4. Exercises#
Exercise 4.1
Modify the implementation of the power method in order to use a termination criteria based on the approximated eigenvalue. That is, terminate when \(|r_k-r_{k-1}|\) or \(|r_k-r_{k-1}|/|r_k|\) is within the desired tolerance.
Exercise 4.2
Modify the implementation of the inverse power method in order to use a termination criteria based on the approximated eigenvector. That is, terminate when the norm of \(x_k-x_{k-1}\) is within the desired tolerance.
Exercise 4.3
Modify the implementation of the power method in order to use a termination criteria based on the approximated eigenpair. That is, terminate when \(||Ax_k-\lambda_kx_k||\) is within the desired tolerance.