Using Mathematica to Explor Gerschgorin Circles and Cassini Ovals

by
Robert Piziak
Danny W. Turner
Department of Mathematics
Baylor University
Waco, Texas 76798-7328

Mathematica in Education
Vol.3 No.3
Summer 1994
(c) TELOS/Springer-Verlag

Introduction

In this article, we use the graphics power of Mathematica to explore the geometric approach to the estimation and location of the eigenvalues of a real or complex matrix.  It is hoped that instructors in Linear Algebra and Numerical Analysis courses will be able to use this paper to implement the real time graphing power of Mathematica to give students a way to visualize how the eigenvalues of various kinds of matrices are distributed.  While our focus is on what might be of use in the classroom, there are some practical issues at hand here as well.  One might ask why bother with rough geometric estimations when eigenvalues and even eigenvectors are just the push of a button away.  Most of our students now carry TI-85 calculators and indeed this is truely the case.  However, those familiar with current research in numerical linear algebra know the story is not quite so simple. Though the sophomore linear algebra texts typically teach the characteristic polynomial as the way to compute eigenvalues, it is well known that for large matrices, this polynomial is highly unstable and essentially useless for finding eigenvalues.  For small to medium size matrices, packages are available with subroutines for solving eigenvalue problems(usually using the QR algorithm).  For large problems, there are a number of different iterative methods that are used(see [Saa91]).  For some of these approaches, initial estimates of the eigenvalues are useful.  Examples include shift-invert Krylov subspace methods and preconditioning methods (see [MoS86]).  It is interesting that rough estimates of eigenvalues are also needed for some techniques for solving large systems of linear equations(see [Ash91] for polynomial preconditioning and [HaY81] for Chebychev iteration).   The eigenvalue estimates can be obtained from the geometric approach we are about to describe and would be tedious or even impossible to obtain without the help of a system like Mathematica.  Since this article is for a general audience, we begin with theory concerning Gerschgorin circles and illustrate with a series of examples.  Then results of Brauer are stated and a number of examples involving Cassini ovals are presented.

Gerschgorin's Circle Theorem

The  Russian mathematician S. Gerschgorin published in [Ger31] an interesting way to estimate the eigenvalues of a square matrix with real or complex entries.  Simply stated, he showed that the eigenvalues of such a matrix must lie in the union of certain disks of the complex plane.  The proof of this result is remarkably easy, relying only on the triangle inequality.
  More precisely, let A = [aij] be an n-by-n matrix over the reals or complexes.  The radii of our circles (or disks) are constructed by summing the absolute values (magnitudes) of the entries of each row of A omitting the diagonal element.  Symbolically, the ith radius ri is 
               n
          ri = · |aij|  .
              j=1
              j\= i 
              
         
Let Di be the (closed) disk in the complex plane C centered at aii with radius ri.  This is the ith Gerschgorin disk; in symbols                 
 
                   Di  =  {z in C |  |z - aii| ² ri}. 
   
The boundaries of these disks are the circles
     
                   Ci = {z in C |  |z - aii| = ri}, 
     
called Gerschgorin circles.

Theorem:Gerschgorin's Circle Theorem.

Let A be an n-by-n matrix over the reals or complexes.  Then each eigenvalue of A is contained in at least one of the disks Di.  In other words, if b1,b2,...,bp are the distinct eigenvalues of A, then the set {b1,b2,...,bp} is contained in the union of D1,D2,...,Dn.  Also, if the ith component of an eigenvector has maximal modulus, then the corresponding eigenvalue must lie in the ith disk.  Moreover, if a set of k disks does not intersect the remaining n-k disks, exactly k eigenvalues must be contained in the set.  In particular, if one disk is disjoint from all the others, that disk contains exactly one eigenvalue.

The proof requires no advanced material and may be found in several of the references at the end of this paper [see for example [Hag88], [Ste73], or [SS90]].

Let's set up some simple Mathematica functions to facilitate our exploration of Gerschgorin circles.  We will display the circles rather than the disks.  The first function gcircles finds the circles(centers and radii) with a choice of color.  The second function listgcircles generates a table of the circles so that we may see the centers and radii. The last function displayer can be used to display graphics objects in a convenient form.  We take advantage of the package Colors.

