Jacobi雅克比算法计算特征向量-全网最简单

一、算法原理

• 算法涉及数据：
  矩阵V：存储特征向量。
  矩阵A：表示需要求特征向量的实对称矩阵。
• 算法过程：
• （1）初始化V为对角矩阵，即主对角线的元素是1，其他元素都为0。
• （2）在A的非主对角线元素中，找到绝对值最大的元素 Apq ，即使用p、q代表A中绝对值最大元素的下标。
• （3）计算tan2φ、cosφ、sinφ以及矩阵U。其计算公式为：
  
  
   我们遍历数组可以求得最大值Apq，即可得到下标值p和q。上式中apq即Apq，而aqq和qpp分别为App和Aqq，我们有p、q、A所以这些量都很好取得。对于矩阵U，四个特殊位置的值为三角函数值，其他对角线地方为1，非对角线地方为0。只有三角函数值的计算方法我们不知道了，这将在代码中展示。
• （4）Jacobi算法是一种迭代算法，每一次迭代会求出矩阵B，然后使用B作为新的A即目标矩阵继续进行迭代。矩阵B的计算方法如下：
  
   使用矩阵B更新A后，我们使用特征向量矩阵V乘以矩阵U得到新的特征向量，使用这个新的特征向量更新V。
• （5）遍历矩阵A求出非对角线元素中的最大值k，如果k小于算法设定的阈值e，则停止计算，否则重复执行2-5。
• （6）当停止计算时，得到特征值L和特征向量V。可以根据特征值的大小顺序重新排列矩阵的特征值和特征向量。

二、算法实现

