特征多项式学习笔记

Dementor · 2022-01-17 16:16:28 · 题解

特征多项式

本篇文章仅仅是根据自己的理解写的学习笔记

特征多项式是矩阵的一个知识点，在线性代数中有重要的应用

目标时间复杂度 O (n ^ 3)

总体过程分为 2 步走：

求上海森堡矩阵
递推求特征多项式

定义

对于一个 n \times n 的矩阵 A ，定义它的特征多项式为：

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

其中 I_n 是一个 n 阶的单位矩阵，最后的 p_A (x) 是一个 n 次多项式

求特征多项式

暴力可以直接取 0,1,...,n 一共 n + 1 个值代入 x 求行列式，再用多项式插值求出特征多项式，复杂度是 O (n ^ 4) 的，在这里不细讲

如果一个矩阵次对角线下方的位置全部都是 0 那么把其称为 “上海森堡矩阵”

那么接下来有一个算法一般叫做 “海森堡算法” 可以在 O (n ^ 3) 的时间内递推出特征多项式

将一个普通矩阵等价变换成我们需要的上海森堡矩阵的方法是 类高斯消元

高斯消元大家都非常熟悉了，但是我们需要在这基础上做一些调整， 这非常重要 ，千万不能遗漏

结论

对于初等可逆矩阵 P ，有：

\det (xI_n - A) = \det (xI_n - PAP^{-1})

证明：

\det (xI_n - PAP^{-1}) = \det (xPI_nP^{-1} - PAP^{-1})

= \det ((xPI_n - PA)P^{-1})

= \det ((PxI_n - PA)P^{-1})

= \det (P(xI_n - A)P^{-1})

= \det (P) \times \det (xI_n - A) \times \det (P^{-1})

= \det (xI_n - A)

我们称 B = PAP^{-1} 为 A 的相似矩阵，于是我们得到了一个表述：

相似矩阵的特征多项式相同

由此，我们也加深了对 “特征” 二字的理解

更改高斯消元

一般的高斯消元可以理解为令一个普通矩阵 A 变成 PA ，这个 P 即为一个初等矩阵

我们发现，此时的 A' 并不与原来的 A 相似，但它们也就相差了一个 P^{-1}

高斯消元用到的三个操作：

交换两行
将一行 \times k 加到另一行上
将一行 \times k

经过验证，这三个操作对应的初等矩阵 P 所需要多乘的 P^{-1} 实际上就是对列做一次一模一样的操作

同时发现，这样的操作是不会影响我们消成上海森堡矩阵的过程的

至此，我们完成了 第一步

递推

记 p_i(x) 为保留 A[1..i][1..i] 时的特征多项式

此时 A 是一个与读入矩阵相似的上海森堡矩阵

首先可以列出初值：

p_0(x) = \{1\}

p_1(x) = \{-A[1][1], 1\}

这里直接默认列表位置为对应多项式幂次的系数，没有写出来的不存在（即那个幂次系数为 0）

对 p_2 (x) = \det (xI_2 - A_2) 直接展开

p_2 (x) = (x - A[2][2]) p_1 (x) - A[2][1] \times A[1][2] p_0 (x)

对 p_3 (x) 、甚至是 p_4 (x) 展开，计算，归纳之后

发现可以写成递推式：

p_i (x) = (x - A[i][i])p_{i-1}(x) - \sum\limits_{m=1}^{i-1}A[i-m][i](\prod\limits_{j=i-m+1}^{i}A[j][j - 1])p_{i-m-1}(x)

这个算法的名字就是“海森堡算法”

至此，我们完成了 第二步

时间复杂度

高斯消元的复杂度 O (n ^ 3)

海森堡算法递推的复杂度是 O (n ^ 3)

最后的答案即为 p_n (x)

综上，时间复杂度是 O (n ^ 3)

达到了我们的目标

代码实现

从整个算法步骤知道代码其实不复杂，一步一步写就能成功实现算法

code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 510;
const ll mod = 998244353;
ll n, a[N][N], p[N][N];
ll qpow (ll a, ll b) {
    ll ans = 1ll;
    while (b) {
        if (b & 1) ans = ans * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ans;
}
ll inv (ll x) {return qpow (x, mod - 2);}
void Hesb (ll a[N][N], int n) {
    for (register int i=1;i<n;i++) {
        int p = -1;
        for (register int j=i+1;j<=n;j++) {
            if (a[j][i]) {
                p = j;
                break ;
            }
        }
        if (p == -1) continue ;
        else {
            swap (a[p], a[i + 1]);
            for (register int j=1;j<=n;j++) swap (a[j][p], a[j][i + 1]);
        }
        for (register int j=i+2;j<=n;j++) {
            ll tmp = a[j][i] * inv (a[i + 1][i]) % mod;
            for (register int k=1;k<=n;k++) a[j][k] = (a[j][k] + mod - tmp * a[i + 1][k] % mod) % mod;
            for (register int k=1;k<=n;k++) a[k][i + 1] = (a[k][i + 1] + tmp * a[k][j] % mod) % mod;
        }
    }
}
ll neg (ll x) {
    return mod - x;
}
int main () {
    cin >> n;
    for (register int i=1;i<=n;i++) {
        for (register int j=1;j<=n;j++) {
            scanf ("%lld\n", &a[i][j]);
        }
    }
    Hesb (a, n);
    p[0][0] = 1;
    for (register int i=1;i<=n;i++) {
        for (register int m=1;m<i;m++) {
            ll tmp = a[m][i];
            for (register int j=m;j<i;j++) tmp = tmp * a[j + 1][j] % mod;
            for (register int j=0;j<=n;j++) p[i][j] = (p[i][j] + neg (tmp * p[m - 1][j] % mod)) % mod;
        }
        for (register int j=1;j<=n;j++) p[i][j] = (p[i][j] + p[i - 1][j - 1]) % mod;
        for (register int j=0;j<=n;j++) p[i][j] = (p[i][j] + neg (a[i][i] * p[i - 1][j] % mod)) % mod;
    }
    for (register int i=0;i<=n;i++) printf ("%lld ", p[n][i]); puts ("");
}

可以通过洛谷上的模板题：P7776 特征多项式

参考文献

oi-wiki 特征多项式

特征多项式 学习笔记