矩阵 学习笔记
引入
矩阵的引入来自线性方程组，将其左边每一项的系数和右边的常数抽象出来就是矩阵。
一般用 \([]\) 或 \(()\) 表示矩阵。
定义
主对角线
对于一个矩阵 \(A\)，\(A_{i,i}\) 为该矩阵主对角线。
单位矩阵
一般用 \(I\) 表示，它的主对角线是 \(1\)，其余的是 \(0\)。
全零矩阵
对于一个矩阵，如果其所有元素均为 \(0\)，则称其为全零矩阵，记作 \(0\)。
同型矩阵
两个矩阵，行数与列数对应相同，称为同型矩阵。
转置矩阵
如果矩阵 \(A\) 和矩阵 \(B\) 满足：
- 矩阵 \(A\) 的行数等于矩阵 \(B\) 的列数，矩阵 \(A\) 的列数等于矩阵 \(B\) 的行数。
- 对于 \(\forall i,j, A_{i,j} = B_{j,i}\)。
则称 \(A\) 是 \(B\) 的转置矩阵，记作 \(A = B^T\)。
方阵
一个行数和列数相等的矩阵。
一个大小为 \(n\) 的矩阵叫作 \(n\) 阶矩阵或 \(n\) 阶方阵。
三角矩阵
如果方阵主对角线左下方的元素均为 \(0\)，称为上三角矩阵。如果方阵主对角线右上方的元素均为 \(0\)，称为下三角矩阵。
两个上（下）三角矩阵的乘积仍然是上（下）三角矩阵。如果对角线元素均非 \(0\)，则上（下）三角矩阵可逆，逆也是上（下）三角矩阵。
单位三角矩阵
如果上三角矩阵 \(A\) 的对角线全为 \(1\)，则称 \(A\) 是单位上三角矩阵。如果下三角矩阵 \(A\) 的对角线全为 \(1\)，则称 \(A\) 是单位下三角矩阵。
两个单位上（下）三角矩阵的乘积仍然是单位上（下）三角矩阵，单位上（下）三角矩阵的逆也是单位上（下）三角矩阵。
运算
线性运算
对两个同型矩阵可以作加（减）法，得到的答案是对应元素的和（差）。
矩阵的转置
即将一个矩阵变成它的转置矩阵。
矩阵乘法
如果矩阵 \(A\) 的行数等于矩阵 \(B\) 的列数，列数等于矩阵 \(B\) 行数，则这两个矩阵可以相乘，相乘的积矩阵 \(C\)，满足 \(C_{i,j} = \sum A_{i,k} \times B_{k,j}\)，即矩阵 \(C\) 的第 \(i\) 行第 \(j\) 列等于矩阵 \(A\) 的第 \(i\) 行乘上矩阵 \(B\) 的第 \(j\) 列。
可以发现，矩阵 \(C\) 的行数等于矩阵 \(A\) 的行数，列数等于矩阵 \(B\) 的列数。
因此矩阵乘法不满足交换律。
矩阵乘法运算律
- 乘法结合律，\(A \times (B \times C) = (A \times B) \times C\)。
- 乘法分配律，\(A \times (B + C) = A \times B + A \times C\)。
矩阵的逆
如果 \(A \times B = I\)，那么 \(B\) 是 \(A\) 的逆矩阵，记作 \(A^{-1}\)。
注意：逆矩阵不一定存在，但如果逆矩阵存在则是唯一的，至于怎么求，可以用高斯消元。
行列式
行列式。
实现
可以发现矩阵可以使用二维数组进行模拟。
显然地，线性运算时间复杂度 \(\mathcal{O}(n ^ 2)\)，乘法运算时间复杂度 \(\mathcal{O}(n^3)\)。
对于普通的矩阵这两个复杂度是不能优化的。
常数优化
当然进行矩阵乘法时，可以通过改变循环顺序加速高速缓存的访问优化部分常数（还挺猛），对于小矩阵也可以通过循环展开实现。
struct Matrix{
int mat[N][N];
Matrix(){
memset(mat, 0, sizeof(mat));
}
int * operator [] (const int & x){
return mat[x];
}
Matrix operator * (const Matrix & x){
Matrix ans;
for (int k = 1; k <= m; k++) {
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= x.m; j++) {
//multi
}
}
}
return ans;
}
};
矩阵乘法的运用
矩阵快速幂
过程
因为矩阵有结合律，所以自然可以使用快速幂进行优化矩乘。
实现
Matrix qpow(int b, Matrix & mat){//mat 是转移矩阵
Matrix ans;
//prework
for (; b; b >>= 1, mat = mat * mat) {
if (b & 1) {
ans = ans * mat;// matrix multi
}
}
return ans;
}
矩阵加速递推
一个经典的矩阵快速幂问题，斐波拉契数列。
在斐波拉契数列中，
很显然的可以构造矩阵，
实现
此处求斐波拉契第 \(n\) 项模 \(10 ^ 9 + 7\) 的结果的部分代码。
i64 n;
Matrix qpow(i64 b){
Matrix mat, ans;
ans[1][1] = 1, ans[1][2] = 0;
ans[2][1] = 0, ans[2][2] = 1;
mat[1][1] = 1, mat[1][2] = 1;
mat[2][1] = 1, mat[2][2] = 0;
for (; b; b >>= 1, mat = mat * mat) {
if (b & 1) {
ans = ans * mat;
}
}
return ans;
}
int main(){
scanf("%lld", &n);
Matrix ans = qpow(n - 1);
printf("%lld", ans[1][1]);
return 0;
}
线性常系数齐次递推式
当然，有时会在递推式中给定一些常数。
比如：
\(F_i = F_{i - 1} + F_{i - 2} + i ^ 2 + i\)。
\(i\) 的递推可以由 \(i - 1 + 1\) 得到，\(i ^ 2\) 可以由 \((i - 1) ^ 2 + i ^ 2 - (i - 1) ^ 2 = 2i - 1 + (i - 1)^2\) 得到。
所以维护 \(i^0,i^1, i^2\) 即可，就可以构造如下矩阵实现状态转移。
矩阵表修改
THUSCH2017 大魔法师
对于每个修改都可以通过矩阵表示，然后整体使用线段树维护即可。
比如操作一，可以构造一下矩阵实现：
操作四，可以构造一下矩阵实现：
参考代码
#include <cstdio>
#include <cstring>
#include <algorithm>
using i64 = long long ;
const int mod = 998244353 ;
const int N = 3e5 + 5 ;
struct Matrix{
i64 mat[4][4];
i64* operator [] (const int & _x) {return mat[_x];}
Matrix operator * (const Matrix & _x) {
Matrix ans;
memset(ans, 0, sizeof(ans));
for (int k = 0; k < 4; k++) {
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
ans[i][j] = (ans[i][j] + mat[i][k] * _x.mat[k][j]) % mod;
}
}
}
return ans;
}
Matrix operator + (const Matrix & _x) {
Matrix ans;
for(int i = 0; i < 4; i++) {
for(int j = 0; j < 4; j++) {
ans.mat[i][j] = (ans.mat[i][j] + mat[i][j]) % mod;
}
}
return _y;
}
};
Matrix magic[N];
#define lson(x) x << 1
#define rson(x) x << 1 | 1
struct Segment{
int l , r;
Matrix now;
Matrix lazy;
}s[N<<2];
struct SegmentTree{
void pushup(int k){
s[k].now = s[lson(k)].now + s[rson(k)].now;
}
void build(int k, int l, int r){
for(int i = 0; i < 4; i++) s[k].lazy[i][i] = 1;
s[k].l = l , s[k].r = r;
if(l == r) {
s[k].now = magic[l];
return;
}
int mid = l + r >> 1;
build(lson(k), l, mid);
build(rson(k), mid + 1, r);
pushup(k);
}
void pushdown(int k){
s[lson(k)].now = s[lson(k)].now * s[k].lazy;
s[rson(k)].now = s[rson(k)].now * s[k].lazy;
s[lson(k)].lazy = s[lson(k)].lazy * s[k].lazy;
s[rson(k)].lazy = s[rson(k)].lazy * s[k].lazy;
memset(s[k].lazy.mat, 0, sizeof(s[k].lazy.mat));
for(int i = 0; i < 4; i++)
s[k].lazy[i][i] = 1;
}
void modify(int k, int l, int r,const Matrix & tag){
if(l <= s[k].l && s[k].r <= r) {
s[k].now = s[k].now * tag;
s[k].lazy = s[k].lazy * tag;
return;
}
pushdown(k);
int mid = s[k].l + s[k].r >> 1;
if (l <= mid) modify(lson(k), l, r, tag);
if (r > mid) modify(rson(k), l, r, tag);
pushup(k);
}
Matrix query(int k, int l, int r){
if(l <= s[k].l && s[k].r <= r)
return s[k].now;
pushdown(k);
int mid = s[k].l + s[k].r >> 1;
Matrix val;
memset(val.mat , 0, sizeof(val.mat));
if (l <= mid) val = val + query(lson(k), l , r);
if (r > mid) val = val + query(rson(k), l , r);
return val;
}
}t;
#undef lson
#undef rson
int n,m;
int main(){
scanf("%d",&n);
for(int i = 1; i <= n; i++)
for(int j = 0; j < 3; j++)
scanf("%lld",&magic[i][0][j]) , magic[i][0][3] = 1;
t.build(1 , 1 , n);
scanf("%d",&m);
for(int i = 1; i <= m; i++) {
int op , l , r ;
i64 val;
scanf("%d %d %d",&op,&l,&r);
Matrix tag;
memset(tag.mat, 0, sizeof(tag.mat));
for(int j = 0; j < 4; j++) tag[j][j] = 1;
if(op == 1) {
tag[1][0] = 1;
t.modify(1,l,r,tag);
} else if(op == 2) {
tag[2][1] = 1;
t.modify(1,l,r,tag);
} else if(op == 3) {
tag[0][2] = 1;
t.modify(1,l,r,tag);
} else if(op == 4) {
scanf("%lld",&val);
tag[3][0] = val;
t.modify(1,l,r,tag);
} else if(op == 5) {
scanf("%lld",&val);
tag[1][1] = val;
t.modify(1,l,r,tag);
} else if(op == 6) {
scanf("%lld",&val);
tag[2][2] = 0;
tag[3][2] = val;
t.modify(1,l,r,tag);
} else {
Matrix ans = t.query(1,l,r);
for(int j = 0; j < 3; j++)
printf("%lld ",ans[0][j]);
puts("");
}
}
return 0;
}
矩阵统计定长路数量
TJOI2017 可乐
简要题意
有一个机器人，他初始在 \(1\) 号节点，它每秒可以去到与当前节点相邻的点或者原地不动或者自爆，问 \(t\) 秒钟，机器人一共有多少种行走方案。
节点总个数 \(n \land 1 \le n \le 30\)，边数 \(m \land 1 \le m \le 30\)
思路
可以构造出一个 dp 式子，\(f_{i,j,k}\) 表示时间 \(i\) 机器人从节点 \(j\) 走到节点 \(k\) 的方案数。
至于自爆，可以使机器人前往一个 \(0\) 号节点且 \(0\) 号节点出度为 \(0\) 入度为 \(n\)。
原地不动就直接连接自环即可。
因此，就有状态转移：
\(f_{i,j,k} = \sum f_{i - 1,j,p} + f_{i - 1, p, k}\)。
这个 dp 显然可以使用矩阵加速，实现 \(\mathcal{O}(n ^ 3 \log t)\) 解决问题。
最后答案 \(ans = \sum\limits_{i = 0} ^ n f_{t, 1, i}\)。
参考代码
#include <cstring>
#include <stdio.h>
using u32 = unsigned int ;
using i64 = long long ;
using u64 = unsigned long long ;
const int mod = 2017;
const int N = 35 ;
int n, m;
i64 t, ans;
struct Matrix {
i64 mat[N][N];
Matrix(){
memset(mat, 0, sizeof(mat));
}
i64 * operator [] (const int & x){return mat[x];}
Matrix operator * (const Matrix & x) {
Matrix ans;
for (int k = 0; k <= n; k++) {
for (int i = 0; i <= n; i++) {
for (int j = 0; j <= n; j++) {
ans[i][j] = (ans[i][j] + mat[i][k] * x.mat[k][j] % mod) % mod;
}
}
}
return ans;
}
}mat;
void qpow(i64 k){
Matrix Ans;
for (int i = 0; i <= n; i++) {
Ans[i][i] = 1;
}
for (; k; k >>= 1, mat = mat * mat) {
if (k & 1) {
Ans = Ans * mat;
}
}
for (int i = 0; i <= n; i++) {
ans = (ans + Ans[1][i]) % mod;
}
}
int main(){
scanf("%d %d", &n, &m);
for (int i = 1; i <= m; i++) {
int u, v;
scanf("%d %d", &u, &v);
mat[u][v] = mat[v][u] = 1;
}
scanf("%lld", &t);
for (int i = 0; i <= n; i++) {
mat[i][i] = mat[i][0] = 1;
}
qpow(t);
printf("%lld", ans);
return 0;
}
推广
此方法可以推广，即只需要有一个矩阵表示时间 \(0\) （或者行走了 \(1\) 条边时）时的方案数，然后通过矩阵乘法加速转移就可以完成实现行走 \(k\) 条边时，\(i\) 到 \(j\) 的方案数。
矩阵求定长最短路
问题
给定一个图，问从一号节点经过 \(k(1 \le k \le 10^9)\) 条边与 \(n(1 \le n \le 50)\) 节点的最短路。
思路
显然可以设状态 \(f_{k,i,j}\) 表示经过 \(k\) 条边从 \(i\) 节点到 \(j\) 节点的最短路。
显然这是 Floyd 算法，因此有状态装移 \(f_{k,i,j} = \min\{f_{k - 1, i, p} + f_{k - 1, p, j}\}\)。
和上面的一样可以构造矩阵，实现 \(\mathcal{O}(n ^ 3 \log k)\) 解决。
注意：
因为此题的转移是通过类似 Floyd 算法转移的，因此矩阵乘法需要改为 Floyd。
限长最短路/路径计数
路径计数
问题1
问一个图中，长度小于等于 \(k\) 的路径有多少条。
做法
可以通过添加一个虚点记录答案解决这个问题，并在图中 \((v, v^{'})\) 和 \((v^{'}, v^{'})\) 这两条边即可。
当然也可以通过构造矩阵时，用代数方式直接添加新的变量统计答案即可。
最后使用矩阵快速幂求解即可。
问题2
问一个图中，经过的边数小于等于 \(k\) 的最短路。
做法
每个点向自己连边并且此边边权为 \(0\) 即可。
最后也是用矩阵快速幂迭代 \(k\) 此即可。
