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

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**

- Consider the Toeplitz matrix
Find its eigenvalues by hand. Check your results with Maple.

- 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 commandA3 := 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.- 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. - 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.
- 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 . - Find an expression sequence for . Are the
eigenvalues sensitive to changes in
**a**?

- By using the routines

Wed Jan 31 11:33:59 EST 1996