特征多项式优化线性递推

发表于 | 更新于 | 分类于 | 评论数： | 阅读次数：

欢迎大家去看叉姐的论文啊~

欢迎大家去看这篇文章的姊妹篇：利用简单线性表示优化线性递推

假设大家都学过矩阵快速幂了，下面主要讲将特征多项式黑科技。

温习：矩阵快速幂

对于一个 $k$ 次的线性递推:

$ h_n = a_1 h_{n-1} + a_2 h_{n-2} + ... + a_k h_{n-k} , \forall\ n \in \mathbb{N}\ ,\ n > k $

我们构造转移矩阵

$ A = \begin{bmatrix} a_1 & a_2 & \cdots & a_{k-2} & a_{k-1} & a_k\\ 1 & 0 & \cdots & 0 & 0 & 0\\ 0 & 1 & \cdots & 0 & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & 0 & 0\\ 0 & 0 & \cdots & 0 & 1 & 0 \end{bmatrix}_{k \times k}$

和初始向量

$x = \begin{bmatrix} h_{k-1} \\ h_{k-2} \\ \vdots \\ h_2 \\ h_1 \\ h_0 \end{bmatrix}_{k \times 1}$

易见

$Ax = \begin{bmatrix} a_1 & a_2 & \cdots & a_{k-2} & a_{k-1} & a_k\\ 1 & 0 & \cdots & 0 & 0 & 0\\ 0 & 1 & \cdots & 0 & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & 0 & 0\\ 0 & 0 & \cdots & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} h_{k-1} \\ h_{k-2} \\ \vdots \\ h_2 \\ h_1 \\ h_0 \end{bmatrix} = \begin{bmatrix} h_k \\ h_{k-1} \\ \vdots \\ h_3 \\ h_2 \\ h_1 \end{bmatrix} $

显然这样可以顺次推下去，于是我们只要得到 $A^{n-k+1}x$ 就可以得到 $ h_n $ ，由于一次矩阵乘法复杂度为 $ O(k^3) $ (听说可以更小一点，但我不会。。。) 矩阵乘法又满足快速幂性质，可以类比普通乘法快速幂得到一个 $O(k^3 \log n)$ 的优秀做法。

但是这个复杂度依然不够优秀，我们把眼光放到线性代数上。

线性代数

我在这章将直接给出结论。在本文的结尾会补充一些线性代数的内容。

特征值和特征向量

若非零向量 $x$ 有 $Ax = \lambda x\ $ , 则称 $x$ 为 $A$ 的一个特征向量 ， $\lambda$ 是其对应的特征值。

特征多项式

我们记矩阵 $A$ 的特征多项式 $f(\lambda)=|A - \lambda I|$。 令 $f(\lambda)=0$ 解得的 $\lambda$ 即为矩阵 $A$ 的特征值。

化零多项式

对于一个多项式 $ f(x) = p_0 x^0 + p_1 x^1 + … + p_n x^n $ 和一个矩阵 $ A $ , 若有 $f(A) = p_0 A^0 + p_1 A^1 + … + p_n A^n = 0$ ，则称 $f(A)$ 是矩阵 $A$ 的一个化零多项式。

Caylay-Camilton 定理

凯莱-哈密尔顿定理十分简洁。

Caylay-Camilton 定理: 一个矩阵的特征多项式是它的化零多项式

即若有矩阵 $A$ 的特征多项式 $f(\lambda)= p_0 \lambda^0 + p_1 \lambda^1 + … + p_n \lambda^n $ ，则 $f(A) = p_0 A^0 + p_1 A^1 + … + p_n A^n = 0$ 。

特征多项式优化线性递推

线性递推的小结论

对于线性递推 $a_{n+k}=c_1a_{n+k-1}+c_2a_{n+k-2}+…+c_ka_n$ ，它的转移矩阵的特征多项式是

$A^k=c_1A^{k-1}+c_2A^{k-2}+...+c_{n-1}A+c_nI$

优化开始

我们其实发现矩阵快速幂求得就是 $A^{n-k+1}x$ 。我们考虑将 $A^x$ 降次。

其实我们发现 $A^x=c_1A^{x-1}+c_2A^{x-2}+…+c_kA^{x-k}$ ，就可以把 $A^x$ 降次为 $A^{x-1},A^{x-2},…,A^{x-k}$ ，如果我们按着这样降次下去，就可以将 $A^{n-k+1}$ 降次到 $I,A,A^2,…,A^{k-1}$ 的线性组合。这样做的时间复杂度是 $O(nk)$ ，这显然是不够优秀的。

不过，如果我们得到 $A^n=p_0I+p_1A+p_2A^2+…+p_{k-1}A^{k-1}$ 和 $A^m=q_0I+q_1A+q_2A^2+…+q_{k-1}A^{k-1}$ ，就有 $A^{n+m}=A^nA^m$ $=(p_0I+p_1A+p_2A^2+…+p_{k-1}A^{k-1})(q_0I+q_1A+q_2A^2+…+q_{k-1}A^{k-1})$ $=\sum_{i=0}^{2k-2} \sum_{j=0}^i p_jq_{i-j} A^i$

我们在这里得到一个次数为 $2k-2$ 的多项式，我们用上面降次的方法将多项式降到最高次数 $k-1$ 的多项式。

如果我们暴力计算多项式乘积的话，多项式乘积的时间复杂度是 $O(k^2)$ 的。如果我们暴力降次的话，时间复杂度是 $O(k^2)$ 。于是一次乘法的总复杂度为 $O(k^2)$ 。利用快速幂，我们可以做到 $O(k^2 \log n)$ 。

