前言
在联合省选day2t3中,存在一种使用行列式求导来计算生成树边权和的方法。需要计算出每个位置的代数余子式
常见的做法均是套用
A
×
A
∗
=
∣
A
∣
×
I
A\times A^*=|A|\times I
A×A∗=∣A∣×I的等式求解,这种做法在模意义下矩阵
r
a
n
k
(
A
)
=
n
−
1
rank(A)=n-1
rank(A)=n−1时不能正确得出伴随矩阵
以下给出正确的求解伴随矩阵的方法
代数余子式
给 n n n阶方阵 A = ( a i , j ) A=(a_{i,j}) A=(ai,j),定义 a i , j a_{i,j} ai,j的余子式 M i , j M_{i,j} Mi,j为将 i i i行 j j j列划去后的行列式, a i , j a_{i,j} ai,j的代数余子式 A i , j A_{i,j} Ai,j为 ( − 1 ) i + j M i , j (-1)^{i+j}M_{i,j} (−1)i+jMi,j
伴随矩阵
若
A
i
,
j
A_{i,j}
Ai,j为代数余子式矩阵,我们定义
A
∗
=
(
A
i
j
)
T
A^*=(A_{ij})^{\mathrm T}
A∗=(Aij)T
A
∗
A^*
A∗为矩阵
A
A
A的伴随矩阵
我们存在一个结论,即
A
×
A
∗
=
A
∗
×
A
=
∣
A
∣
×
I
n
A\times A^*=A^*\times A=|A|\times I_n
A×A∗=A∗×A=∣A∣×In
其中
I
n
I_n
In为
n
n
n阶单位矩阵,
∣
A
∣
|A|
∣A∣为矩阵
A
A
A的行列式
计算
若 r a n k ( A ) = n rank(A)=n rank(A)=n,则 A − 1 A^{-1} A−1是良定义的,可以直接用 A ∗ = ∣ A ∣ × A − 1 A^*=|A|\times A^{-1} A∗=∣A∣×A−1来进行计算。
若 r a n k ( A ) ≤ n − 2 rank(A)\leq n-2 rank(A)≤n−2,不难发现任意 n − 1 n-1 n−1阶子式行列式均为 0 0 0,故 A ∗ = O A^*=O A∗=O,其中 O O O为空矩阵
我们只需要特殊考虑 r a n k ( A ) = n − 1 rank(A)=n-1 rank(A)=n−1的情况
不妨设 A A A的列向量组为 c 1 , c 2 ⋯ c n c_1,c_2\cdots c_n c1,c2⋯cn,由其线性相关,可以得到存在一组不全为 0 0 0的系数 q 1 , q 2 ⋯ q n q_1,q_2\cdots q_n q1,q2⋯qn,满足 ∑ q i c i = 0 \sum q_ic_i=0 ∑qici=0。不妨假设 q c ≠ 0 q_c\not =0 qc=0
考虑同行两余子式 M r , i M_{r,i} Mr,i与 M r , c M_{r,c} Mr,c的关系(在此后的表述中, M r , i M_{r,i} Mr,i可能为矩阵,也可能为该矩阵的行列式。视语境而定)。不失一般性地令 i < c i<c i<c,对于 i > c i>c i>c的情况是同理推导的
若我们将 M r , i M_{r,i} Mr,i的第 c − 1 c-1 c−1列向前置换至第 i i i列,同时原本第 [ i , c − 2 ] [i,c-2] [i,c−2]列向后循环移位。我们能够得到与 M r , c M_{r,c} Mr,c十分相似的一个矩阵。唯一的不同在于 M r , i M_{r,i} Mr,i的第 i i i列为原本的第 c c c列去掉行 r r r,而 M r , c M_{r,c} Mr,c的第 i i i列为原本的第 i i i列去掉行 r r r
我们构造新矩阵
M
′
M'
M′,将
M
r
,
i
M_{r,i}
Mr,i(循环移位后)的第
i
i
i列乘上
q
c
q_c
qc,将
M
r
,
c
M_{r,c}
Mr,c的第
i
i
i列乘上
q
i
q_i
qi,令
M
′
M'
M′为这两个矩阵的和。这里两个矩阵的和的计算方式是这样的:将唯一不同的两列求和,其他所有的列照写。
由行列式的性质,不难发现有
∣
M
′
∣
=
(
−
1
)
c
−
i
−
1
q
c
∣
M
r
,
i
∣
+
q
i
∣
M
r
,
c
∣
|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'
M′的其他列乘上原本对应的列
c
c
c的
q
c
q_c
qc,再加到
M
′
M'
M′的第
i
i
i列上,行列式不变。而第
i
i
i列全为
0
0
0,此时可以证明
∣
M
′
∣
=
0
|M'|=0
∣M′∣=0,因而有
(
−
1
)
c
−
i
q
c
∣
M
r
,
i
∣
=
q
i
∣
M
r
,
c
∣
∣
M
r
,
i
∣
=
q
i
q
c
(
−
1
)
c
−
i
∣
M
r
,
c
∣
(
−
1
)
r
+
i
∣
M
r
,
i
∣
=
q
i
q
c
(
−
1
)
c
−
i
+
r
+
i
∣
M
r
,
c
∣
A
r
,
i
=
q
i
q
c
A
r
,
c
(-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∣=qcqi(−1)c−i∣Mr,c∣(−1)r+i∣Mr,i∣=qcqi(−1)c−i+r+i∣Mr,c∣Ar,i=qcqiAr,c
通过对
i
>
c
i>c
i>c的情况加以分析,不难发现,上述结论仍然成立。故有
A
r
,
i
=
q
i
q
c
A
r
,
c
A_{r,i}=\frac{q_i}{q_c}A_{r,c}
Ar,i=qcqiAr,c
对于行向量组的分析,也是相似的。
对于一组行向量
p
=
(
p
1
p
2
⋯
p
n
)
p=(p_1p_2\cdots p_n)
p=(p1p2⋯pn),列向量
q
=
(
q
1
q
2
⋯
q
n
)
q=(q_1q_2\cdots q_n)
q=(q1q2⋯qn),我们求出这样的两个向量使得
p
A
=
0
,
A
q
=
0
pA=0,Aq=0
pA=0,Aq=0。不妨设
p
r
,
q
c
≠
0
p_r,q_c\not =0
pr,qc=0,就有
A
i
,
j
=
p
i
q
j
p
r
q
c
A
r
,
c
A_{i,j}=\frac{p_iq_j}{p_rq_c}A_{r,c}
Ai,j=prqcpiqjAr,c
求解向量
p
,
q
p,q
p,q,代数余子式
A
r
,
c
A_{r,c}
Ar,c均可以在
Θ
(
n
3
)
\Theta(n^3)
Θ(n3)的时间内用高斯消元法求解。综上所述,求解一个矩阵的伴随矩阵的时间复杂度为
Θ
(
n
3
)
\Theta(n^3)
Θ(n3)
代码
调用Mat::get_adjoint_matrix(a,n)即可求解
n
n
n阶矩阵
a
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;
}
}
}