关于伴随矩阵/矩阵初等变换/动态矩阵求逆的一些想法

起因是某道 BEST 定理题需要求出矩阵修改一个位置后所有主对角线上的（代数）余子式。也就是求出删除 \(i\) 行 \(i\) 列后的伴随矩阵。而伴随矩阵可以用逆矩阵计算，问题又变成了删除一行一列再加入一行一列的矩阵求逆。（因为 Laplace 矩阵不满秩而所有的主子式满秩，所以你不能求出原矩阵的逆矩阵，只能先求出原矩阵删除 \(i\) 行 \(i\) 列的逆矩阵再删再加）。

环一周场切了此题，他的方法利用了 Laplace 矩阵的性质似乎更加好。我自己研究了很久终于研究出一种做法，很有可能做捞了不过也记录一下。

关于伴随矩阵

设 \(n\) 阶方阵 \(A\) 每个位置的代数余子式构成的矩阵为 \(M\)。

行列式的 Laplace 展开讲的是对于一个行列式我们可以将其值表示为第 \(i\) 行所有代数余子式的线性组合，即： \(|A|=\sum_{j=1}^n A_{i,j}M_{i,j}\)

特别地，将 \(A\) 中的一行复制到令一行上去，考察这个有两行完全相等的行列式，其值一定是 0，我们把它按照其中一行展开，得到了这样的结果：\(\forall k\neq i,\sum_{j=1} A_{i,j}M_{k,j}=0\)。

这是一个类似点乘的形式，我们转置 \(M\) 矩阵，设 \(A^*=M^T\)，由上面的结果就有了 \(AA^*=|A|I_n\)，我们也得到了计算非奇异矩阵所有位置的代数余子式的方法：\(A^*=|A|A^{-1}\)。

这个结果有什么用呢？考虑到你单点修改一个行列式的取值，由 Laplace 展开的结果，行列式值的变化量正好就是这个位置变化量的 \(M_{i,j}\) 倍，所以求出了伴随矩阵，我们就能够刻画单点修改行列式的结果！

关于矩阵的初等变换

初等行/列变换也许是我们 OIer 最先接触的线代内容了。早在学习高斯消元的时候我们就有所掌握。后续我们学习矩阵求逆/行列式求值的时候也在使用基于初等变换的方式。但是我们真的理解透了它的线代意义吗？为什么它可以胜任各种矩阵操作呢？

首先一个误区：初等行列变换真的都是线性变换吗？答案其实是否定的！

比如说给一列加上另一列的 \(k\) 倍，我们可以用类似这样的矩阵 \(\begin{pmatrix}1&0 &k\\ 0&1 &0\\0&0 &1\\\end{pmatrix}\) 描述。但是给一行加上另一行的 \(k\) 倍呢？

答案：左乘刚刚那个矩阵！由于矩阵的左右乘不尽相同，当你顺序进行一系列的初等变换后，描述成矩阵乘法角度相当于左乘一个矩阵再右乘一个矩阵。

其实所有基于初等行/列变换的算法关键在于其可以执行这样一个操作：给定一个矩阵 \(A\)，将其分解成若干个初等行/列变换矩阵的乘积。而左/右乘一个初等行/列变换矩阵的代价是很小的（可以做到 \(O(n)\)），这样，就为我们绝大多数的线性代数算法提供了基石。

关于动态求逆矩阵

现在回到原来的问题上来，我们知道了初等变换矩阵拥有左乘/右乘的区别，绝大多数线性代数算法只会把矩阵分解成单纯左乘或是单纯右乘，但是执行动态求逆操作时我们发现左右乘的区别至关重要！（原来我没有深入理解初等变换导致我直接按照初等变换的顺序一股脑把矩阵乘上去调了 INF 年）

考察经过一次初等变换之后逆矩阵会如何变化，假设该初等变换相当于右乘一个矩阵 \(T\)，那么 \((AT)^{-1}=T^{-1}A^{-1}\)，即：

• 交换原矩阵的 \(i,j\) 行 \(\iff\) 交换逆矩阵的 \(i,j\) 列

• 原矩阵的 \(i\) 行乘 \(k\) \(\iff\) 逆矩阵的 \(i\) 列乘 \(k^{-1}\)

• 原矩阵的 \(i\) 行加上 \(j\) 行的 \(k\) 倍 \(\iff\) 逆矩阵的 \(i\) 列减去 \(j\) 列的 \(k\) 倍

