Using Optimization to Extract Roots of Real Coefficient Polynomials

[Namir]

Given the polynomail P(x) with real coefficients and order n. This polynomials can be factorized by extracting a quadrtic polynomial q(x). This step yields the reduced polynomial B(x) and the remainder mononial r(x).

Code:

P(x) = q(x)B(x) + r(x)

Where:

Code:

P(x) = x^n + a(1)*x^(n-1) + ... + a(n-1)*x + a(n)
q(x) = x^2 + c(1)*x + c(2)
r(x) = d(1)*x + d(2)
B(x) = x^(n-2) + b(1)*x^(n-3) + ... + b(n-3)*x + b(n-2)

When q(x) is an exact factor of P(x), r(x)=0, that is d(1) and d(2) are zero.

The Lin-Bairstow algorithm factors out q(x) repeatedly to obtain real and imaginary roots calculated from the quadratic polynomial q(x).

I devised another algorithm that extracts the quadratic polynomial q(x) by optimization:

Given tthe coefficients a(1) to a(n), the tolerance toler, and the optimized function Fx(c(1),c(2)) that returns the value of sqrt(d(1)^2 + d(2)^2). Hereis a generic version of the pseudo-code for the Optimization-Based Quadratic Reducrion (OBQR):

1. Initialize variables.
2. While order>=2.
2.1 Initialize c(1) and c(2).
2.2 Minimize function Fx by using an algorithm that changes the values of c(1) and c(2) to yield a value for sqrt(d(1)^2 + d(2)^2) that is close to zero. Use the tolerance value and possibly a maximum iteration limit to stop the optimization process.
2.3 Extract x^2 + c(1)*x + c(2) from P(x) and obtain B(x).
2.4 Calculates the roots of q(x), such that roots(i,1)=real_part and roots(i,2) = imaginary_part (or 0 for real roots).
2.5 Set the coefficients of P(x) to be those of B(x).
3. If order = 1, calculate the last root as -a(1).
4. Return the roots.

The optimizing function divides P(x) by q(x) and obtains the coefficients for the simple polynomial r(x). The function then returns sqrt(d(1)^2 + d(2)^2).

You have a wide variety of choices for the optimizing algorithms. You can choose legacy deterministic algorithms or the newer evolutionary algorithms that rely on random numbers. The latter algorithms yield slightly different results each time you use them.

The efficiency and accuracy of the quadratic polynomial extraction depends on the algorithm you choose, the implementation (whether it is yours, one you find in some toolbox library, or part of the programming language you are using) and the code efficiency of the implementation.


Here is a Matlab version that uses the built-in fminsearch() function (which implements the Nelder-Mead Simplex method) which works very well.

Code:

function [xroots,nroots] = slb_fminsearch(polyCoeffs, x)
% Obtains the roots of a polynomial using obtaimizatio0n-based quadratic
% reduction.
  global polcoeff
  
  if nargin<3, x=[0 0]; end

  rng('shuffle')
  polcoeff = polyCoeffs/polyCoeffs(1);
  order=length(polcoeff)-1;
  xroots=zeros(order,2);
  nroots=0;
  while order>=2   
    x = fminsearch(@optimFun, x);
    xaug = [1 x]; % augment x
    [r1,i1,r2,i2]=quadratic(xaug);
    xroots(nroots+1,1)=r1;
    xroots(nroots+1,2)=i1;
    xroots(nroots+2,1)=r2;
    xroots(nroots+2,2)=i2; 
    nroots = nroots + 2;
    order = order - 2;
    [q,~] = deconv(polcoeff, xaug);
    polcoeff = q;      
  end
  
  if order == 1
    xroots(nroots+1,1)=-polcoeff(2)/polcoeff(1);
    xroots(nroots+1,2)=0;
    nroots = nroots + 1;
  end
end

