Accuracy of Computing Determinants

Specific Domains Numerics
1

While trying out the R packages JuliaCall and XRJulia for calling Julia functions from R, I got interested in quite accurately determining the determinant of the 12-by-12 Hilbert matrix.

    julia> A = [1/(i+j-1) for i in 1:12, j in 1:12];
    julia> det(A)
    2.732614032613825e-78

Let’s compare this with some other systems I can put my hands on.

    | System      | Value        | Comment  |
    |-------------|--------------|----------|
    | R           | 2.550555e-78 | Float64  |
    | Julia       | 2.732614e-78 |    "     |
    | MATLAB      | 2.858986e-78 |    "     |
    | Octave      | 2.683444e-78 |    "     |
    | Maxima      | 2.637781e-78 | Rational |
    | Pari-GP     | 2.637781e-78 |    "     |
    | Maple       | 2.637781e-78 |    "     |
    | Mathematica | 2.637781e-78 |    "     |

Is this good enough? Pari-GP gets it exact even if we define the matrix in floating point, probably because it calculates in Float128 by default (I know we could do this in Julia, too). Though the first four systems use Float64, the differences seem to be significant. Are there systematic differences in the way determinants are computed?

2

Floating point can be even weirder than that. At that level of accuracy, it’s architecture-dependent. For example:

julia> A = [1/(i+j-1) for i in 1:12, j in 1:12]
12×12 Array{Float64,2}:
 1.0        0.5        0.333333   …  0.1        0.0909091  0.0833333
 0.5        0.333333   0.25          0.0909091  0.0833333  0.0769231
 0.333333   0.25       0.2           0.0833333  0.0769231  0.0714286
 0.25       0.2        0.166667      0.0769231  0.0714286  0.0666667
 0.2        0.166667   0.142857      0.0714286  0.0666667  0.0625
 0.166667   0.142857   0.125      …  0.0666667  0.0625     0.0588235
 0.142857   0.125      0.111111      0.0625     0.0588235  0.0555556
 0.125      0.111111   0.1           0.0588235  0.0555556  0.0526316
 0.111111   0.1        0.0909091     0.0555556  0.0526316  0.05
 0.1        0.0909091  0.0833333     0.0526316  0.05       0.047619
 0.0909091  0.0833333  0.0769231  …  0.05       0.047619   0.0454545
 0.0833333  0.0769231  0.0714286     0.047619   0.0454545  0.0434783

julia> det(A)
2.550554736789249e-78

with just a Julia generic binary on Windows 10. The compiler can optimize differently on different architectures which can change the result at this level!!! I heard about weird stuff like this, but first encountered it when doing some precision debugging of DelayDiffEq.jl where different OS’s would give slightly different answers:

The moral of the story is, 64-bit floating points are not trustworthy at this level, no matter what you use. Different BLAS/LAPACK backends can do different ordering of the summations (and the compiler can sometimes change it, and SIMD can affect this, etc), and some will be better/worse on a given problem + architecture combo, but it can change depending on the problem you’re looking to solve. The only way to really make this accurate is embrace the fact that Julia gives you generic algorithms over abstract types and compute it using something that is efficient yet accurate enough for your problem. For example, Julia can do it with rationals as well:

julia> A = [Rational{BigInt}(1/(i+j-1)) for i in 1//1:12//1, j in 1//1:12//1]
12×12 Array{Rational{BigInt},2}:
 1//1   1//2   1//3   1//4   1//5   …  1//8   1//9   1//10  1//11  1//12
 1//2   1//3   1//4   1//5   1//6      1//9   1//10  1//11  1//12  1//13
 1//3   1//4   1//5   1//6   1//7      1//10  1//11  1//12  1//13  1//14
 1//4   1//5   1//6   1//7   1//8      1//11  1//12  1//13  1//14  1//15
 1//5   1//6   1//7   1//8   1//9      1//12  1//13  1//14  1//15  1//16
 1//6   1//7   1//8   1//9   1//10  …  1//13  1//14  1//15  1//16  1//17
 1//7   1//8   1//9   1//10  1//11     1//14  1//15  1//16  1//17  1//18
 1//8   1//9   1//10  1//11  1//12     1//15  1//16  1//17  1//18  1//19
 1//9   1//10  1//11  1//12  1//13     1//16  1//17  1//18  1//19  1//20
 1//10  1//11  1//12  1//13  1//14     1//17  1//18  1//19  1//20  1//21
 1//11  1//12  1//13  1//14  1//15  …  1//18  1//19  1//20  1//21  1//22
 1//12  1//13  1//14  1//15  1//16     1//19  1//20  1//21  1//22  1//23

julia> float(det(A))
2.637780651253547321325265140355620571956761490130113311468588530215014561601783e-78
8 Likes
3

Notice that there are two different matrices considered here. There is the actual Hilbert matrix

A = Rational{BigInt}[1//(i+j-1) for i in 1:12, j in 1:12]

and the floating point representation of this matrix Float64.(A). These two matrices have slightly different determinants even when computed exactly

julia> Float64(det(A))
2.6377806512535473e-78

julia> Float64(det(Rational{BigInt}.(Float64.(A))))
2.687225581661903e-78

Furthermore, since the matrix is so ill conditioned, different ways of computing the determinant in floating point arithmetic give quite deifferent results.

LU

julia> det(lufact(Float64.(A)))
2.550554736789249e-78

Bunch-Kaufman using upper and lower triangle for storage

julia> det(Symmetric(Float64.(A), :L))
2.8930243994501074e-78

julia> det(Symmetric(Float64.(A), :U))
2.557691999626182e-78

Cholesky using uppler and lower triangle for storage

julia> det(cholfact(Symmetric(Float64.(A), :L)))
2.6834441543359152e-78

julia> det(cholfact(Symmetric(Float64.(A), :U)))
2.8372335715728185e-78
7 Likes
4

The discussion here is very interesting. In https://github.com/JuliaLang/julia/issues/24568 we had a long discussions on the problems/merits of -ffast-math. My understanding from it is that if you don’t use -ffast-math, or other value changing optimizations, and Julia on purpose does not do that, then you get predictable results from Julia. That is the motivation to keep the pi sum benchmark without -ffast-math.

However, from the discussion here it seems that Julia’s determinant is not predictable and depends on the platform. Isn’t that inconsistent with the conclusion that we have reached in https://github.com/JuliaLang/julia/issues/24568?

Because if it cannot be predicted what exact floating point number you get for a well defined problem (even if ill-conditioned), then we might as well give up on predictability like that and enable -ffast-math for the pi sum benchmark I would think.