計算任意一個圖生成樹的個數——Kirchhoff 的Matrix Tree 方法Java實作

計算任意一個圖的生成樹的個數，是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;
	}

}