浅谈矩阵在OI中的应用

本文作者太菜了，不会用专业的数学术语表达，在本文中用相对通俗的语言进行讲解。

定义

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

可以想象为一个 $n\times m$ 的二维数组，可以表示为：

A = [\begin{matrix} a_{1, 1} & a_{1, 2} & \dots & a_{1, m} \\ a_{2, 1} & ⋱ & a_{2, m} \\ ⋮ & ⋱ \\ a_{n, 1} & a_{n, m - 1} & a_{n, m} \end{matrix}]

$A= \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,m} \\ a_{2,1} & \ddots & & a_{2,m} \\ \vdots & & \ddots & \\ a_{n,1} & & a_{n,m-1} & a_{n,m} \\ \end{bmatrix}$

如果 $n=m$ ，那么 $A$ 也被称为 $n$ 阶矩阵。

单位矩阵 $I$ 表示 $I_{i,i}=1$ ，其余都为 $0$ 的矩阵。

基本运算

加减法表示两个大小相同的矩阵每个对于位置相加。

即

A + B = [\begin{matrix} a_{1, 1} + b_{1, 1} & a_{1, 2} + b_{1, 2} & \dots & a_{1, m} + b_{1, m} \\ a_{2, 1} + b_{2, 1} & ⋱ & a_{2, m} + b_{2, m} \\ ⋮ & ⋱ \\ a_{n, 1} + b_{n, 1} & a_{n, m - 1} + b_{n, m - 1} & a_{n, m} + b_{n, m} \end{matrix}]

$A+B= \begin{bmatrix} a_{1,1}+b_{1,1} & a_{1,2}+b_{1,2} & \cdots & a_{1,m}+b_{1,m} \\ a_{2,1}+b_{2,1} & \ddots & & a_{2,m}+b_{2,m} \\ \vdots & & \ddots & \\ a_{n,1}+b_{n,1} & & a_{n,m-1}+b_{n,m-1} & a_{n,m}+b_{n,m} \\ \end{bmatrix}$

一个矩阵乘上一个数，表示每个位置都乘上这个数。

即

A \times b = [\begin{matrix} a_{1, 1} \times b & a_{1, 2} \times b & \dots & a_{1, m} \times b \\ a_{2, 1} \times b & ⋱ & a_{2, m} \times b \\ ⋮ & ⋱ \\ a_{n, 1} \times b & a_{n, m - 1} \times b & a_{n, m} \times b \end{matrix}]

$A\times b= \begin{bmatrix} a_{1,1}\times b & a_{1,2}\times b & \cdots & a_{1,m}\times b \\ a_{2,1}\times b & \ddots & & a_{2,m}\times b \\ \vdots & & \ddots & \\ a_{n,1}\times b & & a_{n,m-1}\times b & a_{n,m}\times b \\ \end{bmatrix}$

乘法

两个矩阵相乘仅当第一个矩阵 $A$ 的列数和另一个矩阵 $B$ 的行数相等时才能乘。如 $A$ 是 $m\times n$ 矩阵和 $B$ 是 $n\times p$ 矩阵，它们的乘积 $C$ 是一个 $m\times p$ 矩阵。

乘法满足

C_{i, j} = \sum_{k = 1}^{n} A_{i, k} \times B_{k, j}

$C_{i,j}=\sum_{k=1}^n A_{i,k}\times B_{k,j}$

计该运算为 $C=AB$ 。

也就是 $C_{i,j}$ 表示 $A$ 的第 $i$ 行， $B$ 的第 $j$ 列的每一个元素乘起来的和。

Matr operator *(const Matr &x,const Matr &y){
	Matr z;
	memset(z.cnt,0,sizeof(z.cnt));
	for(int k=1;k<=n;++k)
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				z.cnt[i][j]=(z.cnt[i][j]+x.cnt[i][k]*1ll*y.cnt[k][j]%mod)%mod;
	return z;
}

