矩阵加速学习笔记

发表于 更新于 分类于 阅读次数： 本文字数： 8.4k
矩阵用于加速递推可以配合矩阵快速幂，把复杂度降到 log 级别，非常适合用来加速计算。

一些题目中，我们会遇到一些类似于 \(f_i = f_{i-1} + k\) 的递推式，如果我们求第 \(n\) 项，当 \(n\) 不大时暴力计算确实可行，但是一旦 \(n\) 达到类似 \(10^{18}\) 的范围时，我们是一定需要优化的，通过矩阵加速递推，我们可以大大降低复杂度，在一些 dp 题非常好用。

下面以一些例题来介绍矩阵加速：

矩阵加速递推

例题

P1962 斐波那契数列

题意

给定 \(f_1 = f_2 = 1\)，求满足 \(f_{i+2} = f_{i+1} + f_i\) 的第 \(n\) 项 模 \(10^{9} + 7\) 的值。\(n< 2^{63}\)。

题解

由于\(\begin{bmatrix}f_{i+1} & f_{i}\end{bmatrix}\begin{bmatrix}1 & 1\\1& 0\end{bmatrix} = \begin{bmatrix}f_{i+2}&f_i\end{bmatrix}\)，我们求的 \(f_{n} = \begin{bmatrix}1&1\end{bmatrix}\begin{bmatrix}1&1\\1&0\end{bmatrix}^{n-2}(n\ge3)\)，矩阵的幂采用快速幂处理，那么时间复杂度为 \(O(2^3\log n)\)。这里需要注意的是，每乘一次，我们求的值下标就会加一，但是 \(f_1,f_2\) 是已知的，因此乘一次矩阵得到的就是 \(f_3\) 了，\(f_n\) 只需要乘 \(n - 2\)次。

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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <bits/stdc++.h>
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;
using u128 = unsigned __int128;

constexpr int mod = 1e9 + 7;
struct Matrix {
    vector<vector<i64>> a;
    Matrix() { a.assign(3, vector<i64>(3)); }
    Matrix operator*(const Matrix b) const {
        Matrix res;
        for (int i = 1; i <= 2; i++) {
            for (int j = 1; j <= 2; j++) {
                for (int k = 1; k <= 2; k++) {
                    res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
                }
            }
        }
        return res;
    }
    void inita() { a[1][1] = a[1][2] = 1; }
    void initb() { a[1][1] = a[2][1] = a[1][2] = 1; }
};

Matrix qpow(Matrix a, i64 b) {
    Matrix res;
    for (int i = 1; i <= 2; i++) {
        res.a[i][i] = 1;
    }
    while (b) {
        if (b % 2 == 1) {
            res = res * a;
        }
        a = a * a;
        b /= 2;
    }
    return res;
}

void solve() {
    i64 n;
    cin >> n;
    if (n <= 2) {
        cout << 1 << "\n";
        return;
    }
    Matrix a, b;
    a.inita(), b.initb();
    a = a * qpow(b, n - 2);
    cout << a.a[1][1] << "\n";
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    int t = 1;
    // cin >> t;
    while (t--) {
        solve();
    }
    return 0;
}

P1707 刷题比赛

题意

在第 \(1\) 天 nodgd,Coicoi,Nicole 都只做了 \(1\) 道题。

在第 \(2\) 天 nodgd,Coicoi,Nicole 都只做了 \(3\) 道题。

1、nodgd 同学第 \(k+2\) 天刷题数量
\[a_{k+2}=pa_{k+1}+qa_k+b_{k+1}+c_{k+1}+rk^2+tk+1\]

2、Ciocio 同学第 \(k+2\) 天刷题数量
\[b_{k+2}=ub_{k+1}+vb_k+a_{k+1}+c_{k+1}+w^k\]

3、Nicole 同学第 \(k+2\) 天刷题数量
\[c_{k+2} = xc_{k+1}+yc_k + a_{k+1} + b_{k+1} + z^k+k+2\]

