特征多项式与常系数齐次线性递推

定义

对于 $n$ 阶矩阵，若存在非零列向量 $x$ 和数 $\lambda$ 满足 $Ax=\lambda x$，则称 $\lambda$ 和 $x$ 为一组对应的特征值和特征向量。

在确定了特征值之后，可以得到对应的无穷多个解。

$$\\ $$

求解特征值和特征向量

求解特征值和特征向量：

容易发现，$\lambda$ 是一个特征值，只需要满足 $Ax=\lambda x$ 有解，以 $x$ 为元容易列出方程，其常数项为均 $0$，系数矩阵为。

$$\left( \begin{array}{cccc}\lambda-A_{1,1}& -A_{1,2}& -A_{1,3}& \cdots & -A_{1,n}\\ -A_{2,1}&\lambda- A_{2,2} & -A_{2,3}& \cdots & -A_{2,n}\\ -A_{3,1} & -A_{3,2} & \lambda-A_{3,3} & \cdots & -A_{3,n}\\ \vdots& \vdots &\vdots &\vdots &\vdots \\ -A_{n,1} & -A_{n,2} & -A_{n,3} &\cdots & \lambda-A_{n,n}\end{array}\right) =\lambda I-A$$

其中$I$是单位矩阵。

这个方程有非零解的充要条件是：$det(\lambda I-A)=0$ （因为如果不为 $0$，则矩阵满秩，所有向量线性无关，无法得到0向量）

而 $det(\lambda I-A)$ 是一个 $n$ 次多项式 $p(\lambda)$，称为特征多项式，所有的特征值 $\lambda$ 就是 $p(\lambda)$ 的根。

$$\\ $$

应用

加速矩阵乘法：

由 $Ax=\lambda x$，迭代该式可以得到 $A^nx=\lambda^nx$ 。

特殊矩阵的特征值

上三角矩阵

$$\lambda I-A=\left( \begin{array}{cccc}\lambda-A_{1,1}& -A_{1,2}& -A_{1,3}& \cdots & -A_{1,n}\\ 0 & \lambda-A_{2,2} & -A_{2,3}& \cdots & -A_{2,n}\\ 0 & 0 & \lambda-A_{3,3} & \cdots & -A_{3,n}\\ \vdots& \vdots &\vdots &\vdots &\vdots \\ 0 & 0 & 0 &\cdots & \lambda-A_{n,n}\end{array}\right)$$

带入行列式即可知道 $\displaystyle det(\lambda I-A)=\prod (\lambda -A_{i,i})$ 。

也就是说，主对角线上所有的 $A_{i,i}$ 都是 $det(\lambda I-A)=0$ 的根。

$$\\ $$

零化多项式

对于一个矩阵 $A$，它的一个零化多项式 $\varphi(\lambda)$ 是满足 $\varphi(A)=0$ 的多项式，定义域包含矩阵。

最小多项式：次数最低的零化多项式。

$\text{Cayley-Hamilton}$ 定理

特征多项式：$p(\lambda)=|\lambda I-A|$，$\lambda$ 定义域不止是 $\mathbb{ R }$，还可以是矩阵。

$\text{Cayley-Hamilton}$ 定理指出：矩阵的特征多项式也是它的零化多项式。

即令：

$$\varphi(\lambda)=det(\lambda I-A)=\lambda^n+a_1\lambda^{n-1}+\cdots+a_{n-1}\lambda +a_n$$

则有：

$$\varphi(A)=A^n+a_1A^{n-1}+\cdots+a_{n-1}A +a_n=O$$

证明：

将 $\varphi(\lambda)$ 改写为：

$$\varphi(\lambda)=(\lambda-\lambda_1)(\lambda-\lambda_2)\cdots(\lambda-\lambda_n)$$

由定理：任意的 $n$ 阶矩阵都能相似为上三角矩阵 可知，存在可逆矩阵 $P$，使得：

$$PAP^{-1}= \left( \begin{array}{cccc} \lambda_1&*&\cdots&*\\ &\lambda_2&\ddots&\vdots\\ &&\ddots&*\\ &&&\lambda_n \end{array}\right)$$

