矩阵特征值和特征向量定义
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;i<tempm->dshape.shape[2];i++){ //在矩阵非对角线上找到最大的元素
for(j=0;j<tempm->dshape.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;i<tempm->dshape.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;j<tempm->dshape.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;i<resultm->dshape.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;i<resultm->dshape.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 <stdio.h>
#include <malloc.h>
#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
