矩阵乘法学习笔记

【前言】

矩阵乘法通常用于加速线性递推式，能够在较为优秀的时间里完成状态的转移。

题目特征通常有：

1. 单次转移规则简单但转移次数很多。
2. 转移规则满足结合律。

当然还可以扩展到图上应用，之后都会提到。

【前置芝士】

1.【矩阵】

一个 \(n\times m\) 的矩阵可以看成是一个 \(n\times m\) 的二维数组。

它同样有着和实数一样的某些运算。

2.【矩阵加法】

只有规模相同的矩阵才能相加减，规则为对应位置相加减，即若 \(C=A+B\)，则有：

\[C_{i,j}=A_{i,j}+B_{i,j} \]

减法同理。

3.【矩阵乘法】

较为复杂，但同样可以理解。

假设 \(A\) 是 \(n\times m\) 的矩阵，\(B\) 是 \(m\times p\) 的矩阵，\(C=A\times B\)，则有：

\[C_{i,j}=\sum_{k=1}^m A_{i,k}\times B_{k,j} \]

显然，只有第一个矩阵的列数等于第二个矩阵的行数时，两个矩阵才可以相乘。

值得思考的是为什么矩阵乘法的定义这么复杂，换句话说，为什么不指定更简单的运算规则达到同样的效果呢？

其实，从矩阵的本质出发，我们很容易找到答案。

矩阵的本质是简化的线性方程组，每一个位置代表了对应未知数的系数。

可以尝试随便写几个线性方程组，简化为矩阵后进行相乘，发现结果确实是上述规则。（具体可以看这里）

同时根据这一点，我们也可以证明矩阵乘法是否符合交换律、结合律、分配率。

结果是：矩阵乘法满足结合律和分配率，不满足交换律。

4.【单位矩阵】

类似于普通乘法中的 \(1\)，设单位矩阵为 \(S\)，有 \(A\times S=A\)。

其构造方式为：

\[S=\begin{vmatrix}1&0&0&0&...&0\\0&1&0&0&...&0\\0&0&1&0&...&0\\0&0&0&1&...&0\\...&...&...&...&1&...\\0&0&0&0&...&1\end{vmatrix} \]

即主对角线上的数字为 \(1\)，其余为 \(0\)。

不难验证其正确性。

【主要思想】

1.【从斐波那契说起】

题目放这里

著名的斐波那契数列有 \(F_i=F_{i-1}+F_{i-2}\) 的递推式，求第 \(n\) 项的时间复杂度为 \(O(n)\)。

但是如果结合矩阵乘法，会有意想不到的优化。

2.【构造状态矩阵】

要得到 \(F_i\)，显然只需要知道 \(F_{i-1},F_{i-2}\) 即可，所以每次只需要保留最后两位。

设 \(F(n)\) 表示一个 \(1\times 2\) 的矩阵，\(F(n)=[Fib_n\ Fib_{n+1}]\)。

我们希望根据 \(F(n-1)=[Fib_{n-1}\ Fib_{n}]\) 得到 \(F(n)\)，这需要另一个矩阵的辅助。

3.【构造转移矩阵】

不难发现我们需要将 \(F(n-1)\) 的第一项变为 \(Fib_{n-1}+Fib_{n}\)，第二项变为 \(Fid_{n}\)。

所以不难构造出一个转移矩阵：

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

根据上述乘法原则，不难发现 \(F(n)=F(n-1)\times A\)。

4.【矩阵快速幂加速递推】

假设我们需要第 \(F(n)\)，设 \(F(0)=[0\ 1]\)，显然 \(F(n)=F(0)\times A^n\)。

所有具有乘法结合律的乘方运算都可以用快速幂来加速，矩阵同理。

不仅如此，快速幂是贯穿矩阵学习的重要思想，请不了解的同学寻找相关资料学习。

5.【Code】

代码大概长酱紫：

void Mul(int f[2],int a[2][2]){
    int c[2];
    memset(c, 0, sizeof(c));
    for(int j=0;j<2;j++)
        for(int k=0;k<2;k++)
            c[j] = (c[j] + 1LL * f[k] * a[k][j] % MOD) % MOD;
    memcpy(f, c, sizeof(c));
}

