ACM

矩阵相关

矩阵学习记录

定义

矩阵是一个按照长方阵列排列的复数或实数集合。

形如:

$$ A = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix} $$

矩阵的计算

可以类比实数，同样有加减乘除。

矩阵的加法

矩阵的加法一般在两个矩阵为同型矩阵时使用。

$$ A = \begin{bmatrix} 4 & 3 \\ 4 & 0 \\ 1 & 5 \end{bmatrix} ,\quad B = \begin{bmatrix} 1 & 5 \\ 0 & 1 \\ 1 & 3 \end{bmatrix} ,\quad A + B = \begin{bmatrix} 4+1 & 3+5 \\ 4+0 & 0+1 \\ 1+1 & 5+3 \end{bmatrix} = \begin{bmatrix} 5 & 8 \\ 4 & 1 \\ 2 & 8 \end{bmatrix} $$

同时，交换律也适用于矩阵加法。

$$ A + B = B + A $$$$ (A + B) + C = A + (B + C) $$

矩阵的减法

减法与加法类似。

$$ A - B = \begin{bmatrix} 4 & 3 \ 4 & 0 \ 1 & 5 \end{bmatrix}

  • \begin{bmatrix} 1 & 5 \ 0 & 1 \ 1 & 3 \end{bmatrix} = \begin{bmatrix} 4-1 & 3-5 \ 4-0 & 0-1 \ 1-1 & 5-3 \end{bmatrix} = \begin{bmatrix} 3 & -2 \ 4 & -1 \ 0 & 2 \end{bmatrix} $$

矩阵的乘法

两个矩阵能够相乘，必要条件是前一个矩阵的列数等于后一个矩阵的行数。

A 是一个 a * c 的矩阵，B 是一个 c * b 的矩阵，A * B 得到的结果是 a * b 的矩阵。

$$ \begin{aligned} A &= \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix} ,\quad B = \begin{bmatrix} 7 & 8 \\ 9 & 1 \\ 2 & 3 \end{bmatrix} \\ \\ AB &= \begin{bmatrix} (1 \cdot 7 + 2 \cdot 9 + 3 \cdot 2) & (1 \cdot 8 + 2 \cdot 1 + 3 \cdot 3) \\ (4 \cdot 7 + 5 \cdot 9 + 6 \cdot 2) & (4 \cdot 8 + 5 \cdot 1 + 6 \cdot 3) \end{bmatrix} \\ \\ &= \begin{bmatrix} (7 + 18 + 6) & (8 + 2 + 9) \\ (28 + 45 + 12) & (32 + 5 + 18) \end{bmatrix} \\ \\ &= \begin{bmatrix} 31 & 19 \\ 85 & 55 \end{bmatrix} \end{aligned} $$

同时，矩阵乘法满足

$$(AB)C = a(BC)$$

$$(A + B)C = AC + BC$$

$$C(A + B) = CA + CB$$

实数中，有任意数乘以 1 还等于他本身这个概念，在矩阵中，这个 1 即单位矩阵。

单位矩阵左上角到右下角都为 1，其余都是 0。

A 矩阵是一个 n * m 的矩阵，他的单位矩阵大小为 m * 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

struct Matrix {
    vector<vector<ll>> mat;
    Matrix(int m, bool type) {
        mat.assign(m + 1, vector<ll>(m + 1, 0));
        if (type) {
            for (int i = 1; i <= m; i++) {
                mat[i][i] = 1;
            }
        }
    }
};

Matrix operator* (const Matrix &a, const Matrix &b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 0);
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= m; j++) {
            for (int k = 1; k <= m; k++) {
                res.mat[i][j] = (res.mat[i][j] + a.mat[i][k] * b.mat[k][j] % mod) % mod;
            }
        }
    }
    return res;
}

Matrix tpow(Matrix a, ll b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 1);
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

例题

P3390 【模板】矩阵快速幂

照搬上面快速幂代码即可。

 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

#include <bits/stdc++.h>
using namespace std;
#define endl "\n"
typedef long long ll;

const int N = 1e5 + 7, mod = 1e9 + 7;

struct Matrix {
    vector<vector<ll>> mat;
    Matrix(int m, bool type) {
        mat.assign(m + 1, vector<ll>(m + 1, 0));
        if (type) {
            for (int i = 1; i <= m; i++) {
                mat[i][i] = 1;
            }
        }
    }
};

Matrix operator* (const Matrix &a, const Matrix &b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 0);
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= m; j++) {
            for (int k = 1; k <= m; k++) {
                res.mat[i][j] = (res.mat[i][j] + a.mat[i][k] * b.mat[k][j] % mod) % mod;
            }
        }
    }
    return res;
}

Matrix tpow(Matrix a, ll b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 1);
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

void solve() {
    ll n, k;
    cin >> n >> k;
    Matrix a(n, 0);
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            cin >> a.mat[i][j];
        }
    }
    if (k == 0) {
        Matrix res(n, 1);
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                cout << res.mat[i][j] << " ";
            }
            cout << endl;
        }
        return;
    }
    Matrix res = tpow(a, k);
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            cout << res.mat[i][j] << " ";
        }
        cout << endl;
    }
}

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

