計算任意一個圖生成樹的個數——Kirchhoff 的Matrix Tree 方法Java實現
阿新 • • 發佈:2019-02-08
計算任意一個圖的生成樹的個數,是Kirchhoff提出的理論,通常稱為Matrix Tree Theorem,原理很簡單:
Let G be a graph with V(G)={v1,v2,...,vn},let A={aij}be the adjacentcy matrix of G,and let C={cij}be the n*n matrix, where cij=deg vi if i=j; cij=-aij if i!=j; Then the number of spanning trees of G is the vlaue of any cofactor(餘子式) of C
但是需要計算行列式的值,這裡需要點小技巧,所以這裡選擇了LUP分解法來計算。
package trees; import matrix.DeterminantCalculator; /** * 計算任意一個圖的生成樹的個數,是Kirchhoff提出的理論,通常稱為Matrix Tree Theorem * Let G be a graph with V(G)={v1,v2,...,vn},let A={aij}be the adjacentcy matrix of G, * and let C={cij}be the n*n matrix, where cij=deg vi if i=j; cij=-aij if i!=j; Then the * number of spanning trees of G is the vlaue of any cofactor(餘子式) of C * @author xhw * */ public class NumberOfSpanningTree { /** * @param args */ public static void main(String[] args) { double a[][]={{0,1,1,0}, {1,0,1,0}, {1,1,0,1}, {0,0,1,0}}; int n=numberOfSpanningTree(a); System.out.println("numberOfSpanningTree:"+n); } public static int numberOfSpanningTree(double[][] a) { double c[][]=generateKirchhoffMatrix(a); double confactor[][]=new double[c.length-1][c.length-1]; for(int i=1;i<c.length;i++) { for(int j=1;j<c.length;j++) { confactor[i-1][j-1]=c[i][j]; //System.out.print(c[i][j]+" "); } //System.out.println(); } DeterminantCalculator dc=new DeterminantCalculator(); int n=(int)dc.det(confactor); return n; } /** * C={cij}be the n*n matrix, where cij=deg vi if i=j; cij=-aij if i!=j * @param a * @return */ public static double[][] generateKirchhoffMatrix(double[][] a) { int length=a.length; double c[][]=new double[length][length]; for(int i=0;i<length;i++) { int deg=0; for(int j=0;j<length;j++) { deg+=a[i][j]; c[i][j]=-a[i][j]; } c[i][i]=deg; } return c; } } package matrix; /** * 矩陣和可以用來快速地計算矩陣的行列式,因為det(A) = det(L) det(U),而三角矩陣的行列式就是對角線元素的乘積。如果要求L 是單位三角矩陣,Uii的乘積 * 那麼同樣的方法也可以應用於LUP分解,只需乘上P的行列式,即相應置換的符號差。 * @author xhw * */ public class DeterminantCalculator { private LUPDecomposition lu; /** * @param args */ public static void main(String[] args) { // TODO Auto-generated method stub } public DeterminantCalculator() { this.lu=new LUPDecomposition(); } public double det(double a[][]) { a=lu.decomposition(a); if(a==null) return 0; double d=1; for(int i=0;i<a.length;i++) { d=d*a[i][i]; } int n=lu.getExchangeTimes(); if(n%2==0) return d; else return -d; } } package matrix; /** * LUP分解 * P是置換矩陣,L是下三角矩陣,並且對角線為1,U是上三角矩陣;P的出現主要是為了選主元,在選取Ade非對角線元素中選主元避免除數為0, * 除此以外,除數的值也不能過小,否則導致計算中數值不穩定,因此所選的主元是個較大的值。詳細參見演算法導論P461 * @author xhw * */ public class LUPDecomposition{ private int p[]; private int exchangeTimes=0; /** * @param args */ public static void main(String[] args) { // TODO Auto-generated method stub } public double[][] decomposition(double a[][]) { int length=a.length; //p是置換矩陣,p[i]=j說明P的第i行第j列為1 p=new int [length]; for(int i=0;i<length;i++) { p[i]=i; } for(int k=0;k<length;k++) { double max=0; int maxK=0; for(int i=k;i<length;i++) { if(Math.abs(a[i][k])>max) { max=Math.abs(a[i][k]); maxK=i; } } if(max==0) { System.out.println("singular matrix"); return null; } if(k!=maxK) { //交換k和maxk行 exchange(p,k,maxK); exchangeTimes++; for(int i=0;i<length;i++) { double temp=a[k][i]; a[k][i]=a[maxK][i]; a[maxK][i]=temp; } } //“原地”計算LU,矩陣a的上半部分為U,下半部分為L for(int i=k+1;i<length;i++) { a[i][k]=a[i][k]/a[k][k]; for(int j=k+1;j<length;j++) { a[i][j]=a[i][j]-a[i][k]*a[k][j]; } } } return a; } public void exchange(int p[],int k,int maxK) { int temp=p[k]; p[k]=p[maxK]; p[maxK]=temp; } public int[] getP() { return p; } public void setP(int[] p) { this.p = p; } public int getExchangeTimes() { return exchangeTimes; } public void setExchangeTimes(int exchangeTimes) { this.exchangeTimes = exchangeTimes; } }