Illustrating Cayley-Hamilton with Python

Posted on by John

If you take a square matrix M, subtract x from the elements on the diagonal, and take the determinant, you get a polynomial in x called the characteristic polynomial of M. For example, let

M = \left[ \begin{matrix} 5 & -2 \\ 1 & \phantom{-}2 \end{matrix} \right]

Then

\left| \begin{matrix} 5-x & -2 \\ 1 & 2-x \end{matrix} \right| = x^2 - 7x + 12

The characteristic equation is the equation that sets the characteristic polynomial to zero. The roots of this polynomial are eigenvalues of the matrix.

The Cayley-Hamilton theorem says that if you take the original matrix and stick it into the polynomial, you’ll get the zero matrix.

\left[ \begin{matrix} 5 & -2 \\ 1 & \phantom{-}2 \end{matrix} \right]^2 - 7\left[ \begin{matrix} 5 & -2 \\ 1 & \phantom{-}2 \end{matrix} \right] + 12\left[ \begin{matrix} 1 & 0 \\ 0 & 1\end{matrix} \right] = \left[ \begin{matrix} 0 & 0 \\ 0 & 0\end{matrix} \right]

In brief, a matrix satisfies its own characteristic equation. Note that for this to hold we interpret constants, like 12 and 0, as corresponding multiples of the identity matrix.

You could verify the Cayley-Hamilton theorem in Python using scipy.linalg.funm to compute a polynomial function of a matrix.

>>> from scipy import array
>>> from scipy.linalg import funm
>>> m = array([[5, -2], [1, 2]])
>>> funm(m, lambda x: x**2 - 7*x + 12)

This returns a zero matrix.

I imagine funm is factoring M into something like PDP−1 where D is a diagonal matrix. Then

f(M) = P f(D) P−1.

This is because f can be applied to a diagonal matrix by simply applying f to each diagonal entry independently. You could use this to prove the Cayley-Hamilton theorem for diagonalizable matrices.

Related posts

4 thoughts on “Illustrating Cayley-Hamilton with Python”

  1. Saurish Chakrabarty

    Nice…

    The dot product function seems to use integers when they can (you get zero exactly).
    >>> np.dot(m,m)-7*m+12*np.eye(2,dtype=int)
    If needed, one could rewrite a version of funm for polynomials using dot products.

  3. Max

    Also possible: m@m – 7*m + 12*np.eye(2). The function evaluating a polynomial P = [a0, …, an] = a0 + a1*X + …+an*X**n for a matrix could be
    def matpol(P,M):
    Mk=np.eye(np.ndim(M));
    S=P[0]*Mk
    for k in range(1,len(P)):
    Mk = Mk @ M
    S+=P[k]*Mk
    return S

    With this, matpol([12,-7,1],m) yields the zero array.

Comments are closed.