跳转至

矩阵

本文转载（或修改）自 OI-Wiki

本文介绍线性代数中一个非常重要的内容——矩阵（Matrix），主要讲解矩阵的性质、运算，以及矩阵乘法的一些应用。

线性变换

视频讲解（来源：3Blue1Brown）

所谓"变换"，可以认为是一种运动，也可以认为是一个输入输出均为向量的函数。

我们从运动的观念来看，我们知道平面直角坐标系 \(xOy\) 的基底是\(\boldsymbol{i} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}\)和\(\boldsymbol{j} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}\)

我们将这个坐标平面旋转、缩放、翻转，这个过程可以认为是"变换"

而如果一个变换满足： - 原点依旧是原点，位置不变 - 直线变换后依旧是直线，且相互平行的直线变换后依旧平行 那么这种变换被称为线性变换

根据线性变换的特点，我们发现，如果向量 \(\boldsymbol{a} = x\boldsymbol{i} + y\boldsymbol{j}\)，那么线性变换之后仍然有 \(\boldsymbol{a} = x\boldsymbol{i} + y\boldsymbol{j}\)

因此我们可以通过对基向量 \(\boldsymbol{i}\) 和 \(\boldsymbol{j}\) 的变换来描述线性变换

设基向量\(\boldsymbol{i} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}\)在变换后为\(\begin{bmatrix} a \\ c \end{bmatrix}\)，\(\boldsymbol{j} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}\)在变换后为\(\begin{bmatrix} b \\ d \end{bmatrix}\)，我们将这两个变换后的向量记到一起，即

\[\boldsymbol{A} = \begin{bmatrix} a & b \\ c & d \end{bmatrix}\]

这个就是矩阵

一个向量\(\boldsymbol{p} = \begin{bmatrix} x \\ y \end{bmatrix}\)在矩阵上的变换认为是向量与矩阵的乘法，记得到的向量为 \(\boldsymbol{q}\)，那么有

\[\boldsymbol{q} = \boldsymbol{A}\boldsymbol{p} = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = x \begin{bmatrix} a \\ c \end{bmatrix} + y \begin{bmatrix} b \\ d \end{bmatrix} = \begin{bmatrix} ax + by \\ cx + dy \end{bmatrix}\]

根据这个定义，线性变换也可以推广到更高维

定义

对于矩阵 \(A\)，主对角线是指 \(A_{i,i}\) 的元素。

一般用 \(I\) 来表示单位矩阵，就是主对角线上为 \(1\)，其余位置为 \(0\)。

同型矩阵

两个矩阵，行数与列数对应相同，称为同型矩阵。

方阵

行数等于列数的矩阵称为方阵。方阵是一种特殊的矩阵。对于「 \(n\) 阶矩阵」的习惯表述，实际上讲的是 \(n\) 阶方阵。阶数相同的方阵为同型矩阵。

研究方程组、向量组、矩阵的秩的时候，使用一般的矩阵。研究特征值和特征向量、二次型的时候，使用方阵。

主对角线

方阵中行数等于列数的元素构成主对角线。

对称矩阵

如果方阵的元素关于主对角线对称，即对于任意的 \(i\) 和 \(j\)， \(i\) 行 \(j\) 列的元素与 \(j\) 行 \(i\) 列的元素相等，则将方阵称为对称矩阵。

对角矩阵

主对角线之外的元素均为 \(0\) 的方阵称为对角矩阵，一般记作：

\[\operatorname{diag}\{\lambda_1,\cdots,\lambda_n\}\]

式中的 \(\lambda_1,\cdots,\lambda_n\) 是主对角线上的元素。

对角矩阵是对称矩阵。

如果对角矩阵的元素均为 \(1\)，称为单位矩阵，记为 \(I\)。只要乘法可以进行，无论形状，任何矩阵乘单位矩阵仍然保持不变。

三角矩阵

