矩陣的特征值和特征向量的雅克比算法C/C++實現

矩陣的特征值和特征向量是線性代數以及矩陣論中很重要的一個概念。在遙感領域也是經經常使用到。比方多光譜以及高光譜圖像的主成分分析要求解波段間協方差矩陣或者相關系數矩陣的特征值和特征向量。

依據普通線性代數中的概念，特征值和特征向量能夠用傳統的方法求得，可是實際項目中一般都是用數值分析的方法來計算，這里介紹一下雅可比迭代法求解特征值和特征向量。

雅克比方法用於求實對稱陣的所有特征值、特征向量。

對於實對稱陣 A，必有正交陣 U。使

U ^TA U = D。

當中 D 是對角陣，其主對角線元 l_i 是 A 的特征值. 正交陣 U 的第 j 列是 A 的屬於 l_i 的特征向量。

原理：Jacobi 方法用平面旋轉對矩陣 A 做類似變換，化A 為對角陣，進而求出特征值與特征向量。

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

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

easy驗證 U_pq是正交陣。

對於向量x，U_pq x 相當於把坐標軸Ox_p和 Ox_q 於所在的平面內旋轉角度 j .

變換過程: 在保證類似條件下，使主對角線外元素趨於零！

記 n 階方陣A = [a_ij],  對 A 做以下的變換：

A₁= U_pq^TAU_pq,                         

          A₁ 仍然是實對稱陣，由於，U_pq^T =U_pq^-1，知A₁與 A 的特征值同樣。

 

前面說了雅可比是一種迭代算法。所以每一步迭代時，須要求出旋轉后新的矩陣，那么新的矩陣元素怎樣求，這里給出詳細公式例如以下：

由上面的一組公式能夠看到:

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

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

(3) j是旋轉角度。能夠由以下的公式計算：

歸納能夠得到雅可比迭代法求解矩陣特征值和特征向量的詳細過程例如以下：

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

(2) 在A的非主對角線元素中，找到絕對值最大元素 a_{pq 。}

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

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

(5) 若當前迭代前的矩陣A的非主對角線元素中最大值小於給定的閾值e時。停止計算；否則, 令A = A₁ , 反復運行(2) ~ (5)。 停止計算時。得到特征值 l_i≈(A₁) _ij ，i,j= 1,2,…,n.以及特征向量V。

(6) 這一步可選。

依據特征值的大小從大到小的順序又一次排列矩陣的特征值和特征向量。

 

到如今為止，每一步的計算過程都十分清楚了，寫出代碼也就不是難事了，詳細代碼例如以下：

/**
* @brief 求實對稱矩陣的特征值及特征向量的雅克比法 
* 利用雅格比(Jacobi)方法求實對稱矩陣的所有特征值及特征向量 
* @param pMatrix				長度為n*n的數組。存放實對稱矩陣
* @param nDim					矩陣的階數 
* @param pdblVects				長度為n*n的數組，返回特征向量(按列存儲) 
* @param dbEps					精度要求 
* @param nJt					整型變量。控制最大迭代次數 
* @param pdbEigenValues			特征值數組
* @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;		//迭代次數
	while(1)
	{
		//在pMatrix的非對角線上找到最大元素
		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;	//p  
				int w = i*nDim + nCol;	//q
				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;	//p
				int w = nCol*nDim + j;	//q
				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;		//p   
			int w = i*nDim + nCol;		//q
			dbMax = pdblVects[u]; 
			pdblVects[u] = pdblVects[w]*dbSinTheta + dbMax*dbCosTheta; 
			pdblVects[w] = pdblVects[w]*dbCosTheta - dbMax*dbSinTheta; 
		} 

	}

	//對特征值進行排序以及又一次排列特征向量,特征值即pMatrix主對角線上的元素
	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;
}