幂迭代和逆幂迭代求最大最小特征值

参考链接 https://wenku.baidu.com/view/ee7ecbeca98271fe910ef9fc.html?from=search

幂迭代算法：

逆幂迭代算法：

实际在使用中A可以先进行LU分解

无论是幂迭代，还是逆幂迭代，只能求出最大和最小特征值与对应的特征向量。

具体问题

python实现

import numpy as np
import matplotlib.pyplot as plt

eigenval1_ = []
eigenval2_ = []
eigenvec1_ = []
eigenvec2_ = []
cnt_ = []
x0 = np.array([1 / np.sqrt(2), 1 / np.sqrt(2)])
x0.resize((2, 1))
dis = 10 ** -12


def normalized_power_iteration(A):
    v0 = x0
    distance = 1
    times = 0.0

    while distance > dis:
        v1 = np.dot(A, v0)
        v1 = v1 / np.linalg.norm(v1)
        distance = np.linalg.norm(v1 - v0)
        v0 = v1
        times += 1.0
    eig = np.linalg.norm(np.dot(A, v1)) / np.linalg.norm(v1)
    return (v1, eig, times)


def normalized_inverse_power_iteration(A):
    # 2x2矩阵 进行LU分解，在这里进行手动的LU分解求逆，因为题目找那个说了是一个2x2的矩阵，实际上当矩阵的维数比较高的时候，应该用LU分解然后解线性方程组来求解。
    L = np.array([1, 0, A[1, 0] / A[0, 0], 1])
    L.resize((2, 2))
    LI = np.array([1, 0, -L[1, 0], 1])
    LI.resize((2, 2))

    U = np.array([A[0, 0], A[0, 1], 0, A[1, 1] - A[0, 1] * A[1, 0] / A[0, 0]])
    U.resize((2, 2))
    UI = np.array([1 / U[0, 0], -U[0, 1] / (U[0, 0] * U[1, 1]), 0, 1 / U[1, 1]])
    UI.resize((2, 2))
    lu = np.dot(UI, LI)

    v0 = x0
    distance = 1
    times = 0

    while distance > dis:
        u1 = np.dot(lu, v0)
        m1 = np.linalg.norm(u1)
        v1 = u1 / m1
        distance = np.linalg.norm(v1 - v0)
        v0 = v1
        times += 1
    eig = np.linalg.norm(np.dot(lu, v1)) / np.linalg.norm(v1)
    return (v1, 1 / eig, times)


for i in np.arange(0, n):
    (vec, val, round1) = normalized_power_iteration(As[i, :, :])
    eigenval1_.append(val)
    eigenvec1_.append(vec[0, 0])
    eigenvec1_.append(vec[1, 0])
    cnt_.append(round1)

    (vec, val, round2) = normalized_inverse_power_iteration(As[i, :, :])
    eigenvec2_.append(vec[0, 0])
    eigenvec2_.append(vec[1, 0])
    eigenval2_.append(val)

eigenval1 = np.array(eigenval1_)
eigenvec1 = np.array(eigenvec1_)
eigenvec1.resize((n, 2))

cnt = np.array(cnt_)

eigenval2 = np.array(eigenval2_)
eigenvec2 = np.array(eigenvec2_)
eigenvec2.resize((n, 2))

plt.xlabel("bilv")
plt.ylabel("round times")
plt.title("result")
plt.show()