矩阵加速(数列)

发表于 | 分类于
本文字数： 4.1k

基本概念

矩阵即由$n$行$m$列的数字(也可以嵌套矩阵)构成的一个方阵。如:$\begin{bmatrix}1&2\\3&4\end{bmatrix}$，$\begin{bmatrix}a\\b\end{bmatrix}$等都是矩阵。

有一种特殊的矩阵叫做单位矩阵，它相当于数字乘法中的$1$。单位矩阵的行列数相同，且每一个满足$i=j$的数字$a_{ij}$的值都为1，其他的数字都为$0$。如$\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}$就是一个$3$行$3$列的单位矩阵。

矩阵和数字、矩阵和矩阵之间可以进行一些运算。

本文将介绍矩阵的几种基本的运算，并利用矩阵来用$O(logn)$的复杂度求出$k$阶常系数线性递推关系数列的第$n$项。

矩阵运算

矩阵与数字的运算

直接用矩阵的每一个数字加减乘除矩阵外的数字即可。

矩阵加减法

设矩阵$A$有$a$行$b$列，矩阵$B$有$c$行$d$列，则仅当满足$a=c$且$b=d$时，两个矩阵才可以相加减。

设运算得到的矩阵为$C$，则$c_{ij} = a_{ij} \pm b_{ij}$。这里的$\pm$对应了矩阵的加减。

矩阵的加法满足交换律和结合律。

矩阵乘法

设矩阵$A$有$a$行$b$列，矩阵$B$有$c$行$d$列，则仅当满足$b=c$时，两个矩阵才可以相乘。

设运算得到的矩阵为$C$，则$c_{ij} = \sum\limits_{k = 1}^{b}a_{ik}\times b_{kj}$。C为$a\times d$大小的矩阵。

矩阵乘法不满足交换律，但满足结合律。

矩阵快速幂

由于矩阵乘法满足结合律，我们就可以和数字的快速幂一样，用$O(logk)$的复杂度求出矩阵的$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
struct matrix {
    int n, m;
    long long data[3][3];
    matrix() { memset(data, 0, sizeof(data)); }
    matrix(int nv, int mv) {
        n = nv, m = mv;
        memset(data, 0, sizeof(data));
    }
    void build() {             //此函数将自身矩阵初始化为单位矩阵
        for (int i = 1; i <= n; i++) data[i][i] = 1;
    }
    matrix operator * (const matrix &a) const {
        matrix x = matrix(n, a.m);
        for (int i = 1; i <= x.n; i++)
            for (int j = 1; j <= x.m; j++)
                for (int k = 1; k <= m; k++)
                    x.data[i][j] = ((data[i][k] * a.data[k][j]) % mod + x.data[i][j]) % mod;
        return x;
    }
};
matrix fastpow(matrix base, long long k) {         //矩阵快速幂
    matrix ans = matrix(base.n, base.m);
    ans.build();
    while (k) {
        if (k & 1) ans = ans * base;
        base = base * base;
        k >>= 1;
    }
    return ans;
}

矩阵加速

线性递推的玄学优化版本，将$O(n)$复杂度的线性递推优化为$O(logn)$。这里之所以说是优化，是因为矩阵加速的实质也是递推，只是形式不同。

引例

我们先来看一个非常著名的数列:斐波那契数列。在这个数列中， $F[1]=F[2]=1$ ， $F[n]=F[n-1]+F[n-2] (n\ge 3)$ 。现在要求出数列的第$n$项， $n$ 在 $long long$ 范围内。对应的洛谷题目

由于$n$比较大，用$O(n)$的递推会超时。我们考虑用矩阵来表示这里的递推关系，即

这里的$A$是另一个矩阵，它只能由我们推导出来。

我们由$F[n]=1\times F[n-1]+1\times F[n-2]$，$F[n-1]=1\times F[n-1]+0\times F[n-2]$得出

我们可以发现，每进行一次$A\times \begin{bmatrix}F[n-1]\\F[n-2]\end{bmatrix}$，我们当前递推到的下标值都会加1，因此我们只需要算出$A^{n-2}\times \begin{bmatrix}F[n-1]\\F[n-2]\end{bmatrix}$对应的矩阵，取这个矩阵中的$F[n]$值，就得到了问题的答案。我们可以先用矩阵快速幂算出$A^{n-2}$，再和$\begin{bmatrix}F[n-1]\\F[n-2]\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
#include<cstdio>
#include<cstring>
using namespace std;

const int MOD = 1e9 + 7;
long long n;
struct matrix {
    int n, m;
    long long data[3][3];
    matrix() { memset(data, 0, sizeof(data)); }
    matrix(int nv, int mv) {
        n = nv, m = mv;
        memset(data, 0, sizeof(data));
    }
    void build() {
        for (int i = 1; i <= n; i++) data[i][i] = 1;
    }
    matrix operator * (const matrix &a) const {
        matrix x = matrix(n, a.m);
        for (int i = 1; i <= x.n; i++)
            for (int j = 1; j <= x.m; j++)
                for (int k = 1; k <= m; k++)
                    x.data[i][j] = ((data[i][k] * a.data[k][j]) % MOD + x.data[i][j]) % MOD;
        return x;
    }
}A;

matrix fastpow(matrix a, long long k) {
    matrix ans = matrix(a.n, a.m), base = a;
    ans.build();
    while (k) {
        if (k & 1) ans = ans * base;
        base = base * base;
        k >>= 1;
    }
    return ans;
}

int main() {
    A.n = A.m = 2;
    A.data[1][1] = A.data[1][2] = A.data[2][1] = 1;
    scanf("%lld", &n);
    if (n <= 2) printf("%d\n", 1);
    else {
        matrix base = matrix(2, 1);
        base.data[1][1] = base.data[2][1] = 1;
        matrix ans = fastpow(A, n - 2) * base;
        printf("%lld\n", ans.data[1][1]);
    }
    return 0;
}