矩阵的特征值和特征向量的雅克比算法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;
}