-對於正定的對稱矩陣，奇異值等於特征值，奇異向量等於特征向量。在這種情況下用奇異值分解就把特征值和特征向量求出來了。但是只要是方陣，它就有特征值和特征向量，對於一般的方陣，特征值和特征向量怎么求呢（當然我指的是數值求法）？這就要用本文即將介紹的“冪法”。
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平移反冪法
使用位移反冪法首先需要提供一個好的起始點，這個點要接近一個特征向量，然后我們的位移反冪法才能夠以更高的精度算出這個特征向量。QM和Given'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;
}