幂法

含义

幂法是计算矩阵的按模最大的特征值和相应特征向量的一种向量迭代法

适用于大型稀疏矩阵

基础知识

求矩阵A特征值问题等价于求A的特征方程。

∣λE−A∣=[λ−a11−a12⋯−a1n−a21λ−a22⋯−a2n⋯⋯⋯⋯λ−an1−an2⋯λ−ann]=0|\lambda E-A| = \left[ \begin{matrix} \lambda-a_{11} & -a_{12} & \cdots & -a_{1n} \\ -a_{21} & \lambda-a_{22} & \cdots & -a_{2n} \\ \cdots & \cdots & \cdots & \cdots \\ \lambda-a_{n1} & -a_{n2} & \cdots & \lambda-a_{nn} \\ \end{matrix} \right] = 0∣λE−A∣=⎣⎢⎢⎡λ−a11−a21⋯λ−an1−a12λ−a22⋯−an2⋯⋯⋯⋯−a1n−a2n⋯λ−ann⎦⎥⎥⎤=0
求A的属于特征值\lambda的特征向量等价求
(λE−A)x=0(\lambda E -A)x =0(λE−A)x=0 的非零解
设λ为A∈Rx∗n的特征值，x称为A的与特征值λ相对应的一个特征向量，即Ax=λx,（x≠0）,则有：设\lambda 为A\in R^{x*n}的特征值，x称为A的与特征值\lambda 相对应的一个特征向量，即Ax = \lambda x,（x\neq0）,则有：设λ为A∈Rx∗n的特征值，x称为A的与特征值λ相对应的一个特征向量，即Ax=λx,（x̸=0）,则有：
（１）cx（c≠0）也是A的与特征值λ相对应的一个特征向量,即Acx=λcxcx（c\neq0）也是A的与特征值\lambda 相对应的一个特征向量,即Acx = \lambda cxcx（c̸=0）也是A的与特征值λ相对应的一个特征向量,即Acx=λcx
（２）λk为Ak的特征值，Akx=λkx\lambda ^k为A^k的特征值，A^kx=\lambda ^kxλk为Ak的特征值，Akx=λkx

思想

设n阶矩阵A的特征值和特征向量为
λ1,λ2,...,λnv1,v2,...vn\lambda_1,\lambda_2,...,\lambda_n　　　 v_1,v_2,...v_nλ1,λ2,...,λn　　　v1,v2,...vn
v1,v2,...vn线性无关，且v_1,v_2,...v_n线性无关，且v1,v2,...vn线性无关，且
∣λ1∣>∣λ2∣>...>∣λn∣|\lambda_1|>|\lambda_2|>...>|\lambda_n|∣λ1∣>∣λ2∣>...>∣λn∣
称λ1为A的按模最大特征值（主特征值）称\lambda_1为A的按模最大特征值（主特征值）称λ1为A的按模最大特征值（主特征值）

任取非零向量x0=(x1(0),xx(0),⋯ ,xn(0))x_0 =(x_1^{(0)},x_x^{(0)},\cdots,x_n^{(0)})x0=(x1(0),xx(0),⋯,xn(0))
特征向量构成n维空间的一组基.任取的初始向量X(0)由它们的线性组合给出
x0=a1v1+a2v2+…+anvnx_0 =a_1v_1+a_2v_2+…+a_nv_nx0=a1v1+a2v2+…+anvn
假设a1≠0,由此构造向量序列xk=Axk−1，k=1,2,3...假设a_1\neq0,由此构造向量序列x_k = Ax_{k-1}，k=1,2,3...假设a1̸=0,由此构造向量序列xk=Axk−1，k=1,2,3...
其中，x1=Ax0=a1Av1+a2Av2+...+anAvn=x_1=Ax_0=a_1Av_1+a_2Av_2+...+a_nAv_n=x1=Ax0=a1Av1+a2Av2+...+anAvn=
＝a1λ1v1+a2λ2v2+...+anλnvn＝a_1\lambda_1v_1+a_2\lambda_2v_2+...+a_n\lambda_nv_n　　　　　　　　　　　　　　　　　　　　　＝a1λ1v1+a2λ2v2+...+anλnvn