/**
* @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)
{
	// 第一步：初始化特征向量矩阵pbdEigenValues即V
	for(int i = 0; i < nDim; i ++) 
	{   
		// 对角线为1，其他为0
		pdblVects[i*nDim+i] = 1.0f; 
		for(int j = 0; j < nDim; j ++) 
		{ 
			if(i != j)   
				pdblVects[i*nDim+j]=0.0f; 
				// 为了方便特征向量数组的形式为1维
				// 其实换算很简单，就是行数*长度+宽度
				// 因此不要觉得复杂，直接看成pblVects[i][j]即可
		} 
	} 
 	
	int nCount = 0;		// 记录迭代的次数
	while(1)
	{
		// 第二步：对矩阵pMatrix即矩阵A，求出最大索引nRow、nCol即p、q
		double dbMax = pMatrix[1];
		int nRow = 0;
		int nCol = 1; // 因为求的是非对角线上最大值，所有从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;				
 
 		// 迭代次数累计1
		nCount++;
 
 		// 有了p、q，取得App、Apq、Aqq
		double dbApp = pMatrix[nRow*nDim+nRow];
		double dbApq = pMatrix[nRow*nDim+nCol];
		double dbAqq = pMatrix[nCol*nDim+nCol];
 
		// 第三步：计算旋转角度φ,其实很简单
		// 能计算出tan2φ，那么φ=0.5arctan(tan2φ)
		double dbAngle = 0.5*atan2(-2*dbApq,dbAqq-dbApp);
		// 有了φ计算sinφ、cosφ、sin2φ、cos2φ
		double dbSinTheta = sin(dbAngle);
		double dbCosTheta = cos(dbAngle);
		double dbSin2Theta = sin(2*dbAngle);
		double dbCos2Theta = cos(2*dbAngle);
		
		// 第四步：计算B，并且使用B更新A

		// 首先可以得到Bpp、Bqq、Bpq，其实Bqp=Bpq，因此可以先更新这四个值
		// 根据第四步图中式1
		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];
 
 		// 更新A中其他值(即行和列不等于p或q的位置的值)
		// 根据式3，Bip=cosφ*Aip+sinφ*Aiq、Biq=-sinφ*Aip+cosφ*Aiq
		
		// 其中i不等于p或q，因此根据式2将计算出B中第p和q列元素值
		for(int i = 0; i < nDim; i ++) 
		{ 
			// 列不等于p或q
			if((i!=nCol) && (i!=nRow)) 
			{ 
				int u = i*nDim + nRow;	// 多次索引，每次都计算有耗费  
				int w = i*nDim + nCol;	// 所以提前算出u和w
				// u索引到第i行第p列，w索引到第i行第q列
				// 因此Au=Aip，Aw=Aiq
				
				// 	Au即Aip，因此dbMax此处就是Aip					
				dbMax = pMatrix[u]; 
				
				// Bu=Bip=cosφ*Aip+sinφ*Aiq，为了工整将sin写在前面，其实表达式的值是正确的
				pMatrix[u]= pMatrix[w]*dbSinTheta + dbMax*dbCosTheta; 
				
				// 考虑计算顺序，整个循环从上到下计算更新A矩阵的第p列和第q列
				// 当计算Bip时需要Aip和Aiq，Aip和Aiq是未被更新的原始数据，因此Bip计算正确
				// 当计算Biq时需要Aip和Aiq，上一步已经更新了Aip为Bip，使用保存原始Aip的dpMax代替即可
				pMatrix[w]= pMatrix[w]*dbCosTheta - dbMax*dbSinTheta; 
			} 
		} 
 
 		// 同理更新Bpj和Bqj，即更新第p和q行。其他的Bij=Aij，因此不需要计算了。
		for (int j = 0; j < nDim; j ++)
		{
			if((j!=nCol) && (j!=nRow)) 
			{ 
				int u = nRow*nDim + j;	// 第p行第j个元素的索引
				int w = nCol*nDim + j;	// 第q行第j个元素的索引
				dbMax = pMatrix[u]; 
				pMatrix[u]= pMatrix[w]*dbSinTheta + dbMax*dbCosTheta; 
				pMatrix[w]= pMatrix[w]*dbCosTheta - dbMax*dbSinTheta; 
			} 
		}
 
		// 第四步：当更新完A后，需要更新存储特征向量的矩阵V
		// V=VU，V最初完全是单位矩阵，U比起单位矩阵多了4个特殊元素值
		// 如果U也是单位矩阵E，那么VU=EE=E，即更新后V不会发生变化
		// 那么现在考虑U多出E的4个特殊元素值影响，VU是V行乘U列，U列4个特殊元素处于第p和第q列
		// 因此我们使用V的每一行乘以U的第p列和第p列即可
		// 记VU=S，那么需要计算的是Sip和Siq
		for(int i = 0; i < nDim; i ++) 
		{ 
			int u = i*nDim + nRow;		 // 类似上文中计算B第p和q列元素值
			int w = i*nDim + nCol;		 // Au代表Aip,Aw代表Aiq
			
			// 记录原始数据Vip
			dbMax = pdblVects[u]; 		
			
			// 根据矩阵乘法Sip等于V矩阵第i行乘以U矩阵第p列
			// 由于U矩阵第p列的两个特殊元素为Upp和Uqp
			// 因此Sip=Vip*Upp+Viq*Uqp=Vu*cosφ+Vw*sinφ，和下面的式子对的上
			pdblVects[u] = pdblVects[w]*dbSinTheta + dbMax*dbCosTheta; 
			// 同理Siq=Vip*Upq+Viq*Uqq=Vu*-sinφ+Vw*cosφ，也是一致的
			pdblVects[w] = pdblVects[w]*dbCosTheta - dbMax*dbSinTheta;
			// 为什么Sip=Vip*Upp+Viq*Uqp?首先Sip表示V第i行乘以U第p列
			// 因此表达式中V元素的行数必为i，而U元素的列数必为p，将未知数表达为x
			// 得到Sip=Vip*Uxp+Viq*uxp。为什么只考虑Vip和viq，因为此行中只有它们可乘到U中特殊元素
			// 因为Vip表示行第p个元素，因此自然乘以U列中第p个元素，即行数为p，所以此元素为Upp
		} 
 
	}
 
	// 对特征值进行排序以及重新排列特征向量,特征值即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 ) ); 
	} 
 
 	// 安装特征值大小顺序，将特征向量存储于pdbTmpVec中
	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)
	{
		// 存储第j个特征向量为矩阵中的一列
		for (int i = 0; i < nDim; i ++)
		{
			// 存储为一列因此用i*nDim+j。
			// iter->second表示特征向量在V中的列序号，因此这里是i*nDim + iter->second。
			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];
		// 如果特征向量和值小于0
		if(dSumVec<0)
		{
			// 将特征向量反向
			for(int j = 0;j < nDim; j ++)
				pdbTmpVec[j * nDim + i] *= -1;
		}
	}
 
 	// 将pdbTmpVec中数据拷贝给pdblVects即V
	memcpy(pdblVects,pdbTmpVec,sizeof(double)*nDim*nDim);
	// 释放new的空间
	delete []pdbTmpVec;
 
	return 1;
}

三、算法原理

• 参考资料。