eigen.java,v

包含了模式识别中常用的一些分类器设计算法

	double s = 0.0;
	double r = 0.0;
	double q = 0.0;
	double p = 0.0; 
	double anorm = 0.0;
	double tmp = 0.0;
	double a[][] = null;

	int hqr_max_iterations = 30;

	n = M.col;
	
	a = new double[n][n];
	
	for( i = 0; i < n; i++ )
	    for( j = 0; j < n; j++ )
		a[i][j] = (double)M.Elem[i][j];
	
	anorm = Math.abs(a[0][0]);
	for (i = 2; i <= n; i++)
	    for (j = (i - 1); j <= n; j++)
		anorm += Math.abs(a[i - 1][j - 1]);
	nn = n;
	t = 0.0;
	while (nn >= 1) 
	{
	    its = 0;
	    do {
		for (l = nn; l >= 2; l--) 
		{
		    s = Math.abs(a[l - 2][l - 2]) + Math.abs(a[l - 1][l - 1]);
		    if (s == 0.0) 
			s = anorm;
		    if ((double)(Math.abs(a[l - 1][l - 2]) + s) == s) 
			break;
		}
		x = a[nn - 1][nn - 1];
		if (l == nn) 
		{
		    wr[nn - 1] = x + t;
		    wi[nn - 1] = 0.0; 
		    nn--;
		} 
		else 
		{
		    y = a[nn - 2][nn - 2];
		    w = a[nn - 1][nn - 2] * a[nn - 2][nn - 1];
		    if (l == (nn-1)) 
		    {
			p = 0.5 * (y - x);
			q = p * p + w;
			z = Math.sqrt((double)Math.abs(q));
			x += t;
			if (q >= 0.0) 
			{
			    if(p > 0) 
			    {
				tmp = Math.abs(z);
			    } 
			    else
			    {
				tmp = -Math.abs(z);
			    }	    
			    z = p + tmp;
			    wr[nn - 2] = wr[nn - 1] = x + z;
			    if (z != 0) 
				wr[nn - 1] = x - w / z;
			    wi[nn - 2] =wi[nn - 1] = 0.0;
			} 
			else 
			{
			    wr[nn - 2] = wr[nn - 1] = x + p;
			    wi[nn - 2] = -(wi[nn - 1] = z);
			}
			nn -= 2;
		    } 
		    else 
		    {
			if (its == hqr_max_iterations) 
			{
			    //System.out.println("Too many iterations in HQR"); 
			    return; 
			}
			if (its == 10 || its == 20) 
			{
			    t += x;
			    for (i = 1; i <= nn; i++) 
				a[i - 1][i - 1] -= x;
			    s = Math.abs(a[nn - 1][nn - 2]) + Math.abs(a[nn - 2][nn - 3]);
			    y = x = 0.75 * s;
			    w = -0.4375 * s * s;
			}
			++its;
			for (m = (nn - 2); m >= l; m--) 
			{
			    z = a[m - 1][m - 1];
			    r = x - z;
			    s = y - z;
			    p = (r * s - w) / a[m][m - 1] + a[m - 1][m];
			    q = a[m][m] - z - r - s;
			    r = a[m + 1][m];
			    s = Math.abs(p) + Math.abs(q) + Math.abs(r);
			    p /= s;
			    q /= s;
			    r /= s;
			    if (m == l) 
				break;
			    u = Math.abs(a[m - 1][m - 2]) * (Math.abs(q) + Math.abs(r));
			    v = Math.abs(p) * (Math.abs(a[m - 2][m - 2]) +
					   Math.abs(z) + Math.abs(a[m][m]));
			    if ((double)(u + v) == v) 
				break;
			}
			for (i = m + 2; i <= nn; i++) 
			{
			    a[i - 1][i - 3] = 0.0;
			    if  (i != (m+2)) a[i - 1][i - 4] = 0.0;
			}
			for (k = m; k <= nn - 1; k++) 
			{
			    if (k != m) 
			    {
				p = a[k - 1][k - 2];
				q = a[k][k - 2];
				r = 0.0;
				if (k != (nn-1)) 
				    r = a[k + 1][k - 2];
				x = Math.abs(p) + Math.abs(q) + Math.abs(r);
				if (x != 0) 
				{
				    p /= x;
				    q /= x;
				    r /= x;
				}
			    }
			    if(p > 0) 
			    {
				tmp = Math.abs(Math.sqrt((double)(p * p + q * q + r * r)));
			    } 
			    else 
			    {
				tmp = -Math.abs(Math.sqrt((double)(p * p + q * q + r * r)));
			    }	  
			    s = tmp;
			    if (s != 0) 
			    {
				if (k == m) 
				{
				    if (l != m)
					a[k - 1][k - 2] = -a[k - 1][k - 2];
				} 
				else
				    a[k - 1][k - 2] = -s * x;
				p += s;
				x = p / s;
				y = q / s;
				z = r / s;
				q /= p;
				r /= p;
				for (j = k; j <= nn; j++) 
				{
				    p = a[k - 1][j - 1] + q * a[k][j - 1];
				    if (k != (nn-1)) 
				    {
					p += r * a[k + 1][j - 1];
					a[k + 1][j - 1] -= p * z;
				    }
				    a[k][j - 1] -= p * y;
				    a[k - 1][j - 1] -= p * x;
				}
				mmin = nn< k + 3 ? nn : k + 3;
				for (i = l; i <= mmin; i++) 
				{
				    p = x * a[i - 1][k - 1] + y * a[i - 1][k];
				    if (k != (nn-1)) 
				    {
					p += z * a[i - 1][k + 1];
					a[i - 1][k + 1] -= p * r;
				    }
				    a[i - 1][k] -= p * q;
				    a[i - 1][k - 1] -= p;
				}
			    }
			}
		    }
		}
	    } while (l < nn - 1);
	}
	
	for( i = 0; i < n; i++ )
	    for( j = 0; j < n; j++ )
		M.Elem[i][j] = (double)a[i][j];
    }

    /**
     * this routine reduces a matrix to hessenberg's form by the elimination
     * method. see numerical recipes, the art of scientific computing, second 
     * edition, cambridge university press. pages 484 - 486.
     *  
     * @@param   M input matrix
     *
     */
    public static void elmhes(Matrix M) 
    {

	int n = 0;
	int m = 0;
	int j = 0;
	int i = 0;

	double y = 0.0;
	double x = 0.0;
	double tmp = 0.0;
	double a[][] = null;

	n = M.col;
	
	a = new double[n][n];
	
	for( i = 0; i < n; i++ )
	    for( j = 0; j < n; j++ )
		a[i][j] = (double)M.Elem[i][j];
	
	for (m = 2; m < n; m++) 
	{
	    x = 0.0;
	    i = m;
	    for (j = m;j <= n; j++) 
	    {
		if (Math.abs(a[j - 1][m - 2]) > Math.abs(x)) 
		{
		    x = a[j - 1][m - 2];
		    i = j;
		}
	    }
	    if (i != m) 
	    {
		for (j = m - 1; j <= n; j++) 
		{
		    tmp = a[i - 1][j - 1]; 
		    a[i - 1][j - 1] = a[m - 1][j - 1]; 
		    a[m - 1][j - 1] = tmp;
		}
		for (j = 1; j <= n; j++) 
		{
		    tmp = a[j - 1][i - 1]; 
		    a[j - 1][i - 1] = a[j - 1][m - 1]; 
		    a[j - 1][m - 1] = tmp;
		}
	    }
	    if (x != 0) 
	    {
		for (i = m + 1; i <= n; i++) 
		{
		    y = a[i - 1][m - 2];
		    if (y != 0) 
		    {
			y /= x;
			a[i - 1][m - 2] = y;
			for (j = m; j <= n; j++)
			    a[i - 1][j - 1] -= y * a[m - 1][j - 1];
			for (j = 1; j <= n; j++)
			    a[j - 1][m - 1] += y * a[j - 1][i - 1];
		    }
		}
	    }
	}
	
	for( i = 0; i < n; i++ )
	    for( j = 0; j < n; j++ )
		M.Elem[i][j] = (double)a[i][j];
    }
    
    /**
     * given a matrix this routine replaces it by a balanced matrix with
     * identical eigenvalues. see numerical recipes, the art of scientific
     * computing, second edition, cambridge university press. pages 482 - 484.
     *
     * @@param M input matrix
     *
     */
    public static void balanc(Matrix M) 
    {
	
	int n = 0;
	int last = 0;
	int j = 0;
	int i = 0;
	
	double s = 0.0;
	double r = 0.0;
	double g = 0.0;
	double f = 0.0;
	double c = 0.0;
	double sqrdx = 0.0;
	double a[][] = null;
	double radix = 2.0;
	
	n = M.col;
	
	a = new double[n][n];
	
	for( i = 0; i < n; i++ )
	    for( j = 0; j < n; j++ )
		a[i][j] = (double)M.Elem[i][j];
	
	sqrdx = radix * radix;
	last = 0;
	while (last == 0) 
	{
	    last = 1;
	    for (i = 1; i <= n; i++) 
	    {
		r = c = 0.0;
		for (j = 1; j <= n; j++)
		    if (j != i) 
		    {
			c += Math.abs(a[j - 1][i - 1]);
			r += Math.abs(a[i - 1][j - 1]);
		    }
		if ((c != 0) && (r != 0)) 
		{
		    g = r / radix;
		    f = 1.0;
		    s = c + r;
		    while (c < g) 
		    {
			f *= radix;
			c *= sqrdx;
		    }
		    g= r * radix;
		    while (c > g) 
		    {
			f /= radix;
			c /= sqrdx;
		    }
		    if ((c + r) / f < 0.95 * s) 
		    {
			last = 0;
			g = 1.0 / f;
			for (j = 1; j <= n; j++) 
			    a[i - 1][j - 1] *= g;
			for (j = 1; j <= n; j++) 
			    a[j - 1][i - 1] *= f;
		    }
		}
	    }
	}
	
	for( i = 0; i < n; i++ )
	    for( j = 0; j < n; j++ )
		M.Elem[i][j] = (double)a[i][j];
    }
    
    /**
     * this method computes the eigen values of a symmetric matrix
     *
     * @@param   T input matrix
     *
     * @@return  eigen values of the input matrix
     *
     */
    public static double[] compEigenVal(Matrix T) 
    {
	
	// declare local variables
	//
	double wr[] = null;
	double wi[] = null;
	
	// allocating memory to store the eigen values
	//
	wr = new double[T.row];
	wi = new double[T.row];
	
	// computing the eigen values
	//
	if( T.row == 1 ) 
	{
	    wr[0] = T.Elem[0][0];
	    wi[0] = 0.0;
	}
	else 
	{
	    balanc( T );
	    elmhes( T );
	    hqr( T, wr, wi );
	}
	
	// return the eigen values (real part only)
	//
	return wr;
    }

    /**    
     * this method sorts the eigen values in decreasing order of importance
     *
     * @@param   wr real eigen values
     *
     * @@return  sorted eigen values in increasing order
     */
    public static double[] sortEigen(double wr[]) 
    {

	// declare local variables
	//
	int i = 0;
	int j = 0;
	int ind = 0;
	double tmp = 0.0;
	double largest = 0.0;

	// sorting the eigen values (real part) in decreasing order
	//
	for(i=0; i < wr.length; i++) 
	{
	    ind = i;
	    largest = wr[i];    
	    for(j=i+1; j < wr.length; j++) 
	    {
		if(wr[j] > largest) 
		{
		    ind = j;
		    largest = wr[j];
		}
	    }
	    if(i != j) 
	    {     
		tmp = wr[i];
		wr[i] = wr[ind];
		wr[ind] = tmp;
	    }
	}
	
	// return the sorted array
	//
	return wr;
    }
}
@


