【学习笔记】矩阵乘法 / 矩阵快速幂的应用-编程知识

【学习笔记】矩阵乘法 / 矩阵快速幂的应用

news/2024/12/23 7:55:28/文章来源:https://www.cnblogs.com/Ayaka-0928/p/18622973

Basic Tips

不同角度看待线性方程组：矩阵

对于一个方程组：
$
\begin{cases}
4x_1 + 16x_2 = 114 \\
3x_1 + 3x_2 = 514 \\
7x_1 + 3x_2 = 416
\end{cases}
$

利用一种打包处理的思想将其写成矩阵的形式：
$\begin{bmatrix}x_1 & x_2\end{bmatrix} \cdot \begin{bmatrix}4 & 16 \\\\ 3 & 3 \\\\ 7 & 3 \end{bmatrix} = \begin{bmatrix}114 \\\\ 514 \\\\ 416 \end{bmatrix}$

矩阵的存储一般使用二维数组。
矩阵数乘 / 矩阵加减

对于一个矩阵 $A$，让它乘上一个数 $x$，结果便是 $A$ 中每个位置的元素乘上 $x$。

矩阵加减只能在大小相同的两个矩阵之间进行，例如大小相同两个矩阵 $A$ 跟 $B$ 相加 / 减，得到的结果矩阵即为两矩阵对应位置上的元素相加 / 减。
矩阵乘法 / 矩阵快速幂

对于大小为 $X \times K$ 的矩阵 $A$ 和大小为 $K \times Y$ 的矩阵 $B$，记它们的乘积矩阵为 $C$，可得 $C_{i,j} = \displaystyle\sum^{K}_{k=1}{A_{i,k} \times B_{k,j}}$，注意，矩阵乘法并不满足 交换律，矩阵乘法可以快速的处理一些 $ax + by = c$ 的关系，需要关注枚举顺序，矩阵快速幂通常用来加速递推、优化 dp 转移，在图论中根据邻接矩阵可以求限制路径长度的路径方案数，复杂度取决于转移矩阵的大小，一般是 $O(\log n)$ 带常数。
矩阵的运算律

矩阵并不满足一般的交换律，即 $A \times B$ 不等同于 $B \times A$，结果可能不一样，矩阵满足结合和左右分配律，即 $A \cdot (BC) = (AB) \cdot C$。

矩阵快速幂的实现

相较于数的快速幂，矩阵快速幂要初始化为单位矩阵。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod = 1e9 + 7;
int n,k;struct node{int mat[105][105];node(){memset(mat,0,sizeof(mat));}inline void init(){for(int i = 1;i <= n;i ++) mat[i][i] = 1;}
}x; // 结构体分装会更方便node operator *(const node &x,const node &y)
{node z;for(int k = 1;k <= n;k ++){for(int i = 1;i <= n;i ++){for(int j = 1;j <= n;j ++)z.mat[i][j] = (z.mat[i][j] + x.mat[i][k] * y.mat[k][j] % mod) % mod;}}return z;
} // 重定义乘号运算inline node qpow(node x,int k)
{node ret;ret.init();while(k){if(k & 1) ret = ret * x;x = x * x;k >>= 1;}return ret;
} // 初始化单位矩阵signed main()
{scanf("%lld %lld",&n,&k);for(int i = 1;i <= n;i ++){for(int j = 1;j <= n;j ++)scanf("%lld",&x.mat[i][j]);}node Ans = qpow(x,k);for(int i = 1;i <= n;i ++){for(int j = 1;j <= n;j ++)printf("%lld%c",Ans.mat[i][j]," \n"[j == n]);}return 0;
}

What's More

单位矩阵 $I$：除对角线元素为 $1$ 外其余元素皆为 $0$。

小细节：实现矩阵比较方便的是使用结构体分装，通过构造函数进行初始化操作，部分存在结合律的运算也可以通过变形的矩阵乘法实现，所以对于题目中出现的 $01$ 异或、gcd、min、max 等均可以往这方面思考。

Examples

P1962 斐波那契数列

对于 $n \leq 2^{63}$ 的数据，正常的递推肯定会爆，因此我们来观察转移式子：

$ \large \begin{aligned} & f_n = 1 \times f_{n - 1} + 1 \times f_{n - 2} \\ & f_{n - 1} = 1 \times f_{n - 1} + 0 \times f_{n - 2} \end{aligned} $

这就是上文所提到的线性方程组的形式，写成矩阵：

$ \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} $
$\cdot$
$ \begin{bmatrix} f_{n - 1} \\ f_{n - 2} \end{bmatrix} $
$=$
$ \begin{bmatrix} f_n \\ f_{n - 1} \end{bmatrix} $

