Skip to content

Characteristic Polynomial

This section focuses only on square matrices, i.e., the linear transformation corresponding to matrix \(A\) maps \(n\) vectors to \(n\) vectors.

In practical problems, we often need to consider repeatedly applying the same transformation. Describing it solely as "the linear transformation corresponding to matrix \(A\) transforms the identity matrix \(I\) into \(A\)" can be abstract. In such cases, the best approach is to find "fixed points," i.e., parts that remain unchanged under the transformation.

However, in reality, the linear transformation corresponding to matrix \(A\) may have no fixed points. Thus, we settle for finding collinear parts or parts that undergo simple deformations.

Eigenvalues and Eigenvectors

Under the linear transformation corresponding to matrix \(A\), some vectors only change in scale but not in direction.

Let \(V\) be a linear space over field \(F\), and \(T\) be a linear transformation on \(V\). If there exists \(\lambda \in F\) and a nonzero vector \(\xi \in V\) such that:

\[ T\xi=\lambda\xi \]

then \(\lambda\) is called an eigenvalue of \(T\), and \(\xi\) is called an eigenvector of \(T\) belonging to eigenvalue \(\lambda\).

Eigenvectors lie on the same line and maintain their direction under the linear transformation (compression to zero is also considered as maintaining direction). Eigenvectors are not unique; any vector collinear with an eigenvector is also an eigenvector. However, the zero vector is not considered an eigenvector, as vectors with direction are naturally nonzero. The eigenvalue of an eigenvector is the scaling factor.

In practical applications, for eigenvectors with the same eigenvalue, a set of basis vectors is usually chosen to represent them all.

Let \(\alpha_1,\alpha_2,\cdots,\alpha_n\) be a basis of \(V\), and let \(A\) be the matrix of \(T\) under this basis, i.e.:

\[ T(\alpha_1,\alpha_2,\cdots,\alpha_n)=(\alpha_1,\alpha_2,\cdots,\alpha_n)A \]

Let \(\lambda_0\) be an eigenvalue of \(T\), and \(\xi\) be an eigenvector of \(T\) belonging to \(\lambda_0\), with a nonzero vector \(X\) satisfying:

\[ \xi=(\alpha_1,\alpha_2,\cdots,\alpha_n)X \]

Then:

\[ T\xi=\lambda_0\xi \]
\[ T(\alpha_1,\alpha_2,\cdots,\alpha_n)X=\lambda_0(\alpha_1,\alpha_2,\cdots,\alpha_n)X \]
\[ (\alpha_1,\alpha_2,\cdots,\alpha_n)AX=\lambda_0(\alpha_1,\alpha_2,\cdots,\alpha_n)X \]
\[ AX=\lambda_0X \]
\[ (A-\lambda_0I)X=0 \]

Thus, the corresponding determinant is also \(0\).

Characteristic Polynomial

Consider an \(n\times n\) matrix \(A\), where \(n\geq 0\land n\in\mathbb{Z}\). Let \(\lambda\) be a parameter. The matrix \(\lambda I-A\) is called the characteristic matrix of \(A\).

The determinant of the characteristic matrix is called the characteristic polynomial of \(A\). It expands to an \(n\)-degree polynomial whose roots are the eigenvalues of \(A\), denoted as \(p_A(\lambda)\):

\[ p_A(\lambda)=\det(\lambda I_n-A)=\begin{vmatrix} \lambda-a_{11} & -a_{12} & \cdots & -a_{1n} \\ -a_{21} & \lambda-a_{22} & \cdots & -a_{2n} \\ \vdots & \vdots & & \vdots \\ -a_{n1} & -a_{n2} & \cdots & \lambda-a_{nn} \\ \end{vmatrix} \]

where \(I_n\) is an \(n\times n\) identity matrix. Some definitions use \(p_A(\lambda)=\det(A-\lambda I_n)\), which differs from our definition only by a sign \((-1)^n\). However, our definition ensures \(p_A(\lambda)\) is monic, while the alternative is monic only when \(n\) is even. Note that the determinant of a \(0\times 0\) matrix is well-defined as \(1\).

The nonzero solution vector \(X\) corresponding to \((\lambda_0 I-A)X=0\) is called an eigenvector of \(A\) belonging to \(\lambda_0\).

The linear transformation \(T\) has eigenvalue \(\lambda_0\) if and only if matrix \(A\) has eigenvalue \(\lambda_0\).

The linear transformation \(T\) has eigenvector \(\xi\) if and only if matrix \(A\) has eigenvector \(X\), where:

