文章目录

矩阵
- __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实现矩阵