python QR分解求解线性方程组和矩阵本征值和本征向量

下面的代码提供了两个函数

• solve_linear_equ, 利用QR分解求解线性方程组，输入是一个二维的非奇异的系数方阵和一个常数array，输出是该线性方程组的解
• eigen, 输入是一个实方阵(写的时候没有考虑复数的情况，我也不清楚如果有复数的特征值会发生什么事)，输出是有两个，一个是一维的array，对应该矩阵所有的特征值，一个是二维的array，对应相应的特征向量。

import numpy as np


# coefficients of givens matrix
def Givens(alpha: np.float64, beta: np.float64):
    eta = np.sqrt(alpha * alpha + beta * beta)
    if eta != np.float64(0):
        c = alpha / eta
        s = beta / eta
    else:
        c = 1
        s = 0
    return c, s, eta


# modify hess_mat in place, return Q
# Q = G1^T * G2^T * ... * Gn^T
def upperHessenbergQR(hess_mat: np.array):
    size = np.size(hess_mat, 0)
    Q = np.identity(size)
    for i in range(size - 1):
        c, s, eta = Givens(hess_mat[i][i], hess_mat[i + 1][i])
        if eta != 0:
            hess_mat[i][i] = eta
            hess_mat[i + 1][i] = 0
            tmp_1 = c * hess_mat[i, i + 1:] + s * hess_mat[i + 1, i + 1:]
            tmp_2 = -s * hess_mat[i, i + 1:] + c * hess_mat[i + 1, i + 1:]
            hess_mat[i, i + 1:] = tmp_1
            hess_mat[i + 1, i + 1:] = tmp_2

            tmp_3 = np.matmul(Q[:, i:i + 2], np.array([[c, -s], [s, c]]))
            Q[:, i:i + 2] = tmp_3
    return Q


def HouseHolder(u: np.array):
    gamma = np.sqrt(np.sum(u * u))
    w = np.zeros(np.size(u), dtype=np.float64)
    if gamma != 0 and gamma != u[0]:
        detla = np.sqrt(gamma * (gamma - u[0]))
        w[0] = (u[0] - gamma) / detla
        w[1:] = u[1:] / detla
    else:
        gamma = 0
    return gamma, w


def QRHouseHolder(mat: np.array):
    size = np.size(mat, 0)
    Q = np.identity(size, dtype=np.float64)
    for i in range(size - 1):
        gamma, u = HouseHolder(mat[i:, i])
        temp = np.identity(len(u), dtype=np.float64) - np.matmul(u.reshape(-1, 1), u.reshape(1, -1))
        H = np.identity(size, dtype=np.float64)
        H[i:, i:] = temp
        if gamma != 0:
            Q = np.matmul(H, Q)
            v = np.matmul(u, mat[i:, i + 1:])
            mat[i][i] = gamma
            mat[i + 1:, i] = 0
            for j in range(i, size):
                mat[j, i + 1:] = mat[j, i + 1:] - u[j - i] * v
        else:
            mat[i + 1:, i] = 0
    return Q.T


# solve linear equation by QR decomposition
def solve_linear_equ(mat_copy: np.array, constants: np.array):
    R = mat_copy.copy()
    Q = QRHouseHolder(R)
    right = np.matmul(Q.T, constants.reshape(-1, 1)).reshape(-1)
    size = np.size(constants)
    roots = np.zeros(size)
    roots[size - 1] = 1 if R[size - 1][size - 1] == 0 else right[size - 1] / R[size - 1][size - 1]

    for i in range(size - 1, -1, -1):
        tmp = right[i] - np.sum(R[i, i + 1:] * roots[i + 1:])
        roots[i] = 1 if R[i][i] == 0 else tmp / R[i][i]
    return roots


# transform the matrix into hessenberg matrix
def HouseholderReduction(mat: np.array):
    size = np.size(mat, 0)
    for i in range(size - 2):
        gamma, u = HouseHolder(mat[i + 1:, i])
        if gamma != 0:
            v = np.matmul(u, mat[i + 1:, i + 1:])
            mat[i + 1][i] = gamma
            mat[i + 2:, i] = 0
            for j in range(i + 1, size):
                mat[j, i + 1:] = mat[j, i + 1:] - u[j - i - 1] * v
            v = np.matmul(mat[:, i + 1:], u)
            for j in range(i + 1, size):
                mat[:, j] = mat[:, j] - u[j - i - 1] * v
        else:
            mat[i + 2:, i] = 0


"""
def HouseholderReduction2(mat: np.array):
    size = np.size(mat, 0)
    First = 1
    for i in range(size - 2):
        gamma, u = HouseHolder(mat[i + 1:, i])
        if gamma != 0:
            temp = np.identity(len(u)) - u.reshape(-1, 1) * u.reshape(1, -1)
            H = np.identity(size)
            H[i + 1:, i + 1:] = temp
            if First == 0:
                print(np.matmul(H, mat))
            mat = np.matmul(np.matmul(H, mat), H)
            mat[i + 2:, i] = 0
            if First == 0:
                print(mat)
                First = 1
    return mat
"""


def check_converge(mat: np.array, ord: int):
    return np.abs(mat[ord][ord - 1]) < 1.0e-10


