矩阵特征值计算的python实现(1)——幂法与反幂法

\qquad 已知 $n\times n$ 矩阵 $A=[a_{ij}]$，矩阵 $A$ 的特征多项式定义为：

$\phi (\lambda )=\det (\lambda \mathbf{I}-A)=\left|\begin{array}{cccc}\lambda -a_{11}& -a_{12}& \cdots & -a_{1n}\\ -a_{21}& \lambda -a_{22}& \cdots & -a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ -a_{n1}& -a_{n2}& \cdots & \lambda -a_{nn}\end{array}\right|$

\qquad 矩阵 $A$ 的特征方程定义为：

$\phi (\lambda )=\det (\lambda \mathbf{I}-A)=0$

$(1)$ 矩阵 $A$ 的特征方程的解，即为特征值$\text{(Eigenvalue)}$。

\qquad 　　矩阵 $A$ 的所有特征值的集合，记为 $\lambda (A)$。

$(2)$ 若 $\lambda$ 为矩阵 $A$ 的特征值，齐次方程 $(\lambda \mathbf{I}-A)\mathbf{x}=0$ 的非零解 $\mathbf{x}$ 为矩阵 $A$ 的对应于特征值 $\lambda$ 的特征向量$\text{(Eigenvector)}$。
\qquad

1. 矩阵特征值界的估计

\qquad 考虑 $n\times n$ 矩阵 $A=[a_{ij}]$，令：

