跳转至

特徵多項式

特徵的這部分只研究方陣，即矩陣 \(A\) 對應的線性變換將 \(n\) 個向量映射到 \(n\) 個向量。

由於在實際問題中，經常要考慮連續進行重複的變換，如果只用「矩陣 \(A\) 對應的線性變換將單位陣 \(I\) 變換為 \(A\)」的描述，就會很抽象。此時最好的辦法是找「不動點」，即變換當中不動的部分。

然而事實上，矩陣 \(A\) 對應的線性變換很可能沒有不動點，於是退而求其次，尋找共線或者類似於簡單變形的部分。

特徵值與特徵向量

在矩陣 \(A\) 對應的線性變換作用下，一些向量的方向不改變，只是伸縮了。

設 \(V\) 是 \(F\) 上的線性空間，\(T\) 是 \(V\) 上的線性變換。若存在 \(F\) 中的 \(\lambda\) 與 \(V\) 中的 非零向量 \(\xi\)，使得：

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

則稱 \(\lambda\) 為 \(T\) 的一個 特徵值，而 \(\xi\) 為 \(T\) 的 屬於特徵值 \(\lambda\) 的一個特徵向量。

特徵向量在同一直線上，在線性變換作用下保持方向不改變（壓縮到零也認為是方向不改變）。特徵向量不唯一，與特徵向量共線的向量都是特徵向量，但是規定零向量不是特徵向量，擁有方向的向量自然是非零向量。特徵向量的特徵值就是它伸縮的倍數。

在實際應用中，一般對於擁有相同特徵值的特徵向量，會選取一組基作為它們全體的代表。

設 \(\alpha_1,\alpha_2,\cdots,\alpha_n\) 是 \(V\) 的一組基，\(T\) 在這組基下的矩陣為 \(A\)，即：

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

設 \(\lambda_0\) 是 \(T\) 的一個特徵值，\(\xi\) 為 \(T\) 的屬於特徵值 \(\lambda_0\) 的一個特徵向量，且有非零向量 \(X\) 滿足：

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

於是有：

\[ 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 \]

所以相應的行列式也為 \(0\)。

特徵多項式

考慮一個 \(n\times n\) 的矩陣 \(A\)，其中 \(n\geq 0\land n\in\mathbb{Z}\)。設 \(\lambda\) 為一個參量，矩陣 \(\lambda I-A\) 稱為 \(A\) 的 特徵矩陣。

特徵矩陣的行列式稱為 \(A\) 的 特徵多項式，展開為一個 \(n\) 次多項式，根為 \(A\) 的特徵值，記為 \(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} \]

其中 \(I_n\) 為一個 \(n\times n\) 的單位矩陣。一些地方會定義為 \(p_A(\lambda)=\det(A-\lambda I_n)\) 與我們的定義僅相差了一個符號 \((-1)^n\)，但採用這種定義得到的 \(p_A(\lambda)\) 一定為首一多項式，而另外的定義則僅當 \(n\) 為偶數時才是首一多項式。需要注意的是 \(0\times 0\) 的矩陣行列式為 \(1\) 是良定義的。

相應於 \((\lambda_0 I-A)X=0\) 的非零解向量 \(X\)，稱為 \(A\) 的屬於 \(\lambda_0\) 的特徵向量。

線性變換 \(T\) 有特徵值 \(\lambda_0\) 等價於矩陣 \(A\) 有特徵值 \(\lambda_0\)。

線性變換 \(T\) 有特徵向量 \(\xi\) 等價於矩陣 \(A\) 有特徵向量 \(X\)，其中有：

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

根據代數基本定理，特徵多項式可以分解為：

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

稱 \(d_i\) 為特徵值 \(\lambda_i\) 的 代數重數。全體代數重數的和為空間維數 \(n\)。

求解矩陣的全部特徵值及特徵向量

分為以下步驟：

  • 計算行列式 \(|\lambda I-A|\)。
  • 求出多項式 \(f(\lambda)=|\lambda I-A|\) 在域 \(F\) 中的全部根，即 \(A\) 的特徵值。
  • 對 \(A\) 的每個特徵值 \(\lambda\)，解齊次線性方程組 \((\lambda I-A)X=0\)，求出它的一組基礎解系 \(X_1,\cdots,X_t\)，則 \(A\) 的屬於 \(\lambda\) 的全部特徵向量為：
