eigen.java,v

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

a24 7
// class: Eigen
//
// this class takes in square matrix and computes the eigen values and
// its corresponding eigen vectors for the matrix
//
class Eigen {

d30 1
a30 1

d64 20
a83 17
    // method: ludcmp
    //
    // arguments: 
    //    int n: dimensions of the input square matrix
    //     double[][] a: input matrix
    //    int[] indx: row permutations effected by partial pivoting
    //    Double d: +/- 1, if row interchanges was even or odd respectively
    //
    // return: none
    //
    // 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 are of scientific
    // computing, second edition, cambridge university press. pages 46 - 47.
    //
    public static void ludcmp(int n, int indx[], double a[][], Double d) {
d98 6
a103 4
	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;
d108 1
a108 1
	    vv[i]=1.0/big;
d110 14
a123 10
	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++)
d125 5
a129 4
		a[i][j]=sum;
		if ( (dum=vv[i]*Math.abs(sum)) >= big) {
		    big=dum;
		    imax=i;
d132 7
a138 5
	    if (j != imax) {
		for (k=0;k<n;k++) {
		    dum=a[imax][k];
		    a[imax][k]=a[j][k];
		    a[j][k]=dum;
d141 1
a141 1
		vv[imax]=vv[j];
d143 8
a150 5
	    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;
d155 17
a171 15
    // method: lubksb
    //
    // arguments: 
    //    int n: dimensions of the input square matrix
    //     double[][] a: input matrix
    //    int[] indx: row permutations effected by partial pivoting
    //    double[] b: right hand side vector B of A.X = B
    //
    // return: none
    //
    // given a matrix a [1..n][1..n], this routine solves the equation
    // A.X = B.  see numerical recipes, the are of scientific
    // computing, second edition, cambridge university press. page 47.
    //
    static void lubksb(double a[][], double b[], int n, int indx[]) {
d179 5
a183 4
	for (i=0;i<n;i++) {
	    ip=indx[i];
	    sum=b[ip];
	    b[ip]=b[i];
d185 13
a197 8
		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];
d200 17
a216 16
    
    // method: linearSolve
    //
    // arguments: 
    //    Matrix M: input matrix
    //    double eigVal: eigen value
    //    double[] x: random vector
    //    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 are of scientific computing, second edition,
    // cambridge university press. pages 493 - 495
    //
d218 2
a219 1
				   double y[], double x[]) {
d245 2
a246 2
	ludcmp(n,indx,a,d);
	lubksb(a,bVec,n,indx);
d251 15
a265 13
    // method: calcEigVec
    //
    // arguments: 
    //    double eigVal: eigen value
    //    double[] eigVec: eigen vectors
    //
    // return: none
    //
    // 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
    //
    public static void normEigVec(double eigVal, double eigVec[]) {
d269 2
a270 1
	for(int i=0; i < eigVec.length; i++) {
d275 15
a289 13
    // method: calcEigVec
    //
    // arguments: 
    //    Matrix M: input matrix
    //    double eigVal: eigen value
    //    double[] eigVec: eigen vectors
    //
    // return: none
    //
    // given a matrix 'M', an Eigenvalue 'eigVal', calculate the
    // the Eigenvector 'eigVec'
    //
    public static void calcEigVec(Matrix M, double eigVal, double eigVec[]) {
d324 2
a325 1
	do {
d331 2
a332 1
	    if( len_y > maxGF ) {
d343 2
a344 1
	for( i = 0; i < n; i++ ) {
d355 2
a356 1
	    else {
d358 2
a359 1
		if( decrease < MIN_DECREASE ) {
d361 4
a364 2
		    sum_nom = 0.0; sum_denom = 0.0;
		    for( i = 0; i < n; i++ ) {
d372 2
a373 1
		    if( deltaL == 0.0 ) {
d376 2
a377 1
		    else {
d381 2
a382 1
		if( done == 0 ) {
d387 2
a388 1
		    for( i = 0; i < n; i++ ) {
d397 2
a398 1
	if( betterEigVal == 1 ) {
d403 13
a415 11
    // method: random_unit_vector
    //
    // arguments: 
    //    double[] x: input vector
    //    int dim: vector dimension
    //
    // return: none
    //
    // this routine generates a random vector of length dim with unit size
    //
    public static void random_unit_vector(double x[], int dim) {
d422 2
a423 1
	for( i = 0; i < dim; i++ ) {
d432 12
a443 9
    // method: myrandom
    //
    // arguments: none
    // return   : random number between 0 and 1
    //
    // this routine uses knuth's subtractive algorithm to generate
    // uniform random numbers in (0, 1)
    //
    public static double myrandom() {
d451 2
a452 1
	if ((rand_init < 0) || (temp_init == 0)) {
d461 2
a462 1
	    for (i = 0; i < RAND_ARRAY - 1; i++) {
d466 2
a467 1
		if (mk < 0) {
d473 4
a476 2
	    for (k = 0; k < 5; k ++) {
		for (i = 0; i < RAND_ARRAY; i++) {
d479 2
a480 1
		    if (m_array[i] < 0) {
d491 2
a492 1
	if (++next_i == RAND_ARRAY) {
d495 2
a496 1
	if (++next_ip == RAND_ARRAY) {
d500 2
a501 1
	if (mj < 0) {
d511 13
a523 11
    // method: norm_vector
    //
    // arguments: 
    //    double[] V: input vector
    //    int dim: vector dimension
    //
    // return: none
    //
    // this routine takes in a vector and normalizes the vector
    //
    public static void norm_vector(double V[], int dim) {
d536 13
a548 11
    // method: vector_len
    //
    // arguments: 
    //    double[] V: input vector
    //    int dim: vector dimension
    //
    // return: none
    //
    // this routine returns the length of the vector
    // 
     public static double vector_len(double V[], int dim) {
d553 2
a554 1
	for( i = 0; i < dim; i++ ) {
d560 14
a573 12
    // method: copy_FV
    //
    // arguments: 
    //    double[] src: source vector
    //    double[] dest: destination vector
    //    int dim: vector dimension
    //
    // return: none
    //
    // this routine copies one vectors contents to another
    // 
    public static void copy_FV(double src[], double dest[], int dim) {
d577 2
a578 1
	for( i = 0; i < dim; i++ ) {
d583 12
a594 11
    // method: Euclidian_Distance
    //
    // arguments: 
    //    double[] S1: first vector
    //    double[] S2: second vector
    //    int dim: vector dimension
    //
    // return: none
    //
    // this function computes the euclidean distance between two vwctors
    //
d596 2
a597 1
					     double S2[], int dim) {
d603 2
a604 1
	for( i = 0; i < dim; i++ ) {
d611 14
a624 13
    // method: hqr
    //
    // arguments: 
    //    Matrix M: input matrix
    //    double[] wr: real part of the eigen vectors
    //    double[] w1: immaginary part of the eigen vectors
    //
    // return: none
    //
    // this routine finds all eignvalues of an upper hessenberg matrix. see 
    // numerical recipes, the are of scientific computing, second edition, 
    // cambridge university press. pages 491 - 493.
    //
d662 9
a670 8
	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;
d672 24
a695 16
		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));
d697 4
a700 2
			if (q >= 0.0) {
			    if(p > 0) {
d702 3
a704 1
			    } else {
d707 10
a716 7
			    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);
d719 5
a723 2
		    } else {
			if (its == hqr_max_iterations) {
d727 2
a728 1
			if (its == 10 || its == 20) {
d730 5
a734 4
			    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;
d737 9
a745 8
			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);
d749 7
a755 5
			    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;
d757 4
a760 3
			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;
d762 12
a773 8
			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) {
d779 7
a785 4
			    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)));
d787 5
a791 3
			    s=tmp;
			    if (s != 0) {
				if (k == m) {
d793 4
a796 3
					a[k-1][k-2] = -a[k-1][k-2];
				} else
				    a[k-1][k-2] = -s*x;
d798 3
a800 3
				x=p/s;
				y=q/s;
				z=r/s;
d803 7
a809 5
				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;
d811 2
a812 2
				    a[k][j-1] -= p*y;
				    a[k-1][j-1] -= p*x;
d814 8
a821 6
				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;
d823 2
a824 2
				    a[i-1][k] -= p*q;
				    a[i-1][k-1] -= p;
d830 1
a830 1
	    } while (l < nn-1);
d838 14
a851 12
    // method: elmhes
    //
    // arguments: 
    //    Matrix M: input matrix
    //
    // return: none
    //
    // this routine reduces a matrix to hessenberg's form by the elimination
    // method. see numerical recipes, the are of scientific computing, second 
    // edition, cambridge university press. pages 484 - 486.
    //  
    public static void elmhes(Matrix M) {
d871 34
a904 24
	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) {
d906 5
a910 5
			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];
d921 15
a935 12
    // method: balanc
    //
    // arguments: 
    //    Matrix M: input matrix
    //
    // return: none
    //
    // given a matrix this routine replaces it by a balanced matrix with
    // identical eigenvalues. see numerical recipes, the are of scientific
    // computing, second edition, cambridge university press. pages 482 - 484.
    //
    public static void balanc(Matrix M) {
d941 1
a941 1

d959 21
a979 16
	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) {
d983 3
a985 2
		    g=r*radix;
		    while (c>g) {
d989 8
a996 5
		    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;
d1007 12
a1018 10
    // method: compEigenVal
    //
    // arguments: 
    //    Matrix T: input matrix
    //
    // return   : eigen values of the input matrix
    //
    // this method computes the eigen values of a symmetric matrix
    //
    public static double[] compEigenVal(Matrix T) {
d1024 1
a1024 1

d1029 1
a1029 1

d1032 2
a1033 1
	if( T.row == 1 ) {
d1037 2
a1038 1
	else {
d1043 1
a1043 1

d1049 12
a1060 10
    // method: sortEigen
    //
    // arguments: 
    //    double[] wr: real eigen values
    //
    // return   : sorted eigen values in increasing order
    //
    // this method sorts the eigen values in decreasing order of importance
    //
    public static double[] sortEigen(double wr[]) {
d1072 2
a1073 1
	for(i=0; i < wr.length; i++) {
d1076 4
a1079 2
	    for(j=i+1; j < wr.length; j++) {
		if(wr[j] > largest) {
d1084 2
a1085 1
	    if(i != j) {     
d1091 1
a1091 1

a1096 3

// 
// end of file
@


1.1
log
@Initial revision
@
text
@d95 1
a95 1
		System.out.println("Singular matrix in routine LUDCMP");
d631 1
a631 1
			    System.out.println("Too many iterations in HQR"); 
@