void Mulself(int a[2][2]){
    int c[2][2];
    memset(c, 0, sizeof(c));
    for(int i=0; i<2; i++)
        for(int j=0; j<2; j++)
            for(int k=0; k<2; k++)
                c[i][j] = (c[i][j] + 1LL * a[i][k] * a[k][j] % MOD) % MOD;
    memcpy(a, c, sizeof(c));
}

int f[2]={0,1};
int a[2][2]={{0,1},{1,1}};
for(; n; n >>= 1){
    if(n & 1) Mul(f, a);
    Mulself(a);
}
printf("%lld\n", f[0]);

6.【小结】

通过这道题，我们对矩阵乘法加速递推有了初步了解。

一般来说，此类题目有以下特点：

1. 可以抽象出一个长度为 \(n\) 的一维向量，该向量在每个单位时间内只发生一次变化。
2. 变化形式为线性递推。（加上若干数或乘上一个数）
3. 该递推式规则本身保持不变。（根据这个构造转移矩阵）
4. 递推轮数很多但向量长度 \(n\) 不大。

那么我们可以考虑矩阵乘法优化，这个长度为 \(n\) 的矩阵称为状态矩阵，上面出现过的 \(A\) 称为转移矩阵。

解题的关键就是构造出合适的转移矩阵和状态矩阵。

转移矩阵一般的构造规则是：

如果状态矩阵中的第 \(x\) 个数对下一时间单位的第 \(y\) 个数有影响，那么将 \(A_{x,y}\) 赋为恰当的值。

时间复杂度为 \(O(n^3\log T)\)，\(T\) 为递推总轮数。

【经典例题】

【典例一】

【模板】矩阵加速（数列）

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

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

求 \(a\) 数列的第 \(n\) 项。

不难想到推导过程。

可以定义状态矩阵：\(F(n)=[a_n,a_{n-1},a_{n-2}]\)。

显然 \(F(n+1)=[a_n+a_{n-2},a_n,a_{n-1}]\)。

所以有：

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

主要代码如下：

void Mul(LL f[3], LL a[3][3]){
    LL c[3];
    memset(c, 0, sizeof(c));
    for(int j=0; j<3; j++)
        for(int k=0; k<3; k++)
            c[j] = (c[j] + (f[k] * a[k][j]) % MOD) % MOD;
    memcpy(f, c, sizeof(c));
}

void Mulself(LL a[3][3]){
    LL c[3][3];
    memset(c, 0, sizeof(c));
    for(int i=0; i<3; i++)
        for(int j=0; j<3; j++)
            for(int k=0; k<3; k++)
                c[i][j] = (c[i][j] + (a[i][k] * a[k][j]) % MOD) % MOD;
    memcpy(a, c, sizeof(c));
}

int main(){
    int T=read();
    while(T--){
        n=read();
        if(n<=3) {puts("1"); continue;}
        LL f[3]={1,1,1};
        LL a[3][3]={{1,1,0},{0,0,1},{1,0,0}};
        n -= 3;
        for(; n; n >>= 1){
            if(n & 1) Mul(f, a);
            Mulself(a);
        }
        printf("%lld\n", f[0]);
    }
    return 0;
}

【典例二】

GTY's birthday gift

给定大小为 \(n\) 的正整数可重集 \(S\)，每次操作可以加入数 \(a+b(a,b\in S)\)，求 \(k\) 次操作后 \(S\) 的和的最大值。

显然每次贪心的寻找最大的两个数相加，关键是求和。

这需要状态矩阵加上一维表示和。

令 \(F(n)=[f(n-1),f(n),S(n)=\sum_{i=1}^n f(i)]\)（\(f(n),f(n-1)\) 分别表示当前集合内的最大值和次大值）

显然 \(F(n+1)=[f(n),f(n-1)+f(n),S(n)+f(n-1)+f(n)]\)。

则有：

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

主要代码如下：

