【学习笔记】矩阵乘法 / 矩阵快速幂的应用

Basic Tips

• 不同角度看待线性方程组：矩阵

  对于一个方程组：

  $$\begin{cases} 4x_1 + 16x_2 = 114 \\ 3x_1 + 3x_2 = 514 \\ 7x_1 + 3x_2 = 416 \\ \end{cases}$$

  利用一种打包处理的思想将其写成矩阵的形式：
  $\begin{bmatrix} x_1 & x_2 \end{bmatrix} \cdot \begin{bmatrix}4 & 16 \\\\ 3 & 3 \\\\ 7 & 3 \end{bmatrix} = \begin{bmatrix}114 \\\\ 514 \\\\ 416 \end{bmatrix}$

  矩阵的存储一般使用二维数组。

• 矩阵数乘 / 矩阵加减

  对于一个矩阵 $A$，让它乘上一个数 $x$，结果便是 $A$ 中每个位置的元素乘上 $x$。

  矩阵加减只能在大小相同的两个矩阵之间进行，例如大小相同两个矩阵 $A$ 跟 $B$ 相加 / 减，得到的结果矩阵即为两矩阵对应位置上的元素相加 / 减。

• 矩阵乘法 / 矩阵快速幂

  对于大小为 $X \times K$ 的矩阵 $A$ 和大小为 $K \times Y$ 的矩阵 $B$，记它们的乘积矩阵为 $C$，可得 $C_{i,j} = \displaystyle\sum^{K}_{k=1}{A_{i,k} \times B_{k,j}}$，注意，矩阵乘法并不满足 交换律，矩阵乘法可以快速的处理一些 $ax + by = c$ 的关系，需要关注枚举顺序，矩阵快速幂通常用来加速递推、优化 dp 转移，在图论中根据邻接矩阵可以求限制路径长度的路径方案数，复杂度取决于转移矩阵的大小，一般是 $O(\log n)$ 带常数。

• 矩阵的运算律

  矩阵并不满足一般的交换律，即 $A \times B$ 不等同于 $B \times A$，结果可能不一样，矩阵满足结合和左右分配律，即 $A \cdot (BC) = (AB) \cdot C$。

• 矩阵快速幂的实现

  相较于数的快速幂，矩阵快速幂要初始化为单位矩阵。

  #include<bits/stdc++.h>
  #define int long long
  using namespace std;
  const int mod = 1e9 + 7;
  int n,k;
  
  struct node{
      int mat[105][105];
      node(){memset(mat,0,sizeof(mat));}
      inline void init(){for(int i = 1;i <= n;i ++) mat[i][i] = 1;}
  }x; // 结构体分装会更方便
  
  node operator *(const node &x,const node &y)
  {
      node z;
      for(int k = 1;k <= n;k ++)
      {
          for(int i = 1;i <= n;i ++)
          {
              for(int j = 1;j <= n;j ++)
                  z.mat[i][j] = (z.mat[i][j] + x.mat[i][k] * y.mat[k][j] % mod) % mod;
          }
      }
      return z;
  } // 重定义乘号运算
  
  inline node qpow(node x,int k)
  {
      node ret;
      ret.init();
      while(k)
      {
          if(k & 1) ret = ret * x;
          x = x * x;
          k >>= 1;
      }
      return ret;
  } // 初始化单位矩阵
  
  signed main()
  {
      scanf("%lld %lld",&n,&k);
      for(int i = 1;i <= n;i ++)
      {
          for(int j = 1;j <= n;j ++)
              scanf("%lld",&x.mat[i][j]);
      }
      node Ans = qpow(x,k);
      for(int i = 1;i <= n;i ++)
      {
          for(int j = 1;j <= n;j ++)
              printf("%lld%c",Ans.mat[i][j]," \n"[j == n]);
      }
      return 0;
  }
  