x2=Ax1=a1λ1Av1+a2λ2Av2+...+anλnAvn=x_2=Ax_1=a_1\lambda_1Av_1+a_2\lambda_2Av_2+...+a_n\lambda_nAv_n=x2=Ax1=a1λ1Av1+a2λ2Av2+...+anλnAvn=
＝a1λ12v1+a2λ22v2+...+anλn2vn＝a_1\lambda_1^2v_1+a_2\lambda_2^2v_2+...+a_n\lambda_n^2v_n　　　　　　　　　　　　　　　　　　　　　＝a1λ12v1+a2λ22v2+...+anλn2vn

得到向量序列：
x1＝a1λ1v1+a2λ2v2+...+anλnvnx_1＝a_1\lambda_1v_1+a_2\lambda_2v_2+...+a_n\lambda_nv_nx1＝a1λ1v1+a2λ2v2+...+anλnvn
x2=a1λ12v1+a2λ22v2+...+anλn2vnx_2= a_1\lambda_1^2v_1+a_2\lambda_2^2v_2+...+a_n\lambda_n^2v_nx2=a1λ12v1+a2λ22v2+...+anλn2vn
…\dots　　　　　　…
xk=a1λ1kv1+a2λ2kv2+...+anλnkvnx_k= a_1\lambda_1^kv_1+a_2\lambda_2^kv_2+...+a_n\lambda_n^kv_nxk=a1λ1kv1+a2λ2kv2+...+anλnkvn

下面按模最大特征值λ1是单根的情况讨论:
xk=a1λ1kv1+a2λ2kv2+...+anλnkvn=λ1k[a1v1+a2(λ2λ1)kv2+...+an(λnλ1)2vn]=λ1k(a1v1+εk)x_k= a_1\lambda_1^kv_1+a_2\lambda_2^kv_2+...+a_n\lambda_n^kv_n=\lambda_1^k[a_1v_1 +a_2(\frac{\lambda_2}{\lambda_1})^kv_2+...+a_n(\frac{\lambda_n}{\lambda_1})^2v_n]=\lambda_1^k(a_1v_1+\varepsilon_k)xk=a1λ1kv1+a2λ2kv2+...+anλnkvn=λ1k[a1v1+a2(λ1λ2)kv2+...+an(λ1λn)2vn]=λ1k(a1v1+εk)

εk=a2(λ2λ1)kv2+...+an(λnλ1)2vn\varepsilon_k=a_2(\frac{\lambda_2}{\lambda_1})^kv_2+...+a_n(\frac{\lambda_n}{\lambda_1})^2v_nεk=a2(λ1λ2)kv2+...+an(λ1λn)2vn
由于∣λ1∣>∣λ2∣>...>∣λn∣故∣λjλ1∣<1(j=2,3,...n),故：|\lambda_1|>|\lambda_2|>...>|\lambda_n|　　故|\frac{\lambda_j}{\lambda_1}|<1(j=2,3,...n),故：∣λ1∣>∣λ2∣>...>∣λn∣　　故∣λ1λj∣<1(j=2,3,...n),故：
lim⁡k→+∞εk=0则有lim⁡k→+∞1λ1kxk=lim⁡k→+∞λ1k(a1v1+εk)λ1k=a1v1\lim_{k\rightarrow+\infty}\varepsilon_k=0　　则有\lim_{k\rightarrow+\infty}\frac{1}{\lambda_1^k}x_k=\lim_{k\rightarrow+\infty}\frac{\lambda_1^k(a_1v_1+\varepsilon_k)}{\lambda_1^k}=a_1v_1k→+∞limεk=0　　则有k→+∞limλ1k1xk=k→+∞limλ1kλ1k(a1v1+εk)=a1v1

lim⁡k→+∞1λ1kxk=a1v1就是λ1对应的特征向量\lim_{k\rightarrow+\infty}\frac{1}{\lambda_1^k}x_k=a_1v_1 就是\lambda_1对应的特征向量k→+∞limλ1k1xk=a1v1就是λ1对应的特征向量
用向量xk+1的第i个分量与向量xk的第i个分量之比的极限：用向量x_{k+1}的第i个分量与向量x_k的第i个分量之比的极限：用向量xk+1的第i个分量与向量xk的第i个分量之比的极限：
lim⁡k→+∞(xk+1)i(xk)i=lim⁡k→+∞λ1k+1(a1v1+εk+1)iλ1k(a1v1+εk)i=lim⁡k→+∞λ1[a1(v1)i+(εk+1)i]a1(v1)i+(εk)i=λ1\lim_{k\rightarrow+\infty}\frac{(x_{k+1})_i}{(x_k)_i}=\lim_{k\rightarrow+\infty}\frac{\lambda_1^{k+1}(a_1v_1+\varepsilon_{k+1})_i}{\lambda_1^k(a_1v_1+\varepsilon_k)_i}=\lim_{k\rightarrow+\infty}\frac{\lambda_1[a_1(v_1)_i+(\varepsilon_{k+1})_i]}{a_1(v_1)_i+(\varepsilon_k)_i}=\lambda_1k→+∞lim(xk)i(xk+1)i=k→+∞limλ1k(a1v1+εk)iλ1k+1(a1v1+εk+1)i=k→+∞lima1(v1)i+(εk)iλ1[a1(v1)i+(εk+1)i]=λ1

得到两个极限：
lim⁡k→+∞1λ1kxk=a1v1lim⁡k→+∞(xk+1)i(xk)i=λ1\lim_{k\rightarrow+\infty}\frac{1}{\lambda_1^{k}}x_k=a_1v_1　　　\lim_{k\rightarrow+\infty}\frac{(x_{k+1})_i}{(x_k)_i}=\lambda_1k→+∞limλ1k1xk=a1v1　　　k→+∞lim(xk)i(xk+1)i=λ1

因此当k趋近无穷时，就能近似得到按模最大的特征值和对应的特征向量：
λ1≈(xk+1)i(xk)ixk≈λ1ka1v1\lambda_1 \approx \frac{(x_{k+1})_i}{(x_k)_i} 　　　x_k\approx\lambda_1^k a_1v_1λ1≈(xk)i(xk+1)i　　　xk≈λ1ka1v1