function [r1,i1,r2,i2]=quadratic(coeff)
  discr=coeff(2)^2-4*coeff(1)*coeff(3);
  if discr>0
      i1=0;
      i2=0;
      r1=(-coeff(2)+sqrt(discr))/2/coeff(1);
      r2=(-coeff(2)-sqrt(discr))/2/coeff(1);
  else
      r1 =-coeff(2)/2/coeff(1);
      r2 = r1;
      i1 = sqrt(abs(discr))/2/coeff(1);
      i2 = -i1;
  end
end

function y = optimFun(x)
  global polcoeff
  xaug = [1 x];
  [~,r] = deconv(polcoeff,xaug);
  y = norm(r);
end

Here i an example in Matlab. I am using the Matlab function roots() to obtain accurate roots for comparison.

Code:

>> C1=[1 -2 3 -4 5 -6 7 -8 9 -10];
>> roots(C1)

ans =

   1.3386 + 0.0000i
   1.0970 + 0.7570i
   1.0970 - 0.7570i
   0.4657 + 1.2295i
   0.4657 - 1.2295i
  -0.3103 + 1.2423i
  -0.3103 - 1.2423i
  -0.9217 + 0.7964i
  -0.9217 - 0.7964i

>> [xroots,nroots] = slb_fminsearch(C1,[0 0])

xroots =

    0.4657    1.2295
    0.4657   -1.2295
   -0.9217    0.7964
   -0.9217   -0.7964
   -0.3103    1.2423
   -0.3103   -1.2423
    1.0970    0.7570
    1.0970   -0.7570
    1.3386         0


nroots =

     9

Here is a Matlab function that uses the evolutionary particle swarm optimization function (that is part of Matlab as of Matlab 2018a):

Code:

function [xroots,nroots] = slb(polyCoeffs,toler,maxIter,maxRetry)
%SLB Summary of this function goes here
%   Detailed explanation goes here
  global polcoeff
  
  if nargin<2, toler=1e-10; end
  if nargin<3, maxIter=10000; end
  if nargin<4, maxRetry=1; end

  rng('shuffle')
  polcoeff = polyCoeffs/polyCoeffs(1);
  order=length(polcoeff)-1;
  xroots=zeros(order,2);
  nroots=0;
  while order>=2   
    lb = [-20 -20];
    ub = [20 20];
    nvars = 2;
    options = optimoptions('particleswarm', 'MaxIterations', maxIter);
    options = optimoptions('particleswarm', 'FunctionTolerance', toler);
    numRetry = maxRetry;
    while numRetry>0
      [x,~,exitflag] = particleswarm(@optimFun,nvars,lb,ub,options);
      if exitflag==1, break; end
      numRetry = numRetry - 1;
    end
    xaug = [1 x]; % augment x
    [r1,i1,r2,i2]=quadratic(xaug);
    xroots(nroots+1,1)=r1;
    xroots(nroots+1,2)=i1;
    xroots(nroots+2,1)=r2;
    xroots(nroots+2,2)=i2; 
    nroots = nroots + 2;
    order = order - 2;
    [q,~] = deconv(polcoeff, xaug);
    polcoeff = q;      
    if exitflag~=1
      fprintf('Exit flag is %i ', exitflag)
      fprintf('for roots (%f,i%f) and (%f,i%f)\n', r1, i1, r2, i2);
    end 
  end

  if order == 1
    xroots(nroots+1,1)=-polcoeff(2)/polcoeff(1);
    xroots(nroots+1,2)=0;
    nroots = nroots + 1;
  end
end

function [r1,i1,r2,i2]=quadratic(coeff)
  discr=coeff(2)^2-4*coeff(1)*coeff(3);
  if discr>0
      i1=0;
      i2=0;
      r1=(-coeff(2)+sqrt(discr))/2/coeff(1);
      r2=(-coeff(2)-sqrt(discr))/2/coeff(1);
  else
      r1 =-coeff(2)/2/coeff(1);
      r2 = r1;
      i1 = sqrt(abs(discr))/2/coeff(1);
      i2 = -i1;
  end
end

