矩阵列主元三角分解

一，列主元三角分解定理

如果A为非奇异矩阵，则存在排列矩阵P使 PA=LU,其中L为下三角矩阵，U为上三角矩阵，即A=P^-1LU

二，列主元三角分解Python代码

# 自己原创
def pivot_lu_decomposition(coefficient_matrix: np.ndarray, right_hand_side_vector: np.ndarray):
    """
    实现方程Ax=b系数矩阵A的pivoted LU decomposition
    :param coefficient_matrix: 初始系数矩阵A
    :param right_hand_side_vector: 初始常数列向量b
    :return: 排列矩阵P，单位下三角矩阵L，上三角矩阵U，常数列向量b
    """
    # first step: evaluate the lu decomposition condition
    rows, columns = coefficient_matrix.shape
    if rows == columns:  # judge if it is a square matrix
        for k in range(rows):  # 判断各阶顺序主子式是否为0
            if det(coefficient_matrix[:k + 1, :k + 1]) == 0:
                raise Exception("cannot generate LU decomposition")
        else:
            # pivot LU decomposition
            lower_triangular_matrix = np.eye(rows)
            permutation_matrix = np.eye(rows)
            for k in range(rows - 1):
                max_index = np.argmax(abs(coefficient_matrix[k:, k]))
                coefficient_matrix[[k, max_index + k], :] = coefficient_matrix[[max_index + k, k], :]
                right_hand_side_vector[[k, max_index + k], :] = right_hand_side_vector[[max_index + k, k], :]
                permutation_matrix[[k, max_index + k], :] = permutation_matrix[[max_index + k, k], :]
                for row in range(k + 1, rows):
                    multiplier = coefficient_matrix[row, k] / coefficient_matrix[k, k]
                    coefficient_matrix[row, k:] += -multiplier * (coefficient_matrix[k, k:])
                    right_hand_side_vector[row] += -multiplier * right_hand_side_vector[k]
                    lower_triangular_matrix[row, k] = multiplier
    else:
        raise Exception("ERROR:please pass a square matrix.")

    return permutation_matrix, lower_triangular_matrix, coefficient_matrix, right_hand_side_vector