python實(shí)現(xiàn)矩陣的示例代碼

更新時(shí)間：2023年07月07日 10:34:51 作者：星火流明

本文主要介紹了python實(shí)現(xiàn)矩陣的示例代碼，文中通過示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧

矩陣

使用python構(gòu)建一個(gè)類，模擬矩陣，可以進(jìn)行各種矩陣的計(jì)算，與各種方便的用法

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一個(gè)包含所有元素的列表，shape指定形狀，默認(rèn)為列向量，dtype是數(shù)據(jù)類型
        # 使用一個(gè)數(shù)組模擬矩陣，通過操作這個(gè)數(shù)組完成矩陣的運(yùn)算
        self.shape = (len(matrix), 1)
        if shape:
            self.shape = shape
        self.array = array(dtype, matrix)

getitem

由于矩陣是一個(gè)二維數(shù)組，應(yīng)當(dāng)支持諸如matrix[1, 2]，matrix[1:3, 2]，matrix[1:3, 2:4]之類的取值
所以我們需要使用slice類的indice方法實(shí)現(xiàn)__getitem__，并支持切片

    def __getitem__(self, item: 'a index of two dimensions'):
        # 使用slice類的indices方法，實(shí)現(xiàn)二維切片
        rows, cols = item
        # 下面對(duì)傳入的指針或者切片進(jìn)行處理，使其可以統(tǒng)一處理
        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)))  # 新矩陣的形狀
        # 通過遍歷按照順序?qū)⒃丶尤胄碌木仃?
        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:
        # 若是數(shù)則返回?cái)?shù)
            return res[0]
        # 若是矩陣則返回矩陣
        return Matrix(res, shape)

由于要支持切片，所以需要在方法中新建一個(gè)矩陣用于返回值

setitem

由于沒有序列協(xié)議的支持，我們需要自己實(shí)現(xiàn)__setitem__

    def __setitem__(self, key: 'a index or slice of two dimensions', value):
        # 使用slice類的indices方法，實(shí)現(xiàn)二維切片的廣播賦值
        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):
            # 對(duì)于傳入的值是矩陣，則需要判斷形狀
            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是一個(gè)數(shù)，這里的邏輯基本與__getitem__相同，實(shí)現(xiàn)了廣播。

而若傳入的是一個(gè)矩陣，則需要判斷形狀，對(duì)對(duì)應(yīng)的元素進(jìn)行賦值，這是為了方便LU分解。

reshape

reshape用于改變形狀，對(duì)于上面的實(shí)現(xiàn)方法，只需要改變matrix.shape就可以了
注意改變前后的總元素?cái)?shù)應(yīng)當(dāng)一致

    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

實(shí)現(xiàn)__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

對(duì)于加法與乘法的支持，這里的乘法是元素的乘法，不是矩陣的乘法
同樣的，實(shí)現(xiàn)廣播

? ? def __add__(self, other):
? ? ? ? shape = self.shape
? ? ? ? res = zeros(shape) ?# 創(chuàng)建一個(gè)新的零矩陣，用于返回
? ? ? ? if isinstance(other, Matrix):
? ? ? ? ? ? # 實(shí)現(xiàn)同樣形狀的矩陣元素之間的加法
? ? ? ? ? ? if self.shape != other.shape:
? ? ? ? ? ? ? ? # 如果矩陣的形狀對(duì)不上，就返回錯(cuò)誤
? ? ? ? ? ? ? ? raise ShapeError
? ? ? ? ? ? for i in range(shape[0]):
? ? ? ? ? ? ? ? for j in range(shape[1]):
? ? ? ? ? ? ? ? ? ? res[i, j] = self[i, j] + other[i, j]
? ? ? ? else:
? ? ? ? ? ? # 實(shí)現(xiàn)廣播
? ? ? ? ? ? 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) ?# 創(chuàng)建一個(gè)新的零矩陣，用于返回
? ? ? ? if isinstance(other, Matrix):
? ? ? ? ? ? # 實(shí)現(xiàn)同樣形狀的矩陣元素之間的乘法
? ? ? ? ? ? if self.shape != other.shape:
? ? ? ? ? ? ? ? # 如果矩陣的形狀對(duì)不上，就返回錯(cuò)誤
? ? ? ? ? ? ? ? raise ShapeError
? ? ? ? ? ? for i in range(shape[0]):
? ? ? ? ? ? ? ? for j in range(shape[1]):
? ? ? ? ? ? ? ? ? ? res[i, j] = self[i, j] * other[i, j]
? ? ? ? else:
? ? ? ? ? ? # 實(shí)現(xiàn)廣播
? ? ? ? ? ? for i in range(shape[0]):
? ? ? ? ? ? ? ? for j in range(shape[1]):
? ? ? ? ? ? ? ? ? ? res[i, j] = self[i, j] * other
? ? ? ? return res

