矩阵乘法初涉

数学

本文简单介绍了矩阵乘法的基本概念和一些有趣的应用.

作者水平有限,错误在所难免,欢迎大家批评指正!

矩阵快速幂是什么东西?

矩阵快速幂用于优化 常线性齐次递推 , 可以把一个线性时间复杂度的算法优化到 $\log$ 级别.需要注意,要优化的递推式必须是 常线性齐次递推.

常线性齐次递推(线性常系数齐次递推关系)

线性常系数齐次递推关系的定义(百度百科)
$a_n+c_1a_{n-1}+c_2a_{n-2}+...+c_ka_{n-k}=0$
$a_0=d_0,a_1=d_1,...,a_{k-1}=d_{k-1}$
若$c_1,c_2,...,c_k,d_0,d_1,...,d_{k-1}$都是常数，则称之为 $k$ 阶的线性常系数齐次递推关系.

从定义可以看出,线性常系数齐次递推关系就是突出了 常系数 和 线性 两个特点.

我们判断一个递推式是否可以用矩阵快速幂主要就是根据这两个特点来判断的,这个问题在实际的应用中需要注意,常常容易忽略.尤其是常系数这个特点.

矩阵快速幂的一般做法

我们以广为人知的斐波那契数列作为例子, 众所周知, 斐波那契数列的递推式是 $F_n = F_{n-1} + F_{n-2} (n \geq 2)$, 同时我们令 $F_0 = 0, F_1 = 1$. 可以看出,斐波那契数列的递推式是常线性齐次递推. 那么我们就可以考虑用矩阵快速幂优化它.

怎么优化呢?一般对于一个 $k$ 阶的常线性齐次递推, 我们需要构造一个 $k \cdot k$ 的矩阵. 这个可能不是很好理解, 换句话说, 为了从第 $n-1$ 项推到第 $n$ 项, 我们需要储存多少信息, 就需要多少大的方阵. 这有点像动态规划状态的设计.

那么我们考虑斐波那契数列的递推式, 如果要求 $F_n$ , 我们需要知道什么呢?

显然, 我们需要知道 $F_{n-1}$ 和 $F_{n-2}$, 也就是说, 我们需要存储 $2$ 个信息, 所以我们需要构造一个 $2 \times 2$ 的矩阵. 具体应该如何构造呢? 我一般是这样做的:

$\left[\begin{array}{cc}? & ? \\ ? & ? \end{array} \right]^n$$\left[\begin{array}{cc}F_{n-2} \\ F_{n-1} \end{array} \right]^n$$\left[\begin{array}{cc}F_{n-1} \\ F_{n} \end{array} \right]^n$

然后根据矩阵乘法的相关性质, 填入值使等式成立:

$\left[\begin{array}{cc}0 & 1 \\ 1 & 1 \end{array} \right]^n$$\left[\begin{array}{cc}F_{n-2} \\ F_{n-1} \end{array} \right]^n$$\left[\begin{array}{cc}F_{n-1} \\ F_{n} \end{array} \right]^n$

一般这样就可以了.

另外一种形式

其实斐波那契数列还可以用一个方阵的幂表示, 就像这个方阵:$\left[\begin{array}{cc}0 & 1 \\ 1 & 1 \end{array} \right]^n$ 的左下角就是 $F_n$ 的值, 有些矩阵乘法可以用方阵的幂表示, 原因是其初值的特殊性导致初值方阵可以被省略.

一些题目

A. 后四位

求斐波那契数列第 $n$ 项 $F_n$ 的后四位(不足四位输出所有位).

这是一道模板题, 其实就是求 $F_n \ mod \ 10000$;

