Possible bug in determinant calculation


#1

I found a (not poorly conditioned) integer matrix that julia computes the determinant incorrectly for.

using LinearAlgebra
A=[1 1 0 0 1 0 0 0; 1 0 1 0 0 1 0 0; 1 0 0 1 0 0 1 0; 0 1 1 1 0 0 0 0; 0 1 0 0 0 0 1 1; 0 0 1 0 1 0 0 1; 0 0 0 1 1 1 0 0; 0 0 0 0 1 1 0 1];

julia> det(A)
96.0

julia> det(float.(A))
6.0

The second value is correct. The LU-decomposition when it is treated as an Int64 array has one enormous and one tiny value. This does not happen when I first convert it to float.

The example, is hardly unique:

n=8; for it=1:10^6;  R=rand(0:1,n,n); if abs(det(R)-det(float.(R)))>0.01; println(R); end; end

should generate a few examples in a matter of seconds.


#2

I think the generic path (which is called for Matrix{Int}) does the calculation without pivoting. What you are seeing is simply ill-conditioned matrix LU w/o pivoting. Perhaps the pivoting path could be called when LinearAlgebra.lutype(T) <: AbstractFloat, which holds for Int.


#3

Filed


#4

I’d like to make a request that det returns a number that has the same type as the entries in the matrix. So if the matrix is populated by Ints then the determinant should also be an Int.


#5

It’s a little tricky. If the elements are Int64 and big do you want the determinant modulo 2^64? Because in general the answer can be bigger and then it can’t be the same type.


#6

Yeah, I think I’d still want an Int64 result even if it’s “wrong” in the same way that multiplication or addition could result in wrong answers. It’d be consistent with the rest of Julia:

julia> 2^62
4611686018427387904

julia> 2^63
-9223372036854775808

#7

Note that you can get the right answer in this case by using rationals:

julia> A=[1 1 0 0 1 0 0 0; 1 0 1 0 0 1 0 0; 1 0 0 1 0 0 1 0; 0 1 1 1 0 0 0 0; 0 1 0 0 0 0 1 1; 0 0 1 0 1 0 0 1; 0 0 0 1 1 1 0 0; 0 0 0 0 1 1 0 1];

julia> det(A//1)
6//1

Though in the general case you would likely see an overflow at some point.


#8

You can do this without overflowing by using Rational{BigInt}, although this can get expensive quickly. I don’t think it should use BigInt unless you explicitly request this.

There are specialized algorithms for integer determinants (e.g. Dodgeson condensation), but in general doing computational number theory efficiently is probably outside the scope of the standard LinearAlgebra library and is more something for packages like Nemo.jl.