跳到正文
LESSON

7.3.13 特征多项式

我们考虑一个 n × n n\times n n × n 的矩阵 A A A ，其中 n ≥ 0 ∧ n ∈ Z n\geq 0\land n\in\mathbb{Z} n ≥ 0 ∧ n ∈ Z 。它的特征多项式记为 p A ( x ) p_A(x) p A ​ ( x ) 其中

阅读正文 →

特征多项式的定义

我们考虑一个 n×nn\times n 的矩阵 AA，其中 n≥0∧n∈Zn\geq 0\land n\in\mathbb{Z}。它的特征多项式记为 pA(x)p_A(x) 其中

pA(x)=det⁡(xIn−A) p_A(x)=\det(xI_n-A)

其中 InI_n 为一个 n×nn\times n 的单位矩阵。一些地方会定义为 pA(x)=det⁡(A−xIn)p_A(x)=\det(A-xI_n) 与我们的定义仅相差了一个符号 (−1)n(-1)^n，但我们采用这种定义得到的 pA(x)p_A(x) 一定为首一多项式，而另外的定义则仅当 nn 为偶数时才是首一多项式。需要注意的是 0×00\times 0 的矩阵行列式为 11 是良定义的。

特征多项式的求法

若 n×nn\times n 的矩阵 AA 为上三角矩阵如

A=[a1,1a1,2⋯a1,na2,2⋯a2,n⋱⋮an,n] 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}

那么

pA(x)=det⁡(xIn−A)=[x−a1,1−a1,2⋯−a1,nx−a2,2⋯−a2,n⋱⋮x−an,n]=∏i=1n(x−ai,i) \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}

可轻松求得，下三角矩阵也是类似的。但如果 AA 不属于这两种矩阵，我们需要使用相似变换来使得矩阵变为容易求得特征多项式的形式。

相似变换

对于 n×nn\times n 的矩阵 AA 和 BB，当存在 n×nn\times n 的可逆矩阵 PP 满足

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

我们说 AA 和 BB 相似，记变换 A↦P−1APA\mapsto P^{-1}AP 为相似变换。且我们说 AA 和 P−1APP^{-1}AP 有相同的特征多项式。

考虑

det⁡(xIn−P−1AP)=det⁡(xP−1InP−P−1AP)=det⁡(P−1xInP−P−1AP)=det⁡(P−1)⋅det⁡(P)⋅det⁡(xIn−A)=det⁡(xIn−A)=pA(x) \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}

得证，且我们发现 A↦PAP−1A\mapsto PAP^{-1} 也是一样的。另外 pA(0)=(−1)n⋅det⁡(A)p_A(0)=(-1)^n\cdot \det(A)，因为 pA(0)=det⁡(−1⋅InA)=det⁡(−1⋅In)⋅det⁡(A)p_A(0)=\det(-1\cdot I_nA)=\det(-1\cdot I_n)\cdot \det(A) 故 det⁡(A)=det⁡(P−1AP)\det(A)=\det(P^{-1}AP)。

使用高斯消元进行相似变换

对 n×nn\times n 的矩阵 BB 进行高斯消元行变换的三个基本操作有

  • 将 BB 的第 ii，jj 行互换：B↦PijBB\mapsto P_{ij}B。
  • 将 BB 的第 ii 行乘 k≠0k\neq 0：B↦Di(k)BB\mapsto D_i(k)B。
  • 将 BB 的第 jj 行的 kk 倍加到第 ii 行：B↦Tij(k)BB\mapsto T_{ij}(k)B。

对于 1≤i<j≤n1\leq i\lt j\leq n 这几个 n×nn\times n 的初等矩阵分别为

Pij=[Ii−101Ij−i−110In−j] P_{ij}= \begin{bmatrix} I_{i-1}&&&&\\ &0&&1&\\ &&I_{j-i-1}&&\\ &1&&0&\\ &&&&I_{n-j} \end{bmatrix}

其中两个 11 分别在 PijP_{ij} 的第 ii 和 jj 行，注意 Pij−1=PijP_{ij}^{-1}=P_{ij}。

Di(k)=[Ii−1kIn−i] D_i(k)= \begin{bmatrix} I_{i-1}&&\\ &k&\\ &&I_{n-i} \end{bmatrix}

其中 kk 在 Di(k)D_i(k) 的第 ii 行第 ii 列，注意 Di(k)−1=Di(k−1)D_i(k)^{-1}=D_i(k^{-1})。

Tij(k)=[Ii−11k⋱1In−j] T_{ij}(k)= \begin{bmatrix} I_{i-1}&&&&\\ &1&&k&\\ &&\ddots &&\\ &&&1&\\ &&&&I_{n-j} \end{bmatrix}

其中 kk 在 Tij(k)T_{ij}(k) 的第 ii 行第 jj 列，注意 Tij(k)−1=Tij(−k)T_{ij}(k)^{-1}=T_{ij}(-k)。

若我们记 EijE_{ij} 为第 ii 行第 jj 列的元素为 11、其余为零的 n×nn\times n 矩阵，那么

  • Pij=In−Eii−Ejj+Eij+EjiP_{ij}=I_n-E_{ii}-E_{jj}+E_{ij}+E_{ji}
  • Di(k)=In+(k−1)EiiD_i(k)=I_n+(k-1)E_{ii}
  • Tij(k)=In+kEijT_{ij}(k)=I_n+kE_{ij}

