Polynomial Roots With Multiplicities

Specific Domains Numerics
1

When using the roots function from the Polynomials package, the values found for multiple roots can be quite inaccurate. Just as an example, here is the polynomial with zeros at 1 and -1, each of multiplicity 4.

julia> using Pynomials
julia> p = Poly([1, 0, -4, 0, 6, 0, -4, 0, 1])
Poly(1 - 4*x^2 + 6*x^4 - 4*x^6 + x^8)
julia> Polynomials.roots(p)
8-element Array{Complex{Float64},1}:
  -1.000063312800363 + 6.33123618252134e-5im 
  -1.000063312800363 - 6.33123618252134e-5im 
 -0.9999366871996364 + 6.33132485906425e-5im 
 -0.9999366871996364 - 6.33132485906425e-5im 
  1.0000782706569455 + 0.0im                 
  0.9999999979609331 + 7.826862289848257e-5im
  0.9999999979609331 - 7.826862289848257e-5im
  0.9999217334211853 + 0.0im

I understand that this inaccuracy is introduced by the eigenvalue function that finds all these roots in one go. Each of these roots can be refined by applying, e.g., Newton’s or Laguerre’s method.

I thought, package PolynomialRoots would provide such a feature but, at least in the way I use it, returns a similarly inaccurate solution.

julia> using PolynomialRoots
julia> p = complex([1.0, 0, -4, 0, 6, 0, -4, 0, 1]);
julia> PolynomialRoots.roots(p)
8-element Array{Complex{Float64},1}:
  -1.0000703710890786 + 6.890751979805258e-9im 
   0.9999999994387893 + 8.914437556930486e-5im 
  -1.0000000085541791 - 7.037275037868812e-5im 
   0.9999999926924392 - 8.91443743617057e-5im  
   1.0000891483097958 - 3.3737785193170825e-9im
   -0.999929625585618 - 6.892300583852365e-9im 
   0.9999108595589759 + 3.3725709200610236e-9im
  -0.9999999947711241 + 7.037275192729228e-5im

Is there a way – or a Julia function – to “polish” these roots all in once and without utilizing multiple-precision numbers?

2

PolynomialRoots.jl has a “polish” keyword argument, see the docstring. I can’t try it now to see if it does the trick.

Edit: just tried, it doesn’t help much. Multiple precision instead would give much better results:

julia> Complex{Float64}.(PolynomialRoots.roots(big.(p)))
8-element Array{Complex{Float64},1}:
                 1.0 - 3.2604121516820557e-21im
                -1.0 + 3.2718015945471428e-21im
                 1.0 + 6.556134361696995e-17im 
                -1.0 - 6.555274705703271e-17im 
                 1.0 - 6.556134361696995e-17im 
 -0.9999999999999999 - 3.2718015945471435e-21im
  0.9999999999999999 + 3.2604121516820564e-21im
                -1.0 + 6.555274705703271e-17im

but I know that you would like to avoid this

3

If you know you have multiplicities and you have real coefficients, the PolynomialZeros package has a multroot function that can work for this example:

using Polynomials, PolynomialsZeros
julia> p = Poly([1, 0, -4, 0, 6, 0, -4, 0, 1.0])
Poly(1.0 - 4.0*x^2 + 6.0*x^4 - 4.0*x^6 + 1.0*x^8)

julia> PolynomialZeros.MultRoot.multroot(p)
[ Info: The multiplicity count may be in error--the initial guess for the roots failed to improve when refined along the pejorative manifold.
([-1.0, 1.0], [4, 4])
1 Like
4

Real coefficients only would be okay.
I don’t know in advance whether I have multiple roots.
Guess you mean PolynomialZeros, not PolynomialsZeros.

Unfortunately, I get an error when using it the way you described:

julia> using Polynomials, PolynomialZeros
julia> p = Poly([1, 0, -4, 0, 6, 0, -4, 0, 1.0])
Poly(1.0 - 4.0*x^2 + 6.0*x^4 - 4.0*x^6 + 1.0*x^8)

julia> PolynomialZeros.MultRoot.multroot(p)
ERROR: UndefVarError: diagm not defined
Stacktrace:
 [1] #pejroot#5(::Nothing, ::Float64, ::Int64, ::Function, 
 ::Poly{Float64}, ::Array{Float64,1}, ::Array{Int64,1}) at 
 [...].julia/packages/PolynomialZeros/vjnjt/src/agcd/multroot.jl:81
 [2] pejroot(::Poly{Float64}, ::Array{Float64,1}, ::Array{Int64,1}) at 
 [...].julia/packages/PolynomialZeros/vjnjt/src/agcd/multroot.jl:76
 [3] #multroot#10(::Float64, ::Float64, ::Float64, ::Float64, ::Function, 
 ::Poly{Float64}) at 
 [...].julia/packages/PolynomialZeros/vjnjt/src/agcd/multroot.jl:213
 [4] multroot(::Poly{Float64}) at /home/hwb/.julia/packages
 /PolynomialZeros/vjnjt/src/agcd/multroot.jl:157
 [5] top-level scope at none:0

Why should diagm not be defined – do I need to load another package?
This is Julia 1.1.0 on Manjaro Linux.
Thanks for your help.

5

Sorry, my guess is you need to use master. I just pushed many changes I had locally. The multroot algorithm doesn’t perform well with high degree polynomials with no multiplicities, though should be okay for degree 8-15.

6

I beueve that it is generally not possible to calculate multiplicities without exact arithmetic (eg. Factoring polynomials with rational coefficients).

7

@dpsanders I am not so much interested in exact multiplicities, but in more accurate root finding. There is a lot of research going on in this area, see for example the articles mentioned in the PolynomialZeros package. In Matlab there is the poly_root function that gives better result than roots for multiple roots.

In another post here on Julia Discourse I asked for Julia implementations of “Rank-revealing QR/SVD Decompositions”. It appears these play an important role in determining the multiplicites of polynomial roots.

8

For refining all roots at once, there’s the Aberth method and the Durand-Kerner method. I don’t know of an implementation of either in Julia, though I think I coded one of them once.

The Aberth article I linked states (with a reference) that these both converge linearly at a multiple root.

9

I know the Aberth method (1967), the Jenkins-Traub method (1972) that is implemented in R in the polyroot function, or the (Weierstrass-)Durand-Kerner method (1966) that I once implemented myself. In my experiences, none of them is really much better than computing the eigenvalues of the companion matrix.

There are more modern approaches, such as Zeng’s MULTROOT toolbox (exploited in the PolynomialZeros package), or “rank-revealing decomposition” methods. I would just like to hear of experiences or implementations of such methods in Julia (or other scientific computing software).

@j_verzani – Thanks for pointing out your package. I have not yet downloaded the ‘master’ version, still the references were very helpful, giving some insight into recent advances in polynomial root solving.