实对称矩阵特征值求解算法：Jacobi行循环法

72 阅读 0 评论 48 点赞

我是靠谱客的博主优美小甜瓜，最近开发中收集的这篇文章主要介绍实对称矩阵特征值求解算法：Jacobi行循环法，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

做过程故障诊断，需要计算含小于100个变量的数据的主成分，下面这个代码还是合适的，特此备份。经测试，当矩阵A阶数n=100时，耗费时间在0.2s左右。

int jacobi_loop(double** A, double** V, double* eigsv, double epsl, int maxt, int n){
/************************************************************************
* 作者: taoshaohui from Qingdao University of Science & Technology, 2010/08/17
* 函数名称: jacobi_loop
* 函数功能: Jacobi行循环法求取n阶实对称矩阵A的所有特征值和特征向量V(:,i)
* 输入参数:

A: n阶实对称矩阵, 算法结束后其第j个对角线元素为A的第j个特征值

      n:    A的阶数
      V:    n阶单位矩阵, 算法结束后其第j列存储对应第j个特征值的特征向量
     eigsv:按从大到小顺序排列的特征值向量
     epsl: 迭代收敛标准, 当A所有非对角线元素绝对值皆小于epsl时算法成功
     maxt: 算法最大迭代次数

* 输出参数: V, eigsv
* 返回值 : 整数success, 0迭代未收敛/1迭代成功
* 注意    : 本程序中，2×2Givens旋转矩阵为: [cn, sn; -sn, cn]
* 文献    : 《矩阵计算》, G. H. Golub & C. F. Van Loan, 2002, 袁亚湘译, P 494.
************************************************************************/
int success=0; // 函数返回值
double tao, t, cn, sn; // 临时变量
double maxa;
for (int it=0; it<maxt; it++)
{
  maxa=0;
  for (int p=0; p<n-1; p++)
  {
   for (int q=p+1; q<n; q++)
   {
    if (fabs(A[p][q])>maxa) // 记录非对角线元素最大值
    {
     maxa=fabs(A[p][q]);
    }

    if (fabs(A[p][q])>epsl) // 非对角线元素非0时才执行Jacobi变换
    {
     // 计算Givens旋转矩阵的重要元素:cos(theta), 即cn, sin(theta), 即sn
     tao=0.5*(A[q][q]-A[p][p])/A[p][q];

     if (tao>=0) // t为方程 t^2 + 2*t*tao - 1 = 0的根, 取绝对值较小的根为t
     {
      t=-tao+sqrt(1+tao*tao);
     }
     else
     {
      t=-tao-sqrt(1+tao*tao);
     }
     cn=1/sqrt(1+t*t);
     sn=t*cn;

     // Givens旋转矩阵之转置左乘A, 即更新A的p行和q行
     for (int j=0; j<n; j++)
     {
      double apj=A[p][j];
      double aqj=A[q][j];
      A[p][j]=cn*apj-sn*aqj;
      A[q][j]=sn*apj+cn*aqj;
     }

     // Givens旋转矩阵右乘A, 即更新A的p列和q列
     for (int i=0; i<n; i++)
     {
      double aip=A[i][p];
      double aiq=A[i][q];
      A[i][p]=cn*aip-sn*aiq;
      A[i][q]=sn*aip+cn*aiq;
     }

     // 更新特征向量存储矩阵V, V=J0×J1×J2...×Jit, 所以每次只更新V的p, q两列
     for (int i=0; i<n; i++)
     {
      double vip=V[i][p];
      double viq=V[i][q];
      V[i][p]=cn*vip-sn*viq;
      V[i][q]=sn*vip+cn*viq;
     }
    }

}

  if (maxa<epsl) // 非对角线元素已小于收敛标准，迭代结束
  {
   // 特征值向量
   for (int j=0; j<n; j++)
   {
    eigsv[j]=A[j][j];
   }
   // 对特征值向量从大到小进行排序, 并调整特征向量顺序 (直接插入法)

double* tmp=new double[n];
   for (int j=1; j<n; j++)
   {
    int i=j;
    double a=eigsv[j];
        for (int k=0; k<n; k++)
    {
     tmp[k]=V[k][j];
    }
    while(i>0 && eigsv[i-1]<a){
     eigsv[i]=eigsv[i-1];
     for (int k=0; k<n; k++)
     {
      V[k][i]=V[k][i-1];
     }
     i--;
    }
    eigsv[i]=a;
    for (int k=0; k<n; k++)
    {
     V[k][i]=tmp[k];
    }
   }