雅克比公式求矩阵最大特征值（C语言版）

妖狐艹你老母 2022-05-16 09:56 231阅读 0赞

由于上次发表的自适应滤波器LMS算法涉及到要寻求矩阵最大特征值，因此一直寻思着找个C语言的雅克比算法，下面的代码也是自己试过了，望各位有问题可以提问，谢谢！！！

矩阵的特征值和特征向量是线性代数以及矩阵论中非常重要的一个概念。在遥感领域也是经常用到，比如多光谱以及高光谱图像的主成分分析要求解波段间协方差矩阵或者相关系数矩阵的特征值和特征向量。  
根据普通线性代数中的概念，特征值和特征向量可以用传统的方法求得，但是实际项目中一般都是用数值分析的方法来计算，这里介绍一下雅可比迭代法求解特征值和特征向量。  
雅克比方法用于求实对称阵的全部特征值、特征向量。  
对于实对称阵 A，必有正交阵 U，使  
U TA U = D。  
其中 D 是对角阵，其主对角线元 li 是 A 的特征值. 正交阵 U 的第 j 列是 A 的属于 li 的特征向量。  
原理：Jacobi 方法用平面旋转对矩阵 A 做相似变换，化A 为对角阵，进而求出特征值与特征向量。

***jacobi.c***
    
    #include "JACOBI.h"
    
    //求实对称矩阵的特征值及特征向量的雅格比法 
    //利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 
    //返回值小于0表示超过迭代jt次仍未达到精度要求 
    //返回值大于0表示正常返回 
    //a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值 
    //n-矩阵的阶数 
    //v-长度为n*n的数组，返回特征向量(按列存储) 
    //eps-控制精度要求 
    //jt-整型变量，控制最大迭代次数 
    
    double eejcb(double an[][F_COUNT],int n,double eps,int jt) 
    { 
        static int i,j,p,q,u,w,t,s,l; 
        static double fm,cn,sn,omega,x,y,d; 
        l=1;   //迭代计数器 
    
        while (1) 
        { 
            fm = 0.0; 
            /*找出非对角线元素的最大值*/ 
            for (i=0; i<=n-1; i++) 
            { 
                for (j=0; j<=n-1; j++) 
                { 
                    d=fabs(an[i][j]); 
                    if ((i!=j)&&(d>fm)) 
                    { 
                        fm=d; 
                        p=i; 
                        q=j; 
                    } 
                }           
            } 
    
            /*精度及迭代要求是否超过*/     
            if (fm<eps) 
            { 
                break; 
            } 
            if (l>jt) 
            { 
                break;
            } 
    
            l=l+1; 
            u=p*n+q; 
            w=p*n+p; 
            t=q*n+p; 
            s=q*n+q; 
            x=-an[p][q]; 
            y=(an[q][q]-an[p][p])/2.0; 
            omega=x/sqrt(x*x+y*y); 
            if (y<0.0) 
            { 
                omega=-omega; 
            } 
    
            sn=1.0+sqrt(1.0-omega*omega); 
            sn=omega/sqrt(2.0*sn); 
            cn=sqrt(1.0-sn*sn); 
            fm=an[p][p]; 
            an[p][p]=fm*cn*cn+an[q][q]*sn*sn+an[p][q]*omega; 
            an[q][q]=fm*sn*sn+an[q][q]*cn*cn-an[p][q]*omega; 
            an[p][q]=0.0; 
            an[q][p]=0.0; 
    
            for (j=0; j<=n-1; j++) 
            { 
                if ((j!=p)&&(j!=q)) 
                { 
                    u=p*n+j; 
                    w=q*n+j; 
                    fm=an[p][j]; 
                    an[p][j] = fm*cn+an[q][j]*sn; 
                    an[q][j] =-fm*sn+an[q][j]*cn; 
                } 
            } 
    
            for (i=0; i<=n-1; i++) 
            { 
                if ((i!=p)&&(i!=q)) 
                { 
                    u=i*n+p; 
                    w=i*n+q; 
                    fm = an[i][p]; 
                    an[i][p] = fm*cn+an[i][q]*sn; 
                    an[i][q] = -fm*sn+an[i][q]*cn; 
                } 
            }       
        } 
    
        /*计算得到最大特征值，即在an矩阵中的主对角线中寻找最大值*/    
        fm = 0.0;
        for (i=0; i<=n-1; i++)
        {
            if (an[i][i] > fm)
            {
                fm = an[i][i];
            }
        }
    
        //mu = 1/fm;  //求最大特征值的倒数 
    
        return 1/fm; 
    }

JACOBI.h
    
    #ifndef _JACOBI_H_
    #define _JACOBI_H_
    
    #include "LMS.h"
    
    
    double eejcb(double an[][F_COUNT],int n,double eps,int jt);
    
    #endif