eigen.java,v

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

head	1.4;
access;
symbols;
locks; strict;
comment	@# @;


1.4
date	2005.05.23.22.34.06;	author rirwin;	state Exp;
branches;
next	1.3;

1.3
date	2005.03.11.05.03.17;	author patil;	state Exp;
branches;
next	1.2;

1.2
date	2005.01.20.02.41.42;	author patil;	state Exp;
branches;
next	1.1;

1.1
date	2004.12.28.00.04.32;	author patil;	state Exp;
branches;
next	;


desc
@@


1.4
log
@Fixed Javadoc errors.
@
text
@/**
 * file: Eigen.java
 *
 * last edited: Ryan Irwin
 *
 */

// import necessary java libraries
//

//These imports are not needed - Phil T. 6-23-03
//import java.awt.*;
//import java.applet.*;
//import javax.swing.*;
//import java.lang.Math.*;

/**
 * this class takes in square matrix and computes the eigen values and
 * its corresponding eigen vectors for the matrix
 *
 */
class Eigen 
{
    
    // *********************************************************************
    //
    // declare global variables and components
    //
    // *********************************************************************
    
    // declare constants
    //
    public static final int MAXRAND = (int)2147483647;
    public static final double INFINITY = (double)1000000000000.0;
    public static final double MIN_INITIAL_GROWTH = (double)1000.0;
    public static final int MAX_INITIAL_ITER = (int)10;
    public static final double EPSILON = (double)0.0000001;
    public static final int MAX_ITER = (int)8;
    public static final int MAX_L_UPDATE = (int)5;
    public static final double MIN_DECREASE = (double)0.01;
    public static final double TINY = (double)1.0e-20;
    public static final int MY_MAXINT = (int)32000;
    public static final int RAND_INIT = (int)27;
    public static final int RAND_BIG = (int)1000000000;
    public static final int RAND_SEED = (int)161803398;
    public static final int RAND_ARRAY = (int)56;
    public static final int RAND_CONST = (int)31;
    public static final double RAND_SCALE = (double)(1.0 / RAND_BIG);
    
    // declare global variables
    //
    public static int rand_init = RAND_INIT;
    public static int temp_init = (int)0;
    public static int m_array[] = new int [RAND_ARRAY];
    public static int next_i = (int)0;
    public static int next_ip = (int)0;

    // *********************************************************************
    //
    // declare class methods
    //
    // *********************************************************************
    
    /**
     * given a matrix a[1..n][1..n], this routine replaces it by the
     * LU decomposition of a row wise permutation of itself. this
     * routine is combined with lubksb to solve linear equations to
     * invert a matrix. see numerical recipes, the art of scientific
     * computing, second edition, cambridge university press. pages 46 - 47.
     *
     * @@param   n dimensions of the input square matrix
     * @@param   indx row permutations effected by partial pivoting
     * @@param   a input matrix
     * @@param   d +/- 1, if row interchanges was even or 
     *          odd respectively
     *
     */
    public static void ludcmp(int n, int indx[], double a[][], Double d) 
    {
	
	int i = 0;
	int imax = 0;
	int j = 0;
	int k = 0;
	double big = 0.0;
	double dum = 0.0;
	double sum = 0.0;
	double temp = 0.0;
	double vv[] = null;
	
	vv = new double[n];
	
	d = new Double(1.0);
	for (i = 0; i < n; i++) 
	{
	    big = 0.0;
	    for (j = 0; j < n; j++)
		if ((temp = Math.abs(a[i][j])) > big) 
		    big = temp;
	    if (big == 0.0) { 
		//System.out.println("Singular matrix in routine LUDCMP");
		return;
	    }
	    vv[i] = 1.0 / big;
	}
	for (j = 0; j < n; j++) 
	{
	    for (i = 0; i < j; i++) 
	    {
		sum = a[i][j];
		for (k = 0; k < i; k++) 
		    sum -= a[i][k] * a[k][j];
		a[i][j] = sum;
	    }
	    big = 0.0;
	    for (i = j; i < n; i++) 
	    {
		sum = a[i][j];
		for (k = 0; k < j; k++)
		    sum -= a[i][k]*a[k][j];
		a[i][j] = sum;
		if ( (dum = vv[i] * Math.abs(sum)) >= big) 
		{
		    big = dum;
		    imax = i;
		}
	    }
	    if (j != imax) 
	    {
		for (k = 0; k < n; k++) 
		{
		    dum = a[imax][k];
		    a[imax][k] = a[j][k];
		    a[j][k] = dum;
		}
		d = new Double(-d.doubleValue());
		vv[imax] = vv[j];
	    }
	    indx[j] = imax;
	    if (a[j][j] == 0.0) 
		a[j][j]=TINY;
	    if (j != n-1) 
	    {
		dum = 1.0 / (a[j][j]);
		for (i = j + 1; i < n; i++) 
		    a[i][j] *= dum;
	    }
	}
    }

