roots has issue finding the root of a polynomial with degree 14

[ahmetumutdurmus]

Hi,

I am trying to write a general function to compute the Internal Rate of Return of a bond. This is basically a root finding problem and I was trying to use PolynomialRoots.jl as a root finder. The roots function works fine up to a degree of 13 yet is not able to find the "correct" roots when the degree of the polynomial is increased to 14. For example 1.03 is a root of both polynomials p1 and p2 where:

julia> p1 = [103, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, -100]
14-element Array{Int64,1}:
  103
    3
    3
    3
    3
    3
    3
    3
    3
    3
    3
    3
    3
 -100

julia> p2 = [103, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, -100]
15-element Array{Int64,1}:
  103
    3
    3
    3
    3
    3
    3
    3
    3
    3
    3
    3
    3
    3
 -100

Yet when roots is executed for both polynomials the following outputs are obtained:

julia> roots(p1)
13-element Array{Complex{Float64},1}:
      1.03-2.13163e-16im
 -0.354605-0.935016im
 -0.748511+0.663123im
  0.120537+0.992709im
  0.885456-0.464723im
 -0.748511-0.663123im
  0.885456+0.464723im
 -0.354605+0.935016im
 -0.970942-0.239316im
  0.120537-0.992709im
 -0.970942+0.239316im
  0.568065-0.822984im
  0.568065+0.822984im

julia> roots(p2)
14-element Array{Complex{Float64},1}:
  0.697266-0.419327im
 0.0780983+0.692255im
 -0.705236+0.0995877im
 -0.255796-0.709392im
  0.891654-0.00461915im
 -0.532087-0.516335im
 -0.292549+0.555004im
  0.684213+0.400383im
  0.421208-0.663104im
 -0.692196-0.226128im
 0.0834591-0.763406im
 -0.568103+0.389205im
  0.405608+0.623887im
 -0.185538+0.54199im

Notice that 1.03 is included in the first array yet not the second array. Yet 1.03 is a known root of both polynomials. So rootsmisses this root for the second polynomial.

Thanks in advance.

[giordano]

Thank you so much for the report. Indeed, the second polynomial evaluates at about 102 for each of the found roots for p2, while the first one correctly evaluates to about 0:

julia> PolynomialRoots.evalpoly(roots(p1), p1)
13-element Array{Complex{Float64},1}:
  -9.094947017729282e-13 + 3.32913226551714e-13im  
  -1.564224092715191e-13 + 4.1190877941585115e-13im
 -6.9784461093520134e-15 + 5.654123852771096e-14im 
   3.490439179678138e-16 + 1.1991155346654323e-13im
   -7.69978360342577e-14 + 2.641645392280729e-13im 
 -1.4634125542698607e-14 + 1.884707950923698e-13im 
 -1.5838622493146062e-13 + 6.604113480701821e-14im 
  1.5813312362319058e-13 - 7.972427988693891e-14im 
   6.561477091702649e-14 + 2.7887217117352817e-13im
 -1.5589116955448176e-13 - 3.526810396074801e-15im 
   6.561477091702649e-14 - 2.7887217117352827e-13im
 -1.9345164540469515e-15 + 2.3390608302340083e-14im
 -1.9345164540469515e-15 - 2.3390608302340083e-14im

julia> PolynomialRoots.evalpoly(roots(p2), p2)
14-element Array{Complex{Float64},1}:
 102.08947146280127 + 0.9151682784614544im
  102.0894714628013 + 0.9151682784614517im
  102.0894714628013 + 0.9151682784614563im
  102.0894714628013 + 0.9151682784614509im
  102.0894714628013 + 0.915168278461436im 
  102.0894714628013 + 0.9151682784614503im
  102.0894714628013 + 0.9151682784614514im
  102.0894714628013 + 0.9151682784614511im
 102.08947146280128 + 0.9151682784614481im
  102.0894714628013 + 0.9151682784614529im
 102.08947146280131 + 0.9151682784614431im
  102.0894714628013 + 0.9151682784614512im
  102.0894714628013 + 0.9151682784614517im
  102.0894714628013 + 0.915168278461451im

You can get much better results by polishing the result:

julia> roots(p2, polish = true)
14-element Array{Complex{Float64},1}:
   0.9009688679024191 - 0.4338837391175581im    
   0.6234898018587336 + 0.7818314824680298im    
  -0.9009688679024191 - 0.4338837391175581im    
  -0.2225209339563144 - 0.9749279121818236im    
                 1.03 - 4.0389678347315804e-28im
  -0.6234898018587336 - 0.7818314824680298im    
  -0.9009688679024191 - 0.4338837391175581im    
   0.9009688679024191 + 0.4338837391175581im    
   0.6234898018587335 - 0.7818314824680298im    
  -0.6234898018587335 - 0.7818314824680298im    
   0.2225209339563144 - 0.9749279121818236im    
                 -1.0 + 0.0im                   
   0.6234898018587336 + 0.7818314824680298im    
 -0.18553787475062206 - 0.5419896481974021im    

