关于矩阵的特征值和特征向量求解，大部分的数学运算库都进行了提供，下面是使用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] ,等于fmfor(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);elseEigenVectors(iRow, ii) = TmpEigenVectors(iRow, index);}iter++;}return true; }求解最核心的方法是雅可比方法。关于雅可比方法网上有很多资料，这里不再进行说明。
参考资料：
[1]
[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]
版权声明：本文为博主原创文章，未经博主允许不得转载。
今天不想走，明天就要跑了。
