Unitary Diagonalization and Quadratic Forms

General Diagonalization

Let's review the basic idea of diagonalization by way of an example we have already worked on.  Here is the Leslie model we examined.  Recall that the species of animal we were looking at lived 3 years,  with survival rates  b1 = 0.5,  b2 = 0.25,  and number of offspring for first year animals as a1 = 0,  second year animals  at a2 = 4, and third year animals at a3  = 3.  The population distribution from one year to the next is given by 
	
		  xn+1 = A . xn, 
		  
where xn is the vector of animals sorted into age groups, and A is given below:

A = {{0, 4, 3}, {0.5, 0, 0}, {0, 0, 0.25}}

The idea is to simplify the  matrix equation above by making a simple linear change of variables from variables x to variables  x' with P as the ªchange of basis matrixº, which has  to be invertible, otherwise we couldn't solve for one set of variables in terms of the other:
		x  =  P. x' 
		
This means that   xn =  P . x'n   and  xn+1 =  P . x'n+1.
Now substitute into the first matrix equation and get
	P . x'n+1  =  A . P . x'n,    i.e.,
	     x'n+1  =  P -1 . A . P . x'n.
So what we are really interested in is the 

DIAGONALIZATION THEOREM:  The   n´n  matrix  A  is  diagonalizable,  i.e., there is an invertible matrix  P  such that P -1 . A . P  =  D  is a diagonal matrix, if and only if  A  has a set of  n linearly independent eigenvectors  v1, v2, ... , vn,  in which case we may choose  P = [v1, v2, ..., vn].

Let's test it out on our matrix A:

P = Transpose[Eigenvectors[A]]

Inverse[P].A.P

Unitary Diagonalization

Before we see the crowning glory of this aspect of linear algebra, it is necessary to introduce a bit of terminology:  a square matrix A is called normal if  A . A* = A* . A .
Examples of such matrices are:  any Hermitian matrix A.  Remember that this is a square matrix  A  such that  A = A*.  In this case, the definition of normal becomes  simply
A . A  =  A . A .  Another example is any unitary matrix U.  For in this case,  U* = U-1, which certainly commutes with U under multiplication. (U . U* = I = U* . U)
	Also recall that we say that a matrix is unitarily diagonalizable if there exists a unitary matrix U (orthogonal if A is real)  such that  U* . A . U is a diagonal matrix.  Now for the theorem that  makes this definition worthwhile:

THEOREM:  The square matrix A is unitarily diagonalizable  if and only if A is a normal matrix.

	By the way, part of this is easy to see and gives you an idea where someone might have gotten the idea for the definition of ªnormalityº:  namely, its not hard to check that if  A  is unitarily diagonalizable by unitary  U,  then  A is normal.  The other direction requires more work.

Of course, finding this matrix U is another matter altogether.  The procedure goes something like this:  first, find all the eigenvalues of the matrix  A,  counting multiplicities.  Next,  find an orthonormal basis for each eigenspace.   In case of multiple eigenvalues,  this will involve the Gram-Schmidt Algorithm.  Now simply stack all these vectors together to form the matrix  U;  the theory guarantees us that we will be able to find  a full set of n of these, where  A  is n´n, and that they will automatically form an orthonormal set.  Fortunately, in Mathematica we have some tools at our disposal.  The matrix that results from the following discussion is really an example matrix, lifted out of your textbook, page 321.  MatrixForm it to see that is real symmetric.

a = {{4,2,2},{2,4,2},{2,2,4}}

MatrixForm[a]

Classification of Quadratic Forms

Now here is how the matrix is going to appear in a rather different context:  consider the problem of classifying the quadric surface given by the equation 
		
		4 ( x2  +  y2  + z2  + x y + x z + y z)  = 1.

This is a 3-D version of the classification of conics that you did in Calculus II.  Here is how we will handle such problems.  We treat the left hand side as a ªquadratic formº, that is, a function of the form
		Q(x, y, z)  =  [x, y, z]* . A . [x, y, z].
		
(Check it out,  this is a scalar function of x, y and z.)  But we don't want just any 3´3 matrix A.  We only want to deal with symmetric matrices (Hermitian if we happen to be dealing with complex numbers).  Here's how to find a symmetric  A:  find any  B whatsoever that works,  that is,  the above formula gives the desired quadratic form.  For example,  with the form above, we could use the following  B:

b = {{4,0,0},{4,4,0},{4,4,4}}

