雅可比算法求矩阵的特征值和特征向量

目的

求一个实对称矩阵的所有特征值和特征向量。

前置知识

对于一个实对称矩阵\(A\),必存在对角阵\(D\)和正交阵\(U\)满足\[D=U^TAU\]\(D\)的对角线元素为\(A\)的特征值,\(U\)的列向量为\(A\)的特征向量。

定义\(n\)阶旋转矩阵\[G(p,q,\theta)=
\begin{bmatrix}
1 & & & & & \cdots& & & & & 0\ &\ddots & & & & & & & & &\ & &1 & & & & & & & &\ & & &\cos\theta & & & &-\sin\theta & & &\ & & & &1 & &0 & & & &\ & & & & &\ddots & & & & &\ & & & &0 & &1 & & & &\ & & &\sin\theta & & & &\cos\theta & & &\ & & & & & & & &1 & &\ & & & & & & & & &\ddots &\0& & & & & & & &0 & &1
\end{bmatrix}
\]即在单位矩阵的基础上,修改\(a_{pp}=a_{qq}=\cos\theta,a_{qp}=-a_{pq}=\sin\theta\)

对于\(n\)阶向量\(\alpha\),\(\alpha\cdot G(p,q,\theta)\)的几何意义是把\(\alpha\)在与第\(p\)维坐标轴和第\(q\)维坐标轴平行的平面内旋转角度\(\theta\),并且旋转后的模长保持不变。

算法原理

大概思路使通过旋转变换使非对角线上的元素不断变小,最后得到与原矩阵相似的对角矩阵。

每次找到矩阵\(A\)绝对值最大的的非对角线元素,设为\(a_{pq}\),令\(U=G(p,q,\theta)\),将\(A\)变换为\(U^TAU\)

变换后的值为

通过令\(b_{p,q}=0\)解得\[\theta=\frac{1}{2}\arctan\frac{2a_{pq}}{a_{qq}-a_{pp}}\]特别地当\(a_{qq}=a_{pp}\)时\(\theta=\frac{\pi}{4}\)

注意到旋转操作并不会改变每个行向量或列向量的模长,即矩阵\(A\)的F-范数\(||A||_F=\sqrt{\sum_i\sum_ja_{ij}^2}\)是不变的,并且通过计算可以得出\[b_{ip}^2+b_{iq}^2=a_{ip}^2+a_{iq}^2\]从而可以得知非对角线元素的平方和变小,对角线上元素的平方和增大,故非主对角线上元素的平方和收敛。

算法流程

(1)令矩阵\(T=E\),即初始化单位矩阵

(2)找到\(A\)中绝对值最大的非对角选元素\(a_{pq}\)

(3)找到对应的角度\(\theta\),构造矩阵\(U=G(p,q,\theta)\)

(4)令\(A=U^TAU,T=TU\)

(5)不停地重复(2)到(4),直到\(a_{pq}<\epsilon\)或迭代次数超过某个限定值,则\(A\)的对角线元素近似等于\(A\)的特征值,\(T\)的列向量为\(A\)的特征向量

代码

#include<bits/stdc++.h>
using namespace std;

const int N=1005;
const double eps=1e-5;
const int lim=100;

int n,id[N];
double key[N],mat[N][N],EigVal[N],EigVec[N][N],tmpEigVec[N][N];

bool cmpEigVal(int x,int y)
{
    return key[x]>key[y];
}

