概率期望中高斯消元的几种用法

前置知识：高斯消元法
博主理解浅显，只能膜piao别人的总结
戳别人家的题解
咳咳……还是简单介绍两句
它可以用 $O(n^3)$ 的复杂度解出 n 元方程组
表示方法：矩阵
tips：一般情况下高斯消元可能出现无解、无穷解的情况，我的做法里面没有判断，由于矩阵对角线上不会出现0。

概率与期望：

概率：发生的可能性
期望：概率的加权平均数（表示对权值的一个预期值）
eg. 某图中从起点经过 i 步到达终点的可能性为p[i]，则到达终点的期望步数为$\sum p[i]\ast i$

对于期望的小理解：
其实很简单，在物理中，我们学习过重心的概念，可以把一个物体当作一个质点来考虑，对于期望其实有着类似的道理。期望与权值本身是不同的，但是期望表达的是权值的一个大概分布位置，权值可能有很多很多个，但是期望只有一个。同样可以说，重心就是物质质量的期望分布位置。

期望的好处在于，当我们求的东西（如期望）与物体本身权值无关时，我们可以用期望代替计算这个物体的权值。（想想那句话，当研究对象与物体的形状（质量分布）无关时，我们可以把它当作质点处理）
于是我们可以借助物理的知识来比较好的理解期望！

用一道题来结合期望与高斯消元

[Hnoi2013]游走（高斯消元）

不解释题意差评
vjudge
bzoj
还是解释一下吧
无向图，随机游走，从1号点走到$n$号点 ，给每条边编号，经过代价即为编号，求最小期望代价。

假如能够求出每条边的期望经过次数，剩下的部分就很水了（排序编号直接算）。

关键是求每条边的期望经过次数。
发现：除了终点以外，其余的点只要被经过1次，就必然会出去一次，出去的这一次由于是随机游走，则有可能对这个点连接的任意一条边造成贡献。
设一条边 $(u,v)$连接了不是终点的两个点$u、v$
$siz[u]$表示u连接了$siz[u]$条边
$f[u]$表示$u$点的期望经过次数。
(这里的经过指一定要从那个点出来的经过)
$$E(u,v)=\frac{f[u]}{siz[u]}+\frac{f[v]}{siz[v]}$$
对于终点：直接$f[n]=0$就好了

关键是求每个点的期望经过次数。
对于每个点，它们都有可能从任何一个相邻点转移过来
因此有 $f[u]=\sum \frac{f[v]}{siz[v]}(v与u相邻)$
此时想要单独解出任何一个点都是不行的，但是可以发现我们得到了$n-1$个（除去终点）方程组，再加上f[n]=0其实一共是n个，也就是解的出来。
构造矩阵进行高斯消元即可。

（丢出一只又奇怪又慢的代码，$Init$打的之丑陋，泪奔.jpg）
fc数组就是构造出来的方程矩阵，$Init$是把它传进去……
主要是这样子的话就可以连着一个$struct$直接Ctrl c下来粘到另一份代码上面搞高斯消元了（懒人操作）
（如果数组大小不同$Init$里x+=N；那里需要修改）
主函数专心构造$fc$矩阵就好啦，把它丢进去就不管了。

