# Accuracy of Computing Determinants

**URL:** <https://discourse.julialang.org/t/accuracy-of-computing-determinants/6782>\
**Category:** Numerics\
**Created:** [October 30, 2017, 6:04pm UTC](https://discourse.julialang.org/t/accuracy-of-computing-determinants/6782 "2017-10-30T18:04:14Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Hans\_W\_Borchers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hans_w_borchers/32/8481_2.png) [@Hans\_W\_Borchers](https://discourse.julialang.org/u/Hans_W_Borchers)\
**Post date:** [October 30, 2017, 6:04pm UTC](https://discourse.julialang.org/t/accuracy-of-computing-determinants/6782/1 "2017-10-30T18:04:14Z")

</div>

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

```julia
    | 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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 30, 2017, 6:14pm UTC](https://discourse.julialang.org/t/accuracy-of-computing-determinants/6782/2 "2017-10-30T18:14:22Z")

</div>

> [@Hans\_W\_Borchers](#):
>
> Though the first four systems use Float64, the differences seem to be significant. Are there systematic differences in the way determinants are computed?

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

```julia
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:

> <https://github.com/SciML/DelayDiffEq.jl/pull/22>
>
> opening this so it's easier to track the printing since AppVeyor doesn't run non…-PRs by default.

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

```

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [October 30, 2017, 7:35pm UTC](https://discourse.julialang.org/t/accuracy-of-computing-determinants/6782/3 "2017-10-30T19:35:26Z")

</div>

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

```julia-auto
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-auto
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-auto
julia> det(lufact(Float64.(A)))
2.550554736789249e-78

```

#### Bunch-Kaufman using upper and lower triangle for storage

```julia-auto
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-auto
julia> det(cholfact(Symmetric(Float64.(A), :L)))
2.6834441543359152e-78

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

```

---

<div class="post-metadata">

**Author:** ![certik](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/certik/32/2651_2.png) [@certik](https://discourse.julialang.org/u/certik)\
**Post date:** [November 19, 2017, 1:13am UTC](https://discourse.julialang.org/t/accuracy-of-computing-determinants/6782/4 "2017-11-19T01:13:54Z")

</div>

The discussion here is very interesting. In [https://github.com/JuliaLang/julia/issues/24568](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](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.