See what happens in the following cell:

Simplify[{x,y,z} . b . {x,y,z}]

But rather than use B, use the average of B and its transpose as A.  This matrix will give exacly the same quadratic form as B does and it has the virtue of being symmetric (It is always true that ((B + B*)/2)* = (B + B*)/2.)  Let's check it out.  Notice, by the way, when working with vectors in Mathematica, you don't have to transpose the vector to do multiplication on the left with the matrix B.  Mathematica will do this automatically, and in fact, doesn't even know what the Transpose of a vector is.

a = (b + Transpose[b])/2

Now make sure this gives the same quadratic form:

Simplify[{x,y,z} . a . {x,y,z}]

O.K., now let's get down to business.  First find the eigenvalues of A.  Normally, it would be wise to turn A into a numerical matrix first, since Mathematica handles these in any size.  But this matrix A is sufficiently simple that  Mathematica will  be able to handle it in exact arithmetic.

Eigenvalues[a]

Now let's get the eigenspaces for each eigenvalue.

n = NullSpace[a-8 IdentityMatrix[3]]

This gives us one eigenvector.  Let's name it.

v1 = n[[1]]

n = NullSpace[a-2 IdentityMatrix[3]]

Now name these two eigenvectors

v2 = n[[1]]

v3 = n[[2]]

O.K., O.K., so we worked too hard.  I just wanted to remind you of the definition of the eigenspace.  Here's the lazy (but fast) way of doing the whole thing:

Eigensystem[a]

Notice that the eigenvectors belonging to 2 are linearly independent, but not orthogonal.  Let's turn all these v's into an orthonormal set.  Start with v1.

u1 = v1/Sqrt[v1.v1]

u2 = v2/Sqrt[v2.v2]

u3p = v3 - (v3 . u2) u2

u3 = u3p/Sqrt[u3p . u3p]

Finally, let's put them all together to get the desired orthogonal matrix u.

u = Transpose[{u1, u2, u3}]

Now put it to the tests.

Transpose[u] . u

d = Transpose[u] . a . u

Great!  Now let's see what the significance of this is for the quadratic form we were studying, namely, Q = {x,y,z} . a . {x,y,z}.  Make an orthogonal change of variables given by 

	{x,y,z} =  u . {x', y', z'} ,

so that  

	Q = {x,y,z} . a . {x,y,z} = {x', y', z'} . Transpose[u] . a . u . {x', y', z'}

and see what we get:

Simplify[u . {xp, yp, zp} . a . u . {xp, yp, zp}]

Now we can see that there is an orthonormal basis of R3 which give the axes of the surface Q(x,y,z) = 1, which is nothing more than an ellipsoid with major axis of length 8 and minor axes of length 2 each.  We could even plot it with some more work with Mathematica.

For Hermitian matrices, there is a sneaky shortcut in Mathematica.  This shortcut is based on ideas we don't develop in 314H (maybe we should) about singular values of a matrix. I won't go into detail about it, but try this command:

SingularValues[a]

Here Mathematica choked on symbolic matrices.  So redo it with numbers:

{uu, dd, vv} = SingularValues[N[a]]

Now check this out:

uu . Transpose[uu]

uu . a . Transpose[uu]

Exercises

(1)  For each of the following matrices An, determine whether they are diagonalizable, and if 
       so, find an invertible matrix Sn and a diagonal matrix Dn such that An = Sn-1DnSn:
       
                        æ5   1   0ö
       (a) A1  =  ÷ 0   5   1ç
                        è0   0   5ø
                        
                        æ 1   5   2ö
       (b) A2 =   ç 0   3  -1÷    
                        è 0   0   2ø
                        
                        æ 1  0  2  3 -1ö
                       ÷  3  7  5  8  2 ç 
       (c)  A3 =  ç-1  0  5  2  7 ç
                        ç 4  3  2  7 -1÷
                        è 6  9 -1  3 -8ø

(2)  Repeat the above exercise with Sn being unitary.

(3)  Classify the following quadratic forms (i.e. determine their general shape).
      (a) 3 x2 + 5 x y - 8 y2 = 7
      (b) x y = -1
      (c)  7 x2 - x y - 9 y2 + 7 x - 8 y = 2  (Hint : first make a substitution x' = x + h,
                y' = y + k)
      (d)  7 x2 + 8 y2 - z2 + x y - x z + 3 y z = 0