数值分析--基本矩阵类

为了方便后面数值分析的学习，先定义一些矩阵的基本操作

Matrix.py

import copy

class Matrix:
    def __init__(self, matrix):
        # 列表类型检查
        assert isinstance(matrix, list), "参数必须为二维列表！"
        for row in matrix:
            assert isinstance(row, list), "参数必须为二维列表！"
        for row in matrix:
            assert len(row) == len(matrix[0]), "参数第二维度不一致！"
        for row in matrix:
            for each in row:
                assert isinstance(each, (int, float, complex)), "元素类型只能为int, float或complex！"
        self.__matrix = matrix
        self.__row = len(matrix)
        self.__col = len(matrix[0])

    # 矩阵间加法 +
    def __add__(self, other):
        assert isinstance(other, Matrix), "加数必须是 Matrix 类"
        if self.__row != other.__row or self.__col != other.__col:
                raise ValueError("两矩阵维度不一致！")
        result = Matrix( [[0 for i in range(self.__col)]  for i in range( self.__row ) ] )
        for i in range(self.__row):
            for j in range(self.__col):
                result.__matrix[i][j] = self.__matrix[i][j] + other.__matrix[i][j]
        return result

    # 矩阵间减法 -
    def __sub__(self, other):
        assert isinstance(other, Matrix), "加数必须是 Matrix 类"
        if self.__row != other.__row or self.__col != other.__col:
            raise ValueError("两矩阵维度不一致！")
        result = Matrix([[0 for i in range(self.__col)] for i in range(self.__row)])
        for i in range(self.__row):
            for j in range(self.__col):
                result.__matrix[i][j] = self.__matrix[i][j] - other.__matrix[i][j]
        return result
    # 矩阵乘矩阵 *
    def __mul__(self, other):
        assert isinstance(other, Matrix), "右乘数必须是 Matrix 类"
        if self.__col != other.__row:
            raise ValueError("两矩阵维度不一致！")
        result = Matrix([[0 for i in range(other.__col)] for i in range(self.__row)])
        for i in range(self.__row):
            for j in range(other.__col):
                for m in range(self.__col):
                    result.__matrix[i][j] += self.__matrix[i][m] * other.__matrix[m][j]
        return result
    # 数与矩阵乘 *
    def __rmul__(self, other):
        assert isinstance(other, (int, float, complex)), "左乘系数类型错误"
        result = Matrix([[0 for i in range(self.__col)] for i in range(self.__row)])
        for i in range(self.__row):
            for j in range(self.__col):
                result.__matrix[i][j] = self.__matrix[i][j] * other
        return result
    # 矩阵除以数
    def __truediv__(self, other):
        assert isinstance(other, (int, float, complex)), "被除数系数类型错误"
        result = Matrix([[0 for i in range(self.__col)] for i in range(self.__row)])
        for i in range(self.__row):
            for j in range(self.__col):
                result.__matrix[i][j] = self.__matrix[i][j] / other
        return result
    # 数除以矩阵
    def __rtruediv__(self, other):
        assert isinstance(other, (int, float, complex)), "除数系数类型错误"
        result = Matrix([[0 for i in range(self.__col)] for i in range(self.__row)])
        for i in range(self.__row):
            for j in range(self.__col):
                result.__matrix[i][j] = other / self.__matrix[i][j]
        return result
    # 矩阵与矩阵元素是否相等
    def __eq__(self, other):
        assert isinstance(other, Matrix), "必须是 Matrix 类"
        for i in range(self.__row):
            for j in range(self.__col):
                if self.__matrix[i][j] != other.__matrix[i][j]:
                    return False
        return True
    # 矩阵字符化
    def __str__(self):
        return '\n'.join(list(map(str, self.__matrix)))
    # 行迭代器
    def __iter__(self):
        return self.__matrix.__iter__()
    # 元素获取，看齐numpy
    def __getitem__(self, item):
        if isinstance(item, int):
            return self.__matrix[item]
        if isinstance(item, tuple):
            row, col = item[0], item[1]
            if isinstance(row, int) and isinstance(col, int):
                return self.__matrix[row][col]
            if isinstance(row, slice) and isinstance(col, int):
                return Matrix( [ [ line[col] ] for line in self.__matrix[row]] )
            if isinstance(row, int) and isinstance(col, slice):
                return Matrix( [self.__matrix[row][col]] )
            if isinstance(row, slice) and isinstance(col, slice):
                return Matrix( [ line[col] for line in self.__matrix[row] ] )
    # len()
    def __len__(self):
        return (self.__row)*(self.__col)

    # 析构函数
    def __del__(self):
        del self.__matrix
    # def __copy__(self):
    # 元素和
    def sum(self):
        value = 0
        for i in range(self.__row):
            for j in range(self.__col):
                value += self.__matrix[i][j]
        return value

    # 某行某列的代数余子式
    def cofacor(self, row, col):
        assert self.__row == self.__col, "必须是方阵"
        mat_1 = copy.deepcopy(self)
        temp_mat = []
        # 删除对应行和列
        for row_index, line in enumerate(mat_1):
            if row_index == row:
                continue
            line.pop(col)
            temp_mat.append(line)
        return pow(-1, row+col+2) * Matrix(temp_mat).det()

    # 定义法迭代法行列式
    def det(self):
        assert  self.__row == self.__col, "必须是方阵"
        if self.__row == 1 and self.__col == 1:
            return self.__matrix[0][0]
        if self.__row == 2 and self.__col == 2:
            return self.__matrix[0][0]*self.__matrix[1][1] - self.__matrix[0][1]*self.__matrix[1][0]
        # 取出第一行与之代数余子式相乘
        value = 0
        for col_index, m in enumerate(self.__matrix[0]):
            value += m * copy.deepcopy(self).cofacor(0, col_index)
        return value
    # 当前矩阵的伴随矩阵
    def adj(self):
        assert self.__row == self.__col, "必须是方阵"
        result = [[0 for i in range(self.__col)] for i in range(self.__row)]
        for i in range(self.__row):
            for j in range(self.__col):
                result[i][j] = self.cofacor(j,i)
        return Matrix(result)
    # 矩阵的逆
    def inv(self):
        assert self.__row == self.__col, "必须是方阵"
        det = self.det()
        assert det != 0, "矩阵必须非奇异"
        return self.adj() / det
    # 矩阵转置
    def transpose(self):
        return Matrix(list(map(list, zip(*self.__matrix))))

    # 矩阵的秩
    def rank(self):
        def col_sort(mat, col_index = 0):
            assert isinstance(mat, Matrix), "对象必须为矩阵"
            temp_mat = Matrix( sorted(mat.__matrix, key=lambda x:abs(x[col_index]), reverse=True) )
            return temp_mat
        def cal_col_k(mat, col_index=0):
            assert isinstance(mat, Matrix), "对象必须为矩阵"
            k = [[0] for i in range(mat.__row)]
            if abs(mat.__matrix[0][col_index]) < 1e-9:
                return 0, None
            for i in range(1,mat.__row):
                k[i][0] = mat.__matrix[i][col_index] / mat.__matrix[0][col_index]
            return 1, Matrix(k)
        def elimination(mat, k):
            assert isinstance(mat, Matrix), "对象必须为矩阵"
            temp_mat = [[0 for col in range(mat.__col)] for j in range(mat.__row) ]
            for i in range(mat.__row):
                for j in range(mat.__col):
                    temp_mat[i][j] = mat.__matrix[i][j] -( mat.__matrix[0][j] * k[i][0])
            return Matrix(temp_mat)
        
        mat = copy.deepcopy(self)
        # 列排序
        mat = col_sort(mat, col_index=0)
        if mat.__col == 1 or mat.__row == 1:
                return int(bool(abs(mat.sum())))
        # 列比例
        ret, K = cal_col_k(mat, col_index=0)
        if ret == 0:
            return mat[1:,1:].rank()
        # 消源
        mat = elimination(mat, K)
        return 1 + mat[1:,1:].rank()

