__author__ = 'Lenovo'
import pprint
from pprint import pprint
from math import sqrt
import numpy as np
np.set_printoptions(suppress=True, linewidth=160, precision=10)
import scipy
from scipy import linalg
import pprint
def sign(x):
return cmp(x, 0)
def norm(x):
return sum(i ** 2 for i in x) ** 0.5
def house(x):
n = len(x)
mu = norm(x)
v = x
if mu != 0:
beta = x[0] + sign(x[0]) * mu
v = v / beta
v[0] = 1
return v
def rowHouse(a, v):
m = len(v)
beta = -2.0 / sum(i ** 2 for i in v)
vm = np.array([v]).transpose()
w = beta * np.dot(a.transpose(), vm)
p = np.dot(vm, w.transpose())
pa = a + p
return pa
def swap(a, r, k):
temp = a[:, k].copy()
a[:, k] = a[:, r]
a[:, r] = temp
def swapElems(c, r, k):
temp = c[k]
c[k] = c[r]
c[r] = temp
def qrColumnPivoting(permutation, a):
m = len(a)
n = len(a[0])
c = np.zeros(n)
for i in range(0, n):
c[i] = sum(j**2 for j in a[:,i])
r = 0
t = c.max()
k = c.argmax()
while t > 0:
permutation[r] = k
permutation[k] = r
swap(a, r, k)
swapElems(c, r, k)
v = house(a[r:m, r])
a[r:m, r:n] = rowHouse(a[r:m, r:n], v)
a[r+1:m, r] = v[1:m]
for i in range(r+1, n):
c[i] -= a[r, i]**2
if r < n-1:
t = c[r+1:].max()
k = c[r+1:].argmax()
else:
t = 0
r += 1
return permutation, a, r
def makeQ(a, r):
m = len(a)
q = np.identity(m)
for j in range(r-1, -1, -1):
v = np.hstack(([1], a[j+1:m, j]))
q[j:m, j:m] = rowHouse(q[j:m, j:m], v)
return q
def cleanA(a):
for i in range(0, len(a[0])):
a[i+1:len(a), i] = np.zeros(len(a) - (i+1))
return a
def fod(a, r):
rr = a[:r].transpose()
n = len(rr)
u = np.identity(n-r+1)
for i in range(0, r):
v = house(rr[r-1:n, i])
rr[r-1:n, i:r] = rowHouse(rr[r-1:n, i:r], v)
vt = np.array([v]).transpose()
vsq = np.dot(vt, np.array([v]))
new_u = np.identity(n-r+1) - 2/sum(j ** 2 for j in v)*vsq
u = np.dot(new_u, u)
#u = np.vstack((np.hstack((np.identity(r-1), np.zeros((r-1, n - (r-1))))), np.hstack((np.zeros((m-r+1, n - (m-r+1))), u))))
z1 = np.zeros((r-1, n - (r-1)))
z2 = np.zeros((n - (r-1), r-1))
big_new_u = np.vstack((np.hstack((np.identity(r-1), z1)), np.hstack((z2, new_u))))
return big_new_u
a = np.array([[1,2,3], [1,5,6], [1,8,9], [1,11,12]], dtype=float)
m = len(a)
n = len(a[0])
x = np.array([[1],[1],[1]])
b = np.dot(a, x)
p, a, rank = qrColumnPivoting(np.arange(n), a)
q = makeQ(a, rank)
a = cleanA(a)
u = fod(a, rank)
print "A:"
pprint.pprint(a)
print "P:"
pprint.pprint(p)
print "Q:"
pprint.pprint(q)
print "Q * Qt:"
pprint.pprint(np.dot(q, q.transpose()))
print "U:"
pprint.pprint(u)
print "U * Ut:"
pprint.pprint(np.dot(u, u.transpose()))
a1a2 = np.dot(q.transpose(), b)
pprint.pprint(a1a2)
y = np.vstack((np.dot(a[:rank, :rank].transpose(), a1a2[:rank]), np.zeros(n-rank)))
x = np.dot(u, y)
print "X:"
pprint.pprint(x)