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

2024-09-03 12:18

本文主要是介绍矩阵的特征值和特征向量的雅克比算法C/C++实现,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

矩阵的特征值和特征向量是线性代数以及矩阵论中非常重要的一个概念。在遥感领域也是经常用到,比如多光谱以及高光谱图像的主成分分析要求解波段间协方差矩阵或者相关系数矩阵的特征值和特征向量。

根据普通线性代数中的概念,特征值和特征向量可以用传统的方法求得,但是实际项目中一般都是用数值分析的方法来计算,这里介绍一下雅可比迭代法求解特征值和特征向量。

雅克比方法用于求实对称阵的全部特征值、特征向量。

对于实对称阵 A,必有正交阵 U,使

TA U = D

其中 是对角阵,其主对角线元 li 是 的特征值. 正交阵 的第 列是 的属于 li 的特征向量。

原理:Jacobi 方法用平面旋转对矩阵 做相似变换,化为对角阵,进而求出特征值与特征向量。

既然用到了旋转,这里就介绍一下旋转矩阵。

对于  q,下面定义的 n 阶矩阵Upq 平面旋转矩阵。

容易验证 Upq是正交阵。对于向量xUpq x 相当于把坐标轴Oxp和 Oxq 于所在的平面内旋转角度 j .

变换过程: 在保证相似条件下,使主对角线外元素趋于零!

记 n 阶方阵= [aij],  对 A 做下面的变换:

A1UpqTAUpq,                         

          A1 仍然是实对称阵,因为,UpqT =Upq-1,知A1与 的特征值相同。

 

前面说了雅可比是一种迭代算法,所以每一步迭代时,需要求出旋转后新的矩阵,那么新的矩阵元素如何求,这里给出具体公式如下:

由上面的一组公式可以看到:

(1)矩阵A1 的第p 行、列与第 q 行、列中的元素发生了变化,其它行、列中的元素不变。

(2)p、q分别是前一次的迭代矩阵A的非主对角线上绝对值最大元素的行列号

(3) j是旋转角度,可以由下面的公式计算:

归纳可以得到雅可比迭代法求解矩阵特征值和特征向量的具体步骤如下:

(1) 初始化特征向量为对角阵V,即主对角线的元素都是1.其他元素为0。

(2) 在A的非主对角线元素中,找到绝对值最大元素 apq 。

(3) 用式(3.14)计算tan2j,求 cosj, sin及矩阵Upq .

(4) 用公式(1)-(4)求A1;用当前特征向量矩阵V乘以矩阵Upq得到当前的特征向量V。

(5) 若当前迭代前的矩阵A的非主对角线元素中最大值小于给定的阈值e时,停止计算;否则, 令A = A1 , 重复执行(2) ~ (5)。 停止计算时,得到特征值 li≈(A1) ij ,i,j= 1,2,…,n.以及特征向量V。

(6) 这一步可选。根据特征值的大小从大到小的顺序重新排列矩阵的特征值和特征向量。

 

到现在为止,每一步的计算过程都十分清楚了,写出代码也就不是难事了,具体代码如下:

 

