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.
Search
Member List
Calendar
Help
/
Using Optimization to Extract Roots of Real Coefficient Polynomials
