Algorithm for finding eigenvectors from Hubbard & Hubbard's Linear Algebra book

eigenvectors.py

eps = 1e-6
max_newton_iters = 1000

def vec(arr):
    return [[a] for a in arr]

def matmul(A, B):
    assert len(A[0]) == len(B)
    r = [[0 for _ in range(len(B[0]))] for _ in range(len(A))]

    for i in range(len(A)):
        for j in range(len(B[0])):
            for k in range(len(A[0])):
                r[i][j] += A[i][k] * B[k][j]
    return r

def gaussian_elimination(matrix):
    n = len(matrix[0])
    m = len(matrix)
    for c in range(n):
        if c >= m: continue
        pivot = matrix[c][c]
        if pivot == 0: continue
        for c2 in range(n):
            matrix[c][c2] /= pivot
        assert(matrix[c][c] == 1)
        for row in range(c+1, m):
            fact = matrix[row][c]
            for col in range(n):
                matrix[row][col] -= fact * matrix[c][col]
    for c in range(n):
        if c >= m: continue
        for row in range(0, c):
            fact = matrix[row][c]
            for col in range(n):
                matrix[row][col] -= fact * matrix[c][col]
    return matrix

def eval_poly(coeffs, t):
    x = 1
    r = 0
    for c in coeffs:
        r += c*x
        x = t*x
    return r

def eval_diff(coeffs, t):
    x = 1
    r = 0
    for n, c in enumerate(coeffs[1:]):
        r += c*(n+1)*x
        x = t*x
    return r

def newtons_method(func, diff):
    x = 0
    for _ in range(max_newton_iters):
        yn = func(x)
        if abs(yn - 0) < eps:
            return x
        m = diff(x)

        if m == 0: # avoid divide by zero
            x = random.randrange(-10, 10)
            continue

        x_new = x - yn/m
        if abs(x - x_new) < eps:
            return x_new
        x = x_new
    raise Exception("[ERR] max_newton_iters exceeded.")

# divides coeff / (x - l)
def poly_divide(coeffs, l):
    r = [0 for _ in coeffs]
    sub = 0
    for i, c in reversed(list(enumerate(coeffs))):
        q = c - sub
        # (c x^i)/(x - l) => c x^(i - 1)
        r[i - 1] = q
        sub = -q * l
    if abs(r[-1]) > eps:
        print("[WARN] Polynomial division error greater than eps: ", r[-1])
    return r[:-1]

def get_all_roots(coeffs):
    roots = []
    for _ in coeffs[:-1]:
        poly = lambda t: eval_poly(coeffs, t)
        diff = lambda t: eval_diff(coeffs, t)
        r0 = newtons_method(poly, diff)
        coeffs = poly_divide(coeffs, r0)
        roots.append(r0)
    return roots

def calc_eigenvectors(eigenvalues, coeffs, powers_mat):
    r = []
    for l in eigenvalues:
        q = poly_divide(coeffs, l)
        q.append(0)
        eigenvector = matmul(powers_mat, vec(q))
        r.append(eigenvector)
    return r

def eigenanalysis(A):
    w = vec([int(i == 0) for i in range(len(A[0]))])

    powers = [w]
    for i in range(1, len(A) + 1):
        powers.append(matmul(A, powers[i - 1]))
    powers_mat = [[powers[j][i][0] for j in range(len(A) + 1)] for i in range(len(A))]
    reduced = gaussian_elimination(powers_mat)
    print("Row reduction finished")

    coeffs = [-reduced[i][-1] for i in range(len(powers_mat))]
    coeffs.append(1)

    eigenvalues = get_all_roots(coeffs)
    powers_mat = [[powers[j][i][0] for j in range(len(A) + 1)] for i in range(len(A))]
    eigenvectors = calc_eigenvectors(eigenvalues, coeffs, powers_mat)
    return eigenvalues, eigenvectors

def check_eigevectors(A, eigenvalues, eigenvectors):
    no_warn = True
    for l, v in zip(eigenvalues, eigenvectors):
        v1 = matmul(A, v)
        v2 = [l * vi[0] for vi in v]
        err = max(a[0] - b for a, b in zip(v1, v2))
        if err > eps:
            print("[WARN] check_eigevectors, err is ", err)
            no_warn = False
    if no_warn: print("check_eigevectors passed.")

def random_symetric_matrix(n, min_range=-1, max_range=1):
    A = [[random.uniform(min_range, max_range) for _ in range(n)] for _ in range(n)]
    for i in range(len(A)):
        for j in range(len(A[0])):
            if i < j:
                A[i][j] = A[j][i]
    return A


import random
# A = [[1, -1, 0], [-1, 2, -1], [0, -1, 1]]
A = random_symetric_matrix(5)
print(A)

eigenvalues, eigenvectors = eigenanalysis(A)
print("Eigenvalues: ", ', '.join(map(lambda l: "{:.2f}".format(l), eigenvalues)))
print("Eigenvectors:")
for v in eigenvectors:
    print('[', ', '.join(map(lambda vi: "{:.2f}".format(vi[0]), v)), ']')
check_eigevectors(A, eigenvalues, eigenvectors)