使用MTL库求解矩阵特征值和特征向量

关于矩阵的特征值和特征向量求解，大部分的数学运算库都进行了提供，下面是使用MTL库的接口进行封装。

#include <mtl/matrix.h>
#include <mtl/mtl.h>
#include <mtl/dense1D.h>
#include <mtl/utils.h>
#include <mtl/lu.h>

/*! 对角阵 */
typedef mtl::matrix <double, mtl::diagonal<>, mtl::dense<>, mtl::row_major>::type DiagMatrix;
/*! 对称阵 */
typedef mtl::matrix <double, mtl::symmetric<mtl::lower>, mtl::packed<mtl::external>, mtl::column_major>::type SymmMatrix;
/*! 标准矩阵，数值存在矩阵中 */
typedef mtl::matrix <double, mtl::rectangle<>, mtl::dense<>, mtl::row_major>::type Matrix;
/*! 标准矩阵，数值存在外部 */
typedef mtl::matrix <double, mtl::rectangle<>, mtl::dense<mtl::external>, mtl::row_major>::type ExtMatrix;
/*! 标准向量，数值存在矩阵中 */
typedef mtl::dense1D<double> Vector;
/*! 标准向量，数值存在外部 */
typedef mtl::external_vec<double> ExtVector;

/**
* @brief 求实对称矩阵的特征值及特征向量
* @param MMatrix	所求的矩阵
* @param EigenValues	矩阵的特征值
* @param EigenVectors	矩阵的特征向量（按照列来存），如果特征值的序号为1，那么对应的特征向量为第1列
* @param Percent	矩阵的特征值百分比
* @param AccPercent	矩阵的特征值累积百分比
* @param deps		累计误差
* @return 返回值小于0表示超过迭代jt次仍未达到精度要求，返回值大于0表示正常返回 
*/ 
bool GetMatrixEigen(Matrix MMatrix, Vector &EigenValues, Matrix &EigenVectors, Vector *Percent = NULL, Vector *AccPercent = NULL, double deps = 0.00000001);

下面是求解矩阵特征值和特征向量的实现：

/**
* @brief 求实对称矩阵的特征值及特征向量的雅格比法 
* 利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 
* @param a		长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值 
* @param n		矩阵的阶数 
* @param v		长度为n*n的数组，返回特征向量(按列存储) 
* @param eps	控制精度要求 
* @param jt		整型变量，控制最大迭代次数 
* @return 返回false表示超过迭代jt次仍未达到精度要求，返回true表示正常返回 
*/ 
bool Eejcb(double a[], int n, double v[], double eps, int jt)
{ 
	int i,j,p,q,u,w,t,s,l; 
	double fm,cn,sn,omega,x,y,d; 

	l=1; 
	//初始化特征向量矩阵使其全为0 
	for(i=0; i<=n-1; i++) 
	{   
		v[i*n+i] = 1.0; 
		for(j=0; j<=n-1; j++) 
		{ 
			if(i != j)   
				v[i*n+j]=0.0; 
		} 
	} 

	while(true)   //循环 
	{   
		fm = 0.0; 
		for(i=0; i<=n-1; i++)   //找出,矩阵a(特征值),中除对角线外其他元素的最大绝对值 
		{ 
			//这个最大值是位于a[p][q] ,等于fm 
			for(j=0; j<=n-1; j++) 
			{ 
				d = fabs(a[i*n+j]); 

				if((i!=j) && (d> fm)) 
				{ 
					fm = d;   
					p = i;   
					q = j; 
				} 
			} 
		} 

		if(fm < eps)     //精度复合要求 
			return true; //正常返回 

		if(l > jt)       //迭代次数太多 
			return false;//失败返回 

		l ++;       //   迭代计数器 
		u = p*n + q; 
		w = p*n + p;   
		t = q*n + p;   
		s = q*n + q; 
		x = -a[u]; 
		y = (a[s]-a[w])/2.0;		//x y的求法不同 
		omega = x/sqrt(x*x+y*y);	//sin2θ 

		//tan2θ=x/y = -2.0*a[u]/(a[s]-a[w]) 
		if(y < 0.0) 
			omega=-omega; 

		sn = 1.0 + sqrt(1.0-omega*omega);   
		sn = omega /sqrt(2.0*sn);		//sinθ 
		cn = sqrt(1.0-sn*sn);			//cosθ 

		fm = a[w];   //   变换前的a[w]   a[p][p] 
		a[w] = fm*cn*cn + a[s]*sn*sn + a[u]*omega; 
		a[s] = fm*sn*sn + a[s]*cn*cn - a[u]*omega; 
		a[u] = 0.0; 
		a[t] = 0.0; 

		//   以下是旋转矩阵,旋转了了p行,q行,p列,q列 
		//   但是四个特殊点没有旋转(这四个点在上述语句中发生了变化) 
		//   其他不在这些行和列的点也没变 
		//   旋转矩阵,旋转p行和q行 
		for(j=0; j<=n-1; j++) 
		{ 
			if((j!=p) && (j!=q)) 
			{ 
				u = p*n + j; 
				w = q*n + j; 
				fm = a[u]; 
				a[u] = a[w]*sn + fm*cn; 
				a[w] = a[w]*cn - fm*sn; 
			} 
		} 

		//旋转矩阵,旋转p列和q列 
		for(i=0; i<=n-1; i++) 
		{ 
			if((i!=p) && (i!=q)) 
			{ 
				u = i*n + p;   
				w = i*n + q; 
				fm = a[u]; 
				a[u]= a[w]*sn + fm*cn; 
				a[w]= a[w]*cn - fm*sn; 
			} 
		} 

		//记录旋转矩阵特征向量 
		for(i=0; i<=n-1; i++) 
		{ 
			u = i*n + p;   
			w = i*n + q; 
			fm = v[u]; 
			v[u] =v[w]*sn + fm*cn; 
			v[w] =v[w]*cn - fm*sn; 
		} 
	} 

	return true; 
}

