cds-numerical-methods/Week 4/8 Eigenvalues and Eigenvectors.ipynb

CDS: Numerical Methods Assignments¶

• See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.

• Solutions must be submitted via the Jupyter Hub.

• Make sure you fill in any place that says YOUR CODE HERE or "YOUR ANSWER HERE".

Submission¶

1. Name all team members in the the cell below
2. make sure everything runs as expected
3. restart the kernel (in the menubar, select Kernel$\rightarrow$Restart)
4. run all cells (in the menubar, select Cell$\rightarrow$Run All)
5. Check all outputs (Out[*]) for errors and resolve them if necessary
6. submit your solutions in time (before the deadline)

team_members = "Kees van Kempen, Koen Vendrig"

Eigenvalues and Eigenvectors¶

In the following you will implement your own eigenvalue / eigenvector calculation routines based on the inverse power method and the iterated QR decomposition.

In [ ]:

import numpy as np
import math as m

Task 1¶

We start by implementing the inverse power method, which calculates the eigenvector corresponding to an eigenvalue which is closest to a given parameter $\sigma$. In detail, you should implement a Python function $\text{inversePower(A, sigma, eps)}$, which takes as input the $n \times n$ square matrix $A$, the parameter $\sigma$, as well as some accuracy $\varepsilon$. It should return the eigenvector $\mathbf{v}$ (corresponding to the eigenvalue which is closest to $\sigma$) and the number of needed iteration steps.

To do so, implement the following algorithm. Start by setting up the needed input:

\begin{align} B &= \left( A - \sigma \mathbf{1} \right)^{-1} \\ \mathbf{b}^{(0)} &= (1,1,1,...) \end{align}

where $\mathbf{b}^{(0)}$ is a vector with $n$ entries. Then repeat the following and increase $k$ each iteration until the error $e$ is smaller than $\varepsilon$:

\begin{align} \mathbf{b}^{(k)} &= B \cdot \mathbf{b}^{(k-1)} \\ \mathbf{b}^{(k)} &= \frac{\mathbf{b}^{(k)}}{|\mathbf{b}^{(k)}|} \\ e &= \sqrt{ \sum_{i=0}^n \left(|b_i^{(k-1)}| - |b_i^{(k)}|\right)^2 } \end{align}

Return the last vector $\mathbf{b}^{(k)}$ and the number of needed iterations $k$.

Verify your implementation by calculating all the eigenvectors of the matrix below and comparing them to the ones calculated by $\text{numpy.linalg.eig()}$. Then implement a unit test for your function.

\begin{align*} A = \begin{pmatrix} 3 & 2 & 1\\ 2 & 3 & 2\\ 1 & 2 & 3 \end{pmatrix}. \end{align*}

In [ ]:

def inversePower(A, sigma, eps):
    """
    Estimates the eigenvectors of matrix A by the inverse power method.

    Args:
        A: an n x n matrix
        sigma: an initial guess for an eigenvector
        eps: the desired accuracy

    Returns:
        A tuple (b, k) is returned, containing:
            b: the eigenvector b corresponding to the eigenvalue
               closests to sigma after k iterations
            k: the number of iterations done
            
    See also:
        https://www.sciencedirect.com/topics/mathematics/inverse-power-method
    """
    
    # Does https://johnfoster.pge.utexas.edu/numerical-methods-book/LinearAlgebra_EigenProblem1.html help?
    
    # A should be n x n.
    n = len(A)
    assert len(A.shape) == 2 and A.shape[0] == A.shape[1]
    
    # Setup some initial values.
    #B = np.linalg.inv(A - sigma*np.ones(n))
    B = np.linalg.inv(A - sigma*np.eye(n))
    #B = 1/(A - sigma*np.ones(n))
    #B = 1/(A - sigma*np.eye(n))
    b = np.ones(n)
    k = 0
    e = 0
    
    while e > eps or k == 0:
        b_prev = b.copy()
        k += 1
        
        b = B @ b
        b /= np.sqrt(b @ b) # although b = np.linalg.norm(b) could be used
        e = np.sqrt(np.sum( (np.abs(b_prev) - np.abs(b))**2) )
    
    return b, k

In [ ]:

# Use this cell for your own testing ...

A = np.array([[3, 2, 1], [2, 3, 2], [1, 2, 3]])

sigma_list = [.03, 6., 1.99999989999999]
#sigma = 0.03

for sigma in sigma_list:
    print("For sigma =", sigma)
    b, k = inversePower(A, sigma, 1e-6)
    print("b =", b)
    lam = np.dot(A, b) / b
    print("lam =", lam)
    print()
    print()

print(np.linalg.eig(A)[0])

In [ ]:

def test_inversePower():
    A = np.array([[3, 2, 1], [2, 3, 2], [1, 2, 3]])
        
test_inversePower()

Task 2¶

Next you will need to implement the tri-diagonalization scheme following Householder. To this end implement a Python function $\text{tridiagonalize(A)}$ which takes a symmetric matrix $A$ as input and returns a tridiagonal matrix $T$ of the same dimension. Therefore, your algorithm should execute the following steps:

First use an assertion to make sure $A$ is symmetric. Then let $k$ run from $0$ to $n-1$ and repeat the following:

\begin{align} q &= \sqrt{ \sum_{j=k+1}^n \left(A_{j,k}\right)^2 } \\ \alpha &= -\operatorname{sgn}(A_{k+1,k}) \cdot q \\ r &= \sqrt{ \frac{ \alpha^2 - A_{k+1,k} \cdot \alpha }{2} } \\ \mathbf{v} &= \mathbf{0} \qquad \text{... vector of dimension } n \\ v_{k+1} &= \frac{A_{k+1,k} - \alpha}{2r} \\ v_{k+j} &= \frac{A_{k+j,k}}{2r} \quad \text{for } j=2,3,\dots,n \\ P &= \mathbf{1} - 2 \mathbf{v}\mathbf{v}^T \\ A &= P \cdot A \cdot P \end{align}

