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* T*A 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* :

*A*1**= ***U**pqT***AU**PQ,

*a*1 is still a real symmetric array, because,*upqt *=*upq*-1, known as *a*1 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 *a*1, 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 tan2*J*, to find Cos*j, *Sin*j* and Matrix *upq.*

**(4)** Use the formula (1)-(4) to find *A*1, 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* = *a*1 is repeated (2) ~ (5). When you stop the calculation, you get the eigenvalues *Li*≈ (*A*1) *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