1.3
log
@Comments changed to Java Documentation Style.
@
text
@d4 2
a17 2
 * class: Eigen
 *
a64 10
     * method: ludcmp
     *
     * @@param   int n: dimensions of the input square matrix
     * @@param   double[][] a: input matrix
     * @@param   int[] indx: row permutations effected by partial pivoting
     * @@param   Double d: +/- 1, if row interchanges was even or 
     *          odd respectively
     *
     * @@return  none
     *
d71 6
a151 9
     * method: lubksb
     *
     * @@param   int n: dimensions of the input square matrix
     * @@param   double[][] a: input matrix
     * @@param   int[] indx: row permutations effected by partial pivoting
     * @@param   double[] b: right hand side vector B of A.X = B
     *
     * @@return  none
     *
d156 5
d194 4
a197 1
     * method: linearSolve
d201 1
a202 8
     * @@param   double[] y: random vector
     *
     * @@return  none
     *
     * this routine uses the ludcmp and lubksb routines to find eigenvectors 
     * corresponding to a given eigen value using the inverse iteration method.
     * see numerical recipes, the art of scientific computing, second edition,
     * cambridge university press. pages 493 - 495
a239 7
     * method: calcEigVec
     *
     * @@param   double eigVal: eigen value
     * @@param   double[] eigVec: eigen vectors
     *
     * @@return  none
     *