直接矩阵快速幂转移即可. 这道题的代码作为模板加了注释.

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
int n;
struct Matrix{ //矩阵类型
    int mat[3][3];
    Matrix(){ //默认构造函数(用于初始化矩阵)
        memset(mat, 0, sizeof(mat));
    }
    Matrix operator *(const Matrix &b){ //重载乘法
        Matrix Res;
        for (int i = 1; i <= 2; ++i)
            for (int j = 1; j <= 2; ++j)
                for (int k = 1; k <= 2; ++k)
                    Res.mat[i][j] = (Res.mat[i][j] + mat[i][k] * b.mat[k][j]) % 10000;
        return Res;
    }
};
Matrix KEY(){ //函数，返回系数方阵
    Matrix Res;
    Res.mat[1][2] = 1;
    Res.mat[2][1] = 1;
    Res.mat[2][2] = 1;
    return Res;
}
Matrix E(){ //函数，返回单位矩阵E
    Matrix Res;
    Res.mat[1][1] = Res.mat[2][2] = 1;
    return Res;
}
Matrix Qpow(Matrix x, int k){ //矩阵快速幂
    Matrix Res = E(), T = x;
    while (k){
        if (k & 1) Res = Res * T; //二进制拆分
        T = T * T;
        k >>= 1;
    }
    return Res;
}
int main(){
    while (scanf("%d", &n) == 1){
        if (n == -1) break;
        Matrix x = KEY();
        if (n == 0) printf("0\n");
        else if (n == 1) printf("1\n");
        else{
            x = Qpow(x, n-1);
            printf("%d\n", x.mat[2][2]);
        }
    }
    return 0;
}

B. Matrix Power Series

给出一个 $n \times n$ 的矩阵 $A$, 求 $A + A^2 + \dots + A^k$

考虑分治解决这个问题, 用类似快速幂的方法, 讨论 $k$ 的奇偶性.

如果 $k$ 是偶数, 那么 $A + A^2 + \dots + A^k = (A + A^2 + \dots + A^{k \over 2}) + (A + A^2 + \dots + A^{k \over 2}) A^{k \over 2}$, 递归解决此问题的同时用快速幂计算 $A^k$, 时间复杂度可以做到 $O(n^3 \log k)$. 具体的实现方法可以看代码. 当然如果要优化常数可以写成非递归的版本.

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
int n, k, m, mat[35][35];
inline int mo(int x){
    return x >= m ? x - m : x;
}
struct Matrix{
    int mat[35][35];
    Matrix(){
        memset(mat, 0, sizeof(mat));
    }
    Matrix(int a[35][35], int n){ //用二维数组初始化矩阵
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                mat[i][j] = a[i][j];
    }
    Matrix operator +(const Matrix &b){ //重载加法
        Matrix Res;
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                Res.mat[i][j] = mo(mat[i][j] + b.mat[i][j]);
        return Res;
    }
    Matrix operator *(const Matrix &b){
        Matrix Res;
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                for (int k = 1; k <= n; ++k)
                    Res.mat[i][j] = (Res.mat[i][j] + mat[i][k] * b.mat[k][j]) % m;
        return Res;
    }
    void out(){ //调试输出
        for (int i = 1; i <= n; ++i){
            for (int j = 1; j < n; ++j)
                printf("%d ", mat[i][j]);
            printf("%d\n", mat[i][n]);
        }
    }
}x;
Matrix E(){
    Matrix Res;
    for (int i = 1; i <= n; ++i)
        Res.mat[i][i] = 1;
    return Res;
}
Matrix ZERO(){ //返回 0 矩阵
    Matrix Res;
    return Res;
}
typedef pair<Matrix,Matrix> Pair;
Pair solve(int k){
    if (k == 0) return make_pair(ZERO(), E());
    Pair T = solve(k/2);
    if (k % 2 == 0) return make_pair(T.first + T.first * T.second, T.second * T.second);
    else return make_pair(T.first + T.first * T.second + T.second * T.second * x, T.second * T.second * x);
}
int main(){
    scanf("%d%d%d", &n, &k, &m);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            scanf("%d", &mat[i][j]);
    Matrix T(mat, n);
    x = T;
    solve(k).first.out();
    return 0;
}

C. 计算

给出 $n$, $F_i$ 表示斐波那契数列的第 $n$ 项, 计算 $\Sigma_{i=1}^n F_i \cdot C_n^i$

我们可以发现这个式子很像二项式定理

同时 $\left[\begin{array}{cc}0 & 1 \\ 1 & 1 \end{array} \right]^n$ 的左下角就是 $F_n$ 的值, 所以代入可以得到 $\Sigma_{i=1}^n \left[\begin{array}{cc}0 & 1 \\ 1 & 1 \end{array} \right]^n \cdot C_n^i$, 即 $(\left[\begin{array}{cc}0 & 1 \\ 1 & 1 \end{array} \right] + E)^n$

