# Verify Cayley-Hamilton Theorem

Consider the matrix:

A = {{-1, -4, -2}, {0, 1, 1}, {-6, -12, 2}}


Now, the characteristic polynomial is found:

p[λ_] = CharacteristicPolynomial[A, λ]


Now, how can I substitute the matrix A into the polynomial to verify that the answer is the zero matrix?

• Maybe p[lambda] /. lambda^k_. -> MatrixPower[A,k]. Commented Apr 4, 2015 at 19:42
• Did not seem to work: In[147]:= p[[Lambda]] Out[147]= -2 + [Lambda] + 2 [Lambda]^2 - [Lambda]^3 In[149]:= p[lambda] /. lambda^k_. -> MatrixPower[A, k] Out[149]= {{0, -2, -2}, {-2, 0, -2}, {-2, -2, 0}} Commented Apr 4, 2015 at 20:42
• Right. I missed a step. You need to remove the constant term c, then replace it with c*IdentityMatrix[...]. Commented Apr 4, 2015 at 22:13

To verify that matrix is a zero of its characteristic polynomial, The Characteristic polynomial of the matrix is found, then evaluated for the matrix. The result should be the zero matrix.

Clear[x]
a = {{-1, -4, -2}, {0, 1, 1}, {-6, -12, 2}};
n = Length[a];
p = CharacteristicPolynomial[a, x];
(Sum[ Coefficient[p, x, i] MatrixPower[a, i], {i, 0, Exponent[p, x]}]) // MatrixForm


another way:

cl = CoefficientList[p, x];
Sum[MatrixPower[a, j - 1] cl[[j]], {j, 1, Length[cl]}]


• As usual, an excellent answer. I think that the second part will be the easiest for my students to understand. Thanks for you help. Commented Apr 4, 2015 at 21:32
• Here's the Horner way to do it: Fold[#1.a + #2 IdentityMatrix[n] &, ConstantArray[0, {n, n}], Reverse[CoefficientList[CharacteristicPolynomial[a, x], x]]]. Commented May 2, 2015 at 16:31

Mathematica actually has a function purpose-built for the operation you're looking for.

MatrixFunction[f, m] gives the matrix generated by the scalar function f at the matrix argument m.

In your case, MatrixFunction[p, A] will return the 3-by-3 zero matrix, as desired.

Here is a more "adventurous" way to prove the Cayley-Hamilton theorem that in my opinion has a lot of educational value because it re-derives the characteristic polynomial while at the same time verifying the theorem:

linearCombination =
First@Solve[
A.A.A == α IdentityMatrix[3] + β A + γ A.A, {α, β, γ}]

(* ==> {α -> -2, β -> 1, γ -> 2} *)

(α IdentityMatrix[3] + β A + γ A.A - A.A.A == 0 A) /. linearCombination

(* ==> True *)


In words, I ask Mathematica to check if we can express $A^3$ as a linear combination of lower powers of $A$. Indeed, it finds the required coefficients $\alpha$, $\beta$ and $\gamma$, thereby giving a constructive proof that the claim is true. The Cayley-Hamilton theorem is now verified (in this example) by checking that the matrix polynomial I just found has as its roots exactly the eigenvalues of $A$:

Table[(α + β a + γ a^2 - a^3 == 0) /.
linearCombination, {a, Eigenvalues[A]}]

(* ==> {True, True, True} *)


Therefore, the above polynomial is (by definition) the characteristic polynomial. So here I turned the traditional steps of the proof upside down.

• thank you for this neat answer+1:) Commented Apr 5, 2015 at 7:48

This is a somewhat high-brow way of showing the Cayley-Hamilton theorem, through the power of holomorphic functional calculus.

As I mentioned in this answer, one of the standard ways to define a matrix function is through a Cauchy-like construction:

$$f(\mathbf A) = \frac{1}{2\pi i} \oint_\gamma f(z)\, (z \mathbf I- \mathbf A)^{-1}\,\mathrm dz$$

where the closed contour $\gamma$ encloses the eigenvalues of $\mathbf A$. Here is how to use this definition for a Mathematica demonstration of Cayley-Hamilton.

First, define the matrix and the corresponding characteristic polynomial:

A = {{-1, -4, -2}, {0, 1, 1}, {-6, -12, 2}};
p[λ_] = CharacteristicPolynomial[A, λ];


To construct the integrand, do this:

n = Length[A];
ip = Simplify[p[z] Inverse[z IdentityMatrix[n] - A]];


Now, to evaluate the contour integral, we can use the Cauchy integral theorem (and relatedly, the residue theorem). Let's look at what ip looks like:

ip
{{-14 + 3 z - z^2, 4 (-8 + z), 2 (1 + z)},
{6, 14 + z - z^2, -1 - z},
{6 (-1 + z), 12 (-1 + z), 1 - z^2}}


Could it be?

And @@ Flatten[Map[PolynomialQ[#, z] &, ip, {2}]]
True


GASP! All the matrix entries are polynomials!

This means that the conditions of the integral theorem are satisfied (polynomials certainly are holomorphic functions!), and the matrix contour integral should evaluate to the zero matrix. If you still want to try things out in Mathematica anyway, we can use the syntax for Integrate[] that allows piecewise linear paths in the complex plane:

Integrate[ip, {z, 1, I, -1, -I, 1}]/(2 Pi I)
{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}


and we're done.

See also this math.SE question and this paper.