Implementation of the Jacobian algorithm of the Matrix's eigenvalues and eigenvectors

The eigenvalues and eigenvectors of matrices are a very important concept in linear algebra and matrix theory. It is also frequently used in the field of remote sensing. For example, the principal component analysis of multispectral and hyperspectral images requires the solution of the covariance matrix between the bands or the eigenvalues and eigenvectors of the correlation coefficient matrix.

According to the concept of ordinary linear algebra, eigenvalues and eigenvectors can be obtained by traditional methods, but the actual project is usually calculated by numerical analysis method, and the Jacobian iterative method is introduced to solve the eigenvalues and eigenvectors.

The Jacobian method is used for all eigenvalues and eigenvectors of a realistic symmetric array.

For real symmetric matrix a, there must be orthogonal matrix U. Make

U TA U = D.

The D is a diagonal array, and the main diagonal element Li is a characteristic value. The J column of the orthogonal array U is the characteristic vector of a which belongs to Li .

principle : The Jacobi method uses plane rotation to make a similar transformation to matrix a , and a is a diagonal array, and then the eigenvalues and eigenvectors are obtained.

Now that the rotation is used, here is a description of the rotation matrix.

For P ≠ Q, the following defined N -order matrix UPQ is a planar rotation matrix.

Easy Verify that the upq is an orthogonal array.

For vector x,upq x corresponds to the rotation angle J of the axis Oxp and oxq in the plane in which it is located.

Transformation process: In order to ensure similar conditions, so that the main diagonal outside the element to zero!

Remember n -order square matrix a = [AIJ], make the following transformations to A :

A1= UpqTAUPQ,

a1 is still a real symmetric array, because,upqt =upq-1, it is known that a1 is the same as the eigenvalues of a .

It says earlier that Jacobian is an iterative algorithm. So when each iteration of the step requires a new matrix to be rotated, then the new matrix element is asked, here is a detailed formula such as the following:

It can be seen by a set of formulas above:

(1) The p line, column and line Q of the Matrix a1, the elements in the column change, the other rows, the elements in the column are not changed.

(2) p, q each is the line number of the absolute maximum element on the non-main diagonal of the previous iteration matrix A

(3) J is the angle of rotation. Can be calculated by the following formula:

The detailed process of solving matrix eigenvalues and eigenvectors with Jacobian iterative method is summarized as follows:

(1) the initialization characteristic vector is diagonal array v. That is, the elements of the main diagonal are all 1. Other elements are 0.

(2) in the non-main diagonal element of a , find the absolute maximum element APQ.

(3) using formula (3.14) to calculate tan2J, to find Cosj, Sinj and Matrix upq.

(4) Use the formula (1)-(4) to find A1, and the current eigenvector v by multiplying the current eigenvector matrix V by the matrix UPQ .

(5) if the maximum value of the non-main diagonal element of matrix A before the current iteration is less than the given threshold E. Stop the calculation; otherwise, make a = a1, run repeatedly (2) ~ (5). When the calculation is stopped. Get eigenvalues li≈ (A1) ij , i,j=,...,N. and the Eigenvector v.

(6) This step is optional.

The eigenvalues and eigenvectors of the Matrix are arranged again according to the size of the eigenvalues from large to small.

So far, every step of the calculation process is very clear, writing code is not difficult, detailed code such as the following:

/*** @brief The eigenvalues of the realistic symmetry matrix and the Jacobian method of the eigenvector * Use the vertices (Jacobi) method to seek all the eigenvalues and eigenvectors of the symmetric matrix * @param pmatrix an array of n*n length. Holds the real symmetric matrix * @param the order of the Ndim matrix * @param an array of length n*n, returns the eigenvector (stored by column) * @param dbeps Precision requirements * @param NJT integer variable. Controls the maximum number of iterations * @param pdbeigenvalues array of eigenvalues * @return */bool Cpcaalg::jacbicor (double * Pmatrix,int nDim, double *pdblvects, Doub  Le *pdbeigenvalues, double dbeps,int nJt) {for (int i = 0; i < NDim; i + +) {Pdblvects[i*ndim+i] = 1.0f; for (int j = 0; J < NDim; J + +) {if (i! = j) pdblvects[i*ndim+j]=0.0f;}} int ncount = 0;//Iteration number while (1) {//finds the largest element on a non-diagonal pMatrix, double Dbmax = pmatrix[1];int Nrow = 0;int Ncol = 1;for (int i = 0; i < NDim; i + +)//line {for (int j = 0; J < NDim; J + +)//column {double d = fabs (Pmatrix[i*ndim+j]); if ((I!=j) && (d> Dbmax)) {D   Bmax = D;   Nrow = i; Ncol = j;  }}}if (Dbmax < dbeps)//precision meets requirements break; if (ncount > NJt)//number of iterations exceeds the limit break;ncount++;d ouble dbapp = pmatrix[nrow*ndim+nrow];d ouble DBAPQ = Pmatrix[nrow*ndim +ncol];d ouble dbaqq = pmatrix[ncol*ndim+ncol];//Calculate rotation angle Double dbangle = 0.5*atan2 ( -2*dbapq,dbaqq-dbapp);d ouble Dbsintheta = sin (dbangle);d ouble Dbcostheta = cos (dbangle);d ouble Dbsin2theta = sin (2*dbangle);d ouble dbcos2theta = cos (2*dbangle);p Matrix[nrow*ndim+nrow] = dbApp* Dbcostheta*dbcostheta + Dbaqq*dbsintheta*dbsintheta + 2*dbapq*dbcostheta*dbsintheta;pmatrix[ncol*ndim+ncol] = dbApp* Dbsintheta*dbsintheta + Dbaqq*dbcostheta*dbcostheta-2*dbapq*dbcostheta*dbsintheta;pmatrix[nrow*ndim+ncol] = 0.5* ( Dbaqq-dbapp) *dbsin2theta + dbapq*dbcos2theta;pmatrix[ncol*ndim+nrow] = pmatrix[nrow*ndim+ncol];for (int i = 0; i < NDim; i + +) {if ((I!=ncol) && (I!=nrow)) {int u = i*ndim + nrow;//p int w = i*ndim + Ncol;//qdbmax = pmatrix[u]; pmatr ix[u]= Pmatrix[w]*dbsintheta + Dbmax*dbcostheta; Pmatrix[w]= Pmatrix[w]*dbcostheta-dbmax*dbsintheta; }} for (int j = 0; J < NDim; J + +) {if ((J!=ncol) && (J!=nrow)) {int u = nrow*ndim + J;//pint w = Ncol*ndim + j ;//qdbmax = Pmatrix[u]; pmatrix[u]= Pmatrix[w]*dbsintheta + DbMAx*dbcostheta; Pmatrix[w]= Pmatrix[w]*dbcostheta-dbmax*dbsintheta; The}//computes the eigenvectors for (int i = 0; i < NDim; i + +) {int u = i*ndim + nrow;//p int w = i*ndim + Ncol;//qdbmax = Pdblvects[u]; Pdblvects[u] = Pdblvects[w]*dbsintheta + Dbmax*dbcostheta; PDBLVECTS[W] = Pdblvects[w]*dbcostheta-dbmax*dbsintheta; }//to sort the eigenvalues and to arrange the eigenvectors again, the eigenvalues are pmatrix the elements on the main diagonal std::map<double,int> mapeigen;for (int i = 0; i < NDim; i + +) {pdbe Igenvalues[i] = Pmatrix[i*ndim+i];mapeigen.insert (Make_pair (Pdbeigenvalues[i],i));} Double *pdbtmpvec = new Double[ndim*ndim];std::map<double,int>::reverse_iterator iter = Mapeigen.rbegin (); for ( int j = 0; Iter! = Mapeigen.rend (), J < NDim; ++ITER,++J) {for (int i = 0; i < NDim; i + +) {Pdbtmpvec[i*ndim+j] = Pdblvects[i*ndim + Iter->second];} Eigenvalues are arranged again pdbeigenvalues[j] = Iter->first;} Set positive sign for (int i = 0; i < NDim; i + +) {Double Dsumvec = 0;for (int j = 0; J < NDim; J + +) Dsumvec + = Pdbtmpvec[j * NDi M + i];if (dsumvec<0) {for (int j = 0;j < NDim; J + +) Pdbtmpvec[j * nDim + i] *=-1;}} memcpy (pdblvects,pdbtmpvec,sizeof (double) *ndim*ndim);d elete []pdbtmpvec;return 1;}

Implementation of the Jacobian algorithm of the Matrix's eigenvalues and eigenvectors