function y = optimFun(x)
  global polcoeff
  xaug = [1 x];
  [~,r] = deconv(polcoeff,xaug);
  y = norm(r);
end

Using the above function with C1 gives:

Code:

>> roots(C1)

ans =

   1.3386 + 0.0000i
   1.0970 + 0.7570i
   1.0970 - 0.7570i
   0.4657 + 1.2295i
   0.4657 - 1.2295i
  -0.3103 + 1.2423i
  -0.3103 - 1.2423i
  -0.9217 + 0.7964i
  -0.9217 - 0.7964i

>> [xroots,nroots] = slb(C1,1e-20,1000,50)
Optimization ended: relative change in the objective value 
over the last OPTIONS.MaxStallIterations iterations is less than OPTIONS.FunctionTolerance.
Optimization ended: relative change in the objective value 
over the last OPTIONS.MaxStallIterations iterations is less than OPTIONS.FunctionTolerance.
Optimization ended: relative change in the objective value 
over the last OPTIONS.MaxStallIterations iterations is less than OPTIONS.FunctionTolerance.
Optimization ended: relative change in the objective value 
over the last OPTIONS.MaxStallIterations iterations is less than OPTIONS.FunctionTolerance.

xroots =

   -0.9217    0.7964
   -0.9217   -0.7964
    0.4657    1.2295
    0.4657   -1.2295
   -0.3103    1.2423
   -0.3103   -1.2423
    1.0970    0.7570
    1.0970   -0.7570
    1.3386         0


nroots =

     9

Here is a Matlab function that uses a simple-search algorithm for optimizing.

Code:

function [xroots,nroots] = slb_star1(polyCoeffs,SearchLen, MinSearchLen, MaxIter)
% Obtains the roots of a polynomial using obtaimizatio0n-based quadratic
% reduction.
  global polcoeff
  
  if nargin<3, Xinit=[0 0]; end
  if nargin<2, toler=1e-10; end

  rng('shuffle')
  polcoeff = polyCoeffs/polyCoeffs(1);
  order=length(polcoeff)-1;
  xroots=zeros(order,2);
  nroots=0;
  while order>=2    
    x = star(SearchLen, MinSearchLen, MaxIter);
    xaug = [1 x]; % augment x
    [r1,i1,r2,i2]=quadratic(xaug);
    xroots(nroots+1,1)=r1;
    xroots(nroots+1,2)=i1;
    xroots(nroots+2,1)=r2;
    xroots(nroots+2,2)=i2; 
    nroots = nroots + 2;
    order = order - 2;
    [q,~] = deconv(polcoeff, xaug);
    polcoeff = q;      
  end
  
  if order == 1
    xroots(nroots+1,1)=-polcoeff(2)/polcoeff(1);
    xroots(nroots+1,2)=0;
    nroots = nroots + 1;
  end
end

function X = star(SearchLen, MinSearchLen, MaxIter)
% Function bfgs performs multivariate optimization using the
% Broyden-Fletcher-Goldfarb-Shanno method.
%
% Input
% SearchLen - vector for initial seach lengths.
% MinSearchLen - vector of minimum search leangth.
% MaxIter - maximum number of iterations
%
% Output
%
% X - array of optimized variables
%

  N = 2;
  X = [0 0];
  Star = [0 0.5 1;
          0 1 2;
          0 2 3;
          0.5 0.5 4;
          1 1 5;
          2 2 6;
          0.5 0 7;
          1 0 8;
          2 0 9;
          -0.5 0.5 10;
          -1 1 11;
          -2 2 12;
          -0.5 0 13;
          -1 0 14
          -2 0 15;
          -0.5 -0.5 16;
          -1 -1 17;
          -2 -2 18;
          -0.5 0 19;
          -1 0 20;
          -2 0 21];
        
  [nrows,~] = size(Star);
  Xs = X;
  for iter=1:MaxIter
    f0 = optimFun(X);
    Xbest = X;
    irow = 0;
    for i=1:nrows
      for j=1:2
        Xs(j) = X(j) + Star(i,j)*SearchLen(j);
      end
      fx = optimFun(Xs);
      if fx < f0
        Xbest = Xs;
        f0 = fx;
        irow = Star(i,3);    
      end
    end
    
    % No improvement?
    if irow==0
      for i=1:2
        if SearchLen(i) > MinSearchLen(i)
          SearchLen(i) = SearchLen(i)/2;
        end
      end
    else
      X = Xbest;
    end
    
    if SearchLen(1) < MinSearchLen(1) && SearchLen(2) < MinSearchLen(2), break; end
  end
  X = Xbest;
