4.3 Dodgson算法

醒过来摸鱼

已于 2023-05-19 14:29:05 修改

阅读量586

点赞数

分类专栏：线性代数【更新中】文章标签：算法线性代数矩阵

于 2023-01-22 14:22:44 首次发布

本文链接：https://blog.csdn.net/m0_66201040/article/details/128748796

版权

线性代数【更新中】专栏收录该内容

71 篇文章 38 订阅

订阅专栏

文章目录

算法过程
Python实现

算法过程

Dodgson不仅仅是一个数学家，其实他更著名的一个身份是作家，写出了著名的儿童文学作品《爱丽丝梦游仙境》（Alice in Wonderland）。他发明文学作品时用笔名Lewis Carroll，发表数学论文时用Dodgson，所以很难把这两个名字联系起来。
闲话不说了，先说说他的算法过程。总体来说，算法过程是一个类似于Chió算法的过程，把矩阵不断缩小，直到缩小到 $1\times 1$ 的矩阵。但是细节与Chió算法有很大不同,Dodgson算法是每次缩小两阶。算法具体是如下过程：

先通过行倍加，让中心块interior没有0，中心块是一个 $(n-2)\times (n-2)$ 的矩阵，这个中心块是原矩阵的子矩阵就叫做 $S$ 吧。
把每相邻四个元素的行列式计算出来，组成一个 $(n-1)\times (n-1)$ 的矩阵，记作矩阵 $B$ 。
还没完，把矩阵B的每相邻四个元素的行列式计算出来，这样组成一个 $(n-2)\times (n-2)$ 的矩阵，记作矩阵 $C$ ，把 $C$ 的每个元素除于A的中心块的对应位置元素，得到矩阵 $D$ 。
如果D是 $1\times 1$ 的矩阵，那么原矩阵的行列式就等于 $D$ 的行列式。如果不是，再进行一轮迭代，就计算B的行列式，但是需要跳过第2步的计算，用D作为第2步计算的结果。

我以一个 $5\times 5$ 例子讲述这个过程：
$\begin{pmatrix}-1 & 1 & -1 & 2 & 2\\ 1 & -1 & 2 & 2 & 3\\ 1 & 1 & 2 & -1 & -3\\ 2 & -1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1 & 1\\ \end{pmatrix}\\ B= \begin{pmatrix}0 & 1 & -6 & 2\\ 2 & -4 & -6 & -3\\ -3 & 3 & 3 & 2\\ 3 & -2 & 0 & 0\\ \end{pmatrix}\\ C= \begin{pmatrix}-2 & -30 & 30\\ -6 & 6 & -3\\ -3 & 6 & 0\\ \end{pmatrix}\\ D= \begin{pmatrix}2 & -15 & 15\\ -6 & 3 & 3\\ 3 & 6 & 0\\ \end{pmatrix}\\$
得到的结果D不是 $1\times 1$ 的矩阵，那么计算B的行列式，但是把D作为第二步的运算结果：
$\begin{pmatrix}0 & 1 & -6 & 2\\ 2 & -4 & -6 & -3\\ -3 & 3 & 3 & 2\\ 3 & -2 & 0 & 0\\ \end{pmatrix}\\ B= \begin{pmatrix}2 & -15 & 15\\ -6 & 3 & 3\\ 3 & 6 & 0\\ \end{pmatrix}\\ C= \begin{pmatrix}-84 & -90\\ -45 & -18\\ \end{pmatrix}\\ D= \begin{pmatrix}21 & 15\\ -15 & -6\\ \end{pmatrix}$
$D$ 是一个 $2\times 2$ 的矩阵还不行，还得缩小，再来：
$\begin{pmatrix}2 & -15 & 15\\ -6 & 3 & 3\\ 3 & 6 & 0\\ \end{pmatrix}\\ B= \begin{pmatrix}21 & 15\\ -15 & -6\\ \end{pmatrix}\\ C= \begin{pmatrix}99\\ \end{pmatrix}\\ D= \begin{pmatrix}33\\ \end{pmatrix}$
最终结果出来了，是33.

Python实现

我只贴相关核心代码,我的代码中没有判断中心块是否含有0，不过这是为了学习算法，而不是用于实际生产项目，所以我就没继续完善：

    def interior(self):
        # Dodgson算法的中心块子矩阵
        n = len(self.__vectors)
        array = [[0 for _ in range(n - 2)] for _ in range(n - 2)]
        # i 代表行 j代表列
        for i in range(0, n - 2):
            for j in range(0, n - 2):
                # 如果column = 0
                array[j][i] = self.__vectors[j + 1][i + 1]
        return Matrix(array)

    def dodgson_condensation(self):
        n = len(self.__vectors)
        array = [[0 for _ in range(n - 1)] for _ in range(n - 1)]
        for i in range(0, n - 1):
            for j in range(0, n - 1):
                # 每个元素是2x2矩阵的行列式
                m = Matrix([[self.__vectors[i][j], self.__vectors[i][j + 1]],
                            [self.__vectors[i + 1][j], self.__vectors[i + 1][j + 1]]])
                array[i][j] = m.determinant2x2()
        return Matrix(array)

    def dodgson(self, b=None):
        n = len(self.__vectors)
        if n == 1:
            return self.__vectors[0][0]
        if n == 2 and b is None:
            return self.determinant2x2()

        s = self.interior()
        if b is None:
            b = self.dodgson_condensation()
        c = b.dodgson_condensation()
        d = c.divide_elements(s)

        if len(d.__vectors) == 1:
            return d.dodgson()
        return b.dodgson(d)

    def determinant2x2(self):
        return self.__vectors[0][0] * self.__vectors[1][1] - self.__vectors[0][1] * self.__vectors[1][0]

    def divide_elements(self, other):
        array = copy.deepcopy(self.__vectors)
        for i, vector in enumerate(array):
            for j, e in enumerate(vector):
                vector[j] = vector[j] / other.__vectors[i][j]
        return Matrix(array)