#include<cstdio>
#include<vector>
#include<algorithm>
using namespace std;
typedef long long LL;
const int INF=0x3f3f3f3f;
const int N=505;
const int M=N*N/2;
#define eps 1e-7
int read(){
	bool f=0;int x=0;char c=getchar();
	while(c<'0'||'9'<c){if(c=='-')f=1;c=getchar();}
	while('0'<=c&&c<='9') x=(x<<3)+(x<<1)+(c^48),c=getchar();
	return !f?x:-x;
}
struct GAUSS {
	int n; double *a[N];
	double ans[N];
	void Init(int _n,double *x) {
		n=_n;
		for(int i=0;i<=n;i++) {
			a[i]=x;
			x+=N;
		}
	}
	inline double Abs(double x){return x>0?x:-x;}
	bool Guass(){
		double tmp;
		for(int i=1;i<=n;i++) {
			int p=i;
			for(int j=i+1;j<=n;j++)
				if(Abs(a[j][i])>Abs(a[p][i]))
					p=j;
			swap(a[p],a[i]);
			if(Abs(a[i][i])<eps)
				return false;
			for(int j=1;j<=n;j++)
				if(j!=i) {
					tmp=a[j][i]/a[i][i];
					a[j][i]=0;
					for(int k=i+1;k<=n+1;k++)
						a[j][k]-=tmp*a[i][k];
				}
		}
		for(int i=1;i<=n;i++)
			ans[i]=a[i][n+1]/a[i][i];
		return true;
	}
}GE;
int n,m,siz[N];
vector<int>G[N];
double fc[N][N];
int s[M],t[M];
double edge[M];
int main(){
	n=read(); m=read();
	for(int i=1;i<=m;i++) {
		s[i]=read(),t[i]=read();
		G[s[i]].push_back(t[i]);
		G[t[i]].push_back(s[i]);
		siz[s[i]]++,siz[t[i]]++;
	}
	for(int u=1;u<n;u++) {
		for(int i=0;i<siz[u];i++)
			fc[u][G[u][i]]=-1.0/siz[G[u][i]];
		fc[u][u]=1;
	}
	siz[n]=1; fc[n][n]=1;
	fc[1][n+1]=1;
	GE.Init(n,&fc[0][0]);
	GE.Guass();
	for(int i=1;i<=m;i++)
		edge[i]=GE.ans[s[i]]/siz[s[i]]+GE.ans[t[i]]/siz[t[i]];
	sort(edge+1,edge+m+1);
	double tot=0;
	for(int i=1;i<=m;i++)
		tot+=edge[i]*(m-i+1);
	printf("%.3lf\n",tot);
	return 0;
}

扯明白了一道题，可以再扯一道

[CF963E]Circles of Waiting（带状矩阵）

从坐标轴上$(0,0)$开始，已知向上左下右的概率是$p1,p2,p3,p4$（满足$p1+p2+p3+p4=1$）
问期望走多少步后与原点距离大于r $(x^2+y^2>r^2)$
$(r<=50)$
$$f[i]=p1\ast (f[上]+1)+p2\ast (f[左]+1)+p3\ast (f[下]+1)+p4\ast (f[右]+1)$$
$$f[i]=p1\ast f[上]+p2\ast f[左]+p3\ast f[下]+p4\ast f[右]+1$$

假设半径为$r$ 的圆内有n个点，那么必然有$n$个方程
但是半径为50的圆内有7000+的点，$O(n^3)$的级别太大

将所有点按照横纵坐标排序编号，发现任意一个点和它的转移点的编号差为最大为$2r$，这时候构造出来的方程矩阵为：
第$i$行表示编号为$i$的点列出来的方程，此处 $d=2r$，灰色部分都不可能会有值，因为第i行的任何一个元与 $i$ 的编号差不会大于$d$
这样的矩阵被称为 带状矩阵（常见于网格图！）

面对一个带状矩阵，首先明确我们最终目的是消除矩阵内的元素到只留下对角线，特别注意不能破坏矩阵的带状性质（交换两行）。
先考虑消除红色部分：
枚举每一行，使用这一行的去消以后的行：
可以发现当前行只需要枚举$d$个（因为这一行后面都是0了）
还可以发现往后只需要枚举$d$行（因为再下面的这一列也没了）

这样子做下来是 $O(d^{2}n)$
变成这样：

此时第$n$行已经只有一个元，得到解，把它回代到前面的$d$个方程里面去即可得到每一个解，这个过程$O(dn)$
细节是常数项那一列不在带状里面，运算需要单独写出来。

