矩陣的特徵值和特徵向量的雅克比演算法C C 實現

2021-09-30 13:35:39 字數 4344 閱讀 3916

矩陣的特徵值和特徵向量是線性代數以及矩陣論中非常重要的乙個概念。在遙感領域也是經常用到,比如多光譜以及高光譜影象的主成分分析要求解波段間協方差矩陣或者相關係數矩陣的特徵值和特徵向量。

根據普通線性代數中的概念,特徵值和特徵向量可以用傳統的方法求得,但是實際專案中一般都是用數值分析的方法來計算,這裡介紹一下雅可比迭代法求解特徵值和特徵向量。

雅克比方法用於求實對稱陣的全部特徵值、特徵向量。

對於實對稱陣 a,必有正交陣 u,使

u ta u = d。

其中 d 是對角陣,其主對角線元 li 是 a 的特徵值. 正交陣 u 的第 j 列是 a 的屬於 li 的特徵向量。

原理:jacobi 方法用平面旋轉對矩陣 a 做相似變換,化a 為對角陣,進而求出特徵值與特徵向量。

既然用到了旋轉,這裡就介紹一下旋轉矩陣。

對於 p ≠q,下面定義的n 階矩陣upq 是平面旋轉矩陣。

容易驗證 upq是正交陣。對於向量x,upq x 相當於把座標軸oxp和 oxq 於所在的平面內旋轉角度 j .

變換過程: 在保證相似條件下,使主對角線外元素趨於零!

記 n 階方陣a = [aij],  對 a 做下面的變換:

a1= upqtaupq,                         

a1 仍然是實對稱陣,因為,upqt =upq-1,知a1與 a 的特徵值相同。

前面說了雅可比是一種迭代演算法,所以每一步迭代時,需要求出旋轉後新的矩陣,那麼新的矩陣元素如何求,這裡給出具體公式如下:

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

(1)矩陣a1 的第p 行、列與第 q 行、列中的元素發生了變化,其它行、列中的元素不變。

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

(3)j是旋轉角度,可以由下面的公式計算:

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

(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) 這一步可選。根據特徵值的大小從大到小的順序重新排列矩陣的特徵值和特徵向量。

到現在為止,每一步的計算過程都十分清楚了,寫出**也就不是難事了,具體**如下:

[cpp]view plain

copy

/*** @brief 求實對稱矩陣的特徵值及特徵向量的雅克比法 

* 利用雅格比(jacobi)方法求實對稱矩陣的全部特徵值及特徵向量 

* @param pmatrix                長度為n*n的陣列,存放實對稱矩陣

* @param ndim                   矩陣的階數 

* @param pdblvects              長度為n*n的陣列,返回特徵向量(按列儲存) 

* @param dbeps                  精度要求 

* @param njt                    整型變數,控制最大迭代次數 

* @param pdbeigenvalues         特徵值陣列

* @return 

*/boolcpcaalg::jacbicor(double* pmatrix,intndim,double*pdblvects,double*pdbeigenvalues,doubledbeps,intnjt)  

}   

intncount = 0;     //迭代次數

while(1)  

}  }if(dbmax 

break;    

if(ncount > njt)       //迭代次數超過限制

break;  

ncount++;  

doubledbapq = pmatrix[nrow*ndim+ncol];  

doubledbaqq = pmatrix[ncol*ndim+ncol];  

//計算旋轉角度

doubledbsintheta = sin(dbangle);  

doubledbcostheta = cos(dbangle);  

doubledbsin2theta = sin(2*dbangle);  

doubledbcos2theta = cos(2*dbangle);  

dbaqq*dbsintheta*dbsintheta + 2*dbapq*dbcostheta*dbsintheta;  

dbaqq*dbcostheta*dbcostheta - 2*dbapq*dbcostheta*dbsintheta;  

pmatrix[ncol*ndim+nrow] = pmatrix[nrow*ndim+ncol];  

for(inti = 0; i 

}   

for(intj = 0; j 

}  //計算特徵向量

for(inti = 0; i 

}  //對特徵值進行排序以及重新排列特徵向量,特徵值即pmatrix主對角線上的元素

std::map<double,int> mapeigen;  

for(inti = 0; i 

double*pdbtmpvec =newdouble[ndim*ndim];  

std::map<double,int>::reverse_iterator iter = mapeigen.rbegin();  

for(intj = 0; iter != mapeigen.rend(),j 

//特徵值重新排列

pdbeigenvalues[j] = iter->first;  

}  //設定正負號

for(inti = 0; i 

}  memcpy(pdblvects,pdbtmpvec,sizeof(double)*ndim*ndim);  

deletepdbtmpvec;  

return1;  

}  

矩陣的特徵值和特徵向量的雅克比演算法C C 實現

矩陣的特徵值和特徵向量是線性代數以及矩陣論中非常重要的乙個概念。在遙感領域也是經常用到,比如多光譜以及高光譜影象的主成分分析要求解波段間協方差矩陣或者相關係數矩陣的特徵值和特徵向量。根據普通線性代數中的概念,特徵值和特徵向量可以用傳統的方法求得,但是實際專案中一般都是用數值分析的方法來計算,這裡介紹...

特徵值 特殊矩陣的特徵值和特徵向量

特徵值與特徵向量 2 前 言 1 今天我們來討論一類特殊矩陣的特徵值和特徵向量。秩1 矩陣的性質希望同學們還沒有完全遺忘,正好通過今天的內容帶著大家複習下。2 i 雖然今天的矩陣不是抽象矩陣,但是想通過定義法求特徵值較為麻煩。這裡我們需要做乙個轉換 ax 0有非零解說明0是a的特徵值。ii 這裡我們...

矩陣的特徵值和特徵向量

定義 線性變換是指乙個n維列向量被左乘乙個n階矩陣後得到另乙個n維列向量,它是同維向量空間中的把乙個向量線性對映成了另乙個向量。即 y ax y,x rna aij a aij n n 如果對於數 存在乙個n維零列向量x 即x rn且x 0 使得ax x 則稱數 為矩陣a的乙個特徵值,x為矩陣a對應...