Needs["Graphics`Colors`"]

gcircles[A_,color_]:= Module[{absA=Abs[A],r},
   r = Apply[Plus,Transpose[absA] - absA *
                IdentityMatrix[Length[A]]];
  Table[{color, Circle[{Re[A[[i,i]]], Im[A[[i,i]]]}, r[[i]]]},
  {i, 1, Length[A]}]]
  
listgcircles[A_]:= Table[gcircles[A,Red][[i,2]],
  {i, Length[A]}] //TableForm
  
displayer[graphs_]:= Show[graphs,
  AspectRatio -> Automatic,Axes->Automatic,
  AxesOrigin -> {0,0},
  DisplayFunction -> $DisplayFunction];

Let's look at a simple example.  Note the necessity of using Graphics with the gcircles function so the Show command in displayer will operate on graphics objects.

A={{3,0,-1,-1/4,1/4},{0,5,1/2,0,1},
   {-1/4,0,6,1/4,1/2},{0,-1,1/2,-3,1/4},
   {1/6,-1/6,1/3,1/3,4}};
MatrixForm[A]
listgcircles[A]
p1=Graphics[gcircles[A,Red]];
displayer[p1]

There is an isolated circle here and our theory says it must contain a single eigenvalue.  We can illustrate the theorem by having Mathematica estimate and plot the eigenvalues of A.  The function eigenpoints creates the points for the eigenvalues with a choice of point thickness.

eigenpoints[A_,thick_]:= 
   Module[{evals = Eigenvalues[N[A]]},
   		Table[{{PointSize[thick],
   		        Point[{Re[evals[[i]]], Im[evals[[i]]]}]}},
        {i, 1, Length[A]}]];

Now let's repeat the example above with eigenvalues superposed.    Again note the use of Graphics with eigenpoints before we can display the output.

Print["Eigenvalues of A: ", TableForm[Eigenvalues[N[A]]]]
p2  =Graphics[eigenpoints[A, .02]];
displayer[{p1, p2}];

As further illustration, we can compute the characteristic polynomial of A and superimpose its graph on the picture above.  We will also include an inset of the situation at a different scale to get a better view of the complete graph of the polynomial.  Here is a function to find the characteristic polynomial of a (square) matrix.

charpoly[x_,A_]:= 
	Det[x*IdentityMatrix[Length[A]]-A] /;
  		MatrixQ[A] && Apply[Equal, Dimensions[A]]

We now compute the characteristic polynomial of the matrix in the example given above and superpose all the graphics, including the inset.

p[x_]=charpoly[x,A]

Show[p1,p2,
		Plot[p[x],{x,-6,12},
  			DisplayFunction->Identity,
  			PlotRange->{-4,4}],
  	DisplayFunction->$DisplayFunction,
  	AspectRatio->Automatic,Axes->Automatic,
  	Prolog -> Rectangle[{6.5,0},{12,5},
  			{{Plot[p[x]/500, {x,-6,8},
  				PlotRange->{-4,4},
  				PlotLabel->"p(x)/500",
  				DisplayFunction->Identity,
  				Frame->True]},
  	p1,p2}]]

Our functions work on complex matrices just as well.  We will use a randomly generated 6-by-6 complex matrix to illustrate.

A=Table[Random[Complex,{-10-10I,10+10I},2],{6},{6}];
MatrixForm[A]
Print["Eigenvalues of A: ",TableForm[Eigenvalues[N[A]]]]
p1=Graphics[gcircles[A,Green]];
p2=Graphics[eigenpoints[A,.02]];
displayer[{p1,p2}]

Note that Gerschgorin's theorem does not tell us in which disks eigenvalues actually lie.  Indeed, an interesting fact that can be easily illustrated is that some disks may be devoid of eigenvalues.  Consider the following example.

A={{4,-1},{2,2}};
MatrixForm[A]
Print["Eigenvalues of A: ",TableForm[Eigenvalues[N[A]]]]
p1=Graphics[gcircles[A,Black]];
p2=Graphics[eigenpoints[A,.03]];
displayer[{p1,p2}]

