CODEWITHSUNDEEP
pythonalgorithmmathlinear-algebra
pycoder
pycoder

Reputation: 547

Precise determinant of integer NxN matrix

Determinant definition has only additions, subtractions and multiplications. So a determinant of a matrix with integer elements must be integer.

However numpy.linalg.det() returns a "slightly off" floating-point number:

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

It gets worse for a larger matrix:

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

And it's wrong! I'm sure the correct answer is:

>>> 320 * 1024**7
377789318629571617095680

Of course, for a big matrix it may be a rather long integer. But python has long integers built in.

How can I get an exact integer value of the determinant instead of approximate floating point value?

Upvotes: 7

Views: 2161

Answers (2)

pycoder
pycoder

Reputation: 547

A simple practical way to calculate determinant of an integet matrix is 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]

This algorithm is reasonably fast (O(N³) complexity).

And it's an integer preserving algorithm. It does have a division. But as long as all the elements of M are integer, all intermediate calculations would be integer too (the division remainder will be zero).

As a bonus the same code works for fractions/floating-point/complex elements if you drop the assert line and replace the integer division // with a regular division /.

PS: Another alternative is to use sympy instead of numpy:

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

But somewhy that is MUCH slower than the above det() function.

# 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]

Summary:

  • det(M) is 5+ times slower than numpy.linalg.det(M),
  • det(M) is ~50 times faster than sympy.Matrix(M).det()

It becomes even faster without the assert line.

Upvotes: 15

kaya3
kaya3

Reputation: 51037

@pycoder's answer is the preferred solution; for comparison, I wrote a Gaussian elimination function using the Fraction class which allows exact arithmetic of rational numbers. It is about 11 times slower than the Bareiss algorithm on the same benchmark.

from fractions import Fraction

def det(matrix):
    matrix = [[Fraction(x, 1) for x in row] for row in matrix]
    n = len(matrix)
    d, sign = 1, 1
    for i in range(n):
        if matrix[i][i] == 0:
            j = next((j for j in range(i + 1, n) if matrix[j][i] != 0), None)
            if j is None:
                return 0
            matrix[i], matrix[j] = matrix[j], matrix[i]
            sign = -sign
        d *= matrix[i][i]
        for j in range(i + 1, n):
            factor = matrix[j][i] / matrix[i][i]
            for k in range(i + 1, n):
                matrix[j][k] -= factor * matrix[i][k]
    return int(d) * sign

Upvotes: 6

Related Questions