python 矩阵运算

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

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

#!/usr/bin/env python
#coding = utf-8
'''
Author: Yang XU
E-mail: xuy1202@gmail.com
'''

import copy
import itertools
import types

class 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的行列式会被求成一个很小的数，有时间再解决