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