（以上的字母 \(p,q,r,t,u,v,w,x,y,z\) 都是给定的常数，并保证是正整数）

求第 \(n\) 天每个人的刷题数量。 （\(n\le10^{16}\)） 由于结果很大，输出结果 \(\bmod \space m\) 的值即可。

题解

是一个递推关系式，我们首先确定每一轮会变化的值：\(a_{k+1},b_{k+1},c_{k+1},a_k,b_k,c_k,k,k^2,w^k,z^k\)，因此我们的答案矩阵需要包括它们，再考虑每个数从上一轮的转移，\(a,b,c\) 按照公式即可，但是注意还有一个常数 1，我们也加入矩阵中。

比较麻烦的是 \(k,k^2,w^k,z^k\)，我们比较两次的结果：

\((k+1) = (k) + 1\)，因此通过 \(1\) 和 \(k\) 转移

\((k+1)^2 = (k^2) + 2(k) + 1\)，通过 \(k^2,k,1\)，需要注意 \(k\) 的系数

\((w^{k+1}) = (w^k)w\)，通过 \(w^k\) 转移，系数为 \(w\)。\(z^k\) 同理。

那么，我们就可以构造出转移矩阵：

image-20250328165954718

还需要注意的一点是，这道题的模数可能比较大，相乘会爆 long long，因此进行乘法时可以采用类似快速幂的方法，把乘法拆成加法运算（龟速乘），时间复杂度增加一个 log，也可以使用 \(O(1)\) 的方法：https://oi-wiki.org/math/binary-exponentiation/#%E6%A8%A1%E6%84%8F%E4%B9%89%E4%B8%8B%E5%A4%A7%E6%95%B4%E6%95%B0%E4%B9%98%E6%B3%95。

以下代码复杂度为 \(O(11^3\log n\log m)\)

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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include <bits/stdc++.h>
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;
using u128 = unsigned __int128;

i64 mod;
i64 p, q, r, t, u, v, w, x, y, z, b, c;

i64 qmul(i64 a, i64 b) {
    i64 res = 0;
    while (b) {
        if (b % 2 == 1) {
            res = (res + a) % mod;
        }
        a = (a + a) % mod;
        b /= 2;
    }
    return res;
}

// i64 qmul(i64 a, i64 b) {
//     unsigned long long c =
//         (unsigned long long)a * b -
//         (unsigned long long)((long double)a / mod * b + 0.5L) * mod;
//     if (c < mod) return c;
//     return c + mod;
// }

struct Matrix {
    vector<vector<i64>> a;
    Matrix() { a.assign(12, vector<i64>(12)); }
    Matrix operator*(const Matrix b) const {
        Matrix res;
        for (int i = 1; i <= 11; i++) {
            for (int j = 1; j <= 11; j++) {
                for (int k = 1; k <= 11; k++) {
                    res.a[i][j] =
                        (res.a[i][j] + qmul(a[i][k], b.a[k][j])) % mod;
                }
            }
        }
        return res;
    }
    void inita() {
        a[1][1] = a[1][2] = a[1][3] = 3, a[1][4] = a[1][5] = a[1][6] = 1;
        a[1][7] = a[1][8] = 1;
        a[1][9] = w;
        a[1][10] = z;
        a[1][11] = 1;
    }
    void initb() {
        a[1][1] = p, a[2][1] = 1, a[3][1] = 1, a[4][1] = q, a[7][1] = r,
        a[8][1] = t, a[11][1] = 1;
        a[1][2] = 1, a[2][2] = u, a[3][2] = 1, a[5][2] = v, a[9][2] = 1;
        a[1][3] = 1, a[2][3] = 1, a[3][3] = x, a[6][3] = y, a[8][3] = 1,
        a[10][3] = 1, a[11][3] = 2;
        a[1][4] = 1, a[2][5] = 1, a[3][6] = 1;
        a[7][7] = a[11][7] = 1;
        a[8][7] = 2;
        a[8][8] = a[11][8] = 1;
        a[9][9] = w;
        a[10][10] = z;
        a[11][11] = 1;
    }
};

