import randomimport copyEPS = 0.00001class MatrixException( Exception ):passclass Matrix( object ):def __init__( self, rows, cols, values_list = None, description = None ):self.rows = rowsself.cols = colsself.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] = 0if description:self.description = descriptiondef __getitem__( self, index ):#以后写if isinstance( index, str ):passelse: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@propertydef T( self ):return self.transpose()@staticmethoddef 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] * factorres = Matrix.format_matrix( res )return resdef __sub__( self, other ):try:if self.rows != other.rows or self.cols != other.cols:raise MatrixExceptionM = 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 Mexcept MatrixException:print "two matrix must be the same size."def __rmul__( self, other ):try:if not isinstance( other, float ):raise TypeErrorreturn 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 MatrixExceptionM = 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 Mexcept MatrixException:print "self.cols must be equal to other.rows"else:pass@staticmethoddef format_matrix( M ):global EPSfor r in xrange( M.rows ):for c in xrange( M.cols ):if abs( M[r][c] ) < EPS:M[r][c] = 0.0else:M[r][c] = round( M[r][c], 7 )return Mdef 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 Mdef 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 * QR = H * Rreturn Q , Rclass 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 = sizedef __sub__( self, other ):try:if other.rows != other.cols:raise TypeErrorself.matrix = [ [ 1 if c == r else 0 for c in xrange( other.rows ) ] for r in xrange( other.rows ) ]self.rows = self.cols = other.rowsreturn 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 Icut_index = mirror_base.index( 1 )sign = -1 if vector[cut_index] < 0 else 1sigma = sign * ( cross_product( vector, vector ) ) ** 0.5U = 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 Hdef 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 = resglobal Isign = -1 if vector[cut_index] < 0 else 1sigma = sign * ( cross_product( vector[cut_index:],vector[cut_index:] ) ) ** 0.5mirror = 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 HI = IdentityMatrix()M1 = Matrix( 5, 5 )M1.random_create()Q, R = M1.QR()print "Q:\n", Qprint "R:\n", R
版权声明：本文为博主原创文章，未经博主允许不得转载。
，每个人在他的人生发轫之初，总有一段时光，
