冪法（指數迭代法）

　　冪法是通過迭代來計算矩陣的主特征值（按模最大的特征值）與其對應特征向量的方法，適合於用於大型稀疏矩陣。

基本定義

　　設$A = (a_{ij})\in R^{n\times n}$，其特征值為$\lambda_i$，對應特征向量$x_i(i=1,...,n)$，即$Ax_i = \lambda_i x_i(i=1,...,n)$，且$\{x_1,...,x_n\}$線性無關。

　　任取一個非零向量$v_0\in R^{n}$，且$v_0\ne 0$，構造一個關於矩陣$A$的乘冪的向量序列：

$v_k = Av_{k-1}=A^2v_{k-2}=A^3v_{k-3}=...=A^kv_{0}$

　　稱$v_k$為迭代向量。

　　設特征值$\lambda_i$的前$r$個為絕對值最大的特征值（ppt中分為$\lambda_1$強占優和非強占優，感覺沒必要），即有：

$|\lambda_1|=|\lambda_2|=...=|\lambda_r|>|\lambda_{r+1}|\ge...\ge|\lambda_n|$

　　由於$\{x_1,...,x_n\}$ 線性無關，所以構成$R^n$的一個基，於是$v_0$能被表達為：

$v_0=\sum\limits_{i=1}^{n}\alpha_i x_i$（且設$\alpha_1...\alpha_r$非全零）

 　　由$Ax_i = \lambda_i x_i$：

$v_k = Av_{k-1}=...=A^kv_{0}=\sum\limits_{i=1}^{n}A^k\alpha_i x_i=\sum\limits_{i=1}^{n}\lambda_i^k\alpha_i x_i=\lambda_1^k(\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_k)$

　　其中：

$\varepsilon_k = \sum\limits_{i=r+1}^{n}(\frac{\lambda_i}{\lambda_1})^k\alpha_ix_i$ 

　　因為$\lambda_1$最大，所以有$|\frac{\lambda_i}{\lambda_1}|<1　(i=r+1,...,n)$，從而有：

$\lim\limits_{k\to \infty}(\frac{\lambda_i}{\lambda_1})^k=0　(i=r+1,...,n)$

　　所以有：

$\lim\limits_{k\to \infty}\varepsilon_k = 0$

$\lim\limits_{k\to \infty}v_k = \lim\limits_{k\to \infty}\lambda_1^k(\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_k) = \lim\limits_{k\to \infty}\lambda_1^k(\sum\limits_{i=1}^{r}\alpha_ix_i)$

　　因為在上式中$(\sum\limits_{i=1}^{r}\alpha_ix_i)$是固定項，可以看出，迭代到后期，$v_{k+1}$和$v_k$的各個元素有固定比值$\lambda_1$，即：

$\lim\limits_{k\to \infty}\frac{(v_{k+1})_i}{(v_{k})_i} = \lambda_1$

　　這樣，收斂到主特征值后，還可另外計算它對應的一個特征向量（其實就是構成$v_0$的前$r$項之和，而且只能算一個）：

$\lim\limits_{k\to \infty}\frac{v_{k}}{\lambda_1^k} = \sum\limits_{i=1}^{r}\alpha_ix_i$

　　其中收斂速度由比值$|\frac{\lambda_{r+1}}{\lambda_1}|$決定，越小收斂越快。

冪法改進

　　以上定義中，隨着不斷的迭代，如果$|\lambda_1|>1$，$v_k$會無限大；如果$|\lambda_1|<1$，$v_k$則會趨於0。在計算機中就會導致溢出而出錯。當然如果$|\frac{\lambda_{r+1}}{\lambda_1}|$很小，使得$\varepsilon_k$能快速收斂為0，並且$|\lambda_1|$不會太大或太小以至於在$\varepsilon_k$收斂前$v_k$就溢出了，那用上述方法就行了。

　　但是通常不會碰到那么好的運氣，所以就在每次迭代的時候就要對$v_k$進行“規范化”，防止它溢出。規范化就是將向量除以一個標量，使向量的長度為1，而向量的長度有多種度量法，比如$||\cdot||_2、||\cdot||_\infty$等。一般使用$||\cdot||_\infty$，就是向量各個分量的絕對值的最大值（用別的也一樣，但是這個復雜度最小，計算機最好算）。

　　為了便於理解，將一次迭代分為兩步，一步乘矩陣$A$，一步除以向量長度進行規范化：

