Wrong results

[saschatimme]

I did some test computations and found the following troubling result:

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.

[giordano]

I think this is going haywire because of the magnitude difference:

julia> using PolynomialRoots

julia> d = 20; c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1];

julia> PolynomialRoots.evalpoly(roots(c), c)
20-element Array{Complex{Float64},1}:
 9.081074789250304e16 + 1.2232077690076664e17im
 9.081074789250294e16 + 1.2232077690076717e17im
 9.081074789250203e16 + 1.2232077690076765e17im
 9.081074789250253e16 + 1.2232077690076682e17im
 9.081074789250288e16 + 1.2232077690076715e17im
 9.081074789250278e16 + 1.2232077690076686e17im
 9.081074789250277e16 + 1.223207769007676e17im 
 9.081074789250288e16 + 1.2232077690076754e17im
 9.081074789250253e16 + 1.2232077690076669e17im
 9.081074789250283e16 + 1.2232077690076715e17im
 9.081074789250288e16 + 1.2232077690076699e17im
 9.081074789250286e16 + 1.223207769007674e17im 
 9.081074789250362e16 + 1.2232077690076762e17im
 9.081074789250314e16 + 1.2232077690076688e17im
 9.081074789250256e16 + 1.2232077690076763e17im
 9.081074789250288e16 + 1.2232077690076726e17im
 9.081074789250272e16 + 1.2232077690076722e17im
 9.081074789250307e16 + 1.2232077690076648e17im
  9.08107478925015e16 + 1.2232077690076675e17im
 9.081074789250309e16 + 1.2232077690076734e17im

It's interesting that all the points are being evaluated more or less at the same, wildly wrong, value 🤔

When I don't scale the central elements, I get more reasonable results:

julia> d = 20; c = [randn(ComplexF64); randn(ComplexF64, d - 1); 1];

julia> PolynomialRoots.evalpoly(roots(c), c)
20-element Array{Complex{Float64},1}:
   8.326672684688674e-16 - 7.549516567451064e-15im 
   8.881784197001252e-15 + 4.9960036108132044e-15im
  3.3306690738754696e-15 + 7.105427357601002e-15im 
  1.8189894035458565e-11 + 0.0im                   
   9.103828801926284e-15 + 3.552713678800501e-15im 
   4.440892098500626e-16 + 3.6637359812630166e-15im
  1.0658141036401503e-14 + 0.0im                   
  1.4210854715202004e-14 + 0.0im                   
  1.3655743202889425e-14 - 5.495603971894525e-15im 
 -1.5987211554602254e-14 - 1.687538997430238e-14im 
  -9.947598300641403e-14 + 5.684341886080802e-14im 
    7.66053886991358e-15 - 1.1102230246251565e-15im
  1.3322676295501878e-15 + 1.2434497875801753e-14im
   7.105427357601002e-15 - 3.552713678800501e-15im 
  1.2212453270876722e-14 - 1.7763568394002505e-15im
                     0.0 - 9.992007221626409e-15im 
   5.329070518200751e-15 - 3.552713678800501e-15im 
  -8.881784197001252e-16 + 5.329070518200751e-15im 
  2.6645352591003757e-15 - 3.552713678800501e-15im 
   4.440892098500626e-16 + 1.942890293094024e-16im

I need to check whether the original Fortran implementation has the same troubles.

[saschatimme]

Scaling the central values is exactly the point of the example :)
The example is motivated by the results from this article: https://arxiv.org/abs/2001.05281

[giordano]

Scaling the central values is exactly the point of the example :)

Sure, I just wanted to make sure the root finder didn't get broken recently 🙂

[nsajko]

This has something to do with how the machine epsilon is used:

julia> using DynamicPolynomials, 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;

julia> c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1];

julia> for nbits in union(20:1:28, (64, 256))
         setprecision(nbits)
         local max_err = maximum(abs.(c - poly(ComplexF64.(roots(Complex{BigFloat}.(c), polish = false)))))
         println(nbits, "  ", max_err)
       end
20  1.519972640161657e-5
21  8.309140910263302e-6
22  3.878784367264271e-6
23  3.317700077628581e-6
24  1.4608485903698948e-6
25  8.50493687573117e-7
26  3.488836894278648e16
27  7.314329793362672e14
28  5.7399276787994944e17
64  4.896172571983926e15
256  4.896172572182062e15

julia> c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1];

julia> for nbits in union(20:1:28, (64, 256))
         setprecision(nbits)
         local max_err = maximum(abs.(c - poly(ComplexF64.(roots(Complex{BigFloat}.(c), polish = false)))))
         println(nbits, "  ", max_err)
       end
20  7.367320484632937e-6
21  7.998143075577152e-6
22  3.3275014752343082e-6
23  1.6213547856398292e-6
24  7.281698361264757e-7
25  3.37870495686925e12
26  9.920784674292804e11
27  1.391496433864446e13
28  6.00896607233545e12
64  4.534357255620861e12
256  4.534357256375943e12

julia> c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1];

julia> for nbits in union(20:1:28, (64, 256))
         setprecision(nbits)
         local max_err = maximum(abs.(c - poly(ComplexF64.(roots(Complex{BigFloat}.(c), polish = false)))))
         println(nbits, "  ", max_err)
       end
20  3.5433584824734573e-6
21  1.583409787159591e-6
22  1.356164689199754e-6
23  32341.518738578667
24  805614.4762540918
25  51793.74408068894
26  37838.74576482569
27  43181.31659319327
28  24163.354164083048
64  16709.59795418854
256  16709.597954278048

