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

xAlec · 2024-12-18 20:11:23 · 算法·理论

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 & 3 & 7 \\ 16 & 3 & 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} f_{n - 1} \\ f_{n - 2} \end{bmatrix} \cdot \begin{bmatrix} 1 & 1 \\ 1 & 0 \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，将矩阵乘法转换成 Floyd 的转移，答案即为 e^N_{s,e}，注意处理重边。

#include <bits/stdc++.h>
// #define int long long

using 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(n^3 + q)。

#include <bits/stdc++.h>
#define int unsigned int

using 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;
}

Summary

其实只要找到转移的式子，列出矩阵，其它的就可以自行实现了，尤其是 dp 方程，在列出暴力转移的情况下可以考虑矩阵优化，同时部分数据结构也可以维护矩阵。

trick：由于图论中用矩阵快速幂处理邻接矩阵需要保证边权 w \in \{0,1\}，因此在边权的范围较小可以考虑拆点分层图。

Others

P1306 斐波那契公约数重点是斐波那契相关定理的证明。

P2151 [SDOI2009] HH去散步点转边的 trick。

P4159 [SCOI2009] 迷路分层图即拆点。

P6569 [NOI Online #3 提高组] 魔法值矩乘变形成异或和。

UVA Power of Matrix 矩阵套矩阵。

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

Basic Tips

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

矩阵数乘 / 矩阵加减

矩阵乘法 / 矩阵快速幂

矩阵的运算律