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

2024-09-03 12:18
文章标签 算法 c++ 实现 矩阵 特征值 特征向量 雅克

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

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

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

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

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

U TA U = D。

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

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

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

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

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

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

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

A1= UpqTAUpq,                         

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

 

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

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

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

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

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

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

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

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

(3) 用式(3.14)计算tan2j，求 cosj, sinj 及矩阵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。

相关文章

不懂推荐算法也能设计推荐系统

不懂推荐算法也能设计推荐系统

本文以商业化应用推荐为例，告诉我们不懂推荐算法的产品，也能从产品侧出发， 设计出一款不错的推荐系统。 相信很多新手产品，看到算法二字，多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法（注：传统是指科班出身的产品都会接触过）。但对于推荐算法，多数产品对着网上搜到的资源，都会无从下手。特别当某些推荐算法 和 “AI”扯上关系后，更是加大了理解的难度。 但，不了解推荐算法，就无法做推荐系
阅读更多...
hdu1043（八数码问题，广搜 + hash（实现状态压缩） ）

hdu1043（八数码问题，广搜 + hash（实现状态压缩） ）

利用康拓展开将一个排列映射成一个自然数，然后就变成了普通的广搜题。 #include<iostream>#include<algorithm>#include<string>#include<stack>#include<queue>#include<map>#include<stdio.h>#include<stdlib.h>#include<ctype.h>#inclu
阅读更多...
康拓展开（hash算法中会用到）

康拓展开（hash算法中会用到）

康拓展开是一个全排列到一个自然数的双射（也就是某个全排列与某个自然数一一对应） 公式： X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! 其中,a[i]为整数,并且0<=a[i]<i,1<=i<=n。（a[i]在不同应用中的含义不同）； 典型应用： 计算当前排列在所有由小到大全排列中的顺序，也就是说求当前排列是第
阅读更多...
【C++ Primer Plus习题】13.4

【C++ Primer Plus习题】13.4

大家好,这里是国中之林! ❥前些天发现了一个巨牛的人工智能学习网站，通俗易懂，风趣幽默，忍不住分享一下给大家。点击跳转到网站。有兴趣的可以点点进去看看← 问题: 解答: main.cpp #include <iostream>#include "port.h"int main() {Port p1;Port p2("Abc", "Bcc", 30);std::cout <<
阅读更多...
csu 1446 Problem J Modified LCS （扩展欧几里得算法的简单应用）

csu 1446 Problem J Modified LCS （扩展欧几里得算法的简单应用）

这是一道扩展欧几里得算法的简单应用题，这题是在湖南多校训练赛中队友ac的一道题，在比赛之后请教了队友，然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意：把题意转变下就变成了：求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点：求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个
阅读更多...
C++包装器

C++包装器

包装器 在 C++ 中，“包装器”通常指的是一种设计模式或编程技巧，用于封装其他代码或对象，使其更易于使用、管理或扩展。包装器的概念在编程中非常普遍，可以用于函数、类、库等多个方面。下面是几个常见的 “包装器” 类型： 1. 函数包装器 函数包装器用于封装一个或多个函数，使其接口更统一或更便于调用。例如，std::function 是一个通用的函数包装器，它可以存储任意可调用对象（函数、函数
阅读更多...
综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

LntonAIServer视频质量诊断功能中的抖动检测是一个专门针对视频稳定性进行分析的功能。抖动通常是指视频帧之间的不必要运动，这种运动可能是由于摄像机的移动、传输中的错误或编解码问题导致的。抖动检测对于确保视频内容的平滑性和观看体验至关重要。 优势 1. 提高图像质量 - 清晰度提升：减少抖动，提高图像的清晰度和细节表现力，使得监控画面更加真实可信。 - 细节增强：在低光条件下，抖
阅读更多...
C++11第三弹：lambda表达式 | 新的类功能 | 模板的可变参数

C++11第三弹：lambda表达式 | 新的类功能 | 模板的可变参数

🌈个人主页： 南桥几晴秋 🌈C++专栏： 南桥谈C++ 🌈C语言专栏： C语言学习系列 🌈Linux学习专栏： 南桥谈Linux 🌈数据结构学习专栏： 数据结构杂谈 🌈数据库学习专栏： 南桥谈MySQL 🌈Qt学习专栏： 南桥谈Qt 🌈菜鸡代码练习： 练习随想记录 🌈git学习： 南桥谈Git 🌈🌈🌈🌈🌈🌈🌈🌈🌈🌈🌈🌈🌈�
阅读更多...
【数据结构】——原来排序算法搞懂这些就行，轻松拿捏

【数据结构】——原来排序算法搞懂这些就行，轻松拿捏

前言：快速排序的实现最重要的是找基准值，下面让我们来了解如何实现找基准值 基准值的注释：在快排的过程中，每一次我们要取一个元素作为枢纽值，以这个数字来将序列划分为两部分。 在此我们采用三数取中法，也就是取左端、中间、右端三个数，然后进行排序，将中间数作为枢纽值。 快速排序实现主框架： //快速排序 void QuickSort(int* arr, int left, int rig
阅读更多...
【C++】_list常用方法解析及模拟实现

【C++】_list常用方法解析及模拟实现

相信自己的力量，只要对自己始终保持信心，尽自己最大努力去完成任何事，就算事情最终结果是失败了，努力了也不留遗憾。💓💓💓 目录   ✨说在前面 🍋知识点一：什么是list？ •🌰1.list的定义 •🌰2.list的基本特性 •🌰3.常用接口介绍 🍋知识点二：list常用接口 •🌰1.默认成员函数 🔥构造函数（⭐） 🔥析构函数 •🌰2.list对象
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2