matmul

matmul矩陣乘法，運(yùn)算符為@

    def __matmul__(self, other):
        # 實(shí)現(xiàn)矩陣的乘法
        if self.shape[1] != other.shape[0]:
            # 對(duì)形狀進(jìn)行判斷
            raise ShapeError
        if self.shape[0] == 1 and other.shape[1] == 1:
            # 行向量與列向量的乘積，就是它們的數(shù)量積
            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]):
                # 將兩個(gè)矩陣分別按行向量與列向量分塊，然后相乘
                try:
                    # 由于切片返回的可能是數(shù)，而數(shù)不支持'@'運(yùn)算符，所以使用異常處理語句
                    res.append(self[i, :] @ other[:, j])
                except TypeError:
                    res.append(self[i, :] * other[:, j])
        return Matrix(res, shape)

將矩陣分成向量進(jìn)行矩陣乘法

LU分解

用屬性self._lu，構(gòu)建一個(gè)新的矩陣，作為L(zhǎng)U分解表。self._lu并不會(huì)在初始化時(shí)創(chuàng)建，而是在需要用到LU分解表時(shí)計(jì)算。
同時(shí)，我們維護(hù)一個(gè)self.changed屬性，用來判斷在需要用到LU分解表時(shí)是否需要重新進(jìn)行LU分解

? ? def __init__(self, matrix: 'a list of one dimension', shape: 'a tuple of shape' = None,
? ? ? ? ? ? ? ? ?dtype: 'data type code' = "d"):
? ? ? ? # matrix一個(gè)包含所有元素的列表，shape指定形狀默認(rèn)為列向量，dtype是數(shù)據(jù)類型
? ? ? ? # 使用一個(gè)數(shù)組模擬矩陣，通過操作這個(gè)數(shù)組完成矩陣的運(yùn)算
? ? ? ? self.shape = (len(matrix), 1)
? ? ? ? if shape:
? ? ? ? ? ? self.shape = shape
? ? ? ? self.array = array(dtype, matrix)
? ? ? ? self._changed = True
? ? ? ? self._primary = list(range(shape[0])) # 只進(jìn)行行交換

顯然，當(dāng)我們修改了矩陣的元素后就需要重新進(jìn)行LU分解，重寫 setitem ，在開始時(shí)修改self.changed屬性

    def __setitem__(self, key: 'a index or slice of two dimensions', value):
        # 使用slice類的indices方法，實(shí)現(xiàn)二維切片的廣播賦值
        self._change = True
        ...

在lu分解中需要選擇主元，用屬性self._primary儲(chǔ)存當(dāng)前矩陣的主元表，因此我們要重寫 getitem 與 setitem

...
row = self._primary[row]  # 通過主元表進(jìn)行賦值，將行換為
index = row * self.shape[1] + col
...

下面，我們來實(shí)現(xiàn)LU分解

? ? def _primary_update(self, k):
? ? ? ? # 選擇絕對(duì)值最大的數(shù)作為主元
? ? ? ? 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) # 新建一個(gè)矩陣儲(chǔ)存LU分解表
? ? ? ? rows, cols = self.shape
? ? ? ? _lu = self._lu
? ? ? ? step = min(rows, cols)
? ? ? ? for k in range(step):
? ? ? ? ? ? if _lu[k, k] == 0:
? ? ? ? ? ? ? ? # 如果當(dāng)前對(duì)角元素為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]

轉(zhuǎn)置

用一個(gè)方法實(shí)現(xiàn)轉(zhuǎn)置，而不是維護(hù)一個(gè)屬性

    def trans(self):
        shape = self.shape[::-1]
        res = zeros(shape) # 創(chuàng)建一個(gè)零矩陣用于返回
        for i in range(shape[0]):
            for j in range(shape[1]):
                res[i, j] = self[j, i]
        return res

利用LU分解求行列式

原矩陣的行列式就是L與U的對(duì)角元素的乘積

    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))

到此這篇關(guān)于python實(shí)現(xiàn)矩陣的示例代碼的文章就介紹到這了,更多相關(guān)python 矩陣內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

您可能感興趣的文章: