矩阵

定义

矩阵是一个按照长方阵列排列的复数或实数集合。

一个矩阵应当形如:

A=[a1,1 a1,2 a1,na2,1 a2,2 a2,n   am,1  am,2   am,n]A= \begin{bmatrix} a_{1,1}\ &a_{1,2}\ & \cdots &a_{1,n}\\ a_{2,1}\ &a_{2,2}\ & \cdots &a_{2,n}\\ \vdots\ &\vdots\ &\ddots\ &\vdots\\ a_{m,1}\ &\ a_{m,2}\ &\cdots\ &\ a_{m,n} \end{bmatrix}

n=mn=m 时,这个矩阵就变成了方阵。

特别的,我们定义 A0A^0 为单位矩阵

I=[1 0 00 1 0   0  0   1]I= \begin{bmatrix} 1\ &0\ & \cdots &0\\ 0\ &1\ & \cdots &0\\ \vdots\ &\vdots\ &\ddots\ &\vdots\\ 0\ &\ 0\ &\cdots\ &\ 1 \end{bmatrix}


(为了源代码好看,就这样吧。)

运算

矩阵中一般定义有矩阵乘法和矩阵加法。

由于日常主要使用的是矩阵乘法,所以这里只讲矩阵乘法。

矩阵乘法

在进行矩阵乘法的时候,需要满足一个前提条件,即矩阵 AA 的行数等于矩阵 BB 的列数。

两个大小分别为 m×nm\times nn×pn\times p 的矩阵 A,BA,B 相乘后是一个大小为 m×pm\times p 的矩阵 CC,满足:Ci,j=k=1nai,kbk,j(1im,1jp)C_{i,j}=\sum\limits_{k=1}^na_{i, k}\cdot b_{k,j}(1\le i\le m, 1\le j\le p)

另外,对于单位矩阵 II,有 A×I=AA\times I=A

矩阵乘法满足分配律和结合律,但不满足交换律。即:

  • A(BC)=(AB)CA(BC)=(AB)C
  • (A+B)C=AC+BC(A+B)C=AC+BC
  • ABBAAB\ne BA

洛谷 B2105 矩阵乘法

模板题,按照定义完成即可。

实例代码:

#include <bits/stdc++.h>

using namespace std;
const int MAXN = 105;
int n, m, k, a[MAXN][MAXN], b[MAXN][MAXN], c[MAXN][MAXN];
signed main(){
  ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
  cin >> n >> m >> k;
  for (int i = 1; i <= n; i++){
    for (int j = 1; j <= m; j++){
      cin >> a[i][j];
    }
  }
  for (int i = 1; i <= m; i++){
    for (int j = 1; j <= k; j++){
      cin >> b[i][j];
    }
  }
  for (int i = 1; i <= n; i++){
    for (int j = 1; j <= k; j++){
      for (int p = 1; p <= m; p++){
        c[i][j] += a[i][p] * b[p][j];
      }
      cout << c[i][j] << " \n"[j == k];
    }
  }
  return 0;
}

矩阵快速幂

由于矩阵乘法满足结合律,即 A(BC)=(AB)CA(BC)=(AB)C,那么我们就可以利用类似快速幂的思想来求。

直接放代码(题目链接):

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

using namespace std;
const int  Mod = 1e9 + 7;
int n, k;
vector<vector<int>> a(105, vector<int>(105)), ans(105, vector<int>(105));
vector<vector<int>> mul(vector<vector<int>> a, vector<vector<int>> b){
  vector<vector<int>> ans(105, vector<int>(105));
  for (int i = 1; i <= n; i++){
    for (int j = 1; j <= n; j++){
      for (int k = 1; k <= n; k++){
        (ans[i][j] += a[i][k] * b[k][j] % Mod) %= Mod;
      }
    }
  }
  return ans;
}
void power(vector<vector<int>> a, int b){
  for (int i = 1; i <= n; i++){
    ans[i][i] = 1;
  }
  for (; b; b >>= 1, a = mul(a, a)){
    if (b & 1){
      ans = mul(ans, a);
    }
  }
  for (int i = 1; i <= n; i++){
    for (int j = 1; j <= n; j++){
      cout << ans[i][j] << " \n"[j == n];
    }
  }
}
signed main(){
  ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
  cin >> n >> k;
  for (int i = 1; i <= n; i++){
    for (int j = 1; j <= n; j++){
      cin >> a[i][j];
    }
  }
  power(a, k);
  return 0;
}

矩阵加速

我们直接来看例题:洛谷 P1939 矩阵加速(数列)

题意:已知数列 aa,且满足

