Cayley-Hamilton and Generalized Eigenspace Projections

The following 3 x 3 matrix will be used to illustrate the results.  It has eigenvalues 3, -2, -2

A={{-1,-1,0},{1,-3,0},{-4,-1,3}};
A//MatrixForm
p = CharacteristicPolynomial[A,z]//Factor

The following function takes a polynomial in the variable z and calculates the corresponding matrix polynomial with z replaced by the matrix A

MakeMatrixPoly[f_,A_]:=Module[{coeffs,i,n},
coeffs=CoefficientList[f,z];
n=Length[coeffs];
Inner[Times,coeffs,Table[MatrixPower[A,i],{i,0,n-1}]]]

The Cayley-Hamilton theorem says that if we do this with the characteristic polynomial we get the zero matrix.

MakeMatrixPoly[p,A]

The factors g1, g2 of p corresponding to the distinct eigenvalues are used to construct the projection matrices onto the generalized eigenspaces.  Here we set up the equation g1*h1+g2+h2 = 1 and solve it to find the polynomials h1, h2.

g1=(z-3);g2=(z+2)^2;
h1=(a*z+b);h2=c;
g1*h1+g2*h2//Expand
cl = CoefficientList[%,z]
abcsol = Solve[cl=={1,0,0},{a,b,c}][[1]]
h1=h1/.abcsol;h2=h2/.abcsol;
{h1,h2,g1*h1+g2*h2//Expand}

Projection Matrices onto the generalized eigenspaces are the matrix polynomials corresponding to the ordinary polynomials g1*h1 and g2*h2

{g1*h1//Expand,g2*h2//Expand}
proj1=MakeMatrixPoly[g1*h1,A];proj2=MakeMatrixPoly[g2*h2,A];
{proj1//MatrixForm,proj2//MatrixForm}

Here we test the equations that the projections should satisfy: P1+P2 = I, Pi^2=Pi and Pi Pj = 0.

{proj1+proj2//MatrixForm,
proj1.proj1//MatrixForm,
proj1.proj2//MatrixForm}

Jordan Normal Form and Real Canonical Form

By introducing a basis consisting of chains of generalized eigenvectors, we can put our matrix into Jordan form.
The matrix has only two independent eigenvectors.  We will use the one for eigenvalue 3 as part of a basis.

evects = Eigenvectors[A]
v1=evects[[1]]

For the rest of the Jordan basis we need a chain of gev’s.  So we need to find a gev of order 2.  The columns of the projection matrix proj1 will work.   The first column is {1,0,1}  We multiply this by A -(-2) I to get the second vector in the chain (it will be an eigenvector).

v3= {1,0,1};
v2= (A-(-2)*IdentityMatrix[3]).v3

The change of basis matrix has columns v1, v2, v3 and we get Jordan form.

P= {v1,v2,v3}//Transpose;
B= Inverse[P].A.P;
{P//MatrixForm,B//MatrixForm}

Here is another example, with a repeated complex conjugate pair of eigenvalues.

A = {{-3,-5,6,-1},{2,3,-5,2},{0,2,-3,2},{-2,-2,1,-1}};
A//MatrixForm
p= CharacteristicPolynomial[A,z]//Factor
Eigenvalues[A]
Eigenvectors[A]

First we will find the Jordan form by selecting a basis consisting of chains of eigenvectors.  It suffices to find a chain for the eigenvalue  -1+2*i since the other chain will be the complex conjugate.  The simplest way is to solve the equation (A - lambda I)^2 v = 0.

A1=  A-(-1+2*I)IdentityMatrix[4];
A1.A1//Simplify
gevs = NullSpace[%]

The first vector in the gev list is actually an eigenvector.  Let’s take the other vector as the vector v2 of the chain and calculate v1= A1.v2.  Finally we can use the conjugates as the other two vectors of the new basis.  We get the Jordan form.

v2=gevs[[2]];v1=A1.v2; v3= Conjugate[v1];v4=Conjugate[v2];
P={v1,v2,v3,v4}//Transpose;
B= Inverse[P].A.P;
{P//MatrixForm,B//MatrixForm}

To get the real canonical form, we use the real and imaginary parts of these vectors instead.

u1= Re[v1]; w1=Im[v1];u2= Re[v2];w2=Im[v2];
Q= {u1,w1,u2,w2}//Transpose;
Breal= Inverse[Q].A.Q;
{Q//MatrixForm,Breal//MatrixForm}