int main(){
    while(~scanf("%d %d", &n, &k)){
        LL ans=0;
        for(int i=1; i<=n; i++) 
            p[i]=read(),ans=(ans+p[i]) % MOD;
        sort(p+1, p+n+1);
        LL f[3]={p[n-1], p[n], ans};
        LL a[3][3]={{0, 1, 1},{1, 1, 1},{0, 0, 1}};
        for(; k; k >>= 1){
            if(k & 1) Mul(f, a);
            Mulself(a);
        }
        printf("%lld\n", f[2]);
    }
    return 0;
}

【典例三】

Power of Matrix

给定矩阵 \(A\) 和 \(k\)，求 \(A+A^2+A^3+...+A^k\)。

关于矩阵加法的定义已经给出，不再赘述。

直接计算的复杂度为 \(O(n^3k)\)，显然不可取。

考虑分治，每次将序列二分并分治下去，计 \(F(k)=A+A^2+A^3+...+A^k\)。

显然对于每个 \(k\)：

\[F(k)=\begin{cases} F(\frac{k}{2})+A^{\frac{k}{2}}\times F(\frac{k}{2}) & k\mod 2=0\\ F(\frac{k}{2})+A^{\frac{k}{2}}\times F(\frac{k}{2})+A^k & k\mod 2=1 \end{cases} \]

将矩阵乘法在结构体中重载，代码还挺好写的，之后都采用这种方式，而且代码不再重复给出。

struct Matrix{
    int a[50][50];
    int MOD = 10; 

    Matrix() {memset(a, 0, sizeof(a));}

    Matrix operator * (const Matrix &b){
        Matrix c;
        for(int i=0; i<n; i++)
            for(int j=0; j<n; j++)
                for(int k=0; k<n; k++)
                    c.a[i][j] = (c.a[i][j] + a[i][k] * b.a[k][j] % MOD) % MOD;
        return c;
    }

    Matrix operator + (const Matrix &b){
        Matrix c;
        for(int i=0; i<n; i++)
            for(int j=0; j<n; j++)
                c.a[i][j] = (a[i][j] + b.a[i][j]) % MOD;
        return c;
    }

    Matrix operator ^ (int k){
        Matrix c, b;
        for(int i=0; i<n; i++) c.a[i][i]=1;//这里用到了单位矩阵
        memcpy(b.a, a, sizeof(a));
        for(; k; k >>= 1){
            if(k & 1) c = c * b;
            b = b * b;
        }
        return c;
    }
} A;

Matrix dfs(int k){
    if(k == 1) return A;
    Matrix P = dfs(k >> 1);
    P = P + P * (A ^ (k >> 1));
    if(k & 1) P = P + (A ^ k);
    return P;
}

【典例四】

石头游戏

太长了，自己看

同样将问题转化为一维向量。

定义长度为 \(n\times m+1\) 的状态矩阵 \(F\)，下标为 \(0\) 到 \(n\times m\)。

记 \(num(i,j)=(i-1)\times m+j\)。

用 \(F[num(i,j)]\) 记录格子 \((i,j)\) 中的石头数量。

特别的，令 \(F[0]=1\)，把它当做得到石头的窗口，原因后文介绍。

设 \(F(k)\) 为 \(k\) 秒后的状态矩阵，显然 \(F(0)=[1,0,0,...,0]\)。

可惜本题的转移矩阵是不断变化的，而且每个格子还不一样。

但是可以发现 \(lcm(1,2,3,4,5,6)=60\)，所以我们可以记录下前 \(60\) 个状态矩阵。

设 \(A_k\) 表示第 \(k\) 个，构造方法为：

1. 若第 \(k\) 秒 \((i,j)\) 的操作字符为 \(N\)，则令 \(A_k[num(i,j),num(i-1,j)]=1\)。其余同理。
2. 若第 \(k\) 秒 \((i,j)\) 的操作字符为数字 \(x\)，则令 \(A_k[0,num(i,j)]=x\)。
3. 令 \(A_k[0,0]=1\)。
4. 令其余为 \(0\)。

注意 \(F[0]=1\)，显然以上矩阵可以胜任我们的状态转移工作。

设 \(A=\Pi_{i=1}^{60} A_i,t=60\times k+b\)，那么 \(F(t)=F(0)\times A^k\times \Pi_{i=1}^b A_i\)，那么：

