Lindsey–Fox algorithm

The Lindsey–Fox algorithm uses the FFT (fast Fourier transform) to very efficiently conduct a grid search in the complex plane to find accurate approximations to the N roots (zeros) of an Nth-degree polynomial. The power of this grid search allows a new polynomial factoring strategy that has proven to be very effective for a certain class of polynomials. This algorithm was conceived of by Pat Lindsey and implemented by Jim Fox in a package of computer programs created to factor high-degree polynomials. It was originally designed and has been further developed to be particularly suited to polynomials with real, random coefficients. In that form, it has proven to be very successful by factoring thousands of polynomials of degrees from one thousand to hundreds of thousand as well as several of degree one million and one each of degree two million and four million. In addition to handling very high degree polynomials, it is accurate, fast, uses minimum memory, and is programmed in the widely available language, Matlab. There are practical applications, often cases where the coefficients are samples of some natural signal such as speech or seismic signals, where the algorithm is appropriate and useful. However, it is certainly possible to create special, ill-conditioned polynomials that it cannot factor, even low degree ones. The basic ideas of the algorithm were first published by Lindsey and Fox in 1992[1] and reprinted in 1996.[2]  After further development, other papers were published in 2003[3][4] and an on-line booklet.[5]  The program was made available to the public in March 2004 on the Rice University web site.[6]  A more robust version-2 was released in March 2006 and updated later in the year.

The strategy implemented in the Lindsey–Fox algorithm to factor polynomials is organized in three stages. The first evaluates the polynomial over a grid on the complex plane and conducts a direct search for potential zeros. The second stage takes these potential zeros and “polishes” them by applying Laguerre's method to bring them close to the actual zeros of the polynomial. The third stage multiplies these zeros together or “unfactors” them to create a polynomial that is verified against the original. If some of the zeros were not found, the original polynomial is “deflated” by dividing it by the polynomial created from the found zeros. This quotient polynomial will generally be of low order and can be factored by conventional methods with the additional, new zeros added to the set of those first found. If there are still missing zeros, the deflation is carried out until all are found or the whole program needs to be restarted with a finer grid. This system has proven to be fast, accurate, and robust on the class of polynomials with real, random coefficients and other similar, well-conditioned polynomials.

Stage one is the reason this algorithm is so efficient and is what sets it apart from most other factoring algorithms. Because the FFT (fast Fourier transform) is used to evaluate the polynomial, a fast evaluation over a dense grid in the complex plane is possible. In order to use the FFT, the grid is structured in polar coordinates. In the first phase of this stage, a grid is designed with concentric circles of a particular radius intersected by a set of radial lines. The positions and spacing of the radial lines and the circles are chosen to give a grid that will hopefully separate the expected roots. Because the zeros of a polynomial with random coefficients have a fairly uniform angular distribution and are clustered close to the unit circle, it is possible to design an evaluation grid that fits the expected root density very well. In the second phase of this stage, the polynomial is evaluated at the nodes of the grid using the fast Fourier transform (FFT). It is because of the extraordinary efficiency and accuracy of the FFT that a direct evaluation is possible. In the third phase of this first stage, a search is conducted over all of the 3 by 3 node cells formed in the grid. For each 3 by 3 cell (see Figure below), if the value of the polynomial at the center node of the cell (the "x") is less than the values at all 8 of the nodes on the edges of the cell (the "o's"), the center is designated a candidate zero. This rule is based on the “Minimum Modulus Theorem” which states that if a relative minimum of the absolute value of an analytic function over an open region exists, the minimum must be a zero of the function. Finally, this set of prospective zeros is passed to the second stage. The number is usually slightly larger than the degree because some were found twice or mistakes were made. The number could be less if some zeros were missed.

Figure 1: Section of the polar coordinate grid showing a 3-node by 3-node cell

Stage two is more traditional than the other two. It “polishes” each of the prospective zeros found by the grid search. The first phase consists of applying an iterative algorithm to improve the accuracy of the location found by the grid search. In earlier versions of the program, Newton's method was used but analysis and experiment showed that Laguerre's method was both more robust and more accurate. Even though it required more calculation than Newton's method for each iteration, it converged in fewer iterations. The second phase of the second stage checks for duplications. A “fuzzy” uniqueness test is applied to each zero to eliminate any cases where on two or more prospective zeros, iterations converged to the same zero. If the number of unique, polished zeros is less than the degree of the polynomial, deflation later will be necessary. If the number is greater, some error has occurred. This stage consumes the largest part of the execution time of the total factorization, but it is crucial to the final accuracy of the roots. One of the two criteria for success in factoring a polynomial is that each root must have been successfully polished against the original polynomial.

