计算矩阵行列式的函数

Question

我想让我的函数计算输入矩阵 A 的行列式，使用行归约将 A 转换为阶梯形式，之后行列式应该只是 A 的对角线的乘积。我可以假设 A 是 n x n np.array

这是我已有的代码：

def determinant(A):
    A = np.matrix.copy(A)
    row_switches = 0
    #  Reduce A to echelon form
    for col in range(A.shape[1]):
        if find_non_zero(A, col) != col:
            #  Switch rows
            A[[find_non_zero(A, col), col], :] = A[[col, find_non_zero(A, col)], :]
            row_switches += 1
        #  Make all 0's below "pivot"
        for row in range(col+1, A.shape[0]):
            factor = A[row, col] / A[col, col]
            A[row, :] = A[row, :] - factor * A[col, :]

    return A.diagonal().prod() * (-1) ** row_switches


# Find first non-zero value starting from diagonal element
def find_non_zero(A, n):
    row = n
    while row < A.shape[0] and A[row, n] == 0:
        row += 1
    return row

然后我将我的结果与 np.linalg.det(A) 进行比较。对于低于 50x50（2.8e-08 差异）的浮点数随机矩阵，差异是可控的，但在 70x70 之后，差异平均在 1000 到 10'000 之间。这可能是什么原因？

我的代码的另一个问题是对于整数矩阵 A = np.random.randint(low=-1000,high=1000,size=(25, 25))，差异更加疯狂：

1820098560（我的）与 1.0853429659737294e+81（numpy）

Answer 1

整数数组有两个问题，您可以通过将函数的第一行更改为 A = np.matrix(A, dtype=float).

来解决它们

您可能会溢出并完全放弃您的结果。

>>> np.arange(1, 10).prod() # correct
362880
>>> np.arange(1, 20).prod() # incorrect
109641728
>>> np.arange(1, 20, dtype=float).prod() # correct
1.21645100408832e+17

无论第 A[row, :] = A[row, :] - factor * A[col, :] 行中 rhs 的结果是什么，它都将被转换为整数。

>>> a = np.zeros((3,), dtype=int)
>>> a[0] = 2.4
>>> a
array([2, 0, 0])

至于浮点数组的不准确性，您不得不接受它们，因为浮点运算的精度有限。当对角线的乘积给你一个像 6.59842495617676e+17 这样的数字而 numpy 给你 6.598424956176783e+17 时，你可以看到结果非常接近。但是它们只能表示这么多数字，当数字非常大时，最后几位的差异实际上意味着 1000s 的差异。这只会变得更糟，你的矩阵越大，结果你的数字就越大。但就相对差异而言，即 (your_method - numpy) / numpy，无论您使用的数字大小如何，它都相当不错。

算法的稳定性

关于你的 factor 值的一点，当它比 wikipedia 非常小时：

One possible problem is numerical instability, caused by the possibility of dividing by very small numbers. If, for example, the leading coefficient of one of the rows is very close to zero, then to row-reduce the matrix, one would need to divide by that number. This means that any error existed for the number that was close to zero would be amplified. Gaussian elimination is numerically stable for diagonally dominant or positive-definite matrices. For general matrices, Gaussian elimination is usually considered to be stable, when using partial pivoting, even though there are examples of stable matrices for which it is unstable.[11]

[snip]

This algorithm differs slightly from the one discussed earlier, by choosing a pivot with largest absolute value. Such a partial pivoting may be required if, at the pivot place, the entry of the matrix is zero. In any case, choosing the largest possible absolute value of the pivot improves the numerical stability of the algorithm, when floating point is used for representing numbers.

如果重要的话，numpy 使用 LAPACK's LU decomposition algorithm 它实现了 Sivan Toledo 的递归 LU 算法的迭代版本。

计算矩阵行列式的函数

Function for calculating the determinant of a matrix

python

numpy

determinants

算法的稳定性