【线性代数】矩阵乘法与线性DP优化

【线性代数】矩阵乘法与线性DP优化


前言#

现在有一道题目如下:

输入一个整数 n (n1018n \leq 10^{18}), 求第 n 个斐波那契数。

众所周知斐波那契数列的递推公式是:fi=fi1+fi2f_i = f_{i-1} + f_{i-2}。通过O(N)O(N)的时间递推可以在 1s 内求出前 1e8 的斐波那契数,但是题目的范围是 1e18,这要怎么办呢?这时候就要引出我们今天要学习的内容——矩阵与矩阵乘法了。

概念#

什么是矩阵(matrix)#

下面是一个2×22 \times 2的矩阵,其中abcd,是里面的元素,矩阵里的元素可以是数字符号或者数学式。

[abcd]\begin{bmatrix}{a} & {b} \\ {c} & {d}\end{bmatrix}

关于矩阵的其他内容我们不再延申,你现在只要知道矩阵是这么样的一个东西就可以了。

矩阵可以用字母代表,那么矩阵  An×mA_{n \times m}  本质上是一个  n  行  m  列的二维数组。

矩阵乘法#

矩阵之间可以相乘,并且满足结合律与分配律——不满足交换律,在n×m(nm)n \times m(n \ne m)这种不是正方形的矩阵中,交换前后两个矩阵相乘会导致结果矩阵的形状不同,我们会在后面解释。

矩阵相乘时,相乘矩阵的宽高必须有一个相同,否则无法配对。

比如下面这个式子,把矩阵A×B=CA \times B = C,可以这样写:

Cn×s=An×m×B_m×sC*{n\times s}=A*{n\times m}\:\times B\_{m\times s}

那么它实际上就等于:

Cij=Aik×B_kj(1km)C*{ij}=\sum A*{ik}\times B\_{kj}\left(1\:\le\:k\:\le\:m\right)

其中 k 的遍历过程也就是相乘矩阵的宽高必须有一个相同的原因,否则匹配不了。

写成代码形式,就和弗洛伊德很像,可以把弗洛伊德看成是魔改的矩阵乘法:

CPP
1
2
3
4
for(int i = 1; i <= an; i++)
    for(int k = 1; k <= am; k++)
        for(int j = 1; j <= bm; j++)
            c[i][j] = c[i][j] + a[i][k] * b[k][j];

用矩阵优化线性 DP#

回到前言我们说的题目,这里回顾一下:

输入一个整数 n (n1018n \leq 10^{18}), 求第 n 个斐波那契数。

我们每一次转移可以用一个矩阵表示:

[fi1,fi2]×[1110]=[fi,fi1]\left[{f_{i-1}},\:f_{i-2}\right]\times\left[\begin{matrix}{1} & {1} \\ {1} & {0}\end{matrix}\right]=\left[{f_i},\:f_{i-1}\right]

因为是线性的,所以我们很容易发现,每一转移其实就是当前的状态矩阵乘上我们的[1110]\left[\begin{matrix}{1} & {1} \\ {1} & {0}\end{matrix}\right],那么每转移一次,i 就加一,我们先处理出f1,f2f_1,f_2,那么通过n2n-2次转移,我们就可以得到fnf_n的值。

延伸#

所有类似于这样的线性 DP 都可以用矩阵来转移,比如fi=fi1+fi3+fi4f_i = f_{i-1} + f_{i-3} + f_{i-4},这种转移,我们可以构造下面这一个转移矩阵:

[1000010010101001]\begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1\end{bmatrix}

然后使用快速幂可以把时间复杂度从O(N)O(N)优化到O(k3logN)O(k^3 \log N) ,其中 k 是状态是数量,也就是转移矩阵的边长。

下面是实例代码,可以解决fi=fi1+fi3+fi4f_i = f_{i-1} + f_{i-3} + f_{i-4}求解fnf_n并且n1018n \leq 10^{18}的情况。

这里因为乘的顺序不一样(交换了),矩阵不能使用交换律,所以转移矩阵稍有不同。(实际上一般情况下习惯把系数放在第一行,也就是代码中转移矩阵的样子)

CPP
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
#include <iostream>
#include <cstdio>
#include <vector>
using namespace std;

typedef long long ll;
typedef vector<vector<ll>> matrix;

matrix multiply(matrix a, matrix b, ll mod)
{
    ll n = a.size();
    matrix c(n, vector<ll>(n));
    for (ll i = 0; i < n; i++)
        for (ll j = 0; j < n; j++)
            for (ll k = 0; k < n; k++)
                c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % mod;
    return c;
}

matrix matrix_pow(matrix a, ll n, ll mod)
{
    ll m = a.size();
    matrix res(m, vector<ll>(m));
    for (ll i = 0; i < m; i++)
        res[i][i] = 1;
    while (n)
    {
        if (n & 1)
            res = multiply(res, a, mod);
        a = multiply(a, a, mod);
        n >>= 1;
    }
    return res;
}

int main()
{
    matrix a = {{1, 0, 1, 1},
                {1, 0, 0, 0},
                {0, 1, 0, 0},
                {0, 0, 1, 0}}; // 转移矩阵
    ll n;
    scanf("%lld", &n); // 2, 4, 6, 9
    matrix ans = matrix_pow(a, n - 4, 10007);
    printf("%lld", (ll)(2 * ans[0][3] + 4 * ans[0][2] + 6 * ans[0][1] + 9 * ans[0][0]) % 10007);
    return 0;
}