Observe the disk with no eigenvalues!  However, a minor change in A
will produce no empty circles.

A={{4,-1},{-1,2}};
MatrixForm[A]
Print["Eigenvalues of A: ",TableForm[Eigenvalues[N[A]]]]
p1=Graphics[gcircles[A,Violet]];
p2=Graphics[eigenpoints[A,.025]];
displayer[{p1,p2}]

Is it a coincidence that there are no empty circles when A is symmetric?  The reader is invited to show(investigate) that real symmetric(complex hermitian) matrices cannot have empty Gerschgorin circles.

Next we observe that column sums work just as well as row sums in Gerschgorin's theorem.  With the notation as above,construct new radii
          n
     cj = ·|aij|         and new disks
          i=1
          i\= j
          
     D'j = {z in C |  |z - ajj| ² cj}.
     
  The conclusions of Gerschgorin's theorem apply equally to these disks.  Indeed, since det(A-bI) = det(AT-bI), A and AT share the same eigenvalues where A is of course square.  Since the row disks for AT are the same as the column disks for A, our claim follows.   In fact one can sometimes sharpen the estimates for an eigenvalue b by computing both the row and the column disks.  We illustrate this below.

A={{1,.2},{.1,2}};
MatrixForm[A]
p1=Graphics[gcircles[A,Orange]]
displayer[p1]

We can argue that the real matrix A must have real eigenvalues since the two circles above are disjoint.  Recall, complex eigenvalues occur in complex conjugate pairs when A is real.  Just by looking at the scale we see we can bound the eigenvalues by .8²b1²1.2 and 1.9²b2²2.1.  Next we compute the column disks and superpose the graphs to see if we can sharpen these bounds.

B=Transpose[A];
MatrixForm[B]
p2=Graphics[gcircles[B,Magenta]]
displayer[p2]

displayer[{p1,p2}]

Now we see from graph p2 that .9²b1²1.1 and 1.8²b2²2.2 so looking at the combined results above, we may conclude .9²b1²1.1 and 1.9²b2²2.1.

  There are many interesting issues yet to be explored concerning Gerschgorin circles.  For example, we can ask how to use the Gerschgorin circles to give a bound on the spectral radius of a matrix.  For comparison purposes we define the function spectralradius to compute the spectral radius of a matrix using Mathematica's built-in eigenvalue routine.  Next we define  the function specradbound that finds the maximum over Gerschgorin circles of distance of center from the origin plus the Gerschgorin radius of this circle to obtain a computationally easy rough bound on the spectral radius.

spectralradius[A_]:= Max[Abs[Eigenvalues[N[A]]]];

specradbound[circs_]:= Max[Table[
  Sqrt[circs[[1,i,1]][[1]]^2+circs[[1,i,1]][[2]]^2]+
  circs[[1,i,2]],{i,1,Length[circs[[1]]]}]]

Again we illustrate with an example.

A={{1,2,3},{2,-40*I,5*I},{10,-20,30}};
MatrixForm[A]
displayer[Graphics[{gcircles[A,Red],eigenpoints[A,.03]}]]
Print["Eigenvalues of A: ",TableForm[Eigenvalues[N[A]]]]
Print["Spectral Radius of A: ",spectralradius[A]]
Print["Upper Bound for the Spectral Radius of A: ",
  specradbound[listgcircles[A]]]

The reader will find it quite interesting to modify (one at a time) entries in the matrix A above and note the effect on the resulting circles and eigenvalues.  With some experimentation one can discover how to manipulate entries to move the circles and eigenvalues in a predictable fashion!

Our next illustration looks at the problem of getting better bounds on eigenvalues by shrinking disks via conjugation with a diagonal matrix.  The function conjmatrix computes DAD-1 for a given matrix A where D is a diagonal matrix taken with strictly positive diagonal entries.  Note that DAD-1 and A are similar and thus have the same eigenvalues.

conjmatrix[A_,v_]:=
  DiagonalMatrix[v].A.Inverse[DiagonalMatrix[v]]

Once again, a good example is called for.

A={{1,.1,.1},{.1,1,.1},{.1,.1,2}};
MatrixForm[A]
listgcircles[A]
displayer[Graphics[gcircles[A,Red]]]

