Pejorative Manifolds of Polynomials and Matrices, part 1

In an unpublished 1972 technical report "Conserving confluence curbs ill-condition," Velvel Kahan coined the descriptive term pejorative manifold. In case you don't use it in everyday conversation, pejorative means "expressing contempt or disapproval."

Velvel's report concerns polynomials with multiple roots, which are usually regarded with contempt because they are so badly conditioned. But Velvel's key observation is that although multiple roots are sensitive to arbitrary perturbation, they are insensitive to perturbations that preserve multiplicity. The same observation can be applied to matrix inverses and eigenvalues.

Contents

The manifold

In this example, the pejorative manifold $\mathcal{M}$ is the set of all sixth-degree polynomials which have a zero of multiplicity three at x = 3. These are severe restrictions, of course, and $\mathcal{M}$ is a tiny subset of the set of all polynomials. But if we stay within $\mathcal{M}$, life is not too bad.

Two polynomials

First a pair of polynomials with multiple roots. All of our polynomials have a root of multiplicity three at x = 3 as well as other multiple roots. This one a double root at x = 1 and a simple root at x = 2.

 x = sym('x');

 p1 = (x-3)^3*(x-2)*(x-1)^2;

$$ p1 = {\left(x-1\right)}^2\,\left(x-2\right)\,{\left(x-3\right)}^3 $$

The next one has a double root at x = 2 and a simple root at x = 1.

 p2 = (x-3)^3*(x-2)^2*(x-1);

$$ p2 = \left(x-1\right)\,{\left(x-2\right)}^2\,{\left(x-3\right)}^3 $$

Let's see the coefficients.

  p1 = expand(p1);

$$ p1 = x^6-13\,x^5+68\,x^4-182\,x^3+261\,x^2-189\,x+54 $$

  p2 = expand(p2);

$$ p2 = x^6-14\,x^5+80\,x^4-238\,x^3+387\,x^2-324\,x+108 $$

To summarize, we have two polynomials, p1 and p2. Both have degree six and have a triple root at x = 3. In addition, pk has a double root at x = k .

Confirm

We can solve the two exactly and verify that they have the desired roots.

  z1 = solve(p1)'
  z2 = solve(p2)'
z1 =
[ 1, 1, 2, 3, 3, 3]
z2 =
[ 1, 2, 2, 3, 3, 3]

Convex combination

Now let's take a convex linear combination. I am using 1/3 and 2/3, but I could use any other weights as long as their sum is 1. Any such convex linear combination is still in $\mathcal{M}$.

   p3 = 1/3*p1 + 2/3*p2;

$$ p3 = x^6-\frac{41\,x^5}{3}+76\,x^4-\frac{658\,x^3}{3}+345\,x^2-279\,x+90 $$

This still has the triple root at x = 3. The roots at x = 1 and x = 2 are now simple and the root at x = 5/3 is the convex linear combination of 1 and 2 with the same coefficients, 5/3 = 1/3*1 + 2/3*2.

   z = solve(p3)'
z =
[ 1, 5/3, 2, 3, 3, 3]

Plots

Looking at plots of the three polynomials, you can appreciate how the triple root at 3 is more sensitive than the blue double root at 1 or the green double root at 2, which are, in turn, more sensitive and any of the simple roots.

   plot_polys(p1,p2,p3)

Floating-point haze

Now a small perturbation. Take a convex linear combination with an irrational weight and use vpa, the variable precision floating point arithmetic, set to 25 digits. We are in a floating-point haze with thickness $10^{-25}$ surrounding $\mathcal{M}$.

   digits(25)
   w = 2/(1+sqrt(vpa(5)))
   q = w*p1 + (1-w)*p2;
   z = w*1 + (1-w)*2
w =
0.6180339887498948482045868
z =
1.381966011250105151795413

