Possible bug in determinant calculation

Specific Domains Numerics
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

https://github.com/JuliaLang/julia/issues/30917

6 Likes
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.

1 Like
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.

1 Like
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.

2 Likes
9

Is there any update on how/when this will be fixed? The bug is still present in the latest release (1.1.1). For example:

julia> using LinearAlgebra

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

julia> det(A)
17.999999999999996

julia> det(Float64.(A))
5.999999999999998
10

It is not marked release critical for 1.2, so it is not likely to be fixed in the next release.

There is an open PR for it:

https://github.com/JuliaLang/julia/pull/30953

11

Thanks for the rapid feedback.

From my perspective, the fact that det potentially returns a dramatically incorrect result for integer matrices is an embarrassment for Julia. This is hardly an obscure, rarely used function. I have no thoughts on the best way to make this work properly, but do suggest that something (anything!) be done soon, and allow better solutions in future versions.

12

It’s a bug. There are many other bugs.

There is no need to dramatize the issue; it is unlikely to make it get fixed any sooner.

You already have a workaround, as suggested above, eg

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

julia> det(Rational{BigInt}.(A))
6//1
2 Likes
13

To be fair, some bugs are more important than others.

1 Like
14

Just checked and this is not yet fixed in Julia 1.2.

15

This has just been fixed on master. I’ve marked it for backporting to 1.3 but there is a risk that it will only make it into 1.4.

3 Likes
16

This is good news. Thank you.

I would like to encourage that the return type of det match the element type of the matrix. So if the entries in a matrix are of type Int, then the return value should also be Int. This is consistent with real-valued matrices returning a real determinant, and complex returning complex. Likewise for Rational.

I’m aware there’s an overflow danger, but we have that with all integer arithmetic operations.

julia> A
2Ă—2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
julia> det(A)
-2.0

julia> det(A .+ 0im)
-2.0 - 0.0im

julia> A = [1//1 2//1; 3//1 4//1]
2Ă—2 Array{Rational{Int64},2}:
 1//1  2//1
 3//1  4//1

julia> det(A)
-2//1
17

Why do need this?

I think it is better to leave the exact type unspecified, as long as the call is type stable. This gives wiggle room for future changes. You can always convert explicitly if preferred.

18

I wouldn’t mind that but it requires a new implementation for integer matrices since the current version is based on the LU factorization.

19

I don’t feel too strongly about this issue, but let me at least add a
motivating example:

A=(rand(Int64,2,2))
2x2 Array{Int64,2}:
-7148224804043827887 1550119287025645287
4310643582572096085 2785513071409872787

A[1,1]*A[2,2]-A[1,2]*A[2,1]
5953626061256297424

det(A)
-2.659348538587869e37

A user might reasonably expect these to return the same value.

What I find perhaps more disappointing is that for a BigInt array it
still gives a floating point number and not a BigInt. BTW - this can
be done efficiently. You pick a sequence of primes p_1,p_2,…,p_k
such that their product is larger than the determinant (easy using
Hadamard’s inequality) and compute the determinant in GF(p_i) for each
i and then use the Chinese remainder theorem to compute the remainder
modulo p_1p_2…*p_k which uniquely specifies the determinant.

Mathematically, the determinant is just a polynomial in the entries of
the matrix and thus if you have a matrix whose elements are in a ring
R, the determinant should be an element of R as well. Integers in
Julia are in point of fact elements of a ring (e.g. Z/2^64Z) and as
such it does make sense that the determinant returns the corresponding
ring element. It would seem crazy to a user if they constructed a
matrix of elements of integers modulo 3^40, took its determinant, and
got a floating point number. In a way, Julia currently does the same.

Tamas, I think your question is a good one. In point of fact, the
determinant is a very important mathematical quantity, but not the
most important for numerical analysis. It is certainly not a good way
to test for invertibility. If M is the 500 by 500 identity matrix then
Julia reports:

det(M/5)
0.0

rank(M)
500

What is most important of course is that if Julia chooses to do a
floating point computation that it gives the correct number! I’m
pleased to see this is getting fixed.

Danny

20

My understanding is that the determinant has a very limited role to play in numerical computing (if any) because of its numerical properties. (Some people go further, eg Sheldon “Down with Determinants” Axler even wrote a neat textbook in line with this idea — but that is a digression here). In almost all scenarios I am aware of, one actually wants something else for practical purposes, eg parts of the SVD, the QR, etc.

LinearAlgebra could include all kinds of methods for special cases of the determinant, but since they use algorithms that are distinct from the usual floating point ones, this imposes a maintenance burden.

Before this is undertaken, it would be great to understand

  1. what the uses cases are,
  2. what are the important trade-offs (eg speed versus accuracy, especially if exact results are to be expected),
  3. how these special algorithms scale with matrix size.

Perhaps a dedicated library with an integer_determinant or something like that would be the right place to experiment first.