# Polynomial Roots With Multiplicities

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?

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

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

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.

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.

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

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

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.

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.