    /**
     * given a matrix a [1..n][1..n], this routine solves the equation
     * A.X = B.  see numerical recipes, the art of scientific
     * computing, second edition, cambridge university press. page 47.
     *
     * @@param   a input matrix
     * @@param   b right hand side vector B of A.X = B
     * @@param   n dimensions of the input square matrix
     * @@param   indx row permutations effected by partial pivoting
     *
     */
    static void lubksb(double a[][], double b[], int n, int indx[]) 
    {
	
	int i= 0;
	int ii = 0;
	int ip = 0;
	int j = 0;
	double sum = 0.0;
	
	for (i = 0; i < n; i++) 
	{
	    ip = indx[i];
	    sum = b[ip];
	    b[ip] = b[i];
	    if (ii != 0)
		for (j = ii; j <= i - 1; j++) 
		    sum -= a[i][j]*b[j];
	    else 
		if (sum != 0) 
		    ii = i;
	    b[i] = sum;
	}
	for (i = n - 1; i >= 0; i--) 
	{
	    sum = b[i];
	    for (j = i + 1; j < n; j++) 
		sum -= a[i][j]*b[j];
	    b[i] = sum / a[i][i];
	}
    }

    /**
     * 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
     *
     * @@param   Matrix M: input matrix
     * @@param   double eigVal: eigen value
     * @@param   double[] y: random vector
     * @@param   double[] x: random vector
     *
     */
    public static void linearSolve(Matrix M, double eigVal, 
				   double y[], double x[]) 
    {
	
	Double d = null;
	double a[][] = null;
	double bVec[] = null;

	int n = 0;
	int i = 0;
	int j = 0;
	int indx[] = null;
	
	n = M.row;
	
	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( i = 0; i < n; i++ )
	    a[i][i] -= eigVal;
	
	bVec = new double[n];
	indx = new int[n];
	copy_FV( x, bVec, n );
	
	ludcmp(n, indx, a, d);
	lubksb(a, bVec, n, indx);
	
	copy_FV( bVec, y, n );
    }
   
    /**
     * an Eigenvalue 'eigVal' and the corresponding Eigenvectors this method
     * normalizes the eigen vectors i.e. divides each Eigenvector with the 
     * square root of the corresponding Eigenvalue
     *
     * @@param   eigVal eigen value
     * @@param   eigVec eigen vectors
     *
     */
    public static void normEigVec(double eigVal, double eigVec[]) 
    {
	
	// loop through each eigenvector normalizing them
	//
	for(int i = 0; i < eigVec.length; i++) 
	{
	    eigVec[i] = eigVec[i] / Math.sqrt(eigVal);
	}
    }

