整数n×N矩阵的精确行列式

Precise determinant of integer NxN matrix(整数nN矩阵的精确行列式)
本文介绍了整数n×N矩阵的精确行列式的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

行列式定义只有加法、减法和乘法。因此,具有整数元素的矩阵的行列式必须是整数。

然而numpy.linalg.det()返回一个浮点数:

>>> import numpy
>>> M = [[-1 if i==j else 1 for j in range(7)] for i in range(7)]
>>> numpy.linalg.det(M)
319.99999999999994

矩阵越大,情况越糟:

>>> M = [[-1024 if i==j else 1024 for j in range(7)] for i in range(7)]
>>> numpy.linalg.det(M)
3.777893186295698e+23
>>> "%.0f" % numpy.linalg.det(M)
'377789318629569805156352'

这是错误的!我确信正确答案是:

>>> 320 * 1024**7
377789318629571617095680
当然,对于大型矩阵,它可能是一个相当长的整数。但PYTHON内置了长整数。

如何获取行列式的精确整数值而不是近似浮点值?

推荐答案

计算整数矩阵行列式的一种简单实用的方法是Bareiss algorithm。

def det(M):
    M = [row[:] for row in M] # make a copy to keep original M unmodified
    N, sign, prev = len(M), 1, 1
    for i in range(N-1):
        if M[i][i] == 0: # swap with another row having nonzero i's elem
            swapto = next( (j for j in range(i+1,N) if M[j][i] != 0), None )
            if swapto is None:
                return 0 # all M[*][i] are zero => zero determinant
            M[i], M[swapto], sign = M[swapto], M[i], -sign
        for j in range(i+1,N):
            for k in range(i+1,N):
                assert ( M[j][k] * M[i][i] - M[j][i] * M[i][k] ) % prev == 0
                M[j][k] = ( M[j][k] * M[i][i] - M[j][i] * M[i][k] ) // prev
        prev = M[i][i]
    return sign * M[-1][-1]

此算法相当快(O(N³)复杂性)。

它是一个整数保持算法。它确实有一个部门。但只要M的所有元素都是整数,所有中间计算也都是整数(除法余数将为零)。

额外的好处是,如果您删除assert行,并将整数除法//替换为规则除法/,则相同的代码适用于分数/浮点/复数元素。


ps:另一种替代方法是使用sympy而不是NumPy:

>>> import sympy
>>> sympy.Matrix([ [-1024 if i==j else 1024 for j in range(7)] for i in range(7) ]).det()
377789318629571617095680

但这比上面的det()函数慢得多。

# Performance test: `numpy.linalg.det(M)` vs `det(M)` vs `sympy.Matrix(M).det()`
import timeit
def det(M):
    ...
M = [[-1024 if i==j else 1024 for j in range(7)] for i in range(7)]
print(timeit.repeat("numpy.linalg.det(M)", setup="import numpy; from __main__ import M", number=100, repeat=5))
#: [0.0035009384155273, 0.0033931732177734, 0.0033941268920898, 0.0033800601959229, 0.0033988952636719]
print(timeit.repeat("det(M)", setup="from __main__ import det, M", number=100, repeat=5))
#: [0.0171120166778564, 0.0171020030975342, 0.0171608924865723, 0.0170948505401611, 0.0171010494232178]
print(timeit.repeat("sympy.Matrix(M).det()", setup="import sympy; from __main__ import M", number=100, repeat=5))
#: [0.9561479091644287, 0.9564781188964844, 0.9539868831634521, 0.9536828994750977, 0.9546608924865723]

摘要:

  • det(M)numpy.linalg.det(M)慢5倍以上,
  • det(M)sympy.Matrix(M).det()快50倍

没有assert行会更快。

这篇关于整数n×N矩阵的精确行列式的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持编程学习网!

本站部分内容来源互联网,如果有图片或者内容侵犯您的权益请联系我们删除!

相关文档推荐

Leetcode 234: Palindrome LinkedList(Leetcode 234:回文链接列表)
How do I read an Excel file directly from Dropbox#39;s API using pandas.read_excel()?(如何使用PANDAS.READ_EXCEL()直接从Dropbox的API读取Excel文件?)
subprocess.Popen tries to write to nonexistent pipe(子进程。打开尝试写入不存在的管道)
I want to realize Popen-code from Windows to Linux:(我想实现从Windows到Linux的POpen-code:)
Reading stdout from a subprocess in real time(实时读取子进程中的标准输出)
How to call type safely on a random file in Python?(如何在Python中安全地调用随机文件上的类型?)