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. In the field of remote sensing, for example, the principal component analysis of multispectral and hyperspectral images requires the solution of the covariance matrix of the band 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 the real symmetric matrix a, there must be orthogonal matrix U, so that

U TA U = D.

Where D is a diagonal array, and its 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 similarity 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 n -order matrix UPQ defined below is a planar rotation matrix.

It is easy to verify that 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: Under the guarantee of similar conditions, so that the main diagonal outside the element to zero!

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

A1= UpqTAUPQ,

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

It says that Jacobian is an iterative algorithm, so each iteration requires a new matrix after rotation, so the new matrix elements, here are the specific formula:

You can see it 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 is the number of the largest element on the non-main diagonal of the previous iteration matrix A, respectively.

(3) J is the rotation angle, which can be calculated by the following formula:

The specific steps for solving the eigenvalues and eigenvectors of matrices by using the Jacobian iterative method are summarized as follows:

(1) the initialization feature vector is diagonal array V, that is, the main diagonal element is 1. The other element is 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 in the non-main diagonal element of matrix A before the current iteration is less than the given threshold E, the calculation is stopped; otherwise, a = a1 is repeated (2) ~ (5). When you stop the calculation, you get the eigenvalues Li≈ (A1) ij , i,j=,...,N. and the Eigenvector v.

(6) This step is optional. the eigenvalues and eigenvectors of the matrix are rearranged from large to small according to the size of the characteristic values.

So far, every step of the calculation process is very clear, write code is not difficult, the specific code is as follows:

/*** @brief The eigenvalues of the realistic symmetry matrix and the Jacobian method of the eigenvector * using vertices (Jacobi) method to seek the full eigenvalues and eigenvectors of the symmetric matrix * @param pmatrix the length of the n*n array, storing the real symmetric matrix * @param the order of the Ndim matrix * @param pdblvects length n*n array, return eigenvectors (stored by column) * @param dbeps Precision requirements * @param NJT integer variable, control maximum iterations * @param pdbeigenvalues eigenvalues Array * @retur n */bool Cpcaalg::jacbicor (double * Pmatrix,int nDim, double *pdblvects, double *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 = 0;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 rearrange the eigenvectors, the eigenvalues are pmatrix the elements on the main diagonal std::map<double,int> mapeigen;for (int i = 0; i < NDim; i + +) {Pdbei Genvalues[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];} Eigenvalue rearrangement 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