    /**
     * given a matrix 'M', an Eigenvalue 'eigVal', calculate the
     * the Eigenvector 'eigVec'
     *
     * @@param   M input matrix
     * @@param   eigVal eigen value
     * @@param   eigVec eigen vectors
     *
     */
    public static void calcEigVec(Matrix M, double eigVal, double eigVec[]) 
    {
	
	double lambdaI = 0;
	double lambdaI1 = 0;
	double xi[] = null;
	double xi1[] = null;
	double y[] = null;
	double y_init_max[] = null;
	double deltaL = 0.0;
	double di = 0.0;
	double di1 = 0.0;
	double len_y = 0.0;
	double sum_nom = 0.0;
	double sum_denom = 0.0;
	double maxGF = 0.0;
	double decrease = 0.0;
	
	int n = 0;
	int i = 0;
	int nrIter = 0;
	int nrInitialIter = 0;
	int nrLambdaUpdates = 0;
	int done = 0;
	int betterEigVal = 0;
	
	n = M.row;
	
	xi = new double[n];
	xi1 = new double[n];
	y = new double[n];
	y_init_max = new double[n];
	
	lambdaI = eigVal;
	
	maxGF = -INFINITY;
	do 
	{
	    nrInitialIter++;
	    random_unit_vector( xi, n );
	    linearSolve( M, lambdaI, y, xi );
	    len_y = vector_len( y, n );
	    
	    if( len_y > maxGF ) 
	    {
		maxGF = len_y;
		copy_FV( y, y_init_max, n );
	    }
	}
	while( len_y < MIN_INITIAL_GROWTH && 
	       nrInitialIter < MAX_INITIAL_ITER );
	copy_FV( y_init_max, y, n );
	
	nrIter = 0;
	di = 0.0;
	for( i = 0; i < n; i++ ) 
	{
	    xi1[i] = y[i] / len_y;
	}
	
	do {
	    nrIter++;
	    di1 = Euclidian_Distance( xi1, xi, n );
	    
	    if( di1 < EPSILON || nrIter >= MAX_ITER || 
		nrLambdaUpdates >= MAX_L_UPDATE )
		done = 1;
	    else 
	    {
		decrease = (double) Math.abs( di1 - di);
		if( decrease < MIN_DECREASE ) 
		{
		    nrLambdaUpdates++;
		    sum_nom = 0.0; 
		    sum_denom = 0.0;
		    for( i = 0; i < n; i++ ) 
		    {
			sum_nom += xi[i] * xi[i];
			sum_denom += xi[i] * y[i];
		    }
		    deltaL = sum_nom / (double)Math.abs(sum_denom);
		    lambdaI1 = lambdaI + deltaL;
		    betterEigVal = 1;
		    lambdaI = lambdaI1;
		    if( deltaL == 0.0 ) 
		    {
			done = 1;
		    }
		    else 
		    {
			nrIter = 0;
		    }
		}
		if( done == 0 ) 
		{
		    copy_FV( xi1, xi, n );
		    di = di1;
		    linearSolve( M, lambdaI, y, xi );
		    len_y = vector_len( y, n );
		    for( i = 0; i < n; i++ ) 
		    {
			xi1[i] = y[i] / len_y;
		    }
		}
	    }
	}
	while( done == 0 );
	
	copy_FV( xi1, eigVec, n );
	if( betterEigVal == 1 ) 
	{
	    eigVal = lambdaI;
	}
    }
    
    /**
     * this routine generates a random vector of length dim with unit size
     *
     * @@param   x input vector
     * @@param   dim vector dimension
     *
     */
    public static void random_unit_vector(double x[], int dim) 
    {
	
	int i = 0;
	int rand = 0;
	
	double scale = 0.0;
	
	for( i = 0; i < dim; i++ ) 
	{
	    scale = myrandom();
	    rand = (int)(MY_MAXINT * scale);
	    x[i] = ((double)rand - (double)(MY_MAXINT/2)) 
		/ (double)MY_MAXINT;
	}
	norm_vector( x, dim );
    }
    