• What's More

  单位矩阵 $I$：除对角线元素为 $1$ 外其余元素皆为 $0$。

  小细节：实现矩阵比较方便的是使用结构体分装，通过构造函数进行初始化操作，部分存在结合律的运算也可以通过变形的矩阵乘法实现，所以对于题目中出现的 $01$ 异或、gcd、min、max 等均可以往这方面思考。

Examples

P1962 斐波那契数列

对于 $n \leq 2^{63}$ 的数据，正常的递推肯定会爆，因此我们来观察转移式子：

$\large \begin{aligned} & 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} \end{aligned}$

这就是上文所提到的线性方程组的形式，写成矩阵：

$\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} f_{n - 1} \\ f_{n - 2} \end{bmatrix} = \begin{bmatrix} f_n \\ f_{n - 1} \end{bmatrix}$

求 $f_n$ 相当于转移矩阵 $A^n$ 乘上初始矩阵。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod = 1e9 + 7;
struct node{
	int mz[3][3];
	node operator * (const node& x) const
	{
		node ret;
		for(int i = 1;i <= 2;i ++)
		{
			for(int j = 1;j <= 2;j ++)
			{
				ret.mz[i][j] = 0;
				for(int k = 1;k <= 2;k ++)
					ret.mz[i][j] = ((ret.mz[i][j] + mz[i][k] * x.mz[k][j] % mod) % mod + mod) % mod;
			}
		}
		return ret;
	}
}feb;
int n;
inline node qpow(node x,int p)
{
	node res = x;
	p --;
	while(p)
	{
		if(p & 1)
			res = res * x;
		x = x * x;
		p >>= 1;
	}
	return res;
}
signed main()
{
	cin.tie(0) -> sync_with_stdio(0);
	cout.tie(0) -> sync_with_stdio(0);
	cin >> n;
	feb.mz[1][1] = feb.mz[1][2] = feb.mz[2][1] = 1;
	if(n <= 2)
		cout << "1" << endl;
	else
	{
		node Answer = qpow(feb,n);
		cout << Answer.mz[1][2] << endl;
	}
	return 0;
}

P2886 [USACO07NOV] Cow Relays G

图论上的应用，鉴于 $T \leq 100$，可以邻接矩阵存图，由于是求最短路，可以用 Floyd，发现可以类比矩阵乘法，记邻接矩阵为 $e$，将矩阵乘法转换成 $\text{Floyd}$ 的转移，答案即为 $e^N_{s,e}$，注意处理重边。

#include <bits/stdc++.h>
// #define int long long

using namespace std;
template <typename Tp>
inline void read(Tp &x){
    Tp res = 0,f = 1;
    char ch = getchar();
    while(!isdigit(ch)){if(ch == '-') f = -1;ch = getchar();}
    while(isdigit(ch)){res = (res << 1) + (res << 3) + (ch ^ 48);ch = getchar();}
    x = res * f;
}

const int maxn = 1005,inf = 0x3f3f3f3f;
int n,t,s,e,apr[maxn],idx;

struct Matrix{
    int mat[maxn][maxn];
    Matrix(){memset(mat,inf,sizeof(mat));}
    Matrix operator * (const Matrix &g) const{
        Matrix ret;
        for(int k = 1;k <= idx;k ++){
            for(int i = 1;i <= idx;i ++){
                for(int j = 1;j <= idx;j ++){
                    ret.mat[i][j] = min(ret.mat[i][j],g.mat[i][k] + mat[k][j]);
                }
            }
        }
        return ret;
    }
}Ans;

Matrix mat_qpow(Matrix x,int k){
    Matrix ret;
    ret = x,k --;
    while(k){
        if(k & 1) ret = ret * x;
        x = x * x;
        k >>= 1;
    }
    return ret;
}

