探讨使用 Python 实现计算 n 阶行列式的多种方法

前些日子刚学线性代数，被手算高阶行列式折磨（主要因为当时才只学了定义），于是想用Python来做验算。自然也就想到了Python强大的第三方科学计算库 —— Numpy。经过一通操作，便有了以下的动态：

图片

以下是输入和输出的结果：

1
2
3
4
input:[[ 2  1  4  1],
       [ 3 -1  2  1],
       [ 1  2  3  2],
       [ 5  0  6  2]]
1
output: 6.217248937900884e-15

但是这个行列式的结果应该是0才对啊！！！

观察 Numpy 计算的结果，猜测可能是浮点数导致的误差。那既然现成的库不靠谱，就只能自己造轮子了。当时我用定义法写了一个简单的小程序。后来随着更深入的学习，我也了解到了更多的计算方法，也翻了一些参考资料，所以今天就来讨论一下计算行列式的多种方法。

我简单总结了一下计算行列式的几种方法，具体如下：

图片

按照以上总结的方法，下面就来一一实现。此处只讨论实现方式，不考虑时间复杂度和优化问题。

一、逆序数法

逆序数法总共分为三步：计算逆序数、计算全排列、通过定义计算结果。

我们将这三步分为三个函数来一步步实现。

首先是计算逆序数，最容易想到的就是这个算法：

1
2
3
4
5
6
7
def get_t(rank):
    t = 0
    for i1 in range(len(rank)):
        for i2 in range(i1 + 1, len(rank)):
            if rank[i1] > rank[i2]:
                t += 1
    return t

简单解释一下就是从一串数组中挨个遍历每一个数，然后查找该数前面大于该数的数量并且计数。

第二个函数是获得 n 个数的全排列，可以使用递归的方法来解决：确定第1位，对n-1位进行全排列，确定第二位，对n-2位进行全排列……

1
2
3
4
5
6
7
8
9
10
11
12
def get_permutations(rank):
    length = len(rank)
    permutations = []
    def _permute(index=0):
        if index == length:
            permutations.append(rank[0:length])
        for i in range(index, length):
            rank[i], rank[index] = rank[index], rank[i]
            _permute(index + 1)
            rank[i], rank[index] = rank[index], rank[i]
    _permute()
    return permutations

后来发现还有更简单直接的方法：调用Python的内置库 itertools 来实现：

1
2
3
import itertools
def get_permutations(rank):
    return list(itertools.permutations(rank))

准备好以上两个函数后就可以根据定义进行计算了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def calculate_det_1(det):
    # 行列式为0时直接返回0
    if det == 0:
        return 0
    # 行列式为一阶时直接返回结果
    if len(det) == 1:    
        return det[0]

    n = len(det)
    p = get_permutations(range(1, n + 1))
    result = 0
    # 遍历全排列
    for r in p:    
        t = get_t(r)
        i = 0
        part = (-1) ** t
        # 遍历排列中元素并计算结果
        for c in r:   
            part *= det[i][c - 1]
            i += 1
        result += part
    return result

二、等价转换法

该方法就是通过不断的变换将行列式转换为上三角形行列式，然后计算对角线的乘积。具体的实现方法如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

def calculate_det_2(det):
    n = len(det)
    result = 1
    # 遍历列
    for col in range(n):
        row = col
        result *= det[row][col]
        # 寻找不是0的行
        while det[row][col] == 0 and row < n - 1:
            row += 1
        # 化简行列式
        for i in range(row + 1, n):
            if det[i][col] == 0:
                pass
            else:
                k = - det[i][col] / det[row][col]
                for j in range(col, n):
                    det[i][j] += det[row][j] * k
    return result

三、代数余子式法

该方法就是将行列式不断降阶，直到变为容易计算的低阶行列式。

该方法比较适合用递归来解决，实现代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

def calculate_det_3(det):
    n = len(det)
    # 处理基础情况
    if n == 1:
        return det[0][0]
    if n == 2:
        return det[0][0]*det[1][1] - det[0][1]*det[1][0]
    # 求代数余子式并计算结果
    result = 0
    for i in range(n):
        C = []
        for j in range(n):
            C.append(det[j][1:])
        C.remove(C[i])
        if i % 2 == 1:
            result -= calculate_det_3(C) * det[i][0]
        else:
            result += calculate_det_3(C) * det[i][0]
    return result

以上就是个人对于使用 Python 计算行列式的三种常用方法的总结，个人能力有限，可能存在错误，欢迎各位指正。