【学习笔记】动态规划—矩阵递推加速

【学习笔记】动态规划—矩阵递推加速

---

【大前言】

矩阵优化 (dp) 通常用于线性递推式的 (dp) 优化，能以优异的时间复杂度实现大量的状态转移。

更完整的 (dp) 优化策略：【学习笔记】动态规划—各种 (DP) 优化

---

[QAQ ]

---

一.【题目特征】

((1).) 类似线性递推（划重点，包括有向图上的递推等等）
((2).) 转移次数 (10^9) 左右（雾）
((3).) 决策点较少（常数）

---

[QAQ ]

---

二.【前置芝士】

1.【前言】

首先要清楚矩阵是个什么东西，在对 (dp) 进行优化时通常只会用到矩阵乘法和矩阵加法，其中矩阵乘法最为常见。

---

2.【运算法则】

(1).【矩阵加法】

矩阵加法的一般式：(C_{i,j}=A_{i,j} + B_{i,j})，其中 (A,B) 均为 (N imes N) 的矩阵，矩阵 (A + B) 得到 (N imes N) 的矩阵 (C) 。

式子的含义为：矩阵 (C) 由矩阵 (A,B) 对应位置上数值相加所得。

划重点：矩阵加法满足交换律 。

(2).【矩阵乘法】

矩阵乘法的一般式：(C_{i,j}=sum_{k=1}^{K} (A_{i,k} imes B_{k,j}))，其中 (A) 为 (N imes K) 的矩阵，(B) 为 (K*M) 的矩阵，矩阵 (A imes B) 得到 (N imes M) 的矩阵 (C) 。

式子的含义为：矩阵 (C_{i,j}) 由矩阵 (A) 第 (i) 行上的 (K) 个数与矩阵 (B) 第 (j) 列上的 (K) 分别相乘并求和得到。

划重点：矩阵乘法不满足交换律，满足结合律，满足分配律（在某些特定情况下满足乘法交换律）。

---

[QAQ ]

---

三.【如何优化加速】

1.【前言】

以 斐波那契数列 ([P1962]) 为例，大部分人都会简单的求 (Fibonacci)，其递推式为 (f(n)=f(n-1)+f(n-2) (n geqslant 2) (且) (n in N^{*}))，其中 (f(1)=f(2)=1) 。

这道题按照正常的递推做法可以过 (60) 分，对于大一点的 (n) 就不行了。

由于其递推方程是固定的，决策点只要两个（(n-1) 和 (n-2)），所以可以考虑用矩阵乘法加速。

---

2.【构造答案矩阵和累乘矩阵】

还是以 斐波那契数列 ([P1962]) 为例，

注意决策点数量：两个，可以先尝试使用二维矩阵，如果不行就试着加一维（辅助递推），大概做法如下：

设 (f(n)) 为 (Fibonacci) 数列第 (n) 项。

先构建一个 (1 imes 2) 的答案矩阵 (F(n))：

(F(n) = egin{vmatrix}f(n)&f(n-1)end{vmatrix})

再构造一个 (2 imes 2) 的累乘矩阵 (base)：

(base = egin{vmatrix}a&b\c&dend{vmatrix})，其中 (a,b,c,d) 均为未知数

我们需要满足：(F(n) imes base = F(n+1)) 。
即 (egin{vmatrix}f(n)&f(n-1)end{vmatrix} imes egin{vmatrix}a&b\c&dend{vmatrix} = egin{vmatrix}af(n)+cf(n-1)&bf(n)+df(n-1)end{vmatrix} = egin{vmatrix}f(n+1)&f(n)end{vmatrix})
即 (af(n)+cf(n-1)=f(n+1),bf(n)+df(n-1)=f(n))
由递推式可知：(a=b=c=1,d=0)
即 (base = egin{vmatrix}1&1\1&0end{vmatrix}) 。

实际上在做题的时候不需要这么麻烦，只需要按这种思路去模拟一下 (base) 就出来了。

---

3.【快速幂加速运算】

对 (F(n)) 稍作转换：(F(n) = F(2) imes base imes base...... imes base)，其中 (base) 一共乘了 (n-2) 次。

为什么要写 (F(2)) 呢？在写题时，这是一个致命的细节问题。
根据定义，(F(1)=egin{vmatrix}f(1)&f(0)end{vmatrix})，其中 (f(0)) 无法计算（或者说毫无意义），所以要用 (f(2)) 作为初始矩阵来计算。

由矩阵乘法结合律可知： (F(n)=F(2) imes base^{n-2} = egin{vmatrix}1&1end{vmatrix} imes egin{vmatrix}1&1\1&0end{vmatrix}^{n-2})