signed main(){
    read(n),read(t),read(s),read(e);
    for(int i = 1;i <= t;i ++){
        int w,u,v;
        read(w),read(u),read(v);
        if(!apr[u]) apr[u] = ++ idx;
        if(!apr[v]) apr[v] = ++ idx;
        Ans.mat[apr[u]][apr[v]] = Ans.mat[apr[v]][apr[u]] = w; // 处理重边
    }
    Ans = mat_qpow(Ans,n);
    printf("%d\n",Ans.mat[apr[s]][apr[e]]);
    return 0;
}

P6569 [NOI Online #3 提高组] 魔法值

由于原图的边权均为 $0$ 或 $1$，因此题面中的式子可以变成：

$\begin{aligned} f_{x,i} = f_{1,i - 1} \times e_{1,x} \oplus f_{2,i - 1} \times e_{2,x} \oplus \dots f_{n,i - 1} \times e_{n,x} \end{aligned}$

啊，这不就是异或版矩乘？

记第 $i$ 天的魔法值矩阵为 $f_i$，邻接矩阵为 $e$，由于邻接矩阵有且仅有 $01$，因此它满足矩乘结合律，得到以下式子：

$\begin{aligned} f_{i} = f_{i - 1} \times e \to f_i = f_0 \times e^i \end{aligned}$

如果每次询问都跑一次快速幂一定会爆，因此看到数据范围 $f_i,a_i < 2^{32}$，不难想到类似状压的方法，预处理邻接矩阵的 $2^i$ 即可，复杂度 $O(31(n^3 + q))$

#include <bits/stdc++.h>
#define int unsigned int

using namespace std;
template <typename Tp>
inline void read(Tp &x){
    Tp res = 0,f = 1;
    char ch = getchar();
    while(!isdigit(ch)){if(ch == '-') f = -1;ch = getchar();}
    while(isdigit(ch)){res = (res << 1) + (res << 3) + (ch ^ 48);ch = getchar();}
    x = res * f;
}

const int maxn = 105;
int n,m,q,f0[maxn],x;

struct Matrix{
    int mat[maxn][maxn];
    Matrix(int x = 0){
        memset(mat,0,sizeof(mat));
        for(int i = 1;i <= n;i ++) mat[i][i] = x;
    }
    Matrix operator * (const Matrix &g) const{
        Matrix ret(0);
        for(register int k = 1;k <= n;k ++){
            for(register int i = 1;i <= n;i ++){
                for(register int j = 1;j <= n;j ++){
                    ret.mat[i][j] ^= mat[i][k] * g.mat[k][j];
                }
            }
        }
        return ret;
    }
};

Matrix d[33];

signed main(){
    read(n),read(m),read(q);
    for(register int i = 1;i <= n;i ++) read(f0[i]);
    for(register int i = 1;i <= m;i ++){
        int u,v;
        read(u),read(v);
        d[0].mat[u][v] = d[0].mat[v][u] = 1;
    }
    for(register int i = 1;i < 32;i ++)
        d[i] = d[i - 1] * d[i - 1];
    while(q --){
        read(x);
        Matrix tmp;
        for(register int i = 1;i <= n;i ++)
            tmp.mat[1][i] = f0[i];
        for(register int i = 0;i < 32;i ++){
            if((x >> i) & 1) tmp = tmp * d[i];
        }
        printf("%u\n",tmp.mat[1][1]);
    }
    return 0;
}

Summary

其实只要找到转移的式子，列出矩阵，其它的就可以自行实现了，尤其是 dp 方程，在列出暴力转移的情况下可以考虑矩阵优化，同时部分数据结构也可以维护矩阵。

trick：由于图论中用矩阵快速幂处理邻接矩阵需要保证边权 $\in \{0,1\}$，因此在边权的范围较小可以考虑拆点分层图。

Others

P1306 斐波那契公约数 重点是斐波那契相关定理的证明

P2151 [SDOI2009] HH去散步 点转边的 $\text{trick}$

P4159 [SCOI2009] 迷路 分层图即拆点

P6569 [NOI Online #3 提高组] 魔法值 矩乘变形成异或和

UVA Power of Matrix 矩阵套矩阵