__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) for i in range(0, r): v = house(rr[i:n, i]) #rr2 = np.copy(rr) rr[i:n, i:r] = rowHouse(rr[i:n, i:r], v) if i < n: rr[i+1:n, i] = v[1:n] # vt = np.array([v]).transpose() # vsq = np.dot(np.array([v]).transpose(), np.array([v])) # new_u = np.identity(n) - 2/sum(j ** 2 for j in v)*vsq #temp = np.vstack((np.hstack((np.identity(r-1), z1)), np.hstack((z2, new_u)))) #t2 = np.dot(temp, rr2) #u = np.dot(new_u, u) #big_new_u = np.vstack((np.hstack((np.identity(r-1), z1)), np.hstack((z2, new_u)))) return makeQ(rr, r) def makePMatr(p): 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) print "A:" pprint.pprint(a) q = makeQ(a, rank) print "Q:" pprint.pprint(q) a = cleanA(a) aold = np.copy(a) print "Qr:" pprint.pprint(np.dot(q, a)) ############################################### u = fod(a, rank) a = np.dot(u.transpose(), aold.transpose()) print "U * Aold:" pprint.pprint(aold) pprint.pprint(np.dot(u.transpose(), aold.transpose())) 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())) pmatr = makePMatr(p) 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)