代数余子式和伴随矩阵
代数余子式
给定 \(n\) 阶方阵 \(A=(a_{ij})\),定义 \(a_{ij}\) 的余子式 \(M_{ij}\) 为 \(A\) 划去第 \(i\) 行第 \(j\) 列后的行列式,\(a_{ij}\) 的代数余子式 \(A_{ij}=(−1)^{i+j}M_{ij}\) 。
代数余子式可以用于行列式的求值,比如按第 \(r\) 行展开:
按第 \(c\) 列展开是同理的。
伴随矩阵
定义
在 \((1)\) 式中,如果把 \(a_{rc}\) 中的 \(r\) 替换成 \(i≠r\),则该乘积对应了将 \(A\) 的第 \(i\) 行替换成第 \(r\) 行后的行列式——该行列式有两行相等,所以它等于 \(0\) 。
所以我们有 \(\sum_{c=1}^na_{ic}A_{jc}=\det A \cdot [i=j]\),把它写成矩阵乘法就是:
我们定义 \(A^*=(A_{ij})^{\mathrm T}\) 是 \(A\) 的伴随矩阵。
同样的结论对列也成立,所以 \(AA^*=A^*A=\det A \cdot I_n\) 。
计算
-
如果 \(r(A)=n\),用高斯消元法分别求出 \(\det A\) 和 \(A−1\),\(A^*=\det A \cdot A^{-1}\) 。
-
如果 \(r(A)⩽n−2\),\(A\) 的任意 \(n−1\) 阶子式都为 \(0\),所以 \(A^*=O\) 。
-
如果 \(r(A)=n−1\):
假设 \(A\) 的列向量组为 \(\vec{c}_{1}, \vec{c}_{2}, \ldots, \vec{c}_{n}\) 。
由它们线性相关,可知存在不全为 \(0\) 的实数 \(q_{1}, q_{2}, \ldots, q_{n}\) 满足 \(\sum_{i=1}^{n} q_{i} \vec{c}_{i}=\overrightarrow{0}\) 。不妨设 \(q_{c} \neq 0\) 。
考虑同行两余子式 \(M_{r i}\) 和 \(M_{r c}\) 的关系,不妨讨论 \(i<c\) 的情况,\(i>c\) 类似。不难发现,如果我们将 \(M_{r i}\) 的第 \(i\) 至 \(c-2\) 列向右推移,而第 \(c-1\) 列移动到第 \(i\) 列,唯一的不同点是 \(M_{r i}\) 的第 \(i\) 列为 \(\vec{c}_{c}\) 删去第 \(r\) 行,\(M_{r c}\) 的第 \(i\) 列为 \(\vec{c}_{i}\) 删去第 \(r\) 行。
联系上述线性相关式,可设 \(M^{\prime}\) 表示在 \(M_{r i}\) 进行推移的基础上,把第 \(i\) 列乘以 \(q_{c \prime}\) 并与 \(M_{r c}\) 的第 \(i\) 列的 \(q_{i}\) 倍相加,所得到的 \(n-1\) 阶行列式。
那么 \(M^{\prime}=q_{c}(-1)^{c-i} M_{r i}+q_{i} M_{r c}\) 。
最后,将 \(M^{\prime}\) 的第 \(i\) 列依次加上其对应 “原行列式的第 \(j \notin\{i, c\}\) 列” 的列的 \(q_{j}\) 倍。
这样,\(M^{\prime}\) 的值不变,而第 \(i\) 列全为 \(0\) 。这说明了 \(q_{c}(-1)^{c-i} M_{r i}+q_{i} M_{r c}=M^{\prime}=0\) 。
所以
即
同样的结论对列也成立。
于是对于列向量 \(\vec{p}=\left(p_{1} p_{2} \cdots p_{n}\right)^{\mathrm{T}}\) 和行向量 \(\vec{q}=\left(q_{1} q_{2} \cdots q_{n}\right)\) 使得 \(A\vec{p}=\overrightarrow{0}, \vec{q} A=\overrightarrow{0}\) ,不妨设 \(p_{r}, q_{c} \neq 0\) ,就有 \(A_{i j}=\frac{p_{i} q_{j}}{p_{r} q_{c}} A_{r c}\) 。
\(\vec{p}\) 与 \(\vec{q}\) 都可以通过高斯消元法解方程组解出,而代数余子式 \(A_{r c}\) 也可以通过高斯消元法求出。
综上所述,求一个矩阵的伴随矩阵,时间复杂度为 \(O\left(n^{3}\right)\) 。
应用
通过伴随矩阵,可以求出矩阵所有元素的余子式。
HDU7253 Tree
给定有向图,求每条边各含于多少棵以 \(R\) 为根的所有生成内向树。
solution:
计算含有 \(\langle u, v\rangle\) 的生成内向树个数,相当于把 \(u\) 的出边除了 \(\langle u, v\rangle\) 外全部删去后的生成内向树计数。
这就是把拉普拉斯矩阵的第 \(u\) 行替换成仅有第 \(u\) 项为 \(1\) ,第 \(v\) 项为 \(-1\) 的行,然后划去第 \(R\) 行第 \(R\) 列求余子式。
这等价于,拉普拉斯矩阵去掉第 \(u, R\) 行、第 \(u, R\) 列的余子式,减去去掉第 \(u, R\) 行、第 \(v, R\) 列的余子式。
将拉普拉斯矩阵去掉第 \(R\) 行、第 \(R\) 列作为矩阵 \(A\) ,那么就转化为求 \(A\) 的所有余子式。
这里不用讨论那么多情况,因为 \(\mathrm{r}(A)<n\) 意味着不存在生成内向树。
点击查看代码
#include <bits/stdc++.h>
using namespace std;
int T;
int n, m;
const int md = 1e9 + 7;
inline int pwr(int x, int y) {
int res = 1;
while (y) {
if (y & 1)
res = 1ll * res * x % md;
x = 1ll * x * x % md;
y >>= 1;
}
return res;
}
struct mat {
int a[505][505];
mat(int _x = 0) {
memset(a, 0, sizeof(a));
if (_x)
for (int i = 1; i <= n; i++) a[i][i] = 1;
}
inline auto operator[](int t) { return a[t]; }
inline int Det(int N) {
int res = 1;
for (int i = 1; i <= N; i++) {
for (int j = i + 1; j <= N; j++) {
while (a[j][i]) {
int tmp = a[i][i] / a[j][i];
for (int t = 1; t <= N; t++) a[i][t] = (a[i][t] - 1ll * tmp * a[j][t]) % md;
for (int t = 1; t <= N; t++) swap(a[i][t], a[j][t]);
res = -res;
}
}
}
for (int i = 1; i <= N; i++) res = 1ll * res * a[i][i] % md;
return (res + md) % md;
}
inline mat Inv() {
mat res(1);
for (int i = 1; i <= n; i++) {
if (!a[i][i]) {
int loc = i;
for (int j = 1; j <= n; j++)
if (a[j][i])
loc = j;
if (i != loc) {
for (int j = 1; j <= n; j++) swap(res[i][j], res[loc][j]);
for (int j = 1; j <= n; j++) swap(a[i][j], a[loc][j]);
}
}
for (int j = 1; j <= n; j++) {
if (i != j && a[j][i]) {
int tmp = 1ll * a[j][i] * pwr(a[i][i], md - 2) % md;
for (int t = 1; t <= n; t++) a[j][t] = (a[j][t] - 1ll * tmp * a[i][t]) % md;
for (int t = 1; t <= n; t++)
res[j][t] = (res[j][t] - 1ll * tmp * res[i][t] % md + md) % md;
}
}
}
for (int i = 1; i <= n; i++)
if (!a[i][i]) {
for (int j = 1; j <= n; j++) swap(res[i][j], res[n][j]);
for (int j = 1; j <= n; j++) swap(a[i][j], a[n][j]);
break;
}
return res;
}
int x[505], y[505];
inline mat calc() {
mat tmp = *this;
tmp = tmp.Inv();
for (int i = 1; i <= n; i++) x[i] = tmp[n][i];
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) tmp[i][j] = a[j][i];
}
tmp = tmp.Inv();
for (int i = 1; i <= n; i++) y[i] = tmp[n][i];
int r = 0, c = 0;
for (int i = 1; i <= n; i++)
if (x[i])
r = i;
for (int i = 1; i <= n; i++)
if (y[i])
c = i;
mat res;
for (int i = 1; i <= n; i++) {
if (i == r)
continue;
for (int j = 1; j <= n; j++)
if (j != c)
tmp[i - (i > r)][j - (j > c)] = a[i][j];
}
res[r][c] = tmp.Det(n - 1);
if ((r + c) & 1)
res[r][c] = md - res[r][c];
int T = 1ll * res[r][c] * pwr(1ll * x[r] * y[c] % md, md - 2) % md;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
if (i != r || j != c)
res[i][j] = 1ll * x[i] * y[j] % md * T % md;
}
}
return res;
}
} mp;
inline void solve() {
scanf("%d%d", &n, &m);
mp = mat();
for (int i = 1; i <= m; i++) {
int x, y;
scanf("%d%d", &x, &y);
mp[x][y]--;
mp[y][y]++;
}
mp = mp.calc();
for (int i = 1; i <= n; i++) printf("%d ", mp[i][i]);
puts("");
}
int main() {
scanf("%d", &T);
while (T--) solve();
return 0;
}