AIOJ

This part of the discussion considers only square matrices, that is, the linear transformation corresponding to matrix $A$ maps $n$ vectors to $n$ vectors.

Since practical problems often require repeated transformations, merely saying that the linear transformation corresponding to matrix $A$ transforms the identity matrix $I$ into $A$ is too abstract. The best approach is to look for fixed points, that is, parts of the transformation that remain unchanged.

However, the linear transformation corresponding to matrix $A$ may have no fixed points at all. We therefore look instead for collinear parts or parts resembling a simple transformation.

Eigenvalues and Eigenvectors

Under the linear transformation corresponding to matrix $A$ , the directions of some vectors do not change; only their magnitudes are scaled.

Let $V$ be a vector space over $F$ , and let $T$ be a linear transformation on $V$ . If there exist an element of $F$ , denoted by $\lambda$ , and a nonzero vector in $V$ , denoted by $\xi$ , such that:

T\xi=\lambda\xi

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

Eigenvectors on the same line retain their direction under a linear transformation (compressing them to zero is also regarded as preserving direction). An eigenvector is not unique: every vector collinear with an eigenvector is also an eigenvector, but the zero vector is not an eigenvector by definition, since a vector with a direction is necessarily nonzero. The eigenvalue of an eigenvector is the factor by which it is scaled.

In applications, for eigenvectors with the same eigenvalue, one usually chooses a basis as a representative of the whole set.

Let $\alpha_1,\alpha_2,\cdots,\alpha_n$ be a basis of $V$ , and let the matrix of $T$ in this basis be $A$ , that is:

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

Let $\lambda_0$ be an eigenvalue of $T$ , let $\xi$ be an eigenvector of $T$ belonging to the eigenvalue $\lambda_0$ , and let a nonzero vector $X$ satisfy:

\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

Therefore, 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 is an $n$ -th degree polynomial whose roots are the eigenvalues of $A$ , and is denoted by $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}

Here $I_n$ is an $n\times n$ identity matrix. Some sources define it as $p_A(\lambda)=\det(A-\lambda I_n)$ , which differs from our definition only by the sign $(-1)^n$ . With our definition, $p_A(\lambda)$ is always monic, whereas the other definition 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 of $(\lambda_0 I-A)X=0$ , denoted by $X$ , 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

By 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}

The $d_i$ are called the algebraic multiplicity of the eigenvalue $\lambda_i$ . The sum of all algebraic multiplicities is the dimension $n$ of the space.

Finding All Eigenvalues and Eigenvectors of a Matrix

The process consists of the following steps:

Compute the determinant $|\lambda I-A|$ .
Find all roots of the polynomial $f(\lambda)=|\lambda I-A|$ in the field $F$ ; these are the eigenvalues of $A$ .
For each eigenvalue of $A$ , denoted by $\lambda$ , solve the homogeneous linear system $(\lambda I-A)X=0$ and find a fundamental solution system $X_1,\cdots,X_t$ . Then all eigenvectors of $A$ belonging to $\lambda$ are:

k_1X_1+k_2X_2+\cdots+k_tX_t

The $k_i$ in this expression are not all zero.

The eigenvectors of the linear transformation $T$ belonging to $\lambda$ are:

\xi_i=(\alpha_1,\cdots,\alpha_n)X_i

Therefore, all eigenvectors belonging to $\lambda$ are:

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

The $k_i$ in this expression are not all zero.

Whether eigenvalues and eigenvectors exist depends on the field containing $V$ .

Similarity Transformation

Introduction

If an $n\times n$ matrix $A$ is an upper-triangular matrix of the form

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}

The characteristic polynomial of a lower-triangular matrix can be obtained similarly. If $A$ is neither of these two types, a similarity transformation is needed to put the matrix into a form whose characteristic polynomial is easier to compute.

Definition

For $n\times n$ matrices $A$ and $B$ , if there exists an invertible $n\times n$ matrix $P$ satisfying

B=P^{-1}AP

then matrices $A$ and $B$ are similar, and the transformation $A\mapsto P^{-1}AP$ is called a similarity transformation. The matrices $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}

This proves the claim; the case of $A\mapsto PAP^{-1}$ is the same. 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.

The theorem shows that the characteristic polynomial of the matrix of a linear transformation is independent of the choice of basis and is determined directly by the linear transformation. It can therefore be called the characteristic polynomial of the linear transformation.

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

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

Here $tr A$ is called the trace of $A$ and is the sum of the entries on the main diagonal of $A$ .

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.

Sylvester's Determinant Identity

