求矩阵特征值的QR方法

QR方法求矩阵特征值

import numpy as np


def qr(A):
    # 通过householder变换进行QR分解
    n = A.shape[0]
    Q, R = np.eye(n), A
    # 从第0列到第n-2列
    for i in range(n - 1):
        # 构造变换矩阵
        v = R[i:, i]
        e = np.zeros(n - i)
        e[0] = 1.0
        w = (v + np.sign(v[0]) * np.sqrt(np.sum(v ** 2)) * e).reshape(-1, 1)
        beta = 2 / (w.T.dot(w))
        # 这里需要注意的是，当i=0时，I是空的二维数组
        I0 = np.eye(i)
        H0 = np.eye(n - i) - beta * w.dot(w.T)
        H = np.zeros((n, n))
        H[:i, :i] = I0
        H[i:, i:] = H0
        R = H.dot(R)
        # 直接置零
        R[i + 1:, i] = 0
        Q = Q.dot(H)
    return Q, R


def eigvalues(A, max_iter, eps1, eps2):
    """
    QR方法计算矩阵A的全部特征值
    :param A: 待求矩阵
    :param max_iter: 最大迭代次数
    :param eps1: QR方法迭代停止条件
    :param eps2: 用于判断是否有复数特征值
    :return: 矩阵A的所有特征值
    """
    n = A.shape[0]
    A0 = A
    # 进行迭代
    for _ in range(max_iter):
        Q, R = qr(A0)
        A1 = R.dot(Q)
        if np.linalg.norm(np.diag(A1 - A0), ord=np.inf) < eps1:
            A0 = A1
            break
        A0 = A1
    # 计算特征值
    i = 0
    vals = []
    while i < n:
        # 判断是否存在复数特征值
        if i < n - 1 and np.abs(A0[i + 1, i]) > eps2:
            tmp = np.linalg.eigvals(A0[i:i + 2, i:i + 2])
            vals.extend([tmp[0], tmp[1]])
            i += 2
        else:
            vals.append(A0[i, i])
            i += 1
    return np.array(vals)


if __name__ == '__main__':
    A = np.array(
        [[3, 2, 3, 4],
         [11, 1, 2, 3],
         [2, 8, 9, 1],
         [-4, 2, 9, 11]]
    )
    print(eigvalues(A, max_iter=100, eps1=1e-5, eps2=1e-3))
    # 用numpy验证结果的正确性
    print(np.linalg.eigvals(A))

测试结果：

# 自己编写的算法
[17.59688698+0.j          4.11231656+4.93916298j  4.11231656-4.93916298j
 -1.82152011+0.j        ]
# numpy自带的算法
[-1.82152011+0.j          4.11231656+4.93916298j  4.11231656-4.93916298j
 17.59688698+0.j        ]