迭代序號	乘$A$（$v_i$）	規范化（$u_i$）
初始化	隨機取$v_0=\sum\limits_{i=0}^{n}\alpha_ix_i$	$u_0=\frac{v_0}{||v_0||}$
1	$v_1=Au_0=\frac{Av_0}{\left\|v_0\right\|}$	$u_1=\frac{v_1}{||v_1||}=\frac{\frac{Av_0}{\left\|v_0\right\|}}{\left\|\frac{Av_0}{\left\|v_0\right\|}\right\|}=\frac{Av_0}{\left\|Av_0\right\|}$
2	$v_2=Au_1=\frac{A^2v_0}{\left\|Av_0\right\|}$	$u_2=\frac{v_2}{||v_2||}=\frac{\frac{A^2v_0}{\left\|Av_0\right\|}}{\left\|\frac{A^2v_0}{\left\|Av_0\right\|}\right\|}=\frac{A^2v_0}{\left\|A^2v_0\right\|}$
...	...	...
k	$v_k=Au_{k-1}=\frac{A^kv_0}{\left\|A^{k-1}v_0\right\|}$	$u_k=\frac{v_k}{||v_k||}=\frac{\frac{A^kv_0}{\left\|A^{k-1}v_0\right\|}}{\left\|\frac{A^kv_0}{\left\|A^{k-1}v_0\right\|}\right\|}=\frac{A^kv_0}{\left\|A^kv_0\right\|}$
...	...	...

　　根據迭代，以下計算主特征值和對應的特征向量。

主特征值　　

　　對於向量序列${v_k}$：

$v_k=\frac{A^kv_0}{\left\|A^{k-1}v_0\right\|}=\frac{\lambda_1^k(\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_k)}{\left\|\lambda_1^{k-1}(\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_{k-1})\right\|}$

　　對$v_k$取范數：

$||v_k||=\frac{\left\|\lambda_1^k(\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_k)\right\|}{\left\|\lambda_1^{k-1}(\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_{k-1})\right\|}=|\lambda_1|\frac{\left\|\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_k\right\|}{\left\|\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_{k-1}\right\|}\to\left|\lambda_1\right|, (k\to\infty)$

　　注意！在計算機中並不是直接通過這個極限$k\to\infty$來計算，而是通過上面表格的迭代來實現$k\to\infty$，由於每一步的“規范化”，才能使得計算成為可能。而且如果直接用這個極限算的話就和前面的沒區別了，“規范化”的優勢沒有用起來。

　　另外，因為算范數的時候是先取向量分量的絕對值再算的，所以改進后的方法計算出來的主特征值不可避免帶有絕對值，因此特征值的符號還要再判斷一下（下面判斷）。而改進之前的是不帶的（因為沒有計算范數進行規范化），如果符合條件，改進之前的也可一試。

主特征值對應的一個特征向量

　　觀察向量$u_k$：

$u_k=\frac{A^kv_0}{\left\|A^kv_0\right\|}=\frac{\lambda_1^k(\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_k)}{\left\|\lambda_1^{k}(\sum\limits_{i=1}^{r}\alpha_ix_i+\varepsilon_{k})\right\|}\to{\rm sign}(\lambda_1^k)\frac{\sum\limits_{i=1}^{r}\alpha_ix_i}{\left\|\sum\limits_{i=1}^{r}\alpha_ix_i\right\|},(k\to\infty)$

　　如上，通過$u_k$可獲得一個單位化的主特征值特征向量。因為還乘上了$\lambda_1^k$的符號，所以向量各個分量的符號和迭代的次數$k$有關。當$k$足夠大（使迭代收斂），連續兩次迭代，$u_k$的分量的符號不變時，則可以判斷主特征值符號為正，否則為負，從而解決了上面主特征值符號未知的問題。

　　另外，求出的單位特征向量（如果用L2范數）是主特征值的所有不相關特征向量的線性組合，與設置的$v_0$有關。當然如果主特征值的沒有重根，且$\alpha_1$不為零（為零的話就迭代出第二大的特征值），那規范化后得出的都是同一個特征向量而和$v_0$無關（如果特征值為負，符號可能隨迭代次數改變）。

Python代碼

import numpy as np

#矩陣A
A = np.matrix(
    [[-5.6,2.,0.],
    [3.,-1.,1.],
    [0.,1.,3.]])
L0=2#范數類型
v = np.matrix([[-2.],[-7.],[1.]])#v0
u = v/(np.linalg.norm(v,L0))#u0