    /**
     * this routine uses knuth's subtractive algorithm to generate
     * uniform random numbers in (0, 1)
     *
     * return   random number between 0 and 1
     *
     */
    public static double myrandom() 
    {
	
	int mj = (int)0;
	int mk = (int)0;
	int ii = (int)0;
	int i = (int)0;
	int k = (int)0;
	
	if ((rand_init < 0) || (temp_init == 0)) 
	{
	    
	    temp_init = (int)1;
	    mj = RAND_SEED - (rand_init < 0 ? -rand_init : rand_init);
	    mj %= RAND_BIG;
	    m_array[RAND_ARRAY - 1] = mj;
	    mk = (int)1;
	    
	    ii = (int)0;
	    for (i = 0; i < RAND_ARRAY - 1; i++)
	    {
		ii = (21 * i) % (RAND_ARRAY - 1);
		m_array[ii] = mk;
		mk = mj - mk;
		if (mk < 0) 
		{
		    mk += RAND_BIG;
		}
		mj = m_array[ii];
	    }
	    
	    for (k = 0; k < 5; k ++) 
	    {
		for (i = 0; i < RAND_ARRAY; i++)
		{
		    m_array[i] -= m_array[1 + (i + RAND_CONST - 1) % 
					 (RAND_ARRAY - 1)];
		    if (m_array[i] < 0) 
		    {
			m_array[i] += RAND_BIG;
		    }
		}
	    }
	    
	    next_i = (int)0;
	    next_ip = RAND_CONST;
	    rand_init = (int)1;
	}
	
	if (++next_i == RAND_ARRAY) 
	{
	    next_i = (int)1;
	}
	if (++next_ip == RAND_ARRAY) 
	{
	    next_ip = (int)1;
	}
	mj = m_array[next_i] - m_array[next_ip];
	if (mj < 0) 
	{
	    mj += RAND_BIG;
	}
	m_array[next_i] = mj;
	
	double randnum = (double)mj * RAND_SCALE;
	
	return randnum;
    }

    /**
     * this routine takes in a vector and normalizes the vector
     *
     * @@param   V input vector
     * @@param   dim vector dimension
     *
     */
    public static void norm_vector(double V[], int dim) 
    {

	int i = 0; 
	double len = 0.0;

	len = vector_len( V, dim );
	if( len != 0.0 ) {
	    for( i = 0; i < dim; i++ ) {
		V[i] /= len;
	    }
	}
    }

    /**
     * this routine returns the length of the vector
     *
     * @@param   V input vector
     * @@param   dim vector dimension
     *
     */
     public static double vector_len(double V[], int dim) 
    {
	
	int i = 0;
	double L = 0.0;
	
	for( i = 0; i < dim; i++ ) 
	{
	    L += V[i] * V[i];
	}
	return( (double)Math.sqrt((double)L) );
    }

    /**
     * this routine copies one vectors contents to another
     *
     * @@param   src source vector
     * @@param   dest destination vector
     * @@param   dim vector dimension
     *
     */
    public static void copy_FV(double src[], double dest[], int dim) 
    {
	
	int i = 0;
	
	for( i = 0; i < dim; i++ ) 
	{
	    dest[i] = src[i];
	}
    }
    
    /**
     * this function computes the euclidean distance between two vwctors
     *
     * @@param   S1 first vector
     * @@param   S2 second vector
     * @@param   dim vector dimension
     *
     * @@return Euclidian Distance
     */
     public static double Euclidian_Distance(double S1[], 
					     double S2[], int dim) 
    {
	
	int i = 0;
	double aux = 0.0;
	double D = 0.0;
	
	for( i = 0; i < dim; i++ ) 
	{
	    aux = S1[i] - S2[i];
	    D += aux * aux;
	}
	return( (double)Math.sqrt((double)D) );
    }
    
    /**
     * this routine finds all eignvalues of an upper hessenberg matrix. see 
     * numerical recipes, the art of scientific computing, second edition, 
     * cambridge university press. pages 491 - 493.
     *
     * @@param   M input matrix
     * @@param   wr real part of the eigen vectors
     * @@param   wi immaginary part of the eigen vectors
     *
     */
     public static void hqr(Matrix M, double wr[], double wi[]) {
	
	int n = 0;
	int nn = 0;
	int m = 0;
	int l = 0;
	int k = 0;
	int j = 0;
	int its = 0;
	int i = 0;
	int mmin = 0;

	double z = 0.0;
	double y = 0.0;
	double x = 0.0;
	double w = 0.0;
	double v = 0.0;
	double u = 0.0;
	double t = 0.0;