Stage three has several phases and possible iterations or even restarting. The first phase of the third stage takes all of the unique, polished zeros that were found in the first two stages and multiplies them together into the coefficient form of a candidate polynomial (“unfactors” the zeros). If the degree of this reconstructed polynomial is the same as that of the original polynomial and if the difference in their coefficients is small, the factorization is considered successful. Often, however, several zeros were missed by the grid search and polish processes of stage one and two, or the uniqueness test discarded a legitimate zero (perhaps it is multiple), so the original polynomial is “deflated” (divided) by the reconstructed polynomial and the resulting (low degree) quotient is factored for the missing zeros. If that doesn't find them all, the deflation process is repeated until they are all found. This allows the finding of multiple roots (or very tightly clustered roots), even if some of them were discarded earlier. If, in the unusual case, these further iterations of deflation do not find all of the missing zeros, a new, finer grid is constructed and the whole process started again at stage one. More details on the third stage are in another module.

Multiple order and clustered roots are unusual in random coefficient polynomials. But, if they happen or if factoring an ill-conditioned polynomial is attempted, the roots will be found with the Lindsey–Fox program in most cases but with reduced accuracy. If there are multiple order zeros (Mth order with M not too high), the grid search will find them, but with multiplicity one. The polishing will converge to the multiple order root but not as fast as to a distinct root. In the case of a cluster with Q zeros that fall within a single cell, they are erroneously identified as a single zero and the polishing will converge to only one of them. In some cases, two zeros can be close to each other in adjacent cells and polish to the same point. In all of these cases, after the uniqueness test and deflation, the quotient polynomial will contain a M − 1 order zero and/or Q − 1 zeros clustered together. Each of these zeros will be found after M − 1 or Q − 1 deflations. There can be problems here because Laguerre's polishing algorithm is not as accurate and does not converge as fast for a multiple zero and it may even diverge when applied to tightly clustered zeros. Also, the condition of the quotient polynomial will be poorer when multiple and clustered zeros are involved. If multiple order zeros are extremely far from the unit circle, the special methods for handling multiple roots developed by Zhonggang Zeng are used. Zeng's method is powerful but designed for polynomials of moderate degrees, and hence only used in special cases [6]. References

Successful completion of the factoring of a polynomial requires matching zeros on the complex plane measured by the convergence of Laguerre's algorithm on each of the zeros. It also requires matching the polynomial reconstructed from the found zeros with the original polynomial by measuring the maximum difference in each coefficient.

Because the FFT is such an efficient means of evaluating the polynomial, a very fine grid can be used which will separate all or almost all of the zeros in a reasonable time. And because of the fineness of the grid, the starting point is close to the actual zero and the polishing almost always converges in a small number of steps (convergence is often a serious problem in traditional approaches). And because the searching and polishing is done on the original polynomial, they can be done on each root simultaneously on a parallel architecture computer. Because the searching is done on a 3 by 3 cell of the grid, no more than three concentric circles of the grid need be kept in memory at a time, i.e., it is not necessary to have the entire grid in memory. And, some parallelization of the FFT calculations can be done.

Deflation is often a major source of error or failure in a traditional iterative algorithm. Here, because of the good starting points and the powerful polisher, very few stages of deflation are generally needed and they produce a low order quotient polynomial that is generally well-conditioned. Moreover, to control error, the unfactoring (multiplying the found roots together) is done in the FFT domain (for degree larger than 500) and the deflation is done partly in the FFT domain and partly in the coefficient domain, depending on a combination of stability, error accumulation, and speed factors.

For random coefficient polynomials, the number of zeros missed by the grid search and polish stages ranges from 0 to 10 or occasionally more. In factoring one 2 million degree polynomial, the search and polish stages found all 2 million zeros in one grid search and required no deflation which shows the power of the grid search on this class of polynomial. When deflation is needed, one pass is almost always sufficient. However, if you have a multiple zero or two (or more) very, very closely spaced zeros, the uniqueness test will discard a legitimate zero but it will be found by later deflation. Stage three has enough tests and alternatives to handle almost all possible conditions. But, by the very definition of random coefficients, it is hard to absolutely guarantee success.

The timings of the Lindsey–Fox program and examples of root distributions of polynomials with random coefficients are here.