Our first choice is D = diagonal[1,1,a] with a variable parameter a.

Clear[a]
m1=conjmatrix[A,{1,1,a}];
MatrixForm[m1]
listgcircles[m1]

The algebraic form of the radii clearly explains what is happening.  To make the radius of the circle about two smaller, we decrease the value of the parameter a.  But this in turn has the unfortunate effect of increasing the radii of the disks centered about one, as the graphics will now clearly illustrate.  We will generate a graphics array to make the point but an animation could easily be constructed here as well.

graphs=Table[Flatten[Graphics[gcircles[m1,Red],Axes->Automatic,
  AspectRatio->Automatic,AxesOrigin->{0,0}]],{a,1.,1/4,-1/4}];
displayer[GraphicsArray[Partition[graphs,2]]]

The attentive reader may have noticed that there should be three circles above.  Indeed there are, but let's try a variation in our diagonal matrix so as to bring this out.

Clear[a]
m2=conjmatrix[A,{1,a,a}];
MatrixForm[m2]
listgcircles[m2]
graphs=Table[Flatten[Graphics[gcircles[m2,Red],Axes->Automatic,
  AspectRatio->Automatic,AxesOrigin->{0,0}]],{a,1.,1/4,-1/4}];
displayer[GraphicsArray[Partition[graphs,2]]]

This choice of D now shows us there is a single eigenvalue in the circle centered at 2 and the other two eigenvalues lie in the union of the two circles centered at 1.

If matrix A is not square, it does not have eigenvalues.  However, it does have singular values and we can estimate(or bound) these noting that the singular values of A are the square roots of the eigenvalues of AAT or ATA.  We illustrate with a random rational matrix.

A=Table[Random[Integer,{-10,10}] /
         (Random[Integer,{-10,10}]+11),
       {4},{8}];
B=A.Transpose[A];
MatrixForm[B]
circtable=listgcircles[B]//N
bd=specradbound[circtable];
Print["Eigenvalues of AAT: ",Chop[Eigenvalues[N[B]]]]
Print["Upper Bound for the spectral radius of AAT: ",bd]
Print["Upper Bound for the singular values of A: ",Sqrt[bd]]
  
displayer[Graphics[{gcircles[B,Blue],eigenpoints[B,.03]}]]

Note that from the circles alone we can deduce that the singular values of A lie in the interval from zero to the square root of our bound on the spectral radius of AAT.

Brauer Regions and Cassini Ovals

In [Fro69], Froberg summarizes results of Brauer in [Brau47] who developed another interesting way to graphically "capture" the eigenvalues of an n-by-n matrix A.  Define the row-based Cassini region Ojk as follows:
  
     Ojk = {z in C |  |z - ajj| |z - akk| ² rjrk}.
     
Each Ojk is an example of a Cassini oval(confer [JaB68]) along with its interior.  Brauer established that the union of the n(n-1)/2 Cassini regions Ojk , 1²j²(n-1) and j<k²n, contains the eigenvalues of A.  Similarly, define a column-based Cassini region by

     O'jk = {z in C |  |z - ajj| |z - akk| ² cjck},
     
and note that the union of these regions(with integers j and k as above) also contains the eigenvalues of A.  It follows that the eigenvalues of A must be contained in the intersection of the union of the row-based Cassini regions and the union of the column-based Cassini regions.

    Below is the function cassovals which can be used to generate the desired Cassini ovals.

cassovals[A_,xmin_,xmax_,ymin_,ymax_,color_,npts_]:= 
  Module[
    {z=x+I*y,dimA=Length[A],absA=Abs[A],r},r=
    Apply[Plus,Transpose[absA]-absA*IdentityMatrix[Length[A]]];
    Table[ContourPlot[
    Abs[z-A[[j,j]]]*Abs[z-A[[k,k]]]-r[[j]]*r[[k]],
    {x,xmin,xmax},{y,ymin,ymax}, Contours->{0.},
    AspectRatio->Automatic,Axes->Automatic,
    ContourShading->False,ContourSmoothing->Automatic,
    ContourStyle->{{color}},
    PlotPoints->npts,DisplayFunction->Identity],
    {j,1,dimA-1},{k,j+1,dimA}]
    ];