$\{\begin{array}{l}{r}_i=\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n|a_{ij}|,i=1,2,\cdots ,n\\ \\ D_i=\left\{z||z-a_{ij}|\le r_i,z\in \mathbf{C}\right\}\end{array}$

\qquad 称复平面上以 $a_{ii}$ 为圆心、$r_i$ 为半径的所有圆盘为矩阵 $A$ 的$\text{Gerschgorin}$圆盘。

\qquad
$\text{Gerschgorin圆盘定理}$

$(1)$ 一个 $n\times n$ 矩阵 $A=[a_{ij}]$ 的每一个特征值必定属于下述的某个圆盘中：

$|\lambda -a_{ii}|\le r_i=\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n|a_{ij}|,i=1,2,\cdots ,n$

\qquad 　或者说，矩阵 $A$ 的特征值都在复平面上 $n$ 个圆盘的并集中。

$(2)$ 如果矩阵 $A$ 的 $m$ 个圆盘组成了一个连通的并集 $S$，且 $S$ 与余下 $n-m$ 个圆盘是分离的，那么 $S$ 内恰好包含矩阵 $A$ 的 $m$ 个特征值。

\qquad 　特别地，如果矩阵 $A$ 的一个圆盘 $D_i$ 与其它圆盘是分离的（孤立圆盘），那么圆盘 $D_i$ 中精确地包含矩阵 $A$ 的一个特征值。
\qquad
$\text{例}$
$A=\left[\begin{array}{ccc}4& 1& 0\\ 1& 0& -1\\ 1& 1& -4\end{array}\right]⟹\{\begin{array}{l}{D}_1=|\lambda -4|\le 1\\ D_2=|\lambda |\le 2\\ D_3=|\lambda +4|\le 2\end{array}$

\qquad 由于 $D_1$ 是孤立圆盘，可得到： $3\le {\lambda }_1\le 5$
\qquad 矩阵 $A$ 的三个特征值位于三个圆盘的并集中，即：${\lambda }_i\in D_1\cup D_2\cup D_3$
\qquad

2. 幂法与反幂法

2.1 幂法

\qquad 假设 $n\times n$ 矩阵 $A=[a_{ij}]$ 的特征值 ${\lambda }_i$ 所对应的特征向量为 ${\mathbf{x}}_i$，且特征值的排列满足：

$|{\lambda }_1|>|{\lambda }_2|\ge |{\lambda }_3|\ge \cdots \ge |{\lambda }_n|$

\qquad 其中，${\lambda }_1=\underset{i}{\max }|{\lambda }_i|$ 为主特征值，且为实数。

\qquad 幂法用来计算矩阵的主特征值 ${\lambda }_1$ 及对应特征向量 ${\mathbf{x}}_1$，是一种迭代的方法。
\qquad

• 幂法的基本思路

\qquad 任取一个非零的初始向量 ${\mathbf{v}}^{(0)}$，构造出一个向量序列 { v ( k ) } \{\boldsymbol v^{(k)}\} {v(k)}，也就是：

$\{\begin{array}{l}{\mathbf{v}}^{(1)}=A{\mathbf{v}}^{(0)}\\ {\mathbf{v}}^{(2)}=A{\mathbf{v}}^{(1)}=A^2{\mathbf{v}}^{(0)}\\ \cdots \\ {\mathbf{v}}^{(k+1)}=A{\mathbf{v}}^{(k)}=A^{k+1}{\mathbf{v}}^{(0)}\\ \cdots \end{array}$

\qquad 假设矩阵 $A$ 的 $n$ 个特征向量 { x i } \{\boldsymbol x_i\} {xi​} 线性无关，初始向量 ${\mathbf{v}}^{(0)}$ 就可以表示为：

${\mathbf{v}}^{(0)}={\alpha }_1{\mathbf{x}}_1+{\alpha }_2{\mathbf{x}}_2+\cdots +{\alpha }_n{\mathbf{x}}_n$

\qquad 那么

$\begin{array}{rl}{\mathbf{v}}^{(k)}& =A^k{\mathbf{v}}^{(0)}={\alpha }_1{A}^k{\mathbf{x}}_1+{\alpha }_2{A}^k{\mathbf{x}}_2+\cdots +{\alpha }_n{A}^k{\mathbf{x}}_n\\ & ={\alpha }_1{\lambda }_1^k{\mathbf{x}}_1+{\alpha }_2{\lambda }_2^k{\mathbf{x}}_2+\cdots +{\alpha }_n{\lambda }_n^k{\mathbf{x}}_n\\ & ={\lambda }_1^k\left[{\alpha }_1{\mathbf{x}}_1+\sum _{i=2}^n{\alpha }_i{\left(\frac{{\lambda }_i}{{\lambda }_1}\right)}^k{\mathbf{x}}_i\right]\end{array}$

$(1)$ 若记 ${\mathbf{\epsilon }}_k=\sum _{i=2}^n{\alpha }_i{\left(\frac{{\lambda }_i}{{\lambda }_1}\right)}^k{\mathbf{x}}_i$，那么：${\mathbf{v}}^{(k)}={\lambda }_1^k\left[{\alpha }_1{\mathbf{x}}_1+{\mathbf{\epsilon }}_k\right]$

$(2)$ 由于 $|{\lambda }_1|>|{\lambda }_i|⟹\left|\frac{{\lambda }_i}{{\lambda }_1}\right|<1(i=2,\cdots ,n)$
$⟹\underset{k\to \infty }{\lim }{\mathbf{\epsilon }}_k=0$　（显然，$\left|\frac{{\lambda }_2}{{\lambda }_1}\right|$ 越小，收敛越快）
$⟹\underset{k\to \infty }{\lim }\frac{{\mathbf{v}}^{(k)}}{{\lambda }_1^k}={\alpha }_1{\mathbf{x}}_1$

\qquad\qquad 因此，迭代过程 ${\mathbf{v}}^{(k+1)}=A{\mathbf{v}}^{(k)}$ 持续到 $k$ 充分大时，${\mathbf{v}}^{(k)}\approx {\alpha }_1{\lambda }_1^k{\mathbf{x}}_1$
\qquad
$(3)$ 若 ${\mathbf{v}}_i^k$ 表示 ${\mathbf{v}}^{(k)}$ 的第 $i$ 个分量，那么

$\frac{{\mathbf{v}}_i^{k+1}}{{\mathbf{v}}_i^k}={\lambda }_1\frac{{\alpha }_1{\mathbf{x}}_{1i}+{\mathbf{\epsilon }}_{k+1,i}}{{\alpha }_1{\mathbf{x}}_{1i}+{\mathbf{\epsilon }}_{k,i}}⟹\underset{k\to \infty }{\lim }\frac{{\mathbf{v}}_i^{k+1}}{{\mathbf{v}}_i^k}={\lambda }_1$

\qquad\qquad\qquad 或者记为：$\frac{{\mathbf{v}}_i^{k+1}}{{\mathbf{v}}_i^k}\to {\lambda }_1$　（收敛速度由$\left|\frac{{\lambda }_2}{{\lambda }_1}\right|$ 确定）

\qquad\qquad 这就说明，相邻两个迭代向量的分量的比值收敛到主特征值。

\qquad

• 主特征值的计算

\qquad 任取非零的初始向量 ${\mathbf{v}}^{(0)}={\mathbf{u}}^{(0)}\ne 0$，构造向量序列

{ v ( 1 ) = A u ( 0 ) = A v ( 0 ) , u ( 1 ) = v ( 1 ) max ⁡ { v ( 1 ) } v ( 2 ) = A u ( 1 ) = A 2 v ( 0 ) max ⁡ { A v ( 0 ) } , u ( 2 ) = v ( 2 ) max ⁡ { v ( 2 ) } ⋯ ⋯ v ( k ) = A u ( k − 1 ) = A k v ( 0 ) max ⁡ { A k − 1 v ( 0 ) } , u ( k ) = v ( k ) max ⁡ { v ( k ) } \qquad\qquad\left\{\begin{array}{lc}\textcolor{red}{\boldsymbol v^{(1)}=A\boldsymbol u^{(0)}}=A\boldsymbol v^{(0)},&\boldsymbol u^{(1)}=\dfrac{\boldsymbol v^{(1)}}{\max\{\boldsymbol v^{(1)}\}}\\ \\ \textcolor{red}{\boldsymbol v^{(2)}=A\boldsymbol u^{(1)}}=\dfrac{A^2\boldsymbol v^{(0)}}{\max\{A\boldsymbol v^{(0)}\}},&\boldsymbol u^{(2)}=\dfrac{\boldsymbol v^{(2)}}{\max\{\boldsymbol v^{(2)}\}}\\ \\ \qquad\qquad\cdots&\cdots \\ \\ \textcolor{red}{\boldsymbol v^{(k)}=A\boldsymbol u^{(k-1)}}=\dfrac{A^k\boldsymbol v^{(0)}}{\max\{A^{k-1}\boldsymbol v^{(0)}\}},&\textcolor{mediumblue}{\boldsymbol u^{(k)}=\dfrac{\boldsymbol v^{(k)}}{\max\{\boldsymbol v^{(k)}\}}} \end{array}\right. ⎩ ⎨ ⎧​v(1)=Au(0)=Av(0),v(2)=Au(1)=max{Av(0)}A2v(0)​,⋯v(k)=Au(k−1)=max{Ak−1v(0)}Akv(0)​,​u(1)=max{v(1)}v(1)​u(2)=max{v(2)}v(2)​⋯u(k)=max{v(k)}v(k)​​

\qquad
\qquad 由于 ${\mathbf{v}}^{(k)}\approx {\alpha }_1{\lambda }_1^k{\mathbf{x}}_1$，如果 $|{\lambda }_1|>1$（或者 $|{\lambda }_1|<1$），随着 $k$ 值越来越大，${\lambda }_1^k$ 会趋于无穷（或趋于 $0$），从而造成计算过程溢出。

\qquad 因此，需要每次迭代时，需要将迭代向量 ${\mathbf{v}}^{(k)}$ 规范化，也就是：

u ( k ) = v ( k ) max ⁡ { v ( k ) } \qquad\qquad\textcolor{mediumblue}{\boldsymbol u^{(k)}=\dfrac{\boldsymbol v^{(k)}}{\max\{\boldsymbol v^{(k)}\}}} u(k)=max{v(k)}v(k)​

此处， max ⁡ { v ( k ) } \max\{\boldsymbol v^{(k)}\} max{v(k)} 表示迭代向量 ${\mathbf{v}}^{(k)}$ 具有最大绝对值的分量

\qquad
幂法求主特征值

步骤 $1:$ 令 $k=0$，${\mathbf{v}}^{(0)}={\mathbf{u}}^{(0)}\ne 0$

步骤 $2:$ 计算 ${\mathbf{v}}^{(k)}=A{\mathbf{u}}^{(k-1)}$， μ ( k ) = max ⁡ { v ( k ) } \mu^{(k)}=\max\{\boldsymbol v^{(k)}\} μ(k)=max{v(k)}

步骤 $3:$ 计算 ${\mathbf{u}}^{(k)}=\frac{{\mathbf{v}}^{(k)}}{{\mu }^{(k)}}$（当 $k\to \infty$，那么 ${\mu }^{(k)}\to {\lambda }_1$）

重复执行步骤 $(2)$ 和 $(3)$，直到 $|{\mu }^{(k)}-{\mu }^{(k-1)}|<\epsilon$，得到主特征值的估计 ${\hat{\lambda }}_1={\mu }^{(k)}$

\qquad
$\text{例1}$：用幂法求 $A=\left[\begin{array}{ccc}1& 1& 0.5\\ 1& 1& 0.25\\ 0.5& 0.25& 2\end{array}\right]$ 的主特征值及对应的特征向量。

import numpy as np
def eig_power(A,v0,eps):
    uk = v0
    flag = 1
    val_old = 0
    n = 0
    while flag:
        n = n+1
        vk = A*uk        
        val = vk[np.argmax(np.abs(vk))]        
        uk = vk/val              
        if (np.abs(val-val_old)<eps):
            flag = 0
        val_old = val
        print(np.asarray(uk).flatten(),val)
    print('max eigenvalue:',val)
    print('eigenvector:',np.asarray(uk).flatten())
    print('iteration:',n)
    return val, uk  
    
if __name__ == '__main__':
    A = np.matrix([[1,   1,      0.5],
                   [1,   1,      0.25],
                   [0.5,   0.25,   2.0]], dtype='float')
    v0 = np.matrix([[1],[1],[1]], dtype='float')
    eps = 1e-10
    val,uk = eig_power(A,v0,eps)

计算结果：

max eigenvalue: 2.536526
eigenvector: [0.74822115 0.64966114 1.        ]
iteration: 41

2.2 幂法的加速方法

$\text{1. 原点平移法}$

$(1)$ 当 $\left|\frac{{\lambda }_2}{{\lambda }_1}\right|$ 的值接近于 $1$，也就是 ${\lambda }_1$ 和 ${\lambda }_2$ 的值比较接近，就会导致 $\frac{{\mathbf{v}}_i^{k+1}}{{\mathbf{v}}_i^k}\to {\lambda }_1$ 的收敛速度变慢。

$(2)$ 为了加快收敛速度，可以人为地拉大 ${\lambda }_1$ 和 ${\lambda }_2$ 的距离，也就是让 $\left|\frac{{\lambda }_2}{{\lambda }_1}\right|$ 的值能够尽可能地偏离 $1$，达到加快收敛速度的目的。

\qquad 为此，引入矩阵 $B=A-p\mathbf{I}$，可得到矩阵 $B$ 的特征值为 ${\lambda }_1-p,{\lambda }_2-p,\cdots ,{\lambda }_n-p$。如果能够选择合适的 $p$ 值，使得：

$\left|\frac{{\lambda }_2-p}{{\lambda }_1-p}\right|<\left|\frac{{\lambda }_2}{{\lambda }_1}\right|$

\qquad 此时，对矩阵 $B$ 应用幂法求出 ${\lambda }_1-p$ 的速度，会比直接用矩阵 $A$ 求 ${\lambda }_1$ 更快。因此，选择一个合适的 $p$ 值非常关键。

\qquad
$\text{例2}$：用原点平移法求 $A=\left[\begin{array}{ccc}1& 1& 0.5\\ 1& 1& 0.25\\ 0.5& 0.25& 2\end{array}\right]$ 的主特征值及对应的特征向量。

import numpy as np
def eig_power(A,v0,eps,p):
    uk = v0
    flag = 1
    val_old = 0
    n = 0
    A = A-p*np.eye(len(A))
    while flag:
        n = n+1
        vk = A*uk        
        val = vk[np.argmax(np.abs(vk))]        
        uk = vk/val       
        if (np.abs(val-val_old)<eps):
            flag = 0
        val_old = val
        print(np.asarray(uk).flatten(),val)
    print('max eigenvalue:',val+p)
    print('eigenvector:',np.asarray(uk).flatten())
    print('iteration:',n)
    return val+p, uk  
    
if __name__ == '__main__':
    A = np.matrix([[1,   1,      0.5],
                   [1,   1,      0.25],
                   [0.5,   0.25,   2.0]], dtype='float')
    v0 = np.matrix([[1],[1],[1]], dtype='float')
    eps = 1e-10
    val,uk = eig_power(A,v0,eps,0.75)

计算结果：

max eigenvalue: 2.536526
eigenvector: [0.74822115 0.64966114 1.        ]
iteration: 26			对比例1结果，直接采用幂法需要迭代41次

\qquad
$\text{2. Rayleigh商加速法}$

\qquad 假设 $n\times n$ 对称矩阵 $A$ 的特征值满足：$|{\lambda }_1|>|{\lambda }_2|\ge |{\lambda }_3|\ge \cdots \ge |{\lambda }_n|$，规范化向量 ${\mathbf{u}}^{(k)}$ 的瑞利商可以给出主特征值 ${\lambda }_1$ 更好的近似：

$\frac{(A{\mathbf{u}}^{(k)},{\mathbf{u}}^{(k)})}{({\mathbf{u}}^{(k)},{\mathbf{u}}^{(k)})}={\lambda }_1+O\left({\left(\frac{{\lambda }_2}{{\lambda }_1}\right)}^{2k}\right)$

\qquad
$\text{例3}$：用$\text{Rayleigh}$商加速法求 $A=\left[\begin{array}{ccc}1& 1& 0.5\\ 1& 1& 0.25\\ 0.5& 0.25& 2\end{array}\right]$ 的主特征值及对应的特征向量。

import numpy as np
def eig_power(A,v0,eps):
    uk = v0
    flag = 1
    val_old = 0
    n = 0
    while flag:
        n = n+1
        vk = A*uk
        t = ((A*uk).T*uk)/(uk.T*uk)        
        val = t[0,0]   
        uk = vk/val     
        if (np.abs(val-val_old)<eps):
            flag = 0
        val_old = val   
        print(np.asarray(uk).flatten(),val)
    print('max eigenvalue:',val)
    print('eigenvector:',np.asarray(uk).flatten())
    print('iteration:',n)
    return val, uk 
    
if __name__ == '__main__':
    A = np.matrix([[1,   1,      0.5],
                   [1,   1,      0.25],
                   [0.5,   0.25,   2.0]], dtype='float')
    v0 = np.matrix([[1],[1],[1]], dtype='float')
    eps = 1e-10
    val,uk = eig_power(A,v0,eps)

计算结果：

max eigenvalue: 2.5365258603724308
eigenvector: [0.92496879 0.80312773 1.23621526]
iteration: 20		对比例1和例2结果，直接采用幂法需要迭代41次，p=0.75原点加速法需要26次

2.3 反幂法

\qquad 与幂法的思路一样，反幂法可以计算最小模值的特征值及其特征向量。

\qquad 仍然假设 $n\times n$ 矩阵 $A=[a_{ij}]$ 的特征值 ${\lambda }_i$ 所对应的特征向量为 ${\mathbf{x}}_i$，且特征值的排列满足：

$|{\lambda }_1|\ge |{\lambda }_2|\ge |{\lambda }_3|\ge \cdots \ge |{\lambda }_n|>0$

\qquad 那么，$A^{-1}$ 的特征值满足：

$\left|\frac{1}{{\lambda }_n}\right|\ge \left|\frac{1}{{\lambda }_{n-1}}\right|\ge \left|\frac{1}{{\lambda }_{n-2}}\right|\ge \cdots \ge \left|\frac{1}{{\lambda }_1}\right|$

\qquad 对矩阵 $A^{-1}$ 使用幂法，也就是“反幂法”，可以计算矩阵 $A^{-1}$ 的主特征值 $\frac{1}{{\lambda }_n}$，也就是矩阵 $A$ 的模值最小的特征值。

\qquad
反幂法求最小模值特征值

步骤 $1:$ 令 $k=0$，${\mathbf{v}}^{(0)}={\mathbf{u}}^{(0)}\ne 0$

步骤 $2:$ 计算 ${\mathbf{v}}^{(k)}=A^{-1}{\mathbf{u}}^{(k-1)}$， μ ( k ) = max ⁡ { v ( k ) } \mu^{(k)}=\max\{\boldsymbol v^{(k)}\} μ(k)=max{v(k)}

\qquad 　 其中，${\mathbf{v}}^{(k)}$ 可通过解线性方程组 $A{\mathbf{v}}^{(k)}={\mathbf{u}}^{(k-1)}$ 求得

反幂法同样可以采用“原点加速法”：
采用原点加速法时， ${\mathbf{v}}^{(k)}=(A-p\mathbf{I}{)}^{-1}{\mathbf{u}}^{(k-1)}$
通过解线性方程组 $(A-p\mathbf{I}){\mathbf{v}}^{(k)}={\mathbf{u}}^{(k-1)}$ 求得 ${\mathbf{v}}^{(k)}$

步骤 $3:$ 计算 ${\mathbf{u}}^{(k)}=\frac{{\mathbf{v}}^{(k)}}{{\mu }^{(k)}}$（当 $k\to \infty$，那么 ${\mu }^{(k)}\to \frac{1}{{\lambda }_n}$）

重复执行步骤 $(2)$ 和 $(3)$，直到 $|{\mu }^{(k)}-{\mu }^{(k-1)}|<\epsilon$，得到最小模值特征值的估计 ${\hat{\lambda }}_n=\frac{1}{{\mu }^{(k)}}$

\qquad
$\text{例4}$：用反幂法求 $A=\left[\begin{array}{ccc}2& 1& 0\\ 1& 3& 1\\ 0& 1& 4\end{array}\right]$ 的对应于 $\lambda =1.2679$ 的特征向量。

import numpy as np   
def eig_invpower(A,v0,eps,p=0):
    uk = v0
    flag = 1
    val_old = 0
    n = 0
    if p!=0:
        A = A-p*np.eye(len(A))
    LU,La,Ua,order0,order1 = Doolittle_pivot(np.asarray(A))   # PA=LU
    while flag:
        n = n+1        
        vk = solveLineq(La,Ua,np.asarray(uk)[order1,:])      
        val = vk[np.argmax(np.abs(vk))]        
        uk = np.asmatrix(vk.reshape(len(A),1))/val        
        print(np.asarray(uk).flatten())
        if (np.abs(1/val-val_old)<eps):
            flag = 0
        val_old = 1/val
    print('min eigenvalue:',1/val+p)
    print('eigenvector:',np.asarray(uk).flatten())
    print('iteration:',n)
    return 1/val+p, uk   

def Doolittle_pivot(A):  # A为np.array，而不是np.matrix
    n = len(A)
    LU = A.copy()
    order1 = np.arange(n)
    for r in range(n):        
        ut = LU[:r,r].reshape(r,1)
        si = A[r:,r] - np.sum(ut*LU[r:,:r].T,axis=0)
        ir = np.argmax(np.abs(si))
        if ir!=0:
            LU[[r,r+ir],:] = LU[[r+ir,r],:]
            order1[[r,r+ir]] = order1[[r+ir,r]]        
        lt = LU[r,:r].reshape(r,1)
        LU[r,r:] = LU[r,r:] - np.sum(lt*LU[:r,r:],axis=0)        
        if r==n-1:
            continue
        LU[r+1:,r] = (LU[r+1:,r] - np.sum(ut*LU[r+1:,:r].T,axis=0))/LU[r,r]       
    U = np.triu(LU)
    L = np.tril(LU) - np.diag(np.diag(LU)) + np.eye(n)
    order0 = []
    [order0.insert(i,np.where(order1==i)[0][0]) for i in range(n)]        
    return LU,L,U,order0,order1

def solveLineq(L,U,b):  # b为np.array，而不是np.matrix
    rows = len(b)
    y = np.zeros(rows)
    y[0] = b[0]/L[0,0]
    for k in range(1,rows):
        y[k] = (b[k] - np.sum(L[k,:k]*y[:k]))/L[k,k]      
    x = np.zeros(rows)
    k = rows-1
    x[k] = y[k]/U[k,k] 
    for k in range(rows-2,-1,-1):
        x[k] = (y[k] - np.sum(x[k+1:]*U[k,k+1:]))/U[k,k]    
    return x    
    
if __name__ == '__main__':
    A = np.matrix([[2,   1,   0],
                   [1,   3,   1],
                   [0,   1,   4]], dtype='float')    
    v0 = np.matrix([[1],[1],[1]], dtype='float')
    eps = 1e-10
    val,uk = eig_invpower(A,v0,eps,1.2679)

计算结果：

A:
 [[0.7321 1.     0.    ]
 [1.     1.7321 1.    ]
 [0.     1.     2.7321]]
L:
 [[ 1.          0.          0.        ]
 [ 0.          1.          0.        ]
 [ 0.7321     -0.26807041  1.        ]]
U:
 [[1.         1.7321     1.        ]
 [0.         1.         2.7321    ]
 [0.         0.         0.00029517]]
L*U:
 [[0.7321 1.     0.    ]
 [1.     1.7321 1.    ]
 [0.     1.     2.7321]]
v1: [ 6776.39884879 -4960.0015972   1815.81991772]
u1: [ 1.         -0.73195243  0.26796237]
v2: [ 20327.46459393 -14880.73682922   5446.72771553]
u2: [ 1.         -0.73205081  0.26794919]
v3: [ 20328.33052096 -14881.37077439   5446.95974656]
u3: [ 1.         -0.73205081  0.26794919]
v4: [ 20328.33054145 -14881.3707894    5446.95975206]
u4: [ 1.         -0.73205081  0.26794919]
min eigenvalue: 1.2679491924311228
eigenvector: [ 1.         -0.73205081  0.26794919]
iteration: 4

【注】：选主元的三角分解求解线性方程组，请参考《解线性方程组的python实现(2)——矩阵三角分解法》第$3$部分。
\qquad