首先考虑如何删除一行一列求逆矩阵：

先假设删除的是 \(i\) 行 \(i\) 列，不是可以通过逆矩阵初等变换操作调整成这样。

注意到删除第 \(i\) 行第 \(i\) 列相当于将原矩阵 \(A\) 的第 \(i\) 行第 \(i\) 列置成 0，然后再在令 \(A_{i,i}=1\)，这样求完逆之后，第 \(i\) 行第 \(i\) 列仍然是只有 \((i,i)\) 这个位置有值。把逆矩阵的第 \(i\) 行第 \(i\) 列删掉就行了。

现在的问题是对原矩阵进行若干次初等变换，使得第 \(i\) 行第 \(j\) 列被消成只有一个 \(1\)。

这看起来是困难的，因为这个东西不弱于解线性方程，但是一个好事是我们已经拥有了原矩阵的逆矩阵，也就是我们已经解完了线性方程！

考虑逆矩阵的向量意义：逆矩阵相当于找出了一个原矩阵所有行向量的 \(n\) 个线性组合得出所有只有一个位置有值的单位向量 \((0,\dots,0,1,0,\dots,0)\)。那么我们就知道了一种进行初等行变换的方式使得第 \(i\) 行被清成只有第 \(i\) 位有 1：即按照逆矩阵第 \(i\) 行将所有行向量线性组合到第 \(i\) 行。然后再通过若干次线性变换消掉第 \(i\) 列上的所有值即可。这个过程中只有行变换，也就是只有左乘。

复杂度做到了下界 \(O(n^2)\)。

接下来是如何加入一行一列求逆矩阵：

同样的套路，我们只需要考虑加入的是第 \(n\) 行第 \(n\) 列的情况。我们考虑先给原矩阵填充一个空行空列，\(A_{n,n}\) 填充 1。然后考虑通过若干次初等变换把它变成新矩阵 \(B\)。

这些初等变换咋找呢？我的方法显得有点憨，也许不是最简思路：发现 \(A\to B\) 需要乘上一个矩阵 \(A^{-1}B\)，由于 \(A,B\) 差别不是很大，计算发现 \(A^{-1}B\) 的组成实际上是由一个大小为 \((n-1)\times (n-1)\) 的单位矩阵拼上某一列再拼上新矩阵 \(B\) 的最后一行。

拼上的那一列可以 \(O(n^2)\) 简单计算。

现在我们考虑乘上 \(A^{-1}B\) 对逆矩阵的影响。我们考虑将其分解成初等变换矩阵。怎么做呢？我们考虑通过 \(O(n)\) 次初等变换把消成单位阵，这是容易做的，这些变换中有的是行变换有的是列变换，注意一下最后把它写成初等变换矩阵的顺序就行了。

复杂度也是 \(O(n^2)\)。

代码是那道 BEST 题的代码，其中 \(a\) 是原矩阵，\(f\) 是原逆矩阵，\(g\) 是计算结果。代码实现了删除第 \(i\) 行第 \(i\) 列然后加入第 \(n\) 行第 \(n\) 列。

