Using Optimization to Extract Roots of Real Coefficient Polynomials
07-10-2018, 11:46 AM
Post: #1
 Namir Senior Member Posts: 690 Joined: Dec 2013
Using Optimization to Extract Roots of Real Coefficient Polynomials
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.
07-10-2018, 12:57 PM
Post: #2
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
07-10-2018, 10:40 PM
Post: #3
 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
(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
07-11-2018, 02:28 AM
Post: #4
 ttw Member Posts: 186 Joined: Jun 2014
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
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
07-11-2018, 08:40 AM
Post: #5
 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
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
07-11-2018, 01:44 PM (This post was last modified: 07-11-2018 02:28 PM by Namir.)
Post: #6
 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
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
07-12-2018, 07:22 AM
Post: #7
 ttw Member Posts: 186 Joined: Jun 2014
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
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
07-12-2018, 01:05 PM
Post: #8
 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
(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

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
07-12-2018, 02:05 PM (This post was last modified: 07-12-2018 02:13 PM by Claudio L..)
Post: #9
 Claudio L. Senior Member Posts: 1,671 Joined: Dec 2013
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
(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.
07-12-2018, 07:03 PM
Post: #10
 ttw Member Posts: 186 Joined: Jun 2014
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
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.
07-13-2018, 01:27 AM
Post: #11
 ttw Member Posts: 186 Joined: Jun 2014
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
Another result on Halley's and Chebychev's methods.

https://arxiv.org/pdf/1807.04476.pdf
07-13-2018, 03:21 AM
Post: #12
 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: Using Optimization to Extract Roots of Real Coefficient Polynomials
(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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)