голубков 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],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)