\[ k_1X_1+k_2X_2+\cdots+k_tX_t \]

該表達式中的 \(k_i\) 不全為零。

  • 線性變換 \(T\) 的屬於 \(\lambda\) 的特徵向量為：
\[ \xi_i=(\alpha_1,\cdots,\alpha_n)X_i \]

因此，屬於 \(\lambda\) 的全部特徵向量為：

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

該表達式中的 \(k_i\) 不全為零。

特徵值與特徵向量是否存在，依賴於 \(V\) 所在的域。

相似變換

引入

若 \(n\times n\) 的矩陣 \(A\) 為上三角矩陣如

\[ 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} \]

那麼

\[ \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} \]

可輕鬆求得，下三角矩陣也是類似的。但如果 \(A\) 不屬於這兩種矩陣，則需要使用相似變換，使得矩陣變為容易求得特徵多項式的形式。

定義

對於 \(n\times n\) 的矩陣 \(A\) 和 \(B\)，當存在 \(n\times n\) 的可逆矩陣 \(P\) 滿足

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

則矩陣 \(A\) 和 \(B\) 相似，記變換 \(A\mapsto P^{-1}AP\) 為相似變換。且 \(A\) 和 \(P^{-1}AP\) 有相同的特徵多項式。

考慮

\[ \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\mapsto PAP^{-1}\) 也是一樣的。另外 \(p_A(0)=(-1)^n\cdot \det(A)\)，因為 \(p_A(0)=\det(-1\cdot I_nA)=\det(-1\cdot I_n)\cdot \det(A)\) 故 \(\det(A)=\det(P^{-1}AP)\)。

定理：相似矩陣有相同的特徵多項式及特徵值，反之不然。

定理表明，線性變換的矩陣的特徵多項式與基的選取無關，而直接由線性變換決定，故可稱之為線性變換的特徵多項式。

矩陣 \(A\) 的特徵多項式 \(f(\lambda)=|\lambda I-A|\) 是一個首一的多項式。根據韋達定理，它的 \(n-1\) 次係數為：

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

其中 \(tr A\) 稱為 \(A\) 的跡，為 \(A\) 的主對角線元素之和。

根據韋達定理，特徵多項式的常數項為：

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

定理：相似的矩陣有相同的跡。

換位公式

定理：無論矩陣 \(A\) 和矩陣 \(B\) 是否方陣，只要乘法能進行，則矩陣 \(AB\) 的跡等於矩陣 \(BA\) 的跡。

一種證法是直接展開，即證畢。另一種證法用到換位公式。

定理：設 \(A\) 為 \(m\) 行 \(n\) 列矩陣，設 \(B\) 為 \(n\) 行 \(m\) 列矩陣，則有：

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

該公式表明 \(AB\) 與 \(BA\) 有相同的非零特徵值。

舒爾（Schur）引理

任意的 \(n\) 階矩陣 \(A\) 都相似於一個上三角陣，即存在滿秩陣 \(P\)，使得 \(P^{-1}AP\) 為上三角陣，它的主對角線上元素為 \(A\) 的全部特徵值。

推論：設 \(A\) 的 \(n\) 個特徵值為 \(\lambda_1,\cdots,\lambda_n\)，\(\phi(x)\) 為任一多項式，則矩陣多項式 \(\phi(A)\) 的 \(n\) 個特徵值為：

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

特別地，\(kA\) 的特徵值為 \(k\lambda_1,\cdots,k\lambda_n\)，\(A^m\) 的特徵值為 \({\lambda_1}^m,\cdots,{\lambda_n}^m\)。

使用高斯消元進行相似變換

對 \(n\times n\) 的矩陣 \(B\) 可以進行高斯消元，其基本操作為初等行變換。

在對矩陣使用上述操作（左乘初等矩陣）後再右乘其逆矩陣即相似變換，左乘為行變換，易發現右乘即列變換。

若能將矩陣通過相似變換變為上三角或下三角的形式，那麼可以輕鬆求出其特徵多項式。但若對主對角線上的元素應用變換 \(A\mapsto T_{ij}(k)AT_{ij}(-k)\) 後會導致原本通過 \(A\mapsto T_{ij}(k)A\) 將第 \(i\) 行第 \(j\) 列的元素消為零後右乘 \(T_{ij}(-k)\) 即將 \(A\) 的第 \(i\) 列的 \(-k\) 倍加到第 \(j\) 列這一操作使得之前消為零的元素現在可能不為零，可能不能將其變為上三角或下三角形式。