#迭代函數
def iterate_fun(A,u,final,L):
    i=0
    while i<final: 
        v = A*u
        u = v/(np.linalg.norm(v,L))
        i=i+1
    print("冪法特征值：")
    print(np.linalg.norm(v,L))
    print("numpy特征值：")
    print(np.linalg.eig(A)[0])
    print("冪法特征向量：")
    print(u)
    print("numpy特征向量：")
    print(np.linalg.eig(A)[1])

iterate_fun(A,u,1010,L0)

　　結果截圖：

　　主特征值與numpy算出來絕對值一致。這里用的是L2范數，且主特征值無重根，所以計算出來的特征向量單位化后和numpy主特征值對應的向量除符號外相同。

迭代加速（原點平移法）

　　冪法迭代的速度由比值$\frac{\left|\lambda\right|_{second\,biggest}}{|\lambda_1|}$決定，當比值比較接近1時，迭代的速度就可能比較慢了。迭代加速的方法是給矩陣$A$減去一個數量陣$pI$：（$I$為單位陣）

$B = A-pI$

　　假設$A$的特征值為：

$\lambda_1=\lambda_2=...=\lambda_r>...>\lambda_n$

　　則$B$的特征值就是：

$\lambda_1-p=\lambda_2-p=...=\lambda_r-p>...>\lambda_n-p$

　　而$A,B$的特征向量都是相同的。

　　注意！這里與上面假設不同，這里沒有絕對值了。那么即使經過數量陣相加，特征值的大小關系還是不變的。用$B$迭代出特征值后再加回$p$即是$A$的特征值。

　　對於$B$，我們要：

　　1、保持$\lambda_1-p$的絕對值在$B$中依舊是最大，對$B$進行迭代才能算出最大特征值。

　　2、由於我們想要減小$B$的迭代比值：

$\frac{|\omega|}{|\lambda_1-p|},\omega=\max\limits_k(|\lambda_{k}-p|),(k=r+1,...,n)$

　　實際上，我們可以推斷$\omega$只會取為$\lambda_{r+1}-p$或是$\lambda_{n}-p$（畫個數軸就能判斷），而當：

$p=\frac{\lambda_{second\,biggest}+\lambda_n}{2}=\frac{\lambda_{r+1}+\lambda_n}{2}$

　　數軸上$p$在$\lambda_{r+1}和\lambda_{n}$中間：（圖中假設$r=1$）

　　有$\lambda_{r+1}-p = -(\lambda_{n}-p)$，此時$B$的迭代比值最小，為：

$\left|\frac{\lambda_{r+1}-\lambda_n}{2\lambda_{1}-\lambda_{r+1}-\lambda_{n}}\right|$

　　迭代速度達到最快，而$p$再小或再大都會使迭代速度減慢。  

　　另外，當$p=\frac{\lambda_1+\lambda_{second\,smallest}}{2}$，還可以求$B$（或$A$）的最小特征值。實際上，通過平移，冪法可以求的就是在數軸兩端的特征值。

　　在實際應用中，$A$的特征值並不知道，所以$p$是不能精確計算的，這個方法告訴我們，當發現收斂速度很慢時，我們可以適當地變動一下加速收斂。

加速python代碼

　　代碼加速迭代計算最小特征值（負數，絕對值最大），加速前迭代55次，加速后迭代25次，速度提升約50%。

import numpy as np

p = (-0.20106557+3.27764013)/2#p的大小

#矩陣A
A = np.matrix(
    [[-5.6,2.,0.],
    [3.,-1.,1.],
    [0.,1.,3.]])
A = A-np.eye(3)*p #平移
L0=2#范數類型 
v = np.matrix([[-2.],[-7.],[1.]])#v0
u = v/(np.linalg.norm(v,L0))#u0

#迭代函數
def iterate_fun(A,v,u,L):
    i=0
    former_v = -1
    after_v = -1
    isPositive = 1
    while True: 
        i=i+1 
        former_v = np.linalg.norm(v,L)
        v = A*u
        after_v = np.linalg.norm(v,L)
        if former_v-after_v==0.:
            if u[0]*v[0] < 0:
                isPositive = -1 
            break
        u = v/(np.linalg.norm(v,L))
    print("迭代次數：")
    print(i)
    print("冪法特征值：")
    print(isPositive*np.linalg.norm(v,L)+p)
    print("numpy特征值：")
    print(np.linalg.eig(A+np.eye(3)*p)[0])
    print("冪法特征向量：")
    print(u)
    print("numpy特征向量：")
    print(np.linalg.eig(A+np.eye(3)*p)[1])

iterate_fun(A,v,u,L0)

　　結果截圖：