将 $PAP^{-1}$ 代入 $\varphi(\lambda)$ 得到：

$$\varphi(PAP^{-1})=(PAP^{-1}-\lambda_1 I)(PAP^{-1}-\lambda_2 I)\cdots(PAP^{-1}-\lambda_n I)$$

计算：

$$\left( \begin{array}{cccc} 0&*&\cdots&*\\ &\lambda_2-\lambda_1&\ddots&\vdots\\ &&\ddots&*\\ &&&\lambda_n-\lambda_1 \end{array}\right) \times \left( \begin{array}{cccc} \lambda_1-\lambda_2&*&\cdots&*\\ &0&\ddots&\vdots\\ &&\ddots&*\\ &&&\lambda_n-\lambda_2 \end{array}\right) \times \cdots \left( \begin{array}{cccc} \lambda_1-\lambda_n&*&\cdots&*\\ &\lambda_2-\lambda_n&\ddots&\vdots\\ &&\ddots&*\\ &&&0 \end{array}\right)=O$$

即 $\varphi(PAP^{-1})=P\varphi(A)P^{-1}=O$，故有 $\varphi(A)=O$

$$\\ $$

求解特征多项式

带入 $n$ 个数，求出得 $det(x I_n-A)$，得到 $n$ 个矩阵，通过高斯消元可以 $O(n^3)$ 地求出行列式。

然后可 $O(n^2)$ 拉格朗日插值求出原来的多项式，总复杂度受限于高斯消元，为 $O(n^4)$ 。

$$\\ $$

求解最小多项式

构造矩阵序列 $a_i=A^i$。

求出它的一个线性递推 $r_i$，即：

$$\displaystyle \sum_{j=0}^{m} r_j a_{i-j}=\sum_{j=0}^{m} r_j A^{i-j}=(\sum_{j=0}^m r_{m-j}A^j)\cdot A^{i-m}=0 \\ \displaystyle \therefore \sum_{j=0}^m r_{m-j}A^j=0$$

所以可以由 $r_i$ 翻转得到 $f(\lambda)$ 。

求解 $a_i$ 前 $n$ 项的复杂度受限于矩阵乘法为 $O(n^4)$，求解递推式的复杂度为 $O(n^3)$ 。

考虑到实际求解递推式时，随机生成了两个向量 $u,v$ 。

实际是计算标量序列 $\{uA^iv\}$ 的递推式，所以实际每次求出 $uA^i$ 复杂度应为 $O(n^2)$ 。

求这个递推式需要用到 $a_i$ 前 $2n$ 项，求解复杂度为 $O(n^3)$ 。

因此总复杂度为 $O(n^3)$ 。

(但是如果只是求出来并没有什么用，因为求解方法是随机的，甚至连检查一次保证正确都需要 $O(n^2(n+e))$ 的时间($e$ 为矩阵非 $0$ 位置个数))

$$\\ $$

求解稀疏方程组

设方程系数用矩阵 $A$ 表示，右侧每个方程的常数用向量 $b$ 表示，答案用向量 $x$ 表示，则满足关系式

$Ax=b$，即 $x=A^{-1}b$ 。

求出 $\{A^ib\}$ 线性递推式，反推出 $A^{-1}b$ 即可。

反推方法：

带入线性递推的 $m$ 项，则 $\sum_{i=0}^{m}\limits A^{m-i}b\cdot r_i=0$

两边同乘 $A^{-1}$，得到 $A^{-1}b\cdot r_m +\sum_{i=0}^{m-1}\limits A^{m-i}br_i=0$

求解矩阵 k 次幂

我们要求解 $A^k$，常规做法是直接用快速幂

设矩阵 $A$ 的一个零化多项式是 $f(\lambda)$ ，显然，$A^k$ 可以用一个多项式表示 $A^k=\sum_0^k w_i A^i$。

$\{w_i\}$ 构成了一个 $k+1$ 次多项式 $F_k(x)$。

存在一种合法的表示是 $F_k(x)=x^k$。

$\because f(A)=0\\ \therefore \forall i, f(A)A^i=0\\ $

所以对于任意实数 $T$，$G_k(x)=x^k-Tf(x)$ 也合法，也就是相当于我们要求出 $x^k$ 对于 $f(x)$ 这个 $n+1$ 多项式取模。

