next up previous
Next: Symbolic Computation With Up: Applied Math 511b Notes Previous: Test File for

Eigenvalue Problems.

Consider the problem of computing --- and, for example, plotting --- the eigenvalues of the matrix

If we wrote a Fortran program that accepted a numerical value for b and then computed the eigenvalues of A by using, for example, some standard library package for computing eigenvalues, we would consider this a numerical approach. In Maple, we could use evalf(Eigenvals(A)) to numerically compute the eigenvalues of a numeric matrix A.

If instead we first found an `explicit formula' (whatever that means) for the eigenvalues, in terms of the parameter b, and then wrote a program to evaluate those eigenvalues using the formula, then we would consider this a symbolic approach.

The distinguishing feature between `numerics' and `symbolics' thus seems to be that numerical methods evaluate objects to numbers as soon as possible, while symbolic methods do so as late as possible (this was articulated by J. Rice at ISSAC '92). Deeper thought shows that the distinction is really artificial (a numeral is, after all, a symbol itself), and that there is a whole spectrum of methods available, depending on when quantities are evaluated.

Let us return to the matrix A. We can use Maple to compute explicit expressions for the eigenvalues by finding the characteristic polynomial and then finding its roots, as we were all taught to do in our first linear algebra class.

> with(linalg):
Warning: new definition for   norm
Warning: new definition for   trace
> A := matrix([[ 0, 1], [1, 2*b]]);

> p := charpoly(A,lambda);

> answer := solve(p,lambda);

On the numerical end of the spectrum we can write (for example) a Matlab program to evaluate the eigenvalues given a numerical value of b:

function lambda = eiga(b)
A = [ 0  1
      1  2*b];
lambda = eig(A);
We can use either the Maple output or the Matlab program for subsequent analysis or plotting. For example, we can plot the eigenvalues for e.g. using either representation.

Which approach is `better'? Surprisingly, neither is totally dominant. The symbolic approach really isn't explicit --- it involves taking a square root, an infinite algebraic operation, which we are really just `used to' --- we have come to accept finding square roots as `basic operations'. However, the computer implementation of `square root' is really quite complicated. ( Geometrically, of course, the idea is simple and even has a finite straightedge-and-compass construction --- hence, perhaps, our ingrained sense of `explicitness' can be rationalized).

But in Matlab, taking eigenvalues is a basic operation! Hence, in some sense, the Matlab program is also `explicit'! One just has to `get used to the idea'. So there is no clear philosophical superiority of one approach over the other.

Are there practical superiorities? Yes. The explicit formula is on the one hand more efficient than the numerical method, and on the other hand less stable than the numerical method (for large values of b).

The last point is not as well-known as it deserves to be. We illustrate it by noting that the product of the two eigenvalues is identically -1, for all values of b. When we plot , (where each piece is calculated as displayed), for large b we get the plot in Figure 2.

  
Figure 2: Numerical instability of the quadratic formula.

We clearly see that something goes wrong for b between and --- in fact, the explicit formula is giving incorrect answers when evaluated (the rounding error made when calculating is revealed by catastrophic cancellation of the leading digits --- see [21] for details).

I emphasize that this shows symbolic expressions themselves can be unstable to evaluate --- and when this happens we don't care how fast they evaluate because they give the wrong answer.

In contrast, the Matlab program can be shown to be stable in the backward error sense: eig(A) will compute the exact eigenvalues of a matrix very near to A. Further, one can show that the eigenvalues of A are insensitive to small changes in A --- in numerical analysis parlance, the eigenvalues of A are well-conditioned [47].

Remark. Numerical experiments show that Matlab will in fact produce exact eigenvalues of a matrix of the form

where is less than the machine epsilon and is usually small relative to b if b > 1 and small absolutely if b < 1. This is an example of an experimental `structured backward error result'. Structured backward error analysis is currently a topic of much research interest. See [60] for an introduction.

Can we, by judicious use of the art of scientific computing, create a program with the efficiency of the symbolic expression with the stability of the numerical eigenvalue calculation? The answer is yes, but the cost is the replacement of the explicit expressions for the eigenvalues with a `computation sequence' --- i.e. a small, `straight-line' program! Assuming that b is real, the following (well-known) computation sequence does the trick:

No catastrophic cancellation occurs here to reveal the rounding errors made. This computation sequence (in words, `compute the root with the largest magnitude first, since that can be done accurately in floating-point arithmetic, and then use that one to compute the small root') is efficient and stable. It is, if you like, a hybrid symbolic/numeric `recipe'.

But what about insight? Remember Hamming's famous motto, `The Purpose of Computing is Insight, not Numbers'! Which formulas give the most insight? Well, each contributes something. The quadratic formula has a symmetry that the computation sequence lacks, the computation sequence allows easier computation of the asymptotics of as , and the eigenvalue formulation retains important information on the sensitivity of the problem to changes in the data. Depending on the problem context, any one of the three formulations could give the required insight. Of course, a good analyst would be aware of all three (and perhaps more) and switch between representations as appropriate.

Final Remarks on this Example. Many ideas are discussed in this book using elementary examples. I will rarely even bother to say `the generalization of these ideas to larger examples is obvious'. In this case perhaps the generalizations do need emphasis.

The problem of finding eigenvalues of symmetric or unsymmetric n by n matrices is well studied. The naive approach, the first learned in a linear algebra course and hence the first thing one tries to use in practice, is to first find the characteristic polynomial and then to try to find its roots. If n is less than 5, then there are explicit formulas in terms of radicals, analogous to the quadratic formula (see the excellent book [6] for simple derivations, and see the exercises). These formulas suffer from instability, just as the quadratic formula does [67]. For n=5, there is an `explicit' formula in terms of elliptic functions [38] --- which can also be unstable. There are no formulas in terms of radicals for , and hence numerical methods must be used.

It has been known for more than thirty years that this naive approach of finding eigenvalues by first finding the characteristic polynomial is inherently unstable, and very bad in a practical sense (see the excellent expository essay by Wilkinson in [104]). Just because we can compute the characteristic polynomial doesn't mean that we should do so as a step to find the roots. Sadly, people still use this method (see e.g. [17] as an example of an otherwise adequate, even good, text that advocates this method). Perversely, part of the attraction of symbolic systems is the thought that this horrid approach can be used in practice, with very high precision used to find the roots of the characteristic polynomial to sufficient accuracy --- that is, one tries to `buy more accuracy by paying for more precision'. This fixes the problem, all right, but the cost can be horrendous. See the exercises.

Exercises

  1. Consider the Toeplitz matrix

    Find its eigenvalues by hand. Check your results with Maple.

  2. Consider the n-by-n generalization of the matrix in question 1:

    A convenient way of entering this matrix in Maple is to use the following short procedure:

    A := n -> linalg[toeplitz]([seq(a^k,k=0..n-1)]);
    
    See ?seq and ?toeplitz for a discussion of what these commands do individually. Here, seq constructs the sequence of elements , and this is enclosed in list brackets (square brackets) since toeplitz expects a list as input, from which it constructs the Toeplitz matrix given by that list of elements. Issuing the command
    A3 := A(3);
    
    will reproduce the matrix of question 1, and assign it to the variable A3 for later manipulation. Similarly A(5) generates the 5-by-5 version of this matrix, etc.
    1. By using the routines charpoly and factor, find the eigenvalues of . Is Maple's answer an `explicit' formula or is it really a computation sequence? See ?label for a discussion of how you can control the intermediate quantities %1, %2, etc.
    2. Find the characteristic polynomial of . Factor it. Given the `explicit' formula for fifth-degree polynomials in terms of elliptic functions mentioned in the text, you see that we could find explicit formulas for the roots of this polynomial, but personally I doubt whether that would be useful.
    3. If a=0, then the characteristic polynomial of is . Thus we expect some difficulty in numerically computing the roots of the characteristic polynomial for small a and large n, because near-multiplicity of the roots makes them sensitive to changes in a (or, indeed, roundoff error). For , use realroot (see ?realroot for usage) to isolate all ten real roots of the characteristic polynomial of , to an accuracy of . (The roots are guaranteed to be all real because the matrix is symmetric). Try to use fsolve or some other numerical rootfinder to get the same results. Compare with direct computation of the eigenvalues using evalf(Eigenvals(A10)). Repeat for .
    4. Find an expression sequence for . Are the eigenvalues sensitive to changes in a?




next up previous
Next: Symbolic Computation With Up: Applied Math 511b Notes Previous: Test File for



Robert Corless
Wed Jan 31 11:33:59 EST 1996