Let's try the row-based regions for a complex matrix A.

A={{4+I,-1,3+2I},{2,2*I,I+1},{5-8I,9-7I,3I}};
MatrixForm[A]
p1=cassovals[A,-12,12,-12,12,NavyBlue,10];
displayer[p1]

Now, we compute the eigenvalues of A and display them along with the Cassini ovals and the Gerschgorin circles.  Notice that the ovals are "tighter" than the circles.  The order of the items in the list argument of displayer below is important in order for the final plot to be scaled correctly!!  In particular, the contour graphics should come first.

p2=Graphics[eigenpoints[A,.06]];
p3=Graphics[gcircles[A,Red]];
displayer[{p1,p2,p3}]

Continuing with the A above, it would be interesting to add the column-based collection of Cassini ovals.

B=Transpose[A];
p4=cassovals[B,-11,15,-11,14,LimeGreen,10];
displayer[{p1,p2,p3,p4}]

The idea of conjugation is also applicable to Brauer regions.  Below we reconsider the matrix used in our earlier conjugation example.  Note that the diagonal matrix is different from the previous illustrations and only one value of parameter a is used.  However, the basic idea and results are the same.  When counting regions, remember that an "oval of Cassini" may in fact consist of two distinct ovals.

Clear[a]
A={{1,.1,.1},{.1,1,.1},{.1,.1,2}};
p1=cassovals[A,.8,2.1,-.2,.2,Blue,12];
B=conjmatrix[A,{a,1,a}];
MatrixForm[B]
a=2;
p2=cassovals[B,.8,2.1,-.2,.2,Red,12];
p3=Graphics[eigenpoints[A,.015]];
displayer[{p1,p2,p3}]

To see more clearly  what's going on, we zoom in a little on each cluster of ovals found above.  First, we look more closely at the cluster of small ovals near x position 2.

p4=cassovals[A,1.9,2.1,-.2,.2,Blue,12];
p5=cassovals[B,1.9,2.1,-.2,.2,Red,12];
displayer[{p4,p5,p3}]

Next, let's examine the cluster of ovals near x position 1.

p6=cassovals[A,.8,1.2,-.2,.2,Blue,12];
p7=cassovals[B,.8,1.2,-.2,.2,Red,12];
displayer[{p6,p7,p3}]

The next matrix A is the one-step transition matrix of a 4-state Markov chain.  In the theory of Markov chains, a number of properties of the chain can be deduced from information about the eigenvalues of A.  For this example, the eigenvalue 1 appears only once.  It is not a coincidence that the chain has just one irreducible closed subclass(states 1 and 2)!
  For this example, the Cassini ovals are also circles and the plot has superposed the Gerschgorin circles, the eigenvalues of A, and the Cassini ovals.  This example shows that eigenvalues may fall on the boundaries of the regions.

A={{0,1,0,0},{1,0,0,0},{.1,0,0,.9},{0,.3,.7,0}};
MatrixForm[A]
p1=Graphics[gcircles[A,Green]];
p2=Graphics[eigenpoints[A,.03]];
p3=cassovals[A,-1,1,-1,1,SkyBlue,6];
displayer[{p3,p2,p1}]

Let's look at a very interesting matrix based on an example in [Var65].  Notice how we may use the specradbound function to help determine a reasonable region of (x,y) values for input into our cassovals function.

A={{1,1,0,0},
   {0,Exp[Pi*I/2],1,0},
   {0,0,Exp[Pi*I],1},
   {1,0,0,Exp[3*Pi*I/2]}}//N;
MatrixForm[A]
Print["Upper Bound for the spectral radius of A: ",
  specradbound[listgcircles[A]]]
p4=cassovals[A,-2,2,-2,2,SkyBlue,12];
p5=Graphics[eigenpoints[A,.025]];
p6=Graphics[gcircles[A,SeaGreen]];
displayer[{p4,p5,p6}]

Finally, here are the row-based Cassini regions and eigenvalues for an A that was previously considered.  The Gerschgorin circles are included for reference.  The choice of the value of the argument "npts" in the cassovals function controls PlotPoints in the ContourPlot routine and can have much influence on the zero level contours actually displayed.