代码还是丢出来：

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
#define ID(x,y) id[x+R][y+R]
const int MOD=1e9+7;
const int R=52;

int r,a1,a2,a3,a4;
int id[R*2][R*2],idcnt;
LL p1,p2,p3,p4;
LL Pow(LL x,int p) {
	if(p==0) return 1;
	if(p&1) return x*Pow(x*x%MOD,p>>1)%MOD;
	return Pow(x*x%MOD,p>>1);
}
int fc[8000][8000];
struct GAUSS {
	int n; int *a[8000];
	void Init(int _n,int *x) {
		n=_n;
		for(int i=0;i<=n;i++) {
			a[i]=x;
			x+=8000;
		}
	}
	inline double Abs(double x){return x>0?x:-x;}
	LL Guass(){
		LL tmp;
		int d=2*r;
		for(int i=1;i<=n;i++) {
			//if(Abs(a[i][i])<eps) return false;
			LL inv=Pow(a[i][i],MOD-2);
			for(int j=i+1;j<=i+d&&j<=n;j++)
				if(j!=i&&a[j][i]) {
					tmp=a[j][i]*inv%MOD;
					a[j][i]=0;
					for(int k=i+1;k<=i+d&&k<=n;k++)
						a[j][k]=(a[j][k]-tmp*a[i][k])%MOD;
					a[j][n+1]=(a[j][n+1]-tmp*a[i][n+1])%MOD;
				}
		}
		int goal=ID(0,0);
		for(int i=n;i>=goal;i--) {
			LL inv=Pow(a[i][i],MOD-2);
			int &t=a[i][n+1];
			t=inv*t%MOD;
			for(int j=i-1;j>=goal;j--)
				a[j][n+1]=(a[j][n+1]-(LL)a[j][i]*t)%MOD;
		}
		return a[goal][n+1];
	}
}GE;
int main() {
	scanf("%d %d %d %d %d",&r,&a1,&a2,&a3,&a4);
	LL inv=Pow(((LL)a1+a2+a3+a4)%MOD,MOD-2);
	p1=a1*inv%MOD; p2=a2*inv%MOD; p3=a3*inv%MOD; p4=a4*inv%MOD;
	for(int x=-r;x<=r;x++)
		for(int y=-r;y<=r;y++)
			if(x*x+y*y<=r*r)
				ID(x,y)=++idcnt;
	for(int x=-r;x<=r;x++)
		for(int y=-r;y<=r;y++)
			if(int f=ID(x,y)) {
				fc[f][f]=-1;
				if(int t=ID(x-1,y))
                    fc[f][t]=p1;
                if(int t=ID(x,y-1))
                    fc[f][t]=p2;
                if(int t=ID(x+1,y))
                    fc[f][t]=p3;
                if(int t=ID(x,y+1))
                    fc[f][t]=p4;
                fc[f][idcnt+1]=-1;
			}
	GE.Init(idcnt,&fc[0][0]);
	printf("%d\n",(GE.Guass()+MOD)%MOD);
	return 0;
}

做完题后：说实话确实没有想到去考虑什么无解、无穷解的东西
毕竟……你说一个期望它无解或者无穷解会怎么样
一脸懵逼.jpg

[Topcoder12984]TorusSailing（主元法）

vjudge
Topcoder

在这个东西上面再次开始向右或向下等概率随机游走，问从$(0,0)$走到(gx,gy)的期望步数，$n、m$为行列数
$（n、m<=100）$
（强烈建议化终点为最右下角的点，起点按相对位置平移）

刚看到这个网格图，令$f[i]$表示编号为i的点走到终点的期望步数

$$f[i]=(f[右]+f[下])/2+1$$
我刷的一下掏出了带状矩阵：从左到右从上到下编号，$d=m$
但是光是元的个数就达到了$nm$级别，这个的平方的内存吃不消啊。况且这个矩阵的左下角貌似有毒——它破坏了带状性质后，原先那样$O(d^{2}s)$的做法就不行了，还是得$O(s^3)$，送命了
（其实我在想有没有一种编号方法可以保留带状性质，但我没构造出来……咕了它吧）

