矩阵的特征值与特征向量算法 |
您所在的位置：网站首页 › 求特征向量计算器 › 矩阵的特征值与特征向量算法 |
|
矩阵特征值和特征向量定义
A为n阶矩阵，若数λ和n维非0列向量x满足Ax=λx，那么数λ称为A的特征值，x称为A的对应于特征值λ的特征向量。式Ax=λx也可写成( A-λE)x=0，并且|λE-A|叫做A 的特征多项式。当特征多项式等于0的时候，称为A的特征方程，特征方程是一个齐次线性方程组，求解特征值的过程其实就是求解特征方程的解。 依据普通线性代数中的概念，特征值和特征向量能够用传统的方法求得，可是实际项目中一般都是用数值分析的方法来计算。 这里介绍一下雅可比(Jacobi)迭代法求解特征值和特征向量。 雅可比(Jacobi)迭代法雅克比方法用于求实对称阵的所有特征值、特征向量。Jacobi算法计算简单、稳定性好、精度高、 求得的特征向量正交性好。但当A为稀疏阵时，Givens旋转变换将破坏其稀疏性，且只能适用于 实对称矩阵。 相关知识 矩阵A与相似矩阵 B = P A P-1的特征值相同。 若矩阵Q满足QT Q = I,则称Q为正交矩阵。显然Q-1 = QT,且正交阵的乘积仍为正交阵。 若A为实对称矩阵，则存在正交阵Q，使Q A QT = diag(λ1,λ2,...,λn)，且QT的列是相应的特征向量。 实对称矩阵的特征值均为实数，且存在标准正交的特征向量系。Givens 旋转矩阵R(p,q,θ)是正交阵，其中 Givens 旋转矩阵R 原理： Jacobi 方法用平面旋转对实对称矩阵 A 做一系列旋转相似变换，从而将A约化为对角阵，进而求出特征值与特征向量。R A RT 与A元素之间的关系： 为使R A RT 为对角矩阵，可选择θ为： 当A为n阶实对称矩阵时，设A有非对角元，apq ≠ 0 ，设Givens 旋转矩阵R(p,q,θ)为： 令C = R A RT,则有： 若令C的非对角元素cpq = cqp = 0,则： C与A的元素满足下列关系： 说明经旋转变换C = R A RT后，C的对角线元素平方和比A的对角线元素平方和增加了2apq2。而C的非对角元素平方和比A的非对角元素平方和减少了2apq2。如果不断的变换下去，则最后非对角元素可趋于0，即可通过一系列旋转变换，使A与一对角阵相似。注意：某步化为0的元素在后续的步骤中可能又非0。但只要不断重复化0过程，则当K→∞时，非对角元素必趋于0。 基本思想：将实对称矩阵A经一系列正交相似变换约化为一个近似的对角阵，从而该对角阵的对角元就是A的近似特征值，由各个正交变换矩阵乘积的转置可得对应的特征向量。 旋转阵Rk+1(k=0,1,2,...)的确定： θ的计算 特征向量的计算 说明H的第i列就是A对应λi的标准正交特征向量的近似值。Ak+1的计算 计算步骤 Jacobi法的收敛性 C代码实现基于前面的C实现矩阵数据结构与计算里构造的矩阵数据结构与相应函数,用C实现了雅可比(Jacobi)迭代法求实对称矩阵的特征值与特征向量。 github源码 文件夹为Matrix，主要有两个文件：Matrix.c 、Matrix.h 设要求的矩阵为n阶的实对称矩阵，则相应的参数如下： 设定最多的迭代次数为n*n*30，若迭代次数超过限制则退出迭代。 设定精度要求为1e-10，若精度符合要求，也不再迭代。 计算后得到的结果为n+1 X n 的矩阵对象，其中第一行为特征向量，每一个特征向量对应的下面的剩余的列为其特征向量。相应代码如下： //雅克比(Jacobi)方法实现实对称矩阵的特征值和特征向量的求解 //返回矩阵第一行为特征值，特征值下面的列为对应的特征向量 Matrix *getSymmetricMatrixEigen(Matrix *m){ Matrix *resultm = NULL; Matrix *tempm = NULL; int nCount = 0,i,j; if(isSymmetricMatrix(m) == 0) return NULL; tempm = copyMatrix(m); if(!tempm) return NULL; resultm = creatIdentitySecondOrderMatrix(m->dshape); if(!resultm) return NULL; while(1){ double dbMax = *(tempm->array + 1); int nRow = 0; int nCol = 1; for(i=0;idshape.shape[2];i++){ //在矩阵非对角线上找到最大的元素 for(j=0;jdshape.shape[3];j++){ if(i != j){ double d = fabs(*(tempm->array + i*tempm->dshape.shape[3] + j)); if(d > dbMax){ dbMax = d; nRow = i; nCol = j; } } } } if(dbMax < 1e-10) break; //精度符合要求，不再迭代 if(nCount > tempm->dshape.shape[3] * tempm->dshape.shape[3] * 30) break; //迭代次数超过限制 nCount++; double dbApp = *(tempm->array + nRow*tempm->dshape.shape[3] + nRow); double dbApq = *(tempm->array + nRow*tempm->dshape.shape[3] + nCol); double dbAqq = *(tempm->array + nCol*tempm->dshape.shape[3] + nCol); //计算旋转角度 double dbAngle = 0.5*atan2(-2*dbApq,dbAqq-dbApp); double dbSinTheta = sin(dbAngle); double dbCosTheta = cos(dbAngle); double dbSin2Theta = sin(2*dbAngle); double dbCos2Theta = cos(2*dbAngle); *(tempm->array + nRow*tempm->dshape.shape[3] + nRow) = dbApp*dbCosTheta*dbCosTheta + dbAqq*dbSinTheta*dbSinTheta + 2*dbApq*dbCosTheta*dbSinTheta; *(tempm->array + nCol*tempm->dshape.shape[3] + nCol) = dbApp*dbSinTheta*dbSinTheta + dbAqq*dbCosTheta*dbCosTheta - 2*dbApq*dbCosTheta*dbSinTheta; *(tempm->array + nRow*tempm->dshape.shape[3] + nCol) = 0.5*(dbAqq-dbApp)*dbSin2Theta + dbApq*dbCos2Theta; *(tempm->array + nCol*tempm->dshape.shape[3] + nRow) = *(tempm->array + nRow*tempm->dshape.shape[3] + nCol); for(i=0;idshape.shape[3];i++){ if((i!=nCol)&&(i!=nRow)){ int u = i*tempm->dshape.shape[3] + nRow; // p int w = i*tempm->dshape.shape[3] + nCol; //q dbMax = *(tempm->array + u); *(tempm->array + u) = *(tempm->array + w) * dbSinTheta + dbMax * dbCosTheta; *(tempm->array + w) = *(tempm->array + w) * dbCosTheta - dbMax * dbSinTheta; } } for(j=0;jdshape.shape[3];j++){ if((j!=nCol)&&(j!=nRow)){ int u = nRow*tempm->dshape.shape[3] + j; //p int w = nCol*tempm->dshape.shape[3] + j; //q dbMax = *(tempm->array + u); *(tempm->array + u) = *(tempm->array + w) * dbSinTheta + dbMax * dbCosTheta; *(tempm->array + w) = *(tempm->array + w) * dbCosTheta - dbMax * dbSinTheta; } } //计算特征向量 for(i=0;idshape.shape[3];i++){ int u = i*resultm->dshape.shape[3] + nRow; // p int w = i*resultm->dshape.shape[3] + nCol; //q dbMax = *(resultm->array + u); *(resultm->array + u) = *(resultm->array + w) * dbSinTheta + dbMax*dbCosTheta; *(resultm->array + w) = *(resultm->array + w) * dbCosTheta - dbMax*dbSinTheta; } } Matrix *eigenVal = (Matrix *)malloc(sizeof(Matrix)); if(!eigenVal) return NULL; eigenVal->dshape.shape[0] = 0; eigenVal->dshape.shape[1] = 0; eigenVal->dshape.shape[2] = 0; eigenVal->dshape.shape[3] = tempm->dshape.shape[3]; eigenVal->length = tempm->dshape.shape[3]; eigenVal->size = eigenVal->length; eigenVal->array = (double *)malloc(eigenVal->size*sizeof(double)); if(!eigenVal->array){ free(eigenVal); return NULL; } for(i=0;idshape.shape[3];i++){ *(eigenVal->array + i) = *(tempm->array + i*tempm->dshape.shape[3] + i); } spliceSecondOrderMatrixRow(eigenVal,resultm); destroyMatrix(tempm); destroyMatrix(resultm); return eigenVal; } 实测注意，待求的矩阵必须是实对称矩阵。 testMatrix.c文件： #include #include #include "Matrix.h" int main(void){ Matrix *m = NULL; Matrix *m2 = NULL; int a[]={0,0,3,3}; double data[] = {1,1,1,1,2,10,1,10,100}; Dshape dshape; initDshape(&dshape,a); m = creatMatrixFromDatas(data,9,dshape); printarray(m); printf("\n"); m2 = getSymmetricMatrixEigen(m); printarray(m2); printf("\n"); destroyMatrix(m); destroyMatrix(m2); return 0; }编译： gcc Matrix.c Matrix.h testMatrix.c -o testMatrix执行testMatrix，输出： 运行结果即矩阵m的特征向量为 λ1 = 0.0946051, 对应的特征向量为: [0.70746, -0.703906, 0.063376] T λ2 = 1.8834,对应的特征向量为: [0.706669, 0.703135, -0.0788656] T λ3 = 101.022,对应的特征向量为: [0.0109521, 0.10058, 0.994869] T |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |