QR 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``` ```import numpy as np import math from numpy import dot np.set_printoptions(suppress=True, linewidth=160, precision=15) 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 givens_qr(A): m, n = A.shape Q = np.identity(m) for j in xrange(n): 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) return Q def givens_qr(A): m, n = A.shape Q = np.identity(m) P = np.identity(n) for j in xrange(n): 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) return Q A = np.random.randn(5, 5) #A[-1] = A[-2] #A[:, -1] = A[:, -1] #A = np.loadtxt('matrix.txt') M = A.copy() Q = givens_qr(A) print 'R matrix:' print A 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 'Matrix A:' print M print 'QR:' print dot(Q, A) print 'Difference:' print dot(Q, A) - M ```