\[ \xi=(\alpha_1,\cdots,\alpha_n)X \]

According to the fundamental theorem of algebra, the characteristic polynomial can be factored as:

\[ f(\lambda)=|\lambda I-A|={(\lambda-\lambda_1)}^{d_1}\cdots{(\lambda-\lambda_m)}^{d_m} \]

Here, \(d_i\) is called the algebraic multiplicity of eigenvalue \(\lambda_i\). The sum of all algebraic multiplicities equals the dimension \(n\) of the space.

Finding All Eigenvalues and Eigenvectors of a Matrix

Follow these steps:

  • Compute the determinant \(|\lambda I-A|\).
  • Find all roots of the polynomial \(f(\lambda)=|\lambda I-A|\) in field \(F\), i.e., the eigenvalues of \(A\).
  • For each eigenvalue \(\lambda\) of \(A\), solve the homogeneous linear system \((\lambda I-A)X=0\) to obtain a basis \(X_1,\cdots,X_t\). Then all eigenvectors of \(A\) belonging to \(\lambda\) are:
\[ k_1X_1+k_2X_2+\cdots+k_tX_t \]

where \(k_i\) are not all zero.

  • The eigenvectors of linear transformation \(T\) belonging to \(\lambda\) are:
\[ \xi_i=(\alpha_1,\cdots,\alpha_n)X_i \]

Thus, all eigenvectors belonging to \(\lambda\) are:

\[ k_1\xi_1+k_2\xi_2+\cdots+k_t\xi_t \]

where \(k_i\) are not all zero.

The existence of eigenvalues and eigenvectors depends on the field \(V\) is over.

Similarity Transformation

Introduction

If an \(n\times n\) matrix \(A\) is upper triangular, e.g.:

\[ A= \begin{bmatrix} a_{1,1}&a_{1,2}&\cdots &a_{1,n}\\ &a_{2,2}&\cdots &a_{2,n}\\ &&\ddots &\vdots \\ &&&a_{n,n} \end{bmatrix} \]

Then:

\[ \begin{aligned} p_A(x)&=\det(xI_n-A)\\ &= \begin{bmatrix} x-a_{1,1}&-a_{1,2}&\cdots &-a_{1,n}\\ &x-a_{2,2}&\cdots &-a_{2,n}\\ &&\ddots &\vdots \\ &&&x-a_{n,n} \end{bmatrix} \\ &=\prod_{i=1}^n(x-a_{i,i}) \end{aligned} \]

This is easily computed. Lower triangular matrices are similar. However, if \(A\) is neither, we use similarity transformations to convert it into a form where the characteristic polynomial is easy to compute.

Definition

For \(n\times n\) matrices \(A\) and \(B\), if there exists an invertible \(n\times n\) matrix \(P\) such that:

\[ B=P^{-1}AP \]

then matrices \(A\) and \(B\) are similar. The transformation \(A\mapsto P^{-1}AP\) is called a similarity transformation. Moreover, \(A\) and \(P^{-1}AP\) have the same characteristic polynomial.

Consider:

\[ \begin{aligned} \det(xI_n-P^{-1}AP)&=\det(xP^{-1}I_nP-P^{-1}AP)\\ &=\det(P^{-1}xI_nP-P^{-1}AP)\\ &=\det(P^{-1})\cdot \det(P)\cdot \det(xI_n-A)\\ &=\det(xI_n-A)\\ &=p_A(x) \end{aligned} \]

Thus proven. The same holds for \(A\mapsto PAP^{-1}\). Also, \(p_A(0)=(-1)^n\cdot \det(A)\), because \(p_A(0)=\det(-1\cdot I_nA)=\det(-1\cdot I_n)\cdot \det(A)\), so \(\det(A)=\det(P^{-1}AP)\).

Theorem: Similar matrices have the same characteristic polynomial and eigenvalues, but the converse is not true.

This theorem shows that the characteristic polynomial of a linear transformation's matrix is independent of the choice of basis and is determined solely by the linear transformation. Hence, it can be called the characteristic polynomial of the linear transformation.

The characteristic polynomial \(f(\lambda)=|\lambda I-A|\) of matrix \(A\) is a monic polynomial. By Vieta's formulas, its \((n-1)\)-th coefficient is:

\[ -(\lambda_1+\cdots+\lambda_n)=-(a_{11}+\cdots+a_{nn})=-tr A \]

where \(tr A\) is called the trace of \(A\), the sum of its main diagonal elements.