bool GetMatrixEigen(Matrix MMatrix, Vector &EigenValues, Matrix &EigenVectors, 
					Vector *Percent, Vector *AccPercent, double deps)
{
	int nDem = MMatrix.nrows();
	double *mat = new double[nDem*nDem];	//输入矩阵，计算完成之后保存特征值
	double *eiv = new double[nDem*nDem];	//计算完成之后保存特征向量

	for (int i=0; i<nDem; i++)				//给输入矩阵和输出特征向量的矩阵赋值初始化
	{
		for (int j=0; j<nDem; j++)
		{
			mat[i*nDem+j] = MMatrix(i,j);
			eiv[i*nDem+j] = 0.0;
		}
	}

	bool rel = Eejcb(mat, nDem, eiv, deps, 100);	//计算特征值和特征向量
	if (!rel)
		return false;

	Vector TmpEigenValues(nDem);		//临时特征值
	Matrix TmpEigenVectors(nDem,nDem);	//临时特征向量

	for (int i=0; i<nDem; i++)			//赋值
	{
		for (int j=0; j<nDem; j++)
		{
			TmpEigenVectors(i,j) = eiv[i*nDem+j];

			if (i == j)
				TmpEigenValues[i] = mat[i*nDem+j];
		}
	}
	delete []mat;
	delete []eiv;

	double dSumEigenValue = 0.0;		//特征值总和
	std::map<double, int> mapEigen;
	for (size_t index=0; index<TmpEigenValues.size(); index++)
	{
		mapEigen[TmpEigenValues[index]] = index;
		dSumEigenValue += TmpEigenValues[index];
	}

	//对协方差矩阵的特征值和特征向量排序并计算累计百分比和百分比
	std::map<double, int>::reverse_iterator iter = mapEigen.rbegin();
	for(size_t ii=0; ii<TmpEigenValues.size(); ii++)
	{
		if(iter == mapEigen.rend())
			break;

		EigenValues[ii] = iter->first;
		int index = iter->second;	//获取该特征值对应的特征向量对应的列号

		if(Percent != NULL)	//计算百分比以及累积百分比
		{
			double dTemp = iter->first / dSumEigenValue;
			(*Percent)[ii] = dTemp;
			if(AccPercent != NULL)
			{
				if(ii != 0)
					(*AccPercent)[ii] = (*AccPercent)[ii-1] + dTemp;
				else
					(*AccPercent)[ii] = dTemp;
			}
		}

		double dMaxValue = ABS(TmpEigenVectors(0, index));
		int iNum = 0;
		for(int iRow=0; iRow<nDem; iRow++)	//获取特征向量中绝对值最大的位置
		{
			double dTemp = ABS(TmpEigenVectors(iRow, index));
			if(dMaxValue < dTemp)
			{
				dMaxValue = dTemp;
				iNum = iRow;
			}
		}

		bool bIsPositive = false;	//判断最大的特征向量中的值是否为正
		if(dMaxValue-TmpEigenVectors(iNum, index)<0.000000001)
			bIsPositive = true;

		for(int iRow=0; iRow<nDem; iRow++)	//确保每一个特征向量中的绝对值最大的都为正
		{
			if(!bIsPositive)
				EigenVectors(iRow, ii) = -TmpEigenVectors(iRow, index);
			else
				EigenVectors(iRow, ii) = TmpEigenVectors(iRow, index);
		}

		iter++;
	}
	return true; 
}

求解最核心的方法是雅可比方法。关于雅可比方法网上有很多资料，这里不再进行说明。

参考资料：

[1]http://osl.iu.edu/research/mtl/mtl2.php3

[2]https://zh.wikipedia.org/wiki/%E5%8D%A1%E7%88%BE%C2%B7%E9%9B%85%E5%8F%AF%E6%AF%94

[3]https://zh.wikipedia.org/wiki/%E9%9B%85%E5%8F%AF%E6%AF%94%E7%9F%A9%E9%98%B5

[4]http://boya.xmu.edu.cn/hhal/numchapt3/num_32.files/frame.htm