因此得幂方法的迭代公式：xk=Axk−1k=1,2,...x_k = Ax_{k-1}　　k=1,2,...xk=Axk−1　　k=1,2,...
当k充分大时，有：
{λ1≈(xk+1)i(xk)ixk≈λ1ka1v1\begin{cases} \lambda_1 \approx \frac{(x_{k+1})_i}{(x_k)_i} \\ x_k\approx\lambda_1^k a_1v_1 \end{cases}{λ1≈(xk)i(xk+1)ixk≈λ1ka1v1

存在的问题是：当∣λ1∣>1时，xk中不为0的分量会随着k→∞而趋于∞，计算机计算时会造成溢出，当|\lambda_1|>1时，x_k中不为0的分量会随着k\rightarrow \infty而趋于 \infty，计算机计算时会造成溢出，当∣λ1∣>1时，xk中不为0的分量会随着k→∞而趋于∞，计算机计算时会造成溢出，当|\lambda_1|<1时,x_k中的分量会随着k\rightarrow \infty而趋于0$，

实际计算时，为了避免计算过程中出现绝对值过大或过小的数参加运算，通常在每步迭代时，将向量“归一化”即用的按模最大的分量。
修改计算公式：
{yk=Axk−1mk=max(yk)(k=1,2,...)xk=yk/mk\begin{cases} y_k=Ax_{k-1}\\ m_k=max(y_k) & (k=1,2,...)\\ x_k=y_k/m_k \end{cases}⎩⎪⎨⎪⎧yk=Axk−1mk=max(yk)xk=yk/mk(k=1,2,...)
上式中mk是向量yk中绝对值最大的一个分量，这时xk分量的模最大为1m_k是向量y_k中绝对值最大的一个分量，这时x_k分量的模最大为1mk是向量yk中绝对值最大的一个分量，这时xk分量的模最大为1
当k充分大时，有：
{λ1≈mkv1≈xk\begin{cases} \lambda_1 \approx m_k \\ v_1 \approx x_k \end{cases}{λ1≈mkv1≈xk

实例

求矩阵
[2321034361]的主特征值和其对应的特征向量\left[ \begin{matrix} 2 &3 & 2 \\ 10 & 3 & 4 \\ 3 & 6 & 1 \end{matrix} \right] 的主特征值和其对应的特征向量 ⎣⎡2103336241⎦⎤的主特征值和其对应的特征向量
根据迭代公式：
{yk=Axk−1mk=max(yk)(k=1,2,...)xk=yk/mk\begin{cases} y_k=Ax_{k-1}\\ m_k=max(y_k) & (k=1,2,...)\\ x_k=y_k/m_k \end{cases}⎩⎪⎨⎪⎧yk=Axk−1mk=max(yk)xk=yk/mk(k=1,2,...)
取x0=(0,0,1)T,则有：x_0=(0,0,1)^T,则有：x0=(0,0,1)T,则有：
y1=(2,4,1)T，m1=4，x1=1m1y1=(0.5,1,0.25)Ty_1 = (2,4,1)^T，m_1 = 4，x_1=\frac{1}{m_1}y_1=(0.5,1,0.25)^Ty1=(2,4,1)T，m1=4，x1=m11y1=(0.5,1,0.25)T
直到k=8

Python实现

#-*- coding:utf-8 -*-
#@author：whs
#@time: 2019/3/2111:25import numpy as npdef Solve(mat, max_itrs, min_delta):"""mat 表示矩阵max_itrs 表示最大迭代次数min_delta 表示停止迭代阈值"""itrs_num = 0delta = float('inf')N = np.shape(mat)[0]# 所有分量都为1的列向量x = np.ones(shape=(N, 1))#x = np.array([[0],[0],[1]])while itrs_num < max_itrs and delta > min_delta:itrs_num += 1y = np.dot(mat, x)#print(y)m = y.max()#print("m={0}".format(m))x = y / mprint("***********第{}次迭代*************".format(itrs_num))print("y = ",y)print("m={0}".format(m))print("x^T为：",x)print("===================================")IOS = np.array([[2, 3, 2],[10, 3, 4],[3, 6, 1]], dtype=float)
Solve(IOS, 10, 1e-10)

输出：

***********第1次迭代*************
y =  [[ 7.][17.][10.]]
m=17.0
x^T为： [[0.41176471][1.        ][0.58823529]]
===================================
***********第2次迭代*************
y =  [[5.        ][9.47058824][7.82352941]]
m=9.470588235294118
x^T为： [[0.52795031][1.        ][0.82608696]]
===================================
***********第3次迭代*************
y =  [[ 5.70807453][11.58385093][ 8.40993789]]
m=11.58385093167702
x^T为： [[0.49276139][1.        ][0.72600536]]
===================================
***********第4次迭代*************
y =  [[ 5.43753351][10.83163539][ 8.20428954]]
m=10.831635388739945
x^T为： [[0.50200485][1.        ][0.75743775]]
===================================
***********第5次迭代*************
y =  [[ 5.5188852 ][11.04979951][ 8.2634523 ]]
m=11.049799514875502
x^T为： [[0.49945569][1.        ][0.74783731]]
===================================
***********第6次迭代*************
y =  [[ 5.49458599][10.98590609][ 8.24620437]]
m=10.985906091381928
x^T为： [[0.50014864][1.        ][0.75061668]]
===================================
***********第7次迭代*************
y =  [[ 5.50153064][11.00395312][ 8.2510626 ]]
m=11.003953118800313
x^T为： [[0.49995948][1.        ][0.74982713]]
===================================
***********第8次迭代*************
y =  [[ 5.49957322][10.99890329][ 8.24970556]]
m=10.998903290037243
x^T为： [[0.50001105][1.        ][0.75004801]]
===================================
***********第9次迭代*************
y =  [[ 5.50011813][11.00030258][ 8.25008117]]
m=11.000302582696586
x^T为： [[0.49999699][1.        ][0.74998675]]
===================================
***********第10次迭代*************
y =  [[ 5.49996747][10.99991685][ 8.24997771]]
m=10.999916852432536
x^T为： [[0.50000082][1.        ][0.75000364]]
===================================

幂法求解矩阵特征值及特征向量相关推荐

使用MTL库求解矩阵特征值和特征向量
关于矩阵的特征值和特征向量求解,大部分的数学运算库都进行了提供,下面是使用MTL库的接口进行封装. #include <mtl/matrix.h> #include <mtl/mtl ...
2021-01-07 matlab数值分析矩阵特征值与特征向量的计算改进乘幂法反幂法
matlab数值分析矩阵特征值与特征向量的计算 1改进乘幂法 function [t,y]=eigIPower(A,v0,ep) [tv,ti]=max(abs(v0)); lam0=v0(ti) ...
矩阵特征值和特征向量详细计算过程(转载)
1.矩阵特征值和特征向量定义 A为n阶矩阵,若数λ和n维非0列向量x满足Ax=λx,那么数λ称为A的特征值,x称为A的对应于特征值λ的特征向量.式Ax=λx也可写成( A-λE)x=0,并且|λE-A ...
求矩阵特征值和特征向量
求矩阵特征值和特征向量的一个小程序代码较长,如果不能执行,就是要建立结构体,大家试试吧,希望能用. // // 实对称三对角阵的全部特征值与特征向量的计算 // // 参数: // 1. doubl ...
基于QR分解迭代求解方阵特征值和特征向量
基于QR分解迭代求解方阵特征值和特征向量一.特征值与特征向量求解的难点线性代数的知识告诉我们如果要求一个方阵的特征值,只需要求解如下的特征方程的根即可: f ( λ ) = ( λ − λ 1 ) ...
python numpy逆_Python使用numpy计算矩阵特征值、特征向量与逆矩阵
原标题:Python使用numpy计算矩阵特征值.特征向量与逆矩阵 Python扩展库numpy.linalg的eig()函数可以用来计算矩阵的特征值与特征向量,而numpy.linalg.inv() ...
Python使用numpy计算矩阵特征值、特征向量与逆矩阵
Python扩展库numpy.linalg的eig()函数可以用来计算矩阵的特征值与特征向量,而numpy.linalg.inv()函数用来计算可逆矩阵的逆矩阵. >>> impor ...
pythonnumpy库求特征向量_Python使用numpy计算矩阵特征值、特征向量与逆矩阵
Python扩展库numpy.linalg的eig()函数可以用来计算矩阵的特征值与特征向量,而numpy.linalg.inv()函数用来计算可逆矩阵的逆矩阵. >>> impor ...
QR分解求矩阵特征值、特征向量 C语言
最近在看一个高光谱图像压缩算法,其中涉及到正交变换,计算正交变换时,需要对普通矩阵求其特征向量.想要在网上找一个现成的程序,可能是我百度的能力不强吧,居然真的没找见.好了废话不多说,下面进入正题. 计 ...
机器学习中的数学基础：（1.1）矩阵特征值和特征向量的几何意义
给定一个二维矩阵先求出该矩阵的特征值与特征向量,特征值分别获是:, 对应的特征向量为: (列向量)PS:此处的U是正交矩阵,根据正交矩阵的性质,可以有如果从定义来理解特征向量的化,某一物体经过该矩 ...

幂法求解矩阵特征值及特征向量

文章目录

幂法

含义

基础知识

思想

实例

Python实现

幂法求解矩阵特征值及特征向量相关推荐

最新文章

热门文章