By Vieta's formulas, the constant term of the characteristic polynomial is:

\[ {(-1)}^n|A|={(-1)}^n(\lambda_1\cdots\lambda_n) \]

Theorem: Similar matrices have the same trace.

Commutation Formula

Theorem: Regardless of whether matrices \(A\) and \(B\) are square, as long as multiplication is defined, the trace of \(AB\) equals the trace of \(BA\).

One proof is by direct expansion. Another uses the commutation formula.

Theorem: Let \(A\) be an \(m\times n\) matrix, and \(B\) be an \(n\times m\) matrix. Then:

\[ \lambda^n|\lambda I_m-AB|=\lambda^m|\lambda I_n-BA| \]

This formula shows that \(AB\) and \(BA\) have the same nonzero eigenvalues.

Schur's Lemma

Any \(n\times n\) matrix \(A\) is similar to an upper triangular matrix. That is, there exists a full-rank matrix \(P\) such that \(P^{-1}AP\) is upper triangular, with the main diagonal elements being all eigenvalues of \(A\).

Corollary: Let the \(n\) eigenvalues of \(A\) be \(\lambda_1,\cdots,\lambda_n\), and let \(\phi(x)\) be any polynomial. Then the eigenvalues of the matrix polynomial \(\phi(A)\) are:

\[ \phi(\lambda_1),\cdots,\phi(\lambda_n) \]

In particular, the eigenvalues of \(kA\) are \(k\lambda_1,\cdots,k\lambda_n\), and the eigenvalues of \(A^m\) are \({\lambda_1}^m,\cdots,{\lambda_n}^m\).

Similarity Transformation via Gaussian Elimination

For an \(n\times n\) matrix \(B\), Gaussian elimination can be performed using elementary row operations.

Applying these operations (left-multiplying by elementary matrices) followed by right-multiplying by their inverses constitutes a similarity transformation. Left multiplication corresponds to row operations, and right multiplication corresponds to column operations.

If we can transform the matrix into upper or lower triangular form via similarity transformations, we can easily compute its characteristic polynomial. However, applying the transformation \(A\mapsto T_{ij}(k)AT_{ij}(-k)\) to the main diagonal elements may cause previously zeroed elements to become nonzero again, preventing triangular form.

Later, we will show that applying transformations to the subdiagonal elements still yields a matrix whose characteristic polynomial is easy to compute.

Upper Hessenberg Matrix

For \(n\gt 2\), a matrix of the form:

\[ H= \begin{bmatrix} \alpha_{1}&h_{12}&\dots&\dots&h_{1n}\\ \beta_{2}&\alpha_{2}&h_{23}&\dots &\vdots \\ &\ddots &\ddots & \ddots &\vdots \\ & &\ddots &\ddots & h_{(n-1)n}\\ &&& \beta_{n}& \alpha_{n} \end{bmatrix} \]

is called an upper Hessenberg matrix, where \(\beta\) is the subdiagonal.

Using similarity transformations to eliminate elements below the subdiagonal yields an upper Hessenberg matrix. Computing the characteristic polynomial of an \(n\times n\) upper Hessenberg matrix can be done in \(O(n^3)\) time.

Let \(H_i\) denote the matrix retaining only the first \(i\) rows and columns of \(H\), and let \(p_i(x)=\det(xI_i-H_i)\). Then:

\[ H_0= \begin{bmatrix} \end{bmatrix},\quad p_0(x)=1 \]
\[ H_1= \begin{bmatrix} \alpha_1 \end{bmatrix},\quad p_1(x)=\det(x I_1-H_1)=x -\alpha_1 \]
\[ H_2= \begin{bmatrix} \alpha_1&h_{12}\\ \beta_2&\alpha_2 \end{bmatrix},\quad p_2(x)=\det(xI_2-H_2)=(x-\alpha_2)p_1(x)-\beta_2h_{12}p_0(x) \]

When computing determinants, we typically expand along the row or column with the most zeros. Here, we expand along the last row:

\[ \begin{aligned} p_3(x)&= \det(xI_3-H_3)\\ &=\begin{vmatrix} x-\alpha_1&-h_{12}&-h_{13}\\ -\beta_2&x-\alpha_2&-h_{23}\\ &-\beta_3&x-\alpha_3 \end{vmatrix}\\ &=(x-\alpha_3)\cdot (-1)^{3+3}p_2(x)-\beta_3\cdot (-1)^{3+2} \begin{vmatrix} x-\alpha_1&-h_{13}\\ -\beta_2&-h_{23} \end{vmatrix}\\ &=(x-\alpha_3)p_2(x)-\beta_3(h_{23}p_1(x)+\beta_2h_{13}p_0(x)) \end{aligned} \]

