竞赛算法(8) 线性代数

线性基#

通过有限的线性基, 可以描述无限的线性空间.

算法竞赛中常用的线性基有两种, 实线性空间中的 $\mathbb{R}^n$ 实数线性基和布尔域线性空间 $\mathbb{Z}_2^n$ 中的异或线性基. (布尔域线性空间中加法为异或, 数乘为与).

以异或线性基为例, 我们可以根据一组布尔序列 $X = \{x_1, x_2, ..., x_m\}$ 构造出一组异或线性基 $B=\{b_1, b_2, ..., b_n\}$ . 这组基有如下性质:

任意非空子集的异或和不为0.
对 $X$ 中的任意元素 $x$ , 都可以在 $B$ 中取出若干元素使其异或和为 $x$ .
$B$ 是满足性质1和2的极小向量组.

:::[异或线性基的典型用法]

判断一个数是否能表示成某数集子集的异或和;
求一个数表示成某数集子集异或和的方案;
求某数集的最大/最小/第k大/第k小子集异或和;
求某数在某数集子集异或和中的排名. :::

矩阵快速幂/矩阵加速递推#

矩阵加速递推可以把线性递推式由 $O(n)$ 优化到 $O(k^3\log n)$ . 其中k为递推式的阶数.

如果一个递推式是线性的, 我们可以把它写成矩阵的形式, 一次递推即相当于左乘一次转移矩阵. 我们需要n次递推, 由于矩阵有结合律, 这相当于我们左乘转移矩阵的n次幂. 我们可以先使用快速幂算法求出矩阵的n次幂, 再用初始向量左乘得到结果.

以计算斐波那契数列第n项 $f_n$ 为例.

构造状态向量. $f_i = f_{i - 1} + f_{i - 2}$ . 我们知道他的初始态有两项 $f_0 = 0, f_1 = 1$ . 这暗示我们状态有两项:

F_n = \begin{bmatrix} f_n \\ f_{n-1} \end{bmatrix}

构建转移矩阵:

f_n = 1 \cdot f_{n-1} + 1 \cdot f_{n-2}

f_{n-1} = 1 \cdot f_{n-1} + 0 \cdot f_{n-2}

于是转移矩阵:

M = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}

转移可以写成:

F_i = MF_{i - 1}

做快速幂即可.

1
// 列向量
2
template <class S, int N>
3
struct ColVector {
4
    array<S, N> a;
5
    ColVector(S def = S()) {
6
        fill(a.begin(), a.end(), def);
7
    }
8

9
    S& operator[](int i) {
10
        return a[i];
11
    }
12
    const S& operator[](int i) const {
13
        return a[i];
14
    }
15
};
16

17
// 矩阵
18
template <typename S, int N>
19
struct Matrix {
20
    array<array<S, N>, N> a;
21
    Matrix(S def = S()) {
22
        for (int i = 0; i < N; i++) {
23
            fill(a[i].begin(), a[i].end(), def);
24
        }
25
    }
26

27
    S* operator[](int i) {
28
        return a[i].data();
29
    }
30
    const S* operator[](int i) const {
31
        return a[i].data();
32
    }
33

34
    // 单位矩阵
35
    static Matrix e() {
36
        Matrix ret;
37
        for (int i = 0; i < N; i++) {
38
            ret[i][i] = 1;
39
        }
40
        return ret;
41
    }
42

43
    // 加法
44
    friend Matrix operator+(Matrix lt, const Matrix& rt) {
45
        for (int i = 0; i < N; i++) {
46
            for (int j = 0; j < N; j++) {
47
                lt[i][j] += rt[i][j];
48
            }
49
        }
50
        return lt;
51
    }
52

53
    // 乘法
54
    friend Matrix operator*(const Matrix& lt, const Matrix& rt) {
55
        Matrix ret;
56
        for (int i = 0; i < N; i++) {
57
            for (int k = 0; k < N; k++) {
58
                if (lt[i][k] == S()) {
59
                    continue;
60
                }
61

62
                for (int j = 0; j < N; j++) {
63
                    ret[i][j] += lt[i][k] * rt[k][j];
64
                }
65
            }
66
        }
67
        return ret;
68
    }
69

70
    Matrix pow(ll b) const {
71
        Matrix ret = e();
72
        Matrix t = *this;
73
        while (b > 0) {
74
            if (b & 1) {
75
                ret = ret * t;
76
            }
77

78
            t = t * t;
79
            b >>= 1;
80
        }
81
        return ret;
82
    }
83

84
    ColVector<S, N> operator*(const ColVector<S, N>& vec) {
85
        ColVector<S, N> ret;
86
        for (int i = 0; i < N; i++) {
87
            for (int k = 0; k < N; k++) {
88
                if (a[i][k] == S()) {
89
                    continue;
90
                }
91

92
                ret[i] = ret[i] + (a[i][k] * vec[k]);
93
            }
94
        }
95
        return ret;
96
    }
97
};

