特征向量的数值求法

来源：互联网发布：弱智吧知乎编辑：程序博客网时间：2024/04/30 10:10

-对于正定的对称矩阵，奇异值等于特征值，奇异向量等于特征向量。在这种情况下用奇异值分解就把特征值和特征向量求出来了。但是只要是方阵，它就有特征值和特征向量，对于一般的方阵，特征值和特征向量怎么求呢（当然我指的是数值求法）？这就要用本文即将介绍的“幂法”。

Power Method幂法

Definition 如果λ₁是矩阵A的所有特征值中绝对值最大的那一个，则称λ₁是A的主特征值；与λ₁对应的特征向量v₁是A的主特征向量。

幂法是用来计算方阵的主特征值（即绝对值最大的特征值）和主特征向量的。由此延伸出来的反幂法用来计算在给定点附近的特征值和特征向量（下文把“特征值和特征向量”简称为“特征对”）。

Definition 特征向量V的归一化是指：V的每一个元素除以V中绝对值最大的那个元素。

Theorem(Power Method) 设A是n×n的方阵，有n个不同的特征值，且|λ1| > |λ2| >= |λ3| >= ... >= |λn|。选择一个合适的X₀，序列 $\{X_k=(x_1^{(k)},x_2^{(k)},...,x_n^{(k)})^T\}$ 和{c_k}由下列递归式产生：

$Y_k=AX_k\qquad(1)$

$X_{k+1}=\frac{1}{c_{k+1}}Y_k\qquad(2)$

其中 $c_{k+1}=\underset{1\le{i}\le{n}}{max}\{|x_i^{(k)}|\}$

经过多次迭代后X_k趋于主特征向量V₁，c_k趋于主特征值λ₁。

Remark 如果X₀选取的刚好是一个特征向量，且X₀又不是主特征向量，则X₀需要重新选取。

Speed of Convergence收敛加速

幂法的收敛速度取决于 $\frac{|\lambda_1|}{|\lambda_2|}$ ，也就是说它的收敛速度是线性的。Aitken $\bigtriangleup{^2}$ 加速法可用于任何线性收敛的序列当中，它采用的加速方式是：

用Aitken来加速我们的幂法，X_k的调整公式为：

$\hat{x_j^{(k)}}=x_j^{(k)}-\frac{(x_j^{(k)}-x_j^{(k-1)})^2}{x_j^{(k)}-2x_j^{(k-1)}+x_j^{(k-2)}},x_j^{(k)}\ne1$

Shifted-Inverse Power Method平移反幂法

使用位移反幂法首先需要提供一个好的起始点，这个点要接近一个特征向量，然后我们的位移反幂法才能够以更高的精度算出这个特征向量。QM和Given's method可以用来获得这种初始点，这里不介绍了。在实际情况中，特征值可能是复数，有多个特征相同或很接近，这都会使得计算变得很复杂需要更高级的算法。我们只考虑简单的情况，所有的特征值各不相同。

Theorem (Shifting Eigenvalues) 假设λ和V是A的一个特征对，a是任何常量，那么λ-a，V就是矩阵 $A-aI$ 的一个特征对。

Theorem (Inverse Eigenvalues) 假设λ和V是A的一个特征对， $\lambda\ne{0}$ ，那么 $\frac{1}{\lambda}$ ，V是A^-1的一个特征对。

Theorem (Shifted-Inverse Eigenvalues) 假设λ和V是A的一个特征对， $a\ne\lambda$ ，那么 $\frac{1}{\lambda-a}$ ，V是 $(A-aI)^{-1}$ 的一个特征对。

Theorem (Shifted-Inverse Power Method) 假设A是一个n×n的矩阵，有n个互不相同的特征值λ₁，λ₂，...，λ_n。对于其中之一个特征值λ_j，可以选择一个常数a，使得 $u_1=\frac{1}{\lambda_j-a}$ 是 $(A-aI)^{-1}$ 的主特征值。进一步地，如果选择一个合适的X₀，序列 $\{X_k=(x_1^{(k)},x_2^{(k)},...,x_n^{(k)})^T\}$ 和{c_k}可由下列递归式产生：

