# In numerical computation, how could increasing the number of bits of precision decrease accuracy?

**URL:** <https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600>\
**Category:** Offtopic\
**Created:** [December 13, 2022, 12:56pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600 "2022-12-13T12:56:22Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [December 13, 2022, 12:56pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/1 "2022-12-13T12:56:22Z")

</div>

In this PolynomialRoots.jl Github Issue, accuracy of the result is relatively OK when the MPFR/`BigFloat` precision is at around (just!) twenty bits, but goes sour suddenly when increasing the number of bits and stays that way when increasing the precision afterwards:

> <https://github.com/giordano/PolynomialRoots.jl/issues/19>
>
> I did some test computations and found the following troubling result:
> 
> \`\`\`jul…ia
> julia\> using DynamicPolynomials
> 
> julia\> using PolynomialRoots
> 
> julia\> function poly(roots)
> @polyvar x
> d = length(roots)
> f = prod(x - xi for xi in roots)
> \[coefficient(f, x^i) for i in 0:d\]
> end
> poly (generic function with 1 method)
> 
> julia\> d = 20
> 20
> 
> julia\> c = \[randn(ComplexF64); randn(ComplexF64, d - 1) \* 1e-4; 1\]
> 21-element Array{Complex{Float64},1}:
> 0.6038117420721143 - 0.4692317923460019im    
> 8.677470319492543e-5 - 1.9856488102339103e-5im 
> -2.8145901548439397e-5 + 2.5958816494607353e-6im 
> -9.933810484391282e-5 - 0.0001480730021139855im 
> -8.820204766296887e-5 - 0.0001042484315957653im 
> 7.88021087360911e-5 + 2.5513881910383407e-5im 
> -5.007828702552179e-5 - 4.850423221617593e-7im  
> 3.788521654994467e-5 + 2.0944404532776053e-5im 
> -1.1315308780483735e-5 - 1.8001817822660527e-5im 
> 0.00011608934805185687 + 0.0001503409434982439im 
> 4.125906444101907e-5 + 9.978779625362144e-5im  
> -5.52502458307844e-5 + 1.244622053838063e-5im  
> 3.170867062110031e-5 - 1.9111588344365587e-5im 
> 9.754810243956861e-5 + 6.369172830914271e-5im  
> 7.028073862086189e-5 - 3.3819410763212794e-5im 
> 0.0001238375611115403 - 5.144632426358618e-5im  
> -5.263257681962449e-5 + 6.227546399484231e-5im  
> 4.6728130048484104e-5 - 0.00011791878030739642im
> 0.00010423858969876596 - 0.00011299961775497805im
> -8.557623396259921e-5 - 2.825091927864975e-5im  
> 1.0 + 0.0im                   
> 
> julia\> r = roots(c)
> 20-element Array{Complex{Float64},1}:
> 7.802412795926557 - 4.178757163336667im 
> -1.2054492306931206 + 8.768496484683855im 
> -8.768489889160787 - 1.2054534486770816im
> -4.178754803019869 - 7.80240884258056im  
> 8.768498428644229 + 1.2054564365850402im
> 1.2054578038795405 - 8.768493818879023im 
> 6.385317333146652 + 6.129227118201205im 
> -6.129219730795856 + 6.38531432406573im  
> 3.856072372001663 - 7.966827038663316im 
> -7.9668235583329405 - 3.856068311866256im 
> -7.802404396639403 + 4.1787600200142485im
> 4.178763518989671 + 7.802411631716743im 
> 6.129228142519897 - 6.385311565395753im 
> -6.385308625464333 - 6.129224228027822im 
> 8.71184335999314 - 1.5631567624600453im
> 1.563163615363657 + 8.711841873146625im 
> -3.856063894157322 + 7.966829725325496im 
> 7.966832198295311 + 3.856071278929593im 
> -1.5631549515439616 - 8.711839169108076im 
> -8.711834912718764 + 1.563159707245344im 
> 
> julia\> c - poly(r)
> 21-element Array{Complex{Float64},1}:
> -7.987123734845452e18 + 3.4639346941464863e18im 
> 896.0000867747032 + 1151.999980143512im     
> -16.00002814590155 + 18.00000259588165im     
> -8.000099338104844 - 14.500148073002114im    
> 0.24991179795233703 - 0.18760424843159576im   
> 0.1563288021087361 + 0.17190051388191038im   
> -0.007862578287025522 - 4.850423221617593e-7im  
> -0.00020625540845005532 + 0.000753366279532776im  
> -1.1315308780483735e-5 - 0.00020110728657266052im
> 2.4536613676856873e-5 - 2.513513072050611e-5im  
> 2.6352546265659437e-6 + 1.5593416637776903e-6im 
> -2.6542446384986325e-8 + 1.0805906987477109e-7im 
> -1.7764205750584404e-8 + 6.601467340955743e-9im  
> -9.523320931490147e-10 - 4.682197821015313e-9im  
> 5.2134349463363316e-11 - 3.523188618357654e-10im 
> 5.818546518080239e-12 + 2.4179138238041714e-13im
> 5.263254309923089e-12 + 7.162041142441677e-13im 
> 1.88754143337546e-13 - 2.443320881988925e-13im 
> -4.330127294054076e-14 - 1.3839515471125718e-14im
> -2.160977560089483e-15 - 1.5989373182683647e-15im
> 0.0 + 0.0im          
> \`\`\`
> 
> I am not sure whether this is a problem with the algorithm or the specific implementation.

How could this possibly be the case?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 13, 2022, 1:38pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/2 "2022-12-13T13:38:28Z")

</div>

> [@nsajko](#):
>
> How could this possibly be the case?

Usually it means that you’re not computing what you think you are, and some part of the computation is limited by the precision of your original data.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [December 13, 2022, 1:42pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/3 "2022-12-13T13:42:42Z")

</div>

But we’re not just talking about accuracy failing to improve as precision is increased. The accuracy actually deteriorates dramatically!  
For example, in the first example in my first comment on the issue, at precision of 25 bits, the maximum error is `8.5e-7`, but thereafter the error is huge (for 26 bits of precision it’s `3.5e16`, note that the exponent is _positive_).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 13, 2022, 2:45pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/4 "2022-12-13T14:45:45Z")

</div>

Looking at your code, it seems that you are:

1. generating random `Complex64` polynomial coefficients `c`
2. computing roots `r` by `c -> BigFloat -> roots -> Complex64`
3. computing new polynomial coefficients `c2` by multiplying monomials `(x-r[i])` from the roots `r`
4. comparing `c - c2`

To quote my comment above, _you’re not computing what you think you are_.

`c` and `c2` **should not match** , so I don’t understand the point of this “test”. (As soon as you round the roots `r` to `ComplexF64`, you dramatically change the coefficients of the corresponding polynomial, and then you introduce additional errors from roundoff when you multiply the monomials together.)

(Realize also that the roots of a polynomial are extremely sensitive to the coefficients — the sensitivity [increases exponentially with the degree](https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial) — so it’s not surprising to me that you are seeing a huge mismatch here.)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 13, 2022, 4:04pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/5 "2022-12-13T16:04:38Z")

</div>

A more appropriate test would seem to be:

1. generating random `Complex64` polynomial coefficients `c`
2. convert to `BigFloat` (an exact conversion) and find roots `r = roots(big.(c))`
3. construct new `BigFloat` polynomial `c2` coefficients by multiplying monomials
4. compare to `c`.

This _should_ match if you have sufficiently high precision. However, I find that it gives a wildly wrong error with higher and higher probability as the degree increases:

```julia
using PolynomialRoots, DynamicPolynomials

function poly(roots)
   @polyvar x
   d = length(roots)
   f = prod(x - xi for xi in roots)
   [coefficient(f, x^i) for i in 0:d]
end

d = 5
c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1]
setprecision(8192)
r = roots(big.(c), polish=true)
c2 = ComplexF64.(poly(r))
c2 - c