易验证其逆矩阵。

我们在对矩阵使用上述操作（左乘初等矩阵）后再右乘其逆矩阵即相似变换，左乘为行变换，易发现右乘即列变换。

若我们能将矩阵通过相似变换变为上三角或下三角的形式，那么可以轻松求出其特征多项式。但我们发现若对主对角线上的元素应用变换 A↦Tij(k)ATij(−k)A\mapsto T_{ij}(k)AT_{ij}(-k) 后会导致原本通过 A↦Tij(k)AA\mapsto T_{ij}(k)A 将第 ii 行第 jj 列的元素消为零后右乘 Tij(−k)T_{ij}(-k) 即将 AA 的第 ii 列的 −k-k 倍加到第 jj 列这一操作使得之前消为零的元素现在可能不为零，可能不能将其变为上三角或下三角形式。

后文将说明我们对次对角线上的元素应用变换后得到的矩阵依然可以轻松得到其特征多项式。

上 Hessenberg 矩阵

对于 n>2n\gt 2 的形如

H=[α1h12……h1nβ2α2h23…⋮⋱⋱⋱⋮⋱⋱h(n−1)nβnαn] 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}

的矩阵我们称为上 Hessenberg 矩阵，其中 β\beta 为次对角线。

我们使用相似变换将次对角线以下的元素消为零后即能得到上 Hessenberg 矩阵，而求出一个 n×nn\times n 上 Hessenberg 矩阵的特征多项式则可在 O(n3)O(n^3) 时间完成。

我们记 HiH_i 为只保留 HH 的前 ii 行和前 ii 列的矩阵，记 pi(x)=det⁡(xIi−Hi)p_i(x)=\det(xI_i-H_i) 那么

H0=[],p0(x)=1 H_0= \begin{bmatrix} \end{bmatrix},\quad p_0(x)=1

H1=[α1],p1(x)=det⁡(xI1−H1)=x−α1 H_1= \begin{bmatrix} \alpha_1 \end{bmatrix},\quad p_1(x)=\det(x I_1-H_1)=x -\alpha_1

H2=[α1h12β2α2],p2(x)=det⁡(xI2−H2)=(x−α2)p1(x)−β2h12p0(x) 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)

在计算行列式时我们一般选择按零最多的行或列余子式展开，余子式即删除了当前选择的元素所在行和列之后的矩阵，在这里我们选择按最后一行进行展开，有

p3(x)=det⁡(xI3−H3)=∣x−α1−h12−h13−β2x−α2−h23−β3x−α3∣=(x−α3)⋅(−1)3+3p2(x)−β3⋅(−1)3+2∣x−α1−h13−β2−h23∣=(x−α3)p2(x)−β3(h23p1(x)+β2h13p0(x)) \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}

观察并归纳，对 2≤i≤n2\leq i\leq n 有

pi(x)=(x−αi)pi−1(x)−∑m=1i−1hi−m,i(∏j=i−m+1iβj)pi−m−1(x) 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)

至此完成了整个算法，该算法一般被称为 Hessenberg 算法。

应用

在信息学中我们一般考虑 (Z/mZ)n×n(\mathbb{Z}/m\mathbb{Z})^{n\times n} 上的矩阵，通常 mm 为素数，进行上述相似变换是简单的，当 mm 为合数时，我们可以考虑类似辗转相除的方法来进行。

 #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;
    }

上述 Hessenberg 算法不具有数值的稳定性，所以 Rn×n\mathbb{R}^{n\times n} 上的矩阵在使用前需要其他算法进行调整或改用其他具有数值稳定性的算法。

我们可以将特征多项式与常系数齐次线性递推联系起来，也可结合 Cayley-Hamilton 定理、多项式取模加速一些域上求矩阵幂次的算法。

Cayley-Hamilton 定理指出

pA(A)=An+c1An−1+⋯+cn−1A+cnI=O \begin{aligned} p_A(A)&=A^n+c_1A^{n-1}+\cdots +c_{n-1}A+c_nI\\ &=O \end{aligned}

其中 OO 为 n×nn\times n 的零矩阵，A∈Cn×nA\in\mathbb{C}^{n\times n} 且 pA(x)=xn+∑i=1ncixn−i∈C[x]p_A(x)=x^n+\sum_{i=1}^nc_ix^{n-i}\in\mathbb{C}[x] 为 AA 的特征多项式。

若我们要求 AKA^K 其中 KK 较大，那么可以求出 f(x)=xK mod pA(x)f(x)=x^K\bmod{p_A(x)} 后利用 f(A)=AKf(A)=A^K。

而 deg⁡(f(x))<n\deg(f(x))\lt n 显然。我们令 f(x)=∑i=0n−1fixif(x)=\sum_{i=0}^{n-1}f_ix^i 且 n=kmn=km 那么

fkm−1xkm−1+⋯+f1x+f0=(⋯(fkm−1xk−1+⋯+fk(m−1))xk+fk(m−1)−1xk−1+⋯+fk(m−2))xk+⋯+fk−1xk−1+⋯+f1x+f0 \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}

令 k=nk=\sqrt{n} 可以发现计算 f(A)f(A) 大约需要 O(n)O(\sqrt{n}) 次矩阵与矩阵的乘法。

上一篇牛顿迭代法