矩阵快速幂还有其他trick, 这里简单介绍一些.

邻接矩阵的快速幂#

以下邻接矩阵的定义: adj[v][u]表示从u到v的边, 请注意这和常见的邻接矩阵是反过来的(方便左乘).

首先考虑这样一个问题: 对于给定的简单图, 顶点数n <= 100, 计算由k <= 1e9步组成的路径数.

我们可以很容易想到这样一个递推方法, 设dp[i][j]表示以节点i为结尾的, 由j步组成的路径数, 设c[j][i]为i到j是否有边的指示器(有边为1, 无边为0), 则:

dp[1][j] = c[1][1] * dp[1][j - 1] + c[1][2] * dp[2][j - 1] + ... + c[1][n] * dp[n][j - 1].
dp[2][j] = c[2][1] * dp[1][j - 1] + ...
…
dp[n][j] = c[n][1] * dp[1][n - 1] + ...

我们可以发现这可以写成矩阵形式. 状态向量即为DP = dp[1][j], dp[2][j], ..., dp[n][j], 转移矩阵为c[i][j]组成的矩阵, 我们会发现这个矩阵正是图的邻接矩阵. 所以使用邻接矩阵建图, 再对邻接矩阵做快速幂即可.

(min, +)矩阵乘法(新定义)#

回想floyd算法的过程. dp[i][j] = min(dp[i][t] + dp[t][j]), 我们会发现dp[i][t]的第二维和dp[t][j]的第一维是相同的, 而且每次t需要遍历所有点, 这恰好和矩阵乘法有相似之处(矩阵乘法: dp[i][j] = dp[i][1] * dp[1][j] + dp[i][2] * dp[2][j] + ...), 也就是说原来的乘法相当于新的加法, 原来的加法变成了新的min操作, 可以证明这样的矩阵乘法是满足结合律的, 所以可以使用快速幂优化. 这样的矩阵乘法也可以叫最短路径矩阵乘法, 因为常用来求拥有k条边的最短路径问题.

例题: CF102644F

代码:

1
template <int N>
2
struct Matrix_mp {
3
    array<array<ll, N>, N> a;
4
    Matrix_mp() {
5
        for (int i = 0; i < N; i++) {
6
            fill(a[i].begin(), a[i].end(), INFLL);
7
        }
8
    }
9

10
    ll* operator[](int i) {
11
        return a[i].data();
12
    }
13
    const ll* operator[](int i) const {
14
        return a[i].data();
15
    }
16

17
    static Matrix_mp e() {
18
        Matrix_mp ret;
19
        for (int i = 0; i < N; i++) {
20
            ret[i][i] = 0;
21
        }
22
        return ret;
23
    }
24

25
    friend Matrix_mp operator+(Matrix_mp lt, const Matrix_mp& rt) {
26
        for (int i = 0; i < N; i++) {
27
            for (int j = 0; j < N; j++) {
28
                lt[i][j] = min(lt[i][j], rt[i][j]);
29
            }
30
        }
31
        return lt;
32
    }
33
    friend Matrix_mp operator*(const Matrix_mp& lt, const Matrix_mp& rt) {
34
        Matrix_mp ret;
35
        for (int i = 0; i < N; i++) {
36
            for (int k = 0; k < N; k++) {
37
                if (lt[i][k] >= INFLL) {
38
                    continue;
39
                }
40

41
                for (int j = 0; j < N; j++) {
42
                    ll now;
43
                    if (rt[k][j] >= INFLL) {
44
                        now = INFLL;
45
                    } else {
46
                        now = lt[i][k] + rt[k][j];
47
                    }
48
                    ret[i][j] = min(ret[i][j], now);
49
                }
50
            }
51
        }
52
        return ret;
53
    }
54

55
    Matrix_mp pow(ll b) const {
56
        Matrix_mp ret = e();
57
        Matrix_mp t = *this;
58
        while (b > 0) {
59
            if (b & 1) {
60
                ret = ret * t;
61
            }
62

63
            t = t * t;
64
            b >>= 1;
65
        }
66
        return ret;
67
    }
68

69
    ColVector<ll, N> operator*(const ColVector<ll, N>& vec) {
70
        ColVector<ll, N> ret(INFLL);
71
        for (int i = 0; i < N; i++) {
72
            for (int k = 0; k < N; k++) {
73
                if (a[i][k] >= INFLL) {
74
                    continue;
75
                }
76
                ll now;
77
                if (a[i][k] >= INFLL || vec[k] >= INFLL) {
78
                    now = INFLL;
79
                } else {
80
                    now = a[i][k] + vec[k];
81
                }
82
                ret[i] = min(ret[i], now);
83
            }
84
        }
85
        return ret;
86
    }
87
};