$Y_k=(A-aI)^{-1}X_k\qquad(3)$

$X_{k+1}=\frac{1}{c_{k+1}}Y_k\qquad(4)$

其中 $c_{k+1}=\underset{1\le{i}\le{n}}{max}\{|x_i^{(k)}|\}$

经过多次迭代后X_k趋于 $(A-aI)^{-1}$ 的主特征向量V_j，X_k同时也是A的主特征向量，c_k趋于 $(A-aI)^{-1}$ 的主特征值u₁。最终我们可以求出A的主特征值：

$\lambda_j=\frac{1}{u_1}+a$

在求解(3)式时需要解一个线性方程组 $(A-aI)Y_k=X_k$ ，常用的方法是雅可比迭代和高斯-赛德尔迭代。当然你也可以用 $A^{-1}(A|I)=(I|A^{-1})$ 的方法进行初等行变换来求得矩阵的逆，那样就不用解线性方程组。不过你要衡量哪种方式快一些，而且矩阵的逆不存在怎么办。

高斯-赛德尔迭代公式为：

$x_i^{k+1}=\frac{b_i-\sum_{j=1}^{i-1}{a_{ij}x_j^{k+1}}-\sum_{j=i+1}^n{a_{ij}x_j^k}}{a_{ii}}\qquad(1\le{i}\le{n})$

注意a_ii即系数矩阵主对角线上的元素不能有0，否则需要事先进行行变换，把0移走。

Exercise

用幂法求一个矩阵的主特征对。

取最大迭代次数为50，收敛时误差为0.000001。

matrix.h

