EigenvaluesQRDecomp.py

#-------------------------------------------------------------------------------
# Name:        Numerical Project Assignment #2, Subject #3
# Purpose:     To determine the eigenvalues of a matrix using QR decomposition
#
# Author:      Addison Euhus
#
# Created:     04/30/2014
# Copyright:   (c) Addison Euhus 2014
#-------------------------------------------------------------------------------

from numpy import *
from numpy import random as rd
from pylab import *

#********************************************************

#Real Valued 4x4 Matrix
A = matrix('1 2 3 7; 4 5 6 7; 7 7 8 9; 1 7 -1 9')

#Complex Valued 4x4 Matrix
#A = matrix('-4 2 0 0; -2 -4 0 0; 0 0 2 3; 0 0 -3 2')

#Random Valued Matrix
#A = rd.random((100,100))

#Diagonal Matrix
#A = eye(7,7)

#*******************************************************

#Constant that is used to check for near-zero matrix values
#Smaller epsilon -> more accurate results for small-valued matrices
epsilon = 10**(-12)

#Iterations done before the program checks that the solution is
#a block matrix and thus calls a function to compute complex eigenvalues
iterate = 100

#Transforms parameter matrix A into matrix with Hessenberg form
def Hessenberg(A):
    #Checks that it is nxn, assigns the value to n
    if(A.shape[0] != A.shape[1]):
        return
    else:
        n = A.shape[0]

    #The following for loop iterates on the parameter matrix A
    #and constructs a matrix with Hessenberg form and equivalent
    #eigenvalues of A. Computes x1, u1, p1, P1, and eventually the new A
    for z in range(0,n-2):
        x1 = A[z+1:,z]
        vect = zeros(shape=(1,n-1-z))
        vect[0,0] = 1
        u1 = x1 - linalg.norm(x1)*vect.T
        if(linalg.norm(u1) == 0):
            continue
        p1 = eye(n-1-z) - 2*(u1*u1.T)/linalg.norm(u1)**2
        P1 = zeros(shape=(n,n))
        for x in range(0,z+1):
            P1[x,x] = 1
        for x in range(0,n-z-1):
            P1[x+z+1,z+1:] = p1[x]
        A = P1*A*P1
    return clean(A)


#Cleans A of approximately zero values using epsilon as the criteria
def clean(A):
    #Checks that it is nxn, assigns the value to n
    if(A.shape[0] != A.shape[1]):
        return
    else:
        n = A.shape[0]

    #Checks near-zero values in the matrix
    for x in range(0,n):
        for y in range(0,n):
            if(abs(A[x,y]) < epsilon):
                A[x,y] = 0
    return A

#Checks below the diagonal of A to determine if it is upper triangular
#Returns true if there are non-zero values below the diagonal
def scanLowerTriangle(A):
    #Checks that it is nxn, assigns the value to n
    if(A.shape[0] != A.shape[1]):
        return
    else:
        n = A.shape[0]

    #Checks below the diagonal for non-zero values
    for x in range(1,n):
        for y in range(0,x):
            if(abs(A[x,y]) > epsilon):
                return True
    return False

#Checks whether the parameter matrix A is a block matrix
#Returns true if it is a block matrix
def isBlockMatrix(A):
    #Checks that it is nxn and even, assigns the value to n
    if(A.shape[0] != A.shape[1] or A.shape[0] % 2 != 0):
        return
    else:
        n = A.shape[0]

    #Automatically returns true if 2x2
    if(n == 2):
        return True

    #Checks each block on the diagonal for non-zero values
    for x in range(2,n,2):
        if(abs(A[x,x-2]) > epsilon):
            return False
        if(abs(A[x,x-1]) > epsilon):
            return False
        if(abs(A[x+1,x-2]) > epsilon):
            return False
        if(abs(A[x+1,x-1]) > epsilon):
            return False
    return True

#Computes the eigenvalues of a block matrix , returns eigenvalues
def complexCompute(A):
    #Checks that it is nxn and even, assigns the value to n
    if(A.shape[0] != A.shape[1] or A.shape[0] % 2 != 0):
        return
    else:
        n = A.shape[0]

    #Evaluates the eigenvalues of each 2x2 block on the diagonal
    twoBytwos = []
    for x in range(0,n,2):
        ai = zeros(shape=(2,2))
        ai[0] = A[x,x:x+2]
        ai[1] = A[x+1,x:x+2]
        twoBytwos.append(ai)
    eigenvals = []
    for ai in twoBytwos:
        eigenvals.append(linalg.eig(ai)[0])
    return eigenvals

#Iterates on A using the QR-Algorithm in order to compute the eigenvalues
def eigenSearch(A):
    #Checks that it is nxn, assigns the value to n
    if(A.shape[0] != A.shape[1]):
        return
    else:
        n = A.shape[0]

    #Values used to check for complex eigenvalues/block matrices
    block = False
    q = 0

    #Iterative loop that using QR-Algorithm where A will converge to
    #either a diagonal or block matrix
    while(scanLowerTriangle(A) == True and block == False):
        Q = qr(A)[0]
        R = qr(A)[1]
        A = R*Q
        q = q+1
        if(q >= iterate):
            if(isBlockMatrix(A) == True and scanLowerTriangle(A) == True):
                block = True
            else:
                q = 0

    #Cleans A and then computes the eigenvalues of A by either checking
    #that it is a block matrix, and calling the complexCompute function,
    #or by pulling the eigenvalues directly off of the diagonal of the
    #triangular matrix
    A = clean(A)
    eigenvals = []
    if(block):
        eigenvals = complexCompute(A)
    else:
        for x in range(0,n):
            eigenvals.append(A[x,x])
    return eigenvals, A

#Runs the algorithm, prints the eigenvectors
print eigenSearch(Hessenberg(A))[0]

#Runs the algorithm, prints the decomposed matrix
#print eigenSearch(Hessenberg(A))[1]