Wiener–Hopf and spectral factorization of real polynomials by Newton’s method

We present a new method for factoring a real polynomial into the product of two polynomials which have their zeros inside and outside the unit circle, respectively. The approach is based on solving a nonlinear system by Newton’s method. The Jacobi matrix is a polynomial of a companion matrix and thus a Bezoutian times the inverse of a triangular Toeplitz matrix. After getting deeper understanding of the displacement structure of the Bezoutian, each Newton step can be reduced to a few applications of discrete Fourier or cosine transforms and the LU-decomposition of a conceived linear system with a Cauchy-like matrix. These structural features are employed to design a fast algorithm, which shows excellent numerical behavior. For example, in the case of the extremely ill-conditioned test polynomial of the degree 2n=20000, we obtain the factorization after 22 Newton steps with an error of 3.4·10-8 in the Euclidean norm, the execution time on a laptop being a couple of minutes. The local quadratic convergence of our Newton method is proved and explicit convergence balls are given on the basis of estimates for the Lipschitz constant which occurs in Kantorovich’s theorem. We also design and test two superfast versions of our method, which, for polynomials of degree 20000, typically need about half a minute to produce an initial coefficient vector and then perform the Newton iteration within 1 s.