矩阵加速

我们通常遇到的问题是：某个递推公式（比如斐波那契数列 $f(n) = f(n-1) + f(n-2)$）需要我们求第 $n$ 项，而 $n$ 非常大（例如 $10^{18}$）。线性递推求解的复杂度是 $O(n)$，当 $n$ 巨大时，这显然是不可接受的。矩阵加速的目标，就是将这类问题的复杂度优化到 $O(k^3 \log n)$，其中 $k$ 是递推关系涉及的项数（对于斐波那契数列是2）。

那么什么问题才能使用矩阵加速？

必须是线性递推关系。这意味着每一项都可以表示为前若干项的线性组合（即各项乘以常数再相加）。例如，$f(n) = 2f(n-1) + 3f(n-2) + 5$ 就是一个线性递推。但 $f(n) = f(n-1) \times f(n-2)$ 或 $f(n) = \sqrt{f(n-1)}$ 就不是。

构造转移矩阵是整个算法的核心。

以斐波那契数列为例：$f(n) = f(n-1) + f(n-2)$，$f(1)=1, f(2)=1$。

我们的目标是找到一个矩阵 M，使得下面关系成立

$$\begin{pmatrix} f(n) \\ f(n-1) \end{pmatrix} = M \times \begin{pmatrix} f(n-1) \\ f(n-2) \end{pmatrix}$$

设 $M = \begin{pmatrix} a & b \ c & d \end{pmatrix}$。

根据矩阵乘法法则：

$$\begin{pmatrix} a \cdot f(n-1) + b \cdot f(n-2) \\ c \cdot f(n-1) + d \cdot f(n-2) \end{pmatrix} = \begin{pmatrix} f(n) \\ f(n-1) \end{pmatrix}$$

对比两边的矩阵，我们得到两个等式：

$a \cdot f(n-1) + b \cdot f(n-2) = f(n)$ $c \cdot f(n-1) + d \cdot f(n-2) = f(n-1)$

根据斐波那契定义 $f(n) = 1 \cdot f(n-1) + 1 \cdot f(n-2)$，我们可以轻松确定第一行：$a=1, b=1$。 对于第二行，为了得到 $f(n-1)$，我们只需要 $1 \cdot f(n-1) + 0 \cdot f(n-2)$，所以 $c=1, d=0$。 因此，我们的转移矩阵 $M$ 就是：

$$M = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$$

此时得到公式：$\begin{pmatrix} f(n) \ f(n-1) \end{pmatrix} = M^{n-1} \times \begin{pmatrix} f(1) \ f(0) \end{pmatrix} = \begin{pmatrix} 1 & 1 \ 1 & 0 \end{pmatrix}^{n-1} \times \begin{pmatrix} 1 \ 0 \end{pmatrix}$

求解第 n 个斐波那契数，转化为了求矩阵 $M$ 的 $n-1$ 次幂。

 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

#include <bits/stdc++.h>
using namespace std;
#define endl "\n"
typedef long long ll;

const int N = 1e5 + 7, mod = 1e9 + 7;

struct Matrix {
    vector<vector<ll>> mat;
    Matrix(int m, bool type) {
        mat.assign(m + 1, vector<ll>(m + 1, 0));
        if (type) {
            for (int i = 1; i <= m; i++) {
                mat[i][i] = 1;
            }
        }
    }
};

Matrix operator* (const Matrix &a, const Matrix &b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 0);
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= m; j++) {
            for (int k = 1; k <= m; k++) {
                res.mat[i][j] = (res.mat[i][j] + a.mat[i][k] * b.mat[k][j] % mod) % mod;
            }
        }
    }
    return res;
}

Matrix tpow(Matrix a, ll b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 1);
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

void solve() {
    ll n;
    cin >> n;
    Matrix M(2, 0);
    M.mat[1][1] = 1;
    M.mat[1][2] = 1;
    M.mat[2][1] = 1;
    M.mat[2][2] = 0;
    Matrix res = tpow(M, n - 1);
    cout << res.mat[1][1] << endl;
}

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

例题 1

P1939 矩阵加速（数列）

$$ a_x = \begin{cases} 1 & x \in \{1, 2, 3\} \\ a_{x-1} + a_{x-3} & x \geq 4 \end{cases} $$

因为 $a_n$ 需要用到 $a_{n-3}$，可以基本确定 M 矩阵大小应该为 $3 * 3$

根据原公式得到

$$ \begin{align*} a[n] &= a[n-1] \times 1 + a[n-2] \times 0 + a[n-3] \times 1 \\ a[n-1] &= a[n-1] \times 1 + a[n-2] \times 0 + a[n-3] \times 0 \\ a[n-2] &= a[n-1] \times 0 + a[n-2] \times 1 + a[n-3] \times 0 \end{align*} $$

所以 M 应该为

$$ \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix} $$
 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

#include <bits/stdc++.h>
using namespace std;
#define endl "\n"
typedef long long ll;

const int N = 1e5 + 7, mod = 1e9 + 7;