求 $f_n$ 相当于转移矩阵 $A^n$ 乘上初始矩阵。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod = 1e9 + 7;
struct node{int mz[3][3];node operator * (const node& x) const{node ret;for(int i = 1;i <= 2;i ++){for(int j = 1;j <= 2;j ++){ret.mz[i][j] = 0;for(int k = 1;k <= 2;k ++)ret.mz[i][j] = ((ret.mz[i][j] + mz[i][k] * x.mz[k][j] % mod) % mod + mod) % mod;}}return ret;}
}feb;
int n;
inline node qpow(node x,int p)
{node res = x;p --;while(p){if(p & 1)res = res * x;x = x * x;p >>= 1;}return res;
}
signed main()
{cin.tie(0) -> sync_with_stdio(0);cout.tie(0) -> sync_with_stdio(0);cin >> n;feb.mz[1][1] = feb.mz[1][2] = feb.mz[2][1] = 1;if(n <= 2)cout << "1" << endl;else{node Answer = qpow(feb,n);cout << Answer.mz[1][2] << endl;}return 0;
}

P2886 [USACO07NOV] Cow Relays G

图论上的应用，鉴于 $T \leq 100$，可以邻接矩阵存图，由于是求最短路，可以用 Floyd，发现可以类比矩阵乘法，记邻接矩阵为 $e$，将矩阵乘法转换成 $\text{Floyd}$ 的转移，答案即为 $e^N_{s,e}$，注意处理重边。

#include <bits/stdc++.h>
// #define int long longusing namespace std;
template <typename Tp>
inline void read(Tp &x){Tp res = 0,f = 1;char ch = getchar();while(!isdigit(ch)){if(ch == '-') f = -1;ch = getchar();}while(isdigit(ch)){res = (res << 1) + (res << 3) + (ch ^ 48);ch = getchar();}x = res * f;
}const int maxn = 1005,inf = 0x3f3f3f3f;
int n,t,s,e,apr[maxn],idx;struct Matrix{int mat[maxn][maxn];Matrix(){memset(mat,inf,sizeof(mat));}Matrix operator * (const Matrix &g) const{Matrix ret;for(int k = 1;k <= idx;k ++){for(int i = 1;i <= idx;i ++){for(int j = 1;j <= idx;j ++){ret.mat[i][j] = min(ret.mat[i][j],g.mat[i][k] + mat[k][j]);}}}return ret;}
}Ans;Matrix mat_qpow(Matrix x,int k){Matrix ret;ret = x,k --;while(k){if(k & 1) ret = ret * x;x = x * x;k >>= 1;}return ret;
}signed main(){read(n),read(t),read(s),read(e);for(int i = 1;i <= t;i ++){int w,u,v;read(w),read(u),read(v);if(!apr[u]) apr[u] = ++ idx;if(!apr[v]) apr[v] = ++ idx;Ans.mat[apr[u]][apr[v]] = Ans.mat[apr[v]][apr[u]] = w; // 处理重边}Ans = mat_qpow(Ans,n);printf("%d\n",Ans.mat[apr[s]][apr[e]]);return 0;
}

P6569 [NOI Online #3 提高组] 魔法值

由于原图的边权均为 $0$ 或 $1$，因此题面中的式子可以变成：

$ \begin{aligned} f_{x,i} = f_{1,i - 1} \times e_{1,x} \oplus f_{2,i - 1} \times e_{2,x} \oplus \dots f_{n,i - 1} \times e_{n,x} \end{aligned} $

啊，这不就是异或版矩乘？

记第 $i$ 天的魔法值矩阵为 $f_i$，邻接矩阵为 $e$，由于邻接矩阵有且仅有 $01$，因此它满足矩乘结合律，得到以下式子：

$ \begin{aligned} f_{i} = f_{i - 1} \times e \to f_i = f_0 \times e^i \end{aligned} $

如果每次询问都跑一次快速幂一定会爆，因此看到数据范围 $f_i,a_i < 2^{32}$，不难想到类似状压的方法，预处理邻接矩阵的 $2^i$ 即可，复杂度 $O(31(n^3 + q))$

#include <bits/stdc++.h>
#define int unsigned intusing namespace std;
template <typename Tp>
inline void read(Tp &x){Tp res = 0,f = 1;char ch = getchar();while(!isdigit(ch)){if(ch == '-') f = -1;ch = getchar();}while(isdigit(ch)){res = (res << 1) + (res << 3) + (ch ^ 48);ch = getchar();}x = res * f;
}const int maxn = 105;
int n,m,q,f0[maxn],x;struct Matrix{int mat[maxn][maxn];Matrix(int x = 0){memset(mat,0,sizeof(mat));for(int i = 1;i <= n;i ++) mat[i][i] = x;}Matrix operator * (const Matrix &g) const{Matrix ret(0);for(register int k = 1;k <= n;k ++){for(register int i = 1;i <= n;i ++){for(register int j = 1;j <= n;j ++){ret.mat[i][j] ^= mat[i][k] * g.mat[k][j];}}}return ret;}
};Matrix d[33];signed main(){read(n),read(m),read(q);for(register int i = 1;i <= n;i ++) read(f0[i]);for(register int i = 1;i <= m;i ++){int u,v;read(u),read(v);d[0].mat[u][v] = d[0].mat[v][u] = 1;}for(register int i = 1;i < 32;i ++)d[i] = d[i - 1] * d[i - 1];while(q --){read(x);Matrix tmp;for(register int i = 1;i <= n;i ++)tmp.mat[1][i] = f0[i];for(register int i = 0;i < 32;i ++){if((x >> i) & 1) tmp = tmp * d[i];}printf("%u\n",tmp.mat[1][1]);}return 0;
}