即 (F(n) = egin{vmatrix}1&1end{vmatrix} imes egin{vmatrix}1&1\1&0end{vmatrix}^{n-2} = egin{vmatrix}f(n)&f(n-1)end{vmatrix})

在具体的代码实现中，我们可以将 (F(n)) 视为 (2 imes 2) 的矩阵，多余的部分赋值为 (0)，即 (F(n) = egin{vmatrix}f(n)&f(n-1)\0&0end{vmatrix})，可以发现，不管怎么乘，它还是这样的形式（第二行和 (1 imes 2) 矩阵的变化相同，第二行依然全为 (0)）。
求 (base) 的幂时用一个快速幂，注意快速幂初始值要设为 (F(2))，如果在算完 (base^{n-2}) 后再乘上 (F(2)) 的话，就违背了矩阵乘法不满足交换律的原则。

时间复杂度为 (O(logn))。

---

4.【Code】

//f[1]=f[2]=1,f[n]=f[n-1]+f[n-2]
//[1 1]
//[1 0]
#include<cstring>
#include<cstdio>
#define LL long long
const int P=1e9+7;
struct QAQ{//结构体打包用起来比较方便
    LL a[3][3];
    void CL(){memset(a,0,sizeof(a));a[1][1]=a[2][1]=a[1][2]=1;}
    QAQ operator * (QAQ &x){
        int i,j,k;QAQ ans;memset(ans.a,0,sizeof(ans.a));
        for(i=1;i<3;i++)
            for(j=1;j<3;j++)
                for(k=1;k<3;k++)
                    (ans.a[i][j]+=a[i][k]*x.a[k][j]%P)%=P;
        return ans;
    }
};
LL n;
inline LL sovle(LL k){
    if(k<1)return 1;//特判很重要
    if(!k)return 1;
    QAQ s,x;x.CL();s.a[1][1]=s.a[1][2]=1,s.a[2][1]=s.a[2][2]=0;//初始化F(2)
	while(k){
        if(k&1)s=(s*x);
        x=(x*x);k>>=1;
    }
    return s.a[1][1]%P;
}
int main(){
    scanf("%lld",&n);
    printf("%lld
",sovle(n-2));
}

---

[QAQ ]

---

四.【矩阵加维】

1.【前言】

缺什么补什么，有不确定的信息就先将递推式写出来，然后根据具体情况加维。

下面一共总结了 (3) 种需要加维的情况（可能不全，欢迎补充）：

---

2.【带常数项 k】

递推方程：(f(n)=f(n-1)+f(n-2)+k) 。

求：(f(n))。

常数项不可忽略，应当专门加一维来计算常数。常数项的递推式是最简单的： (k_n=k_{n-1}+0) 。

(F(n) = egin{vmatrix}f(n)&f(n-1)&kend{vmatrix})，(base = egin{vmatrix}1&1&0\1&0&0\1&0&1end{vmatrix}) 。

---

3.【带未知数项 n】

递推方程：(f(n)=f(n-1)+f(n-2)+n) 。

求：(f(n))。

先将未知项的递推式写出来： ((n)=(n-1)+1) ，虽然 (f(n)) 的转移只有四项，但需要加一维来辅助未知项 (n) 的递推。

(F(n) = egin{vmatrix}f(n)&f(n-1)&n&1end{vmatrix})，(base = egin{vmatrix}1&1&0&0\1&0&0&0\1&0&1&0\1&0&1&1end{vmatrix}) 。

---

4.【求和】

递推方程：(f(n)=f(n-1)+f(n-2)) 。

求：(S(n)=sum_{i=1}^{n} f(i)) 。

暴力计算 (n) 次得出 (f(i)) 数组肯定是不行的，但可以尝试将 (S(n)) 放入矩阵跟随着 (f(n)) 一起递推。

先将前缀和的递推式写出来：(S(n)=S(n-1)+f(n)) 。

(F(n) = egin{vmatrix}f(n)&f(n-1)&S(n-1)end{vmatrix})，(base = egin{vmatrix}1&1&1\1&0&0\0&0&1end{vmatrix}) 。

---

[QAQ ]

---

五.【一些经典题目】

1.【连续幂次和】

给定一个 (n imes n) 的矩阵 (A) 和一个整数 (K)，求 (sum_{i=1}^{K} A^i) 。传送门

(1).【两次分治】

对于单个 (A^i) 可以通过 (log) 次转移得到，需要计算多个时就需要再次分治。

原理：(A^1+A^2+A^3......A^n=) (A^1+A^2+A^3...A^{mid}+A^{mid}*(A^1+A^2+A^3...A^{mid})) 。

每次分治时计算一下 (mid)，递归计算 (sum_{i=1}^{mid} A^i)，(log) 次乘法计算 (A^{mid})，另一半可直接得出（当 (K) 为奇数时还需计算 (A^K)）。

时间复杂度：(n^3log^2K) 。

【Code】

#include<iostream>
#include<cstring>
#include<cstdio>
#include<queue>
#define LL long long
#define Re register LL
using namespace std;
const int N=45;
LL n,K,P=10;
inline void in(Re &x){
    Re f=0;x=0;char c=getchar();
    while(c<'0'||c>'9')f|=c=='-',c=getchar();
    while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
    x=f?-x:x;
}
struct Ma{
    LL n,a[N][N];
    Ma(){memset(a,0,sizeof(a));}
    inline void In(){
        for(Re i=1;i<=n;++i)
            for(Re j=1;j<=n;++j)
                in(a[i][j]),a[i][j]%=P;//一定要边读边膜
    }
    inline void Out(){
        for(Re i=1;i<=n;printf("%lld
",a[i][n]),++i)//卡输出。。
            for(Re j=1;j<n;++j)
                printf("%lld ",a[i][j]);
    }
    inline Ma operator*(Ma O)const{
        Ma ans;ans.n=n;
        for(Re i=1;i<=n;++i)
            for(Re j=1;j<=n;++j)
                for(Re k=1;k<=n;++k)
                    (ans.a[i][j]+=a[i][k]*O.a[k][j]%P)%=P;
        return ans;
    }
    inline Ma operator+(Ma O)const{
        Ma ans;ans.n=n;
        for(Re i=1;i<=n;++i)
            for(Re j=1;j<=n;++j)
                (ans.a[i][j]+=(a[i][j]+O.a[i][j])%P)%=P;
        return ans;
    }
    inline void operator+=(Ma O){*this=*this+O;}
    inline void operator*=(Ma O){*this=*this*O;}
}A;
inline Ma mi(Ma x,Re k){
    Ma s=x;--k;//s初始化为一个x，k减1
    while(k){
        if(k&1)s*=x;
        x*=x,k>>=1;
    }
    return s;
}
inline Ma calc(Ma A,Re k){
    if(k==1){return A;}//边界直接返回
    Ma ans,tmp=calc(A,k>>1);//先算出一小段的
    ans=tmp+(tmp*mi(A,k>>1));//算出A^1+A^2+...A^(k/2*2)
    if(k&1)ans+=mi(A,k);//如果k为奇数则再加上一个A^k
    return ans;
}
int main(){
    // freopen("123.txt","r",stdin);
    while(scanf("%lld%lld",&A.n,&K)&&A.n)A.In(),calc(A,K).Out(),puts("");//卡输出。。
}

(2).【倍增】

(To) (be) (continued...)

时间复杂度：(n^3logK) 。

---

2.【有向图中求合法路径方案数】

给出一张 (n) 个点（从 (0) 到 (n-1) 编号） (m) 条边有向图，每次询问求从 (st) 恰好走 (K) 步到达 (ed) 的方案数，重边视作一条路径，每条边可重复走。

(1).【分析】

为方便分析，将点编号都加一，变为 ([1,n]) 。

设 (f_k(i,j)) 表示从点 (i) 到达点 (j) 恰好走 (k) 步的方案数。所以 (f_k(i,j))

当 (K) 较小时可以使用 (bfs+dp)，若存在一条边 (x,j)，则 (f_{k}(i,j)+=f_{k-1}(i,x)*1) 。但如果 (K leqslant 10^9) 就无法解决了。

若用邻接矩阵 (a(i,j)) 来表示点 (i) 与 (j) 之间是否连边，那么根据上述转移式子可得： (f_{k}(i,j)=sum_{x=1}^{n}f_{k-1}(i,x)*a(x,j)) 。

由于询问是不定向的，起点和终点可能是 ([1,n]) 中的任意一个，所以答案矩阵应包含 (n^2) 个信息，其中 (F(k)) 中的第 (i) 行第 (j) 列表示 (f_k(i,j)) ，即用 (F(K)) 表示恰好走 (K) 步的答案矩阵。

所以对于每条边 ((x,y))，使累乘矩阵中第 (x) 行第 (y) 列的数为 (1)，最后直接求一个 (K) 次幂即可。
(F(k) = egin{vmatrix}f_k(1,1)&f_k(1,2)&...&f_k(1,n)\f_k(2,1)&f_k(2,2)&...&f_k(2,n)\...&...&...&...\f_k(n,1)&f_k(n,2)&...&f_k(n,n)end{vmatrix})

(base) 略。

时间复杂度：(O(n^3logK)) 。

【Code】

#include<cstring>
#include<cstdio>
#define Re register int
using namespace std;
const int N=23;
int x,y,n,m,T,K,P=1000;
inline void in(Re &x){
    Re f=0;x=0;char c=getchar();
    while(c<'0'||c>'9')f|=c=='-',c=getchar();
    while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
    x=f?-x:x;
}
struct Ma{
    int n,a[N][N];
    Ma(){memset(a,0,sizeof(a));}
    inline Ma operator*(Ma O)const{
        Ma ans;ans.n=n;Re i,j,k;
        for(i=1;i<=n;++i)
            for(j=1;j<=n;++j)
                for(k=1;k<=n;++k)
                    (ans.a[i][j]+=a[i][k]*O.a[k][j]%P)%=P;
        return ans;
    }
    inline void operator*=(Ma O){*this=*this*O;}
};
inline Ma mi(Ma x,Re k){
    Ma s=x;--k;
    while(k){
        if(k&1)s*=x;
        x*=x,k>>=1;
    }
    return s;
}
int main(){
    // freopen("123.txt","r",stdin);
    while(~scanf("%d%d",&n,&m)&&(n||m)){
        Ma A;A.n=n;
        while(m--)in(x),in(y),A.a[x+1][y+1]=1;
        in(T);
        while(T--){
            in(x),in(y),in(K);
            printf("%d
",K?mi(A,K).a[x+1][y+1]:x==y);//注意：K=0时要特判，不需要走动时输出1否则输出0
        }
    }
}

(2).【扩展 1】

询问改为：求走的步数不超过 (K) 的方案数。其他条件不变。

求一个 (1) 至 (K) 的连续幂次和即可。

设 (S_K(i,j)=sum_{k=1}^{K}f_k(i,j)) 。

前缀和 (S_K(i,j)) 也要记录 (n^2) 个，要尽量缩减矩阵规模的话，可以把他们一层层地包裹在原 (n imes n) 的矩阵外面（自己口胡的，没有尝试过），但代码难度较大，也可以将原矩阵扩大一倍变成 (2n imes 2n)，代码难度较小，即：(F(k) = egin{vmatrix}f_k(1,1)&...&f_k(1,n)&S_k(1,1)&...&S_k(1,n)\...&...&...&...&...&...\f_k(n,1)&...&f_k(n,n)&S_k(n,1)&...&S_k(n,n)\end{vmatrix}) 。

(base) 略。

时间复杂度：(O((2n)^3logK))

(3).【扩展 2】

增加一个限制条件：一共有四种物品，每条边上有相同或不同种类的若干个，每次经过时如果拿走这些物品则会多消耗一步（相当于走两步），求走的步数不超过 (K) 且能将四种物品都拿全的方案数。

先考虑走两步的转移，在答案矩阵中再加入 (n^2) 个 (f_{k-1}) 的信息，其中一种加维方案（空白处填 (0)）：

(F(k) = egin{vmatrix}f_k(1,1)&...&f_k(1,n)&S_k(1,1)&...&S_k(1,n)\...&...&...&...&...&...\f_k(n,1)&...&f_k(n,n)&S_k(n,1)&...&S_k(n,n)\f_k(1,1)&...&f_k(1,n)\...&...&...\f_k(n,1)&...&f_k(n,n)end{vmatrix})

(base) 略。

再考虑物品限制，假设要求最后只能选 (1,2)，那么在构造矩阵时带有 (3,4) 的边都不加进去。

如上所述，跑若干次计算搞一下容斥即可（我太蒻了，容斥这里只有先咕着了）。

(To) (be) (continued...)

---

3.【？】

(To) (be) (continued...)

---

[QAQ ]

---

六.【题目链接】

【简单题】

• 斐波那契数列 ([P1962])
  【标签】递推/矩阵乘法

• 【模板】矩阵加速（数列）([P1939])
  【标签】递推/矩阵乘法

【中档题】

• 小奇的集合 ([BZOJ4547]) ([HDU5171])
  【标签】递推/矩阵乘法
  【题解】(hzwer)

• (Power) (of) (Matrix) ([UVA11149])
  (Matrix) (Power) (Series) ([Poj3233])
  【标签】递推/矩阵乘法/二分

• (How) (many) (ways??) ([Hdu2157])
  【标签】递推/矩阵乘法

• 牛继电器 (Cow) (Relays) ([P2886]) ([Poj3613])
  【标签】递推/矩阵乘法

---

[QAQ ]

---

七.【参考资料】

• 题解 ([P1962]) 斐波那契数列

• 矩阵的十大经典题目，留份做题