A={{3.,0,-1,-1/4,1/4},{0,5,1/2,0,1},
   {-1/4,0,6,1/4,1/2},{0,-1,1/2,-3,1/4},
   {1/6,-1/6,1/3,1/3,4}};
MatrixForm[A]
p4=cassovals[A,-4,7,-1.3,1.3,SkyBlue,15];
p5=Graphics[eigenpoints[A,.025]];
p6=Graphics[gcircles[A,SeaGreen]];
displayer[{p4,p5,p6}]

With these illustrations and examples as a starting point, we hope that the reader is motivated to continue exploring the properties of eigenvalues via these beautiful geometrical ideas.

Acknowledgements

There does not seem to be a single best source for a treatment of Gerschgorin circles.  In this note we have partially followed the treatment given by Hagar.  The results of Brauer were taken from Froberg's summary.  However, the other sources listed in the bibliography were also quite useful and informative in the preparation of this paper.  We especially thank our colleague Ron Morgan for his helpful suggestions.

References

[Ash91] Ashby,Steven F.,"Minimax Polynomial Preconditioning for Hermitian Linear Systems", SIAM J. Matrix Anal. Appl.,Vol 12, No.4, 766-789,October, (1991).

[Brau47] Brauer,A., "Limits for the characteristic roots of matrices II", Duke Math. J. 14, 21-26 (1947).

[Fro69] Froberg,C.E., Introduction to Numerical Analysis, 2nd ed, Addison-Wesley, Massachusetts (1969).

[Ger31] Gerschgorin,S.A.,"Uber die Abgrenzung der Eigenwerte einer Matrix.", Izv. Akad. Nauk SSSR,Ser. Fiz.-Mat. 6,749-754 (1931).

[Hag88] Hagar, William W., Applied Numerical Linear Algebra, Prentice Hall, New Jersey (1988).

[HaY81] Hageman, L.A. and D.M.Young, Applied Iterative Methods, Academic Press, New York (1981).

[JaB68]James,R.C. and E. F. Beckenbach, James & James Mathematics Dictionary, 3rd ed, Multilingual ed, D. Van Nostrand, New Jersey (1968).

[MeV68] Medley,H.I. and R.S. Varga, "On Smallest Isolated Gerschgorin Disks for Eigenvalues III.", Numerische Mathematik 11,361-369 (1968).

[MoS86] Morgan,R.B. and D. S. Scott, "Generalizations of Davidson's Method for Computing Eigenvalues of Sparse Symmetric Matrices", SIAM J. Sci. Stat. Comput., Vol 7, No, 3, 817-825,July, (1986).

[Saa92] Saad, V., Numerical Methods for Large Eigenvalue Problems, Halsted Press, New York, (1992).

[Ste73] Stewart, G.W., Introduction To Matrix Computations,Academic Press,New York (1973).

[SS90] Stewart, G.W. and Ji-guang Sun, Matrix Perturbation Theory, Academic Press, Harcourt Brace Jovanovich, New York (1990).

[Tau49] Taussky, O.,"A Recurring Theorem on Determinants.", American Mathematical Monthly 46, 672-675 (1949).

[Var64] Varga, R.S.,"On Smallest Isolated Gerschgorin Disks for Eigenvalues.", Numerische Mathematik 6, 366-376 (1964).

[Var65] Varga, R.S., "Minimal Gerschgorin Sets.", Pacific Journal of Mathematics 15, 719-729 (1965).

About the Authors

Robert Piziak was born on 17 January 1943 in Hadley, Massachusetts.  He receved his B.A. in Mathematics from Amherst College (1964) and his M.A. and Ph.D. from the University of Massachusetts (1966, 1969).  He has taught at the University of Florida, Gainesville, Centre College, Danville,Ky., and has been at Baylor University since 1981.

Danny W. Turner was born in North Carolina in 1947.  He received both his B.S.(1969) and Ph.D.(1973) from Clemson University.  Since 1973, Turner has been in the mathematics department at Baylor University, where he is currently Professor of Mathematics and Associate Dean of the College of Arts and Sciences. His main teaching and research interests are in probability, statistics and applied mathematics.