PLU分解的矩阵逆答案

【问题标题】：Matrix inverse from PLU decompositionPLU分解的矩阵逆
【发布时间】：2019-09-01 17:53:24
【问题描述】：

我正在尝试从 plu 分解中计算矩阵逆，但我得到的结果是错误的。使用调试器，我找不到一个可以责备的步骤，所以这里是反转函数的代码：

template<class T>
Matrix<T> Matrix<T>::inv() const {
    std::tuple<Matrix<T>, Matrix<T>, Matrix<T>> PLU(this->decomp_PLU());
    Matrix<T> P(std::get<0>(PLU));
    Matrix<T> L(std::get<1>(PLU));
    Matrix<T> U(std::get<2>(PLU));

    if (!this->lines() == this->cols()) { throw std::logic_error("Matrix must be square"); }
    unsigned N = this->lines();

    Matrix<T> Y(N, N);
    for (unsigned i = 0; i < N; i++) {
        Matrix<T> B(Matrix<T>::gen_col(P.transp().data[i])); // B is the ith column vector of P
        Y[i] = Matrix<T>::solve_climb(L, B).transp().data[0]; // Compute LY=P for each column vector

        Matrix<T> conf(L.dot(Matrix<T>::gen_col(Y[i])));
        if (!conf.allclose(B, 1e-9, 1e-9)) {
            throw std::logic_error("WHYY");
        }
    }
    Y = Y.transp();

    Matrix<T> X(N, N);
    for (unsigned i = 0; i < N; i++) {
        Matrix<T> B(Matrix<T>::gen_col(Y.transp().data[i])); // B is the ith column vector of Y
        X[i] = Matrix<T>::solve_descent(U, B).transp().data[0]; // Compute UX=Y for each column vector

        Matrix<T> conf(U.dot(Matrix<T>::gen_col(X[i])));
        if (!conf.allclose(B, 1e-9, 1e-9)) {
            throw std::logic_error("WHYY");
        }
    }
    X = X.transp();

    return X;
}

函数Matrix<T>::gen_col 从其输入生成一个列向量，solve_climb / solve_descent 分别计算求解 AX=B 的列向量，其中 A 是一个三角矩阵（取决于具体情况，下等或上等）

仅供参考，当尝试检查每个向量的计算是否正确时，代码会抛出逻辑错误“WHYY”。

有任何线索说明代码可能出错的地方吗？

谢谢

编辑：完整代码是 here (matrix_def.h) 和 here (matrix.h)

【问题讨论】：

只是想知道，可能想要运行、测试和调试您的代码的人如何能够做到这一点，因为这里提供的唯一东西是inv 函数，没有别的。不包含头文件！ gen_col 和 solve_descent 的代码在哪里？ Y.transp() 是做什么的？ Matrix<T> X(N,N); 是做什么的？我知道，我可以假设您转置一个矩阵并使用这些函数创建一个矩阵，但就目前而言，此代码的 sn-p 不是minimal reproducible example，因此无法推测解决方案，如果并非不可能调试。
你是对的。用完整代码的链接编辑了我的问题！谢谢
这些行似乎比较了两个列向量，绝对和相对允许误差为 1e-9。是不是太小了点？你能评估一下这两者之间的距离吗？

标签： c++ algorithm math matrix matrix-inverse

【解决方案1】：

由于 L 是下三角形，你应该使用solve_descent，而不是solve_climb。 U（三角形上级）也是如此，您需要使用solve_climb。

【讨论】：