# 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.