void Find_Eigen(int n,double (*a)[N],double *EigVal,double (*EigVec)[N])
{
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            EigVec[i][j]=0;
    for (int i=1;i<=n;i++) EigVec[i][i]=1.0;
    int count=0;
    while (1)
    {
        //统计迭代次数
        count++;
        //找绝对值最大的元素
        double mx_val=0;
        int row_id,col_id;
        for (int i=1;i<n;i++)
            for (int j=i+1;j<=n;j++)
                if (fabs(a[i][j])>mx_val) mx_val=fabs(a[i][j]),row_id=i,col_id=j;
        if (mx_val<eps||count>lim) break;
        //进行旋转变换
        int p=row_id,q=col_id;
        double Apq=a[p][q],App=a[p][p],Aqq=a[q][q];
        double theta=0.5*atan2(-2.0*Apq,Aqq-App);
        double sint=sin(theta),cost=cos(theta);
        double sin2t=sin(2.0*theta),cos2t=cos(2.0*theta);
        a[p][p]=App*cost*cost+Aqq*sint*sint+2.0*Apq*cost*sint;
        a[q][q]=App*sint*sint+Aqq*cost*cost-2.0*Apq*cost*sint;
        a[p][q]=a[q][p]=0.5*(Aqq-App)*sin2t+Apq*cos2t;
        for (int i=1;i<=n;i++)
            if (i!=p&&i!=q)
            {
                double u=a[p][i],v=a[q][i];
                a[p][i]=u*cost+v*sint;a[q][i]=v*cost-u*sint;
                u=a[i][p],v=a[i][q];
                a[i][p]=u*cost+v*sint;a[i][q]=v*cost-u*sint;
            }
        //计算特征向量
        for (int i=1;i<=n;i++)
        {
            double u=EigVec[i][p],v=EigVec[i][q];
            EigVec[i][p]=u*cost+v*sint;EigVec[i][q]=v*cost-u*sint;
        }
    }
    //对特征值排序
    for (int i=1;i<=n;i++) id[i]=i,key[i]=a[i][i];
    std::sort(id+1,id+n+1,cmpEigVal);
    for (int i=1;i<=n;i++)
    {
        EigVal[i]=a[id[i]][id[i]];
        for (int j=1;j<=n;j++)
            tmpEigVec[j][i]=EigVec[j][id[i]];
    }
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            EigVec[i][j]=tmpEigVec[i][j];
    //特征向量为列向量
}

int main()
{
    scanf("%d",&n);
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            scanf("%lf",&mat[i][j]);
    Find_Eigen(n,mat,EigVal,EigVec);
    printf("EigenValues = ");
    for (int i=1;i<=n;i++) printf("%lf ",EigVal[i]);
    printf("\nEigenVector =\n");
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            printf("%lf%c",EigVec[i][j],j==n?'\n':' ');
    return 0;
}

原文地址:https://www.cnblogs.com/beginend/p/11749230.html

时间: 2024-10-03 14:46:02

雅可比算法求矩阵的特征值和特征向量的相关文章

利用QR算法求解矩阵的特征值和特征向量

利用QR算法求解矩阵的特征值和特征向量 为了求解一般矩阵(不是那种幼稚到shi的2 x 2矩阵)的特征值. 根据定义的话,很可能需要求解高阶方程... 这明显是个坑...高阶方程你肿么破... 折腾了好久 1.我要求特征值和特征向量. 2.找到一种算法QR分解矩阵求解特征值 3.QR矩阵分解需要Gram-schimidt正交化分解 有一种很明显的感觉,往往在现在很难有 很系统 很深入 的学习某一个学科的某一门知识. 往往学的时候"靠,学这东西有什么用""学了这么久,也不知道怎么用,不想学" 到后

矩阵的特征值和特征向量的雅克比算法C/C++实现

矩阵的特征值和特征向量是线性代数以及矩阵论中很重要的一个概念.在遥感领域也是经经常使用到.比方多光谱以及高光谱图像的主成分分析要求解波段间协方差矩阵或者相关系数矩阵的特征值和特征向量. 依据普通线性代数中的概念,特征值和特征向量能够用传统的方法求得,可是实际项目中一般都是用数值分析的方法来计算,这里介绍一下雅可比迭代法求解特征值和特征向量. 雅克比方法用于求实对称阵的所有特征值.特征向量. 对于实对称阵 A,必有正交阵 U.使 U TA U = D. 当中 D 是对角阵,其主对角线元 li 是

线性代数精华——矩阵的特征值与特征向量

今天和大家聊一个非常重要,在机器学习领域也广泛使用的一个概念--矩阵的特征值与特征向量. 我们先来看它的定义,定义本身很简单,假设我们有一个n阶的矩阵A以及一个实数\(\lambda\),使得我们可以找到一个非零向量x,满足: \[Ax=\lambda x\] 如果能够找到的话,我们就称\(\lambda\)是矩阵A的特征值,非零向量x是矩阵A的特征向量. 几何意义 光从上面的式子其实我们很难看出来什么,但是我们可以结合矩阵变换的几何意义,就会明朗很多. 我们都知道,对于一个n维的向量x来说,如

