首页 > 其他分享> > 数值分析——幂法和逆幂法（Power Method & Inverse Power Method）

数值分析——幂法和逆幂法（Power Method & Inverse Power Method）

2022-01-28 12:34:07 作者：互联网

本系列整理自博主21年秋季学期本科课程数值分析I 的编程作业，内容相对基础，参考书: David Kincaid, Ward Cheney - Numerical Analysis Mathematics of Scientific Computing (2002, Americal Mathematical Society)

算法引入

幂法（Power Method）

辅助部分

归一化幂法

逆幂法（Inverse Power Method）

算法引入

考虑矩阵A有如下两条性质：

（i）有唯一一个模最大的特征值。

（ii）n个特征向量线性无关。

设u1，u2，…，un为对应特征向量， $Au_i=\lambda_iu_i$ ，由线性无关性知{u1,u2,…,un}组成 $\mathbb{C}^n$ 的一组基，任意向量 $x^{(0)}\in \mathbb{C}^n$ ，可表示为 $x^{(0)}=a_1u_1+a_2u_2+...+a_n u_n$

引入迭代： $x^{(k+1)}=Ax^{(k)}$

不妨取 $x^{(0)}=u_1+u_2+...+ u_n$ ，则

$\begin{aligned} x^{(k)}&=A^ku_1+A^ku_2+...+A^ku_n\\&=\lambda_1^ku_1+\lambda_2^ku_2+...+\lambda_n^ku_n\\&=\lambda_1^k[u_1+(\frac{\lambda_2}{\lambda_1})^ku_2+...+(\frac{\lambda_n}{\lambda_1})^ku_n] \end{aligned}$

由于 $(\lambda_j/\lambda_1)^k\rightarrow 0,as\ k\rightarrow \infty$ ，所以 $x^{(k)}=\lambda_1^k[u_1+\varepsilon ^{(k)}]$

设φ为任一线性函数，则 $\varphi (x^{(k)})=\lambda_1^k[\varphi(u_1)+\varphi(\varepsilon ^{(k)})]$ ，可知比值收敛于λ1：

$\frac{\varphi (x^{(k+1)})}{\varphi (x^{(k)})}=\lambda_1\frac{\varphi(u_1)+\varphi(\varepsilon ^{(k+1)})}{\varphi(u_1)+\varphi(\varepsilon ^{(k)})}\rightarrow \lambda_1$

由此计算矩阵A模最大的特征值λ1，这种方法称为幂法

幂法（Power Method）

这里取向量的无穷范数

辅助部分

#include<iostream>
#include<stdio.h>
#include<math.h>
#include<iomanip>
#include<string.h>
using namespace std;
#define M 60 //maximum iteration step
double InftyNorm(double *x,int n) {
	double max = fabs(x[0]);
	for (int i = 1; i < n; i++) {
		if (fabs(x[i]) > max)
			max = fabs(x[i]);
	}
	return max;
}
void multi(double** A, double* x, double* y, int n) {
	//y=Ax
	for (int i = 0; i < n; i++) {
		y[i] = 0;
		for (int j = 0; j < n; j++) {
			y[i] += A[i][j] * x[j];
		}
	}
}
void output(int k, double* x, double r, int n) {
	cout << setw(3) << "k=" <<setw(2)<< k << "   " << "x(" <<setw(2)<< k << ")=(";
	for (int i = 0; i < n; i++) {
		cout <<setiosflags(ios::fixed)<<setprecision(9)<< setw(12) << x[i];
		if (i < n - 1)
			cout << ",";
	}
	cout << ")";
	cout << "   r(" << k - 1 << ")=" << r << endl;
}

归一化幂法

void NormalizedPowerMethod(double** A, double* x,int p, int n) {
	double *y, r = 0, temp = 0, norm = 1;
	y = (double*)malloc(n * sizeof(double));
	cout << "  Normalized Power Method: " << endl;
	cout << setw(3) << "k=" << setw(2) << 0 << "   " << "x(" << setw(2) << 0 << ")=(";
	for (int i = 0; i < n; i++) {
		cout << setiosflags(ios::fixed) << setprecision(9) << setw(13) << x[i];
		if (i < n - 1)
			cout << ",";
	}
	cout << ")" << endl;
	for (int k = 1; k < M; k++) {
		multi(A, x, y, n);
		if (x[p] == 0) {
			cout << "g(x)=0";//g is a linear function, here g(x)=x_p
			return;
		}
		r = y[p] / x[p];
		norm = InftyNorm(y, n);
		if (norm == 0) {
			cout << "||y||=0";
			return;
		}
		for (int i = 0; i < n; i++)
			x[i] = y[i] / norm;
		output(k, x, r, n);
		if (k > 1)
			if(fabs(temp - r) < pow(10, -10))
				break;
		temp = r;
	}
	cout << endl;
	cout << " The leading eigenvalue of A is " << r << ", and the corresponding eigenvector is (";
	for (int i = 0; i < n; i++) {
		cout << setiosflags(ios::fixed) << setprecision(5) << setw(8) << x[i];
		if (i < n - 1)
			cout << ",";
	}
	cout << ")^T" << endl;
}

逆幂法（Inverse Power Method）

引入

定理若λ为A的特征值，且A非奇异，则 $\lambda^{-1}$ 是 $A^{-1}$ 的特征值

结合幂法，可得出计算A的模最小特征值的方法。

假设A的特征值满足：

$|\lambda_1|\geqslant |\lambda_2|\geqslant ...\geqslant |\lambda_{n-1}|> |\lambda_n|>0$

则可知A非奇异， $A^{-1}$ 的特征值为 $\lambda_{j}^{-1}$ ，满足：

$|\lambda_n^{-1}|> |\lambda_{n-1}^{-1}|\geqslant ...\geqslant |\lambda_{1}^{-1}| >0$

计算A模最小特征值问题，化为计算 $A^{-1}$ 模最大特征值的方法，对 $A^{-1}$ 使用幂法即可，迭代形式：

$x^{(k+1)}=A^{-1}x^{(k)}$

注意：由于矩阵A的维数一般较大，不容易对A求逆，采用对A做LU分解的方法。即由 $x^{(k)}$ 计算 $x^{(k+1)}$ 问题，化为求解线性方程组 $Ax^{(k+1)}=LUx^{(k+1)}=x^{(k)}$ 的问题，用前述向前迭代算法（forward subsititution or forward iteration）和向后迭代算法（backward subsititution）求解 $x^{(k+1)}$

代码

辅助部分

void Print(double** A, int n) {//print the matrix A with dimension n
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			cout <<setw(8)<< A[i][j];
			if (j == n - 1) {
				cout <<endl;
			}	
		}
	}
	cout << endl;
}
void DoolittleFactorization(double** A, double** L, double** U, int n) {
	//apply Doolittle factorization on A, LU=A
	for (int k = 0; k < n; k++) {
		if (k == 0) {
			for (int j = 0; j < n; j++) {
				U[0][j] = A[0][j];//first row of U equals to A's first row
				L[j][0] = A[j][0] / U[0][0];//L's first column
			}
		}
		else {
			L[k][k] = 1;
			for (int i = k; i < n; i++) {
				for (int j = 0; j < k; j++) {
					U[k][i] += L[k][j] * U[j][i];
					if (i + 1 < n) {
						L[i + 1][k] += L[i + 1][j] * U[j][k];
					}
				}
				U[k][i] = (A[k][i] - U[k][i]) / L[k][k];
				if (i + 1 < n) {
					L[i + 1][k] = (A[i + 1][k] - L[i + 1][k]) / U[k][k];
				}
			}
		}
	}
	cout << "  apply Doolitte factorization on A：" << endl;
	cout <<setiosflags(ios::fixed) << setprecision(4)<< " L：" << endl;
	Print(L, n);
	cout << endl;
	cout << " U：" << endl;
	Print(U, n);
	cout << endl;
}
void Solve(double** L, double** U, double* b, double* x, int n) {
	//solve linear equations LUx=b
	double* z;
	z = (double*)malloc(n * sizeof(double));
	for (int i = 0; i < n; i++) {//forward subsititution
		z[i] = 0;
		for (int j = 0; j < i; j++) {
			z[i] += L[i][j] * z[j];
		}
		z[i] = (b[i] - z[i]) / L[i][i];
	}
	for (int i = n - 1; i >= 0; i--) {//backward subsititution
		x[i] = 0;
		for (int j = i + 1; j < n; j++) {
			x[i] += U[i][j] * x[j];
		}
		x[i] = (z[i] - x[i]) / U[i][i];
	}
}

逆幂法

}
void NormalizedInversePowerMethod(double** A, double* x, int p, int n) {
	double** L, ** U;
	L = (double**)malloc(n * sizeof(double*));
	U = (double**)malloc(n * sizeof(double*));
	for (int i = 0; i < n; i++){
		L[i]= (double*)malloc(n * sizeof(double));
		U[i] = (double*)malloc(n * sizeof(double));
		for (int j = 0; j < n; j++) {
			L[i][j] = 0;
			U[i][j] = 0;
		}
	}
	DoolittleFactorization(A, L, U, n);
	double *y, r = 0, temp = 0, norm = 1;
	y = (double*)malloc(n * sizeof(double));
	for (int i = 0; i < n; i++)
		y[i] = 0;
	cout << "  Normalized Inverse Power Method: " << endl;
	cout << setiosflags(ios::right);
	cout << setw(3) << "k=" << setw(2) << 0 << "   " << "x(" << setw(2) << 0 << ")=(";
	for (int i = 0; i < n; i++) {
		cout << setiosflags(ios::fixed) << setprecision(9) << setw(13) << x[i];
		if (i < n - 1)
			cout << ",";
	}
	cout << ")" << endl;

	for (int k = 1; k < M; k++) {
		Solve(L, U, x, y, n);//LUy=x
		if (x[p] == 0) {
			cout << "g(x)=0";//g is a linear function, here g(x)=x_p
			return;
		}
		r = y[p] / x[p];
		norm = InftyNorm(y, n);
		if (norm == 0) {
			cout << "||y||=0";
			return;
		}
		for (int i = 0; i < n; i++)
			x[i] = y[i] / norm;
		output(k, x, r, n);
		if (k > 1)
			if (fabs(temp - r) < pow(10, -10)) 
				break;
		temp = r;
	}
	cout << endl;
	cout << " The minimum modulus eigenvalue of A is " << 1 / r << ", and the corresponding eigenvector is (";
	for (int i = 0; i < n; i++) {
		cout << setiosflags(ios::fixed) << setprecision(6) << setw(8) << x[i];
		if (i < n - 1)
			cout << ",";
	}
	cout << ")^T" << endl;
}

标签：5Clambda,Power,int,double,幂法,7C%,20%,Method
来源： https://blog.csdn.net/lk32767/article/details/122698586