矩陣的特征值和特征向量是線性代數(shù)以及矩陣論中非常重要的一個(gè)概念。在遙感領(lǐng)域也是經(jīng)常用到，比如多光譜以及高光譜圖像的主成分分析要求解波段間協(xié)方差矩陣或者相關(guān)系數(shù)矩陣的特征值和特征向量。

根據(jù)普通線性代數(shù)中的概念，特征值和特征向量可以用傳統(tǒng)的方法求得，但是實(shí)際項(xiàng)目中一般都是用數(shù)值分析的方法來(lái)計(jì)算，這里介紹一下雅可比迭代法求解特征值和特征向量。

雅克比方法用于求實(shí)對(duì)稱(chēng)陣的全部特征值、特征向量。

對(duì)于實(shí)對(duì)稱(chēng)陣?A，必有正交陣?U，使

U?^TA U = D。

其中?D?是對(duì)角陣，其主對(duì)角線元?l_i?是?A?的特征值. 正交陣?U?的第?j?列是?A?的屬于?l_i?的特征向量。

原理：Jacobi 方法用平面旋轉(zhuǎn)對(duì)矩陣?A?做相似變換，化A?為對(duì)角陣，進(jìn)而求出特征值與特征向量。

既然用到了旋轉(zhuǎn)，這里就介紹一下旋轉(zhuǎn)矩陣。

對(duì)于?p?≠?q,下面定義的?n?階矩陣U_pq?是平面旋轉(zhuǎn)矩陣。

容易驗(yàn)證?U_pq是正交陣。對(duì)于向量x，U_pq?x?相當(dāng)于把坐標(biāo)軸Ox_p和?Ox_q?于所在的平面內(nèi)旋轉(zhuǎn)角度?j?.

變換過(guò)程: 在保證相似條件下，使主對(duì)角線外元素趨于零！

記?n?階方陣A?= [a_ij],? 對(duì)?A?做下面的變換：

A₁=?U_pq^TAU_pq, ? ? ? ? ? ? ? ? ? ? ? ??

????????? A_1?仍然是實(shí)對(duì)稱(chēng)陣，因?yàn)?#xff0c;U_pq^T?=U_pq^-1，知A₁與?A?的特征值相同。

?

前面說(shuō)了雅可比是一種迭代算法，所以每一步迭代時(shí)，需要求出旋轉(zhuǎn)后新的矩陣，那么新的矩陣元素如何求，這里給出具體公式如下：

由上面的一組公式可以看到:

(1)矩陣A₁?的第p?行、列與第?q?行、列中的元素發(fā)生了變化，其它行、列中的元素不變。

(2)p、q分別是前一次的迭代矩陣A的非主對(duì)角線上絕對(duì)值最大元素的行列號(hào)

(3)?j是旋轉(zhuǎn)角度，可以由下面的公式計(jì)算：

歸納可以得到雅可比迭代法求解矩陣特征值和特征向量的具體步驟如下：

(1)?初始化特征向量為對(duì)角陣V，即主對(duì)角線的元素都是1.其他元素為0。

(2)?在A的非主對(duì)角線元素中，找到絕對(duì)值最大元素?a_pq?。

(3)?用式(3.14)計(jì)算tan2j，求 cosj,?sinj?及矩陣U_{pq .}

(4)?用公式(1)-(4)求A₁；用當(dāng)前特征向量矩陣V乘以矩陣U_pq得到當(dāng)前的特征向量V。

(5)?若當(dāng)前迭代前的矩陣A的非主對(duì)角線元素中最大值小于給定的閾值e時(shí)，停止計(jì)算；否則, 令A(yù)?=?A₁?, 重復(fù)執(zhí)行(2) ~ (5)。 停止計(jì)算時(shí)，得到特征值l_i≈(A₁)?_ij?，i,j= 1,2,…,n.以及特征向量V。

(6) 這一步可選。根據(jù)特征值的大小從大到小的順序重新排列矩陣的特征值和特征向量。

?

到現(xiàn)在為止，每一步的計(jì)算過(guò)程都十分清楚了，寫(xiě)出代碼也就不是難事了，具體代碼如下：

/** * @brief 求實(shí)對(duì)稱(chēng)矩陣的特征值及特征向量的雅克比法 * 利用雅格比(Jacobi)方法求實(shí)對(duì)稱(chēng)矩陣的全部特征值及特征向量 * @param pMatrix 長(zhǎng)度為n*n的數(shù)組，存放實(shí)對(duì)稱(chēng)矩陣 * @param nDim 矩陣的階數(shù) * @param pdblVects 長(zhǎng)度為n*n的數(shù)組，返回特征向量(按列存儲(chǔ)) * @param dbEps 精度要求 * @param nJt 整型變量，控制最大迭代次數(shù) * @param pdbEigenValues 特征值數(shù)組 * @return */ 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; //迭代次數(shù)while(1){//在pMatrix的非對(duì)角線上找到最大元素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) //迭代次數(shù)超過(guò)限制break;nCount++;double dbApp = pMatrix[nRow*nDim+nRow];double dbApq = pMatrix[nRow*nDim+nCol];double dbAqq = pMatrix[nCol*nDim+nCol];//計(jì)算旋轉(zhuǎn)角度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; //p int w = i*nDim + nCol; //qdbMax = 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; //pint w = nCol*nDim + j; //qdbMax = pMatrix[u]; pMatrix[u]= pMatrix[w]*dbSinTheta + dbMax*dbCosTheta; pMatrix[w]= pMatrix[w]*dbCosTheta - dbMax*dbSinTheta; } }//計(jì)算特征向量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; } }//對(duì)特征值進(jìn)行排序以及重新排列特征向量,特征值即pMatrix主對(duì)角線上的元素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;}//設(shè)定正負(fù)號(hào)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; }

from:?http://blog.csdn.net/zhouxuguang236/article/details/40212143

總結(jié)

以上是生活随笔為你收集整理的矩阵的特征值和特征向量的雅克比算法C/C++实现的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。