【學習筆記】動態規划—矩陣遞推加速
【大前言】
矩陣優化 \(dp\) 通常用於線性遞推式的 \(dp\) 優化，能以優異的時間復雜度實現大量的狀態轉移。
更完整的 \(dp\) 優化策略：【學習筆記】動態規划—各種 \(DP\) 優化
一.【題目特征】
\((1).\) 類似線性遞推（划重點，包括有向圖上的遞推等等）
\((2).\) 轉移次數 \(10^9\) 左右（霧）
\((3).\) 決策點較少（常數）
二.【前置芝士】
1.【前言】
首先要清楚矩陣是個什么東西，在對 \(dp\) 進行優化時通常只會用到矩陣乘法和矩陣加法，其中矩陣乘法最為常見。
2.【運算法則】
(1).【矩陣加法】
矩陣加法的一般式：\(C_{i,j}=A_{i,j} + B_{i,j}\)，其中 \(A,B\) 均為 \(N \times N\) 的矩陣，矩陣 \(A + B\) 得到 \(N \times N\) 的矩陣 \(C\) 。
式子的含義為：矩陣 \(C\) 由矩陣 \(A,B\) 對應位置上數值相加所得。
划重點：矩陣加法滿足交換律 。
(2).【矩陣乘法】
矩陣乘法的一般式：\(C_{i,j}=\sum_{k=1}^{K} (A_{i,k} \times B_{k,j})\)，其中 \(A\) 為 \(N \times K\) 的矩陣，\(B\) 為 \(K*M\) 的矩陣，矩陣 \(A \times B\) 得到 \(N \times M\) 的矩陣 \(C\) 。
式子的含義為：矩陣 \(C_{i,j}\) 由矩陣 \(A\) 第 \(i\) 行上的 \(K\) 個數與矩陣 \(B\) 第 \(j\) 列上的 \(K\) 分別相乘並求和得到。
划重點：矩陣乘法不滿足交換律，滿足結合律，滿足分配律（在某些特定情況下滿足乘法交換律）。
三.【如何優化加速】
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 \times 2\) 的答案矩陣 \(F(n)\)：
\(F(n) = \begin{vmatrix}f(n)&f(n-1)\end{vmatrix}\)
再構造一個 \(2 \times 2\) 的累乘矩陣 \(base\)：
\(base = \begin{vmatrix}a&b\\c&d\end{vmatrix}\)，其中 \(a,b,c,d\) 均為未知數
我們需要滿足：\(F(n) \times base = F(n+1)\) 。
即 \(\begin{vmatrix}f(n)&f(n-1)\end{vmatrix} \times \begin{vmatrix}a&b\\c&d\end{vmatrix} = \begin{vmatrix}af(n)+cf(n-1)&bf(n)+df(n-1)\end{vmatrix} = \begin{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 = \begin{vmatrix}1&1\\1&0\end{vmatrix}\) 。
實際上在做題的時候不需要這么麻煩，只需要按這種思路去模擬一下 \(base\) 就出來了。
3.【快速冪加速運算】
對 \(F(n)\) 稍作轉換：\(F(n) = F(2) \times base \times base......\times base\)，其中 \(base\) 一共乘了 \(n-2\) 次。
為什么要寫 \(F(2)\) 呢？在寫題時，這是一個致命的細節問題。
根據定義，\(F(1)=\begin{vmatrix}f(1)&f(0)\end{vmatrix}\)，其中 \(f(0)\) 無法計算（或者說毫無意義），所以要用 \(f(2)\) 作為初始矩陣來計算。
由矩陣乘法結合律可知： \(F(n)=F(2) \times base^{n-2} = \begin{vmatrix}1&1\end{vmatrix} \times \begin{vmatrix}1&1\\1&0\end{vmatrix}^{n-2}\)
即 \(F(n) = \begin{vmatrix}1&1\end{vmatrix} \times \begin{vmatrix}1&1\\1&0\end{vmatrix}^{n-2} = \begin{vmatrix}f(n)&f(n-1)\end{vmatrix}\)
在具體的代碼實現中，我們可以將 \(F(n)\) 視為 \(2 \times 2\) 的矩陣，多余的部分賦值為 \(0\)，即 \(F(n) = \begin{vmatrix}f(n)&f(n-1)\\0&0\end{vmatrix}\)，可以發現，不管怎么乘，它還是這樣的形式（第二行和 \(1 \times 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\n",sovle(n-2));
}
四.【矩陣加維】
1.【前言】
缺什么補什么，有不確定的信息就先將遞推式寫出來，然后根據具體情況加維。
下面一共總結了 \(3\) 種需要加維的情況（可能不全，歡迎補充）：
2.【帶常數項 k】
遞推方程：\(f(n)=f(n-1)+f(n-2)+k\) 。
求：\(f(n)\)。
常數項不可忽略，應當專門加一維來計算常數。常數項的遞推式是最簡單的： \(k_n=k_{n-1}+0\) 。
\(F(n) = \begin{vmatrix}f(n)&f(n-1)&k\end{vmatrix}\)，\(base = \begin{vmatrix}1&1&0\\1&0&0\\1&0&1\end{vmatrix}\) 。
3.【帶未知數項 n】
遞推方程：\(f(n)=f(n-1)+f(n-2)+n\) 。
求：\(f(n)\)。
先將未知項的遞推式寫出來： \((n)=(n-1)+1\) ，雖然 \(f(n)\) 的轉移只有四項，但需要加一維來輔助未知項 \(n\) 的遞推。
\(F(n) = \begin{vmatrix}f(n)&f(n-1)&n&1\end{vmatrix}\)，\(base = \begin{vmatrix}1&1&0&0\\1&0&0&0\\1&0&1&0\\1&0&1&1\end{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) = \begin{vmatrix}f(n)&f(n-1)&S(n-1)\end{vmatrix}\)，\(base = \begin{vmatrix}1&1&1\\1&0&0\\0&0&1\end{vmatrix}\) 。
五.【一些經典題目】
1.【連續冪次和】
給定一個 \(n \times 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\n",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) = \begin{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\n",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 \times n\) 的矩陣外面（自己口胡的，沒有嘗試過），但代碼難度較大，也可以將原矩陣擴大一倍變成 \(2n \times 2n\)，代碼難度較小，即：\(F(k) = \begin{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) = \begin{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...\)
六.【題目鏈接】
【簡單題】
-
斐波那契數列 \([P1962]\)
【標簽】遞推/矩陣乘法 -
【模板】矩陣加速（數列）\([P1939]\)
【標簽】遞推/矩陣乘法
【中檔題】
-
小奇的集合 \([BZOJ4547]\) \([HDU5171]\)
【標簽】遞推/矩陣乘法
【題解】\(hzwer\) -
\(Power\) \(of\) \(Matrix\) \([UVA11149]\)
\(Matrix\) \(Power\) \(Series\) \([Poj3233]\)
【標簽】遞推/矩陣乘法/二分 -
\(How\) \(many\) \(ways??\) \([Hdu2157]\)
【標簽】遞推/矩陣乘法 -
牛繼電器 \(Cow\) \(Relays\) \([P2886]\) \([Poj3613]\)
【標簽】遞推/矩陣乘法