矩阵的特征值和特征向量是线性代数以及矩阵论中非常重要的一个概念。在遥感领域也是经常用到,比如多光谱以及高光谱图像的主成分分析要求解波段间协方差矩阵或者相关系数矩阵的特征值和特征向量。
根据普通线性代数中的概念,特征值和特征向量可以用传统的方法求得,但是实际项目中一般都是用数值分析的方法来计算,这里介绍一下雅可比迭代法求解特征值和特征向量。
雅克比方法用于求实对称阵的全部特征值、特征向量。
对于实对称阵 A,必有正交阵 U,使
U TA U = D。
其中 D 是对角阵,其主对角线元 li 是 A 的特征值. 正交阵 U 的第 j 列是 A 的属于 li 的特征向量。
原理:Jacobi 方法用平面旋转对矩阵 A 做相似变换,化A 为对角阵,进而求出特征值与特征向量。
既然用到了旋转,这里就介绍一下旋转矩阵。
对于 p ≠ q,下面定义的 n 阶矩阵Upq 是平面旋转矩阵。
![](https://img-blog.csdn.net/20141018160835390?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvemhvdXh1Z3VhbmcyMzY=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
容易验证 Upq是正交阵。对于向量x,Upq x 相当于把坐标轴Oxp和 Oxq 于所在的平面内旋转角度 j .
变换过程: 在保证相似条件下,使主对角线外元素趋于零!
记 n 阶方阵A = [aij], 对 A 做下面的变换:
A1= UpqTAUpq,
A1 仍然是实对称阵,因为,UpqT =Upq-1,知A1与 A 的特征值相同。
前面说了雅可比是一种迭代算法,所以每一步迭代时,需要求出旋转后新的矩阵,那么新的矩阵元素如何求,这里给出具体公式如下:
![](https://img-blog.csdn.net/20141018160855156?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvemhvdXh1Z3VhbmcyMzY=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
由上面的一组公式可以看到:
(1)矩阵A1 的第p 行、列与第 q 行、列中的元素发生了变化,其它行、列中的元素不变。
(2)p、q分别是前一次的迭代矩阵A的非主对角线上绝对值最大元素的行列号
(3) j是旋转角度,可以由下面的公式计算:
![](https://img-blog.csdn.net/20141018160907875?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvemhvdXh1Z3VhbmcyMzY=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
归纳可以得到雅可比迭代法求解矩阵特征值和特征向量的具体步骤如下:
(1) 初始化特征向量为对角阵V,即主对角线的元素都是1.其他元素为0。
(2) 在A的非主对角线元素中,找到绝对值最大元素 apq 。
(3) 用式(3.14)计算tan2j,求 cosj, sinj 及矩阵Upq .
(4) 用公式(1)-(4)求A1;用当前特征向量矩阵V乘以矩阵Upq得到当前的特征向量V。
(5) 若当前迭代前的矩阵A的非主对角线元素中最大值小于给定的阈值e时,停止计算;否则, 令A = A1 , 重复执行(2) ~ (5)。 停止计算时,得到特征值li≈(A1) ij ,i,j= 1,2,…,n.以及特征向量V。
(6) 这一步可选。根据特征值的大小从大到小的顺序重新排列矩阵的特征值和特征向量。
到现在为止,每一步的计算过程都十分清楚了,写出代码也就不是难事了,具体代码如下:
-
-
-
-
-
-
-
-
-
-
-
- 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;
- while(1)
- {
-
- double dbMax = pMatrix[1];
- int nRow = 0;
- int nCol = 1;
- for (int i = 0; i < nDim; i ++)
- {
- for (int j = 0; j < nDim; j ++)
- {
- double d = fabs(pMatrix[i*nDim+j]);
-
- if((i!=j) && (d> dbMax))
- {
- dbMax = d;
- nRow = i;
- nCol = j;
- }
- }
- }
-
- if(dbMax < dbEps)
- break;
-
- if(nCount > nJt)
- break;
-
- nCount++;
-
- double dbApp = pMatrix[nRow*nDim+nRow];
- double dbApq = pMatrix[nRow*nDim+nCol];
- double dbAqq = pMatrix[nCol*nDim+nCol];
-
-
- double dbAngle = 0.5*atan2(-2*dbApq,dbAqq-dbApp);
- double dbSinTheta = sin(dbAngle);
- double dbCosTheta = cos(dbAngle);
- double dbSin2Theta = sin(2*dbAngle);
- double dbCos2Theta = cos(2*dbAngle);
-
- pMatrix[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;
- int w = i*nDim + nCol;
- dbMax = pMatrix[u];
- pMatrix[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;
- int w = nCol*nDim + j;
- dbMax = pMatrix[u];
- pMatrix[u]= pMatrix[w]*dbSinTheta + dbMax*dbCosTheta;
- pMatrix[w]= pMatrix[w]*dbCosTheta - dbMax*dbSinTheta;
- }
- }
-
-
- for(int i = 0; i < nDim; i ++)
- {
- int u = i*nDim + nRow;
- int w = i*nDim + nCol;
- dbMax = pdblVects[u];
- pdblVects[u] = pdblVects[w]*dbSinTheta + dbMax*dbCosTheta;
- pdblVects[w] = pdblVects[w]*dbCosTheta - dbMax*dbSinTheta;
- }
-
- }
-
-
- std::map<double,int> mapEigen;
- for(int i = 0; i < nDim; i ++)
- {
- pdbEigenValues[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];
- }
-
-
- pdbEigenValues[j] = iter->first;
- }
-
-
- for(int i = 0; i < nDim; i ++)
- {
- double dSumVec = 0;
- for(int j = 0; j < nDim; j ++)
- dSumVec += pdbTmpVec[j * nDim + i];
- if(dSumVec<0)
- {
- for(int j = 0;j < nDim; j ++)
- pdbTmpVec[j * nDim + i] *= -1;
- }
- }
-
- memcpy(pdblVects,pdbTmpVec,sizeof(double)*nDim*nDim);
- delete []pdbTmpVec;
-
- return 1;
- }