Observing and generalizing, for \(2\leq i\leq n\):

\[ p_i(x)=(x-\alpha_i)p_{i-1}(x)- \sum_{m=1}^{i-1}h_{i-m,i} \left( \prod_{j=i-m+1}^{i}\beta_j \right) p_{i-m-1}(x) \]

This completes the algorithm, commonly known as the Hessenberg algorithm.

Cayley–Hamilton Theorem

For any \(n\times n\) matrix \(A\) with characteristic polynomial \(f(\lambda)=|\lambda I-A|\), it holds that \(f(A)=0\).

For linear transformations \(T\), a parallel result holds: if \(f(\lambda)\) is the characteristic polynomial of \(T\), then \(f(T)\) is the zero transformation.

This theorem implies that for any matrix \(A\), there exists a polynomial that annihilates it.

Minimal Polynomial

Let \(V\) be an \(n\)-dimensional vector space. Since linear transformations correspond to matrices with \(n^2\) elements, all linear transformations form an \(n^2\)-dimensional linear space.

For a specific linear transformation \(T\), from applying it \(0\) times to \(n\) times, there are \(n^2+1\) linear transformations, and their corresponding matrices must be linearly dependent. Thus, there exists a nonzero polynomial \(f\) such that \(f(T)\) is the zero transformation. We say \(T\) satisfies polynomial \(f\). Among all polynomials satisfied by \(T\), there exists one of minimal degree.

The monic polynomial of minimal degree that annihilates matrix \(A\) is called the minimal polynomial of \(A\), denoted \(m_A(\lambda)\).

By the polynomial division algorithm, the minimal polynomial is unique and divides any annihilating polynomial of \(A\). In particular, the minimal polynomial divides the characteristic polynomial.

Theorem: Ignoring multiplicities, the characteristic polynomial \(f(\lambda)\) and the minimal polynomial \(m_A(\lambda)\) of matrix \(A\) have the same roots.

Theorem: Eigenvectors of matrix \(A\) belonging to distinct eigenvalues are linearly independent.

Applications

In informatics, we typically consider matrices over \((\mathbb{Z}/m\mathbb{Z})^{n\times n}\), where \(m\) is often prime. Performing similarity transformations is straightforward in this case. When \(m\) is composite, we can use methods similar to the Euclidean algorithm.

Implementation
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include <cassert>
#include <iostream>
#include <random>
#include <vector>

using Matrix = std::vector<std::vector<int>>;
using i64 = int64_t;

Matrix to_upper_Hessenberg(const Matrix &M, int mod) {
  Matrix H(M);
  int n = H.size();
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      if ((H[i][j] %= mod) < 0) H[i][j] += mod;
    }
  }
  for (int i = 0; i < n - 1; ++i) {
    int pivot = i + 1;
    for (; pivot < n; ++pivot) {
      if (H[pivot][i] != 0) break;
    }
    if (pivot == n) continue;
    if (pivot != i + 1) {
      for (int j = i; j < n; ++j) std::swap(H[i + 1][j], H[pivot][j]);
      for (int j = 0; j < n; ++j) std::swap(H[j][i + 1], H[j][pivot]);
    }
    for (int j = i + 2; j < n; ++j) {
      for (;;) {
        if (H[j][i] == 0) break;
        if (H[i + 1][i] == 0) {
          for (int k = i; k < n; ++k) std::swap(H[i + 1][k], H[j][k]);
          for (int k = 0; k < n; ++k) std::swap(H[k][i + 1], H[k][j]);
          break;
        }
        if (H[j][i] >= H[i + 1][i]) {
          int q = H[j][i] / H[i + 1][i], mq = mod - q;
          for (int k = i; k < n; ++k)
            H[j][k] = (H[j][k] + i64(mq) * H[i + 1][k]) % mod;
          for (int k = 0; k < n; ++k)
            H[k][i + 1] = (H[k][i + 1] + i64(q) * H[k][j]) % mod;
        } else {
          int q = H[i + 1][i] / H[j][i], mq = mod - q;
          for (int k = i; k < n; ++k)
            H[i + 1][k] = (H[i + 1][k] + i64(mq) * H[j][k]) % mod;
          for (int k = 0; k < n; ++k)
            H[k][j] = (H[k][j] + i64(q) * H[k][i + 1]) % mod;
        }
      }
    }
  }
  return H;
}

