c++ - Boost 库，如何从 lu_factorize() 获取行列式？

Question

我正在尝试使用 boost C++ 库计算行列式。我找到了我在下面复制的函数 InvertMatrix() 的代码。每次我计算这个逆时，我也想要行列式。我很清楚如何通过从 LU 分解乘以 U 矩阵的对角线来计算。有一个问题，我能够正确计算行列式，除了符号。根据旋转的不同，我有一半的时间得到的符号不正确。有没有人对如何每次都获得正确的标志提出建议？提前致谢。

template<class T>
bool InvertMatrix(const ublas::matrix<T>& input, ublas::matrix<T>& inverse)
{
 using namespace boost::numeric::ublas;
 typedef permutation_matrix<std::size_t> pmatrix;
 // create a working copy of the input
 matrix<T> A(input);
 // create a permutation matrix for the LU-factorization
 pmatrix pm(A.size1());

 // perform LU-factorization
 int res = lu_factorize(A,pm);
 if( res != 0 ) return false;

这里是我在计算行列式时插入的最佳镜头。

 T determinant = 1;

 for(int i = 0; i < A.size1(); i++)
 {
  determinant *= A(i,i);
 }

结束我的代码部分。

 // create identity matrix of "inverse"
 inverse.assign(ublas::identity_matrix<T>(A.size1()));

 // backsubstitute to get the inverse
 lu_substitute(A, pm, inverse);

 return true;
}

Answer 1

置换矩阵 pm 包含确定符号变化所需的信息:您需要将行列式乘以置换矩阵的行列式。

细读源文件 lu.hpp我们找到一个名为 swap_rows 的函数，它告诉我们如何将置换矩阵应用于矩阵。它很容易修改以产生置换矩阵的行列式(置换的符号)，假设每个实际交换贡献一个因子 -1:

template <typename size_type, typename A>
int determinant(const permutation_matrix<size_type,A>& pm)
{
    int pm_sign=1;
    size_type size=pm.size();
    for (size_type i = 0; i < size; ++i)
        if (i != pm(i))
            pm_sign* = -1; // swap_rows would swap a pair of rows here, so we change sign
    return pm_sign;
}

另一种选择是使用不进行任何旋转的 lu_factorize 和 lu_substitute 方法(请参阅源代码，但基本上删除 pm 在对 lu_factorize 和 lu_substitute 的调用中)。该更改将使您的行列式计算按原样工作。但是要小心:删除主元会使算法的数值稳定性降低。