Theorem: Whether or not matrices $A$ and $B$ are square, as long as their product is defined, the trace of $AB$ equals the trace of $BA$ .

One proof is by direct expansion. Another proof uses Sylvester's determinant identity.

Theorem: Let $A$ be an $m$ -by- $n$ matrix and $B$ an $n$ -by- $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

Every $n$ -th order matrix $A$ is similar to an upper-triangular matrix. In other words, there exists a full-rank matrix $P$ such that $P^{-1}AP$ is upper triangular, with all eigenvalues of $A$ on its main diagonal.

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

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

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

Similarity Transformations Using Gaussian Elimination

Gaussian elimination can be performed on an $n\times n$ matrix $B$ using elementary row operations.

After applying the above operations to a matrix (left-multiplying by elementary matrices), right-multiplying by the inverse elementary matrix gives a similarity transformation. The left multiplication is a row operation, and the right multiplication is correspondingly a column operation.

If a matrix can be transformed into upper- or lower-triangular form by similarity transformations, its characteristic polynomial can be obtained easily. However, applying $A\mapsto T_{ij}(k)AT_{ij}(-k)$ to an element on the main diagonal may cause the following operation to restore a previously eliminated zero: first apply $A\mapsto T_{ij}(k)A$ to eliminate the element in row $i$ , column $j$ , then right-multiply by $T_{ij}(-k)$ , which adds the $A$ -column indexed by $i$ multiplied by $-k$ to column $j$ . Thus the matrix may fail to become upper or lower triangular.

The following section shows that applying transformations to elements on the subdiagonal still produces a matrix whose characteristic polynomial can be computed easily.

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.

We can obtain an upper Hessenberg matrix by using similarity transformations to eliminate all elements below the subdiagonal. The characteristic polynomial of an $n\times n$ upper Hessenberg matrix can then be computed in $O(n^3)$ time.

Let $H_i$ be the matrix obtained from $H$ by retaining only its first $i$ rows and first $i$ columns, 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 a determinant, we usually expand along a row or column with as many zeros as possible. The corresponding minor is the matrix obtained by deleting the row and column containing the selected element. 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 the pattern and generalizing, for $2\leq i\leq n$ we have

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, which is generally called the Hessenberg algorithm.

Cayley–Hamilton Theorem

For any $n$ -th order matrix $A$ with characteristic polynomial $f(\lambda)=|\lambda I-A|$ , we have $f(A)=0$ .

There is a corresponding result for a linear transformation $T$ : if $f(\lambda)$ is the characteristic polynomial of $T$ , then $f(T)$ is the zero transformation.

This theorem implies that every matrix $A$ has a polynomial that annihilates it.

Minimal Polynomial

Let $V$ be an $n$ -dimensional vector space. Since the matrix corresponding to a linear transformation has $n^2$ entries, all linear transformations form an $n^2$ -dimensional vector space.

For a particular linear transformation $T$ , the transformations obtained by applying it from $0$ through $n$ times give $n^2+1$ linear transformations in total, and their corresponding matrices must be linearly dependent. Therefore, there is a nonzero polynomial $f$ such that $f(T)$ is the zero transformation; we say that transformation $T$ satisfies polynomial $f$ . Among all polynomials satisfied by $T$ , one denoted by $f$ has minimum degree.

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

By the Euclidean algorithm for polynomials, the minimal polynomial is unique and divides every polynomial that annihilates $A$ . In particular, the minimal polynomial divides the characteristic polynomial.

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

Theorem: Eigenvectors of matrix $A$ belonging to distinct eigenvalues are linearly independent.

Applications

In informatics, we usually consider matrices over $(\mathbb{Z}/m\mathbb{Z})^{n\times n}$ , with $m$ usually prime, so the similarity transformations above are straightforward. When $m$ is composite, we can use an approach similar to the Euclidean algorithm.

Implementation

#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, so matrices over $\mathbb{R}^{n\times n}$ need to be preprocessed with another algorithm or handled with a different numerically stable algorithm.

The characteristic polynomial can be connected with constant-coefficient homogeneous linear recurrences. Together with the Cayley–Hamilton theorem and polynomial modular reduction, it can also accelerate algorithms for computing matrix powers over some 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}

Here $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 need $A^K$ for a large $K$ , we can first 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}

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

References

Rizwana Rehman, Ilse C.F. Ipsen.La Budde’s Method for Computing Characteristic Polynomials.
Marshall Law.Computing Characteristic Polynomials of Matrices of Structured Polynomials.
Mike Paterson.On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials.