特征向量的数值求法

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

Power Method幂法

Definition 如果λ1是矩阵A的所有特征值中绝对值最大的那一个,则称λ1是A的主特征值;与λ1对应的特征向量v1是A的主特征向量。

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

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

Theorem(Power Method) 设A是n×n的方阵,有n个不同的特征值,且|λ1| > |λ2| >= |λ3| >= ... >= |λn|。选择一个合适的X0,序列和{ck}由下列递归式产生:

其中

经过多次迭代后Xk趋于主特征向量V1,ck趋于主特征值λ1

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

Speed of Convergence收敛加速

幂法的收敛速度取决于,也就是说它的收敛速度是线性的。Aitken 加速法可用于任何线性收敛的序列当中,它采用的加速方式是:

用Aitken来加速我们的幂法,Xk的调整公式为:

Shifted-Inverse Power Method平移反幂法

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

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

Theorem (Inverse Eigenvalues) 假设λ和V是A的一个特征对,,那么,V是A-1的一个特征对。

Theorem (Shifted-Inverse Eigenvalues) 假设λ和V是A的一个特征对,,那么,V是的一个特征对。

Theorem (Shifted-Inverse Power Method) 假设A是一个n×n的矩阵,有n个互不相同的特征值λ1,λ2,...,λn。对于其中之一个特征值λj,可以选择一个常数a,使得的主特征值。进一步地,如果选择一个合适的X0,序列和{ck}可由下列递归式产生:

 

其中

 经过多次迭代后Xk趋于的主特征向量Vj,Xk同时也是A的主特征向量,ck趋于的主特征值u1。最终我们可以求出A的主特征值:

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

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

注意aii即系数矩阵主对角线上的元素不能有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);
	//取初始X
	double **X=getMatrix(ROW,1);
	double X0[ROW]={1,1,1,1,1,1};
	getFromArray(X,ROW,1,X0);
	//初始化c
	double c=0;
	//初始化Y
	double **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);
	//取初始X
	double **X=getMatrix(ROW,1);
	for(i=0;i<ROW;++i)
		X[i][0]=1;
	//初始化c
	double c=0;
	//初始化Y
	double **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;
}
posted @ 2012-08-17 14:44  高性能golang  阅读(17294)  评论(0编辑  收藏  举报