ax={1x{1,2,3}ax1+ax3x4a_x= \begin{cases} 1 & x \in\{1,2,3\} \\ a_{x-1}+a_{x-3} & x \geq 4 \end{cases}

an(1n2×109)a_n(1\le n\le 2\times 10^9)109+710^9+7 取模后的值。

由于 nn 过大,显然不能直接求,只能使用矩阵加速。

我们来进行分析。

根据题意可知以下几点(我们假设有 n5n\ge 5):

  • an=an1+an3a_n=a_{n-1}+a_{n-3}
  • an1=an2+an4a_{n-1}=a_{n-2}+a_{n-4}
  • an+1=an1+an2a_{n+1}=a_{n-1}+a_{n-2}

可以发现,有一些数值对我们的答案有比较重要的贡献,即 an,an1,an2,an3a_n,a_{n-1},a_{n-2},a_{n-3}

由于要求的是 ana_n,那么我们将 an,an1,an2a_n,a_{n-1},a_{n-2} 放一起讨论,an1,an2,an3a_{n-1},a_{n-2},a_{n-3} 放一起讨论。(感觉说的比较抽象,感性理解一下。)

列一个表:

an1a_{n-1} an2a_{n-2} an3a_{n-3}
ana_n
an1a_{n-1}
an2a_{n-2}

首先,有 an=1×an1+0×an2+1×an3a_n=1\times a_{n-1}+0\times a_{n-2}+1\times a_{n-3}

同时也很容易得到:an1=1×an1+0×an2+0×an3a_{n-1}=1\times a_{n-1}+0\times a_{n-2}+0\times a_{n-3}an2=0×an1+1×an2+0×an3a_{n-2}=0\times a_{n-1}+1\times a_{n-2}+0\times a_{n-3}

那么这张表就可以长成这样:

an1a_{n-1} an2a_{n-2} an3a_{n-3}
ana_n 1 0 1
an1a_{n-1} 1 0 0
an2a_{n-2} 0 1 0

那么我们的基础矩阵就出来了:

A=[1 1 01 0 00  1 0]A= \begin{bmatrix} 1\ &1\ & 0\\ 1\ &0\ & 0\\ 0\ &\ 1\ &0 \end{bmatrix}

然后就矩阵快速幂就好了。

代码:

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

using namespace std;
const int Mod = 1e9 + 7;
struct Node{
  int a[4][4];
  Node operator*(const Node &x)const{
    Node ans;
    for (int i = 1; i <= 3; i++){
      for (int j = 1; j <= 3; j++){
        ans.a[i][j] = 0;
        for (int k = 1; k <= 3; k++){
          (ans.a[i][j] += a[i][k] * x.a[k][j] % Mod) %= Mod;
        }
      }
    }
    return ans;
  }
}ans, base;
void Solve(){
  for (int i = 1; i <= 3; i++){
    for (int j = 1; j <= 3; j++){
      base.a[i][j] = ans.a[i][j] = 0;
    }
  }
  base.a[1][3] = base.a[1][1] = base.a[2][1] = base.a[3][2] = 1;
  ans.a[1][1] = ans.a[2][2] = ans.a[3][3] = 1;
  int x;
  cin >> x;
  for (; x; x >>= 1, base = base * base){
    if (x & 1){
      ans = ans * base;
    }
  }
  cout << ans.a[2][1] << '\n';
}
signed main(){
  ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
  int t;
  cin >> t;
  while (t--){
    Solve();
  }
  return 0;
}

例题

P2886 [USACO07NOV] Cow Relays G

题意:给定一张 T(2T100)T(2\le T\le 100) 条边的无向连通图,求从 SSEE 经过 N(1N106)N(1\le N\le 10^6) 条边的最短路长度。

由于出现过的点最多只有 2×T2\times T 个,所以我们可以把这些点离散化下来,下记离散化后共有 nn 个点。

但是由于 NN 很大,无法直接暴力计算答案,所以我们继续分析。

我们注意到矩阵乘法的转移 Ci,j=k=1nai,kbk,j(1im,1jp)C_{i,j}=\sum\limits_{k=1}^na_{i, k}\cdot b_{k,j}(1\le i\le m, 1\le j\le p) 和 Floyd 的转移 dpx,y=min(dpx,y,dpx,k+dpk,y)dp_{x,y}=\min(dp_{x,y},dp_{x,k}+dp{{k,y}}) 的转移很像,所以考虑能不能使用矩阵加速。

我们注意到每做一次 Floyd,经过的边的数量就会翻倍(感性理解)。那么由于所有转移互相独立,具有结合律,那么就可以直接使用矩阵快速幂解决了。

时间复杂度:O(T3logN)O(T^3\log N)

空间复杂度:O(T2+N)O(T^2+N)

代码:

#include<bits/stdc++.h>

using namespace std;
const int MAXN = 1e6 + 5, MAXT = 100 + 5;
int n, t, s, e, cnt, num[MAXN];
struct Node{
  int dis[MAXT][MAXT];
  Node operator*(const Node &x)const{
    Node ans;
    for (int i = 1; i <= cnt; i++){
      for (int j = 1; j <= cnt; j++){
        ans.dis[i][j] = 1e9;
      }
    }
    for (int i = 1; i <= cnt; i++){
      for (int j = 1; j <= cnt; j++){
        for (int k = 1; k <= cnt; k++){
          ans.dis[i][j] = min(ans.dis[i][j], dis[i][k] + x.dis[k][j]);
        }
      }
    }
    return ans;
  }
}dis, ans;
Node power(){
  n--;
  ans = dis;
  for (; n; n >>= 1, dis = dis * dis){
    if (n & 1){
      ans = ans * dis;
    }
  }
  return ans;
}
void Solve(){
  for (int i = 1; i < 105; i++){
    for (int j = 1; j < 105; j++){
      dis.dis[i][j] = 1e9;
    }
  }
  cin >> n >> t >> s >> e;
  while (t--){
    int u, v, w;
    cin >> w >> u >> v;
    if (!num[u]){
      num[u] = ++cnt;
    }
    if (!num[v]){
      num[v] = ++cnt;
    }
    dis.dis[num[u]][num[v]] = dis.dis[num[v]][num[u]] = w;
  }
  cout << power().dis[num[s]][num[e]];
}
signed main(){
  ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
  int t = 1;
  while (t--){
    Solve();
  }
  return 0;
}