end

function [r1,i1,r2,i2]=quadratic(coeff)
  discr=coeff(2)^2-4*coeff(1)*coeff(3);
  if discr>0
      i1=0;
      i2=0;
      r1=(-coeff(2)+sqrt(discr))/2/coeff(1);
      r2=(-coeff(2)-sqrt(discr))/2/coeff(1);
  else
      r1 =-coeff(2)/2/coeff(1);
      r2 = r1;
      i1 = sqrt(abs(discr))/2/coeff(1);
      i2 = -i1;
  end
end

function y = optimFun(x)
  global polcoeff
  xaug = [1 x];
  [~,r] = deconv(polcoeff,xaug);
  y = norm(r);
end

Here is an example in Matlab. I am using the Matlab function roots() to obtain accurate roots for comparison.

Code:

>> roots(C1)

ans =

   1.3386 + 0.0000i
   1.0970 + 0.7570i
   1.0970 - 0.7570i
   0.4657 + 1.2295i
   0.4657 - 1.2295i
  -0.3103 + 1.2423i
  -0.3103 - 1.2423i
  -0.9217 + 0.7964i
  -0.9217 - 0.7964i

>> [xroots,nroots] = slb_star1(C1,[2 10], [1e-9 1e-9], 1000)

xroots =

    0.4657    1.2295
    0.4657   -1.2295
   -0.3103    1.2423
   -0.3103   -1.2423
   -0.9217    0.7964
   -0.9217   -0.7964
    1.0971    0.7570
    1.0971   -0.7570
    1.3386         0


nroots =

     9

The above examples show how and WHEN te method works. There are many cases, depending on the values of the real polynomial coefficients and the arguments for the root-extracting functions when the above methods give a partially correct answer, whereas Lin-Bairstow gives complete set of correct answers.

[Thomas Klemm]

Cf. Bairstow's Method for the HP-11C

[Namir]

(07-10-2018 12:57 PM)Thomas Klemm Wrote:  Cf. Bairstow's Method for the HP-11C

I have seen and read your message. I must say that the Lin-Bairstow method is quite robust. The optimization-based quadratic extraction method is as good as the optimization sub-task.

I do like your code and the fact that you implemented the LB method for the HP-11C. My hats off to you!

Namir

[ttw]

Using optimization has been around for a while. It's worth examining. One idea that has been used is squaring the objective function and noting that zeros are now minima. I think that non-zero minima may correspond to pairs of complex roots but I haven't done the math. It's useful to separate roots so that local methods may work better. Graeffe's root-squaring method is often used. There are still some useful things being done.


https://arxiv.org/pdf/1501.02168.pdf
https://www.jstor.org/stable/2690148?seq...b_contents
https://people.mpi-inf.mpg.de/~msagralo/...tation.pdf
http://eprints.maths.ox.ac.uk/16/1/mekwi.pdf

[Namir]

Thanks for the links. The Jenkins–Traub algorithm seems to be among the popular polynomial roots finding methods. There is a version of real polynomial coefficients and and other one for complex polynomial coefficients.

Namir

[Namir]

I stumbled on the the Durand-Kerner method which locates the roots of a polynomial simultaneously using an iteration scheme that resembles the Gauss-Seidel method for solving systems of linear equations.

The Durand-Kerner method employs complex mathematics throughout it's calculations. The initial guess, all complex numbers, are rather arbitrary. The seem to work for most cases.