std::vector<int> get_charpoly(const Matrix &M, int mod) {
  Matrix H(to_upper_Hessenberg(M, mod));
  int n = H.size();
  std::vector<std::vector<int>> p(n + 1);
  p[0] = {1 % mod};
  for (int i = 1; i <= n; ++i) {
    const std::vector<int> &pi_1 = p[i - 1];
    std::vector<int> &pi = p[i];
    pi.resize(i + 1, 0);
    int v = mod - H[i - 1][i - 1];
    if (v == mod) v -= mod;
    for (int j = 0; j < i; ++j) {
      pi[j] = (pi[j] + i64(v) * pi_1[j]) % mod;
      if ((pi[j + 1] += pi_1[j]) >= mod) pi[j + 1] -= mod;
    }
    int t = 1;
    for (int j = 1; j < i; ++j) {
      t = i64(t) * H[i - j][i - j - 1] % mod;
      int prod = i64(t) * H[i - j - 1][i - 1] % mod;
      if (prod == 0) continue;
      prod = mod - prod;
      for (int k = 0; k <= i - j - 1; ++k)
        pi[k] = (pi[k] + i64(prod) * p[i - j - 1][k]) % mod;
    }
  }
  return p[n];
}

bool verify(const Matrix &M, const std::vector<int> &charpoly, int mod) {
  if (mod == 1) return true;
  int n = M.size();
  std::vector<int> randvec(n), sum(n, 0);
  std::mt19937 gen(std::random_device{}());
  std::uniform_int_distribution<int> dis(1, mod - 1);
  for (int i = 0; i < n; ++i) randvec[i] = dis(gen);
  for (int i = 0; i <= n; ++i) {
    int v = charpoly[i];
    for (int j = 0; j < n; ++j) sum[j] = (sum[j] + i64(v) * randvec[j]) % mod;
    std::vector<int> prod(n, 0);
    for (int j = 0; j < n; ++j) {
      for (int k = 0; k < n; ++k) {
        prod[j] = (prod[j] + i64(M[j][k]) * randvec[k]) % mod;
      }
    }
    randvec.swap(prod);
  }
  for (int i = 0; i < n; ++i)
    if (sum[i] != 0) return false;
  return true;
}

int main() {
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int n, mod;
  std::cin >> n >> mod;
  Matrix M(n, std::vector<int>(n));
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j) std::cin >> M[i][j];
  std::vector<int> charpoly(get_charpoly(M, mod));
  for (int i = 0; i <= n; ++i) std::cout << charpoly[i] << ' ';
  assert(verify(M, charpoly, mod));
  return 0;
}

The Hessenberg algorithm above is not numerically stable. For matrices over \(\mathbb{R}^{n\times n}\), other algorithms are needed for adjustment, or alternative numerically stable algorithms should be used.

We can relate the characteristic polynomial to linear recurrences with constant coefficients. It can also be combined with the Cayley–Hamilton theorem and polynomial modular reduction to accelerate algorithms for computing matrix powers over certain fields.

The Cayley–Hamilton theorem states:

\[ \begin{aligned} p_A(A)&=A^n+c_1A^{n-1}+\cdots +c_{n-1}A+c_nI\\ &=O \end{aligned} \]

where \(O\) is the \(n\times n\) zero matrix, \(A\in\mathbb{C}^{n\times n}\), and \(p_A(x)=x^n+\sum_{i=1}^nc_ix^{n-i}\in\mathbb{C}[x]\) is the characteristic polynomial of \(A\).

If we want to compute \(A^K\) for large \(K\), we can compute \(f(x)=x^K\bmod{p_A(x)}\) and then use \(f(A)=A^K\).

Clearly, \(\deg(f(x))\lt n\). Let \(f(x)=\sum_{i=0}^{n-1}f_ix^i\) and \(n=km\). Then:

\[ \begin{aligned} f_{km-1}x^{km-1}+\cdots +f_1x+f_0&=(\cdots (f_{km-1}x^{k-1}+\cdots +f_{k(m-1)})x^k\\ &+f_{k(m-1)-1}x^{k-1}+\cdots +f_{k(m-2)})x^k\\ &+\cdots\\ &+f_{k-1}x^{k-1}+\cdots +f_1x+f_0 \end{aligned} \]

Setting \(k=\sqrt{n}\), we find that computing \(f(A)\) requires approximately \(O(\sqrt{n})\) matrix-matrix multiplications.

References