高斯消元-行列式-矩阵树定理学习笔记

高斯消元

引入

高斯消元听起来非常高大上,其实就是我们初中学的加减消元、代入消元法的程序化实现罢了。
考虑一个方程组 {A1,1x1+A1,2x2+A1,3x3+...+A1,nxn=C1A2,1x1+A2,2x2+A2,3x3+...+A2,nxn=C2...Ak,1x1+Ak,2x2+Ak,3x3+...+Ak,nxn=Cn 我们要求出它的一组解,或者判定无解或无穷多组解。这就是高斯消元的基础应用。

概念

系数矩阵:上面的方程组中系数A组成的n×n矩阵就是它的系数矩阵。
增广矩阵:系数矩阵的最右边补上一列,表示C1,C2,...,Cn,得到的n×(n+1)矩阵就成为增广矩阵。
初等行列变换
1. 以一个非零的数乘矩阵的某一行(列)
2. 把矩阵的某一行(列)的c倍加到另一行(列),这里c是任意数 3. 互换矩阵中两行(列)的位置

阶梯型矩阵:一个矩阵成为阶梯型矩阵,需满足两个条件:
1. 如果它既有零行,又有非零行,则零行在下,非零行在上。
2. 如果它有非零行,则每个非零行的第一个非零元素所在列号自上而下严格单调上升。
转置:把一个矩阵的行变成列,列变成行,得到的矩阵就是原矩阵的转置。记为AT
如图所示,这就是一个阶梯型矩阵
示例1

三角矩阵:三角矩阵分上三角矩阵和下三角矩阵两种。上三角矩阵的对角线左下方的系数全部为零,下三角矩阵的对角线右上方的系数全部为零。

算法流程

它的基本思想是通过初等行变换,把增广矩阵消成上三角矩阵,此时最后一个非零行必然是ax=b的形式。所以我们已经得到了一个未知数的值。把这个值代回上面的所有方程中,则它的上面一行也变成了ax=b的形式...重复下去,我们就可以得到所有未知数的值了。^_^

代码 (SDOI2006 线性方程组)

改编自 ComeIntoPower,在此对他表示谢意QwQ

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include <cstdio>
#include <cmath>
#include <algorithm>

using std::fabs;
using std::swap;

const int maxn=1e3+10;
const double eps=1e-9;

inline int Gauss_Elimination(double (*A)[maxn],double* f,int n)
{
for (int i=1,c=1,j;i<=n;++i)
{
for (j=c;j<=n && fabs(A[j][i])<eps;++j);
if (j==n+1) continue;
for (int k=1;k<=n+1;++k) swap(A[c][k],A[j][k]);
for (int j=c+1;j<=n;++j)
if (fabs(A[j][i])>eps)
{
double t=A[j][i]/A[c][i];
for (int k=i;k<=n+1;++k)
A[j][k]-=t*A[c][k];
}
++c;
}
bool NoAnswer=false,InfAnswer=false;
for (int i=n;i;--i)
{
bool NoVariables=true;
for (int j=i;j<=n;++j)
if (fabs(A[i][j])>eps) NoVariables=false;
if (NoVariables)
if (fabs(A[i][n+1])>eps) NoAnswer=true;// 0=C,C!=0,无解
else InfAnswer=true;// 0=0,无穷多组解
else
{
for (int j=i+1;j<=n;++j) A[i][n+1]-=A[i][j]*f[j];
f[i]=A[i][n+1]/A[i][i];
}
}
if (NoAnswer) return -1;
return !InfAnswer; //无穷多解返回0,有唯一解返回1
}

int main()
{
static double A[maxn][maxn],f[maxn];
int n;
scanf("%d",&n);
for (int i=1;i<=n;++i)
for (int j=1;j<=n+1;++j)
scanf("%lf",&A[i][j]);
int result=Gauss_Elimination(A,f,n);
if (result^1) return printf("%d\n",result)&0;
for (int i=1;i<=n;++i) printf("x%d=%.2lf\n",i,f[i]);
}

解释:

    • Q: 咋判断无解和无穷多组解鸭?
    • A:当消元完成后,回代之前时,方程组若出现除了常数项之外全零的行,则一定无解或无穷多组解。(常数项为0是无穷多组解,非零则无解)
    • Q:代码中那个变量c的作用是啥呀?为啥不能用i代替呢?
    • A:因为如果把c换成i,当存在一个自由元的时候,可能没法消出来全零行,而是直接跳过了本该消成全零行的那一行,将其放在了原来的位置,导致判断无穷解和无解的时候出偏差。(详见洛谷本题讨论的“Hack+1”篇目).

行列式

概念:

主子式:对于一个n阶矩阵A,选取它的任意i行,将行号记为T[1..i],并同时选取第T[1..n]列,得到的新矩阵即为原矩阵的一个i阶主子式。 新矩阵的行列式也可称为|A|的一个i阶主子式。

行列式:行列式是一个标量。对一个矩阵A来说,它的行列式(记作|A|)的定义是 |A|=P(1)δ(P)k=1nAi,Pk 其中,P取遍1..n的所有排列,δ(P)为P的逆序对数。这个柿子很好理解:就是从每行每列各选一个数,一共n个数, 行列式有以下性质:
1. |A|=|AT| (由定义,比较显然,在此不证) 2. 将一个矩阵的两行互换,行列式变号。
为了证明这个定理,我们需要首先证明一个引理: * 引理:一个排列的两项交换,逆序对改变量为奇数。 * 引理证明:不妨设a<b.这样的交换对位置在a前面和在b后面的数没有影响,所以只需要考虑a,b之间的这个子区间。设交换的两项为a,b.把位置在a.b之间的数分类:设属于(,a)的数量为x,属于(a,b)的数量为y,属于(b,+)的数量为z。区间里交换以前的逆序对数为x+z,交换后为(x+y)+(y+z)+1,差值为2y+1,是奇数。

现在,证明这个定理就十分容易了。

  • 证明:交换矩阵的两行,相当于交换P的两项。逆序对改变量为奇数,所以行列式变号。

  • 推论2.1:有两行相同的矩阵,行列式为0

    • 证明:交换这两行,行列式变号且值不变...
  1. 将某一行乘上k,行列式乘上k.
    • 证明:你在这一行选出的每个数,都乘上了一个k。提取公因数即得。
  2. 两个矩阵如果只有一行不同,则它们的行列式和等于将不同的行相加得到的新矩阵的行列式。
    • 证明:从定义下手。挺显然的吧qwq.
  • 推论4.1:将一行乘上k的值加到另一行上,行列式不变。
    • 证明:设被加的行为x,乘上k加到A上的为y。把这个新矩阵拆成两个矩阵A,B。其中A是原矩阵,Bx行改为yk:则由性质2和3,B的行列式为0。又由性质4,新矩阵的行列式为|A|+|B|,和原来相等
    • 有了这个推论和性质2,我们就可以做高斯消元啦!
  1. 每行每列和均为0的矩阵,行列式为0.
    • 证明:对原矩阵补上一列0,进行高斯消元。显然这个方程组有解——所有未知数都相等.则由高斯消元判断无穷多组解的条件可得,消出的上三角矩阵中,一定有全零行,且常数项也是0。所以该矩阵的行列式为0.因为高消之后的矩阵的行列式只是原行列式乘上±1(只在交换行的时候行列式发生变号),所以原矩阵的行列式也为0.
  2. 对于一个上三角矩阵,它的行列式为对角线上数的乘积。
    • 证明:要在每行选取一个列不重复的数,且这些数乘积非0,则只能选对角线上的n个数。

求解

直接求解要枚举全排列,时间复杂度为O(n!n)。不可接受。使用高斯消元将其化为上三角矩阵,可在O(n3)的时间内求n阶矩阵的行列式。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
inline long double determinant(long double (*A)[maxn],int n)
{
int s=1;
for (int i=1,c=1,j;i<=n;++i)
{
for (j=c;j<=n && fabs(A[j][i])<eps;++j);
if (j==n+1) continue;
s=-s;
for (int k=1;k<=n;++k) swap(A[c][k],A[j][k]);
for (int j=c+1;j<=n;++j)
if (fabs(A[j][i])>eps)
{
long double t=A[j][i]/A[c][i];
for (int k=i;k<=n;++k)
A[j][k]-=A[c][k]*t;
}
++c;
}
long double ans=s;
for (int i=1;i<=n;++i)
ans*=A[i][i];
return ans; // 这里要取绝对值
}

矩阵树定理

概念:

  1. 度数矩阵:定义D为图G=(V,E),(|V|=n,|E|=m)的度数矩阵,则D为一个n×n矩阵,其中Di,i为编号为i的结点的度数.
  2. 邻接矩阵:就是我们通常所说的邻接矩阵,记为G
  3. 基尔霍夫矩阵(拉普拉斯矩阵):定义基尔霍夫矩阵K=DG.

定理:

当邻接矩阵不带边权时(若u,v有边,则G(u,v)=G(v,u)=1),则有

矩阵树定理:一个无重边、自环的图G的生成树个数,等于它的基尔霍夫矩阵任意一个n-1阶主子式的行列式的绝对值。

将邻接矩阵加上边权,得到新的邻接矩阵、度数矩阵和基尔霍夫矩阵。即:允许重边(甚至把边数扩展到R) 此时: D(i,i)=(i,j)Ew((i,j)) G(i,j)=w(i,j) K=DG 推广一下,有:

变元矩阵树定理: 1. 对生成树T定义求其边权之积的函数F(T)=eTw(e) 2. 对于每一棵生成树,求其F函数值的和得到H(G):H(G)=T is a spanning tree of GF(T) 3. 则H(G)等于(带边权的)基尔霍夫矩阵的任意一个n-1阶主子式的行列式的绝对值. 4. 容易发现,当边权都为1的时候,它就是普通的矩阵树定理.

对于有向图来说,有:

有向图的矩阵树定理
定义 Ki,j=w(i,j) Ki,i=j=1nw(i,j) 去掉第i行第i列,则能求出以i为根的外向树的数量(边权积)
同样地,定义 Ki,j=w(i,j) Ki,i=j=1nw(j,i) 则能求内向树的数量(边权积)

(原谅我语文学得不好没法简洁地描述上述定理qwq...) ## 证明: ~~我们采用闭眼证明法...嗯!它是对的!^_^..~~
事实上博主太菜并不会证明..(留坑)

应用:

BZOJ1002(假的)题解戳这里~
BZOJ1016题解戳这里~
BZOJ3524题解戳这里~

引用&鸣谢

  • 高斯消元的模板来自CDQZ ComeIntoPower巨佬
  • 感谢Memset0巨佬Hack掉我原先的高斯消元。
  • 矩阵树定理一篇大量来源于Candy?的博客和ReMoon。(其实就是抄的吧!

对以上巨佬表示感谢。

Powered By Valine
v1.5.2