[cpp] view plain copy

  1. /** 
  2. * @brief 求实对称矩阵的特征值及特征向量的雅克比法  
  3. * 利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量  
  4. * @param pMatrix                长度为n*n的数组,存放实对称矩阵 
  5. * @param nDim                   矩阵的阶数  
  6. * @param pdblVects              长度为n*n的数组,返回特征向量(按列存储)  
  7. * @param dbEps                  精度要求  
  8. * @param nJt                    整型变量,控制最大迭代次数  
  9. * @param pdbEigenValues         特征值数组 
  10. * @return  
  11. */  
  12. bool CPCAAlg::JacbiCor(double * pMatrix,int nDim, double *pdblVects, double *pdbEigenValues, double dbEps,int nJt)  
  13. {  
  14.     for(int i = 0; i < nDim; i ++)   
  15.     {     
  16.         pdblVects[i*nDim+i] = 1.0f;   
  17.         for(int j = 0; j < nDim; j ++)   
  18.         {   
  19.             if(i != j)     
  20.                 pdblVects[i*nDim+j]=0.0f;   
  21.         }   
  22.     }   
  23.   
  24.     int nCount = 0;     //迭代次数  
  25.     while(1)  
  26.     {  
  27.         //在pMatrix的非对角线上找到最大元素  
  28.         double dbMax = pMatrix[1];  
  29.         int nRow = 0;  
  30.         int nCol = 1;  
  31.         for (int i = 0; i < nDim; i ++)          //行  
  32.         {  
  33.             for (int j = 0; j < nDim; j ++)      //列  
  34.             {  
  35.                 double d = fabs(pMatrix[i*nDim+j]);   
  36.   
  37.                 if((i!=j) && (d> dbMax))   
  38.                 {   
  39.                     dbMax = d;     
  40.                     nRow = i;     
  41.                     nCol = j;   
  42.                 }   
  43.             }  
  44.         }  
  45.   
  46.         if(dbMax < dbEps)     //精度符合要求   
  47.             break;    
  48.   
  49.         if(nCount > nJt)       //迭代次数超过限制  
  50.             break;  
  51.   
  52.         nCount++;  
  53.   
  54.         double dbApp = pMatrix[nRow*nDim+nRow];  
  55.         double dbApq = pMatrix[nRow*nDim+nCol];  
  56.         double dbAqq = pMatrix[nCol*nDim+nCol];  
  57.   
  58.         //计算旋转角度  
  59.         double dbAngle = 0.5*atan2(-2*dbApq,dbAqq-dbApp);  
  60.         double dbSinTheta = sin(dbAngle);  
  61.         double dbCosTheta = cos(dbAngle);  
  62.         double dbSin2Theta = sin(2*dbAngle);  
  63.         double dbCos2Theta = cos(2*dbAngle);  
  64.   
  65.         pMatrix[nRow*nDim+nRow] = dbApp*dbCosTheta*dbCosTheta +   
  66.             dbAqq*dbSinTheta*dbSinTheta + 2*dbApq*dbCosTheta*dbSinTheta;  
  67.         pMatrix[nCol*nDim+nCol] = dbApp*dbSinTheta*dbSinTheta +   
  68.             dbAqq*dbCosTheta*dbCosTheta - 2*dbApq*dbCosTheta*dbSinTheta;  
  69.         pMatrix[nRow*nDim+nCol] = 0.5*(dbAqq-dbApp)*dbSin2Theta + dbApq*dbCos2Theta;  
  70.         pMatrix[nCol*nDim+nRow] = pMatrix[nRow*nDim+nCol];  
  71.   
  72.         for(int i = 0; i < nDim; i ++)   
  73.         {   
  74.             if((i!=nCol) && (i!=nRow))   
  75.             {   
  76.                 int u = i*nDim + nRow;  //p    
  77.                 int w = i*nDim + nCol;  //q  
  78.                 dbMax = pMatrix[u];   
  79.                 pMatrix[u]= pMatrix[w]*dbSinTheta + dbMax*dbCosTheta;   
  80.                 pMatrix[w]= pMatrix[w]*dbCosTheta - dbMax*dbSinTheta;   
  81.             }   
  82.         }   
  83.   
  84.         for (int j = 0; j < nDim; j ++)  
  85.         {  
  86.             if((j!=nCol) && (j!=nRow))   
  87.             {   
  88.                 int u = nRow*nDim + j;  //p  
  89.                 int w = nCol*nDim + j;  //q  
  90.                 dbMax = pMatrix[u];   
  91.                 pMatrix[u]= pMatrix[w]*dbSinTheta + dbMax*dbCosTheta;   
  92.                 pMatrix[w]= pMatrix[w]*dbCosTheta - dbMax*dbSinTheta;   
  93.             }   
  94.         }  
  95.   
  96.         //计算特征向量  
  97.         for(int i = 0; i < nDim; i ++)   
  98.         {   
  99.             int u = i*nDim + nRow;      //p     
  100.             int w = i*nDim + nCol;      //q  
  101.             dbMax = pdblVects[u];   
  102.             pdblVects[u] = pdblVects[w]*dbSinTheta + dbMax*dbCosTheta;   
  103.             pdblVects[w] = pdblVects[w]*dbCosTheta - dbMax*dbSinTheta;   
  104.         }   
  105.   
  106.     }  
  107.   
  108.     //对特征值进行排序以及重新排列特征向量,特征值即pMatrix主对角线上的元素  
  109.     std::map<double,int> mapEigen;  
  110.     for(int i = 0; i < nDim; i ++)   
  111.     {     
  112.         pdbEigenValues[i] = pMatrix[i*nDim+i];  
  113.   
  114.         mapEigen.insert(make_pair( pdbEigenValues[i],i ) );  
  115.     }   
  116.   
  117.     double *pdbTmpVec = new double[nDim*nDim];  
  118.     std::map<double,int>::reverse_iterator iter = mapEigen.rbegin();  
  119.     for (int j = 0; iter != mapEigen.rend(),j < nDim; ++iter,++j)  
  120.     {  
  121.         for (int i = 0; i < nDim; i ++)  
  122.         {  
  123.             pdbTmpVec[i*nDim+j] = pdblVects[i*nDim + iter->second];  
  124.         }  
  125.   
  126.         //特征值重新排列  
  127.         pdbEigenValues[j] = iter->first;  
  128.     }  
  129.   
  130.     //设定正负号  
  131.     for(int i = 0; i < nDim; i ++)   
  132.     {  
  133.         double dSumVec = 0;  
  134.         for(int j = 0; j < nDim; j ++)  
  135.             dSumVec += pdbTmpVec[j * nDim + i];  
  136.         if(dSumVec<0)  
  137.         {  
  138.             for(int j = 0;j < nDim; j ++)  
  139.                 pdbTmpVec[j * nDim + i] *= -1;  
  140.         }  
  141.     }  
  142.   
  143.     memcpy(pdblVects,pdbTmpVec,sizeof(double)*nDim*nDim);  
  144.     delete []pdbTmpVec;  
  145.   
  146.     return 1;  
  147. }  