```

90% of the time gives (correctly)

```julia
6-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

```

but if we increase the degree to `d = 6` it is wildy wrong about 85% of the time (≈ 15% of the time it correctly gives zero), giving answers for `c - c2` like

```julia
7-element Vector{ComplexF64}:
 3.9439141771536586e6 - 7.66389580196865e6im
 2.2082351722049564e6 - 2.857913579123853e6im
    469737.1655469684 - 420674.6692220386im
    50125.96985319289 - 30578.538953102627im
    2872.684011890461 - 1098.2968892416864im
     84.4854525989318 - 15.599756742987397im
                  0.0 + 0.0im

```

and if you check with `evalpoly` you’ll find that these aren’t roots at all, and indeed p(r) seems to have a constant(!!) offset:

```julia
julia> ComplexF64.(evalpoly.(r, Ref(c)))
6-element Vector{ComplexF64}:
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im

```

In fact, if you look at the “roots” they are all the same, which is clearly wrong:

```julia
julia> ComplexF64.(r)
6-element Vector{ComplexF64}:
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im

```

Conclusion: PolynomialRoots appears to be returning garbage roots some of the time. (This also falls under the category of “you aren’t computing what you think you are”.)

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [December 13, 2022, 6:19pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/6 "2022-12-13T18:19:34Z")

</div>

> [@stevengj](#):
>
> As soon as you round the roots `r` to `ComplexF64`, you dramatically change the coefficients of the corresponding polynomial, and then you introduce additional errors from roundoff when you multiply the monomials together.

If I understand correctly, this is the main mistake in my experiment? It seems easy to fix though, to prevent rounding I just remove the cast to `ComplexF64`, and to prevent additional errors when multiplying the monomials together I just temporarily increase the `BigFloat` precision?

When I repeat the experiment with those corrections, the results still exhibit the same peculiar pattern.

Using this code:

```julia
using Printf, DynamicPolynomials, PolynomialRoots