显然可以通过类似快速幂的方式倍增求解这个多项式，每次对 $f(x)$ 取模复杂度是 $O(n\log n)$ ，总时间复杂度 $O(n\log m\log n)$ 。

最后得到的 $F(x)$ 是一个 $n$ 次多项式，带入就可以快速求出 $A_k$，可以认为这个复杂度是受限于求解 $A^0,A^1,\cdots,A^{n-1}$ 的 $O(n^4)$ 。

对于元矩阵 $A$ 为稀疏矩阵的情况，设其包含 $e$ 个非零位置，那么求解 $B\cdot A$ 的过程是 $O(n\cdot e)$ 的，求解 $A_0,A^1,\cdots,A^{n-1}$ 的过程，是 $O(n^2e)$ 的。

求解零化多项式的复杂度也是 $O(n^2(n+e))$ 的，因此总复杂度为 $O(n^2(n+e))$。

而一般的矩阵快速幂是 $O(n^3\log k)$ 的，这种方法适用情况非常特殊。

另外，对于并不需要知道整个矩阵的答案，并且 $A^0,A^1,\cdots,A^{n-1}$ 特殊的具体问题，这个方法也十分有效。

求解常系数线性齐次递推

问题是要求数列 $f_i=\sum _{j=1}^{n}a_j\cdot f_{i-j}$ 。

给出 $f_0,f_1,\cdots,f_{n-1}$，求第 $k$ 项的值。

线性递推显然可以用初始向量列与转移矩阵的幂次的乘积表示，即$f_i=(S \cdot A^i)_n$，其中$A$为转移矩阵，$S$为初始向量列，我们求的是第$n$项。

对于$n=4$的情况，我们的转移矩阵。

$$A=\left( \begin{array}{cccc} & & & a_4\\ 1 & & & a_3 \\ & 1& & a_2 \\ & & 1 & a_1\end{array}\right)$$

鉴于它的特殊性，我们可以直接求出它的特征多项式表达式。

由

$$\lambda I_n-A= \left( \begin{array}{cccc} \lambda& & & -a_4\\ -1 & \lambda & & -a_3 \\ & -1& \lambda & -a_2 \\ & & -1 & \lambda-a_1\end{array}\right)$$

带入行列式最暴力的求法

枚举一个排列$p_i$，设排列$p$的逆序对为$f(p)$，$|A|=\sum (-1)^{f(p)} \Pi A_{i,p_i}$

实际上合法的排列只有$n$个，枚举 $p_i=n$

那么：