struct Matrix {
    vector<vector<ll>> mat;
    Matrix(int m, bool type) {
        mat.assign(m + 1, vector<ll>(m + 1, 0));
        if (type) {
            for (int i = 1; i <= m; i++) {
                mat[i][i] = 1;
            }
        }
    }
};

Matrix operator* (const Matrix &a, const Matrix &b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 0);
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= m; j++) {
            for (int k = 1; k <= m; k++) {
                res.mat[i][j] = (res.mat[i][j] + a.mat[i][k] * b.mat[k][j] % mod) % mod;
            }
        }
    }
    return res;
}

Matrix tpow(Matrix a, ll b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 1);
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

Matrix M(3, 0);

void init() {
    M.mat[1][1] = M.mat[2][1] = M.mat[1][3] = M.mat[3][2] = 1;
}

void solve() {
    ll n;
    cin >> n;
    Matrix res = tpow(M, n - 1);
    cout << res.mat[1][1] << endl;
}

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

例题 2

MC0365 魔法链路2

小码妹固定了自己的魔法链路，使其稳定为一个长度为 $n$ 的序列，小码妹可以把自己学会的 $m$ 种魔法放置在魔法链路序列中，一次施法即可将魔法全部输出，获得巨大的威力。也就是说有 $n$ 个位置，每个位置可以放 $m$ 种数字。

一次施法的威力由魔法本身的威力与共鸣程度决定，每种魔法都具有自己的威力值 $a_i$，而魔法共鸣只存在于在魔法链路中相邻的两个魔法之间，可以将其鸣值看成一个矩阵，其中 $g_{ij}$ 表示前一次施放魔法 $i$，这一次施放魔法 $j$ 产生的魔法共鸣值。链路中所有魔法的威力值之和与所有相邻魔法共鸣之和的加和为这次施法的威力。

强大的魔法少女小码妹自然想让自己的施法威力最大，你能帮她算算这个最大值吗？

格式

输入格式: 第一行输入两个整数 $n, m$ ($1 \le n \le 10^6, 1 \le m \le 100$)，含义如题面所示； 第二行输入 $m$ 个整数 $a_i$ ($1 \le a_i \le 10^5$)，含义如题面所示； 接下来输入一个 $m \times m$ 的矩阵，表示魔法共鸣值。矩阵的第 $i$ 行第 $j$ 个数即为 $g_{ij}$ ($1 \le g_{ij} \le 10^5$)。

输出格式: 输出一行一个整数，表示你的答案。

用矩阵来思考，$f_{i,j}$ 就是从 $i$ 转移到 $j$ 的最大代价，每次转移造成的贡献其实是固定的，我们就可以通过一个 M 来计算这 n 次转移后的结果。

 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

#include <bits/stdc++.h>
using namespace std;
#define endl "\n"
typedef long long ll;

const int N = 1e5 + 7, mod = 1e9 + 7, INF = 1e9;

struct Matrix {
    vector<vector<ll>> mat;
    Matrix(int m, int type) {
        mat.assign(m + 1, vector<ll>(m + 1, 0));
        if (type == 1) {
            for (int i = 1; i <= m; i++) {
                mat[i][i] = 1;
            }
        } else if (type == 2) {
            mat.assign(m + 1, vector<ll>(m + 1, -INF));
            for (int i = 1; i <= m; i++) {
                mat[i][i] = 0;
            }
        } else {
            mat.assign(m + 1, vector<ll>(m + 1, -INF));
        }
    }
};

Matrix operator* (const Matrix &a, const Matrix &b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 3);
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= m; j++) {
            for (int k = 1; k <= m; k++) {
                res.mat[i][j] = (max(res.mat[i][j], a.mat[i][k] + b.mat[k][j]));
            }
        }
    }
    return res;
}

Matrix tpow(Matrix a, ll b) {
    int m = a.mat.size() - 1;
    Matrix res(m, 2);
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}


void solve() {
    int n, m;
    cin >> n >> m;
    vector<ll> a(m + 1);
    vector<vector<ll>> g(m + 1, vector<ll>(m + 1, 0));
    Matrix M(m, 0);
    for (int i = 1; i <= m; i++) {
        cin >> a[i];
    }
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= m; j++) {
            cin >> g[i][j];
            M.mat[i][j] = g[i][j] + a[j];
        }
    }
    ll ans = 0;
    if (n == 1) {
        for (int i = 1; i <= m; i++) {
            ans = max(ans, a[i]);
        }
        cout << ans << endl;
        return;
    }
    Matrix res = tpow(M, n - 1);
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= m; j++) {
            ans = max(ans, res.mat[j][i] + a[j]);
        }
    }
    cout << ans << endl;
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    int t = 1;
    // cin >> t;
    while (t--) solve();
    return 0;
}
Matrix
Licensed under CC BY-NC-SA 4.0
最后更新于 Sep 11, 2025 16:27 +0800
© 2024 - 2025 KiritoXD
发表了95篇文章 · 总计15万2千字
永远相信美好的事情即将发生。
使用 Hugo 构建
主题 Stack 由 Jimmy 设计