# QR decomposition (Givens)

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149``` ```import numpy as np import math from numpy import dot from scipy.linalg import norm np.set_printoptions(suppress=True, linewidth=160, precision=15) def transposition(n, i, j): trans = np.identity(n) trans[[i, j]] = trans[[j, i]] return trans def givens(n, i, j, x, y): g_matrix = np.identity(n) square = math.sqrt(x ** 2 + y ** 2) cosinus = 1.0 * x / square sinus = 1.0 * y / square g_matrix[i, i] = cosinus g_matrix[j, j] = cosinus g_matrix[i, j] = sinus g_matrix[j, i] = -sinus #g_matrix[[i, j], [i, j]] = np.array([[math.cos(alpha), math.sin(alpha)],[-math.sin(alpha), math.cos(alpha)]]) return g_matrix def givens_2(x, y): g_matrix = np.zeros([2, 2]) square = math.sqrt(x ** 2 + y ** 2) cosinus = 1.0 * x / square sinus = 1.0 * y / square g_matrix[0, 0] = cosinus g_matrix[1, 1] = cosinus g_matrix[0, 1] = sinus g_matrix[1, 0] = -sinus return g_matrix def max_column_index(A): return np.argmax(np.fromiter((norm(col) for col in A.T), dtype=np.float)) def qr_choise(A, tolerance=1e-16): m, n = A.shape Q = np.identity(m) P = np.identity(n) MIN_SIZE = min(m, n) rank = MIN_SIZE for j in xrange(MIN_SIZE): max_index = max_column_index(A[j:, j:]) + j if norm(A[j:, j:].T[max_index - j]) < tolerance: print 'Break at rank', str(j-1) rank = j - 1 break trans = transposition(n, j, max_index) A[:] = dot(A, trans) P[:] = dot(trans, P) for i in xrange(m - 1, j, -1): givens2 = givens_2(A[j, j], A[i, j]) A[[j, i], j:n] = dot(givens2, A[[j, i], j:n]) Q[[j, i], :] = dot(givens2, Q[[j, i], :]) return Q, P, rank def qr_simple(A): m, n = A.shape Q = np.identity(m) MIN_SIZE = min(m, n) for j in xrange(MIN_SIZE): for i in xrange(m-1, j, -1): givens2 = givens_2(A[j, j], A[i, j]) A[[j, i], j:n] = dot(givens2, A[[j, i], j:n]) Q[[j, i], :] = dot(givens2, Q[[j, i], :]) return Q # From citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.379.8623&rep=rep1&type=pdf def practical_qr(A, tolerance=1e-16): m, n = A.shape Q = np.identity(m) P = np.identity(n) MIN_SIZE = min(m, n) matrix_norm = norm(A) # make column norms norms = np.fromiter((norm(col)**2 for col in A.T), dtype=np.float) for j in xrange(MIN_SIZE): max_index = np.argmax(norms[j:]) if norms[max_index] < tolerance * matrix_norm: print 'Break at rank', j-1 break trans = transposition(n, j, max_index + j) norms[[j, max_index + j]] = norms[[max_index + j, j]] A[:] = dot(A, trans) P[:] = dot(trans, P) for i in xrange(m - 1, j, -1): givens2 = givens_2(A[j, j], A[i, j]) giv = givens(m, i, j, A[j, j], A[i, j]) A[[j, i], j:n] = dot(givens2, A[[j, i], j:n]) Q = dot(Q, giv) norms[j+1:] -= A[j, j+1: n]**2 return Q, P A = np.random.randn(5, 6) A[:, 4] = A[:, 1] + 5*A[:, 2] A[:, 3] = 3*A[:, 1] A[:, 5] = 2*A[:, 1] + 4*A[:, 2] A[:, 0] = 6*A[:, 2] + 5*A[:, 1] print 'matrix rank of A:' print np.linalg.matrix_rank(A) #A = np.loadtxt('matrix.txt') M = A.copy() Q, P, rank = qr_choise(A) R = A.copy() U = qr_simple(A[:rank, :].T) T_full = np.zeros(M.shape) T_full[:rank, :] = A[:rank, :] print 'R matrix with rank ', rank, ':' print R print 'T matrix:' print T_full print 'U:' print U print 'UR^t:' print U.shape, R.T.shape print dot(U.T, R.T) print '*'*80 print 'Q matrix:' print Q print '*'*80 print 'Prove that Q is orthogonal : Q*Q^t' print dot(Q, Q.T) print '*'*80 print 'permutation matrix P:' print P print '*'*80 print 'Matrix A:' print M print 'QR:' print dot(dot(Q, R), P) print 'Difference:' print dot(dot(Q.T, R), P) - M print dot(dot(Q.T, dot(T_full, U)), P) - M ```