如果想要继续优化的话，我们可以利用FFT来代替暴力多项式乘法，复杂度 $O(k \log k)$。利用多项式取模来代替暴力降次，复杂度 $O(k \log k)$ 。于是一次乘法总复杂度为 $O(k \log k)$ 。利用快速幂，我们可以做到 $O(k \log k \log n)$ 。

小提示

  1. 当然了我们不会直接记 $A^x$ ，我们只需要记 $A^x=p_0I+p_1A+p_2A^2+…+p_{k-1}A^{k-1}$ 中的 $\{p_0,p_1,…,p_{k-1}\}$ 就可以了
  2. 我们也我们也不需要直接算出 $I,A,A^2,…,A^{k-1}$ 这样的话的复杂度是 $O(k^3)$ 这是我们不能接受的。事实上，我们只需要计算出 $x,Ax,A^2x,A^3x,…,A^{k-1}x$ 就可以了。时间复杂度是 $O(k^2)$ 。
  3. 我们甚至不需要算 $x,Ax,A^2x,…,A^{k-1}x$。因为我们最后关注的 $a_n$ 其实只是 $A^{n-k+1}x$ 的最后一项，我们只需要考虑 $x,Ax,A^2x,…A^{k-1}x$ 的最后一项就可以了。在这里 $x,Ax,A^2x,…,A^{k-1}x$ 的最后一项正好分别是 $a_0,a_1,…,a_{k-1}$
  4. 总结一下 $2,3$ 就是说：在线性递推中，如果我们算到 $A^{n-k+1}=p_0I+p_1A+p_2A^2+…+p_{k-1}A^{k-1}$ 那么 $a_n=p_0a_0+p_1a_1+…+p_{k-1}a_{k-1}$

例：Project Euler 258 A lagged Fibonacci sequence

题意

定义递推数列 $\{g_n\}$ ，满足

  • 初值: 对于 $n < 2000$,有 $g_n = 1$，
  • 递推式：对于 $n \geq 2000$ ,有 $g_n=g_{n-2000}+g_{n-1999}$

求 $g_{10^{18}} \mod 20092010$。

分析

这有什么好分析的。就是两千阶的线性递推吗。直接上代码得了。

代码

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
#include<bits/stdc++.h> 

using namespace std;

const long long mod=20092010;
const int k=2000;

void mul(long long *p,long long *q,long long *c) //p,q为A^n,A^m降次后的系数，c为递推系数数组，函数算出a^{m+n}降次之后的系数并存在数组p中
{
    long long tmp[2*k];
    memset(tmp,0,sizeof(tmp));
    for(int i=0;i<k;i++)
        for(int j=0;j<k;j++)
            tmp[i+j]=(tmp[i+j]+p[i]*q[j])%mod; //利用可加性将a_{m+n}利用(a_0,a_1,...,a_{2k-2})线性表示，并把系数存在tmp中
    for(int i=2*k-2;i>=k;i--) //利用a_n=c_0a_{n-1}+c_1a_{n-2}+...+c_ka_{n-k} 将 a_i 转化为 a_{i-1},a_{i-2},...,a_{i-k}
        if(tmp[i])
            for(int j=0;j<k;j++) //这样依次转化 a_{2k-2},a_{2k-3},...,a_k 就可以得到我们想要的线性表示
                if(c[j])
                    tmp[i-j-1]=(tmp[i-j-1]+tmp[i]*c[j])%mod;
    for(int i=0;i<k;i++)
        p[i]=tmp[i];
}

long long solve(long long *a,long long *c,long long n) //a为初始值数组 c为系数数组 和上文线性递推定义对应
{
//  if(k==1) return pw(c[0],n-1)*a[0]%mod;
    if(n<k) return a[n]; //如果已知就直接返回
    long long tmp[k],res[k]; //初始化快速幂，tmp存 A^{2^x} ，res存答案
    memset(tmp,0,sizeof(tmp));memset(res,0,sizeof(res));
    tmp[1]=res[0]=1;
    while(n>0)
    {
        if(n&1) mul(res,tmp,c);
        mul(tmp,tmp,c);
        n>>=1;
    }
    long long ans=0;
    for(int i=0;i<k;i++)
        ans=(ans+a[i]*res[i])%mod; //a_n=p_0a_0+p_1a_1+...+p_{k-1}a_{k-1}
    return ans;
}

int main()
{
    long long a[k],c[k];
    fill(a,a+k,1);
    memset(c,0,sizeof(c));
    c[1998]=c[1999]=1;
    cout<<solve(a,c,1e18)<<endl;
    return 0;
}

补充知识

特征向量和特征值

我们考虑矩阵的几何意义。矩阵其实等价于坐标的线性变换。比方说一个 $2\times2$ 的方阵可以通过乘上一个二维平面直角坐标，得到新的二维平面直角坐标。这就是坐标变换。
我们并不喜欢用坐标这个词，这很不线性代数。我们引入向量————有方向的量，来代替坐标。一个坐标对应的向量，长度为坐标到原点的距离，方向为从原点出发，经过坐标的射线的方向。
有趣的是，在坐标变换前后，有些向量方向改变，有些向量方向则不变（或变成其反方向）。
我们称那些坐标变换前后方向不变或者变成其反方向的非零向量叫作特征向量，其前后的长度之比记为其特征值。特征值的正负号取决于方向是与原来相同还是恰好相反。

特征多项式

为了求出所有特征值，我们尝试解方程 $Ax = \lambda x$ 。

移项得到 $(A - \lambda I) x = 0 $ , 由于 $x$ 为非零向量 ， 故 $A-\lambda I$ 一定不可逆（否则等式两边左乘其逆矩阵，得到 $x = 0$) ，故有 $A - \lambda I$ 特征值为 $0$ ，即$|A - \lambda I|=0$。

相关文章
打赏还是要有的，万一真给了呢
0%