Matrix qpow(Matrix a, i64 b) {
    Matrix res;
    for (int i = 1; i <= 11; i++) {
        res.a[i][i] = 1;
    }
    while (b) {
        if (b % 2 == 1) {
            res = res * a;
        }
        a = a * a;
        b /= 2;
    }
    return res;
}

void solve() {
    i64 n;
    cin >> n >> mod;
    cin >> p >> q >> r >> t >> u >> v >> w >> x >> y >> z;
    if (n == 1) {
        cout << "nodgd 1" << "\n";
        cout << "Ciocio 1" << "\n";
        cout << "Nicole 1" << "\n";
        return;
    }
    Matrix a, b;
    a.inita();
    b.initb();
    b = qpow(b, n - 2);
    a = a * b;
    cout << "nodgd " << a.a[1][1] << "\n";
    cout << "Ciocio " << a.a[1][2] << "\n";
    cout << "Nicole " << a.a[1][3] << "\n";
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    int t = 1;
    // cin >> t;
    while (t--) {
        solve();
    }
    return 0;
}

图上 DP 优化

定长路径统计

考虑一个经典问题：给定一个图，边权均为 0 或 1，求从起点到终点的最短路径。

显然可以通过 Dijkstra 来求最短路，时间复杂度 \(O(m\log m)\)，同时由于特殊的边权（0 或 1），可以采用 0-1 BFS 的思想，通过双端队列，遇到 0 加到队首，遇到 1 加到队尾来跑一遍 BFS，可以优化到 \(O(m)\)。

将这个问题扩展到计数上：给定一个图，边权要么为 1，要么不存在，求从起点出发，经过总长度为 \(k\) 的路径到达终点的方案数。（\(k\) 的范围很大，而 \(n\le 500\)）

如果没有括号内范围限制，那么自然可以设 \(dp_{i,j}\) 表示走到 \(i\) ，路径长度为 \(j\) 的方案数，显然有复杂度为 \(O(mk)\) 的做法。

但是由于 \(k\) 的大小限制，容易想到需要矩阵乘法来优化，每走一步就是乘一次矩阵，那么 \(k\) 步其实就是初始矩阵乘上距离矩阵的 \(k\) 次方（初始矩阵代表路径长度为 0 ，只有起点有值），通过矩阵快速幂，可以优化到 \(O(500^3\log k)\)​。如果起点和终点不固定，那么我们的初始化和统计答案会有略微区别，中间的思想是一致的。

定长最短路

现在，我们需要对于每个点对 \((u,v)\)，找到从 \(u\) 到 \(v\) 的恰好包含 \(k\) 条边的最短路的长度。（不一定是简单路径，即路径上的点或者边可能走多次），类似地，我们仍然通过邻接矩阵的形式进行状态转移，已知包含 \(k\) 条边的最短路，那么在乘一次矩阵就是包含 \(k+1\) 条边的最短路。但是这里的 “乘法” 有一定区别，因为是求最短路，所以答案应该是取 \(\min(L_{k}[x,p] + L_1[p,y])\)，这样才表示 \(x\) 到 \(y\) 的最短路，复杂度是相同的，为 \(O(n^3\log k)\)。

限长路径统计

问题：给定一个图，边权要么为 1，要么不存在，求从起点出发，经过总长度小于或等于 \(k\) 的路径到达终点的方案数。（\(k\) 的范围很大，而 \(n\le 500\)）

这里和定长路径统计的区别就是：这里的长度是小于或等于，也就是说如果某条路径提前到达了目标点，只要他能通过别的路径在 \(k\) 步之内走回来，每次回到目标点又能算一次新的路径。