这篇关于矩阵的特征值和特征向量的雅克比算法C/C++实现的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1132936

相关文章

pandas中位数填充空值的实现示例

《pandas中位数填充空值的实现示例》中位数填充是一种简单而有效的方法,用于填充数据集中缺失的值,本文就来介绍一下pandas中位数填充空值的实现,具有一定的参考价值,感兴趣的可以了解一下... 目录什么是中位数填充?为什么选择中位数填充?示例数据结果分析完整代码总结在数据分析和机器学习过程中,处理缺失数

Golang HashMap实现原理解析

《GolangHashMap实现原理解析》HashMap是一种基于哈希表实现的键值对存储结构,它通过哈希函数将键映射到数组的索引位置,支持高效的插入、查找和删除操作,:本文主要介绍GolangH... 目录HashMap是一种基于哈希表实现的键值对存储结构,它通过哈希函数将键映射到数组的索引位置,支持

Pandas使用AdaBoost进行分类的实现

《Pandas使用AdaBoost进行分类的实现》Pandas和AdaBoost分类算法,可以高效地进行数据预处理和分类任务,本文主要介绍了Pandas使用AdaBoost进行分类的实现,具有一定的参... 目录什么是 AdaBoost?使用 AdaBoost 的步骤安装必要的库步骤一:数据准备步骤二:模型

使用Pandas进行均值填充的实现

《使用Pandas进行均值填充的实现》缺失数据(NaN值)是一个常见的问题,我们可以通过多种方法来处理缺失数据,其中一种常用的方法是均值填充,本文主要介绍了使用Pandas进行均值填充的实现,感兴趣的... 目录什么是均值填充?为什么选择均值填充?均值填充的步骤实际代码示例总结在数据分析和处理过程中,缺失数

Java对象转换的实现方式汇总

《Java对象转换的实现方式汇总》:本文主要介绍Java对象转换的多种实现方式,本文通过实例代码给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价值,需要的朋友参考下吧... 目录Java对象转换的多种实现方式1. 手动映射(Manual Mapping)2. Builder模式3. 工具类辅助映

Go语言开发实现查询IP信息的MCP服务器

《Go语言开发实现查询IP信息的MCP服务器》随着MCP的快速普及和广泛应用,MCP服务器也层出不穷,本文将详细介绍如何在Go语言中使用go-mcp库来开发一个查询IP信息的MCP... 目录前言mcp-ip-geo 服务器目录结构说明查询 IP 信息功能实现工具实现工具管理查询单个 IP 信息工具的实现服

SpringBoot基于配置实现短信服务策略的动态切换

《SpringBoot基于配置实现短信服务策略的动态切换》这篇文章主要为大家详细介绍了SpringBoot在接入多个短信服务商(如阿里云、腾讯云、华为云)后,如何根据配置或环境切换使用不同的服务商,需... 目录目标功能示例配置(application.yml)配置类绑定短信发送策略接口示例:阿里云 & 腾

python实现svg图片转换为png和gif

《python实现svg图片转换为png和gif》这篇文章主要为大家详细介绍了python如何实现将svg图片格式转换为png和gif,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录python实现svg图片转换为png和gifpython实现图片格式之间的相互转换延展:基于Py

Python利用ElementTree实现快速解析XML文件

《Python利用ElementTree实现快速解析XML文件》ElementTree是Python标准库的一部分,而且是Python标准库中用于解析和操作XML数据的模块,下面小编就来和大家详细讲讲... 目录一、XML文件解析到底有多重要二、ElementTree快速入门1. 加载XML的两种方式2.

Java的栈与队列实现代码解析

《Java的栈与队列实现代码解析》栈是常见的线性数据结构,栈的特点是以先进后出的形式,后进先出,先进后出,分为栈底和栈顶,栈应用于内存的分配,表达式求值,存储临时的数据和方法的调用等,本文给大家介绍J... 目录栈的概念(Stack)栈的实现代码队列(Queue)模拟实现队列(双链表实现)循环队列(循环数组