#ifndef _MATRIX_H#define _MATRIX_H   #include<assert.h>#include<stdlib.h>#include<stdio.h> //初始化一个二维矩阵double** getMatrix(int rows,int columns){    double **rect=(double**)calloc(rows,sizeof(double*));    int i;    for(i=0;i<rows;++i)        rect[i]=(double*)calloc(columns,sizeof(double));    return rect;} //返回一个单位矩阵double** getIndentityMatrix(int rows){    double** IM=getMatrix(rows,rows);    int i;    for(i=0;i<rows;++i)        IM[i][i]=1.0;    return IM;}   //返回一个矩阵的副本double** copyMatrix(double** matrix,int rows,int columns){    double** rect=getMatrix(rows,columns);    int i,j;    for(i=0;i<rows;++i)        for(j=0;j<columns;++j)            rect[i][j]=matrix[i][j];    return rect;} //从一个一维矩阵得到一个二维矩阵void getFromArray(double** matrix,int rows,int columns,double *arr){    int i,j,k=0;    for(i=0;i<rows;++i){        for(j=0;j<columns;++j){            matrix[i][j]=arr[k++];        }    }} //打印二维矩阵void printMatrix(double** matrix,int rows,int columns){    int i,j;    for(i=0;i<rows;++i){        for(j=0;j<columns;++j){            printf("%-10f\t",matrix[i][j]);        }        printf("\n");    }} //释放二维矩阵void freeMatrix(double** matrix,int rows){    int i;    for(i=0;i<rows;++i)        free(matrix[i]);    free(matrix);}   //获取二维矩阵的某一行double* getRow(double **matrix,int rows,int columns,int index){    assert(index<rows);    double *rect=(double*)calloc(columns,sizeof(double));    int i;    for(i=0;i<columns;++i)        rect[i]=matrix[index][i];    return rect;} //获取二维矩阵的某一列  double* getColumn(double **matrix,int rows,int columns,int index){    assert(index<columns);    double *rect=(double*)calloc(rows,sizeof(double));    int i;    for(i=0;i<rows;++i)        rect[i]=matrix[i][index];    return rect;} //设置二维矩阵的某一列   void setColumn(double **matrix,int rows,int columns,int index,double *arr){    assert(index<columns);    int i;    for(i=0;i<rows;++i)        matrix[i][index]=arr[i];} //交换矩阵的某两列void exchangeColumn(double **matrix,int rows,int columns,int i,int j){    assert(i<columns);    assert(j<columns);    int row;    for(row=0;row<rows;++row){        double tmp=matrix[row][i];        matrix[row][i]=matrix[row][j];        matrix[row][j]=tmp;    }}  //得到矩阵的转置double** getTranspose(double **matrix,int rows,int columns){    double **rect=getMatrix(columns,rows);    int i,j;    for(i=0;i<columns;++i){        for(j=0;j<rows;++j){            rect[i][j]=matrix[j][i];        }    }    return rect;} //计算两向量内积double vectorProduct(double *vector1,double *vector2,int len){    double rect=0.0;    int i;    for(i=0;i<len;++i)        rect+=vector1[i]*vector2[i];    return rect;} //两个矩阵相乘double** matrixProduct(double **matrix1,int rows1,int columns1,double **matrix2,int columns2){    double **rect=getMatrix(rows1,columns2);    int i,j;    for(i=0;i<rows1;++i){        for(j=0;j<columns2;++j){            double *vec1=getRow(matrix1,rows1,columns1,i);            double *vec2=getColumn(matrix2,columns1,columns2,j);            rect[i][j]=vectorProduct(vec1,vec2,columns1);            free(vec1);            free(vec2);        }    }    return rect;}//矩阵和一个数相乘double** dotProduct(double **matrix,int rows,int columns,double a){double **rect=getMatrix(rows,columns);    int i,j;    for(i=0;i<rows;++i){        for(j=0;j<columns;++j){        rect[i][j]=matrix[i][j]*a;        }    }    return rect;}//两个矩阵相加double** matrixAdd(double **matrix1,double **matrix2,int rows,int columns){double **rect=getMatrix(rows,columns);    int i,j;    for(i=0;i<rows;++i){        for(j=0;j<columns;++j){        rect[i][j]=matrix1[i][j]+matrix2[i][j];        }    }    return rect;} //得到某一列元素的平方和double getColumnNorm(double** matrix,int rows,int columns,int index){    assert(index<columns);    double* vector=getColumn(matrix,rows,columns,index);    double norm=vectorProduct(vector,vector,rows);    free(vector);    return norm;} //打印向量void printVector(double* vector,int len){    int i;    for(i=0;i<len;++i)        printf("%-15.8f\t",vector[i]);    printf("\n");}   #endif

power.c

#include"matrix.h"#include<math.h>#define ROW 6#define ITERATION 50#define EPSILON 0.000002//找出矩阵元素绝对值最大者double getMaxEle(double **matrix,int rows,int cols){double rect=0;int i,j;for(i=0;i<rows;++i){for(j=0;j<cols;++j){if(fabs(matrix[i][j])>rect)rect=fabs(matrix[i][j]);}}return rect;}int main(){//给矩阵A赋值double **A=getMatrix(ROW,ROW);double A1[ROW*ROW]={87,270,-12,-49,-276,40,-14,-45,6,10,46,-4,-50,-156,4,25,162,-25,94,294,-5,-47,-306,49,1,1,3,1,0,2,16,48,1,-6,-48,8};getFromArray(A,ROW,ROW,A1);//取初始Xdouble **X=getMatrix(ROW,1);double X0[ROW]={1,1,1,1,1,1};getFromArray(X,ROW,1,X0);//初始化cdouble c=0;//初始化Ydouble **Y=getMatrix(ROW,1);//开始迭代int iteration=0;while(iteration++<ITERATION){Y=matrixProduct(A,ROW,ROW,X,1);c=getMaxEle(X,ROW,1);assert(c>0);double **newX=dotProduct(Y,ROW,1,1/c);int i;//计算前后两次X的差值double epsilon=0.0;for(i=0;i<ROW;++i){epsilon+=(newX[i][0]-X[i][0])*(newX[i][0]-X[i][0]);}freeMatrix(X,ROW);X=copyMatrix(newX,ROW,1);freeMatrix(newX,ROW);if(sqrt(epsilon)<EPSILON){break;}}printf("Iteration: %d\n",iteration-1);printf("Dominant Eigenvalue=%f\n",c);printf("Dominant Eigenvector=[%f",X[0][0]/c);int i;for(i=1;i<ROW;i++)printf(",%f",X[i][0]/c);printf("]\n");freeMatrix(A,ROW);freeMatrix(X,ROW);freeMatrix(Y,ROW);return 0;}