function poly_coefs_from_roots(roots::AbstractVector{<:Complex{T}}) where {T <: AbstractFloat}
  local old_prec = precision(T)
  setprecision(T, 64*old_prec)

  @polyvar x
  local f = prod(x - xi for xi in roots)
  local ret = [coefficient(f, x^i) for i in 0:length(roots)]

  setprecision(T, old_prec)

  ret
end

random_poly_coefs(d) =
  Complex{BigFloat}[randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1]

max_coef_error(coefs) =
  maximum(abs.(coefs - poly_coefs_from_roots(roots(coefs, polish = true))))

approx_sprintf(n::AbstractFloat) =
  @sprintf("%.2e", n)

function experiment(d, mantissa_lengths)
  local coefs = random_poly_coefs(d)

  for nbits in mantissa_lengths
    setprecision(nbits)
    local max_err = max_coef_error(coefs)
    println(nbits, " ", approx_sprintf(max_err))
  end

  nothing
end

```

These are some example results:

```none
julia> experiment(20, 19:31)
19 8.39e-06
20 3.49e-06
21 1.89e-06
22 1.08e-06
23 4.16e-07
24 2.70e-07
25 1.45e-07
26 5.80e+07
27 6.55e+16
28 2.05e+18
29 6.16e+18
30 6.24e+18
31 6.80e+18

julia> experiment(20, 19:31)
19 4.01e-06
20 2.17e-06
21 1.10e-06
22 4.27e-07
23 1.86e-07
24 3.86e+07
25 4.82e-08
26 9.67e+10
27 1.33e+11
28 1.80e+11
29 1.42e+11
30 9.14e+10
31 9.59e+10

julia> experiment(20, 19:31)
19 5.51e-06
20 2.16e-06
21 1.27e-06
22 7.19e-07
23 3.68e-07
24 2.19e-07
25 1.24e-07
26 2.17e+13
27 6.89e+15
28 1.67e+16
29 7.63e+16
30 7.93e+16
31 1.37e+17