由于 \(k\) 极大，我们不可能将长度为 \(0\) 到 \(k\) 的路径加起来，必须另想办法。我们使用拆点的思想，对于每个点 \(v\)，都建立一个虚点 $ v$ 用来统计答案，加入 \((v,v\prime)\) 和 \((v\prime,v\prime)\) 两条边，求起点到 \(v\prime\) 长度为 \(k+1\) 的路径数量即可。这样做的原因是，对于提前到达 \(v\prime\) 的点，由于存在自环，它会不停地走直到长度为 \(k+1\)，而能在 \(k\) 步到达 \(v\) 的路径，都能在 \(k+1\) 步到达 \(v\prime\) 点。所以，我们将限长路径统计转化为了定长路径统计，不过由于每个点拆成两个， \(n\) 的数量翻倍，常数扩大为 8 倍。

限长最短路

对于每个点对 \((u,v)\)，找到从 \(u\) 到 \(v\) 的边数小于或等于 \(k\) 的最短路的长度。

同样的思想，我们对每个点都加入长度为 0 的自环，按照定长最短路的方法求恰好包含 \(k\) 条边的最短路的长度即可。

例题

[SCOI2009] 迷路

题意

给定一个有向图，边权为 1 到 9，问从起点 1 到终点 \(n\) ，长度为 \(k\)​ 的路径数。

\(n\le 10, k\le10^9\)。

题解

与上面四种基本情况的差别在于，本题的边权不一定是 1，我们希望转化为 1 的情况。

一种思路是，对于每一条边，既然边权为 \(1\) 到 \(9\)，我们可以把这条边新增 边权 - 1 个点，然后内部的点单向连接，结尾的点连向目标点。但是，这样做的话，考虑一种最坏情况：每一条边的边权都是 9，那么每个点都向其他 \(n-1\) 个点连边，边上有 8 个点，一共的点的数量就是 \(n + 8n(n-1)\)，总点数接近 1000，此时使用矩阵快速幂，在 \(O(1000^3\log t)\) 的复杂度是不能通过的。

上面是拆边的方法，我们还可以尝试拆点：把每个点拆成 9 个点，形成一条链的形状，分别编号为 \(1,2...,9\)，一旦原来的点向其他点连一条边权为 \(x\) 的边，就把编号为 \(x\) 的边连向目标点的起点。这样形成的图与原来的是等价的，但点数最大只有 90，是可以在 \(O(90^3\log t)\) 的复杂度通过的。注意边权为 0 不能连边！！！

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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#include <bits/stdc++.h>
using namespace std;

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;
using u128 = unsigned __int128;

constexpr int mod = 2009;

struct Matrix {
    vector<vector<i64>> a;
    Matrix() { a.assign(91, vector<i64>(91)); }
    Matrix operator*(const Matrix& b) const {
        Matrix res;
        for (int i = 1; i <= 90; i++) {
            for (int j = 1; j <= 90; j++) {
                for (int k = 1; k <= 90; k++) {
                    res.a[i][j] =
                        (res.a[i][j] + a[i][k] * b.a[k][j] % mod) % mod;
                }
            }
        }
        return res;
    }
};

Matrix qpow(Matrix a, int b) {
    Matrix res;
    for (int i = 1; i <= 90; i++) {
        res.a[i][i] = 1;
    }
    while (b) {
        if (b % 2 == 1) {
            res = res * a;
        }
        a = a * a;
        b /= 2;
    }
    return res;
}
void solve() {
    int n, t;
    cin >> n >> t;
    vector<string> p(n + 1);
    for (int i = 1; i <= n; i++) {
        cin >> p[i];
        p[i] = " " + p[i];
    }
    Matrix a, b;
    a.a[1][1] = 1;
    for (int i = 1; i <= n; i++) {
        for (int k = 1; k < 9; k++) {
            b.a[(i - 1) * 9 + k][(i - 1) * 9 + k + 1] = 1;
        }
        for (int j = 1; j <= n; j++) {
            int x = p[i][j] - '0';
            if (x > 0) {
                b.a[(i - 1) * 9 + x][(j - 1) * 9 + 1] = 1;
            }
        }
    }
    a = a * qpow(b, t);
    cout << a.a[1][(n - 1) * 9 + 1];
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    int t = 1;
    // cin >> t;
    while (t--) {
        solve();
    }
    return 0;
}