python实现矩阵的示例代码

目录

• 矩阵
• init
• getitem
• setitem
• reshape
• repr
• add 与 mul
• matmul
• LU分解
• 转置
• 利用LU分解求行列式
• 利用LU分解解线性方程组

矩阵

使用python构建一个类，模拟矩阵，可以进行各种矩阵的计算，与各种方便的用法

init

from array import array
class Matrix:
    def __init__(self, matrix: 'a list of one dimension', shape: 'a tuple of shape' = None, 
                 dtype: 'data type code' = 'd'):
        # matrix一个包含所有元素的列表，shape指定形状，默认为列向量，dtype是数据类型
        # 使用一个数组模拟矩阵，通过操作这个数组完成矩阵的运算
        self.shape = (len(matrix), 1)
        if shape:
            self.shape = shape
        self.array = array(dtype, matrix)

getitem

由于矩阵是一个二维数组，应当支持诸如matrix[1, 2]，matrix[1:3, 2]，matrix[1:3, 2:4]之类的取值
所以我们需要使用slice类的indice方法实现__getitem__，并支持切片

    def __getitem__(self, item: 'a index of two dimensions'):
        # 使用slice类的indices方法，实现二维切片
        rows, cols = item
        # 下面对传入的指针或者切片进行处理，使其可以统一处理
        if isinstance(rows, slice):
            rows = rows.indices(self.shape[0])
        else:
            rows = (rows, rows + 1)
        if isinstance(cols, slice):
            cols = cols.indices(self.shape[1])
        else:
            cols = (cols, cols + 1)
        res = []
        shape = (len(range(*rows)), len(range(*cols)))  # 新矩阵的形状
        # 通过遍历按照顺序将元素加入新的矩阵
        for row in range(*rows):
            for col in range(*cols):
                index = row * self.shape[1] + col
                res.append(self.array[index])
        if len(res) == 1:
        # 若是数则返回数
            return res[0]
        # 若是矩阵则返回矩阵
        return Matrix(res, shape)

由于要支持切片，所以需要在方法中新建一个矩阵用于返回值

setitem

由于没有序列协议的支持，我们需要自己实现__setitem__

    def __setitem__(self, key: 'a index or slice of two dimensions', value):
        # 使用slice类的indices方法，实现二维切片的广播赋值
        rows, cols = key
        if isinstance(rows, slice):
            rows = rows.indices(self.shape[0])
        else:
            rows = (rows, rows + 1)
        if isinstance(cols, slice):
            cols = cols.indices(self.shape[1])
        else:
            cols = (cols, cols + 1)
        if isinstance(value, Matrix):
            # 对于传入的值是矩阵，则需要判断形状
            if value.shape != (len(range(*rows)), len(range(*cols))):
                raise ShapeError
            # 使用x，y指针取出value中的值赋给矩阵
            x = -1
            for row in range(*rows):
                x += 1
                y = -1
                for col in range(*cols):
                    y += 1
                    index = row * self.shape[1] + col
                    self.array[index] = value[x, y]
        else:
            for row in range(*rows):
                for col in range(*cols):
                    index = row * self.shape[1] + col
                    self.array[index] = value

若传入的value是一个数，这里的逻辑基本与__getitem__相同，实现了广播。

而若传入的是一个矩阵，则需要判断形状，对对应的元素进行赋值，这是为了方便LU分解。

reshape

reshape用于改变形状，对于上面的实现方法，只需要改变matrix.shape就可以了
注意改变前后的总元素数应当一致

    def reshape(self, shape: 'a tuple of shape'):
        if self.shape[0] * self.shape[1] != shape[0] * shape[1]:
            raise ShapeError
        self.shape = shape

repr

实现__repr__方法，较为美观的打印矩阵

    def __repr__(self):
        shape = self.shape
        _array = self.array
        return "[" + ",\n".join(str(list(_array[i * shape[1]:(i + 1) * shape[1]])) for i in range(shape[0])) + "]"

add 与 mul

对于加法与乘法的支持，这里的乘法是元素的乘法，不是矩阵的乘法
同样的，实现广播

    def __add__(self, other):
        shape = self.shape
        res = zeros(shape)  # 创建一个新的零矩阵，用于返回
        if isinstance(other, Matrix):
            # 实现同样形状的矩阵元素之间的加法
            if self.shape != other.shape:
                # 如果矩阵的形状对不上，就返回错误
                raise ShapeError
            for i in range(shape[0]):
                for j in range(shape[1]):
                    res[i, j] = self[i, j] + other[i, j]
        else:
            # 实现广播
            for i in range(shape[0]):
                for j in range(shape[1]):
                    res[i, j] = self[i, j] + other
        return res
    def __mul__(self, other):
        shape = self.shape
        res = zeros(shape)  # 创建一个新的零矩阵，用于返回
        if isinstance(other, Matrix):
            # 实现同样形状的矩阵元素之间的乘法
            if self.shape != other.shape:
                # 如果矩阵的形状对不上，就返回错误
                raise ShapeError
            for i in range(shape[0]):
                for j in range(shape[1]):
                    res[i, j] = self[i, j] * other[i, j]
        else:
            # 实现广播
            for i in range(shape[0]):
                for j in range(shape[1]):
                    res[i, j] = self[i, j] * other
        return res