需要注意的是，矩阵乘法满足结合律和分配率，但不满足交换率。

因此，矩阵乘法是可以用快速幂的。

需要注意的是，对于任意矩阵 $A$ ， $A^0=I$ ，所以可以将 $I$ 定义为初值。

【模板】矩阵快速幂

板子题，直接套一个矩阵乘法和快速幂就好了。

时间复杂度 $O(n^3\log k)$

#include<bits/stdc++.h>
using namespace std;
const long long N=105;
long long n,k,mod=1E9+7;
struct hzy{long long cnt[N][N];}a;
hzy operator *(const hzy &x,const hzy &y){
	hzy z;
	memset(z.cnt,0,sizeof(z.cnt));
	for(long long k=1;k<=n;++k)
		for(long long i=1;i<=n;++i)
			for(long long j=1;j<=n;++j)
				z.cnt[i][j]=(z.cnt[i][j]+x.cnt[i][k]*y.cnt[k][j]%mod)%mod;
	return z;
}
hzy kuai(hzy a,long long b){
	hzy c;
	for(long long i=1;i<=n;i++)c.cnt[i][i]=1;
	for(;b;b>>=1){
		if(b&1)c=c*a;
		a=a*a;
	}
	return c;
}
int main(){
	scanf("%lld%lld",&n,&k);
	for(long long i=1;i<=n;i++)
		for(long long j=1;j<=n;j++)
			scanf("%lld",&a.cnt[i][j]);
	a=kuai(a,k);
	for(long long i=1;i<=n;i++){
		for(long long j=1;j<=n;j++)printf("%lld ",a.cnt[i][j]);
		puts("");
	}
}//之前写的代码，写的实在太丑了

矩阵快速幂有啥用？

矩阵快速幂的一个重要应用就是优化 dp 转移。

来个例子斐波那契数列。

众所周知，斐波那契数列是可以线性递推的，时间复杂度 $O(n)$ 。

尝试使用矩阵乘法进行优化。

我们发现， $fib_i$ 只与 $fib_{i-1}$ 和 $f_{i-2}$ 有关。

那么我们构造一个矩阵 $\begin{bmatrix} fib_{i-2} & fib_{i-1} \\ \end{bmatrix}$ 。

我们通过乘法将它变成 $\begin{bmatrix} fib_{i-1} & fib_{i-2}+fib_{i-1} \\ \end{bmatrix}$

容易得到它乘上的矩阵 $A=\begin{bmatrix} 0 & 1 \\ 1 & 1 \\ \end{bmatrix}$

事实上，我们发现 $A_{i,j}$ 就表示在上个状态中的 $f_{t-1,i}$ 给 $f_{t,j}$ 的贡献就是 $A_{i,j}\times f_{t-1,i}$ ，因此很容易推出。

乘一次，就向后转移一次，根据结合律对后面相同的矩阵进行合并。

容易发现我们要运算的实际上是 $\begin{bmatrix} 1 & 1 \\ \end{bmatrix}\times A^{n-1}$ ，直接矩阵快速幂即可。

时间复杂度 $O(\log n)$ 。

mat a,b;
int main(){
	scanf("%lld",&n);
	a.cnt[1][1]=0;a.cnt[1][2]=1;
	a.cnt[2][1]=1;a.cnt[2][2]=1;
	a=kuai(a,n-1);
	b.cnt[1][1]=1;
	b.cnt[1][2]=1;
	b=b*a;
	printf("%lld",b.cnt[1][1]%mod);
}

关于斐波那契数列的计算，还可以用生成函数优化为 $O(1)$ 。

练习：【模板】矩阵加速。

来个例题[TJOI2017]可乐（数据加强版）

还是类似，考虑原地不动就是一个自环，自爆就见一个虚点，所有的点都连向它。

设 $f_{i,j}$ 表示 $i$ 时刻最后停留在 $j$ 的方案数。