如果方阵主对角线左下方的元素均为 \(0\)，称为上三角矩阵。如果方阵主对角线右上方的元素均为 \(0\)，称为下三角矩阵。

两个上（下）三角矩阵的乘积仍然是上（下）三角矩阵。如果对角线元素均非 \(0\)，则上（下）三角矩阵可逆，逆也是上（下）三角矩阵。

单位三角矩阵

如果上三角矩阵 \(A\) 的对角线全为 \(1\)，则称 \(A\) 是单位上三角矩阵。如果下三角矩阵 \(A\) 的对角线全为 \(1\)，则称 \(A\) 是单位下三角矩阵。

两个单位上（下）三角矩阵的乘积仍然是单位上（下）三角矩阵，单位上（下）三角矩阵的逆也是单位上（下）三角矩阵。

运算

矩阵的线性运算

矩阵的线性运算分为加减法与数乘，它们均为逐个元素进行。只有同型矩阵之间可以对应相加减。

矩阵的转置

矩阵的转置，就是在矩阵的右上角写上转置「T」记号，表示将矩阵的行与列互换。

对称矩阵转置前后保持不变。

矩阵乘法

矩阵的乘法是向量内积的推广。

矩阵相乘只有在第一个矩阵的列数和第二个矩阵的行数相同时才有意义。

从线性变换的角度来说，矩阵乘法是对两个不同的变换求复合后的总变换。

视频讲解（来源：3Blue1Brown）

矩阵乘法的实现

设 \(A\) 为 \(P \times M\) 的矩阵，\(B\) 为 \(M \times Q\) 的矩阵，设矩阵 \(C\) 为矩阵 \(A\) 与 \(B\) 的乘积，

其中矩阵 \(C\) 中的第 \(i\) 行第 \(j\) 列元素可以表示为：

\[ C_{i,j} = \sum_{k=1}^MA_{i,k}B_{k,j} \]

在矩阵乘法中，结果 \(C\) 矩阵的第 \(i\) 行第 \(j\) 列的数，就是由矩阵 \(A\) 第 \(i\) 行 \(M\) 个数与矩阵 \(B\) 第 \(j\) 列 \(M\) 个数分别 相乘再相加 得到的。这里的 相乘再相加，就是向量的内积。乘积矩阵中第 \(i\) 行第 \(j\) 列的数恰好是乘数矩阵 \(A\) 第 \(i\) 个行向量与乘数矩阵 \(B\) 第 \(j\) 个列向量的内积，口诀为 左行右列。

矩阵乘法满足结合律，不满足一般的交换律。

值得注意的是，矩阵乘法是从右往左算的。

利用结合律，矩阵乘法可以利用快速幂的思想来优化，这称之为矩阵快速幂。

在OI中，由于线性递推式可以表示成矩阵乘法的形式，也通常用矩阵快速幂来求线性递推数列的某一项。

常数优化

首先对于比较小的矩阵，可以考虑直接手动展开循环以减小常数。

可以重新排列循环以提高空间局部性，这样的优化不会改变矩阵乘法的时间复杂度，但是会得到常数级别的提升。

方阵的逆

方阵 \(A\) 的逆矩阵 \(P\) 是使得 \(A \times P = I\) 的矩阵。

逆矩阵不一定存在。如果存在，可以使用高斯消元进行求解。

参考代码

一般来说，可以用一个二维数组来模拟矩阵。

实现 
 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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN = 260;
const ll inf = 1e9;
struct Matrix
{
    ll lx, ly; // lx: 矩阵的行数，ly: 矩阵的列数
    ll val[MAXN][MAXN];
    Matrix() : lx(0), ly(0)
    {
        for (int i = 0; i < MAXN; i++)
            for (int j = 0; j < MAXN; j++)
                val[i][j] = -inf;
    }
    inline void init() // 单位矩阵
    {
        for(int i = 1; i <= min(lx, ly); i++)
            val[i][i] = 1;
    }
    inline ll *operator[](const int i) { return val[i]; }
};
inline Matrix operator*(Matrix a, Matrix b)
{
    Matrix c;
    c.lx = a.lx, c.ly = b.ly;
    for (int i = 1; i <= a.lx; i++)
        for (int k = 1; k <= a.ly; k++)
            for (int j = 1; j <= b.ly; j++)
                c[i][j] += a[i][k] * b[k][j];
    return c;
}