matmul

matmul矩阵乘法，运算符为@

    def __matmul__(self, other):
        # 实现矩阵的乘法
        if self.shape[1] != other.shape[0]:
            # 对形状进行判断
            raise ShapeError
        if self.shape[0] == 1 and other.shape[1] == 1:
            # 行向量与列向量的乘积，就是它们的数量积
            length = self.shape[1]
            return sum(self[0, i] * self[i, 0] for i in range(length))
        res = []
        shape = (self.shape[0], other.shape[1])
        for i in range(shape[0]):
            for j in range(shape[1]):
                # 将两个矩阵分别按行向量与列向量分块，然后相乘
                try:
                    # 由于切片返回的可能是数，而数不支持'@'运算符，所以使用异常处理语句
                    res.append(self[i, :] @ other[:, j])
                except TypeError:
                    res.append(self[i, :] * other[:, j])
        return Matrix(res, shape)

将矩阵分成向量进行矩阵乘法

LU分解

用属性self._lu，构建一个新的矩阵，作为LU分解表。self._lu并不会在初始化时创建，而是在需要用到LU分解表时计算。
同时，我们维护一个self.changed属性，用来判断在需要用到LU分解表时是否需要重新进行LU分解

    def __init__(self, matrix: 'a list of one dimension', shape: 'a tuple of shape' = None,
                 dtype: 'data type code' = "d"):
        # matrix一个包含所有元素的列表，shape指定形状默认为列向量，dtype是数据类型
        # 使用一个数组模拟矩阵，通过操作这个数组完成矩阵的运算
        self.shape = (len(matrix), 1)
        if shape:
            self.shape = shape
        self.array = array(dtype, matrix)
        self._changed = True
        self._primary = list(range(shape[0])) # 只进行行交换

显然，当我们修改了矩阵的元素后就需要重新进行LU分解，重写 setitem ，在开始时修改self.changed属性

    def __setitem__(self, key: 'a index or slice of two dimensions', value):
        # 使用slice类的indices方法，实现二维切片的广播赋值
        self._change = True
        ...

在lu分解中需要选择主元，用属性self._primary储存当前矩阵的主元表，因此我们要重写 getitem 与 setitem

...
row = self._primary[row]  # 通过主元表进行赋值，将行换为
index = row * self.shape[1] + col
...

下面，我们来实现LU分解

    def _primary_update(self, k):
        # 选择绝对值最大的数作为主元
        max_val = -777
        max_index = 0
        for i in range(k, self.shape[0]):
            x = abs(self[i, k])
            if x > max_val:
                max_val = x
                max_index = i
        self._primary[k], self._primary[max_index] = self._primary[max_index], self._primary[k]
    def _lu_factorization(self):
        self._lu = Matrix(self.array, self.shape) # 新建一个矩阵储存LU分解表
        rows, cols = self.shape
        _lu = self._lu
        step = min(rows, cols)
        for k in range(step):
            if _lu[k, k] == 0:
                # 如果当前对角元素为0，就需要更换主元
                _lu._primary_update(k)
            if _lu[k, k] == 0:
                # 如果更换主元之后仍然为0，就说明该列全为0,跳过
                break
            x = 1 / _lu[k, k]
            _lu[k + 1:, k] *= x
            for i in range(k + 1, rows):
                for j in range(k + 1, cols):
                    _lu[i, j] = _lu[i, j] - _lu[i, k] * _lu[k, j]

转置

用一个方法实现转置，而不是维护一个属性

    def trans(self):
        shape = self.shape[::-1]
        res = zeros(shape) # 创建一个零矩阵用于返回
        for i in range(shape[0]):
            for j in range(shape[1]):
                res[i, j] = self[j, i]
        return res

利用LU分解求行列式

原矩阵的行列式就是L与U的对角元素的乘积

    def det(self):
        if self.shape[0] != self.shape[1]:
            raise ShapeError
        if self._changed:
            self._lu_factorization()
            self._changed = False
        res = 1
        for i in range(self.shape[0]):
            res *= self[i, i]
        return res

利用LU分解解线性方程组

利用LU分解可以快速地解出线性方程组

    def linear_equation(self, y):
        # 利用LU分解表解方程
        if not self.det():
            # 不考虑扁平化的情况，即使可能有解
            raise DetError
        lu = self._lu
        length = self.shape[1]
        z = [0]*length # 先解 L @ z = y
        for i in range(length):
            z_i = y[i, 0]
            for j in range(i):
                z_i -= z[j] * lu[i, j]
            z[i] = z_i
        x = [0]*length # 再解 U @ x = z
        for i in range(length - 1, -1, -1):
            x_i = z[i]
            for j in range(length - 1, i, -1):
                x_i -= x[j] * lu[i, j]
            x[i] = x_i / lu[i, i]
        return Matrix(x, (length, 1))

到此这篇关于python实现矩阵的示例代码的文章就介绍到这了,更多相关python 矩阵内容请搜索自由互联以前的文章或继续浏览下面的相关文章希望大家以后多多支持自由互联！