介绍了方阵的特征多项式以及利用上Hessenberg矩阵的 \(\mathcal{O}(n^3)\) 求法。
Reference
特征值与特征向量
给定 \(n \times n\) 的方阵 \(\mathbf{T}\)，若存在一个非零列向量 \(\mathbf{v}\) 和数 \(\lambda\) 满足 \(\mathbf{T}\mathbf{v} = \lambda \mathbf{v}\)，则称 \(\lambda\) 为 \(\mathbf{T}\) 的一个特征值，而 \(\mathbf{v}\) 是 \(\mathbf{T}\) 的属于特征值 \(\mathbf{v}\) 的一个特征向量。
也就是说，\(\mathbf{v}\) 在 \(\mathbf{T}\) 的线性变换下，只是伸缩了 \(\lambda\) 倍。
特征多项式
将上面的式子变换得到 \((\lambda \mathbf{I}-\mathbf{T})\mathbf{v}=\mathbf{0}\)，相应地，非 \(\mathbf{0}\) 解 \(\mathbf{v}\) 就是一个特征向量。我们称 \(\lambda \mathbf{I} - \mathbf{T}\) 为 \(T\) 的特征矩阵，特征矩阵的行列式 \(\det(\lambda \mathbf{I} - \mathbf{T})\) 是一个 \(n\) 次多项式，称为 \(\mathbf{T}\) 的特征多项式，其根就是 \(\mathbf{T}\) 的特征值，记做 \(p_{\mathbf{T}}(\lambda)\)。
怎么求一个方阵的特征多项式？代入 \(n+1\) 个 \(\lambda\) 后Lagrange插值即可，时间复杂度 \(\mathcal{O}(n^4)\)。
接下来，我们引入相似矩阵和上Hessenberg矩阵来做到 \(\mathcal{O}(n^3)\) 的时间复杂度。
更优秀的做法
Hessenberg矩阵
海森堡矩阵（Hessenberg matrix）是一个和三角阵很相似的方阵，上海森堡矩阵（upper Hessenberg matrix）的次对角线以下元素全为 \(0\)，下海森堡（lower Hessenberg matrix）矩阵的次对角线以上元素全为 \(0\)，下面是两个例子，左侧为上Hessenberg矩阵，右侧为下Hessenberg矩阵。
\[\begin{bmatrix} 1 & 3 & 4 & 2\\ 7 & 2 & 2 & 1\\ 0 & 1 & 4 & 5\\ 0 & 0 & 3 & 9 \end{bmatrix} \begin{bmatrix} 1 & 6 & 0 & 0\\ 7 & 2 & 2 & 0\\ 3 & 1 & 4 & 5\\ 5 & 8 & 3 & 9 \end{bmatrix} \nonumber \]Hessenberg矩阵对我们求特征多项式有什么用处吗？接下来我们试着求出一个上Hessenberg矩阵的特征多项式，并且在之后用相似变换将求一般方阵的特征多项式转换为求Hessenberg矩阵的特征多项式。
设上Hessenberg矩阵 \(\mathbf{H}\) 为
\[\begin{bmatrix} \alpha_1 & h_{1,2} & h_{1,3} & \cdots & h_{1,n} \\ \beta_2 & \alpha_2 & h_{2,3} & \cdots & h_{2,n} \\ & \ddots & \ddots & \ddots & \vdots \\ & & \ddots & \ddots & h_{n-1,n} \\ & & & \beta_n & \alpha_n \end{bmatrix} \nonumber \]那么其特征多项式就是
\[\begin{vmatrix} \lambda-\alpha_1 & -h_{1,2} & -h_{1,3} & \cdots & -h_{1,n} \\ -\beta_2 & \lambda-\alpha_2 & -h_{2,3} & \cdots & -h_{2,n} \\ & \ddots & \ddots & \ddots & \vdots \\ & & \ddots & \ddots & -h_{n-1,n} \\ & & & -\beta_n & \lambda-\alpha_n \end{vmatrix} \nonumber \]记 \(\mathbf{H}_i\) 表示保留 \(\mathbf{H}\) 的前 \(i\) 行前 \(i\) 列得到的方阵，\(p_i(\lambda)\) 表示 \(\mathbf{H}_i\) 的特征多项式。
按最后一行展开 \(p_{n}(\lambda)\)，得到
\[p_{n}(\lambda) = (\lambda - \alpha_n)p_{n-1}(\lambda) + \beta_n \begin{vmatrix} \lambda-\alpha_1 & -h_{1,2} & \cdots & -h_{1,n-2} & -h_{1,n} \\ -\beta_2 & \lambda-\alpha_2 & \cdots & -h_{2,n-2} & -h_{2,n} \\ & \ddots & \ddots & \vdots & \vdots \\ & & \ddots & \lambda - \alpha_{n-2} & -h_{n-2,n} \\ & & & -\beta_{n-1} & -h_{n-1,n} \end{vmatrix} \nonumber \]再按最后一列展开右面的行列式，对于每一个 \(-h_{n-i,n}\)，其余子式皆形如
\[\begin{vmatrix} & & & & & \\ & [\lambda \mathbf{I}-\mathbf{H}_{n-i-1}] & & & & \\ & & & & & \\ & & -\beta_{n-i+1} & \cdots & & \\ & & & \ddots & \vdots & \vdots \\ & & & & -\beta_{n-2} & \vdots \\ & & & & & -\beta_{n-1} \end{vmatrix} \nonumber \]这个东西很特别，它就是
\[p_{n-i-1}(\lambda) \prod_{j=n-i+1}^{n-1}-\beta_j \nonumber \]把每个 \(-h_{n-i,n}\) 展开的结果加起来，得到
\[\sum_{i=1}^{n-1}(-1)^{1+(n-i)+(n-1)+(i-1)} p_{n-i-1}(\lambda)h_{n-i,n}\prod_{j=n-i+1}^{n-1}\beta_j=-\sum_{i=1}^{n-1}p_{n-i-1}(\lambda)h_{n-i,n}\prod_{j=n-i+1}^{n-1}\beta_j \nonumber \]代回到原式子，我们就有
\[p_{n}(\lambda) = (\lambda - \alpha_n)p_{n-1}(\lambda) -\sum_{i=1}^{n-1}p_{n-i-1}(\lambda)h_{n-i,n}\prod_{j=n-i+1}^{n}\beta_j \nonumber \]可以递推计算 \(p_i(\lambda)\)，时间复杂度 \(\mathcal{O}(n^3)\)。
相似矩阵
对于一个 \(n \times n\) 方阵 \(\mathbf{A}\)，定义其相似矩阵 \(\mathbf{B}\) 是可以表示为 \(\mathbf{P} \times \mathbf{A} \times \mathbf{P}^{-1}\) 的 \(n \times n\) 方阵，其中 \(\mathbf{P}\) 是任意一个 \(n \times n\) 的可逆方阵，此时称 \(\mathbf{A}\) 和 \(\mathbf{B}\) 相似。
相似矩阵有什么特别之处吗？
定理： 两个相似矩阵 \(\mathbf{A}\) 和 \(\mathbf{B}\) 具有相同的特征多项式。
证明：
\[\begin{aligned} \det(\lambda \mathbf{I} - \mathbf{P}\mathbf{A}\mathbf{P}^{-1})&=\det(\mathbf{P}(\lambda \mathbf{I})\mathbf{P}^{-1} - \mathbf{P}\mathbf{A}\mathbf{P}^{-1})\newline &=\det(\mathbf{P})\det(\mathbf{P}^{-1})\det(\lambda \mathbf{I} - \mathbf{A})\newline &=\det(\lambda \mathbf{I} - \mathbf{A}) \end{aligned}\nonumber \]证毕！
运用相似变换直接将一个任意方阵 \(\mathbf{A}\) 变为上三角阵是不现实的：在消元时做第 \(i\) 列和第 \(j\) 列的初等列变换会把之前的东西变乱——但是运用相似变换将一个任意方阵 \(\mathbf{A}\) 变为一个Hessenberg矩阵是容易的！此时消元的初等列变换不会干扰到之前做好的值。
于是就得到了 \(\mathcal{O}(n^3)\) 求任意方阵的特征多项式的做法。
示例代码
const int MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
long long r = 1;
for(; b; b >>= 1, a = a * a % MOD)
if(b & 1) r = r * a % MOD;
return r;
}
namespace solver {
int Mat[N][N];
int p[N][N]; // p[i][k] 表示 p_i(Mat) 中 \lambda^k 的系数
inline vector<int> solve(int n) {
for(int i = 1; i <= n; i ++) {
if(Mat[i + 1][i] == 0) {
for(int j = i + 2; j <= n; j ++)
if(Mat[j][i] != 0) {
swap(Mat[i + 1], Mat[j]);
for(int k = 1; k <= n; k ++)
swap(Mat[k][i + 1], Mat[k][j]);
break;
}
}
if(Mat[i + 1][i] == 0) continue;
for(int j = i + 2; j <= n; j ++) {
long long mul = 1ll * Mat[j][i] * ksm(Mat[i + 1][i], MOD - 2) % MOD;
for(int k = 1; k <= n; k ++)
Mat[j][k] = Minus(Mat[j][k], mul * Mat[i + 1][k] % MOD);
for(int k = 1; k <= n; k ++)
Mat[k][i + 1] = Plus(Mat[k][i + 1], mul * Mat[k][j] % MOD);
}
}
for(int i = 0; i <= n; i ++)
for(int j = 0; j <= n; j ++)
p[i][j] = 0; // 只做一次的话可以不清空
p[0][0] = 1;
for(int i = 1; i <= n; i ++) {
p[i][0] = Minus(0, 1ll * Mat[i][i] * p[i - 1][0] % MOD);
for(int j = 1; j <= i; j ++)
p[i][j] = Minus(p[i - 1][j - 1], 1ll * Mat[i][i] * p[i - 1][j] % MOD);
for(int j = 1; j <= i - 1; j ++) {
long long mul = Mat[i - j][i];
for(int k = i - j + 1; k <= i; k ++)
mul = mul * Mat[k][k - 1] % MOD;
for(int k = 0; k <= i; k ++)
p[i][k] = Minus(p[i][k], p[i - j - 1][k] * mul % MOD);
}
}
vector<int> res(n + 1);
for(int i = 0; i <= n; i ++)
res[i] = p[n][i];
return res;
}
}