julia> PolynomialRoots.evalpoly(roots(p2, polish = true), p2)
14-element Array{Complex{Float64},1}:
  1.9977807247279814e-14 - 3.699515267932936e-14im 
  3.3381016150334446e-15 + 1.1110493609124174e-13im
  1.7163173186238412e-14 + 4.9326870239105815e-14im
  -1.581106331772122e-15 - 6.9272794589055566e-15im
 -1.1368683772161603e-12 + 6.901111370345701e-25im 
  1.4039165064188364e-14 + 8.888394887299339e-14im 
  1.7163173186238412e-14 + 4.9326870239105815e-14im
  1.9977807247279814e-14 + 3.699515267932936e-14im 
   9.947598300641403e-14 + 0.0im                   
   9.947598300641403e-14 + 0.0im                   
 -1.5791961046974125e-14 + 6.9272794589055566e-15im
                     0.0 + 0.0im                   
  3.3381016150334446e-15 + 1.1110493609124174e-13im
      102.08947146279803 - 0.9151682784526727im

It gets correctly 13 out of 14 roots, including 1.03. However, the last one is still much wrong :confused:

[giordano]

If I'm not mistaken, the root with real part -0.22252.... should have multiplicity 2, instead of 1. Thus, the last one (-0.18553.... + ....) should be the same as the fourth one. Probably, errors accumulate during division for polishing the roots and the last root deviates too much from the true value.

Unfortunately, high-order polynomials are very difficult to deal with.

Note that also having a very rough idea of where the roots are would help very much. For example:

julia> r = roots(p2, polish = true)
14-element Array{Complex{Float64},1}:
   0.9009688679024191 - 0.4338837391175581im    
   0.6234898018587336 + 0.7818314824680298im    
  -0.9009688679024191 - 0.4338837391175581im    
  -0.2225209339563144 - 0.9749279121818236im    
                 1.03 - 4.0389678347315804e-28im
  -0.6234898018587336 - 0.7818314824680298im    
  -0.9009688679024191 - 0.4338837391175581im    
   0.9009688679024191 + 0.4338837391175581im    
   0.6234898018587335 - 0.7818314824680298im    
  -0.6234898018587335 - 0.7818314824680298im    
   0.2225209339563144 - 0.9749279121818236im    
                 -1.0 + 0.0im                   
   0.6234898018587336 + 0.7818314824680298im    
 -0.18553787475062206 - 0.5419896481974021im

julia> r[end] = r[4]
-0.2225209339563144 - 0.9749279121818236im

julia> r2 = round.(r) # This is a very rough estimation of the real roots
14-element Array{Complex{Float64},1}:
  1.0 - 0.0im
  1.0 + 1.0im
 -1.0 - 0.0im
 -0.0 - 1.0im
  1.0 - 0.0im
 -1.0 - 1.0im
 -1.0 - 0.0im
  1.0 + 0.0im
  1.0 - 1.0im
 -1.0 - 1.0im
  0.0 - 1.0im
 -1.0 + 0.0im
  1.0 + 1.0im
 -0.0 - 1.0im

julia> roots(p2, r2, polish = true)
14-element Array{Complex{Float64},1}:
  0.6234898018587336 + 0.7818314824680298im   
 -0.2225209339563144 + 0.9749279121818236im   
 -0.6234898018587335 + 0.7818314824680298im   
  0.9009688679024191 - 0.4338837391175581im   
  0.9009688679024191 + 0.4338837391175581im   
 -0.9009688679024191 - 0.4338837391175581im   
 -0.9009688679024191 + 0.4338837391175581im   
                1.03 + 3.944304526105059e-31im
  0.6234898018587335 - 0.7818314824680298im   
 -0.6234898018587336 - 0.7818314824680298im   
  0.2225209339563144 - 0.9749279121818236im   
                -1.0 + 0.0im                  
  0.2225209339563144 + 0.9749279121818236im   
 -0.2225209339563144 - 0.9749279121818236im   

julia> PolynomialRoots.evalpoly(roots(p2, r2, polish = true), p2)
14-element Array{Complex{Float64},1}:
  3.3381016150334446e-15 + 1.1110493609124174e-13im
  -1.581106331772122e-15 + 6.9272794589055566e-15im
   9.947598300641403e-14 + 0.0im                   
  1.9977807247279814e-14 - 3.699515267932936e-14im 
  1.9977807247279814e-14 + 3.699515267932936e-14im 
  1.7163173186238412e-14 + 4.9326870239105815e-14im
  1.7163173186238412e-14 - 4.9326870239105815e-14im
 -1.1368683772161603e-12 - 6.739366572603224e-28im 
   9.947598300641403e-14 + 0.0im                   
  1.4039165064188364e-14 + 8.888394887299339e-14im 
 -1.5791961046974125e-14 + 6.9272794589055566e-15im
                     0.0 + 0.0im                   
 -1.5791961046974125e-14 - 6.9272794589055566e-15im
  -1.581106331772122e-15 - 6.9272794589055566e-15im

[ahmetumutdurmus]

Thanks a lot for your response and help. I guess I should close the issue right?

[giordano]

To be honest, I'm not going to look into the algorithm I've just implemented :slightly_smiling_face:

The README.md already mentions that high-order polynomials are subject to this issues, I don't think there is much that can be done. However, I think I've shown you that with some tricks you can get much better results and my recommendation is to always (when possible) check that the number you get is actually a root of the given polynomial

[ahmetumutdurmus]

Alright thanks again 🙂