後文將説明對次對角線上的元素應用變換後得到的矩陣依然可以輕鬆得到其特徵多項式。

上 Hessenberg 矩陣

對於 \(n\gt 2\) 的形如

\[ 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\times n\) 上 Hessenberg 矩陣的特徵多項式則可在 \(O(n^3)\) 時間完成。

我們記 \(H_i\) 為只保留 \(H\) 的前 \(i\) 行和前 \(i\) 列的矩陣，記 \(p_i(x)=\det(xI_i-H_i)\) 那麼

\[ 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) \]

在計算行列式時我們一般選擇按零最多的行或列餘子式展開，餘子式即刪除了當前選擇的元素所在行和列之後的矩陣，在這裏我們選擇按最後一行進行展開，有

\[ \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\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) \]

至此完成了整個算法，該算法一般被稱為 Hessenberg 算法。

Cayley–Hamilton 定理

對於任意的 \(n\) 階矩陣 \(A\)，特徵多項式為 \(f(\lambda)=|\lambda I-A|\)，則必有 \(f(A)=0\)。

對於線性變換 \(T\) 有平行的結果：如果 \(f(\lambda)\) 為 \(T\) 的特徵多項式，則 \(f(T)\) 為零變換。

由本定理可知，對於任意的矩陣 \(A\)，必有可以使其零化的多項式。

最小多項式

設 \(V\) 是一個 \(n\) 維向量空間，由於線性變換對應的矩陣有 \(n^2\) 個元素，一切線性變換構成 \(n^2\) 維線性空間。

對於一個特定的線性變換 \(T\)，從作用 \(0\) 次到作用 \(n\) 次，總共 \(n^2+1\) 個線性變換，它們對應的矩陣一定線性相關。於是存在非零多項式 \(f\)，使得 \(f(T)\) 為零變換，稱變換 \(T\) 滿足多項式 \(f\)。在 \(T\) 滿足的所有多項式 \(f\) 中，存在次數最低的。

可以將矩陣 \(A\) 零化的最小次數的首一多項式稱為 \(A\) 的最小多項式，記為 \(m_A(\lambda)\)。

根據多項式的輾轉相除法，最小多項式是唯一的，且可整除任一 \(A\) 的零化多項式。特別地，最小多項式整除特徵多項式。

定理：在不計重數的情況下，矩陣 \(A\) 的特徵多項式 \(f(\lambda)\) 與最小多項式 \(m_A(\lambda)\) 有相同的根。

定理：矩陣 \(A\) 的屬於不同特徵值的特徵向量線性無關。

應用

在信息學中我們一般考慮 \((\mathbb{Z}/m\mathbb{Z})^{n\times n}\) 上的矩陣，通常 \(m\) 為素數，進行上述相似變換是簡單的，當 \(m\) 為合數時，我們可以考慮類似輾轉相除的方法來進行。

實現 
  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>

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 算法不具有數值的穩定性，所以 \(\mathbb{R}^{n\times n}\) 上的矩陣在使用前需要其他算法進行調整或改用其他具有數值穩定性的算法。

我們可以將特徵多項式與常係數齊次線性遞推聯繫起來，也可結合 Cayley–Hamilton 定理、多項式取模加速一些域上求矩陣冪次的算法。

Cayley–Hamilton 定理指出

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

其中 \(O\) 為 \(n\times n\) 的零矩陣，\(A\in\mathbb{C}^{n\times n}\) 且 \(p_A(x)=x^n+\sum_{i=1}^nc_ix^{n-i}\in\mathbb{C}[x]\) 為 \(A\) 的特徵多項式。

若我們要求 \(A^K\) 其中 \(K\) 較大，那麼可以求出 \(f(x)=x^K\bmod{p_A(x)}\) 後利用 \(f(A)=A^K\)。

而 \(\deg(f(x))\lt n\) 顯然。我們令 \(f(x)=\sum_{i=0}^{n-1}f_ix^i\) 且 \(n=km\) 那麼

\[ \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=\sqrt{n}\) 可以發現計算 \(f(A)\) 大約需要 \(O(\sqrt{n})\) 次矩陣與矩陣的乘法。

參考文獻

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用