暴力 dp 过不了，考虑每次转移的本质都相同，直接矩阵快速幂。

转移矩阵实际上就是邻接矩阵（根据上面推出转移矩阵的结论）。

然后直接做就好了。

时间复杂度 $O(n^3\log t)$ 。

int main(){
	Matr c;
	memset(c.cnt,0,sizeof(c.cnt));
	int ans=0;
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++){
		int x,y;
		scanf("%d%d",&x,&y);
		c.cnt[x][y]=c.cnt[y][x]=1;
	}
	scanf("%d",&k);
	n++;
	for(int i=1;i<=n;i++)c.cnt[i][i]=1;
	for(int i=1;i<n;i++)c.cnt[i][n]=1;
	c=kuai(c,k);
	for(int i=1;i<=n;i++)(ans+=c.cnt[1][i])%=mod;
	printf("%d",ans);
}

仍然考虑上面的这个问题，如果有 $Q$ 次询问，每次询问时刻 $t$ 时的答案怎么做， $Q\le 100$ 。

如果暴力每次做，时间复杂度 $O(Qn^3\log t)$ ，炸了。

容易发现做一次矩阵乘法的时间复杂度是 $O(n^3)$ 的，但我们发现答案的矩阵是 $1\times n$ 的。

答案的这个矩阵乘上转移矩阵是 $O(n^2)$ 的。

设初始矩阵为 $P$ 。

答案就是 $P\times A^t$ 。

快速幂拆成 $P\times A^{2^{t_0}}\times A^{2^{t_1}}\times A^{2^{t_ 2}} ...\times A^{2^{t_k}}$

根据结合律，变成 $(...(((P\times A^{2^{t_0}})\times A^{2^{t_1}})\times A^{2^{t_ 2}}) ...\times A^{2^{t_k}})$

每次运算都是 $O(n^2)$ 的，开始时先预处理后面要乘的 $A^{2^t}$ 。

时间复杂度 $O(Qn^2\log t+n^3\log t)$ 。

练习：

【XR-1】分块

Hanjo 2

都是推出暴力式子后直接矩阵快速幂的。

定长最短路

考虑矩阵乘法：

Matr operator *(const Matr &x,const Matr &y){
	Matr z;
	memset(z.cnt,0,sizeof(z.cnt));
	for(int k=1;k<=n;++k)
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				z.cnt[i][j]=(z.cnt[i][j]+x.cnt[i][k]*1ll*y.cnt[k][j]%mod)%mod;
	return z;
}

你会发现这个东西和 Floyd 不仅形似而且神似。

改亿改，变成 Floyd 的形式：

Matr operator *(const Matr &x,const Matr &y){
	Matr z;
	memset(z.cnt,0,sizeof(z.cnt));
	for(int k=1;k<=n;++k)
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				z.cnt[i][j]=min(z.cnt[i][j],max(x.cnt[i][k]+y.cnt[k][j]));
	return z;
}

现在变成了我们熟悉的形式，但有什么用呢？

考虑设邻接矩阵为 $A$ ，那么根据前面 dp 的思路， $A^L$ 表示走 $L$ 步后得到的最短路的邻接矩阵。

来个模板[USACO07NOV]Cow Relays G。