测试main.py

from Matrix import Matrix

A = Matrix( [ [1, 2, 2],
              [2, 3, 2],
              [2, 2, 3]])

M = Matrix([    [0],
                [0],
                [2]])

B = Matrix( [ [1, 2, 3],
              [3, 2, -1],
              [2, 2, 2]] )

# 加法
print("----矩阵加法----")
print(A + B)
print("----矩阵减法----")
print(A - B)
print("----矩阵乘法----")
print(A * B)
print("----数乘矩阵----")
print(10 * B)
print("----矩阵除以数----")
print(A / 10)
print("----数除以矩阵（元素）----")
print(10 / B)
print("----矩阵元素访问(副本)----")
print(A[1:, 1:])

print("----代数余子式----")
print(A.cofacor(0, 0))
print("----行列式----")
print(A.det())
print("----伴随矩阵----")
print(A.adj())
print("----矩阵逆----")
print(A.inv())
print("----矩阵转置----")
print(A.transpose())
print("----矩阵秩----")
print(A.rank())


# 输出：
# ----矩阵加法----
# [2, 4, 5]
# [5, 5, 1]
# [4, 4, 5]
# ----矩阵减法----
# [0, 0, -1]
# [-1, 1, 3]
# [0, 0, 1]
# ----矩阵乘法----
# [11, 10, 5]
# [15, 14, 7]
# [14, 14, 10]
# ----数乘矩阵----
# [10, 20, 30]
# [30, 20, -10]
# [20, 20, 20]
# ----矩阵除以数----
# [0.1, 0.2, 0.2]
# [0.2, 0.3, 0.2]
# [0.2, 0.2, 0.3]
# ----数除以矩阵（元素）----
# [10.0, 5.0, 3.3333333333333335]
# [3.3333333333333335, 5.0, -10.0]
# [5.0, 5.0, 5.0]
# ----矩阵元素访问----
# [3, 2]
# [2, 3]
# ----代数余子式----
# 5
# ----行列式----
# -3
# ----伴随矩阵----
# [5, -2, -2]
# [-2, -1, 2]
# [-2, 2, -1]
# ----矩阵逆----
# [-1.6666666666666667, 0.6666666666666666, 0.6666666666666666]
# [0.6666666666666666, 0.3333333333333333, -0.6666666666666666]
# [0.6666666666666666, -0.6666666666666666, 0.3333333333333333]
# ----矩阵转置----
# [1, 2, 2]
# [2, 3, 2]
# [2, 2, 3]
# ----矩阵秩----
# 3