julia> experiment(20, 19:31)
19 4.10e-06
20 2.76e-06
21 1.11e-06
22 5.48e-07
23 2.46e-07
24 5.76e+07
25 1.45e+12
26 2.35e+13
27 4.89e+13
28 2.00e+14
29 4.19e+14
30 8.40e+14
31 1.28e+15

```

There’s a new variation on the established pattern in the second example, but the main pattern is the same as before.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [December 13, 2022, 7:23pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/7 "2022-12-13T19:23:54Z")

</div>

This seems to be the most significant part:

> [@stevengj](#):
>
> Conclusion: PolynomialRoots appears to be returning garbage roots some of the time.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [December 13, 2022, 9:08pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/8 "2022-12-13T21:08:00Z")

</div>

Yeah, but my point with this whole thing is that the accuracy is for some reason actually good when one does `setprecision(20)` during `roots`.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [December 13, 2022, 10:08pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/9 "2022-12-13T22:08:11Z")

</div>

I’ve looked at the code in PolynomialRoots and some functions are very long and convoluted. The sort of which I would hesitate to use without formal proof (at least in space probes).

Luckily, at the cost of run-time, this function can be made robust, by doing `evalpoly` internally on found roots and filtering out glaring problems or even throwing an exception.

This is the usual trade-off between allow quick path and double down on safety, and the best solution IMHO is to have both paths and default on the safe path for ‘new’ users.

Perhaps, if a Julia flag is set, these function can log an INFO or WARN level message to remind of quicker versions.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 13, 2022, 10:24pm UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/10 "2022-12-13T22:24:59Z")

</div>

> [@nsajko](#):
>
> Yeah, but my point with this whole thing is that the accuracy is for some reason actually good when one does `setprecision(20)` during `roots`.

Probably just luck? `setprecision(20)` doesn’t work for my `d=6` example above, but in any case my example shows that there is some randomness in whether the problem occurs.

Clearly, there’s a bug here. It’s usually a distraction, in my experience, to try to perform too much mathematical analysis of a bug. Better to just fix the bug.

I submitted an issue with a simple test case, in which the incorrect results are obtained for `ComplexF64` as well (`BigFloat` is just a distraction here): [sometimes returns bogus duplicate roots · Issue #25 · giordano/PolynomialRoots.jl · GitHub](https://github.com/giordano/PolynomialRoots.jl/issues/25)

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [December 14, 2022, 1:04am UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/11 "2022-12-14T01:04:21Z")

</div>

> [@Dan](#):
>
> I’ve looked at the code in PolynomialRoots and some functions are very long and convoluted. The sort of which I would hesitate to use without formal proof (at least in space probes).

Good that you haven’t looked at the [original Fortran code](http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/) because that’s 50% longer 😂

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [December 14, 2022, 1:08am UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/12 "2022-12-14T01:08:29Z")

</div>

> [@giordano](#):
>
> Good that you haven’t looked at the [original Fortran code](http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/)

But I trust the FORTRAN code… it usually comes from the time of Real Programmers

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [December 14, 2022, 1:14am UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/13 "2022-12-14T01:14:36Z")

</div>

Ah yes, because Julia’s intense testing of old numerical libraries written in Fortran has never turned up any bugs 🐛

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [December 14, 2022, 1:20am UTC](https://discourse.julialang.org/t/in-numerical-computation-how-could-increasing-the-number-of-bits-of-precision-decrease-accuracy/91600/14 "2022-12-14T01:20:32Z")

</div>

In the bogus cases found by Steven above it may be an error in the Julia implementation because the Fortran library [gives the correct result](https://github.com/giordano/PolynomialRoots.jl/issues/25#issuecomment-1350145189) 🥲 (yes, the code is a mess, but the starting point wasn’t much better and a mistake in the translation is almost unavoidable)

But @nsajko started this discussion because was looking into [another issue](https://github.com/giordano/PolynomialRoots.jl/issues/19) where the ability of using arbitrary precision may be beneficial, except the solver goes haywire at some point. It’d be great if that was solved by finding the culprit in the other issue 🙂