对于向量，可以认为是列数ly为 \(1\) 的特殊矩阵，因此向量以及向量乘矩阵都可以利用上面的矩阵完成

矩阵乘法的应用

矩阵加速递推

P1962 斐波那契数列

大家都知道，斐波那契数列是满足如下性质的一个数列：

\[F_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ F_{n-1}+F_{n-2} \space (n\ge 3) \end{aligned}\right.\]

请你求出 \(F_n \bmod 10^9 + 7\) 的值。

对于 \(100\%\) 的数据，\(1\le n < 2^{63}\)。

我们发现，当 \(n\) 很小时，直接递推即可，但是这里的 \(n\) 特别大，\(O(n)\) 也无法接受

能不能利用矩阵的性质加速？

即已知\(\boldsymbol{p} = \begin{bmatrix} F_{n-1} \\ F_{n-2} \end{bmatrix}\)，需要找到一个矩阵 \(\boldsymbol{A}\)，使得\(\boldsymbol{q} = \boldsymbol{A}\boldsymbol{p} = \begin{bmatrix} F_n \\ F_{n-1} \end{bmatrix}\)

考虑待定系数，设\(A = \begin{bmatrix} a & b \\ c & d \end{bmatrix}\)，那么有

\[\begin{bmatrix} F_n \\ F_{n-1} \end{bmatrix} = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \begin{bmatrix} F_{n-1} \\ F_{n-2} \end{bmatrix} = \begin{bmatrix} aF_{n-1} + bF_{n-2} \\ cF_{n-1} + dF_{n-2} \end{bmatrix}\]

对比系数以及递推公式，我们得到\(\begin{cases} a = 1 \\ b = 1 \\ c = 1 \\ d = 0 \end{cases}\)，即

\[\boldsymbol{A} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}\]

因此利用矩阵乘法的结合律，令\(\boldsymbol{p} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}\)，那么

\[\boldsymbol{q} = \boldsymbol{A}^{n-2}\boldsymbol{p} = \begin{bmatrix} F_n \\ F_{n-1} \end{bmatrix}\]

的第一行第一列就是答案。

用矩阵快速幂即可，时间复杂度 \(O(\log n)\)

注意

矩阵乘法不满足交换律，所以一定不能写成 \(\begin{bmatrix}1 \\ 1\end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-2}\) 的第一行第一列元素。另外，对于 \(n \leq 2\) 的情况，直接输出 \(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
#include <bits/stdc++.h>
using namespace std;
#define endl '\n'
typedef long long ll;
const ll mod = 1e9 + 7;
ll n;
inline void add(ll &a, ll b) { a = (a + b) % mod; }
struct matrix
{
    int lx, ly;
    ll val[3][3];
    inline void init()
    {
        for(int i = 1; i <= min(lx, ly); i++)
            val[i][i] = 1;
    }
    inline ll *operator[](const int i) { return val[i]; }
} trans, vec;
inline matrix operator*(matrix a, matrix b)
{
    matrix c = {};
    c.lx = a.lx;
    c.ly = b.ly;
    for (int i = 1; i <= a.lx; i++)
        for (int k = 1; k <= a.ly; k++)
            for (int j = 1; j <= b.ly; j++)
                add(c[i][j], a[i][k] * b[k][j] % mod);
    return c;
}
inline matrix pow(matrix a, ll b)
{
    matrix res = {a.lx, a.ly, {}};
    for (res.init(); b; b >>= 1ll, a = a * a)
        if (b & 1ll)
            res = res * a;
    return res;
}
int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> n;
    if (n <= 2)
        cout << 1, exit(0);
    trans[1][1] = 1;
    trans[1][2] = 1;
    trans[2][1] = 1;
    trans[2][2] = 0;
    trans.lx = trans.ly = 2;
    vec.lx = 1;
    vec.ly = 2;
    vec[1][1] = 1;
    vec[2][1] = 1;
    vec = pow(trans, n - 2) * vec;
    cout << vec[1][1];
    return 0;
}
P1939 矩阵加速（数列）