Here are the coefficients.

   disp(fliplr(coeffs(q))')
                         1.0
 -13.38196601125010515179541
  72.58359213500126182154496
 -203.3900966300058885005431
  309.1277174175132491262221
 -240.5654115187641954923808
  74.62616460750567819695231

Find the roots.

   z = solve(q)
z =
                                      0.9999999999999999999999944
                                       1.381966011250105151795437
                                       1.999999999999999999999981
                                       2.999999961118614061145771
 3.000000019440692969427115 - 0.00000003367226800807395156577582i
 3.000000019440692969427115 + 0.00000003367226800807395156577582i

The simple roots have survived to full accuracy. But the triple root at x = 3 has been split by the perturbation into three complex numbers on a circle centered at x = 3 with radius roughly the cube root of the working precision.

   circle(double(z(4:6)-3))

Where did the cube root of precision come from? We are trying to solve the equation

$$ {\left(x-3\right)}^3\,\sigma(x)=0 $$

where $\sigma(x)$ is not zero near $x$. Using floating point arithmetic with precision $\epsilon$, we instead solve something like

$$ {\left(x-3\right)}^3\,\sigma(x)=\epsilon $$

The solution is

$$ x = 3 + \left( \frac{\epsilon}{\sigma(x)}\right)^{1/3}$$

with the cube root taking on three different complex values.

Venture off the manifold

Subtract a small quantity from the constant term in p1. This moves us out of $\mathcal{M}$.

   p = p1 - sym(1.e-4)
p =
x^6 - 13*x^5 + 68*x^4 - 182*x^3 + 261*x^2 - 189*x + 539999/10000

We can try to find the roots exactly, but the modified polynomial does not factor over the rationals.

   solve(p)
ans =
 root(z^6 - 13*z^5 + 68*z^4 - 182*z^3 + 261*z^2 - 189*z + 539999/10000, z, 1)
 root(z^6 - 13*z^5 + 68*z^4 - 182*z^3 + 261*z^2 - 189*z + 539999/10000, z, 2)
 root(z^6 - 13*z^5 + 68*z^4 - 182*z^3 + 261*z^2 - 189*z + 539999/10000, z, 3)
 root(z^6 - 13*z^5 + 68*z^4 - 182*z^3 + 261*z^2 - 189*z + 539999/10000, z, 4)
 root(z^6 - 13*z^5 + 68*z^4 - 182*z^3 + 261*z^2 - 189*z + 539999/10000, z, 5)
 root(z^6 - 13*z^5 + 68*z^4 - 182*z^3 + 261*z^2 - 189*z + 539999/10000, z, 6)

So, find the roots numerically with vpa. This time all of the roots, even the simple ones, are affected.

   digits(20)
   z = vpa(solve(p))
z =
                           0.99647996935769229447
                            1.0035512830254182854
                            1.9999000099960012995
 2.9856883694814744941 - 0.025815296081878984073i
 2.9856883694814744941 + 0.025815296081878984073i
                            3.0286919986579391324

We constructed p1 to have a double root at x = 1, a simple root at x = 2, and a triple root at x = 3. Their multiplicity determines how much the perturbation affects them.

   e = abs(double(z - z1'))
e =
   0.003520030642308
   0.003551283025418
   0.000099990003999
   0.029516982906352
   0.029516982906352
   0.028691998657939

Double precision

Let's leave the Symbolic Toolbox for our final example. Now polynomials are represented by conventional MATLAB vectors. Construct the convex combination with double precision floating point and then use the traditional MATLAB function roots.

  format long
  q1 = sym2poly(p1)
  q2 = sym2poly(p2)
  q3 = 1/3*q1 + 2/3*q2;
  disp('q3 = ')
  fprintf('%25.15f\n',q3)
  z = roots(q3)
q1 =
     1   -13    68  -182   261  -189    54
q2 =
     1   -14    80  -238   387  -324   108
q3 = 
        1.000000000000000
      -13.666666666666664
       76.000000000000000
     -219.333333333333314
      345.000000000000000
     -279.000000000000000
       90.000000000000000
z =
  3.000049783405092 + 0.000086219822788i
  3.000049783405092 - 0.000086219822788i
  2.999900433188159 + 0.000000000000000i
  2.000000000002659 + 0.000000000000000i
  1.666666666665640 + 0.000000000000000i
  1.000000000000021 + 0.000000000000000i

This is the same behavior we had with 32-digit vpa. The simple roots at x = 1, 5/3, and 2 are computed to nearly full double precision accuracy, while the triple roots have accuracy roughly eps^(1/3), which is about five digits.

   circle(z(1:3)-3)

To be continued

I am making this a two-parter. The next post is about matrix eigenvalues.

   close

References

Zhonggang Zeng, The Tale of Polynomials, PowerPoint presentation. <http://homepages.neiu.edu/~zzeng/neiu.ppt> , 2003.

W. Kahan, Conserving confluence curbs ill-condition. Technical Report 6, Computer Science, University of California, Berkeley, 1972.




Published with MATLAB® R2018b

|

Comments

To leave a comment, please click here to sign in to your MathWorks Account or create a new one.