Here is a Matlab implementation.

Code:

function xroots = durand_kerner(polyCoeffs,toler)
% DURAND_KERNER calculates the roots of a polynomial simultaneously. The
% method is relatively simple and easy to implement. The Durand-Kerner
% algorithm resembles the Gauss–Seidel method for solving simultneous linear
% equations.

  n = length(polyCoeffs) - 1;
  xroots = zeros(n,1);
  for i=1:n
    xroots(i) = complex(0.4,0.9)^i;
  end
  bStop = false;
  while ~bStop
    lastRoots=xroots;
    for i=1:n
      prod = 1;
      for j=1:n
        if i~=j
          prod = prod * (xroots(i) - xroots(j));  
        end
      end
      xroots(i) = xroots(i) - polyval(polyCoeffs,xroots(i))/prod;
    end
    
    bStop = true;
    for i=1:n
      if abs(abs(lastRoots(i)) - abs(xroots(i)))>toler
        bStop = false;
        break;
      end
    end
  end
  
end

Here is a simple example:

Code:

>> C1=[1 -3 3 -5];
>> roots(C1)

ans =

   2.5874 + 0.0000i
   0.2063 + 1.3747i
   0.2063 - 1.3747i

>> xroots = durand_kerner(C1,1e-10)

xroots =

   2.5874 - 0.0000i
   0.2063 + 1.3747i
   0.2063 - 1.3747i

[ttw]

The Durand-Kerner method seems useful. I had not known (or had forgotten as I didn't need to solve this problem) of it until now. I did find a couple of (relatively) new papers on the method.
https://www.math.hu-berlin.de/~cc/cc_hom...inding.pdf

https://arxiv.org/pdf/1806.06280.pdf

[Namir]

(07-12-2018 07:22 AM)ttw Wrote:  The Durand-Kerner method seems useful. I had not known (or had forgotten as I didn't need to solve this problem) of it until now. I did find a couple of (relatively) new papers on the method.
https://www.math.hu-berlin.de/~cc/cc_hom...inding.pdf

https://arxiv.org/pdf/1806.06280.pdf

Thank you for the additional pdf links!

Finding the Durand-Kerner method was a real delightful surprise. I implemented the methods in both flavors--the Gauss-Seidel approach and the Jacobi approach. In the first approach, the updated guess for each root is used in updating the next roots (within each iteration). In the second approach, the updated roots are used in the next iteration.

Namir

[Claudio L.]

(07-12-2018 01:05 PM)Namir Wrote:  Finding the Durand-Kerner method was a real delightful surprise. I implemented the methods in both flavors--the Gauss-Seidel approach and the Jacobi approach. In the first approach, the updated guess for each root is used in updating the next roots (within each iteration). In the second approach, the updated roots are used in the next iteration.

Namir

I'm interested in the speed of this method, is it really faster than picking the roots one by one? When I needed a method to implement in newRPL, I didn't find this one but if it's really faster I'd like to see it implemented.
My other question is: does it have guaranteed convergence? or it may fail in some corner cases?
For newRPL I chose the Laguerre's method only because it will always converge to a root no matter how horrible your guess is. In my case this was important because the command that returns all the roots in a polynomial should work all the time, not just every now and then.
But if this method is always convergent and fast, I'm more than interested to replace it.

[ttw]

Take a look at the literature. For polynomials, the method seems pretty good. It's quadratically convergent to simple roots but only linearly convergent to multiple roots. Methods to work around the multiple root convergence are discussed.

What's good is that the method converges to each root and doesn't seem to need much deflation.


There are other methods. It's quite complicated but worthwhile.

[ttw]

Another result on Halley's and Chebychev's methods.

https://arxiv.org/pdf/1807.04476.pdf

[Namir]

(07-13-2018 01:27 AM)ttw Wrote:  Another result on Halley's and Chebychev's methods.

https://arxiv.org/pdf/1807.04476.pdf

ttw, I am hiring you as a research assistant!!! Many thanks for the links.

Namir