使用python实现矩阵

文章目录

• 矩阵
• • __init__
  • __getitem__
  • __setitem__
  • reshape
  • __repr__
  • __add__ 与 __mul__
  • __matmul__
  • LU分解
  • 转置
  • 利用LU分解求行列式
• 利用LU分解解线性方程组

矩阵

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

init

from array import arrayclass 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))

来源地址：https://blog.csdn.net/A2233776/article/details/129492455