7.3.13 Characteristic polynomial

# Definition of characteristic polynomial

Let's consider a $n\times n$ matrix $A$, among which $n\geq 0\land n\in\mathbb{Z}$. Its characteristic polynomial is written as $p_A(x)$ Among them

$$p_A(x)=\det(xI_n-A)$$

Among them $I_n$ for one $n\times n$ the identity matrix. In some places it is defined as $p_A(x)=\det(A-xI_n)$ differs by only one symbol from our definition $(-1)^n$, but using this definition we get $p_A(x)$ must be a leading polynomial, and other definitions are only if $n$ It is a leading polynomial only when it is an even number. It should be noted that $0\times 0$ The matrix determinant of is $1$ It is well defined.

# How to find characteristic polynomials

If $n\times n$ matrix $A$ is an upper triangular matrix such as

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

It can be easily obtained, and the lower triangular matrix is similar. But if $A$ Does not belong to these two matrices, we need to use similarity transformation to make the matrix into a form that is easy to find the characteristic polynomial.

# similarity transformation

for $n\times n$ matrix $A$ and $B$, when exists $n\times n$ the invertible matrix $P$ satisfy

$$B=P^{-1}AP$$

we say $A$ and $B$ Similar, note transformation $A\mapsto P^{-1}AP$ is a similar transformation. Let us say $A$ and $P^{-1}AP$ have the same characteristic polynomials.

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

proved, and we found $A\mapsto PAP^{-1}$ Same thing. In addition $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)$ Therefore $\det(A)=\det(P^{-1}AP)$。

# Similarity transformation using Gaussian elimination

Yes $n\times n$ matrix $B$ The three basic operations for performing Gaussian elimination row transformation are

• will $B$ of the $i$，$j$ Row swap:$B\mapsto P_{ij}B$。
• will $B$ of the $i$ travel by $k\neq 0$：$B\mapsto D_i(k)B$。
• will $B$ of the $j$ OK $k$ doubled to $i$ Line:$B\mapsto T_{ij}(k)B$。

for $1\leq i\lt j\leq n$ These few $n\times n$ The elementary matrices of are respectively

$$P_{ij}= \begin{bmatrix} I_{i-1}&&&&\\ &0&&1&\\ &&I_{j-i-1}&&\\ &1&&0&\\ &&&&I_{n-j} \end{bmatrix}$$

two of them $1$ respectively in $P_{ij}$ of the $i$ and $j$ OK, pay attention $P_{ij}^{-1}=P_{ij}$。

$$D_i(k)= \begin{bmatrix} I_{i-1}&&\\ &k&\\ &&I_{n-i} \end{bmatrix}$$

Among them $k$ in $D_i(k)$ of the $i$ line of work $i$ column, note $D_i(k)^{-1}=D_i(k^{-1})$。

$$T_{ij}(k)= \begin{bmatrix} I_{i-1}&&&&\\ &1&&k&\\ &&\ddots &&\\ &&&1&\\ &&&&I_{n-j} \end{bmatrix}$$

Among them $k$ in $T_{ij}(k)$ of the $i$ line of work $j$ column, note $T_{ij}(k)^{-1}=T_{ij}(-k)$。

If we remember $E_{ij}$ for the first $i$ line of work $j$ The elements of the column are $1$, the rest is zero $n\times n$ matrix, then

• $P_{ij}=I_n-E_{ii}-E_{jj}+E_{ij}+E_{ji}$
• $D_i(k)=I_n+(k-1)E_{ii}$
• $T_{ij}(k)=I_n+kE_{ij}$

It is easy to verify its inverse matrix.

After we use the above operation on the matrix (left multiplication of the elementary matrix), we then right-multiply its inverse matrix, which is the similarity transformation. The left multiplication is the row transformation, and it is easy to find that the right multiplication is the column transformation.

If we can transform the matrix into an upper triangular or lower triangular form through similarity transformation, then its characteristic polynomial can be easily obtained. But we found that if we apply a transformation to the elements on the main diagonal $A\mapsto T_{ij}(k)AT_{ij}(-k)$ will later cause the original pass $A\mapsto T_{ij}(k)A$ General $i$ line of work $j$ The elements of the column are eliminated to zero and then multiplied by the right $T_{ij}(-k)$ Coming soon $A$ of the $i$ column $-k$ doubled to $j$ The operation of column makes it possible that elements that were previously eliminated to zero may now be non-zero, and may not be transformed into upper or lower triangular form.

It will be explained later that the matrix obtained after applying transformation to the elements on the subdiagonal can still easily obtain its characteristic polynomial.

# Upper Hessenberg matrix

for $n\gt 2$ shaped like

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

The matrix is called the upper Hessenberg matrix, where $\beta$ is a sub-diagonal line.

We can obtain the upper Hessenberg matrix by using a similarity transformation to eliminate the elements below the sub-diagonal to zero, and find a $n\times n$ The characteristic polynomial of the upper Hessenberg matrix can be found in $O(n^3)$ Time is done.

we remember $H_i$ reserved for only $H$ before $i$ line and before $i$ column matrix, note $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 calculating the determinant, we generally choose to expand by the cofactor of the row or column with the most zeros. The cosubform is the matrix after deleting the row and column of the currently selected element. Here we choose to expand by the last row, and there is

$$\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}$$

Observe and summarize, yes $2\leq i\leq n$ Yes

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

# Application

In informatics we generally consider $(\mathbb{Z}/m\mathbb{Z})^{n\times n}$ matrix, usually $m$ is a prime number, it is simple to perform the above similarity transformation, when $m$ When it is a composite number, we can consider a method similar to euclidean division.

 #include <iostream> #include <random> #include <vector> typedef std::vector<std::vector<int>> Matrix; typedef long long i64; 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(0); 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 above Hessenberg algorithm is not numerically stable, so $\mathbb{R}^{n\times n}$ The matrix above needs to be adjusted by other algorithms or switched to other algorithms with numerical stability before use.

We can connect characteristic polynomials with homogeneous linear recursion with constant coefficients, or combine it with the Cayley-Hamilton theorem and polynomial modulo to speed up algorithms for finding matrix powers in some domains.

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

Among them $O$ for $n\times n$ The zero matrix of$A\in\mathbb{C}^{n\times n}$ And $p_A(x)=x^n+\sum_{i=1}^nc_ix^{n-i}\in\mathbb{C}[x]$ for $A$ The characteristic polynomial of .

if we ask $A^K$ Among them $K$ is larger, then we can find $f(x)=x^K\bmod{p_A(x)}$ Later use $f(A)=A^K$。

And $\deg(f(x))\lt n$ Obviously. we order $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}$$

Order $k=\sqrt{n}$ Calculations can be found $f(A)$ It takes about $O(\sqrt{n})$ Submatrix and matrix multiplication.