def solve_eigen_vec(mat: np.array, size: int, ord: int):
    roots = np.zeros(size, dtype=np.float64)
    roots[ord] = 1
    for inx in range(ord - 1, -1, -1):
        v = np.sum(mat[inx, inx + 1:] * roots[inx + 1:])
        coeff = mat[ord][ord] - mat[inx][inx]
        if coeff == 0:
            roots[inx] = 1
        else:
            roots[inx] = v / coeff
    return roots / np.sqrt(np.sum(roots * roots))


def eigen(m: np.array):
    mat = m.copy()
    size = np.size(mat, 0)
    ord = size - 1
    HouseholderReduction(mat)
    mat_copy = mat.copy()
    cnt = 0
    while ord:
        u = mat[ord][ord] if cnt > 100 else 0
        cnt += 1
        mat -= u * np.identity(size, dtype=np.float64)
        q = upperHessenbergQR(mat)
        mat = np.matmul(mat, q) + u * np.identity(size, dtype=np.float64)
        if check_converge(mat, ord):
            mat[ord][ord - 1] = 0
            ord -= 1
    eigen_value = np.array([mat[i][i] for i in range(size)])
    # print(eigen_value)

    ord = size - 1
    ret_q = np.identity(size)
    while ord:
        mat_copy -= eigen_value[ord] * np.identity(size)
        q = upperHessenbergQR(mat_copy)
        ret_q = np.matmul(ret_q, q)
        mat_copy = np.matmul(mat_copy, q) + eigen_value[ord] * np.identity(size)
        if check_converge(mat_copy, ord):
            mat_copy[ord][ord - 1] = 0
            ord -= 1
    eigen_value = np.array([mat_copy[i][i] for i in range(size)])
    eigen_vec = np.array(
        [np.matmul(ret_q, solve_eigen_vec(mat_copy, size, i).reshape(-1, 1)).reshape(-1) for i in
         range(size)])
    """e, v = np.linalg.eig(mat_copy)
    #print(e)
    #print(v)
    eigen_vec = np.array(
        [solve_linear_equ(mat_copy, size, i) for i in
         range(size)])
    #print(eigen_value)
    #print(eigen_vec)"""
    return eigen_value, eigen_vec


"""
def HouseHolder(mat: np.array):
    cols = mat.shape[1]

    Q = np.eye(cols)
    R = np.copy(mat)

    for col in range(cols - 1):
        a = np.linalg.norm(R[col:, col])
        e = np.zeros((cols - col))
        e[0] = 1.0
        num = R[col:, col] - a * e
        den = np.linalg.norm(num)

        u = num / den
        H = np.eye((cols))
        H[col:, col:] = np.eye((cols - col)) - 2 * u.reshape(-1, 1).dot(u.reshape(1, -1))
        R = H.dot(R)

        Q = Q.dot(H)

    return Q, R


def check_converge(mat: np.array, size: int):
    for i in range(1, size):
        for j in range(i - 1):
            if mat[i][j] > 1.0e-14:
                return False
    return True


def eigen(mat: np.array):
    size = np.size(mat, 0)
    while 1:
        q, r = HouseHolder(mat)
        new_mat = r * q
        if check_converge(new_mat, size):
            print(new_mat)
            return np.array([new_mat[i][i] for i in range(size)]), np.array([])
        mat = new_mat


y = np.array([[2, -1, 0, -1],
              [-1, 2, -1, 0],
              [0, -1, 2, -1],
              [-1, 0, -1, 2]])

h = y[1:, 1:]
print(h)
h[0][0] = 100
"""
"""
c, _ = eigen(y)
print(c)
t, _ = np.linalg.eig(y)
print(t)
"""
"""
y = np.array([[1, 2, 3] for i in range(3)])
q, r = HouseHolder(y)
print(q)
print(r)
print(np.matmul(q, r))
"""
# y = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]], dtype=np.float64)
# H, Q = hessenberg(y, calc_q=True)
# H_copy = H.copy()
# print(H)
# q, r = qr(H)
# mat = HouseholderReduction2(y)
# print(mat)
# print(q)
# print(r)
# print(np.matmul(q, r))
# my_q = upperHessenbergQR(H)
# print(my_q)
# print(q)
# print(my_q)
# print(np.matmul(my_q, H))
if __name__ == '__main__':
    y = np.array([[2, -1, 0, -1],
                  [-1, 2, -1, 0],
                  [0, -1, 2, -1],
                  [-1, 0, -1, 2]], dtype=np.float64)
    y = np.array([[1, 1, 1, -2], [1, 1, 0, 3], [1, 2, 1, 2.7], [-1, 3, 2, 0.9]], dtype=np.float64)
    ret = np.array([1.3, 2.7, 1, -1])
    print(y)
    print(solve_linear_equ(y, ret))
    print(np.matmul(np.linalg.inv(y), ret.reshape(-1, 1)).reshape(-1))

注：

• Householder变换矩阵应为$H=I-2v{v}^T$，代码中取的$||v||=\sqrt{2}$，因此吸收掉了这个2
• 计算本征值的方法基于位移版的QR迭代算法，首先把矩阵相似变换为上Hessenberg矩阵，然后反复QR迭代直到收敛。然后解线性方程组，得到最终的上三角矩阵的特征向量，再乘以Q矩阵得到原本的特征向量。