已知一个数列 \(F\)，它满足：

\[F_x= \begin{cases} 1 & x \in\{1,2,3\}\\ F_{x-1}+F_{x-3} & x \geq 4 \end{cases} \]

求 \(F\) 数列的第 \(n\) 项对 \(10^9+7\) 取余的值。

一共询问 \(T\) 次。

对于 \(100\%\) 的数据 \(1 \leq T \leq 100\)， \(1 \leq n \leq 2 \times 10^9\)。

和刚刚类似，由于转移涉及 \(F_{n-1}\) 与 \(F_{n-3}\)，故令递推用的向量\(\boldsymbol{p} = \begin{bmatrix} F_n \\ F_{n-1} \\ F_{n-2} \end{bmatrix}\)，转移矩阵\(\boldsymbol{A} = \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & i \end{bmatrix}\)，那么

\[\begin{bmatrix} F_n \\ F_{n-1} \\ F_{n-2} \end{bmatrix} = \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & i \end{bmatrix} \begin{bmatrix} F_{n-1} \\ F_{n-2} \\ F_{n-3} \end{bmatrix} = \begin{bmatrix} aF_{n-1} & bF_{n-2} & cF_{n-3} \\ dF_{n-1} & eF_{n-2} & fF_{n-3} \\ gF_{n-1} & hF_{n-2} & iF_{n-3} \end{bmatrix}\]

对比系数得到转移矩阵

\[\boldsymbol{A} = \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
67
68
69
#include <bits/stdc++.h>
using namespace std;
#define endl '\n'
typedef long long ll;
const ll mod = 1e9 + 7;
inline void add(ll &a, ll b)
{
    a = (a + b) % mod;
}
struct Matrix
{
    int lx, ly;
    ll val[4][4];
    inline void init()
    {
        for (int i = 1; i <= min(lx, ly); i++)
            val[i][i] = 1;
    }
    inline ll *operator[](const int i)
    {
        return val[i];
    }
} trans, vec;
Matrix operator*(Matrix a, Matrix b)
{
    Matrix c = {a.lx, b.ly, {}};
    for (int i = 1; i <= a.lx; i++)
        for (int k = 1; k <= a.ly; k++)
            for (int j = 1; j <= b.ly; j++)
                add(c[i][j], a[i][k] * b[k][j] % mod);
    return c;
}
Matrix qpow(Matrix a, ll b)
{
    Matrix c = {a.lx, a.ly, {}};
    for (c.init(); b; b >>= 1ll, a = a * a)
        if (b & 1ll)
            c = c * a;
    return c;
}
void init()
{
    vec.lx = 3, vec.ly = 1;
    vec[1][1] = vec[2][1] = vec[3][1] = 1;
    trans.lx = trans.ly = 3;
    trans[1][1] = trans[1][3] = trans[2][1] = trans[3][2] = 1;
}
void solve()
{
    ll n;
    cin >> n;
    if (n <= 3)
    {
        cout << 1 << endl;
        return;
    }
    Matrix tvec = qpow(trans, n - 3) * vec;
    cout << tvec[1][1] << endl;
}
int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    init();
    int T;
    cin >> T;
    while (T--)
        solve();
    return 0;
}
P3216 [HNOI2011] 数学作业

小 C 数学成绩优异，于是老师给小 C 留了一道非常难的数学作业题：

