Householder 变换与 QR 分解

import random
import copy

EPS = 0.00001

class MatrixException( Exception ):
pass

class Matrix( object ):

def __init__( self, rows, cols, values_list = None, description = None ):
self.rows = rows
self.cols = cols
self.matrix = [ [ 0 for c in xrange( cols ) ] for r in xrange( rows ) ]
if values_list:
for r in xrange( rows ):
for c in xrange( cols ):
try:
self.matrix[r][c] = values_list[cols * r + c]
except IndexError:
self.matrix[r][c] = 0
if description:
self.description = description

def __getitem__( self, index ):
#以后写
if isinstance( index, str ):
pass
else:
return self.matrix[index]

def __repr__( self ):
for r in xrange( self.rows ):
print self.matrix[r]
return ''

def transpose( self ):
M = Matrix( self.cols, self.rows )
for c in xrange( self.cols ):
for r in xrange( self.rows ):
M[c][r] = self.matrix[r][c]
return M

@property
def T( self ):
return self.transpose()

@staticmethod
def mul( M, factor ):
res = Matrix( M.rows, M.cols )
for r in xrange( M.rows ):
for c in xrange( M.cols ):
res[r][c] = M[r][c] * factor
res = Matrix.format_matrix( res )
return res

def __sub__( self, other ):
try:
if self.rows != other.rows or self.cols != other.cols:
raise MatrixException
M = Matrix( self.rows, self.cols )
for r in xrange( self.rows ):
for c in xrange( self.cols ):
M[r][c] = self[r][c] - other[r][c]
M = Matrix.format_matrix( M )
return M
except MatrixException:
print "two matrix must be the same size."

def __rmul__( self, other ):
try:
if not isinstance( other, float ):
raise TypeError
return Matrix.mul( self, other )
except TypeError:
print "factor must be float."

def __mul__( self, other ):
if isinstance( other, float ):
return Matrix.mul( self, other )
elif isinstance( other, Matrix ):
try:
if self.cols != other.rows:
raise MatrixException
M = Matrix( self.rows, other.cols )
for r in xrange( self.rows ):
for c in xrange( other.cols ):
M[r][c] = cross_product( self[r], map( lambda x: x[c], other ) )
M = Matrix.format_matrix( M )
return M
except MatrixException:
print "self.cols must be equal to other.rows"
else:
pass

@staticmethod
def format_matrix( M ):
global EPS
for r in xrange( M.rows ):
for c in xrange( M.cols ):
if abs( M[r][c] ) < EPS:
M[r][c] = 0.0
else:
M[r][c] = round( M[r][c], 7 )
return M

def random_create( self, min = 1, max = 20 ):
for r in xrange( self.rows ):
for c in xrange( self.cols ):
self[r][c] = random.randint( min, max )

def get_col( self, col ):
M = Matrix( self.rows, 1 )
for r in xrange( self.rows ):
M[r][0] = self[r][col]
return M

def QR( self ):
R = copy.deepcopy( self )
Q = IdentityMatrix( self.rows )
for c in xrange( self.cols ):
col_matrix = R.get_col( c )
H = create_householder_matrix_by_cleared_to_zero( col_matrix, c )
Q = H * Q
R = H * R
return Q , R

class IdentityMatrix( Matrix ):
def __init__( self, size = 1 ):
self.matrix = [ [ 1 if c == r else 0 for c in xrange( size ) ] for r in xrange( size ) ]
self.rows = self.cols = size

def __sub__( self, other ):
try:
if other.rows != other.cols:
raise TypeError
self.matrix = [ [ 1 if c == r else 0 for c in xrange( other.rows ) ] for r in xrange( other.rows ) ]
self.rows = self.cols = other.rows
return Matrix.__sub__( self, other )
except TypeError:
print "matrix.rows must be equal to matrix.cols."

def cross_product( a, b ):
return reduce( lambda x, y : x + y, map( lambda x: x[0] * x[1], zip( a, b ) ) )

def create_householder_matrix_by_mirror_base( vector, mirror_base ):
global I
cut_index = mirror_base.index( 1 )
sign = -1 if vector[cut_index] < 0 else 1
sigma = sign * ( cross_product( vector, vector ) ) ** 0.5
U = vector[:cut_index] + [sigma + vector[cut_index]] + vector[cut_index + 1:]
U = Matrix( len( U ), 1, U )
rho = 0.5 * ( U.transpose() * U )[0][0]
H = I - ( 1.0 / rho ) * ( U * U.transpose() )
return H

def create_householder_matrix_by_cleared_to_zero( vector, cut_index ):
if isinstance( vector, Matrix ):
res = []
for r in xrange( vector.rows ):
res.append( vector[r][0] )
vector = res
global I
sign = -1 if vector[cut_index] < 0 else 1
sigma = sign * ( cross_product( vector[cut_index:],
vector[cut_index:] ) ) ** 0.5
mirror = vector[:cut_index] + [-sigma] + [0] * ( len( vector ) - cut_index - 1 )
U = map( lambda x: x[0] - x[1], zip( vector, mirror ) )
U = Matrix( len( U ), 1, U )
rho = 0.5 * ( U.transpose() * U )[0][0]
H = I - ( 1.0 / rho ) * ( U * U.transpose() )
return H

I = IdentityMatrix()
M1 = Matrix( 5, 5 )
M1.random_create()
Q, R = M1.QR()
print "Q:\n", Q
print "R:\n", R