d244 3
a259 8
     * method: calcEigVec
     *
     * @@param   Matrix M: input matrix
     * @@param   double eigVal: eigen value
     * @@param   double[] eigVec: eigen vectors
     *
     * @@return  none
     *
d263 4
d384 1
a384 4
     * method: random_unit_vector
     *
     * @@param   double[] x: input vector
     * @@param   int dim: vector dimension
d386 2
a387 3
     * @@return  none
     *
     * this routine generates a random vector of length dim with unit size
d409 2
a410 1
     * method: myrandom
a411 1
     * @@param   none
a413 3
     * this routine uses knuth's subtractive algorithm to generate
     * uniform random numbers in (0, 1)
     *
d485 1
a485 1
     * method: norm_vector
d487 2
a488 6
     * @@param   double[] V: input vector
     * @@param   int dim: vector dimension
     *
     * @@return  none
     *
     * this routine takes in a vector and normalizes the vector
d506 1
a506 4
     * method: vector_len
     *
     * @@param   double[] V: input vector
     * @@param   int dim: vector dimension
d508 2
a509 1
     * @@return  none
a510 2
     * this routine returns the length of the vector
     * 
d526 1
a526 5
     * method: copy_FV
     *
     * @@param   double[] src: source vector
     * @@param   double[] dest: destination vector
     * @@param   int dim: vector dimension