和上面说的一模一样，需要注意的是邻接矩阵中 $A_{i,i}=Inf$ ，不然就可以一直待在原地。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1005;
struct hzy{int cnt[N][N];}a;
int k,m,st,en,b[N],x[N],y[N],z[N],id[N],tot,dis[N][N];
hzy operator*(const hzy x,const hzy y){
    hzy z;
    memset(z.cnt,63,sizeof(z.cnt));
    for(int k=1;k<=tot;k++)
        for(int i=1;i<=tot;i++)
            for(int j=1;j<=tot;j++)z.cnt[i][j]=min(z.cnt[i][j],x.cnt[i][k]+y.cnt[k][j]);
    return z;
}
hzy kuai(hzy a,int b){
    hzy ans;
    memset(ans.cnt,63,sizeof(ans.cnt));
    for(int i=1;i<=tot;i++)ans.cnt[i][i]=0;
    for(;b;b>>=1,a=a*a)if(b&1)ans=ans*a;
    return ans;
}
int main() {
    memset(dis,63,sizeof(dis));
    scanf("%d%d%d%d",&k,&m,&st,&en);
    for(int i=1;i<=m;i++){
        scanf("%d%d%d",&z[i],&x[i],&y[i]);
        b[x[i]]=1,b[y[i]]=1;
    }
    for(int i=1;i<=1000;i++)if(b[i])id[i]=++tot;
    for(int i=1;i<=m;i++)x[i]=id[x[i]],y[i]=id[y[i]];
    for(int i=1;i<=m;i++)
        dis[y[i]][x[i]]=dis[x[i]][y[i]]=min(dis[x[i]][y[i]],z[i]);
    for(int i=1;i<=tot;i++)for(int j=1;j<=tot;j++)a.cnt[i][j]=dis[i][j];
    printf("%d",kuai(a,k).cnt[id[en]][id[st]]);
    return 0;
}

如果修改题目条件，变成走至多 $L$ 步的最短路。

直接每个点连一个边权为 $0$ 的自环，那么就相当于可以一直留在原地，就做完了。

再来个例题Good Vertices。

和前面几乎一样，还是跑 Flody，只不过修改一下就好了, $A_{i,j}$ 表示 $i$ 走 $L$ 步到 $j$ 的最小时刻。

mat operator*(mat a,mat b){
    mat ans;
    for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)ans.cnt[i][j]=1e9;
    for(int k=1;k<=n;k++)
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)ans.cnt[i][j]=min(ans.cnt[i][j],max(a.cnt[i][k],b.cnt[k][j]));
    return ans;
}

还有例题[NOI Online #1 入门组] 魔法。

$n$ 很小，依然是 Flody，然后考虑一个转移矩阵表示用一次魔法后的邻接矩阵。

考虑使用一次魔法的改变即可。

转移矩阵：

	for(int k=1;k<=n;k++)
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++)
				dis[i][j]=min(dis[i][j],dis[i][k]+dis[k][j]);
	for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)p.cnt[i][j]=dis[i][j];
	for(int i=1;i<=m;i++)//考虑每一条边，把他变成负数
		for(int l=1;l<=n;l++)
			for(int r=1;r<=n;r++)
				p.cnt[l][r]=min(p.cnt[l][r],dis[l][x[i]]+dis[y[i]][r]-z[i]);

练习：[NOI2020] 美食家。

维护数列

矩阵乘法还可以和数据结构一起食用。

来个题Fibonacci-ish II。

考虑没有修改操作，那么直接用莫队，这样就只需要考虑加一个或者删掉一个数就好了。

考虑新加入一个数会让后面的数都向后平移，即 $a_x\times fib_y$ 变成了 $a_x\times fib_{y+1}$ 。

那么这个变化就相当于乘上了一个 $fib$ 的矩阵。

那么这个东西就可以直接用线段树维护了。

code。

没想到吧，人傻常数大，这个代码 Time limit exceeded on test 26 了。

考虑优化。

其实不需要延迟标记，直接考虑单点修改，只要可以维护两个区间直接的合并就好了。

原来两个区间的和分别是

x_{1} f i b_{1} + x_{2} f i b_{2} . . . + x_{k} f i b_{k}

$x_1fib_1+x_2fib_2...+x_kfib_k$

x_{k + 1} f i b_{1} + x_{k + 2} f i b_{2} . . . + x_{k + n} f i b_{n}

$x_{k+1}fib_1+x_{k+2}fib_2...+x_{k+n}fib_n$

合并成

x_{1} f i b_{1} + x_{2} f i b_{2} . . . + x_{k + n} f i b_{k + n}