At the end return $A$ as $T$.

Apply your routine to the matrix $A$ defined in task 1 as well as to a few random (but symmetric) matrices of different dimension $n$.

Hint: Use $\text{np.outer()}$ to calculate the matrix $\mathbf{v}\mathbf{v}^T$ needed in the definition of the Housholder transformation matrix $P$.

In [ ]:

def tridiagonalize(A):
    """
    Tridiagonalizes a given matrix following the Householder method.
    Args:
        A: NxN matrix
    Returns:
        The matrix A, but now tridiagonalized
    """
    assert A.shape[0] == A.shape[1] and len(A.shape) == 2
    n = len(A)
    for k in range(n-1):
        q = m.sqrt(np.sum(A[k+1:n,k]**2))
        alpha = -1*np.sign(A[k+1,k])*q
        r = m.sqrt((alpha**2 - A[k+1,k]*alpha)/2)
        v = np.zeros(n)
        v[k+1] = (A[k+1,k]-alpha)/(2*r)
        print(np.shape(v),np.shape(A),np.shape(r),np.shape(A[k+2:n]/(2*r)))
        v[k+2:n] = A[k+2:n,k]/(2*r)
        P = np.identity(n)-(2*np.outer(v,v))
        A = np.dot(np.dot(P,A),P)
    return A

In [ ]:

B = np.array([[3,2,1],[2,3,2],[1,2,3]])

tridiagonalize(B)

Task 3¶

Implement the $QR$ decomposition based diagonalization routine for tri-diagonal matrices $T$ in Python as a function $\text{QREig(T, eps)}$. It should take a tri-diagonal matrix $T$ and some accuracy $\varepsilon$ as inputs and should return all eigenvalues in the form of a vector $\mathbf{d}$. By making use of the $QR$ decomposition as implemented in Numpy's $\text{numpy.linalg.qr()}$ the algorithm is very simple and reads:

Repeat the following until the error $e$ is smaller than $\varepsilon$:

\begin{align} Q \cdot R &= T^{(k)} \qquad \text{... do this decomposition with the help of Numpy!} \\ T^{(k+1)} &= R \cdot Q \\ e &= |\mathbf{d_1}| \end{align}

where $T^{(0)} = T$ and $\mathbf{d_1}$ is the first sub-diagonal of $T^{(k+1)}$ at each iteration step $k$. At the end return the main-diagonal of $T^{(k+1)}$ as $\mathbf{d}$.

Implement a unit test for your function based on the matrix $A$ defined in task 1. You will need to tri-diagonalize it first.

In [ ]:

def QREig(T, eps):
    """
    Follows the method of the QR decomposition based diagonalization routine for tridiagonalized matrices. The matrix T is
    diagonalized, resulting in all diagonal elements being an eigenvalue.
    Args:
        T: a already tridiagonalized matrix
        eps: the desired accuracy
    Returns:
        A one dimensional array with the eigenvalues of the matrix T
    """
    e = eps + 1
    while e > eps:
        Q,R = np.linalg.qr(T)
        T = np.matmul(R,Q)
        e = np.sum(np.absolute(np.diag(T,k=1)))
    #print(np.diag(T))
    return np.diag(T)

In [ ]:

A_tridiag = tridiagonalize(A)
print(A_tridiag)

In [ ]:

def test_QREig():
    """
    Tests the QREig function for the matrix A_tridiag (defined in task 2) and eps = 0.001
    """
    eps = 0.001
    x = QREig(A_tridiag,eps)
    print(x)

test_QREig()

Task 4¶

With the help of $\text{QREig(T, eps)}$ you can now calculate all eigenvalues. Then you can calculate all corresponding eigenvectors with the help of $\text{inversePower(A, sigma, eps)}$, by setting $\sigma$ to approximately the eigenvalues you found (you should add some small random noise to $\sigma$ in order to avoid singularity issues in the inversion needed for the inverse power method).

Apply this combination of functions to calculate all eigenvalues and eigenvectors of the matrix $A$ defined in task 1.

In [ ]:

"""
Eigenvalues and -vectors calculated numerically for matrix A using functions from previous tasks.
"""
eps = 0.001
T = QREig(A_tridiag,eps)
for i in range(len(T)):
    sigma = T[i] + np.random.normal(0,0.01,1)
    x = inversePower(A_tridiag,sigma,eps)
    print(x)

Task 5¶

Test your eigenvalue / eigenvector algorithm for larger random (but symmetric) matrices.

In [ ]:

"""
Completely 'random' and totally not self-made matrices C and D are being tridiagonalized and its eigenvalues and -vectors are
calculated.
"""
C = np.array([[5,4,3,2,1]
            ,[3,4,5,1,2]
            ,[5,1,4,2,3]
            ,[3,2,1,5,1]
            ,[5,4,3,2,1]])
#print(C)
C_tridiag = tridiagonalize(C)
T = QREig(C_tridiag,eps)
for i in range(len(T)):
    sigma = T[i] + np.random.normal(0,0.01,1)
    x = inversePower(C_tridiag,sigma,eps)
    print(x)

D = np.array([[6,4,3,2,1]
            ,[3,3,5,1,2]
            ,[5,1,4,2,2]
            ,[3,1,1,5,6]
            ,[1,4,7,2,3]])
#print(C)
C_tridiag = tridiagonalize(C)
T = QREig(C_tridiag,eps)
for i in range(len(T)):
    sigma = T[i] + np.random.normal(0,0.01,1)
    x = inversePower(C_tridiag,sigma,eps)
    print(x)