>> x=poly([1e-50, 1, 1e50]);
>> roots(x) % The two smaller magnitude roots underflow to 0.
ans =
1.000000000000000e+050
0
0
>> lroots(x) % Finds all the roots
ans =
1.000000000000000e+050
1.000000000000000e+000
9.999999999999999e-051
>> u=lroots([1 0 0 0 -1]).'; % The fourth roots of unity. >> x=poly([1e20*u 1e40*u]);Roots underflows the 4 smaller roots:
>> roots(x) ans = -9.999999999999980e+039 -5.439066599736671e+023 +1.000000000000001e+040i -5.439066599736671e+023 -1.000000000000001e+040i 1.000000000000000e+040 0 0 0 0>> lroots(x) % Does not underflow.
>> sort(lroots(x)) ans = 0 -1.000000000000000e+020i 1.000000000000000e+020 0 +1.000000000000000e+020i -1.000000000000000e+020 0 -1.000000000000000e+040i 1.000000000000000e+040 0 +1.000000000000000e+040i -1.000000000000000e+040
The first attempt used the final correction term when polishing a root as an error estimate. This is almost always an excellent indication of the error. Unfortunately it is sometimes many orders of magnitude too small. When this strategy was tested on the Wilkinson polynomial, it worked every time with one exception. Once, the actual error in the polished position was .0004 but the final correction term was 1e-40. When polishing roots of 1,000-degree polynomials it has sometimes happened that the final correction term was also 1e-40 for a root with magnitude nearly 1. This was surely not a good estimate of the root's accuracy.
An alternative method was needed to estimate the error in polished roots. It is well known that there is a way to esimtate the error when using Horner's method to evaluate a polynomial, f(z). Therefore it is also possible to estimate the error in evaluating the derivative, f'(z). The Newton correction term is -f(z)/f'(z) and since error estimates are obtainable for both the numerator and the denominator, we can obtain an error estimate for the ratio. This is the method that is currently used to obtain error estimates on the roots. It tends to be larger than the final polishing correction.
The current strategy needs more validation and perhaps new ideas. Getting good error estimates of the polished root positions is critically important for ill-conditioned polynomials. All factorization programs first get tentative values for all the roots which are then polished against the original polynomial. It is quite common for two initial positions to (approximately) polish to the same position. This gives a false double root which must be corrected by removing one of the polished values. Function uniq(z,err) takes a fuzzy array of complex values, z, with error estimates, err, and removes duplicates.
For example tentative positions were obtained for the roots of the Wilkinson polynomial. Unfortunately, when they were polished against the original polynomial, two of the polished positions were 11.995 and 12.003. This was a false double root for which the correct value was 12. It is difficult to discern that they are false double roots unless you possess reasonably accurate error estimates for each root.
Error estimates for the roots are also used in the phase unwrapping function, unwrapz. If the distance of a root to the unit circle is less than the error estimate for the root, then the root is assumed to lie exactly on the unit circle. For example,
>> x=[1 zeros(1,14) -1];
>> [z,e]=lroots(x);
For half the roots, abs(z) differs from 1.0 by 2^(-53)=1.11e-16; i.e., there is an error in the least significant bit. Every root has estimated error 1.11e-15 so unwrapz will assume all the roots lie on the unit circle, which in this case is correct.
Of course a root can be extremely close to the unit circle and ill-conditioned so the error estimate is large. It will be assumed to lie on the unit circle which may or may not be true. However it is impossible to tell where it actually lies relative to the unit circle. Assuming it lies on the unit circle averages the possibilities. If it actually lies on one side, it causes one wrap while if it lies on the other side it does not wrap at all. Assuming that it lies on the unit circle means it causes half a wrap.
Error estimates for the roots are also used in iterativeDeflation - see below.
A common cause of failure is to have all the missing roots be close together. For example in the polynomial in degree200.zip, a previous version of lroots failed to find 14 roots and all had arguments between 124.9 and 145.6 degrees. The polynomial associated with these 14 roots was approximately ((z-r)*(z-conj(r)))^7 where r has argument 135 degrees. This is a very ill-conditioned polynomial; small changes in the coefficients make large changes in the roots. Unfortunately there were large errors in the computation of this 14 degree polynomial. There were small errors in the 186 roots that were found but that is not the big proglem. The big problem is instability in doing the deflation. The numerator alternates sign every other sample and has a dynamic range of 1.4e7. The denominator alternates sign every other sample and has a dynamic range of 8.1e8. Their ratio has a dynamic range of 1e3. This implies a loss of several significant digits and this caused further problems because the quotient is ill-conditioned.
This problem is not limited to the situation where the original polynomial had a huge dynamic range and was ill-conditioned. This problem can happen with random coefficient polynomials which are extremely well-conditioned and have a small dynamic range. For instance.
>> x=rand(1,1001);
>> [junk,k]=sort(angle(z));
Split z into two subsets, a and b.
>> a=z([1 7 12 16 21 30 40]);
If we put a and b back together and unfactor, it reconstructs x quite well
>> polycmp(x, unfactor([a;b],1))
Now separately unfactor a and b and polynomially multiply the results. The answer should reconstruct x but it does a terrible job. Polycmp returns values between 0 and 2, with 0 meaning the two polynomials are scalar multiples of each other and 2 meaning they are not remotely close.
>> u=unfactor(a,1);
>> plot(x)
Note the dynamic ranges of x, u and v; especially v.
>> max(abs(x)) / min(abs(x))
>> max(abs(u)) / min(abs(u))
>> max(abs(v)) / min(abs(v))
We are supposed to convolve two functions with dynamic ranges of 10^3 and 10^14 and get a function with a dynamic range of 10^4. The only way to do that is to loose about 10 significant digits of accuracy. Function unfactor works in the Fourier domain. If you instead do the unfactoring in the time domain with Matlab's poly, u is the same but v is significantly different. However their convolution also does not at all reconstruct x.
This is precisely the potential problem with factorization. Since the original polynomial was random coefficient, almost all roots are very close to the unit circle and uniformly distributed around it. Thus every root that is close to the unit circle is also very close to a dozen other roots. Suppose initially we found the 986 roots in b. We unfactor them to get v which has a dynamic range of 1.8e14. We divide that into the original random coefficient polynomial to get a polynomial that is approximately ((z-r)*(z-conj(r)))^7. The division produces a polynomial that is inherently ill-conditioned and the division incurrs a tremendous loss of significant digits. This polynomial is factored and the roots polished against the original polynomial. It is normal for them all to polish to a neaby root that had been found previously. Thus iterativeDeflation will not be able to find any more roots and the factorization will fail.
>> z=lroots(x);
After deflation, seconds =0.641 Return code=0 Error=4.83169e-012
>> z=z(k); % Sort by increasing angle.
>> z=z(find(imag(z)>=0)); % Restrict to the upper half plane.
>> b=z([2:6, 8:11, 13:15, 17:20, 22:29, 31:39, 41:end]);
ans = 1.174615960053416e-013
>> v=unfactor(b,1);
>> c=conv(u,v);
>> polycmp(x,c)
ans = 1.52859548744714
>> plot(c)
>> plot(u)
>> plot(v)
ans = 1.979039103039599e+004
ans = 3.400382863410701e+003
ans = 1.851515824403872e+014
Neither test alone is sufficient. To see why the polish test, by itself, is insufficient, take a correct set of roots and replace one of the roots with a second copy of some other root. This new set will pass the polish test but it has a false double root. To see why the unfactor test by itself is insufficient, look at the polynomial in degree270.zip which can be downloaded from this web site.
>> load degree270.mat
>> [r,p]=broots(coeff); % Factor it using broots.
iterativeDeflation needs to find 99 missing roots.
Return code=12 Error=2.33288e-010.
"R" is the raw roots obtained by deflating one root at a time (and its complex conjugate). "P" is the roots obtained by polishing r, removing 99 false duplicates, and then handing the results to iterativeDeflation to find missing roots.
Let's also get Matlab's factorization.
>> m=roots(coeff);
A plot shows that all 3 are nearly equal for 46 roots (and their complex conjugates) that are near the unit circle. However, away from the unit circle, the polished roots, p, are nearly equal to Matlab's roots, m, but the raw roots, r, are dramatically different.
>> plotroots(p,m,r)
>> legend('Polished','Matlab','Raw')
"R" does not pass the polish test. When r is polished, the result is a set with many false double roots. "P" passes the polish test, the roots do not move much if they are re-polished.
>> plotroots(r, laguerre3(coeff,r)) % r fails the polish test
>> plotroots(p, laguerre3(coeff,p)) % p passes the polish test
Even though m and p are both close to the correct set of roots and r is very far from it, r does the best job of passing the unfactor test. The unfactor test by itself can not be trusted.
>> polycmp(coeff, unfactor(r,1))
ans = 1.1824e-014
>> polycmp(coeff, unfactor(m,1))
ans = 4.2577e-014
>> polycmp(coeff, unfactor(p,1))
ans = 2.9774e-010
In summary, p passes both the polish and the unfactor tests. R passes the unfactor test much better than p, but r completely fails the polish test. In order for a set of putative roots to be declared correct, they should pass both tests.
There are two major commercial suites of mathematical software. The provider of one of them was willing to answer a question and they confirmed that after a factorization they do not polish the roots against the original polynomial. The other provider was not willing to discuss the issue but the calling sequence for their routines are available on the web and their factorization programs do not have an argument which allows the user to turn on/off a post-factorization polish. Therefore they likely do not do a final polish either. Most people using their routines are factoring polynomials of degree < 15 and their factorizations are almost certain to be nearly correct. The roots could be improved by a final polish but for most users this would be overkill. However, as the polynomial in degree270.mat shows, for higher degree polynomials it is essential that steps be taken to validate a factorization - the roots must pass both the polish and the unfactor tests.