浅浅地聊一下矩阵与线性映射及矩阵的特征值与特征向量

都说矩阵其实就是线性映射,你明白不?反正一开始我是不明白的: 线性映射用矩阵表示:(很好明白的) 有两个线性空间,分别为V1与V2, V1的一组基表示为,V2的一组基表示为:(注意哦,维度可以不一样啊,反正就是线性空间啊), 1, 现在呢,有一个从V1到V2的映射F, 它可以把V1中的一组基都映射到线性空间V2中去,所以有: 用矩阵可以表示为: 2,现在我们把在V1中有一个向量A,经过映射F变为了向量B,用公式表示为:                                 所以呢,坐标

线性代数 - 05 矩阵的特征值与特征向量

线性代数 - 05 矩阵的特征值与特征向量 一.特征值与特征向量 二.矩阵的相似与矩阵的对角化 三.实对称矩阵的对角化 1.向量的内积与正交矩阵 2.实对称矩阵的特征值与特征向量 线性代数 - 05 矩阵的特征值与特征向量,码迷,mamicode.com

线性代数之矩阵的特征值与特征向量

数学上,线性变换的特征向量(本征向量)是一个非退化的向量,其方向在该变换下不变.该向量在此变换下缩放的比例称为其特征值(本征值). 一个线性变换通常可以由其特征值和特征向量完全描述.特征空间是相同特征值的特征向量的集合.“特征”一词来自德语的eigen.1904年希尔伯特首先 在这个意义下使用了这个词,更早亥尔姆霍尔兹也在相关意义下使用过该词.eigen一词可翻译为”自身的”.“特定于……的”.“有特征的”.或者“个体 的”.这显示了特征值对于定义特定的线性变换有多重要. 线性变换的特征向量是指

特征值和特征向量的几何意义、计算及其性质(一个变换(或者说矩阵)的特征向量就是这样一种向量,它经过这种特定的变换后保持方向不变,只是进行长度上的伸缩而已)

  对于任意一个矩阵,不同特征值对应的特征向量线性无关. 对于实对称矩阵或埃尔米特矩阵来说,不同特征值对应的特征向量必定正交(相互垂直).   一.特征值和特征向量的几何意义 特征值和特征向量确实有很明确的几何意义,矩阵(既然讨论特征向量的问题,当然是方阵,这里不讨论广义特征向量的概念,就是一般的特征向量)乘以一个向量的结果仍是同维数的一个向量.因此,矩阵乘法对应了一个变换,把一个向量变成同维数的另一个向量. 那么变换的效果是什么呢?这当然与方阵的构造有密切的关系,比如可以取适当的二维方阵,使得

特征值与特征向量的求法

设A为n阶方阵,如果数“ ”和n维列向量x使得关系式 成立,则称 为方阵A的特征值,非零向量x称为A对应于特征值“ ”的特征向量. 详见1.3.5和1.3.6节:特征值分解问题. 例1-89  求矩阵 的特征值和特征向量 解: >>A=[-2  1  1;0  2  0;-4  1  3]; >>[V,D]=eig(A) 结果显示: V = -0.7071   -0.2425    0.3015 0          0    0.9045 -0.7071   -0.9701  

第四章 特征值与特征向量

§4.1  特征值与特征向量 §4.1.1特征值与特征向量的概念及其计算 定义1.  设A是数域P上的一个n阶矩阵,l是一个未知量,       称为A的特征多项式,记 |(l)=| lE-A|,是一个P上的关于 l 的n次多项式,E是单位矩阵. |(l)=| lE-A|=ln+a1ln-1+-+an= 0是一个n次代数方程,称为A的特征方程. 特征方程 |(l)=| lE-A|=0的根 (如:l0) 称为A的特征根(或特征值). n次代数方程在复数域内有且仅有n 个根,而在实数域内不一定有根,