回到题目背景，可以发现一点就是列出来的方程十分简单，而且只跟右下角方向的有关，如果只设下图中的红色部分为元，则所有格子都可以被表示出来


单次表示的复杂度是$O(n+m)$，转移$O(nm)$次，复杂度$O(n^3)$
这样子表示出来之后，对于每个红色格子的转移式来列方程，元的个数被优化到了$O(n)$级别，直接高斯消元也就可以接受了。

#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long double LD;
#define eps 1e-9
const int N=102;

using namespace std;
struct GAUSS {
	int n; LD *a[N<<1];
	LD ans[N<<1];
	void Init(int _n,LD *x) {
		n=_n;
		for(int i=0;i<=n;i++) {
			a[i]=x;
			x+=(N<<1);
		}
	}
	inline LD Abs(LD x){return x>0?x:-x;}
	bool Guass() {
		LD tmp;
		for(int i=1;i<=n;i++) {
			int p=i;
			for(int j=i+1;j<=n;j++)
				if(Abs(a[j][i])>Abs(a[p][i]))
					p=j;
			swap(a[p],a[i]);
			if(Abs(a[i][i])<eps)
				return false;
			for(int j=1;j<=n;j++)
				if(j!=i) {
					tmp=a[j][i]/a[i][i];
					a[j][i]=0;
					for(int k=i+1;k<=n+1;k++)
					a[j][k]-=tmp*a[i][k];
				}
		}
		for(int i=1;i<=n;i++)
			ans[i]=a[i][n+1]/a[i][i];
			return true;
	}
}GE;
int id;
LD ep[N][N][N<<1];
LD fc[N<<1][N<<1];
class TorusSailing {
	public:
	LD expectedTime(int n,int m,int gx,int gy) {
		int sx=n-1-gx,sy=m-1-gy;
		for(int j=m-2;j>=0;j--)
			ep[n-1][j][++id]=1;
		for(int i=n-2;i>=0;i--)
			ep[i][m-1][++id]=1;
		for(int i=n-2;i>=0;i--)
			for(int j=m-2;j>=0;j--) {
				for(int k=1;k<=id+1;k++)
					ep[i][j][k]=(ep[i+1][j][k]+ep[i][j+1][k])/2;
				ep[i][j][id+1]++;
			}
		//ep[n-1][m-1]=(ep[n-1][0]+ep[0][m-1])/2+1
		for(int j=1;j<m;j++) {
			for(int k=1;k<=id+1;k++)
				fc[j][k]=ep[n-1][m-1-j][k]-(ep[n-1][m-1-j+1][k]+ep[0][m-1-j][k])/2;
			fc[j][id+1]=1-fc[j][id+1];
		}
		for(int i=1;i<n;i++) {
			for(int k=1;k<=id+1;k++)
				fc[m-1+i][k]=ep[n-1-i][m-1][k]-(ep[n-1-i][0][k]+ep[n-1-i+1][m-1][k])/2;
			fc[m-1+i][id+1]=1-fc[m-1+i][id+1];
		}
		GE.Init(id,&fc[0][0]);
		GE.Guass();
		LD ans=0;
		for(int k=1;k<=id;k++)
			ans+=GE.ans[k]*ep[sx][sy][k];
		ans+=ep[sx][sy][id+1];
		cout<<ans<<endl;
		return ans;
	}
};

总结

无向图or有向图带环的随机游走：思考转移方向+高斯消元
带状矩阵可用于优化网格图（实际上我觉得可以在编号上做文章）
主元法的本质即把消元的过程手动做了，从而使元的个数减少，往往是根据于比较特殊的转移（不会出现互相有关），最好是每一个位置都可以直接用元来表示，才使用主元法。