用位移反幂法求下列矩阵在a=4.2附近的特征值及对应的特征向量。

取最大迭代次数为50，收敛误差为0.000002。

GS.h (高斯-赛德尔迭代)

#ifndef _GS_H#define _GS_H#include"matrix.h"#include<math.h>double** Gauss_Seidel(double **A,int row,double **B){double **X=getMatrix(row,row);int i;for(i=0;i<row;++i)X[i][0]=1;int iteration=0;while(iteration++<50){double **newX=getMatrix(row,row);int i,j;for(i=0;i<row;++i){double sub1=0;double sub2=0;for(j=0;j<i;++j){sub1+=A[i][j]*newX[j][0];}for(j=i+1;j<row;++j)sub2+=A[i][j]*X[j][0];assert(A[i][i]!=0);newX[i][0]=(B[i][0]-sub1-sub2)/A[i][i];}//计算前后两次X的差值double epsilon=0.0;for(i=0;i<row;++i){epsilon+=(newX[i][0]-X[i][0])*(newX[i][0]-X[i][0]);}freeMatrix(X,row);X=copyMatrix(newX,row,1);freeMatrix(newX,row);if(sqrt(epsilon)<0.000001){break;}}//printf("Iteration:%d\n",iteration-1);return X;}#endif

InversePower.c

#include"GS.h"#define ROW 3#define ITERATION 50#define EPSILON 0.0000001#define ALPHA 4.2//找出矩阵元素绝对值最大者double getMaxEle(double **matrix,int rows,int cols){double rect=0;int i,j;for(i=0;i<rows;++i){for(j=0;j<cols;++j){if(fabs(matrix[i][j])>rect)rect=fabs(matrix[i][j]);}}return rect;}int main(){int i;//给矩阵A赋值double **A=getMatrix(ROW,ROW);double A1[ROW*ROW]={0,11,-5,-2,17,-7,-4,26,-10};getFromArray(A,ROW,ROW,A1);//取初始Xdouble **X=getMatrix(ROW,1);for(i=0;i<ROW;++i)X[i][0]=1;//初始化cdouble c=0;//初始化Ydouble **Y=getMatrix(ROW,1);//开始迭代int iteration=0;while(iteration++<ITERATION){double **I=getIndentityMatrix(ROW);double **SUB=dotProduct(I,ROW,ROW,-1*ALPHA);int i,j;double **D=matrixAdd(A,SUB,ROW,ROW);freeMatrix(Y,ROW);Y=Gauss_Seidel(D,ROW,X);freeMatrix(I,ROW);freeMatrix(SUB,ROW);freeMatrix(D,ROW);c=getMaxEle(X,ROW,1);assert(c>0);double **newX=dotProduct(Y,ROW,1,1/c);//计算前后两次X的差值double epsilon=0.0;for(i=0;i<ROW;++i){epsilon+=(newX[i][0]-X[i][0])*(newX[i][0]-X[i][0]);}freeMatrix(X,ROW);X=copyMatrix(newX,ROW,1);freeMatrix(newX,ROW);if(sqrt(epsilon)<EPSILON){break;}}printf("Iteration: %d\n",iteration-1);printf("Dominant Eigenvalue=%f\n",1/c+ALPHA);printf("Dominant Eigenvector=[%f",X[0][0]);for(i=1;i<ROW;i++)printf(",%f",X[i][0]);printf("]\n");freeMatrix(A,ROW);freeMatrix(X,ROW);freeMatrix(Y,ROW);return 0;}