\[ans=\max_{1\leq i\leq n\times m}\{F(i)\} \]

主要代码如下：

int main(){
    scanf("%d %d %d %d", &n, &m, &t, &act);
    for(int i=1; i<=n; i++){
        scanf("%s", s + 1);
        for(int j=1; j<=m; j++) 
            a[i][j] = s[j] - '0' + 1;
    }
    for(int i=1; i<=act; i++) scanf("%s", p[i]);
    for(int k=1; k<=60; k++){
        for(int i=1; i<=n; i++)
            for(int j=1; j<=m; j++){
                int x = a[i][j], y = now[i][j];
                now[i][j] = (y+1) % strlen(p[x]);
                if(p[x][y] >= '0' && p[x][y] <= '9'){
                    f[k].a[0][Num(i, j)] = p[x][y] - '0';
                    f[k].a[Num(i, j)][Num(i, j)] = 1;
                }
                else if(p[x][y] == 'N' && i > 1) f[k].a[Num(i, j)][Num(i-1, j)] = 1;
                else if(p[x][y] == 'S' && i < n) f[k].a[Num(i, j)][Num(i+1, j)] = 1;
                else if(p[x][y] == 'E' && j < m) f[k].a[Num(i, j)][Num(i, j+1)] = 1;
                else if(p[x][y] == 'W' && j > 1) f[k].a[Num(i, j)][Num(i, j-1)] = 1;
            }
        f[k].a[0][0] = 1;
    }
    Matrix P = f[1];
    for(int i=2; i<=60; i++) P = P * f[i];
    
    LL ans[100];
    memset(ans, 0, sizeof(ans));
    ans[0] = 1;
    
    int k = t / 60;
    for(; k; k >>= 1){
        if(k & 1) Mul(ans, P.a);
        P = P * P;
    }

    k = t % 60;
    for(int i=1; i<=k; i++) Mul(ans, f[i].a);
    LL mx = 0;
    for(int i=1; i<=n*m; i++) mx = max(mx, ans[i]);
    printf("%lld\n",mx);
    return 0;
}

【典例五】

How many ways??

询问有向图上从 \(s\) 点恰好经过 \(k\) 个点到达 \(t\) 点的方案数

经典的图上问题转化为矩阵问题。

首先将图转化为邻接矩阵（这注定了 \(n\) 的范围不会很大，这也是此类题目的标志之一）

令 \(C=A\times A\)，那么 \(C[i][j]=\sum A[i][k]×A[k][j]\) ，实际上就等同于 \(i\) 到 \(j\) 恰好经过两条边的路径数量。

\(k\) 步同理。

其实相当于 Floyd 算法的矩阵乘法实现。

而利用快速幂相当于 Floyd 算法的倍增实现。

代码如下：

Matrix Pow(Matrix a, int k){
    Matrix b;
    for(int i=0; i<n; i++) b.a[i][i] = 1;//单位矩阵
    for(; k; k >>= 1){
        if(k & 1) b = b * a; 
        a = a * a;
    }
    return b;
}

int main(){
    while(~scanf("%d %d", &n, &m) && n){
        Matrix A;
        for(int i=1; i<=m; i++){
            int u=read(), v=read();
            A.a[u][v] = 1;
        }
        int T=read();
        while(T--){
            int u=read(), v=read(), k=read();
            Matrix ans = Pow(A, k);
            printf("%d\n", ans.a[u][v]);
        }
    }
    return 0;
}

【典例六】

Cow Relays G

给定一张 \(T\) 条边的无向连通图，求从 \(S\) 到 \(E\) 经过 \(N\) 条边的最短路长度。

从路径个数到最短路长度，其实没变多少。

矩阵乘法时将乘法变为取 min，其实这更像 Floyd 了。

主要思想依然是用矩阵来优化。

同时注意离散化，主要代码如下：

struct Matrix{
	int a[210][210];
	Matrix() {memset(a, 0x3f, sizeof(a));}
	
