# Polynomial Roots With Multiplicities

**URL:** <https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309>\
**Category:** Numerics\
**Created:** [May 17, 2019, 8:25am UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309 "2019-05-17T08:25:16Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Hans\_W\_Borchers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hans_w_borchers/32/8481_2.png) [@Hans\_W\_Borchers](https://discourse.julialang.org/u/Hans_W_Borchers)\
**Post date:** [May 17, 2019, 8:25am UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/1 "2019-05-17T08:25:17Z")

</div>

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?

---

<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:** [May 17, 2019, 9:03am UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/2 "2019-05-17T09:03:48Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [May 17, 2019, 10:26am UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/3 "2019-05-17T10:26:20Z")

</div>

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

```julia
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])

```

---

<div class="post-metadata">

**Author:** ![Hans\_W\_Borchers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hans_w_borchers/32/8481_2.png) [@Hans\_W\_Borchers](https://discourse.julialang.org/u/Hans_W_Borchers)\
**Post date:** [May 17, 2019, 11:54am UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/4 "2019-05-17T11:54:02Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [May 17, 2019, 4:43pm UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/5 "2019-05-17T16:43:24Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [May 17, 2019, 7:08pm UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/6 "2019-05-17T19:08:19Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Hans\_W\_Borchers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hans_w_borchers/32/8481_2.png) [@Hans\_W\_Borchers](https://discourse.julialang.org/u/Hans_W_Borchers)\
**Post date:** [May 18, 2019, 4:32pm UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/7 "2019-05-18T16:32:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [May 18, 2019, 4:47pm UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/8 "2019-05-18T16:47:10Z")

</div>

For refining all roots at once, there’s the [Aberth method](https://en.wikipedia.org/wiki/Aberth_method) and the [Durand-Kerner method](https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_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.

---

<div class="post-metadata">

**Author:** ![Hans\_W\_Borchers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hans_w_borchers/32/8481_2.png) [@Hans\_W\_Borchers](https://discourse.julialang.org/u/Hans_W_Borchers)\
**Post date:** [May 18, 2019, 9:12pm UTC](https://discourse.julialang.org/t/polynomial-roots-with-multiplicities/24309/9 "2019-05-18T21:12:54Z")

</div>

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.