$x_{1}fib_1+x_2fib_2...+x_{k+n}fib_{k+n}$

实际上就是右边的乘上 $Fib^k$ ，因为 $n$ 很小，直接预处理就好了。

code。

练习[THUSCH2017] 大魔法师。

动态dp

~~动态动态规划。~~

~~还在学，快了。~~

To be continue。

行列式

定义一个行列式运算 $|A|$ ，也可以表示为 $\det(A)$ ， $A$ 是一个 $n$ 阶矩阵。

det (A) = \sum_{P \in S} (- 1)^{f (P)} \prod_{i = 1}^{n} A_{i, P_{i}}

$\det(A)=\sum _{P\in S} (-1)^{f(P)} \prod _{i=1}^n A_{i,P_i}$

$P$ 表示一个排列， $S$ 表示全排列集合， $f(P)$ 表示 $P$ 的逆序对数量。

也就是每行每列选一个数乘起来，再乘一个容斥系数。

暴力求的时间复杂度为 $O(n!\times n)$

来一些性质：

如果 $A$ 是三角矩阵，那么 $\det(A)=\prod_{i=1}^n A_{i,i}$ ，这个比较显然。
交换 $A$ 的两行或者两列，那么 $\det(A)$ 取相反数。

证明：考虑两行，取得位置为 $P_{i},P_j$ 。

交换之后实际上就是逆序对中交换了两个数。

讨论一下就可以发现逆序对个数的变化一定是个奇数。

所以行列式的值取相反数。

交换两列也差不多。

如果 $A$ 有两行一模一样，那么 $|A|=0$ 。

证明：交换这两行， $|A|=-|A|$ 。

基本变换，将矩阵一行上所有的数乘 $t$ ，加到另一行上， $|A|$ 不变。

这个也比较简单，加过去，通过乘法分配律拆成 $t+1$ 个矩阵，然后发现后面的矩阵都有两行相等，都是 $0$ ，所以 $|A|$ 不变。

发现没有，这个变换和高斯消元一模一样。

可以用高斯消元消为一个上三角，把对角线乘起来就好了。

时间复杂度 $O(n^3)$ 。

#define For(x,l,r) for(int x=l;x<=r;x++)
int det(mat a){
    int ans=1,ff=1;
    For(i,1,n){
        int k=i;
        while(!a.cnt[k][i]&&k<=n)k++;
        if(k!=i){ff*=-1;For(j,i,n)swap(a.cnt[i][j],a.cnt[k][j]);}//交换时记得取相反数
        if(!a.cnt[i][i])return 0;
        ans=1ll*ans*a.cnt[i][i]%mod;
        int Inv=kuai(a.cnt[i][i],mod-2);
        For(j,i+1,n){//普通高斯消元
            if(!a.cnt[j][i])continue;
            int t=1ll*Inv*a.cnt[j][i]%mod;
            For(p,i,n)a.cnt[j][p]=(a.cnt[j][p]-1ll*t*a.cnt[i][p]%mod+mod)%mod;
        }
    }
    return (ans*ff+mod)%mod;
}

然后发现这并不能过【模板】行列式求值。

当 $p$ 非质数时，原来的做法会有问题。

考虑中间那个消元过程，就是一个不断减掉的过程，可以辗转相减。

#define For(i,l,r) for(int i=l;i<=r;i++)
int det(mat a){
    int ans=1,ff=1;
    For(i,1,n)For(j,i+1,n){
        while(a.cnt[i][i]){
            int div=a.cnt[j][i]/a.cnt[i][i];
            For(k,i,n)a.cnt[j][k]=(a.cnt[j][k]-1ll*div*a.cnt[i][k]%mod+mod)%mod;
            ff=-ff,swap(a.cnt[i],a.cnt[j]);
        }
        ff=-ff,swap(a.cnt[i],a.cnt[j]);
    }
    for(int i=1;i<=n;i++)ans=1ll*ans*a.cnt[i][i]%mod;
    return  (ans*ff+mod)%mod;
}

