正确的求解代数余子式的方法

前言

在联合省选day2t3中，存在一种使用行列式求导来计算生成树边权和的方法。需要计算出每个位置的代数余子式
常见的做法均是套用$A\times A^*=|A|\times I$A×A∗=∣A∣×I的等式求解，这种做法在模意义下矩阵$rank(A)=n-1$rank(A)=n−1时不能正确得出伴随矩阵
以下给出正确的求解伴随矩阵的方法

代数余子式

给$n$n阶方阵$A=(a_{i,j})$A=(ai,j​)，定义$a_{i,j}$ai,j​的余子式$M_{i,j}$Mi,j​为将$i$i行$j$j列划去后的行列式，$a_{i,j}$ai,j​的代数余子式$A_{i,j}$Ai,j​为$(-1)^{i+j}M_{i,j}$(−1)i+jMi,j​

伴随矩阵

若$A_{i,j}$Ai,j​为代数余子式矩阵，我们定义
$A^*=(A_{ij})^{\mathrm T}$A∗=(Aij​)T
$A^*$A∗为矩阵$A$A的伴随矩阵

我们存在一个结论，即
$A\times A^*=A^*\times A=|A|\times I_n$A×A∗=A∗×A=∣A∣×In​
其中$I_n$In​为$n$n阶单位矩阵，$|A|$∣A∣为矩阵$A$A的行列式

计算

若$rank(A)=n$rank(A)=n，则$A^{-1}$A−1是良定义的，可以直接用$A^*=|A|\times A^{-1}$A∗=∣A∣×A−1来进行计算。

若$rank(A)\leq n-2$rank(A)≤n−2，不难发现任意$n-1$n−1阶子式行列式均为$0$0，故$A^*=O$A∗=O，其中$O$O为空矩阵

我们只需要特殊考虑$rank(A)=n-1$rank(A)=n−1的情况

不妨设$A$A的列向量组为$c_1,c_2\cdots c_n$c1​,c2​⋯cn​，由其线性相关，可以得到存在一组不全为$0$0的系数$q_1,q_2\cdots q_n$q1​,q2​⋯qn​，满足$\sum q_ic_i=0$∑qi​ci​=0。不妨假设$q_c\not =0$qc​​=0

考虑同行两余子式$M_{r,i}$Mr,i​与$M_{r,c}$Mr,c​的关系（在此后的表述中，$M_{r,i}$Mr,i​可能为矩阵，也可能为该矩阵的行列式。视语境而定）。不失一般性地令$i<c$i<c，对于$i>c$i>c的情况是同理推导的

若我们将$M_{r,i}$Mr,i​的第$c-1$c−1列向前置换至第$i$i列，同时原本第$[i,c-2]$[i,c−2]列向后循环移位。我们能够得到与$M_{r,c}$Mr,c​十分相似的一个矩阵。唯一的不同在于$M_{r,i}$Mr,i​的第$i$i列为原本的第$c$c列去掉行$r$r，而$M_{r,c}$Mr,c​的第$i$i列为原本的第$i$i列去掉行$r$r

我们构造新矩阵$M'$M′，将$M_{r,i}$Mr,i​（循环移位后）的第$i$i列乘上$q_c$qc​，将$M_{r,c}$Mr,c​的第$i$i列乘上$q_i$qi​，令$M'$M′为这两个矩阵的和。这里两个矩阵的和的计算方式是这样的：将唯一不同的两列求和，其他所有的列照写。
由行列式的性质，不难发现有
$|M'|=(-1)^{c-i-1}q_c|M_{r,i}|+q_i|M_{r,c}|$∣M′∣=(−1)c−i−1qc​∣Mr,i​∣+qi​∣Mr,c​∣
若我们将$M'$M′的其他列乘上原本对应的列$c$c的$q_c$qc​，再加到$M'$M′的第$i$i列上，行列式不变。而第$i$i列全为$0$0，此时可以证明$|M'|=0$∣M′∣=0，因而有
$(-1)^{c-i}q_c|M_{r,i}|=q_i|M_{r,c}|\\ |M_{r,i}|=\frac{q_i}{q_c}(-1)^{c-i}|M_{r,c}|\\ (-1)^{r+i}|M_{r,i}|=\frac{q_i}{q_c}(-1)^{c-i+r+i}|M_{r,c}|\\ A_{r,i}=\frac{q_i}{q_c}A_{r,c}$(−1)c−iqc​∣Mr,i​∣=qi​∣Mr,c​∣∣Mr,i​∣=qc​qi​​(−1)c−i∣Mr,c​∣(−1)r+i∣Mr,i​∣=qc​qi​​(−1)c−i+r+i∣Mr,c​∣Ar,i​=qc​qi​​Ar,c​
通过对$i>c$i>c的情况加以分析，不难发现，上述结论仍然成立。故有
$A_{r,i}=\frac{q_i}{q_c}A_{r,c}$Ar,i​=qc​qi​​Ar,c​
对于行向量组的分析，也是相似的。