由于我是打表的, 就不给出代码了 233.

D. [ZJOI2005]沼泽鳄鱼

题目链接: https://www.luogu.org/problemnew/show/P2579

这是一个类似于 $DP$ 优化的题目, 可以发现转移的周期是 $12$, 所以只需要预处理 $12$ 种转移的矩阵即可.

主要思路是对于完整的周期（也就是 $12$ 种转移完整出现的轮回）用快速幂加速，剩余的转移直接一一计算。

关于矩阵乘法 @fengsongquan 有一个常数优化，矩阵乘法是快速幂中主要的时间消耗，所以对于矩阵乘法的优化是切实有意义的。

matrix operator *(const matrix &b){
        matrix res;
        for (int i = 1; i <= n; ++i)
            for (int k = 1; k <= n; ++k)
                if (mat[i][k] != 0){ //优化，%%% @fengsongquan
                    for (int j = 1; j <= n; ++j)
                        res.mat[i][j] = (res.mat[i][j] + mat[i][k] * b.mat[k][j]) % 10000;
                }
        return res;
}

这是优化后的代码，可以看到，其主要优化了对答案没有贡献的部分。事实证明，对于 $n$ 比较大的情况，这种优化是非常有用的，因为矩阵中会有比较多数量的 $0$，这个时候这个优化就比较有用了。

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <cstring>
using namespace std;
int n, m, nfish, st, ed, k;
int mat[55][55];
struct Fish{
    int k, loc[4];
} fi[30];
struct matrix{
    int mat[55][55];
    matrix(){
        memset(mat, 0, sizeof(mat));
    }
    matrix(const int a[55][55]){
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                mat[i][j] = a[i][j];
    }
    matrix operator *(const matrix &b){
        matrix res;
        for (int i = 1; i <= n; ++i)
            for (int k = 1; k <= n; ++k)
                if (mat[i][k] != 0){ //优化，%%% @fengsongquan
                    for (int j = 1; j <= n; ++j)
                        res.mat[i][j] = (res.mat[i][j] + mat[i][k] * b.mat[k][j]) % 10000;
                }
        return res;
    }
} situ[13];
matrix E(){
    matrix res;
    for (int i = 1; i <= n; ++i)
        res.mat[i][i] = 1;
    return res;
}
matrix Qpow(const matrix &x, int k){
    matrix res = E(), t = x;
    while (k){
        if (k & 1) res = res * t;
        t = t * t;
        k >>= 1;
    }
    return res;
}
void clear(matrix &x, int id){
    for (int i = 1; i <= n; ++i)
        x.mat[i][id] = 0;
}
int main(){
    scanf("%d%d%d%d%d", &n, &m, &st, &ed, &k);
    ++st; ++ed; //坐标转换
    memset(mat, 0, sizeof(mat));
    for (int i = 1; i <= m; ++i){
        int x, y; scanf("%d%d", &x, &y);
        ++x; ++y;
        mat[x][y] = mat[y][x] = 1;
    }
    scanf("%d", &nfish);
    for (int i = 1; i <= nfish; ++i){
        scanf("%d", &fi[i].k);
        for (int j = 0; j < fi[i].k; ++j)
            scanf("%d", &fi[i].loc[j]), ++fi[i].loc[j];
    }
    situ[0] = E();
    for (int i = 1; i <= 12; ++i){
        for (int j = 1; j <= n; ++j)
            for (int k = 1; k <= n; ++k)
                situ[i].mat[j][k] = mat[j][k];
        for (int j = 1; j <= nfish; ++j)
            clear(situ[i], fi[j].loc[i%fi[j].k]); //构造转移矩阵
        situ[0] = situ[0] * situ[i];
    }
    matrix ans = Qpow(situ[0], k/12); //批量转移
    for (int i = 1; i <= k%12; ++i)
        ans = ans * situ[i]; //处理剩余的部分
    printf("%d\n", ans.mat[st][ed]);
    return 0;
}