时间复杂度仍然是 $O(n^3)$ 。

范德蒙德行列式

D_{n} = [\begin{matrix} 1 & 1 & \dots & 1 \\ x_{1} & ⋱ & x_{n} \\ ⋮ & ⋱ \\ x_{1}^{n - 1} & x_{n - 1}^{n - 1} & x_{n}^{n - 1} \end{matrix}] = \prod_{1 \leq i < j \leq n} (x_{j} - x_{i})

$D_n=\begin{bmatrix} 1 & 1 & \cdots & 1 \\ x_1 & \ddots & & x_n \\ \vdots & & \ddots & \\ x_1^{n-1} & & x_{n-1}^{n-1} & x_n^{n-1} \\ \end{bmatrix}=\prod_{1\le i<j\le n}(x_j-x_i)$

考虑证明，使用归纳法，对于 $n=2,3$ ，可以直接暴力展开证明。

如果 $n-1$ 成立，考虑对于 $n$ 阶矩阵，把第 $i-1$ 行每个数乘上 $x_1$ ，然后让第 $i$ 行去减掉，倒着做一遍就可以得到：

D_{n} = [\begin{matrix} 1 & 1 & \dots & 1 \\ 0 & ⋱ & x_{n} - x_{1} \\ ⋮ & ⋱ \\ 0 & (x_{n - 1} - x_{1}) x_{n - 1}^{n - 2} & (x_{n} - x_{1}) x_{n}^{n - 2} \end{matrix}] = \prod_{i = 1}^{n - 1} (x_{n} - x_{i}) D_{n - 1} = \prod_{1 \leq i < j \leq n} (x_{j} - x_{i})

$D_n=\begin{bmatrix} 1 & 1 & \cdots & 1 \\ 0 & \ddots & & x_n-x_1 \\ \vdots & & \ddots & \\ 0 & & (x_{n-1}-x_1)x_{n-1}^{n-2} & (x_n-x_1)x_n^{n-2} \\ \end{bmatrix}=\prod_{i=1}^{n-1}(x_n-x_i)D_{n-1}=\prod_{1\le i<j\le n}(x_j-x_i)$

LGV引理

先来个模板：[NOI2021] 路径交点

这个引理主要解决有向无环图中不相交路径求和问题。

已知一张有向无环图。

M = [\begin{matrix} f_{1, 1} & f_{1, 2} & \dots & f_{1, m} \\ f_{2, 1} & ⋱ & f_{2, m} \\ ⋮ & ⋱ \\ f_{n, 1} & f_{n, m - 1} & f_{n, m} \end{matrix}]

$M= \begin{bmatrix} f_{1,1} & f_{1,2} & \cdots & f_{1,m} \\ f_{2,1} & \ddots & & f_{2,m} \\ \vdots & & \ddots & \\ f_{n,1} & & f_{n,m-1} & f_{n,m} \\ \end{bmatrix}$

$f_{i,j}$ 表示 $i$ 走到 $j$ 的方案数。

M =

$M=$

杨氏矩阵

~~还在学，快了。~~

To be continue。

posted @ 2022-03-21 21:22 houzhiyuan 阅读(460) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 9月杂题选做

· 3月杂题选做（part 1）

· 矩阵学习笔记

· 【6】矩阵学习笔记

· 一些些数学的算法笔记

公告

昵称： houzhiyuan
园龄： 2年11个月
粉丝： 15
关注： 6

+加关注

2025年3月

日

一

二

三

四

五

六

houzhiyuan

浅谈矩阵在OI中的应用

定义

基本运算

乘法

定长最短路

维护数列

动态dp

行列式

LGV引理

杨氏矩阵

公告

搜索

常用链接

我的标签

随笔档案

阅读排行榜

评论排行榜

推荐排行榜

最新评论