给定正整数 \(n,m\)，要求计算 \(\text{Concatenate}(n) \bmod \ m\) 的值，其中 \(\text{Concatenate}(n)\) 是将 \(1 \sim n\) 所有正整数 顺序连接起来得到的数。

例如，\(n = 13\)， \(\text{Concatenate}(n) = 12345678910111213\)。小 C 想了大半天终于意识到这是一道不可能手算出来的题目，于是他只好向你求助，希望你能编写一个程序帮他解决这个问题。

对于 \(30\%\) 的数据，\(1\le n \le 10^6\)；

对于 \(100\%\) 的数据，\(1\le n \le 10^{18}\)， \(1\le m \le 10^9\)。

我们记 \(F_n = \text{Concatenate}(n)\)，那么

\[F_n = F_{n-1} \times 10^k + n\]

这里 \(k = \lg n + 1\) 为 \(n\) 的位数

令答案向量为\(\begin{bmatrix} F_n \\ n + 1 \end{bmatrix}\)，尝试构造转移矩阵，发现不可行

于是再令答案向量为\(\begin{bmatrix} F_n \\ n + 1 \\ 1 \end{bmatrix}\)，那么有转移矩阵\(\begin{bmatrix} 10^k & 1 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{bmatrix}\)

与之前不同的是，当 \(n\) 的位数发生变化的时候要更新 \(k\)，因此分段转移即可

实现 
 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
#include <bits/stdc++.h>
using namespace std;
#define endl '\n'
typedef unsigned long long ll;
ll n, mod;
inline void add(ll &a, ll b) { a = (a + b) % mod; }
struct Matrix
{
    int lx, ly;
    ll val[4][4];
    void init()
    {
        for (int i = 1; i <= min(lx, ly); i++)
            val[i][i] = 1;
    }
    ll *operator[](const int x) { return val[x]; }
} trans, vec;
Matrix operator*(Matrix a, Matrix b)
{
    Matrix c = {a.lx, b.ly, {}};
    for (int i = 1; i <= a.lx; i++)
        for (int k = 1; k <= a.ly; k++)
            for (int j = 1; j <= b.ly; j++)
                add(c[i][j], a[i][k] * b[k][j] % mod);
    return c;
}
Matrix qpow(Matrix a, ll b)
{
    Matrix res = {a.lx, a.ly, {}};
    for (res.init(); b; b >>= 1ull, a = a * a)
        if (b & 1ull)
            res = res * a;
    return res;
}
int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> n >> mod;
    vec.lx = 3, vec.ly = 1;
    vec[3][1] = vec[2][1] = 1;
    trans.lx = trans.ly = 3;
    trans[1][2] = trans[2][2] = trans[2][3] = trans[3][3] = 1;
    for (ll k = 10;; k *= 10)
    {
        trans[1][1] = k % mod;
        if (n < k)
        {
            vec = qpow(trans, n - k / 10 + 1) * vec;
            break;
        }
        vec = qpow(trans, k - k / 10) * vec;
    }
    cout << vec[1][1];
    return 0;
}

小练习

P1349 广义斐波那契数列

广义的斐波那契数列是指形如 \(a_n=p\times a_{n-1}+q\times a_{n-2}\) 的数列。

今给定数列的两系数 \(p\) 和 \(q\)，以及数列的最前两项 \(a_1\) 和 \(a_2\)，另给出两个整数 \(n\) 和 \(m\)，试求数列的第 \(n\) 项 \(a_n\) 对 \(m\) 取模后的结果。

对于 \(100\%\) 的数据，\(p,q,a_1,a_2 \in [0,2^{31}-1]\)， \(1\le n,m \le 2^{31}-1\)。

Tip

转移矩阵是

\[\boldsymbol{A} = \begin{bmatrix} p & q \\ 1 & 0 \end{bmatrix}\]

矩阵加速DP

见矩阵加速DP