The idea was to vary the number of bits in BigFloat values (with setprecision) and see how much does it help, but I was in for a surprise!

Turns out that precision is actually relatively good at around 20 bits of precision, but somewhere between 20 and 30 bits (depending on the input value), the results suddenly become completely wrong. Also note that the maximum error stabilizes, at some point it stops depending on the BigFloat precision (but it's still completely wrong at that point).

I think the machine epsilon is basically the most important thing that changes between, e.g., 22 and 23 bits of precision, so the bug is probably related to eps and that stopping criterion thingie.

[nsajko]

Overriding epsilon with a "good" value when calling roots doesn't seem to help though.

[giordano]

I need to check whether the original Fortran implementation has the same troubles.

Following similar approach as #25 (comment):

import Downloads, PolynomialRoots

src = Downloads.download("http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/cmplx_roots_sg.f90")
libcmplxroots = joinpath(tempdir(), "libcmplxroots.so")
run(`gfortran -x f95 $(src) -shared -fPIC -o $(libcmplxroots)`)

function roots!(roots::Vector{ComplexF64}, poly::Vector{ComplexF64},
                degree::Integer, polish::Bool, use_roots::Bool)
    ccall((:cmplx_roots_gen_, libcmplxroots), Ptr{Cvoid},
          (Ptr{ComplexF64}, # roots
           Ptr{ComplexF64}, # poly
           Ptr{Cint}, # degree
           Ptr{Cint}, # polish_roots_after
           Ptr{Cint}),# use_roots_as_starting_points
          roots, poly, Ref{Cint}(degree),
          Ref{Cint}(polish), Ref{Cint}(use_roots))
    return roots
end

function roots(poly::Vector{<:Number}; polish::Bool=false)
    degree = length(poly) - 1
    roots!(Array{ComplexF64}(undef, degree), float(complex(poly)), degree, polish, false)
end

julia> d = 20; c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1];

julia> evalpoly.(roots(c), (c,)) # With upstream library
20-element Vector{ComplexF64}:
 6.795511325941872e34 - 1.5437795762995263e35im
 6.795511325941971e34 - 1.543779576299517e35im
 6.795511325941929e34 - 1.5437795762995083e35im
  6.79551132594187e34 - 1.5437795762995149e35im
 6.795511325941931e34 - 1.5437795762995116e35im
 6.795511325941875e34 - 1.5437795762995158e35im
 6.795511325941925e34 - 1.5437795762995188e35im
 6.795511325941927e34 - 1.5437795762995103e35im
 6.795511325941969e34 - 1.5437795762995164e35im
 6.795511325941953e34 - 1.5437795762995027e35im
  6.79551132594192e34 - 1.543779576299509e35im
 6.795511325941883e34 - 1.54377957629951e35im
 6.795511325941962e34 - 1.5437795762995177e35im
 6.795511325941907e34 - 1.543779576299515e35im
 6.795511325941883e34 - 1.5437795762995092e35im
 6.795511325941816e34 - 1.5437795762995114e35im
 6.795511325941931e34 - 1.5437795762995142e35im
 6.795511325941916e34 - 1.5437795762995162e35im
 6.795511325941931e34 - 1.5437795762995142e35im
 6.795511325942111e34 - 1.5437795762995188e35im

julia> evalpoly.(PolynomialRoots.roots(c), (c,)) # With PolynomialRoots.jl
20-element Vector{ComplexF64}:
  3.423287799691176e15 - 8.646704914995811e15im
  3.423287799691144e15 - 8.646704914995869e15im
  3.423287799691134e15 - 8.646704914995846e15im
  3.423287799691135e15 - 8.646704914995831e15im
  3.423287799691126e15 - 8.64670491499579e15im
  3.423287799691176e15 - 8.646704914995856e15im
 3.4232877996911025e15 - 8.646704914995833e15im
  3.423287799691136e15 - 8.646704914995835e15im
  3.423287799691091e15 - 8.646704914995843e15im
  3.423287799691131e15 - 8.64670491499586e15im
 3.4232877996911195e15 - 8.646704914995844e15im
  3.423287799691112e15 - 8.646704914995826e15im
  3.423287799691128e15 - 8.646704914995828e15im
  3.423287799691126e15 - 8.646704914995835e15im
  3.423287799691089e15 - 8.646704914995832e15im
  3.423287799691184e15 - 8.646704914995863e15im
  3.423287799691118e15 - 8.646704914995827e15im
 3.4232877996911175e15 - 8.64670491499583e15im
  3.423287799691123e15 - 8.646704914995827e15im
  3.423287799691152e15 - 8.646704914995844e15im

They both get it wildly wrong (but results are different, perhaps because of #25.

[giordano]

@nsajko perhaps looking at #25 might be easier since that one has somewhat less pathological values, and if we're lucky solving that issue might solve the issue here with BigFloat 🙂 🤞

If that's of any help, original Fortran source code, which is giving the correct result for #25, is at https://github.com/giordano/cmplx-roots-sg/blob/472b4d25a42e32a30a6817b937881de8516aa2f2/cmplx_roots_sg.f90, I tried to keep the Julia code as faithful to the original one as possible (but there are some limited changes), also comments should be the same, to make it easier to follow the code.

[saschatimme]

I just briefly looked at the Fortran source code but it looks like polynomials of degree larger than 4 get deflated (i.e. you find one zero, and divide the polynomial by (x - {{zero}})). But polynomial division is not numerically stable, so I don't think the algorithm is numerically stable in the current form. There is a classical paper by Wilkinson discussing some of the related problems https://doi.org/10.1093/imamat/8.1.16 .