Determinant of an integer matrix is float

Why isn’t it an int?

julia> a, c, b, d = M = [10 20; 30 40]; a*d - b*c , det(M)
(-200, -200.00000000000003)
1 Like

Nowadays, Julia includes an integer-determinant algorithm (the Bareiss algorithm), but only uses it by default for BigInt matrices (since for fixed-width integer types like the default Int a determinant will quickly overflow unless the matrix is tiny).

julia> using LinearAlgebra

julia> M = BigInt[10 20; 30 40]
2×2 Matrix{BigInt}:
 10  20
 30  40

julia> det(M) # exact BigInt determinant
-200

You can invoke the Bareiss algorithm explicitly for other integer types by calling LinearAlgebra.det_bareiss, but beware that it can easily overflow for larger matrices:

julia> LinearAlgebra.det_bareiss([10 20; 30 40]) # exact Int determinant
-200

For example, with a 100 \times 100 matrix of random ±1 entries:

julia> A = rand([-1,1], 100,100);

julia> det(A) # floating-point approximation: fast and reasonably accurate
7.014043363270618e77

julia> det(BigInt.(A)) # exact
701404336326996708202480358186962679200983525458711992094310103304386652405760

julia> LinearAlgebra.det_bareiss(A) # overflows — bogus answer
746

PS. I just pushed a PR to clarify the det documentation: document exact BigInt determinants by stevengj · Pull Request #53579 · JuliaLang/julia · GitHub

12 Likes

Using the wikipedia link in stevengj’s answer, the following function might be useful:

bareiss_bits(M) = (n = size(M,1); 
  ndigits(maximum(abs,M); base=2)*n+ceil(Int,log2(n)*n/2))

It calculates a bound on the bits needed to calculate det_bareiss and thus, if it is less than 64 or 128, Int64 and Int128 can be used, respectively (not exact, but hopefullly a conservative evaluation and thus safe).

1 Like

The problem/concern I think with that is that it makes the return type of det depend on not just the type but also the values of the arguments.