# голубков 1.1

 ``` 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``` ```__author__ = 'Lenovo' import pprint import scipy import scipy.linalg # SciPy Linear Algebra Library A = scipy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]]) Q, R = scipy.linalg.qr(A) print "A:" pprint.pprint(A) print "Q:" pprint.pprint(Q) print "R:" pprint.pprint(R) from math import sqrt from pprint import pprint def mult_matrix(M, N): """Multiply square matrices of same dimension M and N""" # Converts N into a list of tuples of columns tuple_N = zip(*N) # Nested list comprehension to calculate matrix multiplication return [[sum(el_m * el_n for el_m, el_n in zip(row_m, col_n)) for col_n in tuple_N] for row_m in M] def trans_matrix(M): """Take the transpose of a matrix.""" n = len(M) return [[ M[i][j] for i in range(n)] for j in range(n)] def norm(x): """Return the Euclidean norm of the vector x.""" return sqrt(sum([x_i**2 for x_i in x])) def Q_i(Q_min, i, j, k): """Construct the Q_t matrix by left-top padding the matrix Q with elements from the identity matrix.""" if i < k or j < k: return float(i == j) else: return Q_min[i-k][j-k] def householder(A): """Performs a Householder Reflections based QR Decomposition of the matrix A. The function returns Q, an orthogonal matrix and R, an upper triangular matrix such that A = QR.""" n = len(A) # Set R equal to A, and create Q as a zero matrix of the same size R = A Q = [[0.0] * n for i in xrange(n)] # The Householder procedure for k in range(n-1): # We don't perform the procedure on a 1x1 matrix, so we reduce the index by 1 # Create identity matrix of same size as A I = [[float(i == j) for i in xrange(n)] for j in xrange(n)] # Create the vectors x, e and the scalar alpha # Python does not have a sgn function, so we use cmp instead x = [row[k] for row in R[k:]] e = [row[k] for row in I[k:]] alpha = -cmp(x,0) * norm(x) # Using anonymous functions, we create u and v u = map(lambda p,q: p + alpha * q, x, e) norm_u = norm(u) v = map(lambda p: p/norm_u, u) # Create the Q minor matrix Q_min = [ [float(i==j) - 2.0 * v[i] * v[j] for i in xrange(n-k)] for j in xrange(n-k) ] # "Pad out" the Q minor matrix with elements from the identity Q_t = [[ Q_i(Q_min,i,j,k) for i in xrange(n)] for j in xrange(n)] # If this is the first run through, right multiply by A, # else right multiply by Q if k == 0: Q = Q_t R = mult_matrix(Q_t,A) else: Q = mult_matrix(Q_t,Q) R = mult_matrix(Q_t,R) # Since Q is defined as the product of transposes of Q_t, # we need to take the transpose upon returning it return trans_matrix(Q), R A = [[12, -51, 4], [6, 167, -68], [-4, 24, -41]] Q, R = householder(A) print "A:" pprint(A) print "Q:" pprint(Q) print "R:" pprint(R) ```