要求完成课堂上讲的关于矩阵分解的LU、QR（Gram-Schmidt）、Orthogonal Reduction (Householder reduction 和Givens reduction)和 URV程序实现，要求如下：

1、一个综合程序，根据选择参数的不同，实现不同的矩阵分解；在此基础上，实现Ax=b方程组的求解，以及计算A的行列式；
2、可以用matlab、Python等编写程序，需附上简单的程序说明，比如参数代表什么意思，输入什么，输出什么等等，附上相应的例子；
3、一定是可执行文件，例如 .m文件等,不能是word或者txt文档。附上源代码，不能为直接调用matlab等函数库;

import numpy as np
import scipy
import copy


def det_P(P: np.ndarray):
    """
    计算LU分解中的P矩阵的行列式
    :param P: LU分解中的P矩阵
    :return: P矩阵的行列式
    """
    h, w = P.shape
    t = 0
    for i in range(h):
        if P[i][i] != 1:
            t = t + 1
            for j in range(h):
                if P[i][j] == 1:
                    s = j
            for j in range(i + 1, h):
                if P[j][i] == 1:
                    P[i][i] = 1
                    P[i][j] = 0
                    P[j][i] = 0
                    P[j]▼显示 = 1
    if t % 2 == 0:
        return 1
    else:
        return -1


def det_triangle(matrix: np.ndarray):
    """
    计算三角矩阵的行列式
    :param matrix: 上三角矩阵或下三角矩阵
    :return: 上三角矩阵或下三角矩阵的行列式
    """
    h, w = matrix.shape
    sum = 1
    for i in range(h):
        sum = sum * matrix[i][i]
    return sum


def det(matrix: np.ndarray):
    """
    利用LU分解计算行列式
    :param matrix:
    :return: 行列式
    """
    P, L, U = matrix_decompose(matrix, "LU")
    print("*" * 30 + "计算行列式" + "*" * 30)
    det_p = det_P(P)
    det_l = det_triangle(L)
    det_u = det_triangle(U)
    return np.rint(det_p * det_u * det_l)


