Demo entry 6345731

QR code

   

Submitted by QR_code_python on Feb 05, 2017 at 12:37
Language: Python. Code size: 3.7 kB.

"""
% Santiago Rendon-Valencia
% Cod. 22031680
% ANLA WS2016 - Assignment 7
% FAU Erlangen Nuremberg
% python 2.7
%
"""
import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt


def house(A):
    m, n = A.shape
    W = np.zeros((m, n))
    for k in range(n):
        x = A[k:m, k]
        e = np.zeros(len(x))
        e[0] = 1.0
        v = np.sign(x[0]) * np.linalg.norm(x) * e + x
        v /= np.linalg.norm(v)
        W[k:m, k] = v
        A[k:m, k:n] -= 2.0 * np.outer(v, np.dot(v, A[k:m, k:n]))
    return W, A


def formQ(W):
    m, n = W.shape
    Q = np.eye(m)
    for e in range(n):
        x = Q[:, e]
        for k in range(n - 1, -1, -1):
            v = W[:, k]
            x -= 2.0 * v * np.dot(v, x)
        Q[:, e] = x
    return Q


def tridiag(A):
    m, n = A.shape
    T = np.copy(A).astype(float)
    for k in range(n - 2):
        v = np.copy(T[k + 1:, k])
        v[0] += np.sign(v[0]) * np.linalg.norm(v)
        v /= np.linalg.norm(v)
        T[k + 1:, k:] -= 2.0 * np.outer(v, np.dot(v, T[k + 1:, k:]))
        T[:, k + 1:] -= 2.0 * np.outer(np.dot(T[:, k + 1:], v), v)
    T = (np.transpose(T) + T) / 2.0
    return T


def qralg(T):
    m, n = T.shape
    errors = []
    if m == 1:
        return T, errors  # Check if the matrix is 1x1, eigenvalue is just the number

    error = np.abs(T[m - 1, m - 2])
    errors.append(error)
    while error > 1.0e-12:
        W, R = house(T)
        Q = formQ(W)
        T = np.dot(R, Q)
        T = tridiag(T)
        error = np.abs(T[m - 1, m - 2])
        errors.append(error)
    return T, errors


def wilkinsonShift(B):
    m, n = B.shape
    delta = (B[-1, -1] - B[-2, -2]) / 2.0
    sign = np.sign(delta)
    if delta == 0:
        sign = 1
    B2 = B[-2, -1] * B[-2, -1]
    mu = B[-1, -1] - (sign * B2) / (np.abs(delta) + np.sqrt((delta * delta) + B2))
    return mu


def shiftedqralg(T):
    m, n = T.shape
    errors = []
    if m == 1:
        return T, errors  # Check if the matrix is 1x1, eigenvalue is just the number
    error = np.abs(T[m - 1, m - 2])
    errors.append(error)
    while error > 1.0e-12:
        mu = wilkinsonShift(T)
        shift = mu * np.eye(m)
        W, R = house(T - shift)
        Q = formQ(W)
        T = np.dot(R, Q) + shift
        T = tridiag(T)
        error = np.abs(T[m - 1, m - 2])
        errors.append(error)
    return T, errors


def qralg_driver(A, shift=0):
    eigenvals = []
    T = tridiag(A)
    m, n = T.shape
    globalErrors = []
    if shift:
        print "QR decomposition with Wilkinson shift:"
        name = "qr_with_shift"
    else:
        print "QR decomposition without shift:"
        name = "qr_without_shift"
    for k in range(m):
        if shift:
            T, error = shiftedqralg(T[:m - k, :m - k])
        else:
            T, error = qralg(T[:m - k, :m - k])
        eigenvalue = np.diagonal(T)[-1]
        eigenvals.append(eigenvalue)
        globalErrors.extend(error)
    print "Number of iterations: ", len(globalErrors)

    # 'sawtooth plot'
    plt.figure(1)
    t = np.arange(len(globalErrors))
    plt.semilogy(t, np.ones(len(t)) * 1.0e-12, 'r--', linewidth=0.5, label='criterion')
    plt.semilogy(globalErrors, 'mo')
    plt.semilogy(globalErrors, label='QR')

    plt.legend()
    plt.title('Number of iterations $Vs$ $abs(error)$')
    plt.xlabel('Number of iterations')
    plt.ylabel('$abs(error)$')
    # plt.show()
    plt.figure(1).savefig(name + '.pdf', format='pdf')
    return sorted(eigenvals, reverse=True)


if __name__ == "__main__":
    A = linalg.hilbert(4)
    # A = np.diag(np.arange(15, 0, -1)) + np.ones((15, 15))

    print qralg_driver(A)
    print qralg_driver(A, 1)

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).