一算法原理

雅可比方法用于求解实对称矩阵的特征值和特征向量,对于实对称矩阵AAA,必有正交矩阵UUU,使得UTAU=DU^{T}AU=DUTAU=D.DDD是一个对角阵,主对角线的元素是矩阵AAA的特征值,正交矩阵UUU的每一列对应于属于矩阵DDD的主对角线对应元素的特征向量.
雅可比方法用平面旋转对矩阵AAA做相似变换,化AAA为对角阵,从而求解出特征值和特征向量.旋转矩阵UpqU_p{_q}Upq,是一个单位阵在第ppp行,第ppp列,第qqq行,第qqq列,元素为cosφcos\varphicosφ,第ppp行第qqq列为−sinφ-sin\varphi−sinφ,第qqq行第ppp列为sinφsin\varphisinφ.对于这样的平面旋转矩阵,不难验证其是一种正交矩阵.因此对于向量xxx,UpqxU_p{_q}xUpqx等同于把第ppp个坐标轴和第qqq个坐标轴共同所确定的平面旋转了φ\varphiφ度.记矩阵A1=UpqTAUpqA_1=U_p{_q}^{T}AU_p{_q}A1=UpqTAUpq.因为旋转矩阵是正交阵,因此实际上矩阵A1A_1A1与矩阵AAA是相似的,因此其特征值是相同的.
设矩阵A1A_1A1第iii行,第jjj列的元素为bijb_i{_j}bij,矩阵AAA的第iii行,第jjj列的元素为aija_i{_j}aij(i=0,1,2,...,n−1,j=0,1,2,...,n−1i=0,1,2,...,n-1,j=0,1,2,...,n-1i=0,1,2,...,n−1,j=0,1,2,...,n−1).式(1-1-1)给出了两矩阵元素之间的运算关系.
{bpp=appcos2φ+aqqsin2φ+2apqcosφsinφbqq=appsin2φ+aqqcos2φ−2apqcosφsinφbpq=bqp=12(aqq−app)sin2φ+apqcos2φbpi=apicosφ+aqisinφ,(i≠p,q)bqi=−apisinφ+aqicosφ,(i≠p,q)bjp=ajpcosφ+ajqsinφ,(j≠p,q)bjq=−ajqsinφ+ajqcosφ,(j≠p,q)bij=bji=aij,i≠p,q;j≠p,q(1-1-1)\begin{cases} b_p{_p}=a_p{_p}cos^2\varphi+a_q{_q}sin^2\varphi+2a_p{_q}cos\varphi{sin\varphi}\\ b_q{_q}=a_p{_p}sin^2\varphi+a_q{_q}cos^2\varphi-2a_p{_q}cos\varphi{sin\varphi}\\ b_p{_q}=b_q{_p}=\frac{1}2(a_q{_q}-a_p{_p})sin2\varphi+a_p{_q}cos2\varphi\\ b_p{_i}=a_p{_i}cos\varphi+a_q{_i}sin\varphi,(i\ne{p},q)\\ b_q{_i}=-a_p{_i}sin\varphi+a_q{_i}cos\varphi,(i\ne{p},q)\\ b_j{_p}=a_j{_p}cos\varphi+a_j{_q}sin\varphi,(j\ne{p},q)\\ b_j{_q}=-a_j{_q}sin\varphi+a_j{_q}cos\varphi,(j\ne{p},q)\\ b_i{_j}=b_j{_i}=a_i{_j},i{\ne}p,q;j{\ne}p,q \end{cases} \tag{1-1-1} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧bpp=appcos2φ+aqqsin2φ+2apqcosφsinφbqq=appsin2φ+aqqcos2φ−2apqcosφsinφbpq=bqp=21(aqq−app)sin2φ+apqcos2φbpi=apicosφ+aqisinφ,(i=p,q)bqi=−apisinφ+aqicosφ,(i=p,q)bjp=ajpcosφ+ajqsinφ,(j=p,q)bjq=−ajqsinφ+ajqcosφ,(j=p,q)bij=bji=aij,i=p,q;j=p,q(1-1-1)
其中有两点需要说明:(1) ppp,qqq分别是前一次的迭代矩阵的非主对角线上绝对值最大元素的行列号
(2)φ\varphiφ是旋转角度,可以由式(1-1-2)确定
tan2φ=−2apqaqq−app(1-1-2)tan2\varphi=\frac{-2a_p{_q}}{a_q{_q}-a_p{_p}} \tag{1-1-2} tan2φ=aqq−app−2apq(1-1-2)
归纳得到雅可比方法求解矩阵特征值和特征向量的具体步骤如下:
(1).初始化特征向量为对角阵V,主对角线元素为1,其他元素为0.
(2).在AAA的非主对角线的元素中,找到绝对值最大元素apqa_p{_q}apq.
(3).用式(1-1-2)计算出旋转矩阵
(4).计算矩阵A1A1A1,用当前的矩阵VVV乘旋转矩阵得到当前的特征矩阵VVV.
(5).若当前迭代的矩阵AAA的非主对角线元素最大值小于给定阈值,停止计算,否则执行上述过程.停止计算时,特征值为矩阵AAA的主对角线元素,特征矩阵为矩阵VVV.

二 C语言实现

#include<stdio.h>
#include<stdlib.h>
float** Matrix_Jac_Eig(float **array, int n, float *eig);
int Matrix_Free(float **tmp, int m, int n);
int main(void)
{int n;printf("请输入矩阵维度:\n");scanf("%d", &n);float **array = (float **)malloc(n * sizeof(float *));if (array == NULL){printf("error :申请数组内存空间失败\n");return -1;}for (int i = 0; i < n; i++){array[i] = (float *)malloc(n * sizeof(float));if (array[i] == NULL){printf("error :申请数组子内存空间失败\n");return -1;}}printf("请输入矩阵元素:\n");for (int i = 0; i < n; i++){for (int j = 0; j < n; j++){scanf("%f", &array[i][j]);}}float *eig = (float *)malloc(n * sizeof(float));float **Result = Matrix_Jac_Eig(array, n, eig);printf("特征矩阵元素:\n");for (int i = 0; i < n; i++){for (int j = 0; j < n; j++){printf("%f ", Result[i][j]);}printf("\n");}printf("特征矩阵元素:\n");for (int i = 0; i < n; i++){printf("%f \n", eig[i]);}Matrix_Free(Result, n, n);free(eig);eig = NULL;return 0;
}
float** Matrix_Jac_Eig(float **array, int n, float *eig)
{//先copy一份array在temp_mat中，因为我实在堆区申请的空间,在对其进行处理//的过程中会修改原矩阵的值,因此要存储起来,到最后函数返回的//时候再重新赋值int i, j, flag, k;flag = 0;k = 0;float sum = 0;float **temp_mat = (float **)malloc(n * sizeof(float *));for (i = 0; i < n; i++){temp_mat[i] = (float *)malloc(n * sizeof(float));}for (i = 0; i < n; i++){for (j = 0; j < n; j++){temp_mat[i][j] = array[i][j];}}//判断是否为对称矩阵for (i = 0; i < n; i++){for (j = i; j < n; j++){if (array[i][j] != array[j][i]){flag = 1;break;}}}if (flag == 1){printf("error in Matrix_Eig: 输入并非是对称矩阵:\n");return NULL;}else{//开始执行算法int p, q;float thresh = 0.0000000001;float max = array[0][1];float tan_angle, sin_angle, cos_angle;float **result = (float **)malloc(n * sizeof(float *));if (result == NULL){printf("error in Matrix_Eig:申请空间失败\n");return NULL;}float **result_temp = (float **)malloc(n * sizeof(float *));if (result_temp == NULL){printf("error in Matrix_Eig:申请空间失败\n");return NULL;}float **rot = (float **)malloc(n * sizeof(float *));if (rot == NULL){printf("error in Matrix_Eig:申请空间失败\n");return NULL;}float **mat = (float **)malloc(n * sizeof(float *));if (mat == NULL){printf("error in Matrix_Eig:申请空间失败\n");return NULL;}for (i = 0; i < n; i++){result[i] = (float *)malloc(n * sizeof(float));if (result[i] == NULL){printf("error in Matrix_Eig:申请子空间失败\n");return NULL;}result_temp[i] = (float *)malloc(n * sizeof(float));if (result_temp[i] == NULL){printf("error in Matrix_Eig:申请子空间失败\n");return NULL;}rot[i] = (float *)malloc(n * sizeof(float));if (rot[i] == NULL){printf("error in Matrix_Eig:申请子空间失败\n");return NULL;}mat[i] = (float *)malloc(n * sizeof(float));if (mat[i] == NULL){printf("error in Matrix_Eig:申请子空间失败\n");return NULL;}}for (i = 0; i < n; i++){for (j = 0; j < n; j++){if (i == j){result[i][j] = 1;}else{result[i][j] = 0;}}}for (i = 0; i < n; i++){for (j = 0; j < n; j++){if (i == j){mat[i][j] = 1;}else{mat[i][j] = 0;}}}max = array[0][1];for (i = 0; i < n; i++){for (j = 0; j < n; j++){if (i == j){continue;}else{if (fabs(array[i][j]) >= fabs(max)){max = array[i][j];p = i;q = j;}else{continue;}}}}while (fabs(max) > thresh){if (fabs(max) < thresh){break;}tan_angle = -2 * array[p][q] / (array[q][q] - array[p][p]);sin_angle = sin(0.5*atan(tan_angle));cos_angle = cos(0.5*atan(tan_angle));for (i = 0; i < n; i++){for (j = 0; j < n; j++){if (i == j){mat[i][j] = 1;}else{mat[i][j] = 0;}}}mat[p][p] = cos_angle;mat[q][q] = cos_angle;mat[q][p] = sin_angle;mat[p][q] = -sin_angle;for (i = 0; i < n; i++){for (j = 0; j < n; j++){rot[i][j] = array[i][j];}}for (j = 0; j < n; j++){rot[p][j] = cos_angle*array[p][j] + sin_angle*array[q][j];rot[q][j] = -sin_angle*array[p][j] + cos_angle*array[q][j];rot[j][p] = cos_angle*array[j][p] + sin_angle*array[j][q];rot[j][q] = -sin_angle*array[j][p] + cos_angle*array[j][q];}rot[p][p] = array[p][p] * cos_angle*cos_angle +array[q][q] * sin_angle*sin_angle +2 * array[p][q] * cos_angle*sin_angle;rot[q][q] = array[q][q] * cos_angle*cos_angle +array[p][p] * sin_angle*sin_angle -2 * array[p][q] * cos_angle*sin_angle;rot[p][q] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +array[p][q] * (2 * cos_angle*cos_angle - 1);rot[q][p] = 0.5*(array[q][q] - array[p][p]) * 2 * sin_angle*cos_angle +array[p][q] * (2 * cos_angle*cos_angle - 1);for (i = 0; i < n; i++){for (j = 0; j < n; j++){array[i][j] = rot[i][j];}}max = array[0][1];for (i = 0; i < n; i++){for (j = 0; j < n; j++){if (i == j){continue;}else{if (fabs(array[i][j]) >= fabs(max)){max = array[i][j];p = i;q = j;}else{continue;}}}}for (i = 0; i < n; i++){eig[i] = array[i][i];}for (i = 0; i < n; i++){for (j = 0; j < n; j++){sum = 0;for (k = 0; k < n; k++){sum = sum + result[i][k] * mat[k][j];}result_temp[i][j] = sum;}}for (i = 0; i < n; i++){for (j = 0; j < n; j++){result[i][j] = result_temp[i][j];}}}for (i = 0; i < n; i++){for (j = 0; j < n; j++){array[i][j] = temp_mat[i][j];}}Matrix_Free(result_temp, n, n);Matrix_Free(rot, n, n);Matrix_Free(mat, n, n);Matrix_Free(temp_mat, n, n);return result;}
}
int Matrix_Free(float **tmp, int m, int n)
{int i, j;if (tmp == NULL){return(1);}for (i = 0; i < m; i++){if (tmp[i] != NULL){free(tmp[i]);tmp[i] = NULL;}}if (tmp != NULL){free(tmp);tmp = NULL;}return(0);
}

三结果

实对称矩阵特征值特征向量求解算法C语言实现相关推荐

Eigen 矩阵的特征值特征向量求解（EVD分解）
用Eigen库求解矩阵的特征值和特征向量. 矩阵的特征值特征向量求解(EVD分解) 一.EVD分解原理 1.特征值 2.特征值分解过程二.EVD分解举例三.Eigen库实求解特征值与特征向量 ...
【机器学习】【线性代数 for PCA】矩阵与对角阵相似、一般矩阵的相似对角化、实对称矩阵的相似对角化
Note:PCA主成分分析用到实对称阵的相似对角化,用个文章复习一下相关概念和计算过程. 1.对角矩阵如果一个矩阵满足如下条件,则它就是一个对角阵: (1)是一个方阵 (2)只有对角线元素是非零元素 ...
证明：对于实对称矩阵，不同特征值对应的特征向量相互正交
前言不同特征值对应的特征向量相互正交,是实对称矩阵的一个重要属性,而且从这个属性出发可以证明实对称矩阵的另一个属性:实对称矩阵必可相似对角化.对于一个 n 维矩阵,其可相似对角化的充分且必要条件是- ...
线性代数学习笔记——第七十三讲——实对称矩阵的特征值与特征向量
1. 实对称矩阵的特征值都是实数 2. 实对称矩阵不同特征值的实特征向量相互正交
利用特征值与特征向量求解弹性力学中的主应力与主平面问题
利用特征值与特征向量求解弹性力学中的主应力与主平面问题前言一.二向应力状态 1. 莫尔圆图解法 2. 特征值与特征向量解法二.三向应力状态前言已知物体在任意一点的六个应力分量(σx,σy,σ ...
QR法求解特征值特征向量
一 QR原理理论依据:任意一个非奇异矩阵(满秩的方阵)A都可以分解为一个正交矩阵Q和一个上三角矩阵R的乘积,且当R对角元符号确定时,分解是唯一的.QR分解是一种迭代方法,迭代格式如下: 当Ak基本收 ...
AI理论知识整理（2）-对称矩阵-特征值与特征向量
把一个m×n矩阵的行,列互换得到的n×m矩阵,称为A的转置矩阵,记为A'或ATA^TAT. 矩阵转置的运算律(即性质): 1.(A')'=A 2.(A+B)'=A'+B' 3.(kA)'=kA'(k为 ...
高等数学-线性代数：已知特征值，求解特征空间的特征向量
高等数学-线性代数:已知特征值,求解特征空间的特征向量[练习]
实对称矩阵的性质_浅谈矩阵的相似对角化（一）
森屿瑾年:浅谈线性变换和矩阵之间的关系zhuanlan.zhihu.com 通过前面的讨论,我们引出了线性变换在不同基下的矩阵之间的关系,知道了线性变换在不同基下的矩阵是相似的,进而我们可以通过选取 ...
《数据结构与算法 Python语言描述》读书笔记
已经发布博客 <数据结构与算法 Python语言描述> 读书笔记第二章抽象数据类型和Python类 2.1 抽象数据类型abstract data type:ADT 2.1.1 使用编 ...

实对称矩阵特征值特征向量求解算法C语言实现

一算法原理

二 C语言实现

三结果

实对称矩阵特征值特征向量求解算法C语言实现相关推荐

最新文章

热门文章

实对称矩阵特征值特征向量求解算法C语言实现

一 算法原理

二 C语言实现

三 结果

实对称矩阵特征值特征向量求解算法C语言实现相关推荐

最新文章

热门文章

一算法原理

三结果