对于一组行向量$p=(p_1p_2\cdots p_n)$p=(p1​p2​⋯pn​)，列向量$q=(q_1q_2\cdots q_n)$q=(q1​q2​⋯qn​)，我们求出这样的两个向量使得$pA=0,Aq=0$pA=0,Aq=0。不妨设$p_r,q_c\not =0$pr​,qc​​=0，就有
$A_{i,j}=\frac{p_iq_j}{p_rq_c}A_{r,c}$Ai,j​=pr​qc​pi​qj​​Ar,c​
求解向量$p,q$p,q，代数余子式$A_{r,c}$Ar,c​均可以在$\Theta(n^3)$Θ(n3)的时间内用高斯消元法求解。综上所述，求解一个矩阵的伴随矩阵的时间复杂度为$\Theta(n^3)$Θ(n3)

代码

调用Mat::get_adjoint_matrix(a,n)即可求解$n$n阶矩阵$a$a的伴随矩阵
可以在这里验证你的代码
https://vjudge.net/problem/OpenJ_POJ-C19J

const int MAXN=35;
const int MAXM=152505;
const int mod=998244353;
inline void ad(int &x,int y){x+=y;if(x>=mod)x-=mod;}
inline void dl(int &x,int y){x-=y;if(x<0)x+=mod;}
inline int pow_mod(int a,int b)
{
	int ret=1;
	for(;b;b>>=1,a=1LL*a*a%mod)if(b&1)ret=1LL*ret*a%mod;
	return ret;
}
struct matrix{int m[MAXN][MAXN];matrix(){memset(m,0,sizeof(m));}};
namespace Mat
{
	int A[MAXN][MAXN],B[MAXN][2*MAXN];matrix mat_ret,mat_lin;
	int dat(const matrix &a,int n)
	{
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=a.m[i][j];
		int ret=1;
		for(int i=1;i<=n;i++)
		{
			int u=i;for(int j=i+1;j<=n;j++)if(A[j][i]){u=j;break;}
			if(u!=i){for(int j=i;j<=n;j++)swap(A[u][j],A[i][j]);ret=1LL*ret*(mod-1)%mod;}
			if(!A[i][i])return 0;int Iv=pow_mod(A[i][i],mod-2);
			for(int j=i+1;j<=n;j++)
			{
				int temp=1LL*A[j][i]*Iv%mod;
				for(int k=i;k<=n;k++)dl(A[j][k],1LL*A[i][k]*temp%mod);
			}ret=1LL*ret*A[i][i]%mod;
		}return ret;
	}
	matrix mat_inv(const matrix &a,int n)
	{
		memset(B,0,sizeof(B));
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)B[i][j]=a.m[i][j];
		for(int i=1;i<=n;i++)B[i][i+n]=1;
		for(int i=1;i<=n;i++)
		{
			int u=i;for(int j=i+1;j<=n;j++)if(B[j][i]){u=j;break;}
			for(int j=1;j<=2*n;j++)swap(B[i][j],B[u][j]);
			if(!B[i][i])continue;int Iv=pow_mod(B[i][i],mod-2);
			for(int j=i+1;j<=n;j++)
			{
				int temp=1LL*B[j][i]*Iv%mod;
				for(int k=i;k<=2*n;k++)dl(B[j][k],1LL*B[i][k]*temp%mod);
			}
		}
		for(int i=n;i>=1;i--)
		{
			int Iv=pow_mod(B[i][i],mod-2);
			for(int j=i;j<=2*n;j++)if(B[i][j])B[i][j]=1LL*B[i][j]*Iv%mod;
			for(int j=i-1;j>=1;j--)if(B[j][i])
			{
				int Z=B[j][i];
				for(int k=i;k<=2*n;k++)if(B[i][k])dl(B[j][k],1LL*B[i][k]*Z%mod);
			}
		}
		memset(mat_ret.m,0,sizeof(mat_ret.m));
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)
			mat_ret.m[i][j]=B[i][j+n];
		return mat_ret;
	}
	int calculate_r(const matrix &a,int n)
	{
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=a.m[i][j];int ret=0;
		for(int i=1;i<=n;i++)
		{
			if(!A[i][i])
			{
				int u=-1;
				for(int j=i+1;j<=n;j++)if(A[j][i]){u=j;break;}
				if(u==-1)continue;
				for(int j=i;j<=n;j++)swap(A[i][j],A[u][j]);
			}++ret;int Iv=pow_mod(A[i][i],mod-2);
			for(int j=i+1;j<=n;j++)
			{
				int temp=1LL*A[j][i]*Iv%mod;
				for(int k=i;k<=n;k++)dl(A[j][k],1LL*A[i][k]*temp%mod);
			}
		}return ret;
	}
	matrix get_adjoint_matrix_full(const matrix &a,int n)
	{
		matrix Iv=mat_inv(a,n);int Z=dat(a,n);
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)
			mat_ret.m[j][i]=1LL*Iv.m[i][j]*Z%mod;
		return mat_ret;
	}
	
	
	int us[MAXN][MAXN],Inv[MAXN];
	VI insert(VI Z,int n,int id)
	{
		VI bel(n+1,0);bel[id]=1;bel[0]=-1;
		for(int i=n;i>=1;i--)if(Z[i])
		{
			if(!A[i][i])
			{
				for(int j=i;j>=1;j--)A[i][j]=Z[j];Inv[i]=pow_mod(A[i][i],mod-2);
				for(int j=1;j<=n;j++)us[i][j]=bel[j];
				return bel;
			}
			int temp=1LL*Inv[i]*Z[i]%mod;
			for(int j=i;j>=1;j--)dl(Z[j],1LL*A[i][j]*temp%mod);
			for(int j=1;j<=n;j++)dl(bel[j],1LL*us[i][j]*temp%mod);
		}bel[0]=0;return bel;
	}
	VI get_G(const matrix &a,int n,int opt)
	{
		VI ret;
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=a.m[i][j];
		if(opt)for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)B[j][i]=A[i][j];
		else for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)B[i][j]=A[i][j];
		
		memset(A,0,sizeof(A));memset(us,0,sizeof(us));
		for(int i=1;i<=n;i++)
		{
			VI Z(n+1,0);
			for(int j=1;j<=n;j++)Z[j]=B[j][i];
			ret=insert(Z,n,i);
			if(ret[0]!=-1)return ret;
		}assert(false);
	}
	matrix get_adjoint_matrix_one(const matrix &a,int n)
	{
		VI Q=get_G(a,n,0);
		VI P=get_G(a,n,1);
		int c=-1,r=-1;
		for(int i=1;i<=n;i++)if(Q[i]){c=i;break;}assert(c!=-1);
		for(int i=1;i<=n;i++)if(P[i]){r=i;break;}assert(r!=-1);
		
		for(int i=1;i<=n;i++)if(i!=r)for(int j=1;j<=n;j++)if(j!=c)
			mat_lin.m[i-(i>r)][j-(j>c)]=a.m[i][j];
		int Ivq=pow_mod(Q[c],mod-2),Ivp=pow_mod(P[r],mod-2);
		int Iv=1LL*Ivq*Ivp%mod;
		mat_ret.m[r][c]=dat(mat_lin,n-1);
		if((r+c)&1)mat_ret.m[r][c]=mod-mat_ret.m[r][c];
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)
		{
			int val=1LL*P[i]*Q[j]%mod*Iv%mod*mat_ret.m[r][c]%mod;
			mat_ret.m[i][j]=val;
		}return mat_ret;
	}
	
	matrix get_adjoint_matrix(const matrix &a,int n)
	{
		int Z=calculate_r(a,n);
		if(Z==n)return get_adjoint_matrix_full(a,n);
		else if(Z==n-1)return get_adjoint_matrix_one(a,n);
		else
		{
			memset(mat_ret.m,0,sizeof(mat_ret.m));
			return mat_ret;
		}
	}
}