d528 3
a530 1
     * @@return  none
a531 2
     * this routine copies one vectors contents to another
     * 
d545 1
a545 1
     * method: Euclidian_Distance
d547 3
a549 7
     * @@param   double[] S1: first vector
     * @@param   double[] S2: second vector
     * @@param   int dim: vector dimension
     *
     * @@return  none
     *
     * this function computes the euclidean distance between two vwctors
d551 1
a569 8
     * method: hqr
     *
     * @@param   Matrix M: input matrix
     * @@param   double[] wr: real part of the eigen vectors
     * @@param   double[] w1: immaginary part of the eigen vectors
     *
     * @@return  none
     *
d574 4
a792 6
     * method: elmhes
     *
     * @@param   Matrix M: input matrix
     *
     * @@return  none
     *
d797 2
a871 7
     * method: balanc
     *
     * @@param 
     *    Matrix M: input matrix
     *
     * @@return  none
     *
d876 2
d953 1
a953 1
     * method: compEigenVal
d955 1
a955 1
     * @@param   Matrix T: input matrix
a958 2
     * this method computes the eigen values of a symmetric matrix
     *
d993 1
a993 1
     * method: sortEigen
d995 1
a995 1
     * @@param   double[] wr: real eigen values
a997 3
     *
     * this method sorts the eigen values in decreasing order of importance
     *
@


1.2
log
@Algorithm Complete and debug statements commented
@
text
@d1 4
a4 2
// file: Eigen.java
//
d15 10