#include <cstdio>
#include <cassert>
#include <algorithm>
using namespace std;
typedef long long ll;
const int P=998244353;
const int N=303;
int read(){
	char c=getchar();int x=0;
	while(c<48||c>57) c=getchar();
	do x=x*10+(c^48),c=getchar();
	while(c>=48&&c<=57);
	return x;
}
int qp(int a,int b=P-2){
	int res=1;
	while(b){
		if(b&1) res=(ll)res*a%P;
		a=(ll)a*a%P;b>>=1;
	}
	return res;
}
int n,m;
int fac[N],inv[N],tmp[N];
int adj[N][N],a[N][N];
int w[N][N],f[N][N],g[N][N];
int inver(int x){
	for(int i=1;i<=n;++i) if(i!=x)
		for(int j=1;j<=n;++j) if(j!=x)
			w[i-(i>x)][j-(j>x)]=a[i][j];
	for(int i=1;i<n;++i)
		for(int j=1;j<n;++j) f[i][j]=(i==j);
	int res=1;
	for(int i=1;i<n;++i){
		int p=i;
		while(p<n&&!w[p][i]) ++p;
		if(p==n) return 0;
		if(p!=i){
			for(int j=i;j<n;++j) swap(w[p][j],w[i][j]);
			for(int j=1;j<n;++j) swap(f[p][j],f[i][j]);
			res=P-res;
		}
		res=(ll)res*w[i][i]%P;
		int iv=qp(w[i][i]);
		for(int j=1;j<n;++j){
			if(j>=i) w[i][j]=(ll)w[i][j]*iv%P;
			f[i][j]=(ll)f[i][j]*iv%P;
		}
		for(int j=i+1;j<n;++j){
			int val=w[j][i];
			for(int k=1;k<n;++k){
				if(j>=i) w[j][k]=(w[j][k]+(ll)(P-val)*w[i][k])%P;
				f[j][k]=(f[j][k]+(ll)(P-val)*f[i][k])%P;
			}
		}
	}
	for(int i=n-1;i;--i)
		for(int j=i-1;j;--j){
			for(int k=1;k<=n;++k) f[j][k]=(f[j][k]+(ll)(P-w[j][i])*f[i][k])%P;
			w[j][i]=0;
		}
	return res;
}
int main(){
	n=read();m=read();fac[0]=inv[1]=1;
	for(int i=1;i<=n;++i) fac[i]=(ll)fac[i-1]*i%P;
	for(int i=2;i<=n;++i) inv[i]=(ll)inv[P%i]*(P-P/i)%P;
	for(int i=1;i<=m;++i){
		int c=read(),u=read(),v=read();
		adj[u][v]=c;++a[u][u];--a[u][v];
	}
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n;++j) if(a[i][j]<0) a[i][j]+=P;
	int cur=inver(n);
	for(int i=1;i<=n;++i) cur=(ll)cur*fac[a[i][i]-1]%P;
	ll res=0;
	for(int i=1;i<=n;++i)
		if(a[i][i]>1){
			if(i<n){
				for(int j=1;j<n;++j)
					for(int k=1;k<n;++k) g[j][k]=f[j][k];
				int iv=qp(f[i][i]);
				for(int j=1;j<n;++j) g[j][i]=(ll)g[j][i]*iv%P;
				for(int j=1;j<n;++j) if(i!=j){
					for(int k=1;k<n;++k)
						g[k][j]=(g[k][j]+(ll)(P-f[i][j])*g[k][i])%P;
				}
				for(int j=1;j<n;++j) if(i!=j){
					for(int k=1;k<n;++k)
						g[k][i]=(g[k][i]+(ll)a[j][i]*g[k][j])%P;
				}
				for(int j=1;j<n-1;++j)
					for(int k=1;k<n-1;++k) g[j][k]=g[j+(j>=i)][k+(k>=i)];
				for(int j=1;j<n-1;++j) g[j][n-1]=g[n-1][j]=0;
				g[n-1][n-1]=1;
				for(int j=1;j<n-1;++j){
					tmp[j]=0;
					for(int k=1;k<n-1;++k)
						tmp[j]=(tmp[j]+(ll)g[j][k]*a[k+(k>=i)][n])%P;
				}
				int cur=a[n][n];
				for(int j=1;j<n-1;++j) cur=(cur+(ll)a[n][j+(j>=i)]*(P-tmp[j]))%P;
				iv=qp(cur);
				for(int j=1;j<n-1;++j){
					int val=a[n][j+(j>=i)];
					for(int k=1;k<n;++k) g[n-1][k]=(g[n-1][k]+(ll)(P-val)*g[j][k])%P;
				}
				for(int j=1;j<n-1;++j)
					for(int k=1;k<n;++k) g[j][k]=(g[j][k]+(ll)iv*(P-tmp[j])%P*g[n-1][k])%P;
				for(int j=1;j<n;++j) g[n-1][j]=(ll)g[n-1][j]*iv%P;
			}
			else{
				for(int j=1;j<n;++j)
					for(int k=1;k<n;++k) g[j][k]=f[j][k];
			}
			for(int j=1;j<=n;++j) if(adj[j][i]){
				int c=adj[j][i];
				for(int k=1;k<=n;++k) if(adj[i][k]==c)
					res+=(ll)cur*(P+1-g[k-(k>i)][j-(j>i)])%P*inv[a[i][i]-1]%P;
			}
		}
		else{
			for(int j=1;j<=n;++j) if(adj[j][i]){
				int c=adj[j][i];
				for(int k=1;k<=n;++k) if(adj[i][k]==c) res+=cur;
			}
		}
	printf("%d\n",int(res%P));
	return 0;
}