Newton’s Method for Factoring a Polynomial We present Newton's method here as an example of the use of Horner's method. It takes on a particularly straight forward form when applied to polynomials because a formula for the derivative is available. The iterative form of Newton's method is z k + 1 = z k − f ( z ) f ' ( z ) size 12{z rSub { size 8{k+1} } =z rSub { size 8{k} } - { {f \( z \) } over { { {f}} sup { ' } \( z \) } } } {} where the new estimate of the zero location, zk+1 size 12{z rSub { size 8{k+1} } } {}, is calculated from the old value, zk size 12{z rSub { size 8{k} } } {}. A simple derivation is seen from the approximate definition of the derivative {} f ' ( z k ) ≈ f ( z k + 1 ) − f ( z k ) z k + 1 − z k size 12{ { {f}} sup { ' } \( z rSub { size 8{k} } \) approx { {f \( z rSub { size 8{k+1} } \) - f \( z rSub { size 8{k} } \) } over {z rSub { size 8{k+1} } - z rSub { size 8{k} } } } } {} by setting f(zk+1)=0 size 12{f \( z rSub { size 8{k+1} } \) =0} {} and solving for zk+1 size 12{z rSub { size 8{k+1} } } {}. Another approach is to truncate the Taylor's series expansion f ( z + δ ) = f ( z ) + f ' ( z ) δ + f ' ' ( z ) δ 2 / 2 + ⋯ size 12{f \( z+δ \) =f \( z \) + { {f}} sup { ' } \( z \) δ+ { {f}} sup { '' } \( z \) δ rSup { size 8{2} } /2+ dotsaxis } {} with z=zk size 12{z=z rSub { size 8{k} } } {} and z+δ=zk+1 size 12{z+δ=z rSub { size 8{k+1} } } {}. Indeed, this approach explains why the standard Newton's method does not work well with multiple order zeros since f(z) size 12{f \( z \) } {} and f'(z) size 12{ { {f}} sup { ' } \( z \) } {} are both zero at the multiple order zero. Newton's method requires the value of the polynomial and its derivative at eac iteration and Horner's method provides both efficiently. One approach to finding the zeros of f(z) size 12{f \( z \) } {} could be to find the minima of the square of the magnitude of the polynomial. This could be done using a steepest descent algorithm which, when used with an appropriate step size choice, would result in guaranteed convergence, but usually is very slow. Another approach would be to use Newton’s method which has quadratic convergence but often must be started very close to the solution for convergence to occur. A third method would be to use the Gauss-Newton algorithm which is an approximation of the Newton scheme for a squared error measure. This approach also has quadratic convergence but can be more tolerant of starting values. When applied to factoring a polynomial, all three are similar or the same. Moore [15,16] pointed out that the steepest descent algorithm applied to (10) gives the same directional derivative as a Newton type algorithm. Closer examination shows that all three of the approaches suggested above give the same direction, and both the Newton and the Gauss-Newton algorithms have the same step size and, therefore, are exactly the same algorithm. This is true because f(z) size 12{f \( z \) } {} is analytic and the Cauchy-Riemann conditions cause the directions to be the same. These relationships are interesting to us because they means that Newton methods always give a direction that reduces the squared magnitude, even if the step size, which approximates a parabolic shape, may be in considerable error. However, since the direction is correct, there always exists a sufficiently small step size that will reduce the squared magnitude, and, since there are no minima other than zeros of f(z) size 12{f \( z \) } {}, Newton's method with a suitable step size control will always converge to a zero from any starting point. Near a first order zero of f(z) size 12{f \( z \) } {}, the surface is approximately parabolic and the method then converges quadratically. Newton type algorithms do, however, have trouble with multiple zeros or near a cluster of zeros where the surface is not approximately parabolic. A simple section of Matlab code which implements Horner's method (also called synthetic division or nested evaluation and described in the appendix) to evaluate f(z) size 12{f \( z \) } {} and f'(z) size 12{ { {f}} sup { ' } \( z \) } {} at each iteration and uses them for updating the estimate of the zero location, z size 12{z} {}, is given in Figure 1.

N = length(a)          % Number of poly coefficients
while (abs(f)>e)    % Newton's iterations
  f = a(1);  fp = 0;   % Horner initial values
  for k = 2:N          % Horner's algorithm
    fp = f + fp*z;     % to eval f'(z)
    f = a(k) + f*z;    % and f(z)
  end
  z = z - f/fp;        % Newton's update of z
end

Matlab code for Newton's Method where z size 12{z} {} must be given as an initial estimate and e size 12{e} {} is the value of f(z) size 12{f \( z \) } {} that is assumed to be close enough to zero. For N size 12{N} {} coefficients, a(k) size 12{a \( k \) } {}, each iteration requires approximately 2N size 12{2N} {} multiplications and the Newton's formula requires one division. In a practical algorithm, measures must be taken to prevent dividing by zero or having overflow or divergence. A good method for starting the algorithm must be used. The practical program discussed later is built around this simple combination of Horner's and Newton's methods.