	Matrix operator * (const Matrix &b){
		Matrix c;
		for(int i=1; i<=n; i++)
			for(int j=1; j<=n; j++)
				for(int k=1; k<=n; k++)
					c.a[i][j] = min(c.a[i][j], a[i][k] + b.a[k][j]);
		return c;
	}
};
int main(){
    k=read(), m=read(), s=read(), t=read();
    memset(id, 0, sizeof(id));
    n = 0;
    Matrix A;
    for(int i=1; i<=m; i++){
        int w=read(), u=read(), v=read();
        u = id[u] ? id[u] : id[u] = ++n;
        v = id[v] ? id[v] : id[v] = ++n;
        A.a[u][v] = A.a[v][u] = min(A.a[u][v], w);
    }
    s = id[s], t = id[t];
    Matrix ans = Pow(A, k);
    printf("%d\n", ans.a[s][t]);
    return 0;
}

【典例七】

魔法值

太长了，自己看

看到 \(n\leq 100\) 和这个转移方式，矩阵不难想到。

显然根据异或的性质，异或偶数次等于没有异或。

所以我们只需要统计每个点 \(u\) 对 \(1\) 号点的影响次数，如果影响次数为偶数次就异或上，然后得到答案。

我们可以用一个 \(0/1\) 矩阵每个位置 \((i,j)\) 表示节点 \(i\) 是否对节点 \(j\) 有贡献，利用矩阵快速幂加速计算。

可惜这样我们每次都要乘一次，时间复杂度：\(O(qn^3\log a)\)，好像不太行。

考虑倍增优化。

假设初始矩阵为 \(A\)，我们可以预处理出 \(A,A^2,A^4,..A^{2^{32}}\) 相当于对其进行了二进制拆分。

时间复杂度变为 \(O(n^3\log a+qn^2\log a)\)，好像可以。

主要代码：

struct Matrix{
    LL a[N][N];
	Matrix() {memset(a, 0, sizeof(a));}
	
	Matrix operator * (const Matrix &b){
		Matrix c;
		for(int i=1; i<=n; i++)
			for(int j=1; j<=n; j++)
				for(int k=1; k<=n; k++)
					c.a[i][j] ^= (a[i][k] & b.a[k][j]);
		return c;
	}
}f[35];

void Mul(int k){
    LL tmp[N];
    memset(tmp, 0, sizeof(tmp));
    for(int i=1; i<=n; i++)
        for(int j=1; j<=n; j++)
            if(f[k].a[i][j]) 
                tmp[i] = tmp[i] ^ ans[j];
    memcpy(ans, tmp, sizeof(tmp));
}

int main(){
    n=read(), m=read(), q=read();
    for(int i=1; i<=n; i++) p[i] = read();
    for(int i=1; i<=m; i++){
        int u=read(), v=read();
        f[0].a[u][v] = f[0].a[v][u] = 1;
    }
    for(int i=1; i<=32; i++)
        f[i] = f[i-1] * f[i-1];
    while(q--){
        LL k=read();
        for(int i=1; i<=n; i++) ans[i] = p[i];
        for(LL i=0; i<=31; i++)
            if(k & (1<<i)) Mul(i);
        printf("%lld\n", ans[1]);
    }
    return 0;
}

【总结】

根据以上例题的学习总结，我们完成了矩阵从基本运算到数列递推到图上问题的完美转化。

值得注意的是，根据典例二、四，我们发现在矩阵乘法加速递推遇到求和或常数项时，通常有一下方法：

在状态矩阵中添加一个额外位置，始终储存数列和或常数项 \(1\)，并乘上转移矩阵中的适当系数。

通过这种方法将系数累加到状态矩阵的相应位置。

根据典例六我们发现，利用矩阵乘法求解 Floyd 时并不需要像 Floyd 那样先枚举中间节点。

因为这本质上是矩阵乘法，满足交换律，而普通的 Floyd 算法基于 DP 实现，有严格的阶段性。

至此矩阵乘法入门的全部内容已整理完毕。

引用资料&特别鸣谢：

1. 矩阵乘法优化DP复习。
2. 【学习笔记】动态规划—矩阵递推加速。
3. 《算法竞赛进阶指南》。

完结撒花