def matrix_decompose(matrix: np.ndarray, mode):
    """
    实现矩阵的LU、QR（Gram-Schmidt）、Orthogonal Reduction (Householder reduction 和Givens reduction)和 URV
    :param matrix: 需要分解的矩阵
    :param mode:
    1、LU 实现矩阵的LU分解
    2、Gram-Schmidt 采用Gram-Schmidt的方式实现矩阵的QR分解
    3、Householder reduction 采用Householder reduction的方式实现矩阵的QR分解
    4、Givens reduction 采用Givens reduction的方式实现矩阵的QR分解
    5、URV 实现矩阵的URV分解
    :return:
    若为LU分解，则一次返回P、L、U矩阵，使得PA=LU
    若为QR分解，则返回Q、R矩阵，使得A=QR
    若为URV分解，则返回R，R，U矩阵，使得A=URV^T
    """
    h, w = matrix.shape
    if h < w:
        print("*" * 30 + mode + "(h<w 类型)" + "*" * 30)
    elif h == w:
        print("*" * 30 + mode + "(h=w 类型)" + "*" * 30)
    else:
        print("*" * 30 + mode + "(h>w 类型)" + "*" * 30)
    if mode == 'LU':
        h, w = matrix.shape
        A = copy.deepcopy(matrix)
        P = np.identity(h)
        L = np.zeros((h, w))
        U = np.zeros((w, w))
        min_len = min(A.shape)
        # 高斯消元法
        for i in range(min_len):
            index = i
            for j in range(i + 1, h):
                if np.abs(A[j, i]) > np.abs(A[index, i]):
                    index = j
            temp = copy.deepcopy(P[i, :])
            P[i, :] = P[index, :]
            P[index, :] = temp
            temp = copy.deepcopy(A[i, :])
            A[i, :] = A[index, :]
            A[index, :] = temp
            print("行交换({},{}):
".format(i, index), A)
            for j in range(i + 1, h):
                c = A[j][i] / A[i][i]
                for k in range(i, w):
                    A[j][k] = A[j][k] - c * A[i][k]
                A[j][i] = c
            print("行消去({},{}以下):
{}".format(i, i, A))

        print("*A*
", A)

        for i in range(0, h):
            if i < w:
                L[i, i] = 1
            for j in range(0, min(i, w)):
                L[i, j] = A[i, j]

        for i in range(0, h):
            for j in range(i, w):
                U[i, j] = A[i, j]

        print("*" * 30 + "分解计算完成" + "*" * 30)
        print("P:
", P)
        print("L:
", L)
        print("U:
", U)
        print("*" * 30 + "验证PA=LU" + "*" * 30)
        restored_A = L @ U
        print("P@A:
", P @ matrix)
        print("restored_A:
", np.rint(restored_A))
        return P, L, U
    elif mode == 'Gram-Schmidt':
        u = []
        h, w = matrix.shape
        A = copy.deepcopy(matrix)
        R = np.zeros((w, w))
        Q = np.zeros((h, w))
        q = []
        # 施密特正交化
        for i in range(w):
            u.append(A[:, i].reshape((-1, 1)))
        for i in range(0, w):
            sum = np.zeros((h, 1))
            for j in range(0, i):
                rji = (q[j].T @ u[i])[0][0]
                R[j][i] = rji
                sum += rji * q[j]
            qi = u[i] - sum
            rii = np.linalg.norm(qi)
            R[i][i] = rii
            qi = qi / rii
            q.append(qi)

        print("h:{},w:{}".format(h, w))
        print(Q.shape)
        for i in range(w):
            Q[:, i] = q[i].reshape(h)

        print("*" * 30 + "分解计算完成" + "*" * 30)
        print("Q:
", Q)
        print("R:
", R)

        print("*" * 30 + "验证A=QR" + "*" * 30)

        A_restored = Q @ R
        A_restored[(A_restored < 1e-12) & (A_restored > -1e-12)] = 0
        print("A_restored:
", np.rint(A_restored))
        print("matrix:
", matrix)

        return Q, R

    elif mode == 'Householder reduction':
        # 计算PA=R 其中R为上三角矩阵
        A_h, A_w = matrix.shape
        Plist = []
        R = np.zeros((A_h, A_w))
        l = min(A_h, A_w)
        A = copy.deepcopy(matrix)
        if A_h > A_w:
            l = l + 1
        for i in range(l - 1):
            print("epoch:{}".format(i + 1))
            min_len = min(A.shape)
            print("A的大小为{},A的最短边大小为{}".format(A.shape, min_len))
            A_1 = A[:, 0].reshape(-1, 1)  # A的第一列
            A_1_len = A_1.shape[0]
            print("A1的长度:
", A_1_len)
            print("A的第一列:
", A_1)
            A_1_norm = np.linalg.norm(A_1)  # A的第一列的模
            print("A的第一列的模:
", A_1_norm)
            e1 = np.identity(A_1_len)[:, 0].reshape((-1, 1))
            print("e1:
", e1)
            u = A_1 - A_1_norm * e1  # u=A1-||A1||e1
            print("u=A1-||A1||e1:
", u)
            P_1 = np.identity(A_1_len) - 2 * (u @ u.T / (u.T @ u))  # P1=I-2*(UUT/UTU)
            print("P1=I-2*(UUT/UTU):
", P_1)
            P_ = np.identity(A_h)
            P_[i:, i:] = P_1  # 扩充大小
            print("P_为P左上角添加1:
", P_)
            Plist.append(P_)
            A = P_1 @ A  # 将A的第一列主元下方的元素消去
            A[(A < 1e-12) & (A > -1e-12)] = 0
            print("A:
", A)
            # 计算R矩阵
            R[i, i:] = A[0, :]
            R[i:, i] = A[:, 0]
            print("R:
", R)
            A_ = A[1:, 1:]
            print("A_:
", A_)
            A = A_
            if A_.shape[0] == 1:
                R[i + 1][i + 1:] = A[0][0:]
                R[(R < 1e-12) & (R > -1e-12)] = 0
                print("上三角矩阵R:
", R)

        P = np.identity(A_h)
        for p in Plist:
            P = p @ P
        print("P:
", P)

        Q = np.linalg.inv(P)
        Q[(Q < 1e-12) & (Q > -1e-12)] = 0
        print("*" * 30 + "分解计算完成" + "*" * 30)
        print("Q:
", Q)
        print("R:
", R)

        print("*" * 30 + "验证A=QR" + "*" * 30)

        A_restored = Q @ R
        A_restored[(A_restored < 1e-12) & (A_restored > -1e-12)] = 0
        print("A_restored:
", np.rint(A_restored))
        print("matrix:
", matrix)

        return Q, R

    elif mode == 'Givens reduction':
        h, w = matrix.shape
        Plist = []
        A = copy.deepcopy(matrix)
        # 构造旋转矩阵
        for i in range(0, h):
            for j in range(0, min(i, w)):
                P_ = np.identity(h)
                a = A[j][j]
                b = A[i][j]
                print("a:{},b:{}".format(a, b))
                if a != 0 or b != 0:
                    d = np.sqrt(a ** 2 + b ** 2)
                    c = a / d
                    s = b / d
                    print("c:{},s:{}".format(c, s))
                    P_[i][i] = c
                    P_[j][j] = c
                    P_[i][j] = -s
                    P_[j][i] = s
                Plist.append(P_)
                print("消去({},{})".format(i, j))
                print("旋转矩阵P_:
", P_)
                A = P_ @ A
                A[(A < 1e-12) & (A > -1e-12)] = 0
                print("A:
", A)

        P = np.identity(h)
        for p in Plist:
            P = p @ P
        print("P:
", P)

        Q = np.linalg.inv(P)
        R = A
        Q[(Q < 1e-12) & (Q > -1e-12)] = 0
        print("*" * 30 + "分解计算完成" + "*" * 30)
        print("Q:
", Q)
        print("R:
", R)

        print("*" * 30 + "验证A=QR" + "*" * 30)

        A_restored = Q @ R
        A_restored[(A_restored < 1e-12) & (A_restored > -1e-12)] = 0
        print("A_restored:
", np.rint(A_restored))
        print("matrix:
", matrix)

        return Q, R

    elif mode == 'URV':
        A = matrix
        m, n = A.shape
        r = np.linalg.matrix_rank(A)

        # 奇异值分解
        U, sigma, VT = np.linalg.svd(A, full_matrices=True)

        # 构造C和零矩阵
        C = np.diag(sigma[:r])
        R = np.zeros((m, n))
        R[:r, :r] = C

        # 截断U和V以获得正确的大小
        U = U[:, :m]
        V = VT.T[:, :n]

        print("*" * 30 + "验证A=URV^T" + "*" * 30)
        print((U @ R @ VT))

        return U, R, V.T
    else:
        raise ValueError("mode error")

    if mode in ['Givens reduction', 'Householder reduction', 'Gram-Schmidt']:
        print("*" * 30 + "验证QR分解" + "*" * 30)
        Q_, R_ = np.linalg.qr(matrix)
        print("Q_:
", Q_)
        print("R_:
", R_)

    elif mode == 'LU':
        print("*" * 30 + "验证LU分解" + "*" * 30)
        P_, L_, U_ = scipy.linalg.lu(matrix)
        print("P_:
", P_)
        print("L_:
", L_)
        print("U_:
", U_)


def AXb(A: np.ndarray, b: np.ndarray):
    """
    解方程AX=b
    :param A: 系数矩阵A
    :param b: b
    :return: 待求未知量X
    """
    n = A.shape[0]
    Ab = np.hstack((A, b))
    Rank_A = np.linalg.matrix_rank(A)
    Rank_Ab = np.linalg.matrix_rank(Ab)

    if Rank_A == Rank_Ab and Rank_A == n:
        P, L, U = matrix_decompose(A, "LU")
        print("*" * 30 + "解方程AX=b" + "*" * 30)
        # LUX=Pb
        Pb = P @ b
        print("P@b
", Pb)
        # LY=Pb 计算Y
        y_list = []
        for i in range(n):
            sum = 0
            for j in range(0, i):
                sum += L[i][j] * y_list[j]
            y_list.append(Pb[i][0] - sum)
        print("y_list:
", y_list)
        # UX=Y 计算X
        x_list = []
        for i in range(n - 1, -1, -1):
            sum = 0
            for j in range(n - 1, i, -1):
                sum += U[i][j] * x_list[n - j - 1]
            x_list.append((y_list[i] - sum) / U[i][i])
        x_list.reverse()
        X = np.array(x_list).reshape((-1, 1))
        print("X:
", X)
        print("*" * 30 + "验证AX=b" + "*" * 30)
        print(np.linalg.solve(A, b))
    else:
        print("不存在唯一解")


if __name__ == '__main__':
    mode = "Householder reduction"

    # # h<w 类型
    # matrix1 = np.array(
    #     [[0, -20, 23, 34, 11, 20, 2], [3, 27, 124, 52, 1, 49, 123], [4, 11, 8, 23, 9, 64, 76],
    #      [12, 45, 12, 86, 95, 19, 99]],
    #     dtype=np.float32)
    # print("matrix1:
{}
matrix1.shape:
{}".format(matrix1, matrix1.shape))
    # matrix_decompose(matrix1, mode)

    # # h=w 类型
    # matrix2 = np.array([[0, -20, 23], [3, 27, -4], [4, 11, -2]], dtype=np.float32)
    # print("matrix2:
{}
matrix2.shape:
{}".format(matrix2, matrix2.shape))
    # matrix_decompose(matrix2, mode)

    # # h>w 类型
    # matrix3 = np.array(
    #     [[0, -20, 23], [3, 27, 124], [4, 11, 8], [12, 45, 12], [123, 1, 80], [12, 99, 662]],
    #     dtype=np.float32)
    # b = np.array([[0], [0], [0], [5], [5], [5]])
    # print("matrix3:
{}
matrix3.shape:
{}".format(matrix3, matrix3.shape))
    # matrix_decompose(matrix3, b, mode)

    # LU分解
    matrix4 = np.array(
        [[0, -20, 23, 319, -89], [3, 27, -4, 987, 34], [4, 11, -2, -237, 10], [12, 67, 234, 66, -1099],
         [32, 562, 64, 88, 24]],
        dtype=np.float32)
    P, L, U = matrix_decompose(matrix4, "LU")
    print(det(matrix4))
    print(np.linalg.det(matrix4))

    # matrix6 = np.array([[1, 2, -3, 4], [4, 8, 12, -8], [2, 3, 2, 1], [-3, -1, 1, -4]], dtype=np.float32)
    # b = np.array(([1], [2], [3], [4]), dtype=np.float32)
    # AXb(matrix6, b)