python 矩阵运算

由于自己基本功不扎实且遗忘，上一篇《python实现求行列式的值》成功出错，其计算的有效性只限于2，3维。尽管我对之前所有数学老师的填鸭式教学报以仇恨式的埋怨，但也对自己的挫深表羞愧...

下面脚本修复了之前求行列式的错误，并丰富了其他的矩阵运算的基本内容，包括求常用的乘法及逆矩阵等。

#!/usr/bin/env python#coding = utf-8'''Author: Yang XUE-mail: [e-mail]'''import copyimport itertoolsimport typesclass Matrix(list):    def __init__(self, values):        if not isinstance(values, (types.ListType, Matrix)) or \           not isinstance(values[0], (types.ListType, Matrix)):            raise ValueError('Unsuported values for Matrix: %s'%str(values))        super(Matrix, self).__init__(values)        def shape(self):        i_length = len(self)        j_length = len(self[0])        return (i_length, j_length)        def __str__(self):        i_length, j_length = self.shape()        if i_length > 10:            _M = self[:5] + self[-5:]            line_eclipse = True        else:            _M = self            line_eclipse = False        returnStr = '<Matrix(%s*%s): %s>\n'%(i_length, j_length, id(self))        count = 0        for _list in _M:            count += 1            if j_length > 8:                _tmpList = _list[:4] + ['...'] + _list[-4:]            else:                _tmpList = _list            returnStr += '%s\n'%str(_tmpList)            if line_eclipse and count == 5:                returnStr += '... ...\n'        return returnStr# 转置def transpose(M):    _tmp = zip(*M)    return Matrix([list(_l) for _l in _tmp])# 阵乘def multiplyMatrix(M1, M2):    M1 = Matrix(M1)    M2 = Matrix(M2)    m, x1 = M1.shape()    x2, n = M2.shape()    if x1 != x2:        raise ValueError('Impossibal Multiplication for M(%s,%s) * M(%s,%s)'%(m,x1,x2,n))        productF = lambda item: item[0]*item[1]    M2 = transpose(M2)    returnMatrix = []    for l_list in M1:        _tmpList = []        returnMatrix.append(_tmpList)        for r_list in M2:            value = sum([ productF(item) for item in zip(l_list, r_list)])            if abs(round(value) - value) < 0.00001: value = int(round(value))            _tmpList.append(value)    return Matrix(returnMatrix)# 数乘def multiplyNumber(N, M):    _M = copy.deepcopy(M)    _N = float(N)    _M = [        [ value*_N for value in _l ]        for _l in _M    ]    return Matrix(_M)# 上三角阵def getUpperTriangularMatrix(M):    # 杜尔里特算法（Doolittle algorithm）    M = Matrix(M)    m, n = M.shape()    j_indexer = itertools.cycle(range(n))    minusF = lambda rate, item: item[0] + (rate * item[1])    for i in range(m):        j = j_indexer.next()        base = float(M[i][j])        # 逢对角线元素为0，将本行压入最底        _count = 0        while base == 0 and _count<(m-i):            _count += 1            M.append(M.pop(i))            base = float(M[i][j])        if base == 0: continue        # 初等行变化，将下三角设为0        _i = i        _j = j        while _i < m-1:            _i += 1            _base = float(M[_i][_j])            if _base == 0: continue            rate = -(_base/base)            M[_i] = [minusF(rate, item) for item in zip(M[_i], M[i])]    return M# 行列式def getDeterminant(M):    # 杜尔里特算法（Doolittle algorithm）    M = Matrix(M)    m, n = M.shape()    if m != n:        raise ValueError('Matrix(%s*%s) %s has NO Det.'%(m, n, M))    j_indexer = itertools.cycle(range(n))    minusF = lambda rate, item: item[0] + (rate * item[1])    product = 1    for i in range(m):        j = j_indexer.next()        base = float(M[i][j])        _count = 0        here_to_tail_span = m-i-1        while base == 0 and _count<here_to_tail_span:            product *= (-1)**here_to_tail_span            _count += 1            M.append(M.pop(i))            base = float(M[i][j])        if base == 0: return 0        product *= base        _i = i        _j = j        while _i < m-1:            _i += 1            _base = float(M[_i][_j])            if _base == 0: continue            rate = -(_base/base)            M[_i] = [minusF(rate, item) for item in zip(M[_i], M[i])]    return product# 伴随矩阵def getAdjugateMatrix(M):    length = len(M)    if length == 1: return [[1]]    _returnM = []    for i in range(length):        _tmpList = []        _returnM.append(_tmpList)        for j in range(length):            _M = copy.deepcopy(M)            _M = [                _l for _l in _M                if _M.index(_l) != i            ]            [ _l.pop(j) for _l in _M ]            _Determinant = getDeterminant(_M)            _tmpList.append(((-1)**(i+j))*_Determinant)    return Matrix(transpose(_returnM))# 逆矩阵def getInversedMatrix(M):    # A* / |A|    _Determinant = getDeterminant(M)    if _Determinant == 0:        raise ZeroDivisionError('%s is NOT Non-Singular, has NO InversedMatrix'%str(M))    k = 1.0/_Determinant    _AdjugateMatrix = getAdjugateMatrix(M)    _returnM = multiplyNumber(k, _AdjugateMatrix)    return _returnM

受限于python计算精度，仍有些小问题，比如应该是0的行列式会被求成一个很小的数，有时间再解决