关于矩阵的特征值和特征向量求解，大部分的数学运算库都进行了提供，下面是使用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