$$p_j=\left\{\begin{aligned} j && j<i \\ n && j=i \\ j-1 && j> i\end{aligned}\right.$$

当 $i=n$ 时，$(-1)^{f(p)} \prod A_{i,p_i}=\lambda ^n-a_1\lambda ^{n-1}$

当 $i\neq n$ 时：

$$\because\left\{\begin{array}{cccc}f(p)=n-i\\ \\ \prod A_{i,p_i}=(-1)^{n-i+1}\lambda^i\cdot a_{n-i+1}\end{array}\right.$$

$$\therefore\ (-1)^{f(p)} \prod A_{i,p_i}=-\lambda^i a_{n-i+1}$$

综上，转移矩阵 $A$ 的特征多项式有简单的表达

$$p(\lambda) = |\lambda I_n-A|=\lambda^n-a_1\lambda^{n-1} -a_2\lambda^{n-2} -\cdots -a^n$$

假设有 $f_0$ 这一项(不需要知道是多少)，那么认为初始向量列为 $S=(f_{-(n-1)},f_{-(n-2)},\cdots ,f_{0})$ 。

这个问题，我们要求的是 $S\cdot A^k$ 的第 $n$ 项，不需要知道整个矩阵。

类似求出 $A^k$ 的过程，求出 $F_k(x)\mod p(\lambda)$。

我们要求解 $(S\cdot A^k)_n=\sum_1^{n}[x^i]{F(x)}(S\cdot A^i)_n$，而 $(S\cdot A^i)_n=f_i$ 已知，求出 $F(x)$ 后直接带入即可

需要用到多项式取模，求解这个表达式是 $O(n\log n\log k)$ 的，求完直接带入即可。

模板代码

#4161. Shlw loves matrixI

非 $\text{NTT}$ 模数，$k\le2000$ ，可以 $O(k^2)$ 卷积，总复杂度 $O(k^2\log n)$。

点击查看代码

#include <bits/stdc++.h>
using namespace std;
int n, m;
const int md = 1e9 + 7;
inline int pwr(int x, int y) {
    int res = 1;
    while (y) {
        if (y & 1)
            res = 1ll * res * x % md;
        x = 1ll * x * x % md;
        y >>= 1;
    }
    return res;
}
struct Poly {
    vector<int> dp;
    Poly(int x = 0) { dp.resize(x); }
    inline auto size() { return dp.size(); }
    inline int& operator[](int t) { return dp[t]; }
    inline Poly Reverse() {
        Poly res = *this;
        reverse(res.dp.begin(), res.dp.end());
        return res;
    }
    inline Poly operator-(Poly b) {
        Poly res(max(dp.size(), b.size()));
        for (int i = 0; i < dp.size(); i++) res[i] = (res[i] + dp[i] + md) % md;
        for (int i = 0; i < b.size(); i++) res[i] = (res[i] - b[i] + md) % md;
        return res;
    }
    inline Poly operator*(int x) {
        Poly res = *this;
        for (auto& it : res.dp) it = 1ll * it * x % md;
        return res;
    }
    inline Poly operator*(Poly y) {
        Poly res(dp.size() + y.size() - 1);
        for (int i = 0; i < dp.size(); i++) {
            for (int j = 0; j < y.size(); j++) res[i + j] = (res[i + j] + 1ll * dp[i] * y[j]) % md;
        }
        return res;
    }
    inline Poly Inv(int len) {
        if (len == 1) {
            Poly res(1);
            res[0] = pwr(dp[0], md - 2);
            return res;
        }
        Poly res = this->Inv((len + 1) / 2), tmp = *this;
        tmp.dp.resize(len);
        res = res * 2 - tmp * res * res;
        res.dp.resize(len);
        return res;
    }
    inline Poly operator/(Poly b) {
        int N = dp.size(), M = b.size();
        if (N < M)
            return Poly();
        Poly res = this->Reverse() * (b.Reverse().Inv(N - M + 1));
        res.dp.resize(N - M + 1);
        return res.Reverse();
    }
    inline Poly operator%(Poly b) {
        Poly G = *this / b, res = *this - (G * b);
        res.dp.resize(min(dp.size(), b.size()));
        return res;
    }
};
inline Poly Polypower(Poly x, int y, Poly MD) {
    Poly res(1);
    res[0] = 1;
    while (y) {
        if (y & 1)
            res = res * x % MD;
        x = x * x % MD;
        y >>= 1;
    }
    return res;
}
int a[40005], k;
int main() {
    scanf("%d%d", &k, &n);
    Poly F(n + 1), X(2);
    F[n] = 1;
    X[1] = 1;
    for (int i = n - 1; ~i; i--) scanf("%d", &F[i]), F[i] = (md - F[i]) % md;
    for (int i = 0; i < n; i++) scanf("%d", &a[i]);
    Poly tmp = Polypower(X, k, F);
    int res = 0;
    for (int i = 0; i < n; i++) res = (res + 1ll * a[i] * tmp[i]) % md;
    printf("%d\n", (res + md) % md);

    return 0;
}

P4723 【模板】常系数齐次线性递推

$\text{NTT}$ 模数，$k\le32000$ ，$O(k\log k)$ 卷积，总复杂度 $O(k\log k\log n)$。

点击查看代码

#include <bits/stdc++.h>
using namespace std;
int n, m;
const int md = 998244353, G = 3, Gi = (md + 1) / 3;
int r[1 << 20];
inline int pwr(int x, int y) {
    int res = 1;
    while (y) {
        if (y & 1)
            res = 1ll * res * x % md;
        x = 1ll * x * x % md;
        y >>= 1;
    }
    return res;
}
inline void NTT(vector<int> &dp, int lim, int W) {
    for (int i = 0; i < (1 << lim); i++)
        if (i < r[i])
            swap(dp[i], dp[r[i]]);
    for (int i = 0; i < lim; i++) {
        int w = pwr(W, (md - 1) / (1 << (i + 1)));
        for (int j = 0; j < (1 << lim); j += (1 << (i + 1))) {
            int Pw = 1;
            for (int t = 0; t < (1 << i); t++) {
                int x = dp[j + t], y = 1ll * Pw * dp[j + (1 << i) + t] % md;
                dp[j + t] = (x + y) % md;
                dp[j + (1 << i) + t] = (x - y + md) % md;
                Pw = 1ll * Pw * w % md;
            }
        }
    }
}
struct Poly {
    vector<int> dp;
    Poly(int x = 0) { dp.resize(x); }
    inline auto size() { return dp.size(); }
    inline int &operator[](int t) { return dp[t]; }
    inline Poly Reverse() {
        Poly res = *this;
        reverse(res.dp.begin(), res.dp.end());
        return res;
    }
    inline Poly operator-(Poly b) {
        Poly res(max(dp.size(), b.size()));
        for (int i = 0; i < dp.size(); i++) res[i] = (res[i] + dp[i] + md) % md;
        for (int i = 0; i < b.size(); i++) res[i] = (res[i] - b[i] + md) % md;
        return res;
    }
    inline Poly operator*(int x) {
        Poly res = *this;
        for (auto &it : res.dp) it = 1ll * it * x % md;
        return res;
    }
    inline Poly operator*(Poly y) {
        Poly x = *this, res(dp.size() + y.size() - 1);
        int lim = 0;
        while ((1 << lim) <= res.size()) lim++;
        for (int i = 0; i < (1 << lim); i++) r[i] = (r[i >> 1] >> 1) + ((i & 1) << (lim - 1));
        vector<int> a(1 << lim), b(1 << lim), c(1 << lim);
        for (int i = 0; i < dp.size(); i++) a[i] = dp[i];
        for (int i = 0; i < y.size(); i++) b[i] = y[i];
        NTT(a, lim, G);
        NTT(b, lim, G);
        for (int i = 0; i < (1 << lim); i++) c[i] = 1ll * a[i] * b[i] % md;
        NTT(c, lim, Gi);
        int inv = pwr(1 << lim, md - 2);
        for (int i = 0; i < res.size(); i++) res[i] = 1ll * c[i] * inv % md;
        return res;
    }
    inline Poly Inv(int len) {
        if (len == 1) {
            Poly res(1);
            res[0] = pwr(dp[0], md - 2);
            return res;
        }
        Poly res = this->Inv((len + 1) / 2), tmp = *this;
        tmp.dp.resize(len);
        res = res * 2 - tmp * res * res;
        res.dp.resize(len);
        return res;
    }
    inline Poly operator/(Poly b) {
        int N = dp.size(), M = b.size();
        if (N < M)
            return Poly();
        Poly res = this->Reverse() * (b.Reverse().Inv(N - M + 1));
        res.dp.resize(N - M + 1);
        return res.Reverse();
    }
    inline Poly operator%(Poly b) {
        Poly G = *this / b, res = *this - (G * b);
        res.dp.resize(min(dp.size(), b.size()));
        return res;
    }
};
inline Poly Polypower(Poly x, int y, Poly MD) {
    Poly res(1);
    res[0] = 1;
    while (y) {
        if (y & 1)
            res = res * x % MD;
        x = x * x % MD;
        y >>= 1;
    }
    return res;
}
int a[40005], k;
int main() {
    scanf("%d%d", &k, &n);
    Poly F(n + 1), X(2);
    F[n] = 1;
    X[1] = 1;
    for (int i = n - 1; ~i; i--) scanf("%d", &F[i]), F[i] = (md - F[i]) % md;
    for (int i = 0; i < n; i++) scanf("%d", &a[i]);
    Poly tmp = Polypower(X, k, F);
    int res = 0;
    for (int i = 0; i < n; i++) res = (res + 1ll * a[i] * tmp[i]) % md;
    printf("%d\n", (res + md) % md);

    return 0;
}